1program midas_ensembleH
2 !
3 !:Purpose: Main program for applying the observation operator to an ensemble
4 ! of states to compute the innovations and store them in eob*_HX binary files. The
5 ! ``LETKF`` program can read these binary files instead of computing the innovations.
6 !
7 ! ---
8 !
9 !:Algorithm: The background (a.k.a trial) ensemble members are read and the non-linear
10 ! observation operators are applied to each ensemble member:
11 ! :math:`H(xb_{i})`,
12 ! the innovations are computed:
13 ! :math:`y-H(xb_{i})`,
14 ! and stored in unformatted binary files. ``ensembleH`` has the ability
15 ! to compute innovations for a batch of ensemble members when the member
16 ! indices are sequential.
17 !
18 ! --
19 !
20 !============================================== ==============================================================
21 ! Input and Output Files Description of file
22 !============================================== ==============================================================
23 ! ``flnml`` In - Main namelist file with parameters user may modify
24 ! ``flnml_static`` In - The "static" namelist that should not be modified
25 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
26 ! ``ensemble/$YYYYMMDDHH_006_$NNNN`` In - Background ensemble member files
27 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
28 ! ``obserr`` In - Observation error statistics
29 ! ``$YYYYMMDDHH_000_$NNNN`` Out - Analysis ensemble member files
30 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
31 ! Remainder are files related to radiance obs:
32 ! ``stats_$SENSOR_assim`` In - Satellite radiance observation errors of difference sensors
33 ! ``stats_tovs`` In - Satellite radiance observation errors
34 ! ``stats_tovs_symmetricObsErr`` In - User-defined symmetric TOVS errors for all sky
35 ! ``Cmat_$PLATFORM_$SENSOR.dat`` In - Inter-channel observation-error correlations
36 ! ``ceres_global.std`` In - High-res surface type and water fraction for radiance obs
37 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
38 ! ``rttov_h2o_limits.dat`` In - Min/max humidity limits applied to analysis
39 ! ``ozoneclim98`` In - Ozone climatology
40 !============================================== ==============================================================
41 !
42 ! --
43 !
44 !:Synopsis: Below is a summary of the ``ensembleH`` program calling sequence:
45 !
46 ! - **Initial setups:**
47 !
48 ! - Read the NAMENSEMBLEH namelist and check/modify some values.
49 !
50 ! - Various modules are setup: ``obsFiles_mod``,
51 ! ``gridStateVector_mod``, ``timeCoord_mod`` (and set up dates
52 ! and ``dateStampList`` variables for both trials and
53 ! increments/analyses).
54 !
55 ! - Setup horizontal and vertical grid objects from first
56 ! ensemble member file.
57 !
58 ! - Setup ``obsSpaceData`` object and read observations from
59 ! files: ``inn_setupObs``.
60 !
61 ! - Setup the observation error statistics in ``obsSpaceData``
62 ! object: ``oer_setObsErrors``.
63 !
64 ! - Allocate objects for ``column_mod``.
65 !
66 ! - Allocate and some setup of objects for
67 ! ``ensembleObservations_mod``.
68 !
69 ! - Allocate some objects for ``gridStateVector_mod``.
70 !
71 ! - Allocate ensemble object and read background ensemble members:
72 ! ``ens_readEnsemble``.
73 !
74 ! - **Computation**
75 !
76 ! - Option to read ensemble mean from file (``gio_readFromFile``) or
77 ! compute ensemble mean (``ens_computeMean``)
78 !
79 ! - Loop over background ensemble members, computing innovation for each,
80 ! with resulting ``H(xb)`` being stored in ``ensObs`` objects both for
81 ! original ensemble and, optionally, for the modulated
82 ! ensemble members.
83 !
84 ! - store the local ``ensObs`` object to binary file.
85 !
86 ! --
87 !
88 !
89 !:Options: `List of namelist blocks <../namelists_in_each_program.html#ensembleH>`_
90 ! that can affect the ``ensembleH`` program.
91 !
92 ! * The use of ``ensembleH`` program is controlled by the namelist block
93 ! ``&NAMENSEMBLEH`` read by the ``ensembleH`` program.
94 !
95 ! * Some of the other relevant namelist blocks used to configure the
96 ! ``ensembleH`` are listed in the following table:
97 !
98 !===================== ================== ===============================================================
99 ! Program/Module Namelist Description of what is controlled
100 !===================== ================== ===============================================================
101 ! ``midas_ensembleh`` ``NAMENSEMBLEH`` number of ensemble members, additional parameters to control
102 ! generating modulated members from original ensembles,
103 ! the member number the current batch starts with, number of
104 ! members in the batch, option to read ensemble mean from file.
105 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal resolution of
106 ! the background state.
107 !===================== ================== ===============================================================
108 !
109 use version_mod
110 use midasMpi_mod
111 use fileNames_mod
112 use gridStateVector_mod
113 use gridStateVectorFileIO_mod
114 use columnData_mod
115 use verticalCoord_mod
116 use horizontalCoord_mod
117 use timeCoord_mod
118 use utilities_mod
119 use ramDisk_mod
120 use stateToColumn_mod
121 use obsFiles_mod
122 use obsSpaceData_mod
123 use obsErrors_mod
124 use innovation_mod
125 use ensembleObservations_mod
126 use ensembleStateVector_mod
127 use enkf_mod
128 implicit none
129
130 type(struct_obs), target :: obsSpaceData
131 type(struct_ens) :: ensembleTrl4D
132 type(struct_gsv) :: stateVectorMeanTrl4D
133 type(struct_gsv) :: stateVector4D
134 type(struct_gsv) :: stateVector4Dmod
135 type(struct_gsv) :: stateVectorWithZandP4D
136 type(struct_gsv) :: stateVectorHeightSfc
137 type(struct_columnData) :: column
138
139 type(struct_eob), target :: ensObs
140 type(struct_eob), pointer :: ensObsGain
141
142 type(struct_vco), pointer :: vco_ens => null()
143 type(struct_hco), pointer :: hco_ens => null()
144
145 integer :: get_max_rss, fclos, fnom, fstopc, ierr
146 integer :: memberIndex, nulnam, dateStampFromObs
147 integer :: nEnsGain, eigenVectorIndex, memberIndexInEnsObs, stepIndex
148 integer, allocatable :: dateStampList(:)
149
150 logical :: useModulatedEns, fileExists
151
152 character(len=256) :: ensFileName
153 character(len=256) :: ensMeanFileName
154 character(len=9) :: obsColumnMode
155 character(len=48) :: obsMpiStrategy
156 character(len=48) :: midasMode
157
158 ! namelist variables
159 character(len=256) :: ensPathName ! path of ensemble member files
160 character(len=20) :: obsTimeInterpType ! type of time interpolation to obs time
161 integer :: nEns ! ensemble size
162 integer :: numRetainedEigen ! number of retained eigen modes used for modulated ensemble
163 integer :: fileMemberIndex1 ! first member index in ensemble set to be read
164 integer :: numFullEns ! number of full ensemble set (needed only for modulated ensemble)
165 real(8) :: vLocalize ! vertical localization lengthscale (needed only for modulated ensemble)
166 logical :: readEnsMeanFromFile ! choose to read ens mean from file (when reading subset of members)
167 NAMELIST /NAMENSEMBLEH/nEns, ensPathName, obsTimeInterpType, numRetainedEigen, &
168 vLocalize, fileMemberIndex1, readEnsMeanFromFile, numFullEns
169
170 midasMode = 'analysis'
171 obsColumnMode = 'ENKFMIDAS'
172 obsMpiStrategy = 'LIKESPLITFILES'
173
174 call ver_printNameAndVersion('ensembleH','Program for applying H to ensemble')
175
176 !
177 !- 0. MPI, TMG initialization
178 !
179 call mmpi_initialize
180 call tmg_init(mmpi_myid, 'TMG_INFO')
181
182 call utl_tmg_start(0,'Main')
183 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
184
185 ! Avoid printing lots of stuff to listing for std file I/O
186 ierr = fstopc('MSGLVL','ERRORS',0)
187
188 ! Setup the ramdisk directory (if supplied)
189 call ram_setup
190
191 ! Setting default namelist variable values
192 nEns = 10
193 ensPathName = 'ensemble'
194 obsTimeInterpType = 'LINEAR'
195 numRetainedEigen = 0
196 vLocalize = -1.0D0
197 fileMemberIndex1 = 1
198 readEnsMeanFromFile = .false.
199 numFullEns = 0
200
201 ! Read the namelist
202 nulnam = 0
203 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
204 read(nulnam, nml=namensembleh, iostat=ierr)
205 if (ierr /= 0) call utl_abort('midas-ensembleH: Error reading namelist')
206 if (mmpi_myid == 0) write(*,nml=namensembleh)
207 ierr = fclos(nulnam)
208
209 if (numRetainedEigen < 0) call utl_abort('midas-ensembleH: numRetainedEigen should be ' // &
210 'equal or greater than zero')
211
212 useModulatedEns = (numRetainedEigen > 0)
213 if (useModulatedEns .and. vLocalize <= 0) then
214 call utl_abort('midas-ensembleH: vLocalize should be greater than zero for modulated ens')
215 end if
216
217 if (useModulatedEns .and. numFullEns < nEns) then
218 call utl_abort('midas-ensembleH: For modulated ensembles the number of full ensembles is needed ' // &
219 'with numFullEns >= nEns')
220 end if
221
222 ! Read the observations
223 call obsf_setup(dateStampFromObs, midasMode)
224
225 ! Use the first ensemble member to initialize datestamp and grid
226 call fln_ensFileName(ensFileName, ensPathName, memberIndex_opt=1, &
227 fileMemberIndex1_opt=fileMemberIndex1)
228
229 ! Setup timeCoord module, get datestamp from ensemble member
230 call tim_setup(fileNameForDate_opt = ensFileName)
231 allocate(dateStampList(tim_nstepobs))
232 call tim_getstamplist(dateStampList,tim_nstepobs,tim_getDatestamp())
233
234 write(*,*) 'midas-ensembleH: dateStamp = ',tim_getDateStamp()
235
236 !- Initialize variables of the model states
237 call gsv_setup
238
239 !- Initialize the Ensemble grid
240 if (mmpi_myid == 0) write(*,*) ''
241 if (mmpi_myid == 0) write(*,*) 'midas-ensembleH: Set hco and vco parameters for ensemble grid'
242 call hco_SetupFromFile(hco_ens, ensFileName, ' ', 'ENSFILEGRID')
243 call vco_setupFromFile(vco_ens, ensFileName)
244
245 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
246
247 ! read in the observations
248 call inn_setupObs(obsSpaceData, hco_ens, obsColumnMode, obsMpiStrategy, midasMode)
249
250 ! Initialize obs error covariances
251 call oer_setObsErrors(obsSpaceData, midasMode)
252
253 call col_setup
254 call col_setVco(column, vco_ens)
255 call col_allocate(column, obs_numheader(obsSpaceData))
256 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
257
258 ! Initialize the observation error covariances
259 call oer_setObsErrors(obsSpaceData, midasMode) ! IN
260
261 ! Allocate and initialize eob object for storing HX values
262 call eob_allocate(ensObs, nEns, obs_numBody(obsSpaceData), obsSpaceData, &
263 fileMemberIndex1_opt=fileMemberIndex1)
264 call eob_zero(ensObs)
265 if (useModulatedEns) then
266 nEnsGain = nEns * numRetainedEigen
267 allocate(ensObsGain)
268 call eob_allocate(ensObsGain, nEnsGain, obs_numBody(obsSpaceData), obsSpaceData, &
269 fileMemberIndex1_opt=fileMemberIndex1)
270 call eob_zero(ensObsGain)
271 else
272 ensObsGain => ensObs
273 end if
274
275 ! Set lat, lon, obs values in ensObs
276 call eob_setLatLonObs(ensObs)
277 if (useModulatedEns) call eob_setLatLonObs(ensObsGain)
278
279
280 ! Read the sfc height from ensemble member 1
281 call gsv_allocate(stateVectorHeightSfc, 1, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
282 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
283 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','TT'/))
284 call gio_readFromFile(stateVectorHeightSfc, ensFileName, ' ', ' ', &
285 containsFullField_opt=.true., readHeightSfc_opt=.true.)
286
287 ! Allocate statevector related to ensemble mean
288 call gsv_allocate(stateVectorMeanTrl4D, tim_nstepobs, hco_ens, vco_ens, &
289 dateStamp_opt=tim_getDateStamp(), &
290 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
291 dataKind_opt=4, allocHeightSfc_opt=.true., &
292 allocHeight_opt=.false., allocPressure_opt=.false.)
293 call gsv_zero(stateVectorMeanTrl4D)
294 call gsv_allocate(stateVector4D, tim_nstepobs, hco_ens, vco_ens, &
295 dateStamp_opt=tim_getDateStamp(), &
296 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
297 dataKind_opt=4, allocHeightSfc_opt=.true., &
298 allocHeight_opt=.false., allocPressure_opt=.false.)
299 call gsv_zero(stateVector4D)
300 if (useModulatedEns) then
301 ! same as stateVector4D
302 call gsv_allocate(stateVector4Dmod, tim_nstepobs, hco_ens, vco_ens, &
303 dateStamp_opt=tim_getDateStamp(), &
304 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
305 dataKind_opt=4, allocHeightSfc_opt=.true., &
306 allocHeight_opt=.false., allocPressure_opt=.false.)
307 call gsv_zero(stateVector4Dmod)
308 end if
309
310 ! Allocate statevector for storing state with heights and pressures allocated (for s2c_nl)
311 call gsv_allocate(stateVectorWithZandP4D, tim_nstepobs, hco_ens, vco_ens, &
312 dateStamp_opt=tim_getDateStamp(), &
313 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
314 dataKind_opt=4, allocHeightSfc_opt=.true.)
315 call gsv_zero(stateVectorWithZandP4D)
316
317 ! Allocate ensembles, read the Trl ensemble
318 call utl_tmg_start(2,'--ReadEnsemble')
319 call ens_allocate(ensembleTrl4D, nEns, tim_nstepobs, hco_ens, vco_ens, dateStampList, &
320 fileMemberIndex1_opt=fileMemberIndex1)
321 call ens_readEnsemble(ensembleTrl4D, ensPathName, biPeriodic=.false.)
322 call utl_tmg_stop(2)
323
324 ! Compute ensemble mean and copy to meanTrl stateVectors
325 if (readEnsMeanFromFile) then
326 ! read the mean stateVector and copy to ensembleTrl4D object
327 call fln_ensTrlFileName(ensMeanFileName, '.', tim_getDateStamp())
328 ensMeanFileName = trim(ensMeanFileName) // '_trialmean'
329 inquire(file=trim(ensMeanFileName),exist=fileExists)
330 if (.not. fileExists) then
331 call utl_abort('midas-ensembleH: the ensemble mean file does not exist')
332 end if
333 do stepIndex = 1, tim_nstepobs
334 call gio_readFromFile(stateVectorMeanTrl4D, ensMeanFileName, ' ', ' ', &
335 containsFullField_opt=.true., readHeightSfc_opt=.true., &
336 stepIndex_opt=stepIndex)
337 end do
338 call ens_copyToEnsMean(ensembleTrl4D, stateVectorMeanTrl4D)
339 else
340 call ens_computeMean(ensembleTrl4D)
341 call ens_copyEnsMean(ensembleTrl4D, stateVectorMeanTrl4D)
342 end if
343
344 do memberIndex = 1, nEns
345
346 write(*,*) ''
347 write(*,*) 'midas-ensembleH: apply nonlinear H to ensemble member ', memberIndex
348 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
349
350 ! copy 1 member to a stateVector
351 call ens_copyMember(ensembleTrl4D, stateVector4D, memberIndex)
352 call gsv_copy(stateVector4D, stateVectorWithZandP4D, allowVarMismatch_opt=.true., &
353 beSilent_opt=.true.)
354 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
355
356 ! Compute and set Yb in ensObs
357 call s2c_nl(stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
358 timeInterpType=obsTimeInterpType, dealloc_opt=.false., &
359 beSilent_opt=.true.)
360
361 ! Compute Y-H(X) in OBS_OMP
362 call inn_computeInnovation(column, obsSpaceData, beSilent_opt=.true.)
363
364 ! Copy to ensObs: Y-HX for this member
365 call eob_setYb(ensObs, memberIndex)
366
367 ! Compute and set Yb in ensObsGain
368 do eigenVectorIndex = 1, numRetainedEigen
369 if (mmpi_myid == 0) write(*,*) 'midas-ensembleH: apply nonlinear H to modulated member ', &
370 eigenVectorIndex, '/', numRetainedEigen
371
372 ! modulate the member with eigenvectors of vertical localization matrix
373 call enkf_getModulatedState(stateVector4D, stateVectorMeanTrl4D, &
374 vLocalize, numRetainedEigen, nEns, &
375 eigenVectorIndex, stateVector4Dmod, &
376 beSilent=.true.)
377
378 call gsv_copy(stateVector4Dmod, stateVectorWithZandP4D, allowVarMismatch_opt=.true., &
379 beSilent_opt=.true.)
380 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
381
382 call s2c_nl(stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
383 timeInterpType=obsTimeInterpType, dealloc_opt=.false., &
384 beSilent_opt=.true.)
385
386 ! Compute Y-H(X) in OBS_OMP
387 call inn_computeInnovation(column, obsSpaceData, filterObsAndInitOer_opt=.false., &
388 beSilent_opt=.true.)
389
390 ! Copy to ensObsGain: Y-HX for this member
391 memberIndexInEnsObs = (eigenVectorIndex - 1) * nEns + memberIndex
392 call eob_setYb(ensObsGain, memberIndexInEnsObs)
393 end do ! eigenVectorIndex
394
395 end do
396 call gsv_deallocate(stateVectorWithZandP4D)
397 if (gsv_isAllocated(stateVector4Dmod)) call gsv_deallocate(stateVector4Dmod)
398 call gsv_deallocate(stateVector4D)
399 call gsv_deallocate(stateVectorMeanTrl4D)
400 call gsv_deallocate(stateVectorHeightSfc)
401
402 ! write local ensObs to file
403 call eob_writeToFiles(ensObs, outputFilenamePrefix='eob_HX', writeObsInfo=.true.)
404 if (useModulatedEns) then
405 call eob_writeToFiles(ensObsGain, outputFilenamePrefix='eobGain_HX', writeObsInfo=.false., &
406 numGroupsToDivideMembers_opt=numRetainedEigen, &
407 maxNumMembersPerGroup_opt=numFullEns)
408 end if
409
410 !
411 !- MPI, tmg finalize
412 !
413 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
414 call utl_tmg_stop(0)
415
416 call tmg_terminate(mmpi_myid, 'TMG_INFO')
417 call rpn_comm_finalize(ierr)
418
419 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
420
421 !
422 !- 7. Ending
423 !
424 if (mmpi_myid == 0) write(*,*) ' --------------------------------'
425 if (mmpi_myid == 0) write(*,*) ' MIDAS-ENSEMBLEH ENDS'
426 if (mmpi_myid == 0) write(*,*) ' --------------------------------'
427
428end program midas_ensembleH
