1module ensPostProcess_mod
2 ! MODULE ensPostProcess_mod (prefix='epp' category='1. High-level functionality')
3 !
4 !:Purpose: Various routines that are used to modify or process
5 ! ensembles, usually produced by the LETKF.
6 !
7 use midasMpi_mod
8 use utilities_mod
9 use mathPhysConstants_mod
10 use timeCoord_mod
11 use verticalCoord_mod
12 use horizontalCoord_mod
13 use ensembleStateVector_mod
14 use gridStateVector_mod
15 use gridStateVectorFileIO_mod
16 use interpolation_mod
17 use randomNumber_mod
18 use controlVector_mod
19 use gridVariableTransforms_mod
20 use bMatrix_mod
21 use humidityLimits_mod
22 use varNameList_mod
23 use fileNames_mod
24 use clibInterfaces_mod
25 use calcHeightAndPressure_mod
26 implicit none
27 save
28 private
29
30 ! public procedures
31 public :: epp_postProcess
32
33contains
34
35 !----------------------------------------------------------------------
36 ! epp_postProcess
37 !----------------------------------------------------------------------
38 subroutine epp_postProcess(ensembleTrl, ensembleAnl, &
39 stateVectorHeightSfc, stateVectorCtrlTrl, &
40 writeTrlEnsemble, outputOnlyEnsMean_opt)
41 !
42 !:Purpose: Perform numerous post-processing steps to the ensemble
43 ! produced by the LETKF algorithm.
44 !
45 implicit none
46
47 ! Arguments:
48 type(struct_ens), intent(inout) :: ensembleTrl
49 type(struct_ens), intent(inout) :: ensembleAnl
50 type(struct_gsv), intent(in) :: stateVectorHeightSfc
51 type(struct_gsv), intent(inout) :: stateVectorCtrlTrl
52 logical, intent(in) :: writeTrlEnsemble
53 logical, optional, intent(in) :: outputOnlyEnsMean_opt
54
55 ! Locals:
56 integer :: ierr, nEns, dateStamp, datePrint, timePrint, imode, randomSeedRandomPert
57 integer :: stepIndex, middleStepIndex, nulnam
58 integer, allocatable :: dateStampListInc(:)
59 type(struct_hco), pointer :: hco_ens
60 type(struct_vco), pointer :: vco_ens
61 type(struct_gsv) :: stateVectorMeanAnl, stateVectorMeanAnlRaw, stateVectorMeanTrl
62 type(struct_gsv) :: stateVectorMeanInc
63 type(struct_gsv) :: stateVectorAnalIncMask
64 type(struct_gsv) :: stateVectorStdDevAnl, stateVectorStdDevAnlRaw, stateVectorStdDevAnlPert, stateVectorStdDevTrl
65 type(struct_gsv) :: stateVectorMeanIncSubSample
66 type(struct_gsv) :: stateVectorMeanAnlSubSample
67 type(struct_gsv) :: stateVectorMeanAnlSfcPres
68 type(struct_gsv) :: stateVectorMeanAnlSfcPresMpiGlb
69 type(struct_ens) :: ensembleTrlSubSample
70 type(struct_ens) :: ensembleAnlSubSample
71 type(struct_ens) :: ensembleAnlSubSampleUnPert
72 character(len=12) :: etiket
73 type(struct_ens) :: ensembleAnlInc
74 type(struct_ens) :: ensembleAnlIncSubSample
75 character(len=256) :: outFileName
76 character(len=4), pointer :: varNames(:)
77 character(len=12) :: hInterpolationDegree = 'LINEAR'
78 integer, external :: fnom, fclos, newdate
79 logical :: outputOnlyEnsMean
80
81 ! Namelist variables
82 integer :: randomSeed ! seed used for random perturbation additive inflation
83 logical :: includeYearInSeed ! switch for choosing to include year in default random seed
84 logical :: writeSubSample ! write sub-sample members for initializing medium-range fcsts
85 logical :: writeSubSampleUnPert ! write unperturbed sub-sample members for medium-range fcsts
86 real(8) :: alphaRTPS ! RTPS coefficient (between 0 and 1; 0 means no relaxation)
87 real(8) :: alphaRTPP ! RTPP coefficient (between 0 and 1; 0 means no relaxation)
88 real(8) :: alphaRandomPert ! Random perturbation additive inflation coeff (0->1)
89 real(8) :: alphaRandomPertSubSample ! Random pert. additive inflation coeff for medium-range fcsts
90 logical :: huLimitsBeforeRecenter ! Choose to apply humidity limits before recentering
91 logical :: imposeSaturationLimit ! switch for choosing to impose saturation limit of humidity
92 logical :: imposeRttovHuLimits ! switch for choosing to impose the RTTOV limits on humidity
93 real(8) :: weightRecenter(vco_maxNumLevels) ! weight applied to recentering increment (between 0 and 1; 0 means no recentering)
94 real(8) :: weightRecenterLand ! weight applied to recentering increment for land variables
95 integer :: numMembersToRecenter ! number of members that get recentered on supplied analysis
96 logical :: useOptionTableRecenter ! use values in the optiontable file
97 character(len=8) :: etiket_anl ! etiket for ensemble output files (member number will be appended)
98 character(len=8) :: etiket_inc ! etiket for ensemble output files (member number will be appended)
99 character(len=8) :: etiket_trl ! etiket for ensemble output files (member number will be appended)
100 character(len=12) :: etiket_anlmean ! etiket for mean of analyses and mean of increments files
101 character(len=12) :: etiket_anlrms ! etiket for rms of analyses files
102 character(len=12) :: etiket_anlmean_raw ! etiket for mean of raw analyses
103 character(len=12) :: etiket_anlrms_raw ! etiket for rms of raw analyses
104 character(len=12) :: etiket_anlmeanpert ! etiket for mean of perturbed analyses
105 character(len=12) :: etiket_anlrmspert ! etiket for rms of perturbed analyses
106 character(len=12) :: etiket_trlmean ! etiket for mean of trials
107 character(len=12) :: etiket_trlrms ! etiket for rms of trials
108 integer :: numBits ! number of bits when writing ensemble mean and spread
109 logical :: useAnalIncMask ! mask out the increment on the pilot zone
110 logical :: writeRawAnalStats ! write mean and standard deviation of the raw analysis ensemble
111 logical :: useMemberAsHuRefState ! use each member as reference state for variable transforms
112
113 NAMELIST /namEnsPostProcModule/randomSeed, includeYearInSeed, writeSubSample, writeSubSampleUnPert, &
114 alphaRTPS, alphaRTPP, alphaRandomPert, alphaRandomPertSubSample, &
115 huLimitsBeforeRecenter, imposeSaturationLimit, imposeRttovHuLimits, &
116 weightRecenter, weightRecenterLand, numMembersToRecenter, &
117 useOptionTableRecenter, &
118 etiket_anl, etiket_inc, etiket_trl, etiket_anlmean, etiket_anlrms, &
119 etiket_anlmeanpert, etiket_anlrmspert, &
120 etiket_anlmean_raw, etiket_anlrms_raw, &
121 etiket_trlmean, etiket_trlrms, numBits, useAnalIncMask, &
122 writeRawAnalStats, useMemberAsHuRefState
123
124 if (present(outputOnlyEnsMean_opt)) then
125 outputOnlyEnsMean = outputOnlyEnsMean_opt
126 else
127 outputOnlyEnsMean = .false.
128 end if
129
130 !- Extract the grid definitions and ensemble size
131 if (ens_isAllocated(ensembleTrl)) then
132 hco_ens => ens_getHco(ensembleTrl)
133 vco_ens => ens_getVco(ensembleTrl)
134 nEns = ens_getNumMembers(ensembleTrl)
135 else
136 hco_ens => ens_getHco(ensembleAnl)
137 vco_ens => ens_getVco(ensembleAnl)
138 nEns = ens_getNumMembers(ensembleAnl)
139 end if
140
141 !- Setting default namelist variable values
142 randomSeed = -999
143 includeYearInSeed = .false.
144 writeSubSample = .false.
145 writeSubSampleUnPert = .false.
146 alphaRTPS = 0.0D0
147 alphaRTPP = 0.0D0
148 alphaRandomPert = 0.0D0
149 alphaRandomPertSubSample = -1.0D0
150 huLimitsBeforeRecenter = .true.
151 imposeSaturationLimit = .false.
152 imposeRttovHuLimits = .false.
153 weightRecenter(:) = -1.0D0 ! means no specified values
154 weightRecenterLand = -1.0D0 ! means same recentering for land as other variables
155 numMembersToRecenter = -1 ! means all members recentered by default
156 useOptionTableRecenter = .false.
157 ! For the next 3 variables, the member number will be appended to this string
158 ! for files '${trialdate}_006_${member}'
159 etiket_trl = 'E27_0_0P'
160 ! for files '${analdate}_000_${member}', 'subspace/${analdate}_000_${member}' and 'subspace_unpert/${analdate}_000_${member}'
161 ! the member=0000 will contain the mean analysis
162 etiket_anl = 'E27_0_0P'
163 ! for files '${analdate}_000_inc_${member}' and 'subspace/${analdate}_000_inc_${member}'
164 ! the member=0000 will contain the mean increment
165 etiket_inc = 'E27_0_0P'
166 !
167 etiket_anlmean = 'E27_0_0PAVG' ! for file '${analdate}_000_analmean'
168 etiket_anlrms = 'E27_0_0PRMS' ! for file '${analdate}_000_analrms'
169 etiket_anlmeanpert = 'E27_0_0PAVGP' ! for file '${analdate}_000_analpertmean'
170 etiket_anlrmspert = 'E27_0_0PRMSP' ! for file '${analdate}_000_analpertrms'
171 etiket_anlmean_raw = 'E27_0_0PAVGR' ! for file '${analdate}_000_analmean_raw'
172 etiket_anlrms_raw = 'E27_0_0PRMSR' ! for file '${analdate}_000_analrms_raw'
173 etiket_trlmean = 'E27_0_0PAVG' ! for file '${trialdate}_006_trialmean'
174 etiket_trlrms = 'E27_0_0PRMS' ! for file '${trialdate}_006_trialrms'
175 !
176 numBits = 16
177 useAnalIncMask = .false.
178 writeRawAnalStats = .false.
179 useMemberAsHuRefState = .false.
180
181 !- Read the namelist
182 nulnam = 0
183 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
184 read(nulnam, nml=namEnsPostProcModule, iostat=ierr)
185 if ( ierr /= 0) call utl_abort('epp_postProc: Error reading namelist')
186 if ( mmpi_myid == 0 ) write(*,nml=namEnsPostProcModule)
187 ierr = fclos(nulnam)
188
189 if (alphaRTPS < 0.0D0) alphaRTPS = 0.0D0
190 if (alphaRTPP < 0.0D0) alphaRTPP = 0.0D0
191 if (alphaRandomPert < 0.0D0) alphaRandomPert = 0.0D0
192 if (alphaRandomPertSubSample < 0.0D0) alphaRandomPertSubSample = 0.0D0
193 if (numMembersToRecenter == -1) numMembersToRecenter = nEns ! default behaviour
194
195 if (writeSubSample) then
196 if (.not.(ens_isAllocated(ensembleTrl).and.ens_isAllocated(ensembleAnl))) then
197 call utl_abort('epp_postProc: subSample can only be produced if both Anl and Trl ensembles available')
198 end if
199 end if
200
201 if (writeSubSampleUnPert) then
202 if (.not.ens_isAllocated(ensembleAnl)) then
203 call utl_abort('epp_postProc: subSampleUnPert can only be produced if Anl ensemble available')
204 end if
205 end if
206
207 if (writeRawAnalStats) then
208 if (.not.ens_isAllocated(ensembleAnl)) then
209 call utl_abort('epp_postProc: RawAnalStats can only be produced if Anl ensemble available')
210 end if
211 end if
212
213 if (ens_isAllocated(ensembleTrl)) then
214 !- Allocate and compute ensemble mean Trl
215 call gsv_allocate( stateVectorMeanTrl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
216 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
217 dataKind_opt=4, allocHeightSfc_opt=.true., &
218 hInterpolateDegree_opt = hInterpolationDegree, &
219 allocHeight_opt=.false., allocPressure_opt=.false. )
220 call gsv_zero(stateVectorMeanTrl)
221 call ens_computeMean(ensembleTrl)
222 call ens_copyEnsMean(ensembleTrl, stateVectorMeanTrl)
223
224 !- Allocate and compute ensemble spread stddev Trl
225 call gsv_allocate( stateVectorStdDevTrl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
226 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
227 hInterpolateDegree_opt = hInterpolationDegree, &
228 dataKind_opt=4, allocHeight_opt=.false., allocPressure_opt=.false. )
229 call gsv_zero(stateVectorStdDevTrl)
230 call ens_computeStdDev(ensembleTrl)
231 call ens_copyEnsStdDev(ensembleTrl, stateVectorStdDevTrl)
232 end if
233
234 if (ens_isAllocated(ensembleAnl)) then
235 !- Allocate and compute ensemble mean Anl
236 call gsv_allocate( stateVectorMeanAnl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
237 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
238 dataKind_opt=4, allocHeightSfc_opt=.true., &
239 hInterpolateDegree_opt = hInterpolationDegree, &
240 allocHeight_opt=.false., allocPressure_opt=.false. )
241 call gsv_zero(stateVectorMeanAnl)
242 call ens_computeMean(ensembleAnl)
243 call ens_copyEnsMean(ensembleAnl, stateVectorMeanAnl)
244
245 !- Allocate and compute ensemble spread stddev Anl
246 call gsv_allocate( stateVectorStdDevAnl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
247 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
248 hInterpolateDegree_opt = hInterpolationDegree, &
249 dataKind_opt=4, allocHeight_opt=.false., allocPressure_opt=.false. )
250 call gsv_zero(stateVectorStdDevAnl)
251 call ens_computeStdDev(ensembleAnl)
252 call ens_copyEnsStdDev(ensembleAnl, stateVectorStdDevAnl)
253
254 !- Keep the mean and standard deviation of the raw analysis for outputs, if requested
255 if (writeRawAnalStats) then
256 call gsv_allocate( stateVectorMeanAnlRaw, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
257 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
258 dataKind_opt=4, allocHeightSfc_opt=.true., &
259 hInterpolateDegree_opt = hInterpolationDegree, &
260 allocHeight_opt=.false., allocPressure_opt=.false. )
261 call gsv_allocate( stateVectorStdDevAnlRaw, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
262 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
263 hInterpolateDegree_opt = hInterpolationDegree, &
264 dataKind_opt=4, allocHeight_opt=.false., allocPressure_opt=.false. )
265 call gsv_copy(stateVectorMeanAnl, stateVectorMeanAnlRaw)
266 call gsv_copy(stateVectorStdDevAnl, stateVectorStdDevAnlRaw)
267 end if
268
269 if (ens_isAllocated(ensembleTrl)) then
270 !- Apply RTPP, if requested
271 if (alphaRTPP > 0.0D0) then
272 call epp_RTPP(ensembleAnl, ensembleTrl, stateVectorMeanAnl, &
273 stateVectorMeanTrl, alphaRTPP)
274 ! recompute the analysis spread stddev
275 call ens_computeStdDev(ensembleAnl)
276 call ens_copyEnsStdDev(ensembleAnl, stateVectorStdDevAnl)
277 end if
278
279 !- Apply RTPS, if requested
280 if (alphaRTPS > 0.0D0) then
281 call epp_RTPS(ensembleAnl, stateVectorStdDevAnl, stateVectorStdDevTrl, &
282 stateVectorMeanAnl, alphaRTPS)
283 ! recompute the analysis spread stddev
284 call ens_computeStdDev(ensembleAnl)
285 call ens_copyEnsStdDev(ensembleAnl, stateVectorStdDevAnl)
286 end if
287 end if
288
289 !- Impose limits on humidity *before* recentering, if requested
290 if (huLimitsBeforeRecenter) then
291 if (imposeSaturationLimit .or. imposeRttovHuLimits) then
292 if (mmpi_myid == 0) write(*,*) ''
293 if (mmpi_myid == 0) write(*,*) 'epp_postProcess: limits will be imposed on the humidity of analysis ensemble'
294 if (mmpi_myid == 0 .and. imposeSaturationLimit ) write(*,*) ' -> Saturation Limit'
295 if (mmpi_myid == 0 .and. imposeRttovHuLimits ) write(*,*) ' -> Rttov Limit'
296 if ( imposeSaturationLimit ) call qlim_saturationLimit(ensembleAnl)
297 if ( imposeRttovHuLimits ) call qlim_rttovLimit (ensembleAnl)
298 ! And recompute analysis mean
299 call ens_computeMean(ensembleAnl)
300 call ens_copyEnsMean(ensembleAnl, stateVectorMeanAnl)
301 end if
302 end if
303
304 !- Recenter analysis ensemble on supplied analysis
305 if (any(weightRecenter(:) >= 0.0d0) .or. useOptionTableRecenter) then
306 write(*,*) 'epp_postProcess: Recenter analyses on supplied analysis'
307 call epp_hybridRecentering(ensembleAnl, weightRecenter, weightRecenterLand, &
308 useOptionTableRecenter, numMembersToRecenter)
309 ! And recompute analysis mean
310 call ens_computeMean(ensembleAnl)
311 call ens_copyEnsMean(ensembleAnl, stateVectorMeanAnl)
312 ! And recompute the analysis spread stddev
313 call ens_computeStdDev(ensembleAnl)
314 call ens_copyEnsStdDev(ensembleAnl, stateVectorStdDevAnl)
315 end if
316
317 !- Impose limits on humidity *after* recentering, if requested
318 if (.not.huLimitsBeforeRecenter) then
319 if (imposeSaturationLimit .or. imposeRttovHuLimits) then
320 if (mmpi_myid == 0) write(*,*) ''
321 if (mmpi_myid == 0) write(*,*) 'epp_postProcess: limits will be imposed on the humidity of analysis ensemble'
322 if (mmpi_myid == 0 .and. imposeSaturationLimit ) write(*,*) ' -> Saturation Limit'
323 if (mmpi_myid == 0 .and. imposeRttovHuLimits ) write(*,*) ' -> Rttov Limit'
324 if ( imposeSaturationLimit ) call qlim_saturationLimit(ensembleAnl)
325 if ( imposeRttovHuLimits ) call qlim_rttovLimit (ensembleAnl)
326 ! And recompute analysis mean
327 call ens_computeMean(ensembleAnl)
328 call ens_copyEnsMean(ensembleAnl, stateVectorMeanAnl)
329 end if
330 end if
331
332 !- If SubSample requested, copy sub-sample of analysis and trial members
333 if (writeSubSample) then
334 ! Copy sub-sampled analysis and trial ensemble members
335 call epp_selectSubSample(ensembleAnl, ensembleAnlSubSample, &
336 ensembleTrl, ensembleTrlSubSample)
337
338 ! Create subdirectory for outputting sub sample increments
339 ierr = clib_mkdir_r('subspace')
340
341 ! Allocate stateVectors to store and output sub-sampled ensemble mean analysis
342 call gsv_allocate( stateVectorMeanAnlSubSample, tim_nstepobsinc, &
343 hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
344 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
345 dataKind_opt=4, allocHeightSfc_opt=.true., &
346 hInterpolateDegree_opt = hInterpolationDegree, &
347 allocHeight_opt=.false., allocPressure_opt=.false. )
348 call gsv_zero(stateVectorMeanAnlSubSample)
349
350 end if
351
352 !- If unperturbed SubSample requested, copy sub-sample of analysis members
353 if (writeSubSampleUnPert) then
354 ! Copy sub-sampled analysis ensemble members
355 call epp_selectSubSample(ensembleAnl, ensembleAnlSubSampleUnPert)
356
357 ! Create subdirectory for outputting sub sample members without perturbations
358 ierr = clib_mkdir_r('subspace_unpert')
359
360 end if
361
362 !- Apply random additive inflation, if requested
363 if (alphaRandomPert > 0.0D0) then
364 ! If namelist value is -999, set random seed using the date (as in standard EnKF)
365 if (randomSeed == -999) then
366 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
367 dateStamp = tim_getDateStamp()
368 ierr = newdate(dateStamp, datePrint, timePrint, imode)
369 timePrint = timePrint/1000000
370 datePrint = datePrint*100 + timePrint
371 if (includeYearInSeed) then
372 ! Remove the century, keeping 2 digits of the year
373 randomSeedRandomPert = datePrint - 100000000*(datePrint/100000000)
374 else
375 ! Remove the year and add 9
376 randomSeedRandomPert = 9 + datePrint - 1000000*(datePrint/1000000)
377 end if
378 else
379 randomSeedRandomPert = randomSeed
380 end if
381 write(*,*) 'epp_postProcess: randomSeed for additive inflation set to ', &
382 randomSeedRandomPert
383 if (ens_isAllocated(ensembleTrl)) then
384 call epp_addRandomPert(ensembleAnl, stateVectorMeanTrl, alphaRandomPert, &
385 randomSeedRandomPert, useMemberAsHuRefState)
386 else
387 call epp_addRandomPert(ensembleAnl, stateVectorMeanAnl, alphaRandomPert, &
388 randomSeedRandomPert, useMemberAsHuRefState)
389 end if
390 end if
391
392 !- Recompute the analysis spread stddev after inflation and humidity limits
393 call gsv_allocate( stateVectorStdDevAnlPert, tim_nstepobsinc, hco_ens, vco_ens, &
394 dateStamp_opt=tim_getDateStamp(), &
395 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
396 hInterpolateDegree_opt = hInterpolationDegree, &
397 dataKind_opt=4, allocHeight_opt=.false., &
398 allocPressure_opt=.false. )
399 call ens_computeStdDev(ensembleAnl)
400 call ens_copyEnsStdDev(ensembleAnl, stateVectorStdDevAnlPert)
401
402 !- If SubSample requested, do remaining processing and output of sub-sampled members
403 if (writeSubSample) then
404
405 ! Apply random additive inflation to sub-sampled ensemble, if requested
406 if (alphaRandomPertSubSample > 0.0D0) then
407 ! If namelist value is -999, set random seed using the date (as in standard EnKF)
408 if (randomSeed == -999) then
409 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
410 dateStamp = tim_getDateStamp()
411 ierr = newdate(dateStamp, datePrint, timePrint, imode)
412 timePrint = timePrint/1000000
413 datePrint = datePrint*100 + timePrint
414 if (includeYearInSeed) then
415 ! Remove the century, keeping 2 digits of the year
416 randomSeedRandomPert = datePrint - 100000000*(datePrint/100000000)
417 else
418 ! Remove the year and add 9
419 randomSeedRandomPert = 9 + datePrint - 1000000*(datePrint/1000000)
420 end if
421 else
422 randomSeedRandomPert = randomSeed
423 end if
424 write(*,*) 'epp_postProcess: randomSeed for additive inflation set to ', &
425 randomSeedRandomPert
426 call epp_addRandomPert(ensembleAnlSubSample, stateVectorMeanTrl, &
427 alphaRandomPertSubSample, randomSeedRandomPert, &
428 useMemberAsHuRefState)
429 end if
430
431 ! Compute analysis mean of sub-sampled ensemble
432 call ens_computeMean(ensembleAnlSubSample)
433
434 ! Shift members to have same mean as full ensemble
435 call ens_recenter(ensembleAnlSubSample, stateVectorMeanAnl, &
436 recenteringCoeffScalar_opt=1.0D0)
437
438 ! Re-compute analysis mean of sub-sampled ensemble
439 call ens_computeMean(ensembleAnlSubSample)
440 call ens_copyEnsMean(ensembleAnlSubSample, stateVectorMeanAnlSubSample)
441
442 end if ! writeSubsample
443
444 !- If SubSample requested, do remaining processing and output of sub-sampled members
445 if (writeSubSampleUnPert) then
446
447 ! Compute analysis mean of sub-sampled ensemble
448 call ens_computeMean(ensembleAnlSubSampleUnPert)
449
450 ! Shift members to have same mean as full ensemble
451 call ens_recenter(ensembleAnlSubSampleUnPert, stateVectorMeanAnl, &
452 recenteringCoeffScalar_opt=1.0D0)
453
454 end if
455
456 end if ! ens_isAllocated(ensembleAnl)
457
458 !
459 !- Transform data before computing the increments and writing to files.
460 !
461 if (ens_isAllocated(ensembleAnl)) then
462 call gvt_transform(ensembleAnl,'AllTransformedToModel',allowOverWrite_opt=.true.)
463 call gvt_transform(stateVectorMeanAnl,'AllTransformedToModel',allowOverWrite_opt=.true.)
464 if (writeSubSample) then
465 call gvt_transform(ensembleAnlSubSample,'AllTransformedToModel',allowOverWrite_opt=.true.)
466 call gvt_transform(stateVectorMeanAnlSubSample,'AllTransformedToModel',allowOverWrite_opt=.true.)
467 end if
468 if (writeSubSampleUnPert) then
469 call gvt_transform(ensembleAnlSubSampleUnPert,'AllTransformedToModel',allowOverWrite_opt=.true.)
470 end if
471 end if
472
473 ! When we read ensemble trials we always need to transform them either for incremnets or for writing
474 if (ens_isAllocated(ensembleTrl)) then
475 call gvt_transform(ensembleTrl,'AllTransformedToModel',allowOverWrite_opt=.true.)
476 call gvt_transform(stateVectorCtrlTrl,'AllTransformedToModel',allowOverWrite_opt=.true.)
477 call gvt_transform(stateVectorMeanTrl,'AllTransformedToModel',allowOverWrite_opt=.true.)
478 if (writeSubSample) then
479 call gvt_transform(ensembleTrlSubSample,'AllTransformedToModel',allowOverWrite_opt=.true.)
480 end if
481 end if
482
483 !
484 !- Compute increments
485 !
486 if (ens_isAllocated(ensembleAnl).and.ens_isAllocated(ensembleTrl)) then
487
488 !- Read the analysis mask (in LAM mode only) - N.B. different from land/sea mask!!!
489 if (.not. hco_ens%global .and. useAnalIncMask) then
490 call gio_getMaskLAM(stateVectorAnalIncMask, hco_ens, vco_ens, hInterpolationDegree)
491 end if
492
493 ! initialize the vector for the ensemble increment
494 allocate(dateStampListInc(tim_nstepobsinc))
495 call tim_getstamplist(dateStampListInc,tim_nstepobsinc,tim_getDatestamp())
496 nullify(varNames)
497 call ens_varNamesList(varNames,ensembleTrl)
498 call ens_allocate(ensembleAnlInc, nEns, tim_nstepobsinc, hco_ens, vco_ens, &
499 dateStampListInc, varNames_opt=varNames)
500 call ens_copy(ensembleTrl,ensembleAnlInc)
501 deallocate(varNames)
502
503 ! Compute the ensemble increments
504 call ens_add(ensembleAnl, ensembleAnlInc, scaleFactorInOut_opt=-1.0D0)
505
506 !- Mask the ensemble increments for LAM grid and recompute ensemble analysis
507 if (.not. hco_ens%global .and. useAnalIncMask) then
508 call ens_applyMaskLAM(ensembleAnlInc, stateVectorAnalIncMask)
509 call ens_copy(ensembleAnlInc,ensembleAnl)
510 call ens_add(ensembleTrl,ensembleAnl)
511 end if
512
513 ! Compute mean increment for converted model variables (e.g. VIS and PR)
514 nullify(varNames)
515 call gsv_varNamesList(varNames, stateVectorMeanAnl)
516 call gsv_allocate( stateVectorMeanInc, tim_nstepobsinc, hco_ens, vco_ens, &
517 dateStamp_opt=tim_getDateStamp(), &
518 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
519 hInterpolateDegree_opt = hInterpolationDegree, &
520 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=varNames )
521 call gsv_copy(stateVectorMeanAnl, stateVectorMeanInc)
522 deallocate(varNames)
523
524 call gsv_add(stateVectorCtrlTrl, stateVectorMeanInc, scaleFactor_opt=-1.0D0)
525
526 !- Mask the mean increment for LAM grid and recompute the mean analysis
527 if (.not. hco_ens%global .and. useAnalIncMask) then
528 call gsv_applyMaskLAM(stateVectorMeanInc, stateVectorAnalIncMask)
529 call gsv_copy(stateVectorMeanInc,stateVectorMeanAnl)
530 call gsv_add(stateVectorCtrlTrl,stateVectorMeanAnl)
531 end if
532
533 ! If subsample reqested compute the subsample increments
534 if (writeSubSample) then
535 ! Compute ensemble increments for subsample
536 nullify(varNames)
537 call ens_varNamesList(varNames,ensembleAnlSubSample)
538 call ens_allocate(ensembleAnlIncSubSample, ens_getNumMembers(ensembleAnlSubSample), &
539 tim_nstepobsinc, hco_ens, vco_ens, dateStampListInc, varNames_opt=varNames)
540 call ens_copy(ensembleTrlSubSample,ensembleAnlIncSubSample)
541 deallocate(varNames)
542
543 call ens_add(ensembleAnlSubSample,ensembleAnlIncSubSample,scaleFactorInOut_opt=-1.0D0)
544
545 !- Mask the ensemble increment for LAM grid and recompute the mean analysis
546 if (.not. hco_ens%global .and. useAnalIncMask) then
547 call ens_applyMaskLAM(ensembleAnlIncSubSample, stateVectorAnalIncMask)
548 call ens_copy(ensembleAnlIncSubSample,ensembleAnlSubSample)
549 call ens_add(ensembleTrlSubSample,ensembleAnlSubSample)
550 end if
551
552 ! Compute mean increment with respect to mean of full trial ensemble
553 nullify(varNames)
554 call gsv_varNamesList(varNames, stateVectorMeanAnlSubSample)
555 call gsv_allocate( stateVectorMeanIncSubSample, tim_nstepobsinc, &
556 hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
557 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
558 dataKind_opt=4, allocHeightSfc_opt=.true., &
559 hInterpolateDegree_opt = hInterpolationDegree, &
560! allocHeight_opt=.false., allocPressure_opt=.false., &
561 varNames_opt=varNames )
562 call gsv_copy(stateVectorMeanAnlSubSample, stateVectorMeanIncSubSample)
563 deallocate(varNames)
564
565 call gsv_add(stateVectorCtrlTrl, stateVectorMeanIncSubSample, scaleFactor_opt=-1.0D0)
566
567 !- Mask the mean increment for LAM grid and recompute the mean analysis
568 if (.not. hco_ens%global .and. useAnalIncMask) then
569 call gsv_applyMaskLAM(stateVectorMeanIncSubSample, stateVectorAnalIncMask)
570 call gsv_copy(stateVectorMeanIncSubSample,stateVectorMeanAnlSubSample)
571 call gsv_add(stateVectorCtrlTrl,stateVectorMeanAnlSubSample)
572 end if
573 end if
574
575 end if !- end of increment computation
576
577 !
578 !- Output everything
579 !
580
581 !- Output ens stddev and mean in trialrms, analrms and analpertrms files
582
583 ! determine middle timestep for output of these files
584 middleStepIndex = (tim_nstepobsinc + 1) / 2
585
586 if (ens_isAllocated(ensembleTrl)) then
587 ! output trialmean, trialrms
588 call utl_tmg_start(5,'--WriteEnsMeanRms')
589 call fln_ensTrlFileName(outFileName, '.', tim_getDateStamp())
590 outFileName = trim(outFileName) // '_trialmean'
591 call ens_copyMaskToGsv(ensembleTrl, stateVectorMeanTrl)
592 do stepIndex = 1, tim_nstepobsinc
593 call gio_writeToFile(stateVectorMeanTrl, outFileName, etiket_trlmean, &
594 typvar_opt='P', writeHeightSfc_opt=.false., numBits_opt=numBits, &
595 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
596 end do
597 call fln_ensTrlFileName(outFileName, '.', tim_getDateStamp())
598 outFileName = trim(outFileName) // '_trialrms'
599 call ens_copyMaskToGsv(ensembleTrl, stateVectorStdDevTrl)
600 call gio_writeToFile(stateVectorStdDevTrl, outFileName, etiket_trlrms, &
601 typvar_opt='P', writeHeightSfc_opt=.false., numBits_opt=numBits, &
602 stepIndex_opt=middleStepIndex, containsFullField_opt=.false.)
603 outFileName = trim(outFileName) // '_ascii'
604 call epp_printRmsStats(stateVectorStdDevTrl,outFileName,elapsed=0.0D0,ftype='F',nEns=nEns)
605 call utl_tmg_stop(5)
606
607 ! output the trial ensemble if requested (because it was interpolated)
608 if (writeTrlEnsemble) then
609 call utl_tmg_start(3,'--WriteEnsemble')
610 if (.not. outputOnlyEnsMean) then
611 call ens_writeEnsemble(ensembleTrl, '.', '', etiket_trl, 'P', &
612 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
613 containsFullField_opt=.true.)
614 end if
615 call utl_tmg_stop(3)
616 end if
617 end if
618
619 ! all outputs related to analysis ensemble
620 if (ens_isAllocated(ensembleAnl)) then
621
622 !- Prepare stateVector with only MeanAnl surface pressure and surface height
623 if (gsv_isAllocated(stateVectorHeightSfc)) then
624 call gsv_allocate( stateVectorMeanAnlSfcPres, tim_nstepobsinc, hco_ens, vco_ens, &
625 dateStamp_opt=tim_getDateStamp(), &
626 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
627 hInterpolateDegree_opt = hInterpolationDegree, &
628 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0'/) )
629 call gsv_zero(stateVectorMeanAnlSfcPres)
630 if (mmpi_myid <= (nEns-1)) then
631 call gsv_allocate( stateVectorMeanAnlSfcPresMpiGlb, tim_nstepobsinc, hco_ens, vco_ens, &
632 dateStamp_opt=tim_getDateStamp(), &
633 mpi_local_opt=.false., &
634 hInterpolateDegree_opt = hInterpolationDegree, &
635 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0'/) )
636 call gsv_zero(stateVectorMeanAnlSfcPresMpiGlb)
637 end if
638 call gsv_copy(stateVectorMeanAnl, stateVectorMeanAnlSfcPres, allowVarMismatch_opt=.true.)
639 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorMeanAnlSfcPres)
640 call gsv_transposeTilesToMpiGlobal(stateVectorMeanAnlSfcPresMpiGlb, stateVectorMeanAnlSfcPres)
641 end if
642
643 ! output analmean, analrms
644 call utl_tmg_start(5,'--WriteEnsMeanRms')
645 call fln_ensAnlFileName(outFileName, '.', tim_getDateStamp())
646 outFileName = trim(outFileName) // '_analmean'
647 call ens_copyMaskToGsv(ensembleAnl, stateVectorMeanAnl)
648 do stepIndex = 1, tim_nstepobsinc
649 call gio_writeToFile(stateVectorMeanAnl, outFileName, etiket_anlmean, &
650 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
651 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
652 end do
653 call fln_ensAnlFileName(outFileName, '.', tim_getDateStamp())
654 outFileName = trim(outFileName) // '_analrms'
655 call ens_copyMaskToGsv(ensembleAnl, stateVectorStdDevAnl)
656 call gio_writeToFile(stateVectorStdDevAnl, outFileName, etiket_anlrms, &
657 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
658 stepIndex_opt=middleStepIndex, containsFullField_opt=.false.)
659 outFileName = trim(outFileName) // '_ascii'
660 call epp_printRmsStats(stateVectorStdDevAnl,outFileName,elapsed=0.0D0,ftype='A',nEns=nEns)
661
662 ! output analmean_raw and analrms_raw, if requested
663 if (writeRawAnalStats) then
664 call fln_ensAnlFileName(outFileName, '.', tim_getDateStamp())
665 outFileName = trim(outFileName) // '_analmean_raw'
666 call ens_copyMaskToGsv(ensembleAnl, stateVectorMeanAnlRaw)
667 do stepIndex = 1, tim_nstepobsinc
668 call gio_writeToFile(stateVectorMeanAnlRaw, outFileName, etiket_anlmean_raw, &
669 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
670 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
671 end do
672 call fln_ensAnlFileName(outFileName, '.', tim_getDateStamp())
673 outFileName = trim(outFileName) // '_analrms_raw'
674 call ens_copyMaskToGsv(ensembleAnl, stateVectorStdDevAnl)
675 call gio_writeToFile(stateVectorStdDevAnlRaw, outFileName, etiket_anlrms_raw, &
676 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
677 stepIndex_opt=middleStepIndex, containsFullField_opt=.false.)
678 outFileName = trim(outFileName) // '_ascii'
679 call epp_printRmsStats(stateVectorStdDevAnlRaw,outFileName,elapsed=0.0D0,ftype='A',nEns=nEns)
680 end if ! writeRawAnalStats
681
682 if (alphaRandomPert > 0.0D0) then
683 ! output analpertmean, analpertrms
684 call fln_ensAnlFileName( outFileName, '.', tim_getDateStamp() )
685 outFileName = trim(outFileName) // '_analpertmean'
686 call ens_copyMaskToGsv(ensembleAnl, stateVectorMeanAnl)
687 do stepIndex = 1, tim_nstepobsinc
688 call gio_writeToFile(stateVectorMeanAnl, outFileName, etiket_anlmeanpert, &
689 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
690 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
691 end do
692 call fln_ensAnlFileName( outFileName, '.', tim_getDateStamp() )
693 outFileName = trim(outFileName) // '_analpertrms'
694 call ens_copyMaskToGsv(ensembleAnl, stateVectorStdDevAnlPert)
695 call gio_writeToFile(stateVectorStdDevAnlPert, outFileName, etiket_anlrmspert, &
696 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
697 stepIndex_opt=middleStepIndex, containsFullField_opt=.false.)
698 outFileName = trim(outFileName) // '_ascii'
699 call epp_printRmsStats(stateVectorStdDevAnlPert,outFileName,elapsed=0.0D0,ftype='P',nEns=nEns)
700 end if
701 call utl_tmg_stop(5)
702
703 !- Output the ensemble mean increment (include MeanAnl Psfc) and ensemble increments
704 if (ens_isAllocated(ensembleTrl)) then
705 call utl_tmg_start(5,'--WriteEnsMeanRms')
706 ! output ensemble mean increment
707 call fln_ensAnlFileName( outFileName, '.', tim_getDateStamp(), 0, ensFileNameSuffix_opt='inc' )
708 call ens_copyMaskToGsv(ensembleAnl, stateVectorMeanInc)
709 ! here we assume 4 digits for the ensemble member!!!!
710 etiket = trim(etiket_inc) // '0000'
711 do stepIndex = 1, tim_nstepobsinc
712 call gio_writeToFile(stateVectorMeanInc, outFileName, etiket, &
713 typvar_opt='R', writeHeightSfc_opt=.false., numBits_opt=numBits, &
714 stepIndex_opt=stepIndex, containsFullField_opt=.false.)
715 if (gsv_isAllocated(stateVectorMeanAnlSfcPres)) then
716 call gio_writeToFile(stateVectorMeanAnlSfcPres, outFileName, etiket, &
717 typvar_opt='A', writeHeightSfc_opt=.true., &
718 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
719 end if
720 end do
721 call utl_tmg_stop(5)
722
723 !- Output all ensemble member increments
724 call utl_tmg_start(3,'--WriteEnsemble')
725 if (.not. outputOnlyEnsMean) then
726 call ens_writeEnsemble(ensembleAnlInc, '.', '', etiket_inc, 'R', &
727 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
728 containsFullField_opt=.false., resetTimeParams_opt=.true.)
729 if (gsv_isAllocated(stateVectorMeanAnlSfcPresMpiGlb)) then
730 ! Also write the reference (analysis) surface pressure to increment files
731 call epp_writeToAllMembers(stateVectorMeanAnlSfcPresMpiGlb, nEns, &
732 etiket='ENS_INC', typvar='A', fileNameSuffix='inc', &
733 ensPath='.')
734 end if
735 end if
736 call utl_tmg_stop(3)
737
738 end if ! allocated(ensembleTrl)
739
740 ! output ensemble mean analysis state
741 call utl_tmg_start(5,'--WriteEnsMeanRms')
742 call fln_ensAnlFileName( outFileName, '.', tim_getDateStamp(), 0 )
743 call ens_copyMaskToGsv(ensembleAnl, stateVectorMeanAnl)
744 ! here we assume 4 digits for the ensemble member!!!!
745 etiket = trim(etiket_anl) // '0000'
746 do stepIndex = 1, tim_nstepobsinc
747 call gio_writeToFile(stateVectorMeanAnl, outFileName, etiket, &
748 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
749 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
750 end do
751 call utl_tmg_stop(5)
752
753 !- Output all ensemble member analyses
754 ! convert transformed to model variables for analysis and trial ensembles
755 call utl_tmg_start(3,'--WriteEnsemble')
756 if (.not. outputOnlyEnsMean) then
757 call ens_writeEnsemble(ensembleAnl, '.', '', etiket_anl, 'A', &
758 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
759 containsFullField_opt=.true.)
760 end if
761 call utl_tmg_stop(3)
762
763 !- Output the sub-sampled ensemble analyses and increments
764 if (writeSubSample) then
765 call utl_tmg_start(5,'--WriteEnsMeanRms')
766 ! Output the ensemble mean increment (include MeanAnl Psfc)
767 call fln_ensAnlFileName( outFileName, 'subspace', tim_getDateStamp(), 0, ensFileNameSuffix_opt='inc' )
768 ! here we assume 4 digits for the ensemble member!!!!
769 etiket = trim(etiket_inc) // '0000'
770 do stepIndex = 1, tim_nstepobsinc
771 call gio_writeToFile(stateVectorMeanIncSubSample, outFileName, etiket, &
772 typvar_opt='R', writeHeightSfc_opt=.false., numBits_opt=numBits, &
773 stepIndex_opt=stepIndex, containsFullField_opt=.false.)
774 if (gsv_isAllocated(stateVectorMeanAnlSfcPres)) then
775 call gio_writeToFile(stateVectorMeanAnlSfcPres, outFileName, etiket, &
776 typvar_opt='A', writeHeightSfc_opt=.true., &
777 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
778 end if
779 end do
780
781 ! Output the ensemble mean analysis state
782 call fln_ensAnlFileName( outFileName, 'subspace', tim_getDateStamp(), 0 )
783 ! here we assume 4 digits for the ensemble member!!!!
784 etiket = trim(etiket_anl) // '0000'
785 do stepIndex = 1, tim_nstepobsinc
786 call gio_writeToFile(stateVectorMeanAnlSubSample, outFileName, etiket, &
787 typvar_opt='A', writeHeightSfc_opt=.false., numBits_opt=numBits, &
788 stepIndex_opt=stepIndex, containsFullField_opt=.true.)
789 end do
790 call utl_tmg_stop(5)
791
792 ! Output the sub-sampled analysis ensemble members
793 call utl_tmg_start(3,'--WriteEnsemble')
794 if (.not. outputOnlyEnsMean) then
795 call ens_writeEnsemble(ensembleAnlSubSample, 'subspace', '', etiket_anl, 'A', &
796 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
797 containsFullField_opt=.true.)
798 end if
799 call utl_tmg_stop(3)
800
801 ! Output the sub-sampled ensemble increments (include MeanAnl Psfc)
802 call utl_tmg_start(3,'--WriteEnsemble')
803 if (.not. outputOnlyEnsMean) then
804 call ens_writeEnsemble(ensembleAnlIncSubSample, 'subspace', '', etiket_inc, 'R', &
805 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
806 containsFullField_opt=.false., resetTimeParams_opt=.true.)
807 ! Also write the reference (analysis) surface pressure to increment files
808 if (gsv_isAllocated(stateVectorMeanAnlSfcPresMpiGlb)) then
809 call epp_writeToAllMembers(stateVectorMeanAnlSfcPresMpiGlb, &
810 ens_getNumMembers(ensembleAnlSubSample), &
811 etiket=etiket_inc, typvar='A', fileNameSuffix='inc', &
812 ensPath='subspace')
813 end if
814 end if
815 call utl_tmg_stop(3)
816
817 end if ! writeSubSample
818
819 !- Output the unperturbed sub-sampled ensemble analyses
820 if (writeSubSampleUnPert) then
821
822 ! Output the sub-sampled analysis ensemble members
823 call utl_tmg_start(3,'--WriteEnsemble')
824 if (.not. outputOnlyEnsMean) then
825 call ens_writeEnsemble(ensembleAnlSubSampleUnPert, 'subspace_unpert', '', etiket_anl, 'A', &
826 numBits_opt=16, etiketAppendMemberNumber_opt=.true., &
827 containsFullField_opt=.true.)
828 end if
829 call utl_tmg_stop(3)
830
831 end if
832
833 end if ! ens_isAllocated(ensembleAnl)
834
835 end subroutine epp_postProcess
836
837 !----------------------------------------------------------------------
838 ! epp_writeToAllMembers (private subroutine)
839 !----------------------------------------------------------------------
840 subroutine epp_writeToAllMembers(stateVector, nEns, etiket, typvar, &
841 fileNameSuffix, ensPath)
842 !
843 !:Purpose: Write the contents of the supplied stateVector to all
844 ! ensemble member files in an efficient parallel way.
845 !
846 implicit none
847
848 ! Arguments:
849 type(struct_gsv), intent(in) :: stateVector
850 integer , intent(in) :: nEns
851 character(len=*), intent(in) :: etiket
852 character(len=*), intent(in) :: typvar
853 character(len=*), intent(in) :: fileNameSuffix
854 character(len=*), intent(in) :: ensPath
855
856 ! Locals:
857 integer :: memberIndex, stepIndex, writeFilePE(nEns)
858 character(len=4) :: memberIndexStr
859 character(len=256) :: outFileName
860
861 do memberIndex = 1, nEns
862 writeFilePE(memberIndex) = mod(memberIndex-1, mmpi_nprocs)
863 end do
864
865 do memberIndex = 1, nEns
866
867 if (mmpi_myid == writeFilePE(memberIndex)) then
868
869 call fln_ensAnlFileName( outFileName, ensPath, tim_getDateStamp(), &
870 memberIndex, ensFileNameSuffix_opt=fileNameSuffix )
871 write(memberIndexStr,'(I4.4)') memberIndex
872
873 do stepIndex = 1, tim_nstepobsinc
874 call gio_writeToFile(stateVector, outFileName, &
875 trim(etiket) // memberIndexStr, &
876 typvar_opt=trim(typvar), writeHeightSfc_opt=.true., &
877 stepIndex_opt=stepIndex, containsFullField_opt=.true., &
878 numBits_opt=16)
879 end do
880
881 end if
882
883 end do
884
885 end subroutine epp_writeToAllMembers
886
887 !--------------------------------------------------------------------------
888 ! epp_RTPS
889 !--------------------------------------------------------------------------
890 subroutine epp_RTPS(ensembleAnl, stateVectorStdDevAnl, stateVectorStdDevTrl, &
891 stateVectorMeanAnl, alphaRTPS)
892 !:Purpose: Apply Relaxation To Prior Spread ensemble inflation according
893 ! to the factor alphaRTPS (usually between 0 and 1).
894 implicit none
895
896 ! Arguments:
897 type(struct_ens), intent(inout) :: ensembleAnl
898 type(struct_gsv), intent(in) :: stateVectorStdDevAnl
899 type(struct_gsv), intent(in) :: stateVectorStdDevTrl
900 type(struct_gsv), intent(in) :: stateVectorMeanAnl
901 real(8) , intent(in) :: alphaRTPS
902
903 ! Locals:
904 integer :: varLevIndex, latIndex, lonIndex, stepIndex, memberIndex
905 integer :: nEns, numVarLev, myLonBeg, myLonEnd, myLatBeg, myLatEnd
906 real(8) :: factorRTPS
907 real(4), pointer :: stdDevTrl_ptr_r4(:,:,:,:), stdDevAnl_ptr_r4(:,:,:,:)
908 real(4), pointer :: meanAnl_ptr_r4(:,:,:,:), memberAnl_ptr_r4(:,:,:,:)
909
910 write(*,*) 'epp_RTPS: Starting'
911
912 call gsv_getField(stateVectorStdDevTrl,stdDevTrl_ptr_r4)
913 call gsv_getField(stateVectorStdDevAnl,stdDevAnl_ptr_r4)
914 call gsv_getField(stateVectorMeanAnl,meanAnl_ptr_r4)
915
916 nEns = ens_getNumMembers(ensembleAnl)
917 numVarLev = ens_getNumK(ensembleAnl)
918 call ens_getLatLonBounds(ensembleAnl, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
919 do varLevIndex = 1, numVarLev
920 memberAnl_ptr_r4 => ens_getOneLev_r4(ensembleAnl,varLevIndex)
921 do latIndex = myLatBeg, myLatEnd
922 do lonIndex = myLonBeg, myLonEnd
923 do stepIndex = 1, tim_nstepobsinc
924 ! compute the inflation factor for RTPS
925 if ( stdDevAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) > 0.0 ) then
926 factorRTPS = 1.0D0 + alphaRTPS * &
927 ( stdDevTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) - &
928 stdDevAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) ) / &
929 stdDevAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex)
930 else
931 factorRTPS = 0.0D0
932 end if
933 ! apply the inflation factor to all Anl members (in place)
934 if (factorRTPS > 0.0D0) then
935 do memberIndex = 1, nEns
936 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
937 factorRTPS * &
938 ( memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
939 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) ) + &
940 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex)
941 end do ! memberIndex
942 end if ! factorRTPS > 0
943 end do ! stepIndex
944 end do ! lonIndex
945 end do ! latIndex
946 end do ! varLevIndex
947
948 write(*,*) 'epp_RTPS: Finished'
949
950 end subroutine epp_RTPS
951
952 !--------------------------------------------------------------------------
953 ! epp_RTPP
954 !--------------------------------------------------------------------------
955 subroutine epp_RTPP(ensembleAnl, ensembleTrl, stateVectorMeanAnl, &
956 stateVectorMeanTrl, alphaRTPP)
957 !:Purpose: Apply Relaxation To Prior Perturbation ensemble inflation according
958 ! to the factor alphaRTPP (usually between 0 and 1).
959 implicit none
960
961 ! Arguments:
962 type(struct_ens), intent(inout) :: ensembleAnl
963 type(struct_ens), intent(inout) :: ensembleTrl
964 type(struct_gsv), intent(in) :: stateVectorMeanAnl
965 type(struct_gsv), intent(in) :: stateVectorMeanTrl
966 real(8) , intent(in) :: alphaRTPP
967
968 ! Locals:
969 integer :: varLevIndex, latIndex, lonIndex, stepIndex, memberIndex
970 integer :: nEns, numVarLev, myLonBeg, myLonEnd, myLatBeg, myLatEnd
971 real(4), pointer :: meanAnl_ptr_r4(:,:,:,:), meanTrl_ptr_r4(:,:,:,:)
972 real(4), pointer :: memberAnl_ptr_r4(:,:,:,:), memberTrl_ptr_r4(:,:,:,:)
973
974 write(*,*) 'epp_RTPP: Starting'
975
976 call gsv_getField(stateVectorMeanAnl, meanAnl_ptr_r4)
977 call gsv_getField(stateVectorMeanTrl, meanTrl_ptr_r4)
978
979 nEns = ens_getNumMembers(ensembleAnl)
980 numVarLev = ens_getNumK(ensembleAnl)
981 call ens_getLatLonBounds(ensembleAnl, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
982
983 do varLevIndex = 1, numVarLev
984 memberAnl_ptr_r4 => ens_getOneLev_r4(ensembleAnl,varLevIndex)
985 memberTrl_ptr_r4 => ens_getOneLev_r4(ensembleTrl,varLevIndex)
986 do latIndex = myLatBeg, myLatEnd
987 do lonIndex = myLonBeg, myLonEnd
988 do stepIndex = 1, tim_nstepobsinc
989 ! apply RTPP to all Anl members (in place)
990 ! member_anl = mean_anl + (1-a)*(member_anl-mean_anl) + a*(member_trl-mean_trl)
991 do memberIndex = 1, nEns
992 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
993 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) + &
994 (1.0D0 - alphaRTPP) * &
995 ( memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
996 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) ) + &
997 alphaRTPP * &
998 ( memberTrl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
999 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) )
1000 end do ! memberIndex
1001 end do ! stepIndex
1002 end do ! lonIndex
1003 end do ! latIndex
1004 end do ! varLevIndex
1005
1006 write(*,*) 'epp_RTPP: Finished'
1007
1008 end subroutine epp_RTPP
1009
1010 !--------------------------------------------------------------------------
1011 ! epp_addRandomPert
1012 !--------------------------------------------------------------------------
1013 subroutine epp_addRandomPert(ensembleAnl, stateVectorRefState, alphaRandomPert, &
1014 randomSeed, useMemberAsHuRefState)
1015 !:Purpose: Apply additive inflation using random perturbations from sampling
1016 ! the B matrix as defined by the regular namelist block NAMBHI, NAMBEN, etc.
1017 ! The scale factor alphaRandomPert (usually between 0 and 1) is used
1018 ! to simply multiply the resulting perturbations before adding to the original
1019 ! members. The perturbations have zero ensemble mean.
1020 implicit none
1021
1022 ! Arguments:
1023 type(struct_ens), intent(inout) :: ensembleAnl
1024 type(struct_gsv), intent(in) :: stateVectorRefState
1025 real(8) , intent(in) :: alphaRandomPert
1026 integer , intent(in) :: randomSeed
1027 logical , intent(in) :: useMemberAsHuRefState
1028
1029 ! Locals:
1030 type(struct_gsv) :: stateVectorPerturbation
1031 type(struct_gsv) :: stateVectorPerturbationInterp
1032 type(struct_gsv) :: stateVectorHuRefState
1033 type(struct_gsv) :: stateVectorHuRefStateInterp
1034 type(struct_gsv) :: stateVectorP0Ref
1035 type(struct_vco), pointer :: vco_randomPert, vco_ens
1036 type(struct_hco), pointer :: hco_randomPert, hco_ens, hco_core
1037 character(len=12) :: etiket
1038 real(8), allocatable :: controlVector_mpiglobal(:), controlVector(:)
1039 real(8), allocatable :: perturbationMean(:,:,:)
1040 real(8), pointer :: PsfcRef(:,:,:,:)
1041 real(8), pointer :: perturbation_ptr(:,:,:)
1042 real(4), pointer :: memberAnl_ptr_r4(:,:,:,:)
1043 integer :: cvIndex, memberIndex, varLevIndex, lonIndex, latIndex, stepIndex
1044 integer :: nEns, numVarLev, myLonBeg, myLonEnd, myLatBeg, myLatEnd, varIndex
1045 integer :: middleStepIndex
1046 logical, save :: firstCall = .true.
1047 character(len=4), pointer :: varNamesWithLQ(:)
1048 character(len=4) :: varName
1049
1050 call utl_tmg_start(4,'--AddEnsRandomPert')
1051
1052 ! Determine middle timestep
1053 middleStepIndex = (tim_nstepobsinc + 1) / 2
1054
1055 ! Get ensemble dimensions
1056 nEns = ens_getNumMembers(ensembleAnl)
1057 numVarLev = ens_getNumK(ensembleAnl)
1058 call ens_getLatLonBounds(ensembleAnl, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
1059 vco_ens => ens_getVco(ensembleAnl)
1060 hco_ens => ens_getHco(ensembleAnl)
1061
1062 ! Define the horiz/vertical coordinate for perturbation calculation
1063 nullify(vco_randomPert)
1064 nullify(hco_randomPert)
1065 call hco_setupFromFile(hco_randomPert, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
1066 if ( hco_randomPert % global ) then
1067 etiket = 'BGCK_STDDEV'
1068 else
1069 etiket = 'STDDEV'
1070 end if
1071 call vco_setupFromFile(vco_randomPert, './bgcov', etiket)
1072
1073 ! Special modification to vco_randomPert if a soil variables is present
1074 if ( vco_getNumLev(vco_ens,'OT','I0') > 0 .or. &
1075 vco_getNumLev(vco_ens,'OT','I1') > 0 ) then
1076 write(*,*)
1077 write(*,*) 'epp_addRandomPert: copying number of levels for land surface variables!'
1078 write(*,*)
1079 vco_randomPert%nlev_Other(:) = vco_ens%nlev_Other(:)
1080 end if
1081
1082 ! Get list of variable names in the ensemble and modify if necessary
1083 nullify(varNamesWithLQ)
1084 call ens_varNamesList(varNamesWithLQ,ensembleAnl)
1085 if (useMemberAsHuRefState .and. ens_varExist(ensembleAnl,'HU')) then
1086 ! Replace HU with LQ so we can apply the transform directly in this routine
1087 do varIndex = 1, size(varNamesWithLQ)
1088 if (varNamesWithLQ(varIndex) == 'HU') varNamesWithLQ(varIndex) = 'LQ'
1089 end do
1090 if (mmpi_myid == 0) write(*,*) 'epp_addRandomPert: varNamesWithLQ = ', varNamesWithLQ(:)
1091 end if
1092
1093 hco_core => hco_randomPert
1094 if (firstCall) then
1095 call utl_tmg_stop(4) ! stop counter, since Bmat has it's own counters
1096 call bmat_setup(hco_randomPert, hco_core, vco_randomPert)
1097 firstCall = .false.
1098 call utl_tmg_start(4,'--AddEnsRandomPert')
1099 end if
1100 call gvt_setup(hco_randomPert, hco_core, vco_randomPert)
1101
1102 call rng_setup(abs(randomSeed))
1103
1104 allocate(controlVector(cvm_nvadim))
1105 allocate(controlVector_mpiglobal(cvm_nvadim_mpiglobal))
1106 allocate(perturbationMean(myLonBeg:myLonEnd,myLatBeg:myLatEnd,numVarLev))
1107 perturbationMean(:,:,:) = 0.0d0
1108
1109 ! NOTE: The following stateVectors will include LQ instead of HU when
1110 ! useMemberAsHuRefState=true. This results in the perturbations
1111 ! being provided by bmat_sqrtB in terms of LQ, allowing the conversion
1112 ! from LQ to HU to be done within this subroutine using the HU field
1113 ! of each ensemble member.
1114 call gsv_allocate(stateVectorPerturbation, 1, hco_randomPert, vco_randomPert, &
1115 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
1116 hInterpolateDegree_opt='LINEAR', &
1117 varNames_opt=varNamesWithLQ)
1118 call gsv_allocate(stateVectorPerturbationInterp, 1, hco_ens, vco_ens, &
1119 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
1120 hInterpolateDegree_opt='LINEAR', &
1121 varNames_opt=varNamesWithLQ)
1122 call gsv_getField(stateVectorPerturbationInterp,perturbation_ptr)
1123
1124 ! allocate a minimal reference P0
1125 call gsv_allocate(statevectorP0Ref, 1, hco_ens, vco_randomPert, &
1126 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
1127 hInterpolateDegree_opt='LINEAR', &
1128 varNames_opt=(/'P0'/))
1129 call gsv_getField(statevectorP0Ref, PsfcRef, 'P0')
1130 PsfcRef(:,:,:,:) = 100000.0D0
1131
1132
1133 ! prepare the reference state HU field for transforming LQ to HU perturbations
1134 if (ens_varExist(ensembleAnl,'HU')) then
1135 call gsv_allocate(stateVectorHuRefState, 1, hco_ens, vco_ens, &
1136 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
1137 allocHeightSfc_opt=.true., hInterpolateDegree_opt='LINEAR', &
1138 hExtrapolateDegree_opt='MINIMUM', &
1139 varNames_opt=(/'HU','P0'/) )
1140 if (.not. useMemberAsHuRefState) then
1141 ! use the single provided state as the reference state and interpolate to perturbation grid
1142 call gsv_copy(stateVectorRefState, stateVectorHuRefState, allowVarMismatch_opt=.true.)
1143 call gsv_allocate(stateVectorHuRefStateInterp, 1, hco_randomPert, vco_randomPert, &
1144 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
1145 allocHeightSfc_opt=.true., hInterpolateDegree_opt='LINEAR', &
1146 hExtrapolateDegree_opt='MINIMUM', &
1147 varNames_opt=(/'HU','P0'/) )
1148 call int_interp_gsv(stateVectorHuRefState, stateVectorHuRefStateInterp)
1149 call gsv_deallocate(stateVectorHuRefState)
1150 end if
1151 end if
1152
1153 do memberIndex = 1, nEns
1154
1155 if( mmpi_myid == 0 ) then
1156 write(*,*)
1157 write(*,*) 'Computing random perturbation number= ', memberIndex
1158 end if
1159
1160 ! global vector random control vector (independent of mpi topology)
1161 do cvIndex = 1, cvm_nvadim_mpiglobal
1162 controlVector_mpiglobal(cvIndex) = rng_gaussian()
1163 end do
1164 call bmat_reduceToMPILocal( controlVector, controlVector_mpiglobal )
1165
1166 call utl_tmg_stop(4) ! stop counter, since Bmat has it's own counters
1167 if (ens_varExist(ensembleAnl,'HU') .and. .not.useMemberAsHuRefState) then
1168 ! Use supplied reference state for LQ to HU conversion
1169 call bmat_sqrtB(controlVector, cvm_nvadim, & ! IN
1170 stateVectorPerturbation, & ! OUT
1171 stateVectorRef_opt=stateVectorHuRefStateInterp) ! IN
1172 else
1173 ! No conversion from LQ to HU done in bmat_sqrtB
1174 call bmat_sqrtB(controlVector, cvm_nvadim, & ! IN
1175 stateVectorPerturbation) ! OUT
1176 end if
1177 call utl_tmg_start(4,'--AddEnsRandomPert')
1178
1179 call int_interp_gsv(stateVectorPerturbation, stateVectorPerturbationInterp, &
1180 statevectorRef_opt=statevectorP0Ref)
1181
1182 ! If desired, use member itself as reference state for LQ to HU conversion
1183 if (ens_varExist(ensembleAnl,'HU') .and. useMemberAsHuRefState) then
1184 call ens_copyMember(ensembleAnl, stateVectorHuRefState, memberIndex)
1185 call gvt_transform( stateVectorPerturbationInterp, & ! INOUT
1186 'LQtoHU_tlm', & ! IN
1187 stateVectorRef_opt=stateVectorHuRefState ) ! IN
1188 end if
1189
1190 ! scale the perturbation by the specified factor
1191 call gsv_scale(stateVectorPerturbationInterp, alphaRandomPert)
1192
1193 write(*,*) 'epp_addRandomPert: member ', memberIndex, ', perturbation min/maxval = ', &
1194 minval(perturbation_ptr), maxval(perturbation_ptr)
1195
1196 do varLevIndex = 1, numVarLev
1197 varName = ens_getVarNameFromK(ensembleAnl,varLevIndex)
1198 memberAnl_ptr_r4 => ens_getOneLev_r4(ensembleAnl,varLevIndex)
1199 do latIndex = myLatBeg, myLatEnd
1200 do lonIndex = myLonBeg, myLonEnd
1201
1202 perturbationMean(lonIndex, latIndex, varLevIndex) = &
1203 perturbationMean(lonIndex, latIndex, varLevIndex) + &
1204 perturbation_ptr(lonIndex, latIndex, varLevIndex) / real(nEns, 8)
1205
1206 ! Add perturbation to member
1207 do stepIndex = 1, tim_nstepobsinc
1208 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
1209 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) + &
1210 perturbation_ptr(lonIndex, latIndex, varLevIndex)
1211 end do ! stepIndex
1212
1213 end do ! lonIndex
1214 end do ! latIndex
1215 end do ! varLevIndex
1216
1217 end do ! memberIndex
1218
1219 ! remove the ensemble mean of the perturbations
1220 do varLevIndex = 1, numVarLev
1221 memberAnl_ptr_r4 => ens_getOneLev_r4(ensembleAnl,varLevIndex)
1222 do latIndex = myLatBeg, myLatEnd
1223 do lonIndex = myLonBeg, myLonEnd
1224 do stepIndex = 1, tim_nstepobsinc
1225 do memberIndex = 1, nEns
1226 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
1227 memberAnl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
1228 perturbationMean(lonIndex, latIndex, varLevIndex)
1229 end do ! memberIndex
1230 end do ! stepIndex
1231 end do ! lonIndex
1232 end do ! latIndex
1233 end do ! varLevIndex
1234
1235 deallocate(controlVector)
1236 deallocate(controlVector_mpiglobal)
1237 deallocate(perturbationMean)
1238 call gsv_deallocate(stateVectorPerturbation)
1239 call gsv_deallocate(stateVectorPerturbationInterp)
1240 call gsv_deallocate(stateVectorP0Ref)
1241 if (gsv_isAllocated(stateVectorHuRefStateInterp)) then
1242 call gsv_deallocate(stateVectorHuRefStateInterp)
1243 end if
1244
1245 call utl_tmg_stop(4)
1246
1247 end subroutine epp_addRandomPert
1248
1249 !-----------------------------------------------------------------
1250 ! epp_selectSubSample
1251 !-----------------------------------------------------------------
1252 subroutine epp_selectSubSample(ensembleAnl, ensembleAnlSubSample, &
1253 ensembleTrl_opt, ensembleTrlSubSample_opt)
1254 !:Purpose: Create sub-sampled ensembles of analyses and trials based on
1255 ! the contents of the ascii files 'sampletable' which lists
1256 ! the member indices for the subsample.
1257 implicit none
1258
1259 ! Arguments:
1260 type(struct_ens) , intent(in) :: ensembleAnl
1261 type(struct_ens) , intent(out) :: ensembleAnlSubSample
1262 type(struct_ens), optional, intent(in) :: ensembleTrl_opt
1263 type(struct_ens), optional, intent(out) :: ensembleTrlSubSample_opt
1264
1265 ! Locals:
1266 type(struct_gsv) :: stateVectorMember
1267 integer :: nulFile, ierr, status, numSubSample
1268 integer :: memberIndex, memberIndexSubSample, memberIndexFull
1269 integer :: memberIndexesSubSample(1000), memberIndexesFull(1000)
1270 integer, allocatable :: dateStampListInc(:)
1271 integer, external :: fnom, fclos
1272
1273 numSubSample = 0
1274
1275 nulFile = 0
1276 ierr = fnom(nulFile, './sampletable', 'FTN+SEQ+R/O', 0)
1277 do
1278 read(nulFile,*, IOSTAT=status) memberIndexSubSample, memberIndexFull
1279 if (status < 0) exit
1280
1281 numSubSample = numSubSample + 1
1282 write(*,*) 'epp_selectSubSample: ', memberIndexSubSample, memberIndexFull
1283 memberIndexesSubSample(numSubSample) = memberIndexSubSample
1284 memberIndexesFull(numSubSample) = memberIndexFull
1285 end do
1286 ierr = fclos(nulFile)
1287
1288 write(*,*) 'epp_selectSubSample: number of subSample members = ', numSubSample
1289
1290 allocate(dateStampListInc(tim_nstepobsinc))
1291 call tim_getstamplist(dateStampListInc,tim_nstepobsinc,tim_getDatestamp())
1292
1293 call ens_allocate(ensembleAnlSubSample, numSubSample, tim_nstepobsinc, &
1294 ens_getHco(ensembleAnl), ens_getVco(ensembleAnl), dateStampListInc)
1295 if (present(ensembleTrlSubSample_opt)) then
1296 call ens_allocate(ensembleTrlSubSample_opt, numSubSample, tim_nstepobsinc, &
1297 ens_getHco(ensembleAnl), ens_getVco(ensembleAnl), dateStampListInc)
1298 end if
1299
1300 call gsv_allocate(stateVectorMember, tim_nstepobsinc, &
1301 ens_getHco(ensembleAnl), ens_getVco(ensembleAnl), &
1302 dateStamp_opt=tim_getDateStamp(), &
1303 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
1304 dataKind_opt=4, allocHeightSfc_opt=.false., &
1305 allocHeight_opt=.false., allocPressure_opt=.false. )
1306
1307 do memberIndex = 1, numSubSample
1308 ! copy analysis ensemble member
1309 call ens_copyMember(ensembleAnl, stateVectorMember, &
1310 memberIndexesFull(memberIndex))
1311 call ens_insertMember(ensembleAnlSubSample, stateVectorMember, &
1312 memberIndexesSubSample(memberIndex))
1313
1314 if (present(ensembleTrl_opt) .and. present(ensembleTrlSubSample_opt)) then
1315 ! copy trial ensemble member
1316 call ens_copyMember(ensembleTrl_opt, stateVectorMember, &
1317 memberIndexesFull(memberIndex))
1318 call ens_insertMember(ensembleTrlSubSample_opt, stateVectorMember, &
1319 memberIndexesSubSample(memberIndex))
1320 end if
1321 end do
1322
1323 call gsv_deallocate(stateVectorMember)
1324 deallocate(dateStampListInc)
1325
1326 end subroutine epp_selectSubSample
1327
1328 !-----------------------------------------------------------------
1329 ! epp_hybridRecentering
1330 !-----------------------------------------------------------------
1331 subroutine epp_hybridRecentering(ensembleAnl, weightRecenter, weightRecenterLand, &
1332 useOptionTableRecenter, numMembersToRecenter)
1333 !:Purpose: Modify an ensemble by recentering the members on a state provided
1334 ! in the file "recentering_analysis".
1335 ! The "weightRecenter" and "numMembersToRecenter" are used in the calculation
1336 ! to determine the amount of recentering and how many members it is
1337 ! applied to. Alternatively the information in the "optiontable" file
1338 ! can be used to perform a different amount of recentering on each
1339 ! member.
1340 implicit none
1341
1342 ! Arguments:
1343 type(struct_ens), intent(inout) :: ensembleAnl
1344 real(8) , intent(in) :: weightRecenter(:)
1345 real(8) , intent(in) :: weightRecenterLand
1346 logical , intent(in) :: useOptionTableRecenter
1347 integer , intent(in) :: numMembersToRecenter
1348
1349 ! Locals:
1350 type(struct_gsv) :: stateVectorRecenterAnl
1351 type(struct_hco), pointer :: hco_ens => null()
1352 type(struct_vco), pointer :: vco_ens => null()
1353 character(len=30) :: recenterAnlFileName = 'recentering_analysis'
1354 character(len=20) :: stringArray(100)
1355 character(len=1000) :: textLine
1356 integer :: stepIndex, memberIndex, columnIndex
1357 integer :: numMembers, numColumns, nulFile, status
1358 logical :: recenterAnlFileExists
1359 real(8), allocatable :: weightArray(:,:)
1360 real(8), allocatable :: weightArrayLand(:)
1361 real(8) :: weightFound
1362 integer, external :: fnom, fclos
1363
1364 ! check if recentering analysis file exists
1365 inquire(file=recenterAnlFileName, exist=recenterAnlFileExists)
1366 if (.not. recenterAnlFileExists) then
1367 write(*,*) 'epp_hybridRecentering: RecenterAnlFileName = ', recenterAnlFileName
1368 call utl_abort('epp_hybridRecentering: The recentering analysis file does not exist')
1369 end if
1370
1371 hco_ens => ens_getHco(ensembleAnl)
1372 vco_ens => ens_getVco(ensembleAnl)
1373 numMembers = ens_getNumMembers(ensembleAnl)
1374
1375 allocate( weightArray(vco_maxNumLevels,0:numMembers) )
1376 allocate( weightArrayLand(0:numMembers) )
1377
1378 ! read the optiontable file, if requested
1379 if (useOptionTableRecenter) then
1380 write(*,*) 'epp_hybridRecentering: using optiontable file to specify recentering weights.'
1381 nulFile = 0
1382 status = fnom(nulFile, './optiontable', 'FMT+SEQ+R/O', 0)
1383 read(nulFile,'(a)', IOSTAT=status) textLine
1384 if (status /= 0) then
1385 call utl_abort('epp_hybridRecentering: unable to read optiontable file')
1386 end if
1387 call utl_parseColumns(textLine, numColumns)
1388 if (mmpi_myid==0) write(*,*) 'epp_hybridRecentering: optiontable file has ', numColumns, ' columns.'
1389 rewind(nulFile)
1390 do memberIndex = 0, numMembers
1391 read(nulFile,'(a)') textLine
1392 call utl_parseColumns(textLine, numColumns, stringArray_opt=stringArray)
1393 if (mmpi_myid==0) write(*,*) memberIndex, (stringArray(columnIndex),columnIndex=1,numColumns)
1394 read(stringArray(numColumns),'(f6.3)') weightFound
1395 weightArray(:,memberIndex) = weightFound
1396 if (mmpi_myid==0) write(*,*) 'weightArray = ', weightArray(1,memberIndex)
1397 end do
1398 status = fclos(nulFile)
1399 else
1400 if (any(weightRecenter(1:vco_getNumLev(vco_ens,'MM')) == -1.d0)) then
1401 write(*,*) ''
1402 write(*,*) 'Number of weightRecenter coefficients needed = ', vco_getNumLev(vco_ens,'MM')
1403 write(*,*) 'Provided values : '
1404 write(*,*) weightRecenter(1:vco_getNumLev(vco_ens,'MM'))
1405 call utl_abort('epp_hybridRecentering: A valid weightRecenter coefficient was not provided for all the vertical levels')
1406 end if
1407 do memberIndex = 0, numMembers
1408 weightArray(:,memberIndex) = weightRecenter(:)
1409 end do
1410 end if
1411
1412 do memberIndex = 0, numMembers
1413 weightArrayLand(memberIndex) = weightRecenterLand
1414 end do
1415
1416 ! allocate and read in recentering analysis state
1417 call gsv_allocate( stateVectorRecenterAnl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
1418 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
1419 dataKind_opt=4, allocHeightSfc_opt=.false., &
1420 hInterpolateDegree_opt = 'LINEAR', &
1421 allocHeight_opt=.false., allocPressure_opt=.false. )
1422 call gsv_zero(stateVectorRecenterAnl)
1423
1424 do stepIndex = 1, tim_nstepobsinc
1425 call gio_readFromFile( stateVectorRecenterAnl, recenterAnlFileName, ' ', ' ', &
1426 stepIndex_opt=stepIndex, containsFullField_opt=.true., &
1427 readHeightSfc_opt=.false. )
1428 end do
1429
1430 ! apply recentering
1431 call ens_recenter(ensembleAnl, stateVectorRecenterAnl, &
1432 recenteringCoeff_opt=weightArray(:,1:numMembers), &
1433 numMembersToRecenter_opt=numMembersToRecenter, &
1434 recenteringCoeffLand_opt=weightArrayLand)
1435
1436 call gsv_deallocate(stateVectorRecenterAnl)
1437 if ( allocated(weightArray) ) deallocate(weightArray)
1438
1439 end subroutine epp_hybridRecentering
1440
1441 !-----------------------------------------------------------------
1442 ! epp_printRmsStats
1443 !-----------------------------------------------------------------
1444 subroutine epp_printRmsStats(stateVectorStdDev,fileName,elapsed,ftype,nEns)
1445 !
1446 !:Purpose: Print statistics of a field to an ASCII output file
1447 !
1448 implicit none
1449
1450 ! Arguments:
1451 type(struct_gsv), intent(in) :: stateVectorStdDev
1452 character(len=*), intent(in) :: fileName
1453 real(8), intent(in) :: elapsed
1454 character(len=1), intent(in) :: ftype
1455 integer, intent(in) :: nEns
1456
1457 ! Locals:
1458 real(8), allocatable :: rmsvalue(:)
1459 type(struct_vco), pointer :: vco
1460 type(struct_hco), pointer :: hco
1461 character(len=4), allocatable :: nomvar_v(:)
1462 character(len=4) :: varLevel
1463 real(8), allocatable :: scaleFactor(:)
1464 real(4), allocatable :: pressureOrDepth(:)
1465 integer :: ierr, latIndex, lonIndex, nulFile
1466 integer :: kIndex, kIndexCount, levIndex, numK, nLev_M, kIndexUU, kIndexVV
1467 real(4), pointer :: stdDev_ptr_r4(:,:,:)
1468 real(8) :: pSfc(1,1)
1469 real(8), pointer :: pressures_T(:,:,:), pressures_M(:,:,:)
1470 integer, external :: fnom, fclos
1471 real(8), save, allocatable :: weight(:,:)
1472 logical, save :: firstCall = .true.
1473
1474 vco => gsv_getVco(stateVectorStdDev)
1475 nLev_M = vco_getNumLev(vco,'MM')
1476
1477 numK = gsv_getNumK(stateVectorStdDev)
1478 allocate(nomvar_v(numK))
1479 allocate(pressureOrDepth(numK))
1480 allocate(rmsvalue(numK))
1481 allocate(scaleFactor(numK))
1482
1483 ! compute a 2D weight field used for horizontal averaging
1484 hco => gsv_getHco(stateVectorStdDev)
1485 if (firstCall) then
1486 ! weights are reduced in the overlap region in case of a Yin-Yang grid
1487 allocate(weight(stateVectorStdDev%ni,stateVectorStdDev%nj))
1488 weight(:,:) = 0.0d0
1489 call hco_weight(hco, weight)
1490 firstCall = .false.
1491 end if
1492 ! compute global mean variance accounting for weights
1493 call gsv_getField(stateVectorStdDev,stdDev_ptr_r4)
1494 rmsvalue(:) = 0.0D0
1495 do kIndex = 1, numK
1496 do latIndex = stateVectorStdDev%myLatBeg, stateVectorStdDev%myLatEnd
1497 do lonIndex = stateVectorStdDev%myLonBeg, stateVectorStdDev%myLonEnd
1498 rmsvalue(kIndex) = rmsvalue(kIndex) + &
1499 (stdDev_ptr_r4(lonIndex,latIndex,kIndex)**2) * &
1500 weight(lonIndex,latIndex)
1501 end do
1502 end do
1503 call mmpi_allreduce_sumreal8scalar(rmsvalue(kIndex),'GRID')
1504 rmsvalue(kIndex) = rmsvalue(kIndex)**0.5
1505 end do
1506
1507 if (vco%vgridPresent) then
1508 ! compute pressure for a column where Psfc=1000hPa
1509 pSfc(1,1) = 1000.0D2 !1000 hPa
1510 ! pressure levels
1511 call czp_fetch3DLevels(vco, pSfc, fldM_opt=pressures_M, fldT_opt=pressures_T)
1512
1513 ! set the variable name and pressure for each element of column
1514 do kIndex = 1, numK
1515 levIndex = gsv_getLevFromK(stateVectorStdDev, kIndex)
1516 nomvar_v(kIndex) = gsv_getVarNameFromK(stateVectorStdDev,kIndex)
1517 varLevel = vnl_varLevelFromVarname(nomvar_v(kIndex))
1518 if (varLevel == 'MM') then
1519 pressureOrDepth(kIndex) = pressures_M(1,1,levIndex)/Psfc(1,1)
1520 else if (varLevel == 'TH') then
1521 pressureOrDepth(kIndex) = pressures_T(1,1,levIndex)/Psfc(1,1)
1522 else
1523 pressureOrDepth(kIndex) = 1.0
1524 end if
1525 end do
1526 deallocate(pressures_M)
1527 deallocate(pressures_T)
1528
1529 else if (vco%nLev_depth > 0) then
1530
1531 ! set the variable name and depth for each element of column
1532 do kIndex = 1, numK
1533 nomvar_v(kIndex) = gsv_getVarNameFromK(stateVectorStdDev,kIndex)
1534 if (vnl_varLevelFromVarName(nomvar_v(kIndex)) == 'SS') then
1535 pressureOrDepth(kIndex) = 0.0
1536 else
1537 levIndex = gsv_getLevFromK(stateVectorStdDev, kIndex)
1538 pressureOrDepth(kIndex) = vco%depths(levIndex)
1539 end if
1540 end do
1541
1542 else
1543
1544 call utl_abort('epp_printRmsStats: Unknown type of levels')
1545
1546 end if
1547
1548 ! set the scaleFactor and rmsvalue for each element of column
1549 do kIndex = 1, numK
1550 if ( (nomvar_v(kIndex) == 'UU') .or. (nomvar_v(kIndex) == 'VV') ) then
1551 scaleFactor(kIndex) = MPC_KNOTS_PER_M_PER_S_R8
1552 else if (nomvar_v(kIndex) == 'P0') then
1553 scaleFactor(kIndex) = MPC_MBAR_PER_PA_R8
1554 else
1555 scaleFactor(kIndex) = 1.0d0
1556 end if
1557 rmsvalue(kIndex) = scaleFactor(kIndex) * rmsvalue(kIndex)
1558 end do
1559
1560 ! write to file
1561 if (mmpi_myid == 0) then
1562 write(*,*) 'epp_printRmsStats: Opening ascii output file: ', trim(fileName)
1563 nulFile = 0
1564 ierr = fnom (nulFile, fileName, 'SEQ+R/W', 0)
1565 if (ierr /= 0) then
1566 call utl_abort('epp_printRmsStats: Cannot open ascii output file')
1567 end if
1568
1569 kIndexCount = 0
1570 if ( (nomvar_v(1) == 'UU') .and. (nomvar_v(1+nLev_M) == 'VV') ) then
1571 do levIndex = 1, nLev_M
1572 kIndexCount = kIndexCount + 2
1573 kIndexUU = levIndex
1574 kIndexVV = levIndex + nLev_M
1575 write(nulFile,100) &
1576 elapsed,ftype,nEns,nomvar_v(kIndexUU),pressureOrDepth(kIndexUU),rmsvalue(kIndexUU)
1577 write(nulFile,100) &
1578 elapsed,ftype,nEns,nomvar_v(kIndexVV),pressureOrDepth(kIndexVV),rmsvalue(kIndexVV)
1579 end do
1580 end if
1581 do kIndex = kIndexCount+1, numK
1582 write(nulFile,100) &
1583 elapsed,ftype,nEns,nomvar_v(kIndex),pressureOrDepth(kIndex),rmsvalue(kIndex)
1584 end do
1585 ierr = fclos(nulFile)
1586 end if
1587
1588100 format(f7.2,1x,A1,1x,I5,1x,A4,1x,f12.7,1x,(2E12.5))
1589
1590 deallocate(nomvar_v)
1591 deallocate(pressureOrDepth)
1592 deallocate(scaleFactor)
1593 deallocate(rmsvalue)
1594
1595 end subroutine epp_printRmsStats
1596
1597end module ensPostProcess_mod
