1program midas_randomPert
2 !
3 !:Purpose: Main program for generating an ensemble of random perturbations,
4 ! based on the B matrix (can be homogeneous and isotropic, ensemble-based
5 ! or a weighted sum of covariance matrices).
6 !
7 ! ---
8 !
9 !:Algorithm: Uses a gaussian pseudorandom generator to produce an ensemble
10 ! of control vectors. These control vectors are transformed into
11 ! pertubations in physical space, through the action of the square-root
12 ! of the background-error covariance matrix. The ensemble of perturbations
13 ! can then be debiased and some smoothing can be applied.
14 ! The perturbations thus obtained can then be interpolated and added
15 ! to a provided ensemble mean and/or be mixed with some previous time
16 ! perturbations.
17 !
18 ! ---
19 !
20 !:File I/O:
21 !
22 !============================================== ==============================================================
23 ! Input and Output Files Description of file
24 !============================================== ==============================================================
25 ! ``analysisgrid`` In - File defining grid on which perturbations are generated
26 ! ``targetgrid`` In - File defining destination ensemble grid
27 ! ``bgcov`` In - Static (i.e. NMC) B matrix file for NWP fields
28 ! ``seaIceAnalysis`` In - Sea-ice analysis, to set under-ice perturbations to zero
29 ! ``pert_$YYYYMMDDHH_000_$NNNN`` In - Previous date perturbations, when ensemble not provided
30 ! ``$YYYYMMDDHH_006_$NNNN`` In/Out - Ensemble member files
31 ! ``pert_$YYYYMMDDHH_000_$NNNN`` Out - Output perturbations
32 !============================================== ==============================================================
33 !
34 !
35 !:Synopsis: Below is a summary of the ``randomPert`` program calling sequence:
36 !
37 ! - **Initial setups:**
38 !
39 ! - Initialize the horizontal and vertical grid objects; for the
40 ! generation of the perturbations from ``analysisgrid`` and
41 ! for the ensemble output from ``targetgrid``.
42 !
43 ! - Setup the B matrices: ``bmat_setup``.
44 !
45 ! - Allocate arrays and a stateVector object on the ``analysisgrid``,
46 ! to receive the perturbations and first statistical moments.
47 !
48 ! - Initialize the pseudorandom generator : ``rng_setup``.
49 !
50 ! - **Generation of the ensemble of perturbations**.
51 ! Iterate until the number of required perturbations (``NENS``)
52 ! is generated:
53 !
54 ! - Generate a random control vector from a unit variance gaussian
55 ! distribution: ``rng_gaussian``.
56 !
57 ! - Transform control vector to physical space perturbation using
58 ! ``bmat_sqrtB``.
59 !
60 ! - **Unbias the ensemble of perturbations** if ``REMOVE_MEAN`` is true.
61 !
62 ! - **Smooth variances** if ``smoothVariances`` is true.
63 !
64 ! - **Read ensemble mean** if ``readEnsMean`` is true:
65 !
66 ! - Allocate a stateVector object on the ``targetgrid``, and then
67 ! read the ensemble mean using ``gio_readFromFile``.
68 !
69 ! - if ``setPertZeroUnderIce`` is true, allocate a stateVector
70 ! object and read sea-ice analysis using ``gio_readFromFile``.
71 ! Where the ice fraction is higher than ``iceFractionThreshold``,
72 ! the perturbations are set to zero.
73 !
74 ! - **Read and debias previous date perturbations** if
75 ! ``previousDateFraction`` is higher than zero.
76 !
77 ! - **Final steps**:
78 !
79 ! - Allocate a statevector object on the ``targetgrid``, copy ocean
80 ! mask and interpolate perturbations using ``int_interp_gsv``
81 !
82 ! - If ``previousDateFraction`` higher than zero, average previous
83 ! date perturbations with the ones generated for the present date.
84 !
85 ! - If ``readEnsMean`` is true, add the ensemble mean to the perturbations.
86 !
87 ! - Write the perturbations using ``gio_writeToFile``. If ``readEnsMean``
88 ! is true, write the ensemble mean as well.
89 !
90 ! ---
91 !
92 !:Options: `List of namelist blocks <../namelists_in_each_program.html#randompert>`_
93 ! that can affect the ``randomPert`` program.
94 !
95 ! * ``&NAMENKF`` define specific parameters to configure the ``randomPert``
96 ! program and is read by the program.
97 !
98 ! * The background-error covariance matrix contributions are defined through
99 ! the namelists ``&NAMBHI``, for the homogeneous and isotropic part,
100 ! and ``&NAMBEN`` for the ensemble part. When the matrix is a weighted
101 ! sum, it is important to make sure to define the ``SCALEFACTOR`` values.
102 ! Note that it would be also possible to define contributions from
103 ! ``&NAMBDIFF`` for diffusion correlations and ``&NAMBCHM`` for homogeneous
104 ! and isotropic chemical constituents covariances.
105 !
106
107 use version_mod
108 use midasMpi_mod
109 use ramDisk_mod
110 use controlVector_mod
111 use gridStateVector_mod
112 use gridStateVectorFileIO_mod
113 use bMatrix_mod
114 use interpolation_mod
115 use verticalCoord_mod
116 use horizontalCoord_mod
117 use timeCoord_mod
118 use randomNumber_mod
119 use utilities_mod
120 use gridVariableTransforms_mod
121 implicit none
122
123 type(struct_gsv) :: stateVectorPert, stateVectorPertInterp
124 type(struct_gsv) :: stateVectorEnsMean, stateVectorIce
125
126 type(struct_vco), pointer :: vco_anl => null()
127 type(struct_hco), pointer :: hco_anl => null()
128 type(struct_hco), pointer :: hco_anlcore => null()
129 type(struct_hco), pointer :: hco_target => null()
130 type(struct_hco), pointer :: hco_targetcore => null()
131
132 real(4), pointer :: seaice_ptr(:,:,:)
133 real(8), pointer :: field(:,:,:), fieldInterp(:,:,:)
134
135 integer :: fclos, fnom, fstopc, newdate, dateStamp, datePrevious, dateStampPrevious
136 integer :: imode, ierr
137 integer :: memberIndex, lonIndex, latIndex, cvIndex, levIndex, nkgdim
138 integer :: datePrint, timePrint, nulnam, randomSeed
139 integer :: get_max_rss, n_grid_point, n_grid_point_glb
140
141 integer :: latPerPEa, latPerPEmaxa, myLatBega, myLatEnda
142 integer :: lonPerPEa, lonPerPEmaxa, myLonBega, myLonEnda
143 integer :: latPerPEt, latPerPEmaxt, myLatBegt, myLatEndt
144 integer :: lonPerPEt, lonPerPEmaxt, myLonBegt, myLonEndt
145
146 real(8) :: previousDateFactPrev, previousDateFactNoise
147 real(8), allocatable :: controlVector(:), controlVector_mpiglobal(:)
148 real(8), allocatable :: gdmean(:,:,:)
149 real(4), allocatable :: ensemble_r4(:,:,:,:)
150 real(4), allocatable :: ensemblePreviousDate_r4(:,:,:,:)
151 real(8), allocatable :: avg_pturb_var(:)
152 real(8), allocatable :: avg_pturb_var_glb(:)
153 real(8), allocatable :: pturb_var(:,:,:)
154
155 logical :: targetGridExists
156 character(len=10) :: dateString, datePreviousString
157 character(len=4) :: memberString
158 character(len=25) :: outFileName, inFileName
159 character(len=64) :: ensMeanFileName = 'ensMeanState'
160 character(len=2) :: typvarOut
161
162 ! Namelist variables
163 logical :: remove_mean ! choose to remove mean from perturbations
164 logical :: smoothVariances ! choose to impose horizontally constant perturbation variances
165 logical :: mpiTopoIndependent ! choose to compute random numbers with mpi-topology-independent method
166 logical :: readEnsMean ! choose to read ens mean and add this to the perturbations
167 logical :: setPertZeroUnderIce ! choose to set perturbation to zero under sea ice (for SST)
168 integer :: nens ! number of perturbations to compute
169 integer :: seed ! initial value of the random seed
170 integer :: numBits ! number of bits to use when writing to standard files
171 character(len=12) :: out_etiket ! the 'etiket' to write to standard files
172 real(4) :: iceFractionThreshold ! ice fraction threshold to use in combination with 'setPertZeroUnderIce'
173 real(4) :: previousDateFraction ! relative amount of previous date perturbations to include in current perturbations
174 NAMELIST /NAMENKF/nens, seed, out_etiket, remove_mean, &
175 smoothVariances, mpiTopoIndependent, numBits, &
176 readEnsMean, setPertZeroUnderIce, iceFractionThreshold, &
177 previousDateFraction
178
179 call ver_printNameAndVersion('randomPert','Generation of random perturbations')
180
181 !
182 !- 0. MPI, tmg initialization
183 !
184 call mmpi_initialize
185 call tmg_init(mmpi_myid, 'TMG_INFO')
186
187 call utl_tmg_start(0,'Main')
188 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
189
190 ierr = fstopc('MSGLVL','ERRORS',0)
191
192 !
193 !- 1. Set/Read values for the namelist NAMENKF
194 !
195
196 !- 1.1 Setting default values
197 nens = 10
198 seed = -999 ! If -999, set random seed using the date
199 remove_mean = .true.
200 out_etiket='RANDOM_PERT'
201 smoothVariances = .false.
202 mpiTopoIndependent = .false.
203 numBits = 32
204 readEnsMean = .false.
205 setPertZeroUnderIce = .false.
206 iceFractionThreshold = 0.2
207 previousDateFraction = -1.0
208
209 !- 1.2 Read the namelist
210 nulnam=0
211 ierr=fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
212 read(nulnam, nml=namenkf, iostat=ierr)
213 if(ierr.ne.0) call utl_abort('midas-randomPert: Error reading namelist')
214 if( mmpi_myid == 0 ) write(*,nml=namenkf)
215 ierr=fclos(nulnam)
216
217 if (readEnsMean) then
218 typvarOut = 'A'
219 else
220 typvarOut = 'R'
221 end if
222
223 if (previousDateFraction > 1.0) then
224 call utl_abort('midas-randomPert: previousDateFraction must be less than 1.0 to give stable results')
225 end if
226
227 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
228
229 !
230 !- 2. Initialization
231 !
232
233 ! Setup the ramdisk directory (if supplied)
234 call ram_setup
235
236 !- 2.1 Set the dateStamp, either from env variable or ensMeanState file
237
238 !- Initialize the Temporal grid and the dateStamp from env variable
239 call tim_setup()
240
241 ! If dateStamp not set by env variable, use date from ensMeanState, if available
242 if (tim_getDateStamp() == 0) then
243 if (readEnsMean) then
244 dateStamp = tim_getDatestampFromFile(ensMeanFileName)
245 call tim_setDatestamp(dateStamp)
246 else
247 call utl_abort('midas-randomPert: DateStamp must be set through env variable')
248 end if
249 end if
250 dateStamp = tim_getDateStamp()
251 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
252 ierr = newdate(dateStamp, datePrint, timePrint, imode)
253 write(dateString, '(I10)') datePrint*100 + timePrint/1000000
254 if( mmpi_myid == 0 ) then
255 write(*,*) ' date= ', datePrint, ' time= ', timePrint, ' stamp= ', dateStamp
256 write(*,*) ' dateString = ', dateString
257 end if
258
259 !- 2.2 Initialize variables of the model states
260 call gsv_setup
261
262 !- 2.3 Initialize the horizontal grids from analysisgrid and targetgrid files
263 if (mmpi_myid == 0) write(*,*)
264 if (mmpi_myid == 0) write(*,*) 'Set hco parameters for analysis grid'
265 call hco_setupFromFile(hco_anl, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
266
267 if ( hco_anl % global ) then
268 hco_anlcore => hco_anl
269 else
270 !- Iniatilized the core (Non-Exteded) analysis grid
271 call hco_setupFromFile( hco_anlcore, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
272 end if
273
274 inquire(file='./targetgrid',exist=targetGridExists)
275 if (targetGridExists) then
276 if (mmpi_myid == 0) write(*,*)
277 if (mmpi_myid == 0) write(*,*) 'Set hco parameters for target grid'
278
279 call hco_setupFromFile(hco_target, './targetgrid', 'ANALYSIS', 'Target' ) ! IN
280
281 if ( hco_target % global ) then
282 hco_targetcore => hco_target
283 else
284 !- Iniatilized the core (Non-Exteded) analysis grid
285 call hco_setupFromFile( hco_targetcore, './targetgrid', 'COREGRID', 'TargetCore' ) ! IN
286 end if
287 else
288 ! no targetgrid file, so use analysisgrid grid instead
289 hco_target => hco_anl
290 hco_targetcore => hco_anlcore
291 end if
292
293 call mmpi_setup_latbands(hco_anl % nj, & ! IN
294 latPerPEa, latPerPEmaxa, myLatBega, myLatEnda ) ! OUT
295 call mmpi_setup_lonbands(hco_anl % ni, & ! IN
296 lonPerPEa, lonPerPEmaxa, myLonBega, myLonEnda ) ! OUT
297
298 call mmpi_setup_latbands(hco_target % nj, & ! IN
299 latPerPEt, latPerPEmaxt, myLatBegt, myLatEndt ) ! OUT
300 call mmpi_setup_lonbands(hco_target % ni, & ! IN
301 lonPerPEt, lonPerPEmaxt, myLonBegt, myLonEndt ) ! OUT
302
303 !- 2.4 Initialize the vertical coordinate from the analysisgrid file
304 call vco_setupFromFile(vco_anl, './analysisgrid', ' ')
305
306 !- 2.5 Initialize the B_hi matrix
307 call bmat_setup(hco_anl, hco_anlcore, vco_anl)
308
309 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
310
311 !- 2.6 Initialize the gridded variable transform module
312 call gvt_setup(hco_anl,hco_anlcore,vco_anl)
313 if ( gsv_varExist(varName='HU') ) call gvt_setupRefFromTrialFiles('HU')
314
315 !- 2.7 Set randomSeed
316 if (seed == -999) then
317 ! If "seed" namelist value is -999, set random seed using the date
318 dateStamp = tim_getDateStamp()
319 if (dateStamp == -1) then
320 call utl_abort('midas-randomPert: dateStamp is not set, cannot be used to set random seed')
321 end if
322 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
323 ierr = newdate(dateStamp, datePrint, timePrint, imode)
324 timePrint = timePrint/1000000
325 datePrint = datePrint*100 + timePrint
326 ! Remove the century, keeping 2 digits of the year
327 randomSeed = datePrint - 100000000*(datePrint/100000000)
328 else
329 ! Otherwise, use value from namelist
330 randomSeed = seed
331 end if
332 write(*,*) 'midas-randomPert: randomSeed for set to ', randomSeed
333
334 !
335 !- 3. Memory allocations
336 !
337
338 !- 3.1 Allocate the stateVectorPert
339 call gsv_allocate(stateVectorPert, 1, hco_anl, vco_anl, &
340 dateStamp_opt=dateStamp, mpi_local_opt=.true., &
341 allocHeight_opt=.false., allocPressure_opt=.false., &
342 hInterpolateDegree_opt='LINEAR')
343 nkgdim = stateVectorPert%nk
344 allocate(ensemble_r4(myLonBega:myLonEnda, myLatBega:myLatEnda, nkgdim, nEns))
345
346 !- 3.2 Allocate auxillary variables
347 allocate(gdmean(myLonBega:myLonEnda, myLatBega:myLatEnda, nkgdim))
348 allocate(controlVector(cvm_nvadim))
349 allocate(controlVector_mpiglobal(cvm_nvadim_mpiglobal))
350
351 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
352
353 !
354 !- 4. Compute an ensemble of random perturbations
355 !
356
357 if( mmpi_myid == 0 ) write(*,*) '******************'
358 if( mmpi_myid == 0 ) write(*,*) 'midas-randomPert: COMPUTE the mean of the random ', &
359 'perturbations of all the members'
360
361 if( mpiTopoIndependent ) then
362 call rng_setup(abs(randomSeed))
363 else
364 call rng_setup(abs(randomSeed+mmpi_myid))
365 end if
366
367 gdmean(:,:,:) = 0.0D0
368 call gsv_getField(stateVectorPert,field)
369
370 !- 4.1 Generate a (potentially) biased ensemble
371 do memberIndex = 1, NENS
372
373 if( mmpi_myid == 0 ) write(*,*) '...computing member number= ', memberIndex
374
375 !- 4.1.1 Create a random control vector in spectral space
376
377 if( mpiTopoIndependent ) then
378 !- Global vector (for testing different mpi topology, less efficient)
379 do cvIndex = 1, cvm_nvadim_mpiglobal
380 controlVector_mpiglobal(cvIndex) = rng_gaussian()
381 end do
382 call bmat_reduceToMPILocal( controlVector, controlVector_mpiglobal )
383 else
384 !- Local vector (different randomSeed for each processor, more efficient)
385 do cvIndex = 1, cvm_nvadim
386 controlVector(cvIndex) = rng_gaussian()
387 end do
388 end if
389
390 !- 4.1.2 Transform to control variables in physical space
391 call bmat_sqrtB(controlVector, cvm_nvadim, & ! IN
392 stateVectorPert ) ! OUT
393
394 !- 4.1.3 Copy perturbations to big array and update ensemble sum
395 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
396 do levIndex = 1, nkgdim
397 do latIndex = myLatBega, myLatEnda
398 do lonIndex = myLonBega, myLonEnda
399 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) &
400 = real(field(lonIndex, latIndex, levIndex), 4)
401 gdmean(lonIndex, latIndex, levIndex) &
402 = gdmean(lonIndex, latIndex, levIndex) + field(lonIndex, latIndex, levIndex)
403 end do
404 end do
405 end do
406 !$OMP END PARALLEL DO
407
408 end do
409
410 call gsv_deallocate(stateVectorPert)
411
412 !- 4.2 Remove the ensemble mean
413 if ( REMOVE_MEAN ) then
414
415 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
416 do levIndex = 1, nkgdim
417 do latIndex = myLatBega, myLatEnda
418 do lonIndex = myLonBega, myLonEnda
419 gdmean(lonIndex, latIndex, levIndex) &
420 = gdmean(lonIndex, latIndex, levIndex) / real(NENS, 8)
421 end do
422 end do
423 end do
424 !$OMP END PARALLEL DO
425
426 !$OMP PARALLEL DO PRIVATE (memberIndex, levIndex, latIndex, lonIndex)
427 do memberIndex = 1, NENS
428 do levIndex = 1, nkgdim
429 do latIndex = myLatBega, myLatEnda
430 do lonIndex = myLonBega, myLonEnda
431 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) = &
432 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) - &
433 real(gdmean(lonIndex, latIndex, levIndex), 4)
434 end do
435 end do
436 end do
437 end do
438 !$OMP END PARALLEL DO
439
440 end if
441
442 !- 4.3 Smooth variances to horizontally constant values
443 if ( smoothVariances ) then
444
445 allocate(pturb_var(myLonBega:myLonEnda, myLatBega:myLatEnda, stateVectorPert%nk))
446 allocate(avg_pturb_var(stateVectorPert%nk), avg_pturb_var_glb(stateVectorPert%nk))
447 pturb_var(:,:,:) = 0.0D0
448 avg_pturb_var(:) = 0.0D0
449 avg_pturb_var_glb(:) = 0.0D0
450
451 do memberIndex = 1, NENS
452 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
453 do levIndex = 1, nkgdim
454 do latIndex = myLatBega, myLatEnda
455 do lonIndex = myLonBega, myLonEnda
456 pturb_var(lonIndex, latIndex, levIndex) = pturb_var(lonIndex, latIndex, levIndex) &
457 + real(ensemble_r4(lonIndex, latIndex, levIndex, memberIndex), 8)**2
458 end do
459 end do
460 end do
461 !$OMP END PARALLEL DO
462 end do
463
464 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
465 do levIndex = 1, nkgdim
466 do latIndex = myLatBega, myLatEnda
467 do lonIndex = myLonBega, myLonEnda
468 pturb_var(lonIndex, latIndex, levIndex) &
469 = pturb_var(lonIndex, latIndex, levIndex) / real(NENS, 8)
470 end do
471 end do
472 end do
473 !$OMP END PARALLEL DO
474
475 do latIndex = myLatBega, myLatEnda
476 do lonIndex = myLonBega, myLonEnda
477 do levIndex = 1, nkgdim
478 avg_pturb_var(levIndex) = avg_pturb_var(levIndex) + pturb_var(lonIndex, latIndex, levIndex)
479 end do
480 end do
481 end do
482
483 n_grid_point=(lonPerPEa)*(latPerPEa)
484 call rpn_comm_allreduce(n_grid_point, n_grid_point_glb, 1, &
485 "mpi_double_precision", "mpi_sum", "GRID", ierr)
486 call rpn_comm_allreduce(avg_pturb_var, avg_pturb_var_glb, nkgdim, &
487 "mpi_double_precision", "mpi_sum", "GRID", ierr)
488
489 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
490 do levIndex = 1, nkgdim
491 do latIndex = myLatBega, myLatEnda
492 do lonIndex = myLonBega, myLonEnda
493 if( pturb_var(lonIndex, latIndex, levIndex) > 0.0d0 ) then
494 pturb_var(lonIndex, latIndex, levIndex) = sqrt( pturb_var(lonIndex, latIndex, levIndex))
495 end if
496 end do
497 end do
498 end do
499 !$OMP END PARALLEL DO
500
501 do levIndex = 1, nkgdim
502 if( avg_pturb_var_glb(levIndex) > 0.0d0 ) then
503 avg_pturb_var_glb(levIndex) = sqrt( avg_pturb_var_glb(levIndex) / real(n_grid_point_glb, 8) )
504 end if
505 end do
506
507 !$OMP PARALLEL DO PRIVATE (memberIndex, levIndex, latIndex, lonIndex)
508 do memberIndex = 1, NENS
509 do levIndex = 1, nkgdim
510 if( avg_pturb_var_glb(levIndex) > 0.0d0 ) then
511 do latIndex = myLatBega, myLatEnda
512 do lonIndex = myLonBega, myLonEnda
513 if( pturb_var(lonIndex, latIndex, levIndex) > 0.0d0 ) then
514 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) = &
515 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) * &
516 real(avg_pturb_var_glb(levIndex)/pturb_var(lonIndex, latIndex, levIndex), 4)
517 end if
518 end do
519 end do
520 end if
521 end do
522 end do
523 !$OMP END PARALLEL DO
524
525 end if
526
527 !- 4.4 Read ensemble mean state
528 if ( readEnsMean ) then
529 if (mmpi_myid == 0) write(*,*) 'midas-randomPert: reading ensemble mean state'
530
531 call gsv_allocate(stateVectorEnsMean, 1, hco_target, vco_anl, &
532 dateStamp_opt=dateStamp, mpi_local_opt=.true., &
533 allocHeight_opt=.false., allocPressure_opt=.false., &
534 hInterpolateDegree_opt='LINEAR')
535 call gio_readFromFile(stateVectorEnsMean, ensMeanFileName, ' ', ' ', &
536 containsFullField_opt=.true.)
537
538 end if
539
540 !- 4.5 Set perturbations to zero under ice
541 if ( setPertZeroUnderIce ) then
542
543 ! read sea-ice analysis
544 call gsv_allocate(stateVectorIce, 1, hco_anl, vco_anl, dataKind_opt = 4, &
545 datestamp_opt = -1, mpi_local_opt = .true., &
546 varNames_opt = (/'LG'/), hInterpolateDegree_opt ='LINEAR')
547 call gio_readFromFile(stateVectorIce, './seaIceAnalysis', ' ','A', &
548 unitConversion_opt=.false., containsFullField_opt=.true.)
549 call gsv_getField(stateVectorIce,seaice_ptr)
550
551 ! set perturbations to zero
552 !$OMP PARALLEL DO PRIVATE (memberIndex, levIndex, latIndex, lonIndex)
553 do memberIndex = 1, NENS
554 do levIndex = 1, nkgdim
555 do latIndex = myLatBega, myLatEnda
556 do lonIndex = myLonBega, myLonEnda
557 if (seaice_ptr(lonIndex, latIndex, 1) >= iceFractionThreshold) then
558 ensemble_r4(lonIndex, latIndex, levIndex, memberIndex) = 0.0
559 end if
560 end do
561 end do
562 end do
563 end do
564 !$OMP END PARALLEL DO
565
566 call gsv_deallocate(stateVectorIce)
567
568 end if
569
570 !- 4.6 Add a fraction of previous date perturbations
571 if (previousDateFraction > 0.0) then
572
573 previousDateFactPrev = previousDateFraction
574 ! reduce the amount of random perturbation to maintain steady-state variance
575 previousDateFactNoise = sqrt(1.0d0 - previousDateFraction**2)
576 if (mmpi_myid == 0) then
577 write(*,*) 'midas-randomPert: Scale factor applied to previous date = ', &
578 previousDateFactPrev
579 write(*,*) 'midas-randomPert: Scale factor applied to random perturbation = ', &
580 previousDateFactNoise
581 end if
582
583 ! determine dateStamp of previous date
584 call incdatr(dateStampPrevious, dateStamp, -tim_windowsize)
585 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
586 ierr = newdate(dateStampPrevious, datePrint, timePrint, imode)
587 datePrevious = datePrint*100 + timePrint/1000000
588 write(datePreviousString, '(I10)') datePrevious
589 write(*,*) 'midas-randomPert: previous date, stamp = ', datePrevious, dateStampPrevious
590
591 call gsv_allocate(stateVectorPert, 1, hco_target, vco_anl, &
592 dateStamp_opt=dateStampPrevious, mpi_local_opt=.true., &
593 allocHeight_opt=.false., allocPressure_opt=.false., &
594 hInterpolateDegree_opt='LINEAR')
595 call gsv_getField(stateVectorPert,field)
596
597 deallocate(gdmean)
598 allocate(gdmean(myLonBegt:myLonEndt, myLatBegt:myLatEndt, nkgdim))
599 gdmean(:,:,:) = 0.0D0
600 allocate(ensemblePreviousDate_r4(myLonBegt:myLonEndt, myLatBegt:myLatEndt, nkgdim, nEns))
601
602 do memberIndex = 1, NENS
603
604 ! Read previous date perturbations (or perturbed analyses)
605 write(memberString, '(I4.4)') memberIndex
606 if (readEnsMean) then
607 inFileName = './'//trim(datePreviousString)//'_000_'//trim(memberString)
608 else
609 inFileName = './pert_'//trim(datePreviousString)//'_'//trim(memberString)
610 end if
611 if( mmpi_myid == 0 ) write(*,*) 'midas-randomPert: reading previous date file = ', inFileName
612
613 call gio_readFromFile(stateVectorPert, inFileName, ' ', ' ', &
614 containsFullField_opt=.true.)
615
616 ! Copy to big array and accumulate sum for computing mean
617 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
618 do levIndex = 1, nkgdim
619 do latIndex = myLatBegt, myLatEndt
620 do lonIndex = myLonBegt, myLonEndt
621 ensemblePreviousDate_r4(lonIndex, latIndex, levIndex, memberIndex) = &
622 real(field(lonIndex, latIndex, levIndex), 4)
623 gdmean(lonIndex, latIndex, levIndex) = gdmean(lonIndex, latIndex, levIndex) &
624 + field(lonIndex, latIndex, levIndex)
625 end do
626 end do
627 end do
628 !$OMP END PARALLEL DO
629
630 end do
631
632 ! Finish computing mean and remove it
633 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
634 do levIndex = 1, nkgdim
635 do latIndex = myLatBegt, myLatEndt
636 do lonIndex = myLonBegt, myLonEndt
637 gdmean(lonIndex, latIndex, levIndex) = gdmean(lonIndex, latIndex, levIndex) &
638 / real(NENS, 8)
639 end do
640 end do
641 end do
642 !$OMP END PARALLEL DO
643
644 !$OMP PARALLEL DO PRIVATE (memberIndex, levIndex, latIndex, lonIndex)
645 do memberIndex = 1, NENS
646 do levIndex = 1, nkgdim
647 do latIndex = myLatBegt, myLatEndt
648 do lonIndex = myLonBegt, myLonEndt
649 ensemblePreviousDate_r4(lonIndex, latIndex, levIndex, memberIndex) = &
650 ensemblePreviousDate_r4(lonIndex, latIndex, levIndex, memberIndex) - &
651 real(gdmean(lonIndex, latIndex, levIndex), 4)
652 end do
653 end do
654 end do
655 end do
656 !$OMP END PARALLEL DO
657
658 call gsv_deallocate(stateVectorPert)
659
660 end if
661
662 !- 4.7 Write the perturbations
663 do memberIndex = 1, NENS
664 if( mmpi_myid == 0 ) write(*,*)
665 if( mmpi_myid == 0 ) write(*,*) 'midas-randomPert: pre-processing for writing member number= ', memberIndex
666
667 call gsv_allocate(stateVectorPert, 1, hco_anl, vco_anl, &
668 dateStamp_opt=dateStamp, mpi_local_opt=.true., &
669 allocHeight_opt=.false., allocPressure_opt=.false., &
670 hInterpolateDegree_opt='LINEAR')
671 call gsv_getField(stateVectorPert,field)
672 call gsv_allocate(stateVectorPertInterp, 1, hco_target, vco_anl, &
673 dateStamp_opt=dateStamp, mpi_local_opt=.true., &
674 allocHeight_opt=.false., allocPressure_opt=.false., &
675 hInterpolateDegree_opt='LINEAR')
676
677 ! Copy mask if it exists
678 call gsv_copyMask(stateVectorEnsMean, stateVectorPertInterp)
679
680 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
681 do levIndex = 1, nkgdim
682 do latIndex = myLatBega, myLatEnda
683 do lonIndex = myLonBega, myLonEnda
684 field(lonIndex, latIndex, levIndex) = ensemble_r4(lonIndex, latIndex, levIndex, memberIndex)
685 end do
686 end do
687 end do
688 !$OMP END PARALLEL DO
689
690 ! interpolate perturbation to the target grid
691 call int_interp_gsv(stateVectorPert, stateVectorPertInterp)
692
693 call gsv_getField(stateVectorPertInterp,fieldInterp)
694
695 ! Average current and previous perturbations
696 if (previousDateFraction > 0.0) then
697 !$OMP PARALLEL DO PRIVATE (levIndex, latIndex, lonIndex)
698 do levIndex = 1, nkgdim
699 do latIndex = myLatBegt, myLatEndt
700 do lonIndex = myLonBegt, myLonEndt
701 fieldInterp(lonIndex, latIndex, levIndex) = &
702 previousDateFactPrev * &
703 ensemblePreviousDate_r4(lonIndex, latIndex, levIndex, memberIndex) + &
704 previousDateFactNoise * &
705 fieldInterp(lonIndex, latIndex, levIndex)
706 end do
707 end do
708 end do
709 !$OMP END PARALLEL DO
710 end if
711
712 ! add perturbation to supplied ensemble mean
713 if ( readEnsMean ) then
714 if (mmpi_myid == 0) write(*,*) 'midas-randomPert: adding the ensemble mean and perturbation'
715 call gsv_add(stateVectorEnsMean, stateVectorPertInterp)
716 end if
717
718 ! determine file name and write to file
719 write(memberString, '(I4.4)') memberIndex
720 if (readEnsMean) then
721 outFileName = './'//trim(dateString)//'_000_'//trim(memberString)
722 else
723 outFileName = './pert_'//trim(dateString)//'_'//trim(memberString)
724 end if
725
726 if( mmpi_myid == 0 ) write(*,*) 'midas-randomPert: processing file = ', outFileName
727 stateVectorPertInterp%etiket = 'UNDEFINED'
728 call gio_writeToFile(stateVectorPertInterp, outFileName, out_etiket, & ! IN
729 numBits_opt=numBits, unitConversion_opt=.true., & ! IN
730 containsFullField_opt=readEnsMean, typvar_opt=typvarOut)
731
732 call gsv_deallocate(stateVectorPertInterp)
733 call gsv_deallocate(stateVectorPert)
734
735 end do
736
737 ! Write out ensemble mean state, if it exists
738 if (readEnsMean) then
739 ! determine file name and write to file
740 memberIndex = 0
741 write(memberString, '(I4.4)') memberIndex
742 outFileName = './'//trim(dateString)//'_000_'//trim(memberString)
743
744 if( mmpi_myid == 0 ) write(*,*) 'midas-randomPert: processing file = ', outFileName
745 stateVectorEnsMean%etiket = 'UNDEFINED'
746 call gio_writeToFile(stateVectorEnsMean, outFileName, out_etiket, & ! IN
747 numBits_opt=numBits, unitConversion_opt=.true., & ! IN
748 containsFullField_opt=.true., typvar_opt=typvarOut)
749 end if
750
751 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
752
753 !
754 !- 5. Memory deallocations
755 !
756 deallocate(gdmean)
757 deallocate(ensemble_r4)
758 if (previousDateFraction > 0.0) then
759 deallocate(ensemblePreviousDate_r4)
760 end if
761 deallocate(controlVector)
762 deallocate(controlVector_mpiglobal)
763
764 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
765 call utl_tmg_stop(0)
766
767 !
768 !- 6. MPI, tmg finalize
769 !
770 call tmg_terminate(mmpi_myid, 'TMG_INFO')
771 call rpn_comm_finalize(ierr)
772
773 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
774
775 !
776 !- 7. Ending
777 !
778 if( mmpi_myid == 0 ) write(*,*) '---------------------------------'
779 if( mmpi_myid == 0 ) write(*,*) ' MIDAS-RANDOMPERT ENDS'
780 if( mmpi_myid == 0 ) write(*,*) '---------------------------------'
781
782end program midas_randomPert
