1program midas_letkf
2 !
3 !:Purpose: Main program for the local ensemble transform Kalman filter
4 ! (LETKF). Several different variations of the LETKF algorithm have
5 ! been implemented. Note that the actual calculation of the analyses
6 ! is in the subroutine ``enkf_LETKFanalyses``. Many aspects of this
7 ! program are controlled through the namelist block NAMLETKF.
8 !
9 ! ---
10 !
11 !:Algorithm: The ``letkf`` programs implements several variations of the LETKF
12 ! ensemble data assimilation algorithm. The following variations
13 ! can be chosen through the namelist:
14 !
15 ! - **LETKF**: standard LETKF
16 !
17 ! - **CVLETKF**: LETKF with deterministic approach to cross validation using
18 ! the "gain form" (currently operational)
19 !
20 ! - **CVLETKF-PERTOBS**: LETKF with stochastic approach to cross validation
21 !
22 ! - **LETKF-Gain**: standard LETKF, but using the "gain form"
23 !
24 ! - **LETKF-Gain-ME**: standard LETKF, but using the "modulated ensemble" and "gain
25 ! form" to account for vertical localization
26 !
27 ! - **CVLETKF-ME**: LETKF with deterministic approach to cross validation using
28 ! the "modulated ensemble" and "gain form" to account for
29 ! vertical localization
30 !
31 ! More detail on the LETKF algorithms with and without cross
32 ! validation (deterministic and stochastic) can be found in the
33 ! paper:
34 ! `LETKF with cross validation <https://journals.ametsoc.org/view/journals/mwre/148/6/MWR-D-19-0402.1.xml>`_
35 !
36 ! --
37 !
38 !============================================== ==============================================================
39 ! Input and Output Files (NWP applicaton) Description of file
40 !============================================== ==============================================================
41 ! ``flnml`` In - Main namelist file with parameters user may modify
42 ! ``flnml_static`` In - The "static" namelist that should not be modified
43 ! ``ensemble/$YYYYMMDDHH_006_$NNNN`` In - Background ensemble member files
44 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
45 ! ``obserr`` In - Observation error statistics
46 ! ``$YYYYMMDDHH_000_$NNNN`` Out - Analysis ensemble member files
47 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
48 ! Remainder are files related to radiance obs:
49 ! ``stats_$SENSOR_assim`` In - Satellite radiance observation errors of difference sensors
50 ! ``stats_tovs`` In - Satellite radiance observation errors
51 ! ``stats_tovs_symmetricObsErr`` In - User-defined symmetric TOVS errors for all sky
52 ! ``Cmat_$PLATFORM_$SENSOR.dat`` In - Inter-channel observation-error correlations
53 ! ``ceres_global.std`` In - High-res surface type and water fraction for radiance obs
54 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
55 ! ``rttov_h2o_limits.dat`` In - Min/max humidity limits applied to analysis
56 ! ``ozoneclim98`` In - Ozone climatology
57 !============================================== ==============================================================
58 !
59 ! --
60 !
61 !:Synopsis: Below is a summary of the ``letkf`` program calling sequence:
62 !
63 ! - **Initial setups:**
64 !
65 ! - Read the NAMLETKF namelist and check/modify some values.
66 !
67 ! - Various modules are setup: ``obsFiles_mod``,
68 ! ``gridStateVector_mod``, ``timeCoord_mod`` (and set up dates
69 ! and ``dateStampList`` variables for both trials and
70 ! increments/analyses).
71 !
72 ! - Setup horizontal and vertical grid objects from first
73 ! ensemble member file and determine if this is an NWP or
74 ! ocean application.
75 !
76 ! - Setup ``obsSpaceData`` object and read observations from
77 ! files: ``inn_setupObs``.
78 !
79 ! - Setup the observation error statistics in ``obsSpaceData``
80 ! object: ``oer_setObsErrors``.
81 !
82 ! - Allocate and some setup of objects for
83 ! ``ensembleObservations_mod``.
84 !
85 ! - Allocate objects for ``column_mod`` and
86 ! ``gridStateVector_mod``.
87 !
88 ! - Allocate ensemble object and read trial ensemble:
89 ! ``ens_readEnsemble``.
90 !
91 ! - Optionally, read a deterministic state for recentering the
92 ! trial ensemble before the analysis (for special applications
93 ! using LETKF for deterministic analysis).
94 !
95 ! - Compute ensemble mean: ``ens_computeMean``.
96 !
97 ! - **LETKF computations:**
98 !
99 ! - Loop over trial members, computing innovation for each, with
100 ! resulting :math:`H(xb)` being stored in ``ensObs`` objects both for
101 ! original ensemble and, optionally, for the modulated
102 ! ensemble members.
103 !
104 ! - Compute some additional quantities in ``ensObs`` objects based
105 ! on :math:`H(xb)` values.
106 !
107 ! - Compute :math:`y-H(xb)` for trial ensemble mean.
108 !
109 ! - Additional observation quality control procedures based on
110 ! quantities computed from trial ensemble.
111 !
112 ! - Gather ``ensObs`` quantities from all MPI tasks onto all MPI
113 ! tasks.
114 !
115 ! - Allocate and prepare objects for ``ensembleStateVector_mod``
116 ! used to store the trial and analysis ensembles with
117 ! temporal resolution of analysis (can be same as temporal
118 ! resolution for innovation calculation or only a single time
119 ! step).
120 !
121 ! - Setup information for interpolating weights from coarse grid
122 ! to the full model grid: ``enkf_setupInterpInfo``.
123 !
124 ! - Perform LETKF analysis to compute the analysis ensemble:
125 ! ``enkf_LETKFanalyses``
126 !
127 ! - **Final steps:**
128 !
129 ! - Optionally, compute :math:`H(xa)` for each member and put values in
130 ! ``ensObs`` object for later output to diag SQLite files.
131 !
132 ! - Write out (update) the observation files: ``obsf_writeFiles``.
133 !
134 ! - Compute :math:`y-H(xa)` for ensemble mean analysis, stored in
135 ! ``obsSpaceData``.
136 !
137 ! - Optionally, do post-processing of the analysis ensemble (same
138 ! functionality as the ``ensPostProcess`` program), or just
139 ! write out the raw analysis ensemble (for later processing by
140 ! ``ensPostProcess``).
141 !
142 ! --
143 !
144 !:Options: `List of namelist blocks <../namelists_in_each_program.html#letkf>`_
145 ! that can affect the ``letkf`` program.
146 !
147 ! * Some of the other relevant namelist blocks used to configure the
148 ! letkf analysis are listed in the following table:
149 !
150 !======================== ============== ==============================================================
151 ! Program/Module Namelist Description of what is controlled
152 !======================== ============== ==============================================================
153 ! ``midas_letkf`` ``NAMLETKF`` LETKF algorithm, number of ensemble members and
154 ! additional parameters for controlling the LETKF analysis
155 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal resolution of
156 ! the background state and the analysis
157 !======================== ============== ==============================================================
158 !
159 use version_mod
160 use midasMpi_mod
161 use mathPhysConstants_mod
162 use fileNames_mod
163 use ensembleObservations_mod
164 use ensembleStateVector_mod
165 use gridStateVector_mod
166 use gridStateVectorFileIO_mod
167 use columnData_mod
168 use tovsNL_mod
169 use verticalCoord_mod
170 use horizontalCoord_mod
171 use oceanMask_mod
172 use timeCoord_mod
173 use obsTimeInterp_mod
174 use utilities_mod
175 use ramDisk_mod
176 use stateToColumn_mod
177 use obsFiles_mod
178 use obsSpaceData_mod
179 use obsErrors_mod
180 use obsFilter_mod
181 use innovation_mod
182 use enkf_mod
183 use ensPostProcess_mod
184 implicit none
185
186 type(struct_obs), target :: obsSpaceData
187 type(struct_ens), target :: ensembleTrl4D
188 type(struct_ens), pointer :: ensembleTrl
189 type(struct_ens) :: ensembleAnl
190 type(struct_gsv) :: stateVectorMeanTrl4D
191 type(struct_gsv) :: stateVectorMemberAnl
192 type(struct_gsv) :: stateVectorMeanAnl
193 type(struct_gsv) :: stateVector4D
194 type(struct_gsv) :: stateVector4Dmod
195 type(struct_gsv) :: stateVectorWithZandP4D
196 type(struct_gsv) :: stateVectorHeightSfc
197 type(struct_gsv) :: stateVectorCtrlTrl
198 type(struct_gsv) :: stateVectorRecenter
199 type(struct_columnData) :: column
200
201 type(struct_eob), target :: ensObs, ensObs_mpiglobal
202 type(struct_eob), pointer :: ensObsGain, ensObsGain_mpiglobal
203
204 type(struct_vco), pointer :: vco_ens => null()
205 type(struct_hco), pointer :: hco_ens => null()
206 type(struct_ocm) :: oceanMask
207
208 integer :: memberIndex, middleStepIndex, stepIndex, randomSeedObs
209 integer :: nulnam, dateStampFromObs, ierr
210 integer :: get_max_rss, fclos, fnom, fstopc
211 integer :: nEnsGain, eigenVectorIndex, memberIndexInEnsObs
212 integer, allocatable :: dateStampList(:), dateStampListInc(:)
213
214 character(len=256) :: ensFileName, ctrlFileName, recenterFileName
215 character(len=9) :: obsColumnMode
216 character(len=48) :: obsMpiStrategy
217 character(len=48) :: midasMode
218
219 logical :: nwpFields ! indicates if fields are on momentum and thermo levels
220 logical :: oceanFields ! indicates if fields are on depth levels
221 logical :: useModulatedEns ! using modulated ensembles is requested by setting numRetainedEigen
222
223 ! interpolation information for weights (in enkf_mod)
224 type(struct_enkfInterpInfo) :: wInterpInfo
225
226 real(4), pointer :: field_Psfc(:,:,:,:)
227
228 ! namelist variables
229 character(len=20) :: algorithm ! name of the chosen LETKF algorithm: 'LETKF', 'CVLETKF'
230 logical :: ensPostProcessing ! do all post-processing of analysis ensemble
231 logical :: recenterInputEns ! read a deterministic state to recenter ensemble
232 integer :: numSubEns ! number of sub-ensembles to split the full ensemble
233 character(len=256) :: ensPathName ! absolute or relative path to ensemble directory
234 integer :: nEns ! ensemble size
235 logical :: randomShuffleSubEns ! choose to randomly shuffle members into subensembles
236 logical :: writeLocalEnsObsToFile ! Controls writing the ensObs to file.
237 integer :: maxNumLocalObs ! maximum number of obs in each local volume to assimilate
238 integer :: weightLatLonStep ! separation of lat-lon grid points for weight calculation
239 integer :: numRetainedEigen ! number of retained eigenValues/Vectors of vertical localization matrix
240 logical :: modifyAmsubObsError ! reduce AMSU-B obs error stddev in tropics
241 logical :: backgroundCheck ! apply additional background check using ensemble spread
242 logical :: huberize ! apply huber norm quality control procedure
243 logical :: rejectHighLatIR ! reject all IR observations at high latitudes
244 logical :: rejectRadNearSfc ! reject radiance observations near the surface
245 logical :: ignoreEnsDate ! when reading ensemble, ignore the date
246 logical :: outputOnlyEnsMean ! when writing ensemble, can choose to only write member zero
247 logical :: outputEnsObs ! to write trial and analysis ensemble members in observation space to sqlite
248 logical :: debug ! debug option to print values to the listings.
249 logical :: readEnsObsFromFile ! instead of computing innovations, read ensObs%Yb from file.
250 real(8) :: hLocalize(4) ! horizontal localization radius (in km)
251 real(8) :: hLocalizePressure(3) ! pressures where horizontal localization changes (in hPa)
252 real(8) :: vLocalize ! vertical localization radius (units: ln(Pressure in Pa) or meters)
253 real(8) :: minDistanceToLand ! for ice/ocean DA: minimum distance to land for assimilating obs
254 character(len=20) :: obsTimeInterpType ! type of time interpolation to obs time
255 character(len=20) :: mpiDistribution ! type of mpiDistribution for weight calculation ('ROUNDROBIN' or 'TILES')
256 character(len=12) :: etiket_anl ! etiket for output files
257
258 NAMELIST /NAMLETKF/algorithm, ensPostProcessing, recenterInputEns, nEns, numSubEns, &
259 ensPathName, randomShuffleSubEns, &
260 hLocalize, hLocalizePressure, vLocalize, minDistanceToLand, &
261 maxNumLocalObs, weightLatLonStep, &
262 modifyAmsubObsError, backgroundCheck, huberize, rejectHighLatIR, rejectRadNearSfc, &
263 ignoreEnsDate, outputOnlyEnsMean, outputEnsObs, &
264 obsTimeInterpType, mpiDistribution, etiket_anl, &
265 readEnsObsFromFile, writeLocalEnsObsToFile, &
266 numRetainedEigen, debug
267
268 ! Some high-level configuration settings
269 midasMode = 'analysis'
270 obsColumnMode = 'ENKFMIDAS'
271 obsMpiStrategy = 'LIKESPLITFILES'
272
273 call ver_printNameAndVersion('letkf','Program for Local Ensemble Transform Kalman Filter')
274
275 !
276 !- 0. MPI, TMG and misc. initialization
277 !
278 call mmpi_initialize
279 call tmg_init(mmpi_myid, 'TMG_INFO')
280
281 call utl_tmg_start(0,'Main')
282
283 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
284
285 ! Avoid printing lots of stuff to listing for std file I/O
286 ierr = fstopc('MSGLVL','ERRORS',0)
287
288 ! Setup the ramdisk directory (if supplied)
289 call ram_setup
290
291 !
292 !- 1. Set/Read values for the namelist NAMLETKF
293 !
294
295 !- 1.1 Setting default namelist variable values
296 algorithm = 'LETKF'
297 ensPostProcessing = .false.
298 recenterInputEns = .false.
299 ensPathName = 'ensemble'
300 nEns = 10
301 numSubEns = 2
302 randomShuffleSubEns = .false.
303 maxNumLocalObs = 1000
304 weightLatLonStep = 1
305 modifyAmsubObsError = .false.
306 backgroundCheck = .false.
307 huberize = .false.
308 rejectHighLatIR = .false.
309 rejectRadNearSfc = .false.
310 ignoreEnsDate = .false.
311 outputOnlyEnsMean = .false.
312 outputEnsObs = .false.
313 hLocalize(:) = -1.0D0
314 hLocalizePressure = (/14.0D0, 140.0D0, 400.0D0/)
315 vLocalize = -1.0D0
316 minDistanceToLand = -1.0D0
317 obsTimeInterpType = 'LINEAR'
318 mpiDistribution = 'ROUNDROBIN'
319 etiket_anl = 'ENS_ANL'
320 readEnsObsFromFile = .false.
321 writeLocalEnsObsToFile = .false.
322 numRetainedEigen = 0
323 debug = .false.
324
325 !- 1.2 Read the namelist
326 nulnam = 0
327 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
328 read(nulnam, nml=namletkf, iostat=ierr)
329 if ( ierr /= 0) call utl_abort('midas-letkf: Error reading namelist')
330 if ( mmpi_myid == 0 ) write(*,nml=namletkf)
331 ierr = fclos(nulnam)
332
333 !- 1.3 Some minor modifications of namelist values
334 if (hLocalize(1) > 0.0D0 .and. hLocalize(2) < 0.0D0) then
335 ! if only 1 value given for hLocalize, use it for entire column
336 hLocalize(2:4) = hLocalize(1)
337 if ( mmpi_myid == 0 ) write(*,*) 'midas-letkf: hLocalize(2:4) are modified after reading namelist. ' // &
338 'hLocalize(2:4)=', hLocalize(1)
339 else if ( hLocalize(1) < 0.0D0 ) then
340 call utl_abort('midas-letkf: hLocalize(1) < 0.0D0')
341 end if
342 hLocalize(:) = hLocalize(:) * 1000.0D0 ! convert from km to m
343 hLocalizePressure(:) = log(hLocalizePressure(:) * MPC_PA_PER_MBAR_R8)
344
345 if (minDistanceToLand > 0.0D0) then
346 minDistanceToLand = minDistanceToLand * 1000.0D0 ! convert from km to m
347 end if
348
349 if ( trim(algorithm) /= 'LETKF' .and. &
350 trim(algorithm) /= 'CVLETKF' .and. &
351 trim(algorithm) /= 'CVLETKF-PERTOBS' .and. &
352 trim(algorithm) /= 'LETKF-Gain' .and. &
353 trim(algorithm) /= 'LETKF-Gain-ME' .and. &
354 trim(algorithm) /= 'CVLETKF-ME' ) then
355 call utl_abort('midas-letkf: unknown LETKF algorithm: ' // trim(algorithm))
356 end if
357
358 if ( numRetainedEigen < 0 ) call utl_abort('midas-letkf: numRetainedEigen should be ' // &
359 'equal or greater than zero')
360
361 useModulatedEns = (numRetainedEigen > 0)
362
363 if ( trim(algorithm) == 'LETKF-Gain-ME' .or. trim(algorithm) == 'CVLETKF-ME' ) then
364 if ( .not. useModulatedEns ) call utl_abort('midas-letkf: numRetainedEigen should be ' // &
365 'equal or greater than one for LETKF algorithm: ' // &
366 trim(algorithm))
367 else
368 if ( useModulatedEns ) call utl_abort('midas-letkf: numRetainedEigen should be ' // &
369 'equal to zero for LETKF algorithm: ' // &
370 trim(algorithm))
371 end if
372
373 ! check for NO varying horizontal localization lengthscale in letkf with modulated ensembles.
374 if ( .not. all(hLocalize(2:4) == hLocalize(1)) .and. useModulatedEns ) then
375 call utl_abort('midas-letkf: Varying horizontal localization lengthscales is NOT allowed in ' // &
376 'letkf with modulated ensembles')
377 end if
378
379 !
380 !- 2. Initialization
381 !
382
383 !- 2.1 Setup the observation file names and get dateStamp from obs
384 call obsf_setup(dateStampFromObs, midasMode)
385
386 !- 2.2 Initialize date/time-related info
387
388 ! Setup timeCoord module, set dateStamp from env variable
389 call tim_setup()
390 if (tim_getDateStamp() == 0) then
391 if (dateStampFromObs > 0) then
392 ! use dateStamp from obs if not already set by env variable
393 call tim_setDateStamp(dateStampFromObs)
394 else
395 call utl_abort('midas-letkf: DateStamp was not set')
396 end if
397 end if
398 if (tim_nstepobsinc /= 1 .and. tim_nstepobsinc /= tim_nstepobs) then
399 call utl_abort('midas-letkf: invalid value for namelist variable DSTEPOBSINC. ' // &
400 'Increments can be either 3D or have same number of time steps as trials')
401 end if
402 allocate(dateStampList(tim_nstepobs))
403 call tim_getstamplist(dateStampList,tim_nstepobs,tim_getDatestamp())
404 allocate(dateStampListInc(tim_nstepobsinc))
405 call tim_getstamplist(dateStampListInc,tim_nstepobsinc,tim_getDatestamp())
406
407 write(*,*) 'midas-letkf: analysis dateStamp = ',tim_getDatestamp()
408
409 !- 2.3 Initialize variables of the model states
410 call gsv_setup
411
412 !- 2.4 Initialize the Ensemble grid
413 if (mmpi_myid == 0) write(*,*) ''
414 if (mmpi_myid == 0) write(*,*) 'midas-letkf: Set hco and vco parameters for ensemble grid'
415 call fln_ensFileName( ensFileName, ensPathName, memberIndex_opt=1, &
416 copyToRamDisk_opt=.false. )
417 call hco_SetupFromFile( hco_ens, ensFileName, ' ', 'ENSFILEGRID')
418 call vco_setupFromFile( vco_ens, ensFileName )
419 if (vco_getNumLev(vco_ens, 'MM') /= vco_getNumLev(vco_ens, 'TH')) then
420 call utl_abort('midas-letkf: nLev_M /= nLev_T - currently not supported')
421 end if
422 nwpFields = (vco_getNumLev(vco_ens,'TH') > 0 .or. vco_getNumLev(vco_ens,'MM') > 0)
423 oceanFields = (vco_getNumLev(vco_ens,'DP') > 0)
424 if (.not.nwpFields .and. .not.oceanFields) then
425 call utl_abort('midas-letkf: vertical coordinate does not contain nwp nor ocean fields')
426 end if
427
428 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
429
430 !- 2.5 Read in the observations and other obs-related set up
431
432 ! Read the observations
433 call inn_setupObs( obsSpaceData, hco_ens, obsColumnMode, obsMpiStrategy, midasMode, &
434 obsClean_opt = .false. )
435
436 ! Initialize obs error covariances and set flag using 'util' column of stats_tovs
437 call oer_setObsErrors(obsSpaceData, midasMode, useTovsUtil_opt=.true.) ! IN
438
439 ! Call suprep again to filter out channels according to 'util' column of stats_tovs
440 call filt_suprep(obsSpaceData)
441
442 ! Allocate vectors for storing HX values
443 call eob_allocate(ensObs, nEns, obs_numBody(obsSpaceData), obsSpaceData)
444 if ( outputEnsObs ) allocate(ensObs%Ya_r4(ensObs%numMembers,ensObs%numObs))
445 call eob_zero(ensObs)
446 if ( useModulatedEns ) then
447 nEnsGain = nEns * numRetainedEigen
448 allocate(ensObsGain)
449 call eob_allocate(ensObsGain, nEnsGain, obs_numBody(obsSpaceData), obsSpaceData)
450 call eob_zero(ensObsGain)
451 else
452 ensObsGain => ensObs
453 end if
454
455 ! Set lat, lon, obs values in ensObs
456 call eob_setLatLonObs(ensObs)
457 if ( useModulatedEns ) call eob_setLatLonObs(ensObsGain)
458
459 !- 2.6 Initialize a single columnData object
460 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
461 call col_setup
462 call col_setVco(column, vco_ens)
463 call col_allocate(column, obs_numheader(obsSpaceData), &
464 mpiLocal_opt=.true., setToZero_opt=.true.)
465 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
466
467 !- 2.7 Read the sfc height from ensemble member 1 - only if we are doing NWP
468 if ( nwpFields ) then
469 call gsv_allocate( stateVectorHeightSfc, 1, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
470 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
471 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','TT'/) )
472 call gio_readFromFile( stateVectorHeightSfc, ensFileName, ' ', ' ', &
473 containsFullField_opt=.true., readHeightSfc_opt=.true. )
474 end if
475
476 !- 2.8 Allocate statevector related to ensemble mean
477 call gsv_allocate( stateVectorMeanTrl4D, tim_nstepobs, hco_ens, vco_ens, &
478 dateStamp_opt=tim_getDateStamp(), &
479 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
480 dataKind_opt=4, allocHeightSfc_opt=.true., &
481 allocHeight_opt=.false., allocPressure_opt=.false. )
482 call gsv_zero(stateVectorMeanTrl4D)
483 call gsv_allocate( stateVectorMeanAnl, tim_nstepobsinc, hco_ens, vco_ens, &
484 dateStamp_opt=tim_getDateStamp(), &
485 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
486 dataKind_opt=4, allocHeightSfc_opt=.true., &
487 allocHeight_opt=.false., allocPressure_opt=.false. )
488 call gsv_zero(stateVectorMeanAnl)
489 call gsv_allocate( stateVector4D, tim_nstepobs, hco_ens, vco_ens, &
490 dateStamp_opt=tim_getDateStamp(), &
491 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
492 dataKind_opt=4, allocHeightSfc_opt=.true., &
493 allocHeight_opt=.false., allocPressure_opt=.false. )
494 call gsv_zero(stateVector4D)
495 if ( useModulatedEns ) then
496 ! same as stateVector4D
497 call gsv_allocate( stateVector4Dmod, tim_nstepobs, hco_ens, vco_ens, &
498 dateStamp_opt=tim_getDateStamp(), &
499 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
500 dataKind_opt=4, allocHeightSfc_opt=.true., &
501 allocHeight_opt=.false., allocPressure_opt=.false. )
502 call gsv_zero(stateVector4Dmod)
503 end if
504
505 !- 2.9 Allocate statevector related to an analysis ensemble member
506 call gsv_allocate( stateVectorMemberAnl, tim_nstepobsinc, hco_ens, vco_ens, &
507 dateStamp_opt=tim_getDateStamp(), &
508 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
509 dataKind_opt=4, allocHeightSfc_opt=.true. )
510 call gsv_zero(stateVectorMemberAnl)
511
512 !- 2.10 Allocate statevector for storing state with heights and pressures allocated (for s2c_nl)
513 call gsv_allocate( stateVectorWithZandP4D, tim_nstepobs, hco_ens, vco_ens, &
514 dateStamp_opt=tim_getDateStamp(), &
515 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
516 dataKind_opt=4, allocHeightSfc_opt=.true. )
517 call gsv_zero(stateVectorWithZandP4D)
518
519 !- 2.11 Allocate ensembles, read the Trl ensemble
520 call utl_tmg_start(2,'--ReadEnsemble')
521 call ens_allocate(ensembleTrl4D, nEns, tim_nstepobs, hco_ens, vco_ens, dateStampList)
522 call ens_readEnsemble(ensembleTrl4D, ensPathName, biPeriodic=.false., &
523 ignoreDate_opt=ignoreEnsDate)
524 call utl_tmg_stop(2)
525
526 !- 2.12 If desired, read a deterministic state for recentering the ensemble
527 if (recenterInputEns) then
528 call gsv_allocate( stateVectorRecenter, tim_nstepobs, hco_ens, vco_ens, &
529 dateStamp_opt=tim_getDateStamp(), &
530 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
531 dataKind_opt=4, allocHeightSfc_opt=.false., &
532 allocHeight_opt=.false., allocPressure_opt=.false. )
533 call gsv_zero(stateVectorRecenter)
534 call fln_ensTrlFileName( recenterFileName, './', tim_getDateStamp() )
535 do stepIndex = 1, tim_nstepobs
536 call gio_readFromFile( stateVectorRecenter, recenterFileName, ' ', ' ', &
537 stepIndex_opt=stepIndex, containsFullField_opt=.true., &
538 readHeightSfc_opt=.false. )
539 end do
540 call ens_recenter( ensembleTrl4D, stateVectorRecenter, recenteringCoeffScalar_opt=1.0d0 )
541 call gsv_deallocate( stateVectorRecenter )
542 end if
543
544 !- 2.13 Compute ensemble mean and copy to meanTrl and meanAnl stateVectors
545 call ens_computeMean(ensembleTrl4D)
546 call ens_copyEnsMean(ensembleTrl4D, stateVectorMeanTrl4D)
547 if (tim_nstepobsinc < tim_nstepobs) then
548 call gsv_copy4Dto3D(stateVectorMeanTrl4D, stateVectorMeanAnl)
549 else
550 call gsv_copy(stateVectorMeanTrl4D, stateVectorMeanAnl)
551 end if
552
553 !
554 !- 3. Compute HX values with results in ensObs/ensObsGain
555 !
556
557 !- 3.1 Loop over all members and compute HX for each
558 if ( readEnsObsFromFile ) then
559 call eob_readFromFiles(ensObs, nEns, inputFilenamePrefix='eob_HX', readObsInfo=.true.)
560 if ( useModulatedEns ) then
561 call enkf_setupModulationFactor(stateVectorMeanTrl4D%vco, numRetainedEigen, nEns, vLocalize, &
562 beSilent=.true.)
563
564 ! refresh assimilation flag before reading the files
565 call eob_setAssFlag(ensObsGain)
566 call eob_readFromFiles(ensObsGain, nEnsGain, inputFilenamePrefix='eobGain_HX', readObsInfo=.false.)
567 end if
568 else
569 do memberIndex = 1, nEns
570
571 write(*,*) ''
572 write(*,*) 'midas-letkf: apply nonlinear H to ensemble member ', memberIndex
573 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
574
575 ! copy 1 member to a stateVector
576 call ens_copyMember(ensembleTrl4D, stateVector4D, memberIndex)
577 call gsv_copy(stateVector4D, stateVectorWithZandP4D, allowVarMismatch_opt=.true., &
578 beSilent_opt=.true.)
579
580 ! copy the surface height field
581 if (nwpFields) then
582 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
583 end if
584
585 ! Compute and set Yb in ensObs
586 call s2c_nl( stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
587 timeInterpType=obsTimeInterpType, dealloc_opt=.false., &
588 beSilent_opt=.true. )
589
590 ! Compute Y-H(X) in OBS_OMP
591 call inn_computeInnovation(column, obsSpaceData, beSilent_opt=.true.)
592
593 ! Copy to ensObs: Y-HX for this member
594 call eob_setYb(ensObs, memberIndex)
595
596 ! Compute and set Yb in ensObsGain
597 do eigenVectorIndex = 1, numRetainedEigen
598 if ( mmpi_myid == 0 .and. debug ) then
599 write(*,*) 'midas-letkf: apply nonlinear H to modulated member ', &
600 eigenVectorIndex, '/', numRetainedEigen
601 end if
602
603 ! modulate the member with eigenvectors of vertical localization matrix
604 call enkf_getModulatedState( stateVector4D, stateVectorMeanTrl4D, &
605 vLocalize, numRetainedEigen, nEns, &
606 eigenVectorIndex, stateVector4Dmod, &
607 beSilent=.true. )
608 if ( debug ) then
609 call gsv_getField(stateVector4Dmod,field_Psfc,'P0')
610 write(*,*) 'midas-letkf: max(Psfc)=', maxval(field_Psfc), &
611 ', min(Psfc)=', minval(field_Psfc)
612 end if
613
614 call gsv_copy(stateVector4Dmod, stateVectorWithZandP4D, allowVarMismatch_opt=.true., &
615 beSilent_opt=.true.)
616 if (nwpFields) then
617 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
618 end if
619 call s2c_nl( stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
620 timeInterpType=obsTimeInterpType, dealloc_opt=.false., &
621 beSilent_opt=.true. )
622
623 ! Compute Y-H(X) in OBS_OMP
624 call inn_computeInnovation( column, obsSpaceData, filterObsAndInitOer_opt=.false., &
625 beSilent_opt=.true. )
626
627 ! Copy to ensObsGain: Y-HX for this member
628 memberIndexInEnsObs = (eigenVectorIndex - 1) * nEns + memberIndex
629 call eob_setYb(ensObsGain, memberIndexInEnsObs)
630 end do ! eigenVectorIndex
631
632 end do ! memberIndex
633 end if
634 if ( gsv_isAllocated(stateVector4Dmod) ) call gsv_deallocate(stateVector4Dmod)
635
636 ! write local ensObs to file
637 if (writeLocalEnsObsToFile) then
638 call eob_writeToFiles(ensObs, outputFilenamePrefix='eob_HX', writeObsInfo=.true.)
639 if (useModulatedEns) then
640 call eob_writeToFiles(ensObsGain, outputFilenamePrefix='eobGain_HX', writeObsInfo=.false., &
641 numGroupsToDivideMembers_opt=numRetainedEigen, &
642 maxNumMembersPerGroup_opt=nEns)
643 end if
644 end if
645
646 !- 3.2 Set some additional information in ensObs/ensObsGain and additional quality
647 ! control before finally communicating ensObs/ensObsGain globally
648
649 ! Compute and remove the mean of Yb
650 call eob_calcAndRemoveMeanYb(ensObs)
651 if ( useModulatedEns ) call eob_calcAndRemoveMeanYb(ensObsGain)
652
653 ! Put HPHT in OBS_HPHT, for writing to obs files
654 call eob_setHPHT(ensObs)
655 if ( useModulatedEns ) call eob_setHPHT(ensObsGain)
656
657 ! Compute random observation perturbations
658 if (trim(algorithm) == 'CVLETKF-PERTOBS') then
659 randomSeedObs = 1 + mmpi_myid
660 call eob_calcRandPert(ensObs, randomSeedObs)
661 end if
662
663 ! Apply obs operators to ensemble mean background for several purposes
664 write(*,*) ''
665 write(*,*) 'midas-letkf: apply nonlinear H to ensemble mean background'
666 write(*,*) ''
667 call gsv_copy(stateVectorMeanTrl4D, stateVectorWithZandP4D, allowVarMismatch_opt=.true.)
668 if (nwpFields) then
669 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
670 end if
671 call s2c_nl( stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
672 timeInterpType=obsTimeInterpType, dealloc_opt=.false. )
673 call tvs_allocTransmission(col_getNumLev(column,'TH')) ! this will cause radiative transmission profiles to be stored for use in eob_setVertLocation
674 call inn_computeInnovation(column, obsSpaceData, beSilent_opt=.false.)
675
676 ! Put y-mean(H(X)) in OBS_OMP for writing to obs files (overwrites y-H(mean(X)))
677 call eob_setMeanOMP(ensObs)
678 if ( useModulatedEns ) call eob_setMeanOMP(ensObsGain)
679
680 ! Set y for family of interest to mean(H(x)) if doing simulated observations
681 call eob_setSimObsVal(ensObs)
682
683 ! Set vertical location for all obs for vertical localization (based on ensemble mean pressure and height)
684 if (vLocalize > 0.0d0) then
685 if (nwpFields) then
686 call eob_setTypeVertCoord(ensObs,'logPressure')
687 if ( useModulatedEns ) call eob_setTypeVertCoord(ensObsGain,'logPressure')
688 else if (oceanFields) then
689 call eob_setTypeVertCoord(ensObs,'depth')
690 if ( useModulatedEns ) call eob_setTypeVertCoord(ensObsGain,'depth')
691 end if
692 call eob_setVertLocation(ensObs, column)
693 if ( useModulatedEns ) call eob_setVertLocation(ensObsGain, column)
694 end if
695
696 ! Modify the obs error stddev for AMSUB in the tropics
697 if (modifyAmsubObsError) call enkf_modifyAmsubObsError(obsSpaceData)
698
699 ! Apply a background check (reject limit is set in the routine)
700 if (backgroundCheck) call eob_backgroundCheck(ensObs)
701
702 ! For ice/ocean DA: remove obs that are too close to land
703 if (minDistanceToLand > 0.0D0) then
704 call ens_getMask(ensembleTrl4D,oceanMask)
705 call eob_removeObsNearLand(ensObs, oceanMask, minDistanceToLand)
706 end if
707
708 ! Set values of obs_sigi and obs_sigo before hubernorm modifies obs_oer
709 call eob_setSigiSigo(ensObs)
710
711 ! Apply huber norm quality control procedure (modifies obs_oer)
712 if (huberize) call eob_huberNorm(ensObs)
713
714 !- Reject all IR radiance observation in arctic and antarctic (.i.e |lat|>60. )
715 if (rejectHighLatIR) call enkf_rejectHighLatIR(obsSpaceData)
716
717 ! Reject radiance observations too close to the surface
718 if (rejectRadNearSfc) call eob_rejectRadNearSfc(ensObs)
719
720 ! Compute inverse of obs error variance (done here to use dynamic GPS-RO, GB-GPS based on mean O-P)
721 call eob_setObsErrInv(ensObs)
722 if ( useModulatedEns ) call eob_setObsErrInv(ensObsGain)
723
724 call utl_tmg_start(141,'----Barr')
725 call rpn_comm_barrier('GRID',ierr)
726 call utl_tmg_stop(141)
727
728 ! Clean and globally communicate obs-related data to all mpi tasks
729 call eob_allGather(ensObs,ensObs_mpiglobal)
730 if ( useModulatedEns ) then
731 allocate(ensObsGain_mpiglobal)
732 call eob_allGather(ensObsGain,ensObsGain_mpiglobal)
733 else
734 ensObsGain_mpiglobal => ensObs_mpiglobal
735 end if
736
737 ! Print number of assimilated obs per family to the listing
738 write(*,*) 'oti_timeBinning: After extra filtering done in midas-letkf'
739 call oti_timeBinning(obsSpaceData,tim_nstepobs)
740
741 !
742 !- 4. Final preparations for computing analyses
743 !
744
745 !- 4.1 Copy trial ensemble to nstepobsinc time steps
746 if (tim_nstepobsinc < tim_nstepobs) then
747 allocate(ensembleTrl)
748 call ens_allocate(ensembleTrl, nEns, tim_nstepobsinc, hco_ens, vco_ens, dateStampListInc)
749 call ens_copy4Dto3D(ensembleTrl4D,ensembleTrl)
750 call ens_deallocate(ensembleTrl4D)
751 else
752 ! number of timesteps is the same, so just point to it
753 ensembleTrl => ensembleTrl4D
754 end if
755
756 !- 4.2 Copy trl ensemble to anl ensemble
757 call ens_allocate(ensembleAnl, nEns, tim_nstepobsinc, hco_ens, vco_ens, dateStampListInc)
758 call ens_copy(ensembleTrl,ensembleAnl)
759
760 !- 4.3 Setup for interpolating weights from coarse to full resolution
761 call enkf_setupInterpInfo(wInterpInfo, stateVectorMeanAnl%hco, weightLatLonStep, &
762 stateVectorMeanAnl%myLonBeg, stateVectorMeanAnl%myLonEnd, &
763 stateVectorMeanAnl%myLatBeg, stateVectorMeanAnl%myLatEnd)
764
765 !
766 !- 5. Main calculation of ensemble analyses
767 !
768
769 !- 5.1 Call to perform LETKF
770 call enkf_LETKFanalyses(algorithm, numSubEns, randomShuffleSubEns, &
771 ensembleAnl, ensembleTrl, &
772 ensObs_mpiglobal, ensObsGain_mpiglobal, &
773 stateVectorMeanAnl, &
774 wInterpInfo, maxNumLocalObs, &
775 hLocalize, hLocalizePressure, vLocalize, &
776 mpiDistribution, numRetainedEigen)
777
778 !- 5.2 Loop over all analysis members and compute H(Xa_member) (if output is desired)
779 if ( outputEnsObs ) then
780
781 do memberIndex = 1, nEns
782
783 write(*,*) ''
784 write(*,*) 'midas-letkf: apply nonlinear H to analysis ensemble member ', memberIndex
785 write(*,*) ''
786
787 ! copy 1 member to a stateVector
788 call ens_copyMember(ensembleAnl, stateVectorMemberAnl, memberIndex)
789
790 if (tim_nstepobsinc < tim_nstepobs) then
791 ! ensembleAnl is only 3D, so need to make 4D for s2c_nl
792 middleStepIndex = (tim_nstepobs + 1) / 2
793 call gsv_copy(stateVectorMemberAnl, stateVectorWithZandP4D, allowVarMismatch_opt=.true., stepIndexOut_opt=middleStepIndex)
794 call gsv_3dto4d(stateVectorWithZandP4D)
795 else
796 call gsv_copy(stateVectorMemberAnl, stateVectorWithZandP4D, allowVarMismatch_opt=.true.)
797 end if
798
799 ! copy the surface height field
800 if (nwpFields) then
801 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
802 end if
803
804 call s2c_nl( stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
805 timeInterpType=obsTimeInterpType, dealloc_opt=.false. )
806
807 ! Compute Y-H(Xa) in OBS_OMAM (used instead of OBS_OMA so that obsSpaceData isn't unintentionally modified)
808 call inn_computeInnovation(column, obsSpaceData, destObsColumn_opt=OBS_OMAM, beSilent_opt=.true., &
809 callFiltTopo_opt=.false., callSetErrGpsgb_opt=.false., analysisMode_opt=.false.)
810
811 ! Copy to ensObs: Y-H(Xa_member) for this member
812 call eob_setYa(ensObs, memberIndex, OBS_OMAM)
813
814 end do
815
816 end if
817
818 !- 6. Output obs files with mean OMP and (unrecentered) OMA
819
820 ! Compute Y-H(Xa_mean) in OBS_OMA
821 write(*,*) ''
822 write(*,*) 'midas-letkf: apply nonlinear H to ensemble mean analysis'
823 write(*,*) ''
824 if (tim_nstepobsinc < tim_nstepobs) then
825 ! meanAnl is only 3D, so need to make 4D for s2c_nl
826 middleStepIndex = (tim_nstepobs + 1) / 2
827 call gsv_copy(stateVectorMeanAnl, stateVectorWithZandP4D, allowVarMismatch_opt=.true., stepIndexOut_opt=middleStepIndex)
828 call gsv_3dto4d(stateVectorWithZandP4D)
829 else
830 call gsv_copy(stateVectorMeanAnl, stateVectorWithZandP4D, allowVarMismatch_opt=.true.)
831 end if
832 if (nwpFields) then
833 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
834 end if
835 call s2c_nl( stateVectorWithZandP4D, obsSpaceData, column, hco_ens, &
836 timeInterpType=obsTimeInterpType )
837 call inn_computeInnovation(column, obsSpaceData, destObsColumn_opt=OBS_OMA, beSilent_opt=.false.)
838
839 ! Write (update) observation files.
840 if (outputEnsObs) then
841 call obsf_writeFiles( obsSpaceData, ensObs_opt=ensObs)
842 else
843 call obsf_writeFiles( obsSpaceData)
844 end if
845
846 !- 7. Post processing of the analysis results (if desired) and write everything to files
847 if (ensPostProcessing) then
848 !- Allocate and read the Trl control member (used to compute control member increment for IAU)
849 call gsv_allocate( stateVectorCtrlTrl, tim_nstepobsinc, hco_ens, vco_ens, &
850 dateStamp_opt=tim_getDateStamp(), &
851 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
852 dataKind_opt=4, allocHeightSfc_opt=.true., &
853 allocHeight_opt=.false., allocPressure_opt=.false. )
854 if (recenterInputEns) then
855 !- Use the deterministic trial, if we are recentering the input ensemble
856 ctrlFileName = trim(recenterFileName)
857 else
858 !- Otherwise, use member 0000
859 call fln_ensFileName(ctrlFileName, ensPathName, memberIndex_opt=0, &
860 copyToRamDisk_opt=.false.)
861 end if
862 do stepIndex = 1, tim_nstepobsinc
863 call gio_readFromFile( stateVectorCtrlTrl, ctrlFileName, ' ', ' ', &
864 stepIndex_opt=stepIndex, containsFullField_opt=.true., &
865 readHeightSfc_opt=.false. )
866 end do
867
868 call epp_postProcess(ensembleTrl, ensembleAnl, stateVectorHeightSfc, stateVectorCtrlTrl, &
869 writeTrlEnsemble=.false., outputOnlyEnsMean_opt=outputOnlyEnsMean)
870 else
871 !- Just write the raw analysis ensemble to files
872 if (mmpi_myid == 0) then
873 write(*,*) 'midas-letkf: No ensemble post-processing requested, so just write the raw analysis ensemble'
874 end if
875 call utl_tmg_start(3,'--WriteEnsemble')
876 if (.not. outputOnlyEnsMean) then
877 call ens_writeEnsemble(ensembleAnl, '.', '', etiket_anl, 'A', &
878 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
879 containsFullField_opt=.true.)
880 end if
881 call utl_tmg_stop(3)
882
883 end if
884
885 !
886 !- 8. MPI, tmg finalize
887 !
888 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
889 call utl_tmg_stop(0)
890
891 call tmg_terminate(mmpi_myid, 'TMG_INFO')
892 call rpn_comm_finalize(ierr)
893
894 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
895
896 if ( mmpi_myid == 0 ) write(*,*) ' --------------------------------'
897 if ( mmpi_myid == 0 ) write(*,*) ' MIDAS-LETKF ENDS'
898 if ( mmpi_myid == 0 ) write(*,*) ' --------------------------------'
899
900end program midas_letkf
