1program midas_diagBmatrix
2 !
3 !:Purpose: Program for computing diagnostics of the background-error covariances (**B**) and
4 ! localization (**L**) matrices
5 !
6 ! ---
7 !
8 !:Algorithm: The ``diagBmatrix`` program performs two different types of diagnostics
9 ! based on the namelist options.
10 !
11 ! - **Columns of B and L**: Allow to extract columns of **B** (the background-error covariances matrix)
12 ! which is a useful way to simulate single-observation (a.k.a pseudo-obs) experiments without using
13 ! real observations.
14 ! See e.g. Figs 5 and 6 in <https://doi.org/10.1175/MWR-D-18-0248.1>.
15 ! If localization is used, the corresponding columns of the localization matrix (**L**) will be outputed.
16 ! - **Implied variances**: Compute the variances resulting from an ensemble of random perturbations
17 ! generated from the provided **B** matrix, including hybrid covariances. This allow to evaluate the impact
18 ! of various choices in the design of **B** on the modification to the variances of the ensemble of
19 ! background-error estimates used to derivate **B**, i.e. the "measured" variances.
20 !
21 ! ---
22 !
23 !============================================= ==============================================================
24 ! Input and Output Files Description of file
25 !============================================= ==============================================================
26 ! ``flnml`` In - Main namelist file with parameters user may modify
27 ! ``ensemble/$YYYYMMDDHH_$HHH_$NNNN`` In - Background ensemble member files
28 ! ``bgcov`` In - **B** nmc matrix for the random perturbations
29 ! ``analysis_grid`` In - Horizontal grid file on which the random perturbations
30 ! are generated
31 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
32 ! Only needed if a LQ-HU transform is used in **B**
33 ! ``columnB_$VAR_$YYYYMMDDHH.fst`` Out - Column of **B**
34 ! ``columnBnorm_$VAR_$YYYYMMDDHH.fst`` Out - Column of **B** normalize by the value at the pseudo-obs location
35 ! ``columnL_$YYYYMMDDHH.fst`` Out - Column of **L** (only when ensemble-derived **B** is use)
36 ! ``stddev_$YYYYMMDDHH.fst`` Out - Implied variances
37 !============================================= ==============================================================
38 !
39 ! --
40 !
41 !:Synopsis: Below is a summary of the ``diagBmatrix`` program calling sequence:
42 !
43 ! - **Initial setups:**
44 !
45 ! - Read the NAMDIAG namelist and check/modify some values.
46 !
47 ! - Setup horizontal and vertical grid objects from the provided template file (./analysisgrid)
48 !
49 ! - Various modules are setup: ``gridStateVector_mod``, ``timeCoord_mod`` and ``bmatrix_mod``
50 !
51 ! - **Columns of B and L:**
52 !
53 ! - Attribute time bin of the pseudo-obs based on the namelist option.
54 !
55 ! - For each variables and prescribed pseudo-obs positions:
56 ! 1) create a Dirac delta distribution (zero everywhere except at obs location).
57 ! This is equivalent to impose (O-F / sigma_obs) = 1 at one location.
58 ! 2) apply **B** ^(1/2)^T **B** ^1/2 and output the results into a file
59 ! 3) normalize the results by the value at the pseudo-obs position
60 ! and output to a file
61 !
62 ! - If bMatrixEnsemble_mod is activated, for each prescribed pseudo-obs positions:
63 ! 1) create a Dirac delta distribution
64 ! 2) apply **L** ^(1/2)^T **L** ^1/2 and output the results into a file
65 !
66 ! - **Implied variances:**
67 !
68 ! - For each member of the chosen random ensemble size:
69 ! 1) create a random control vector
70 ! 2) apply **B** ^1/2
71 ! 3) store the resulting 3D state in randomEns
72 !
73 ! - Compute and remove the ensemble mean
74 !
75 ! - Compute the ensemble standard deviations at each grid point and output to a file
76 !
77 ! - Compute the zonal and domain mean standard deviations and output to the same file
78 !
79 !======================== ============ ==============================================================
80 ! Module Namelist Description of what is controlled
81 !======================== ============ ==============================================================
82 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal resolution of
83 ! the background state and increment
84 ! ``bMatrixEnsemble_mod`` ``NAMBEN`` weight and other parameters for ensemble-based **B** matrix
85 ! component
86 ! ``bMatrixHI_mod`` ``NAMBHI`` weight and other parameters for the climatological **B** matrix
87 ! component based on homogeneous-isotropic covariances
88 ! represented in spectral space
89 ! Other **B** matrix modules various weight and other parameters for each type of **B** matrix
90 !======================== ============ ==============================================================
91 !
92
93 use version_mod
94 use midasMpi_mod
95 use controlVector_mod
96 use gridVariableTransforms_mod
97 use varNameList_mod
98 use gridStateVector_mod
99 use gridStateVectorFileIO_mod
100 use ensembleStateVector_mod
101 use bMatrix_mod
102 use bMatrixEnsemble_mod
103 use localization_mod
104 use horizontalCoord_mod
105 use advection_mod
106 use verticalCoord_mod
107 use timeCoord_mod
108 use randomNumber_mod
109 use utilities_mod
110 use ramDisk_mod
111 implicit none
112
113 type(struct_gsv) :: statevector, statevectorEnsAmplitude
114 type(struct_ens) :: ensAmplitude
115 type(struct_hco), pointer :: hco_anl => null()
116 type(struct_hco), pointer :: hco_core => null()
117 type(struct_vco), pointer :: vco_anl => null()
118 type(struct_loc), pointer :: loc => null()
119 type(struct_adv), pointer :: adv_amplitudeAssimWindow
120
121 real(8), pointer :: field4d(:,:,:,:)
122 real(8), pointer :: field3d(:,:,:)
123 real(8), pointer :: ensAmplitude_oneLev(:,:,:,:)
124 real(4), allocatable :: randomEns(:,:,:,:)
125 real(8), allocatable :: mean(:,:,:)
126 real(8), allocatable :: stddev(:,:,:)
127 real(8), allocatable :: stddev_zm(:,:),stddev_zm2(:,:)
128 real(8), allocatable :: stddev_dm(:,:),stddev_dm2(:,:)
129 real(8), allocatable :: zonalMeanStddev(:)
130 real(8), allocatable :: controlVector(:), controlVector_global(:)
131
132 real(8) :: centralValue, centralValueLocal
133
134 integer :: fclos, fnom, fstopc, newdate, get_max_rss
135 integer :: ierr, nsize, iseed, nulnam, nultxt
136 integer :: ensIndex, index, kIndex, nkgdim, levIndex, lonIndex, latIndex
137 integer :: dateTime, datePrint, timePrint, dateStamp, numLoc, numStepAmplitude
138 integer :: nlevs, nlevs2, varIndex, ip3
139 integer :: locIndex, stepIndexInc, nEns, numBensInstance, instanceIndex
140 integer :: amp3dStepIndex, nLonLatPos, lonLatPosIndex
141 integer :: oneobs_timeStepIndex
142
143 integer, parameter :: lonPosIndex = 1
144 integer, parameter :: latPosIndex = 2
145
146 integer :: latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
147 integer :: myLatBeg, myLatEnd
148 integer :: myLonBeg, myLonEnd
149 integer, allocatable :: dateStampList(:)
150
151 character(len=128) :: filename, filenameInc, filenameIncNorm, filenameEnsAmp
152 character(len=10) :: datestr
153 character(len=4) :: varName
154 character(len=1) :: locIndexString
155 character(len=2) :: instanceIndexString
156
157 character(len=4), parameter :: varNameALFAatm(1) = (/ 'ALFA' /)
158 character(len=4), parameter :: varNameALFAsfc(1) = (/ 'ALFS' /)
159 character(len=4) :: varNameALFA(1)
160
161 ! namelist variables
162 integer :: numperturbations ! number of perturbations for randomization estimate of stddev
163 integer :: nrandseed ! initial random seed value
164 integer :: oneobs_levs(100) ! list of level indexes where B matrix columns are computed
165 integer :: oneobs_lonlat(100,2) ! list of lon,lat index pairs where B matrix columns are computed
166 character(len=128) :: oneobs_timeStep ! can be 'first', 'last' or 'middle'
167 character(len=4) :: oneobs_varName ! can be 'all' or a specific variable name (default='all')
168 logical :: writeEnsAmplitude ! choose to write ensemble amplitude fields (for ensemble B)
169 logical :: writeTextStddev ! choose to write stddev to text file in addition to standard file
170 logical :: writePsiChiStddev ! choose to also write stddev of Psi/Chi in addition to UU/VV
171
172 namelist /namdiag/numperturbations, nrandseed, oneobs_levs, oneobs_lonlat, &
173 oneobs_varName, oneobs_timeStep, writeEnsAmplitude, writeTextStddev, writePsiChiStddev
174
175 call ver_printNameAndVersion('diagBmatrix','Diagnositcs of the B matrix')
176
177 ! MPI, tmg initialization
178 call mmpi_initialize
179 call tmg_init(mmpi_myid, 'TMG_INFO')
180 call utl_tmg_start(0,'Main')
181 ierr = fstopc('MSGLVL','ERRORS',0)
182
183 ! Set default values for namelist NAMDIAG parameters
184 numperturbations = -1
185 nrandseed = 1
186 oneobs_varName = 'all'
187 oneobs_timeStep = 'middle'
188 oneobs_levs(:) = -1
189 oneobs_lonlat(:,:)= -1
190 writeEnsAmplitude = .false.
191 writeTextStddev = .false.
192 writePsiChiStddev = .false.
193
194 ! Read the parameters from NAMDIAG
195 nulnam=0
196 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
197 read(nulnam,nml=namdiag,iostat=ierr)
198 if(ierr.ne.0) call utl_abort('midas-diagBmatrix: Error reading namelist')
199 write(*,nml=namdiag)
200 ierr=fclos(nulnam)
201
202 nlevs=0
203 do index = 1, size(oneobs_levs)
204 if (oneobs_levs(index) >= 1) nlevs=nlevs+1
205 end do
206 nLonLatPos=0
207 do index = 1, size(oneobs_lonlat(:,lonPosIndex))
208 if (oneobs_lonlat(index,lonPosIndex) >= 1 .and. oneobs_lonlat(index,latPosIndex) >= 1) nLonLatPos=nLonLatPos+1
209 end do
210
211 ! Top Level Control setup
212 call ram_setup
213
214 !- Initialize the Temporal grid and set dateStamp from env variable
215 call tim_setup()
216 if (tim_getDateStamp() == 0) then
217 call utl_abort('midas-diagBmatrix: date must be set by env variable MIDAS_DATE')
218 end if
219
220 ! Build date-time string from dateStamp
221 dateStamp = tim_getDateStamp()
222 ierr = newdate(dateStamp,datePrint,timePrint,-3)
223 dateTime = datePrint*100 + timePrint/1000000
224 write(datestr,'(i10.10)') dateTime
225 write(*,*)' datePrint= ',datePrint,' timePrint= ',timePrint
226 write(*,*)' date= ',dateTime,' stamp= ',dateStamp
227
228 ! Initialize variables of the model states
229 call gsv_setup
230
231 ! Initialize the Analysis horizontal grid
232 call hco_SetupFromFile( hco_anl,'./analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
233
234 if ( hco_anl % global ) then
235 hco_core => hco_anl
236 else
237 !- Iniatilized the core (Non-Exteded) analysis grid
238 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
239 end if
240
241 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
242
243 ! Initialize the vertical coordinate from the statistics file
244 call vco_SetupFromFile( vco_anl, & ! OUT
245 './analysisgrid') ! IN
246
247 ! Allocate the statevector
248 call gsv_allocate(statevector, tim_nstepobsinc, hco_anl, vco_anl, &
249 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
250 allocHeight_opt=.false., allocPressure_opt=.false.)
251 call gsv_zero(statevector)
252 nkgdim = statevector%nk
253
254 ! Setup the B matrix
255 call bmat_setup(hco_anl,hco_core,vco_anl)
256
257 !- Initialize the gridded variable transform module
258 call gvt_setup(hco_anl,hco_core,vco_anl)
259 if ( gsv_varExist(varName='HU') ) call gvt_setupRefFromTrialFiles('HU')
260
261 ! Setup of the L matrix done in bmat_setup
262 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
263
264 !
265 !==============================================
266 !- Compute columns of B and L matrices
267 !==============================================
268 !
269 if ( nLevs >= 1 .and. nLonLatPos >= 1 ) then
270
271 !
272 !- Compute columns of the B matrix
273 !
274 select case(trim(oneobs_timeStep))
275 case ('first')
276 oneobs_timeStepIndex = 1
277 case ('middle')
278 if (mod(tim_nstepobsinc,2) == 0) then
279 write(*,*)
280 write(*,*) 'odd number of nstepobsinc a required for obs place in the middle of the analysis window '
281 write(*,*) 'tim_nstepobsinc = ', tim_nstepobsinc
282 call utl_abort('midas-diagBmatrix')
283 end if
284 oneobs_timeStepIndex = (tim_nstepobsinc+1)/2
285 case ('last')
286 oneobs_timeStepIndex = tim_nstepobsinc
287 case default
288 write(*,*)
289 write(*,*) 'Unsupported oneobs_timeStep : ', trim(oneobs_timeStep)
290 call utl_abort('midas-diagBmatrix')
291 end select
292
293 allocate(controlVector(cvm_nvadim))
294
295 write(*,*) '********************************************'
296 write(*,*) 'midas-diagBmatrix: Compute columns of B matrix'
297 write(*,*) '********************************************'
298 write(*,*)
299 write(*,*) ' temporal location = ',trim(oneobs_timeStep), oneobs_timeStepIndex
300 write(*,*) ' number of levels = ',nLevs
301 write(*,*) ' number of lon-lat positions = ', nLonLatPos
302
303 do varIndex = 1, vnl_numvarmax
304
305 if ( .not. gsv_varExist(varName=vnl_varNameList(varIndex)) ) cycle
306 if ( trim(oneobs_varName) /= 'all' .and. (trim(oneobs_varName) /= trim(vnl_varNameList(varIndex))) ) cycle
307
308 filenameInc = 'columnB_' // trim(vnl_varNameList(varIndex)) // '_' // datestr // '.fst'
309 filenameIncNorm= 'columnBnorm_' // trim(vnl_varNameList(varIndex)) // '_' // datestr // '.fst'
310 filenameEnsAmp = 'ensAmplitude_' // trim(vnl_varNameList(varIndex)) // '_'
311
312 write(*,*)
313 write(*,*) 'midas-diagBmatrix: simulating a pseudo-observation of ', trim(vnl_varNameList(varIndex))
314
315 if(vnl_varLevelFromVarname(vnl_varNameList(varIndex)).eq.'SF') then
316 nlevs2 = 1
317 else
318 nlevs2 = nlevs
319 end if
320
321 ip3 = 0
322 do levIndex = 1, nlevs2
323 do lonLatPosIndex = 1, nLonLatPos
324
325 latIndex = oneobs_lonlat(lonLatPosIndex,latPosIndex)
326 lonIndex = oneobs_lonlat(lonLatPosIndex,lonPosIndex)
327
328 call gsv_zero(statevector)
329 call gsv_getField(statevector,field4d,vnl_varNameList(varIndex))
330
331 if ( latIndex >= statevector%myLatBeg .and. latIndex <= statevector%myLatEnd .and. &
332 lonIndex >= statevector%myLonBeg .and. lonIndex <= statevector%myLonEnd ) then
333 if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'SF') then
334 field4d(lonIndex,latIndex,1 ,oneobs_timeStepIndex) = 1.0D0
335 else
336 field4d(lonIndex,latIndex,oneobs_levs(levIndex),oneobs_timeStepIndex) = 1.0D0
337 end if
338 end if
339
340 controlVector(:)=0.0d0
341 call bmat_sqrtBT(controlVector,cvm_nvadim,statevector)
342 call bmat_sqrtB (controlVector,cvm_nvadim,statevector)
343
344 write(*,*)'midas-diagBmatrix: writing out the column of B, levIndex,lonIndex,latIndex=',levIndex,lonIndex,latIndex
345
346 ip3 = ip3 + 1
347
348 do stepIndexInc = 1, tim_nstepobsinc
349 call gio_writeToFile(statevector,filenameInc,'1OBS_'//trim(vnl_varNameList(varIndex)), &
350 stepIndex_opt=stepIndexInc, ip3_opt=ip3, unitConversion_opt=.true.)
351 end do
352
353 ! Normalized the result to get correlation-like pattern
354 centralValueLocal = 0.d0
355 if ( latIndex >= statevector%myLatBeg .and. latIndex <= statevector%myLatEnd .and. &
356 lonIndex >= statevector%myLonBeg .and. lonIndex <= statevector%myLonEnd ) then
357 if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)).eq.'SF') then
358 centralValueLocal = field4d(lonIndex,latIndex,1 ,oneobs_timeStepIndex)
359 else
360 centralValueLocal = field4d(lonIndex,latIndex,oneobs_levs(levIndex),oneobs_timeStepIndex)
361 end if
362 end if
363 call rpn_comm_allreduce(centralValueLocal, centralValue, 1, &
364 "MPI_DOUBLE_PRECISION", "MPI_SUM", "GRID", ierr)
365
366 write(*,*) 'midas-diagBmatrix: centralValue found = ', centralValue
367
368 if (centralValue /= 0.d0) then
369 call gsv_scale(statevector,1.d0/centralValue)
370 else
371 call utl_abort('midas-diagBmatrix: central value equals 0!')
372 end if
373
374 do stepIndexInc = 1, tim_nstepobsinc
375 call gio_writeToFile(statevector,filenameIncNorm,'1OBSNRM_'//trim(vnl_varNameList(varIndex)), &
376 stepIndex_opt=stepIndexInc, ip3_opt=ip3, &
377 unitConversion_opt=.false.)
378 end do
379
380 ! Write the ensemble amplitude fields (i.e., the alphas) when Bens is active
381 if (writeEnsAmplitude) call ben_writeAmplitude('./',filenameEnsAmp, ip3) ! IN
382
383 end do
384 end do
385 end do
386
387 deallocate(controlVector)
388
389 !
390 !- Compute columns of the L matrix
391 !
392 numLoc = ben_getNumLoc()
393 numBensInstance = ben_getNumInstance()
394 numStepAmplitude = ben_getNumStepAmplitudeAssimWindow()
395 amp3dStepIndex = ben_getAmp3dStepIndexAssimWindow()
396
397 if (numStepAmplitude > 1) then
398 allocate(datestampList(numStepAmplitude))
399 call tim_getstamplist(dateStampList,numStepAmplitude,tim_getDatestamp())
400 nEns = ben_getnEns()
401 adv_amplitudeAssimWindow => ben_getAmplitudeAssimWindow()
402 else
403 allocate(datestampList(1))
404 datestampList(1) = dateStamp
405 end if
406
407 do instanceIndex = 1, numBensInstance
408 do locIndex = 1, numLoc ! (this loop will be done only when localization is used in B)
409 loc => ben_getLoc(locIndex,instanceIndex_opt=instanceIndex)
410
411 if (loc%vco%Vcode == 5002 .or. loc%vco%Vcode == 5005) then
412 varNameALFA(:) = varNameALFAatm(:)
413 else ! vco_anl%Vcode == -1
414 varNameALFA(:) = varNameALFAsfc(:)
415 end if
416
417 call mmpi_setup_latbands(loc%hco%nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
418 call mmpi_setup_lonbands(loc%hco%ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
419
420 call ens_allocate(ensAmplitude, loc%nEnsOverDimension, numStepAmplitude, loc%hco, loc%vco, &
421 datestampList=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8)
422
423 call gsv_allocate(statevectorEnsAmplitude, numStepAmplitude, loc%hco, loc%vco, &
424 dateStampList_opt=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8, &
425 mpi_local_opt=.true.)
426
427 allocate(controlVector(loc%cvDim))
428
429 write(*,*) '********************************************'
430 write(*,*) 'midas-diagBmatrix: Compute columns of L matrix'
431 write(*,*) '********************************************'
432
433 write(*,*) ' number of levels = ', nlevs
434 write(*,*) ' number of lon-lat positions = ', nLonLatPos
435
436 if (numLoc > 1) then
437 write(locIndexString,'(i1)') locIndex
438 filename = 'columnL_' // trim(loc%locType) // '_' // locIndexString // '_' // datestr // '.fst'
439 else
440 if (numBensInstance > 1) then
441 write(instanceIndexString,'(I2.2)') instanceIndex
442 filename = 'columnL_i' // trim(instanceIndexString) // '_' // trim(loc%locType) // '_' // datestr // '.fst'
443 else
444 filename = 'columnL_' // trim(loc%locType) // '_' // datestr // '.fst'
445 end if
446 end if
447
448 ip3 = 0
449 do levIndex = 1, nlevs2
450 do lonLatPosIndex = 1, nLonLatPos
451
452 latIndex = oneobs_lonlat(lonLatPosIndex,latPosIndex)
453 lonIndex = oneobs_lonlat(lonLatPosIndex,lonPosIndex)
454
455 call ens_zero(ensAmplitude)
456 call gsv_zero(statevectorEnsAmplitude)
457 if ( latIndex >= statevector%myLatBeg .and. latIndex <= statevector%myLatEnd .and. &
458 lonIndex >= statevector%myLonBeg .and. lonIndex <= statevector%myLonEnd ) then
459 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,oneobs_levs(levIndex))
460 ensAmplitude_oneLev(:,oneobs_timeStepIndex,lonIndex,latIndex) = 1.d0
461 end if
462 controlVector(:)=0.0d0
463
464 if (numStepAmplitude > 1) then
465 call adv_ensemble_ad(ensAmplitude, & ! INOUT
466 adv_amplitudeAssimWindow, nEns) ! IN
467 end if
468
469 call loc_LsqrtAd(loc, & ! IN
470 ensAmplitude, & ! IN
471 controlVector, & ! OUT
472 amp3dStepIndex) ! IN
473 call loc_Lsqrt (loc, & ! IN
474 controlVector, & ! IN
475 ensAmplitude, & ! OUT
476 amp3dStepIndex) ! IN
477
478 if (numStepAmplitude > 1) then
479 call adv_ensemble_tl(ensAmplitude, & ! INOUT
480 adv_amplitudeAssimWindow, nEns) ! IN
481 end if
482
483 write(*,*)'midas-diagBmatrix: writing out the column of L, levIndex,lonIndex,latIndex=',levIndex,lonIndex,latIndex
484 call flush(6)
485
486 ip3 = ip3 + 1
487 call ens_copyMember(ensAmplitude, statevectorEnsAmplitude, 1)
488 do stepIndexInc = 1, numStepAmplitude
489 call gio_writeToFile(statevectorEnsAmplitude,filename,'ONEOBS', &
490 stepIndex_opt=stepIndexInc, ip3_opt=ip3,unitConversion_opt=.false.)
491 end do
492
493 end do
494 end do
495
496 deallocate(controlVector)
497 call ens_deallocate(ensAmplitude)
498 call gsv_deallocate(statevectorEnsAmplitude)
499
500 end do ! localization index (if active in B)
501 end do ! Bens instance
502
503 end if ! if any oneobs selected
504
505 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
506
507 !
508 !==============================================
509 !- Compute the stddev from random perturbations
510 !==============================================
511 !
512 if ( numperturbations > 1 ) then
513
514 write(*,*) '********************************************'
515 write(*,*) 'midas-diagBmatrix: Compute the stddev from random perturbations'
516 write(*,*) '********************************************'
517
518 allocate(controlVector(cvm_nvadim))
519 allocate(controlVector_global(cvm_nvadim_mpiglobal))
520
521 call gsv_zero(statevector)
522
523 ! Allocate the randomEns, mean and stddev
524 allocate(randomEns(statevector%myLonBeg:statevector%myLonEnd,statevector%myLatBeg:statevector%myLatEnd,nkgdim,numperturbations))
525 allocate(mean(statevector%myLonBeg:statevector%myLonEnd,statevector%myLatBeg:statevector%myLatEnd,nkgdim))
526 allocate(stddev(statevector%myLonBeg:statevector%myLonEnd,statevector%myLatBeg:statevector%myLatEnd,nkgdim))
527
528 iseed = abs(nrandseed)
529 call rng_setup(iseed)
530
531 call gsv_getField(statevector,field3d)
532
533 !
534 !- Compute the ensemble of random perturbations
535 !
536 do ensIndex = 1, numperturbations
537 write(*,*) 'midas-diagBmatrix: computing member number= ',ensIndex
538 call flush(6)
539
540 !- Global vector (same for each processors)
541 do index = 1, cvm_nvadim_mpiglobal
542 controlVector_global(index) = rng_gaussian()
543 end do
544
545 !- Extract only the subvector for this processor
546 call bmat_reduceToMPILocal(controlVector, & ! OUT
547 controlVector_global ) ! IN
548
549 !- Transform to control variables in physical space
550 call bmat_sqrtB(controlVector,cvm_nvadim,statevector)
551
552 if ( writePsiChiStddev ) call gvt_transform(statevector,'UVtoPsiChi')
553
554 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
555 do kIndex = 1, nkgdim
556 do latIndex = statevector%myLatBeg, statevector%myLatEnd
557 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
558 randomEns(lonIndex,latIndex,kIndex,ensIndex) = field3d(lonIndex,latIndex,kIndex)
559 end do
560 end do
561 end do
562 !$OMP END PARALLEL DO
563
564 end do ! Loop on ens member
565
566 !
567 !- Compute the ensemble mean
568 !
569 mean(:,:,:) = 0.0d0
570 do ensIndex = 1, numperturbations
571 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
572 do kIndex = 1, nkgdim
573 do latIndex = statevector%myLatBeg, statevector%myLatEnd
574 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
575 mean(lonIndex,latIndex,kIndex) = mean(lonIndex,latIndex,kIndex) + randomEns(lonIndex,latIndex,kIndex,ensIndex)
576 end do
577 end do
578 end do
579 !$OMP END PARALLEL DO
580 end do
581
582 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
583 do kIndex = 1, nkgdim
584 do latIndex = statevector%myLatBeg, statevector%myLatEnd
585 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
586 mean(lonIndex,latIndex,kIndex) = mean(lonIndex,latIndex,kIndex)/real(numperturbations,8)
587 end do
588 end do
589 end do
590 !$OMP END PARALLEL DO
591
592 !
593 !- Remove the ensemble mean from the ensemble
594 !
595 !$OMP PARALLEL DO PRIVATE (lonIndex,ensIndex,latIndex,kIndex)
596 do ensIndex = 1, numperturbations
597 do kIndex = 1, nkgdim
598 do latIndex = statevector%myLatBeg, statevector%myLatEnd
599 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
600 randomEns(lonIndex,latIndex,kIndex,ensIndex) = randomEns(lonIndex,latIndex,kIndex,ensIndex) - mean(lonIndex,latIndex,kIndex)
601 end do
602 end do
603 end do
604 end do
605 !$OMP END PARALLEL DO
606 deallocate(mean)
607
608 !
609 !- Compute the ensemble stddev
610 !
611 stddev(:,:,:) = 0.0d0
612
613 do ensIndex = 1, numperturbations
614 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
615 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
616 do latIndex = statevector%myLatBeg, statevector%myLatEnd
617 do kIndex = 1, nkgdim
618 stddev(lonIndex,latIndex,kIndex) = stddev(lonIndex,latIndex,kIndex) + &
619 (randomEns(lonIndex,latIndex,kIndex,ensIndex)**2)/real(numperturbations,8)
620 end do
621 end do
622 end do
623 !$OMP END PARALLEL DO
624 end do
625 deallocate(randomEns)
626
627 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
628 do kIndex = 1, nkgdim
629 do latIndex = statevector%myLatBeg, statevector%myLatEnd
630 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
631 stddev(lonIndex,latIndex,kIndex) = sqrt(stddev(lonIndex,latIndex,kIndex))
632 end do
633 end do
634 end do
635 !$OMP END PARALLEL DO
636
637 !- Insert results in statevector
638 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
639 do kIndex = 1, nkgdim
640 do latIndex = statevector%myLatBeg, statevector%myLatEnd
641 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
642 field3d(lonIndex,latIndex,kIndex) = stddev(lonIndex,latIndex,kIndex)
643 end do
644 end do
645 end do
646 !$OMP END PARALLEL DO
647 deallocate(stddev)
648
649 !- Write to file
650 call gio_writeToFile(statevector,'stddev_' // datestr // '.fst','GD_STDDEV', &
651 unitConversion_opt=.true.)
652
653 !
654 !- Compute the zonal mean std dev
655 !
656 write(*,*) 'midas-diagBmatrix: Compute the zonal mean stddev'
657 call flush(6)
658
659 allocate(stddev_zm(hco_anl%nj,nkgdim))
660 allocate(stddev_zm2(hco_anl%nj,nkgdim))
661 stddev_zm(:,:) = 0.d0
662 stddev_zm2(:,:) = 0.d0
663
664 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
665 do kIndex = 1, nkgdim
666 do latIndex = statevector%myLatBeg, statevector%myLatEnd
667 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
668 stddev_zm(latIndex,kIndex) = stddev_zm(latIndex,kIndex) + (field3d(lonIndex,latIndex,kIndex)**2)/real(hco_anl%ni,8)
669 end do
670 end do
671 end do
672 !$OMP END PARALLEL DO
673
674 nsize = statevector%nj*nkgdim
675 call rpn_comm_allreduce(stddev_zm,stddev_zm2,nsize, &
676 "MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
677
678 !- Insert results in statevector
679 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
680 do kIndex = 1, nkgdim
681 do latIndex = statevector%myLatBeg, statevector%myLatEnd
682 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
683 if(stddev_zm2(latIndex,kIndex)> 0.0d0) then
684 field3d(lonIndex,latIndex,kIndex) = sqrt(stddev_zm2(latIndex,kIndex))
685 else
686 field3d(lonIndex,latIndex,kIndex) = 0.0d0
687 end if
688 end do
689 end do
690 end do
691 !$OMP END PARALLEL DO
692 deallocate(stddev_zm)
693 deallocate(stddev_zm2)
694
695 call gio_writeToFile(statevector,'stddev_' // datestr // '.fst','ZM_STDDEV', &
696 unitConversion_opt=.true.)
697
698 ! Write the zonal mean stddev to a text file, if requested
699 if ( writeTextStddev ) then
700 allocate(zonalMeanStddev(statevector%latPerPE * mmpi_npey))
701 do varIndex = 1, vnl_numvarmax
702 if (.not. gsv_varExist(varName=vnl_varNameList(varIndex)) ) cycle
703
704 write(*,*) 'midas-diagBmatrix: writing zonal mean stddev to text file for variable: ', vnl_varNameList(varIndex)
705 call gsv_getField(statevector,field3d,vnl_varNameList(varIndex))
706
707 varName = vnl_varNameList(varIndex)
708 if ( varName == 'HU' ) varName = 'LQ'
709 filename = 'stddev_ZM_' // trim(varName) // '_' // datestr // '.txt'
710 nultxt = 0
711 if ( mmpi_myid == 0 ) ierr = fnom(nultxt,trim(filename),'FTN',0)
712
713 do levIndex = 1, gsv_getNumLevFromVarName(statevector,vnl_varNameList(varIndex))
714 nsize = statevector%latPerPE
715 call rpn_comm_gather(field3d(1,:,levIndex), nsize, 'mpi_double_precision', &
716 zonalMeanStddev, nsize, 'mpi_double_precision', 0, 'NS', ierr )
717 if ( mmpi_myid == 0 ) then
718 do latIndex = 1, statevector%nj
719 write(nultxt,*) field3d(1,latIndex,levIndex)
720 end do
721 end if
722 end do
723
724 if ( mmpi_myid == 0 ) ierr = fclos(nulnam)
725
726 end do
727 deallocate(zonalMeanStddev)
728
729 end if
730
731 !
732 !- Compute the domain mean std dev
733 !
734 write(*,*) 'midas-diagBmatrix: Compute the domain mean stddev'
735 call flush(6)
736
737 allocate(stddev_dm(hco_anl%ni,nkgdim))
738 allocate(stddev_dm2(hco_anl%ni,nkgdim))
739 stddev_dm(:,:) = 0.d0
740 stddev_dm2(:,:) = 0.d0
741
742 !$OMP PARALLEL DO PRIVATE (latIndex,kIndex)
743 do kIndex = 1, nkgdim
744 do latIndex = statevector%myLatBeg, statevector%myLatEnd
745 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
746 stddev_dm(lonIndex,kIndex) = stddev_dm(lonIndex,kIndex) + (field3d(lonIndex,latIndex,kIndex)**2)/real(hco_anl%nj,8)
747 end do
748 end do
749 end do
750 !$OMP END PARALLEL DO
751
752 nsize = statevector%ni*nkgdim
753 call rpn_comm_allreduce(stddev_dm,stddev_dm2,nsize, &
754 "MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
755
756 !- Insert results in statevector
757 !$OMP PARALLEL DO PRIVATE (lonIndex,latIndex,kIndex)
758 do kIndex = 1, nkgdim
759 do latIndex = statevector%myLatBeg, statevector%myLatEnd
760 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
761 if(stddev_dm2(lonIndex,kIndex)> 0.0d0) then
762 field3d(lonIndex,latIndex,kIndex) = sqrt(stddev_dm2(lonIndex,kIndex))
763 else
764 field3d(lonIndex,latIndex,kIndex) = 0.0d0
765 end if
766 end do
767 end do
768 end do
769 !$OMP END PARALLEL DO
770 deallocate(stddev_dm)
771 deallocate(stddev_dm2)
772 deallocate(controlVector)
773 deallocate(controlVector_global)
774
775 call gio_writeToFile(statevector,'stddev_' // datestr // '.fst','DM_STDDEV', &
776 unitConversion_opt=.true.)
777
778 end if ! if numperturbations.gt.1
779
780 call gsv_deallocate(statevector)
781
782 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
783
784 ! MPI, tmg finalize
785 call utl_tmg_stop(0)
786 call tmg_terminate(mmpi_myid, 'TMG_INFO')
787 call rpn_comm_finalize(ierr)
788
789 write(*,*) ' --------------------------------'
790 write(*,*) ' midas-diagBmatrix ENDS'
791 write(*,*) ' --------------------------------'
792
793end program midas_diagBmatrix
