1program midas_ensPostProcess
2 !
3 !:Purpose: Post-processing program for the local ensemble transform Kalman filter (LETKF).
4 ! Many aspects of this program are controlled throught the namelist block
5 ! namEnsPostProc defined in epp_postProcess.
6 !
7 ! ---
8 !
9 !:Algorithm: The ``ensPostProcess`` program performs several post-processing tasks of
10 ! LETK ensemble analyses and the ensemble trials. The following tasks are
11 ! performed based on the namelist options.
12 !
13 ! - **Mean and std**: Ensemble mean and standard deviation are computed
14 ! upon reading the ensemble members and after the processes which may
15 ! alter the mean or the standard deviation.
16 !
17 ! - **RTPP**: Relaxation to the prior perturbation will relax the analysis
18 ! perturbation to the trial perturbation with a given amount (:math:`\alpha`) in the namelist.
19 ! It can be formulated as:
20 ! :math:`x^a_i \leftarrow \overline{x^a} + (1-\alpha)(\overline{x^a}-x^a_i)+\alpha(\overline{x^b}-x^b_i)`
21 !
22 ! More details for the RTPP can be found in paper:
23 ! `Impacts of Initial Estimate and Observation Availability on Convective-Scale
24 ! Data Assimilation with an Ensemble Kalman Filter
25 ! <https://journals.ametsoc.org/view/journals/mwre/132/5/1520-0493_2004_132_1238_ioieao_2.0.co_2.xml>`_
26 !
27 ! - **RTPS**: Relaxation to the prior spread will relax the analysis ensemble spread (i.e. standard
28 ! deviation) to the trial ensemble spread with a given amount (:math:`\alpha`) in the namelist.
29 ! With the ensemble spread of analysis (:math:`\sigma^a`)
30 ! and trial (:math:`\sigma^b`), it can be formulated as:
31 ! :math:`x^a_i \leftarrow \overline{x^a} + (\overline{x^a}-x^a_i)(1+\alpha(\sigma^b-\sigma^a)/\sigma^a)`
32 !
33 ! More details for the RTPS can be found in paper:
34 ! `Evaluating Methods to Account for System Errors in Ensemble Data Assimilation
35 ! <https://journals.ametsoc.org/view/journals/mwre/140/9/mwr-d-11-00276.1.xml>`_
36 !
37 ! - **Humidity limits**: Impose the saturation limits and RTTOV limits
38 ! on the humidity for each member.
39 !
40 ! - **Recenter**: The ensemble analyses from the LETKF are recentered
41 ! around the EnVar analysis based on the input coefficent.
42 ! It is also called hybrid gain algorithm and more details can be found in
43 ! paper:
44 ! `Using the hybrid gain algorithm to sample data assimilation uncertainty
45 ! <https://rmets.onlinelibrary.wiley.com/doi/10.1002/qj.3426>`_
46 !
47 ! - **Subsample**: Select 20 ensemble members for the medium range forecasts.
48 !
49 ! - **Random perturbation**: Generate homogeneous and isotropic random perturbations
50 ! and add them to the ensemble analyses.
51 !
52 ! - **Analysis increments**: After all the processes above, the analysis incrememnts are
53 ! computed at the central time for all the members including the control member and
54 ! the subsampled 20 members.
55 !
56 ! - **Mask analysis increments**: For LAM grid, the analysis increments in the blending
57 ! area are masked out.
58 !
59 ! --
60 !
61 !============================================= ==============================================================
62 ! Input and Output Files Description of file
63 !============================================= ==============================================================
64 ! ``flnml`` In - Main namelist file with parameters user may modify
65 ! ``ensemble_trial/$YYYYMMDDHH_006_$NNNN`` In - Background ensemble member files
66 ! ``ensemble_anal/$YYYYMMDDHH_000_$NNNN`` In - Analysis ensemble member files
67 ! ``bgcov`` In - Bnmc matrix for the random perturbations
68 ! ``analysis_grid`` In - Horizontal grid file on which the random perturbations
69 ! are generated
70 ! ``$YYYYMMDDHH_recentering_analysis`` In - Analysis from EnVar for recentering the ensemble analyses
71 ! ``analinc_mask`` In - Mask file for masking the analysis increments for LAM grid
72 ! ``$YYYYMMDDHH_000_inc_$NNNN`` Out - Ensemble analysis increments for trials
73 ! ``$YYYYMMDDHH_000_$NNNN`` Out - Ensemble analyses recentered and perturbed
74 ! ``$YYYYMMDDHH_000_$NNNN_raw`` Out - Ensemble analyses from LETKF
75 ! ``$YYYYMMDDHH_000_analmean`` Out - Mean ensemble analysis without perturbation
76 ! ``$YYYYMMDDHH_000_analrms`` Out - Analysis ensemble spread (i.e. standard deviation) without perturbation
77 ! ``$YYYYMMDDHH_000_analrms_ascii`` Out - Global average of ``$YYYYMMDDHH_000_analrms``
78 ! ``$YYYYMMDDHH_000_analpertmean`` Out - Mean ensemble analysis with perturbation and recentering
79 ! ``$YYYYMMDDHH_000_analpertrms`` Out - Analysis ensemble spread with perturbation and recentering
80 ! ``$YYYYMMDDHH_000_analpertrms_ascii`` Out - Global average of ``$YYYYMMDDHH_000_analpertrms``
81 ! ``$YYYYMMDDHH_006_trialmean`` Out - Mean ensemble trial for all time levels
82 ! ``$YYYYMMDDHH_006_trialrms`` Out - Trial ensemble spread at 6H
83 ! ``$YYYYMMDDHH_006_trialrms_ascii`` Out - Global average of ``$YYYYMMDDHH_006_trialrms``
84 ! ``subspace/$YYYYMMDDHH_000_inc_$NNNN`` Out - Subsampled ensemble analysis increments for progs
85 ! ``subspace/$YYYYMMDDHH_000_$NNNN`` Out - Subsampled ensemble analyses recentered and perturbed
86 ! ``subspace_unpert/$YYYYMMDDHH_000_$NNNN`` Out - Unperturbed subsampled ensemble analyses recentered
87 !============================================= ==============================================================
88 !
89 ! --
90 !
91 !:Synopsis: Below is a summary of the ``ensPostProcess`` program calling sequence:
92 !
93 ! - **Initial setups:**
94 !
95 ! - Read the NAMLETKF namelist and check/modify some values.
96 !
97 ! - Various modules are setup: ``gridStateVector_mod``, ``timeCoord_mod`` (and set up dates
98 ! and ``dateStampList`` variables for both trials and increments/analyses).
99 !
100 ! - Setup horizontal and vertical grid objects from first
101 ! ensemble member file.
102 !
103 ! - Allocate varaibles and read ensemble analyses and trials
104 !
105 ! - **Ensemble postprocess:**
106 !
107 ! - Compute ensemble mean and spread (i.e. standard deviation) for analysis and trial
108 !
109 ! - Perform RTPP and RTPS and recompute analysis spread
110 !
111 ! - Recenter ensemble mean and recompute analysis mean and spread
112 !
113 ! - Apply humidity saturation limit and RTTOV HU limit and recompute analysis mean
114 !
115 ! - Subsample ensemble analysis
116 !
117 ! - Add random perturbation to the analyses and compute perturbed analysis spread
118 !
119 ! - Add random perturbation to the subsampled analyses and recenter subsampled mean to the
120 ! global mean.
121 !
122 ! - Compute analysis increments for all analyses and subsampled analyses including the ensemble mean
123 !
124 ! - Mask LAM analysis increments are recompute the ensemble analyses by adding the increments
125 ! to the trials
126 !
127 ! - Writing the outputs.
128 !
129 !:Options: `List of namelist blocks <../namelists_in_each_program.html#ensPostProcess>`_
130 ! that can affect the ``ensPostProcess`` program.
131 !
132 !========================== ========================== ==============================================================
133 ! Program/Module Namelist Description of what is controlled
134 !========================== ========================== ==============================================================
135 ! ``midas_ensPostProcess`` ``NAMENSPOSTPROC`` Number of ensemble members and read, write control of trials
136 ! and analyses and horizontal interpolation degree
137 ! ``enspostprocess_mod`` ``NAMENSPOSTPROCMODULE`` Control the postprocess algorithms
138 ! ``timeCoord_mod`` ``NAMTIME`` Temporal resolution of the background state and the analysis
139 !========================== ========================== ==============================================================
140 !
141
142 use version_mod
143 use midasMpi_mod
144 use fileNames_mod
145 use ensembleStateVector_mod
146 use gridStateVector_mod
147 use gridStateVectorFileIO_mod
148 use verticalCoord_mod
149 use horizontalCoord_mod
150 use timeCoord_mod
151 use utilities_mod
152 use ramDisk_mod
153 use ensPostProcess_mod
154 implicit none
155
156 type(struct_ens) :: ensembleTrl
157 type(struct_ens) :: ensembleAnl
158 type(struct_vco), pointer :: vco_ens => null()
159 type(struct_hco), pointer :: hco_ens => null()
160 type(struct_gsv) :: stateVectorHeightSfc, stateVectorCtrlTrl
161
162 character(len=256) :: ensPathNameAnl = 'ensemble_anal'
163 character(len=256) :: ensPathNameTrl = 'ensemble_trial'
164 character(len=256) :: ensFileName, gridFileName, ctrlFileName
165 integer, allocatable :: dateStampList(:)
166 integer :: ierr, nulnam, stepIndex
167 logical :: targetGridFileExists
168 integer, external :: get_max_rss, fstopc, fnom, fclos
169
170 integer :: nEns ! ensemble size
171 logical :: readTrlEnsemble ! activate reading of trial ensemble
172 logical :: readAnlEnsemble ! activate reading of analysis ensemble
173 logical :: writeTrlEnsemble ! activate writing of the trial ensemble (useful when it's interpolated)
174 character(len=12) :: hInterpolationDegree ! select degree of horizontal interpolation (if needed)
175 NAMELIST /namEnsPostProc/nEns, readTrlEnsemble, readAnlEnsemble, &
176 writeTrlEnsemble, hInterpolationDegree
177
178 call ver_printNameAndVersion('ensPostProcess','Program for post-processing of LETKF analysis ensemble')
179
180 !
181 !- 0. MPI, TMG and misc. initialization
182 !
183 call mmpi_initialize
184 call tmg_init(mmpi_myid, 'TMG_INFO')
185
186 call utl_tmg_start(0,'Main')
187 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
188
189 ! Avoid printing lots of stuff to listing for std file I/O
190 ierr = fstopc('MSGLVL','ERRORS',0)
191
192 ! Setup the ramdisk directory (if supplied)
193 call ram_setup
194
195 !- Setting default namelist variable values
196 nEns = 256
197 readTrlEnsemble = .true.
198 readAnlEnsemble = .true.
199 writeTrlEnsemble = .false.
200 hInterpolationDegree = 'LINEAR' ! or 'CUBIC' or 'NEAREST'
201
202 !- Read the namelist
203 nulnam = 0
204 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
205 read(nulnam, nml=namEnsPostProc, iostat=ierr)
206 if ( ierr /= 0) call utl_abort('midas-ensPostProcess: Error reading namelist')
207 if ( mmpi_myid == 0 ) write(*,nml=namEnsPostProc)
208 ierr = fclos(nulnam)
209
210 if (.not.readTrlEnsemble .and. .not.readAnlEnsemble) then
211 call utl_abort('midas-ensPostProcess: must read either Trial or Analysis ensemble')
212 end if
213
214 if (writeTrlEnsemble .and. .not.readTrlEnsemble) then
215 call utl_abort('midas-ensPostProcess: cannot write Trial ensemble if it is not read')
216 end if
217
218 if (writeTrlEnsemble .and. readAnlEnsemble) then
219 call utl_abort('midas-ensPostProcess: cannot write Trial ensemble when Analysis ensemble is read')
220 end if
221
222 !- 1. Initialize date/time-related info
223
224 ! Setup timeCoord module, set dateStamp with value from trial or analysis ensemble member 1
225 if (readTrlEnsemble) then
226 call fln_ensFileName(ensFileName, ensPathNameTrl, memberIndex_opt=1, &
227 copyToRamDisk_opt=.false.)
228 else
229 call fln_ensFileName(ensFileName, ensPathNameAnl, memberIndex_opt=1, &
230 copyToRamDisk_opt=.false.)
231 end if
232 call tim_setup(fileNameForDate_opt=ensFileName)
233 allocate(dateStampList(tim_nstepobsinc))
234 call tim_getstamplist(dateStampList,tim_nstepobsinc,tim_getDatestamp())
235
236 write(*,*) 'midas-ensPostProcess: analysis dateStamp = ', tim_getDatestamp()
237
238 !- 2. Initialize variables and grids
239
240 !- Initialize variables of the model states
241 call gsv_setup
242
243 !- Initialize the grid from targetgrid file or from trial or analysis ensemble member 1
244 if (mmpi_myid == 0) write(*,*) ''
245 if (mmpi_myid == 0) write(*,*) 'midas-ensPostProcess: Set hco and vco parameters for ensemble grid'
246 inquire(file='targetgrid', exist=targetGridFileExists)
247 if (targetGridFileExists) then
248 gridFileName = 'targetgrid'
249 else if (readTrlEnsemble) then
250 call fln_ensFileName(gridFileName, ensPathNameTrl, memberIndex_opt=1, &
251 copyToRamDisk_opt=.false.)
252 else
253 call fln_ensFileName(gridFileName, ensPathNameAnl, memberIndex_opt=1, &
254 copyToRamDisk_opt=.false.)
255 end if
256 if (mmpi_myid == 0) then
257 write(*,*) 'midas-ensPostProcess: file use to define grid = ', trim(gridFileName)
258 end if
259 call hco_SetupFromFile(hco_ens, gridFileName, ' ', 'ENSFILEGRID')
260 call vco_setupFromFile(vco_ens, gridFileName)
261
262 !- Read the sfc height from trial ensemble member 1 - only if we are doing NWP!
263 if (vco_getNumLev(vco_ens,'TH') > 0 .or. vco_getNumLev(vco_ens,'MM') > 0) then
264 if (readTrlEnsemble) then
265 call fln_ensFileName(ensFileName, ensPathNameTrl, memberIndex_opt=1, &
266 copyToRamDisk_opt=.false.)
267 else
268 call fln_ensFileName(ensFileName, ensPathNameAnl, memberIndex_opt=1, &
269 copyToRamDisk_opt=.false.)
270 end if
271 call gsv_allocate(stateVectorHeightSfc, 1, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
272 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
273 hInterpolateDegree_opt=hInterpolationDegree, &
274 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','TT'/))
275 call gio_readFromFile(stateVectorHeightSfc, ensFileName, ' ', ' ', &
276 containsFullField_opt=.true., readHeightSfc_opt=.true.)
277 end if
278
279 !- 3. Allocate and read ensembles
280
281 call utl_tmg_start(2,'--ReadEnsemble')
282
283 !- Allocate ensembles, read the Anl ensemble
284 if (readAnlEnsemble) then
285 call fln_ensFileName(ensFileName, ensPathNameAnl, resetFileInfo_opt=.true.)
286 call ens_allocate(ensembleAnl, nEns, tim_nstepobsinc, hco_ens, vco_ens, &
287 dateStampList, hInterpolateDegree_opt=hInterpolationDegree)
288 call ens_readEnsemble(ensembleAnl, ensPathNameAnl, biPeriodic=.false.)
289 end if
290
291 !- Allocate ensembles, read the Trl ensemble
292 if (readTrlEnsemble) then
293 call fln_ensFileName(ensFileName, ensPathNameAnl, resetFileInfo_opt=.true.)
294 call ens_allocate(ensembleTrl, nEns, tim_nstepobsinc, hco_ens, vco_ens, &
295 dateStampList, hInterpolateDegree_opt=hInterpolationDegree)
296 call ens_readEnsemble(ensembleTrl, ensPathNameTrl, biPeriodic=.false.)
297
298 end if
299
300 call utl_tmg_stop(2)
301
302 !- Allocate and read the Trl control member
303 if (readTrlEnsemble .and. readAnlEnsemble) then
304 !- Allocate and read control member Trl
305 call gsv_allocate( stateVectorCtrlTrl, tim_nstepobsinc, hco_ens, vco_ens, &
306 dateStamp_opt=tim_getDateStamp(), &
307 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
308 dataKind_opt=4, allocHeightSfc_opt=.true., &
309 hInterpolateDegree_opt=hInterpolationDegree, &
310 allocHeight_opt=.false., allocPressure_opt=.false. )
311 call fln_ensFileName(ctrlFileName, ensPathNameTrl, memberIndex_opt=0, &
312 copyToRamDisk_opt=.false.)
313 do stepIndex = 1, tim_nstepobsinc
314 call gio_readFromFile( stateVectorCtrlTrl, ctrlFileName, ' ', ' ', &
315 stepIndex_opt=stepIndex, containsFullField_opt=.true., &
316 readHeightSfc_opt=.false. )
317 end do
318 end if
319
320 !- 4. Post processing of the analysis results (if desired) and write everything to files
321 call epp_postProcess(ensembleTrl, ensembleAnl, &
322 stateVectorHeightSfc, stateVectorCtrlTrl, &
323 writeTrlEnsemble)
324
325 !
326 !- 5. MPI, tmg finalize
327 !
328 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
329 call utl_tmg_stop(0)
330
331 call tmg_terminate(mmpi_myid, 'TMG_INFO')
332 call rpn_comm_finalize(ierr)
333
334 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
335
336 if ( mmpi_myid == 0 ) write(*,*) ' --------------------------------'
337 if ( mmpi_myid == 0 ) write(*,*) ' MIDAS-ensPostProcess ENDS'
338 if ( mmpi_myid == 0 ) write(*,*) ' --------------------------------'
339
340end program midas_ensPostProcess
