1module calcStatsLam_mod
2 ! MODULE calcStatsLam_mod (prefix='csl' category='1. High-level functionality')
3 !
4 !:Purpose: To compute homogeneous and isotropic background error covariances
5 ! from forecast error estimate in model variable space (limited-area
6 ! version).
7 !
8 use mathPhysConstants_mod
9 use earthConstants_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use lamSpectralTransform_mod
13 use lamAnalysisGridTransforms_mod
14 use horizontalCoord_mod
15 use verticalCoord_mod
16 use localizationFunction_mod
17 use utilities_mod
18 use menetrierDiag_mod
19 use ensemblestatevector_mod
20 use gridVariableTransforms_mod
21 use varNameList_mod
22 use gridBinning_mod
23 use timeCoord_mod
24 use midasMpi_mod
25 use calcHeightAndPressure_mod
26 implicit none
27 save
28 private
29
30 ! Public Subroutines
31 public :: csl_setup, csl_computeBhi
32 public :: csl_toolbox
33
34 type(struct_hco), pointer :: hco_ens => null() ! Ensemble horizontal grid parameters
35 type(struct_hco), pointer :: hco_bhi => null() ! B matrix horizontal grid parameters
36 type(struct_vco), pointer :: vco_bhi => null() ! B matrix vertical grid parameters
37
38 integer,external :: get_max_rss
39
40 integer, parameter :: cv_model = 1
41 integer, parameter :: cv_bhi = 2
42 integer, parameter :: nMaxControlVar = 10
43
44 type :: struct_cv
45 character(len=4) :: NomVar(2)
46 integer :: varLevIndexStart
47 integer :: varLevIndexEnd
48 integer :: nlev
49 character(len=2) :: GridType
50 integer, allocatable :: ip1(:)
51 end type struct_cv
52
53 type :: struct_bhi
54 type(struct_cv) :: controlVariable(nMaxControlVar)
55 character(len=4) :: momentumControlVar(2)
56 integer :: nControlVariable
57 integer :: nVarLev
58 end type struct_bhi
59
60 type(struct_bhi) :: bhi
61
62 integer :: nEns
63 integer :: ip2_ens
64
65 logical :: initialized = .false.
66 logical :: vertLoc, horizLoc, stdDevScaling
67 logical :: ensContainsFullField
68
69 character(len=2) :: ctrlVarHumidity
70
71 type(struct_ens) :: ensPerts
72
73 real(8), pointer :: pressureProfile_M(:), pressureProfile_T(:)
74
75
76 real(8),allocatable :: scaleFactor_M(:), scaleFactor_T(:)
77 real(8) :: scaleFactor_SF
78
79 ! Namelist variables
80 integer :: nTrunc
81 character(len=12) :: WindTransform
82 logical :: NormByStdDev
83 logical :: SetTGtoZero
84 logical :: writeEnsPert
85 character(len=12) :: SpectralWeights
86 real(8) :: vLocalize_wind ! vertical length scale (in units of ln(Pressure))
87 real(8) :: vlocalize_mass ! vertical length scale (in units of ln(Pressure))
88 real(8) :: vlocalize_humidity ! vertical length scale (in units of ln(Pressure))
89 real(8) :: vlocalize_other ! vertical length scale (in units of ln(Pressure))
90 real(8) :: hlocalize_wind ! horizontal length scale (in km)
91 real(8) :: hlocalize_mass ! horizontal length scale (in km)
92 real(8) :: hlocalize_humidity ! horizontal length scale (in km)
93 real(8) :: hlocalize_other ! horizontal length scale (in km)
94 character(len=4) :: correlatedVariables(vnl_numvarmax)
95
96contains
97
98 !--------------------------------------------------------------------------
99 ! csl_setup
100 !--------------------------------------------------------------------------
101 subroutine csl_setup(nEns_in, hco_ens_in, vco_ens_in, ip2_in)
102 !
103 !:Purpose: To initialize this module
104 !
105 implicit none
106
107 ! Arguments:
108 integer, intent(in) :: nEns_in
109 type(struct_vco), pointer, intent(in) :: vco_ens_in
110 type(struct_hco), pointer, intent(in) :: hco_ens_in
111 integer, intent(in) :: ip2_in
112
113 ! Locals:
114 integer :: nulnam, ier
115 integer :: fclos, fnom
116 integer :: varIndex, levIndex, k
117 integer :: numStep
118 integer, allocatable :: dateStampList(:)
119 real(8) :: SurfacePressure
120 character(len=256) :: enspathname
121 logical :: makeBiPeriodic
122 character(len=4), pointer :: controlVarNames(:)
123
124 ! Namelist variables (local)
125 real(8) :: scaleFactor(vco_maxNumLevels)
126 integer :: grd_ext_x
127 integer :: grd_ext_y
128
129 NAMELIST /NAMCALCSTATS_LAM/nTrunc,grd_ext_x,grd_ext_y,WindTransform,NormByStdDev, &
130 SetTGtoZero,writeEnsPert,SpectralWeights, &
131 vLocalize_wind,vlocalize_mass,vlocalize_humidity, &
132 hLocalize_wind,hlocalize_mass,hlocalize_humidity, &
133 hLocalize_other,vlocalize_other, &
134 correlatedVariables,scaleFactor
135
136 write(*,*)
137 write(*,*) 'csl_setup: Starting...'
138 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
139
140 !
141 !- 1. Initialized the info on the ensemble
142 !
143 nEns=nEns_in
144
145 hco_ens => hco_ens_in
146 vco_bhi => vco_ens_in
147
148 if ( vco_bhi%nlev_T == 1 .and. vco_bhi%nlev_T == 1 ) then
149 write(*,*)
150 write(*,*) 'Spatial dimensions = 2D'
151 else
152 write(*,*)
153 write(*,*) 'Spatial dimensions = 3D'
154 end if
155
156 ip2_ens = ip2_in
157
158 !
159 !- 2. Read namelist NAMCALCSTATS_LAM
160 !
161 nTrunc = 75 ! Default value
162 grd_ext_x = 10 ! Default value
163 grd_ext_y = 10 ! Default value
164 WindTransform = 'PsiChi' ! Default value
165 NormByStdDev = .true. ! Default value
166 SetTGtoZero = .false. ! Default value
167 writeEnsPert = .false. ! Default value
168 correlatedVariables = ' '
169 SpectralWeights = 'lst' ! Default value
170 vLocalize_wind = -1.0d0 ! Default value (no vloc)
171 vLocalize_mass = -1.0d0 ! Default value (no vloc)
172 vLocalize_humidity = -1.0d0 ! Default value (no vloc)
173 vLocalize_other = -1.0d0 ! Default value (no vloc)
174 hLocalize_wind = -1.0d0 ! Default value (no hloc)
175 hLocalize_mass = -1.0d0 ! Default value (no hloc)
176 hLocalize_humidity = -1.0d0 ! Default value (no hloc)
177 hLocalize_other = -1.0d0 ! Default value (no hloc)
178 scaleFactor(:) = 1.0d0 ! Default value (no scaling)
179
180 nulnam = 0
181
182 ier=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
183 read (nulnam,nml=namcalcstats_lam)
184 write(* ,nml=namcalcstats_lam)
185 ier=fclos(nulnam)
186
187 write(*,*)
188 write(*,*) 'Truncation = ', nTrunc
189
190 select case(trim(SpectralWeights))
191 case ('lst')
192 write(*,*)
193 write(*,*) 'Using spectral weights from the lamSpectralTransform module'
194 case ('legacy')
195 write(*,*)
196 write(*,*) 'Using spectral weights from the OLD VAR code in LAM mode'
197 case default
198 write(*,*)
199 write(*,*) 'Unknown spectral weights TYPE : ', trim(SpectralWeights)
200 call utl_abort('csl_setup')
201 end select
202
203 !
204 !- 3. Initialized the extended (bi-periodic) grid
205 !
206
207 !- 3.1 Create the extended and non-extended grid prototype file
208 if (mmpi_myid == 0) then
209 call lgt_createLamTemplateGrids('./analysisgrid', hco_ens, vco_bhi, & ! IN
210 grd_ext_x, grd_ext_y) ! IN
211 end if
212 call rpn_comm_barrier("GRID",ier)
213
214 !- 3.2 Setup the Extended B_HI grid
215 call hco_setupFromFile( hco_bhi,'./analysisgrid', 'ANALYSIS', 'BHI' ) ! IN
216
217 !- 3.3 Setup the LAM analysis grid metrics
218 call lgt_setupFromHCO( hco_bhi, hco_ens ) ! IN
219
220 !
221 !- 4. Read the ensemble
222 !
223 numStep = 1
224 allocate(dateStampList(numStep))
225 dateStampList(:) = -1
226 call ens_allocate(ensPerts, nEns, numStep, hco_bhi, vco_bhi, dateStampList, &
227 hco_core_opt=hco_ens)
228
229 ensContainsFullField = .false.
230 ctrlVarHumidity = 'LQ'
231 ensPathName = './ensemble'
232 makeBiPeriodic = .true.
233 call ens_readEnsemble(ensPerts, ensPathName, makeBiPeriodic, &
234 containsFullField_opt=ensContainsFullField)
235
236 if ( ctrlVarHumidity == 'LQ' .and. ens_varExist(ensPerts,'HU') .and. &
237 ensContainsFullField ) then
238 call gvt_transform(ensPerts,'HUtoLQ')
239 end if
240
241 !
242 !- 5. Setup the control variables (model space and B_hi space)
243 !
244 nullify(controlVarNames)
245 call ens_varNamesList(controlVarNames,ensPerts)
246
247 bhi%nControlVariable = size(controlVarNames) !count(mask)
248 write(*,*)
249 write(*,*) 'Number of Control Variables = ', bhi%nControlVariable
250
251 !- 5.1 Set ControlVariable structure
252 bhi%momentumControlVar(:) = 'NULL'
253 do varIndex = 1, bhi%nControlVariable
254
255 !- Set variable name
256 bhi%controlVariable(varIndex)%nomvar(cv_model) = trim(controlVarNames(varIndex))
257 if ( bhi%controlVariable(varIndex)%nomvar(cv_model) == 'UU' ) then
258 if ( trim(WindTransform) == 'PsiChi') then
259 write(*,*)
260 write(*,*) '--- Momentum Control Variables = Psi-Chi ---'
261 bhi%momentumControlVar(1) = 'PP'
262 else if ( trim(WindTransform) == 'VortDiv') then
263 write(*,*)
264 write(*,*) '--- Momentum Control Variables = Vort-Div ---'
265 bhi%momentumControlVar(1) = 'QR'
266 call utl_abort('Momentum Control Variables = Vort-Div not yest availble')
267 else if ( trim(WindTransform) == 'UV') then
268 write(*,*)
269 write(*,*) '--- Momentum Control Variables = U-V ---'
270 bhi%momentumControlVar(1) = 'UU'
271 else
272 write(*,*)
273 write(*,*) 'Wind Transform not available ', trim((WindTransform))
274 call utl_abort('csl_setup')
275 end if
276 bhi%controlVariable(varIndex)%nomvar(cv_bhi) = trim(bhi%momentumControlVar(1))
277 else if ( bhi%controlVariable(varIndex)%nomvar(cv_model) == 'VV' ) then
278 if ( trim(WindTransform) == 'PsiChi') then
279 bhi%momentumControlVar(2) = 'CC'
280 else if ( trim(WindTransform) == 'VortDiv') then
281 bhi%momentumControlVar(2) = 'DD'
282 else if ( trim(WindTransform) == 'UV') then
283 bhi%momentumControlVar(2) = 'VV'
284 else
285 write(*,*)
286 write(*,*) 'Wind Transform not available ', trim((WindTransform))
287 call utl_abort('csl_setup')
288 end if
289 bhi%controlVariable(varIndex)%nomvar(cv_bhi) = bhi%momentumControlVar(2)
290 else
291 bhi%controlVariable(varIndex)%nomvar(cv_bhi) = bhi%controlVariable(varIndex)%nomvar(cv_model)
292 end if
293
294 !- 5.2 Set Level info
295 bhi%controlVariable(varIndex)%GridType = vnl_varLevelFromVarname(bhi%controlVariable(varIndex)%nomvar(cv_model))
296 bhi%controlVariable(varIndex)%nlev = ens_getNumLev(ensPerts,bhi%controlVariable(varIndex)%GridType)
297
298 write(*,*)
299 write(*,*) 'Control Variable Name ', bhi%controlVariable(varIndex)%nomvar(cv_model)
300 write(*,*) ' Number of Levels = ', bhi%controlVariable(varIndex)%nlev
301 write(*,*) ' Type of Levels = ', bhi%controlVariable(varIndex)%GridType
302
303 allocate( bhi%controlVariable(varIndex)%ip1 (bhi%controlVariable(varIndex)%nlev) )
304
305 if (bhi%controlVariable(varIndex)%GridType == 'TH') then
306 bhi%controlVariable(varIndex)%ip1(:) = vco_bhi%ip1_T(:)
307 else if (bhi%controlVariable(varIndex)%GridType == 'MM') then
308 bhi%controlVariable(varIndex)%ip1(:) = vco_bhi%ip1_M(:)
309 else if (bhi%controlVariable(varIndex)%GridType == 'SF') then
310 bhi%controlVariable(varIndex)%ip1(:) = 0
311 end if
312
313 bhi%controlVariable(varIndex)%varLevIndexStart = &
314 ens_getOffsetFromVarName(ensPerts,bhi%controlVariable(varIndex)%nomvar(cv_model)) + 1
315 bhi%controlVariable(varIndex)%varLevIndexEnd = &
316 ens_getOffsetFromVarName(ensPerts,bhi%controlVariable(varIndex)%nomvar(cv_model)) + &
317 bhi%controlVariable(varIndex)%nlev
318
319 end do
320
321 bhi%nVarLev = ens_getNumK(ensPerts)
322
323 !
324 !- 6. Transform u-wind and v-wind to control variables
325 !
326 if (writeEnsPert) then
327 call ens_writeEnsemble(ensPerts, './', 'MODELVAR_', 'MODELVAR', 'E', &
328 containsFullField_opt = ensContainsFullField)
329 end if
330
331 if ( bhi%momentumControlVar(1) == 'PP' .and. bhi%momentumControlVar(2) == 'CC' ) then
332 call gvt_transform(ensPerts,'UVtoPsiChi')
333 call ens_modifyVarName(ensPerts, 'UU', 'PP')
334 call ens_modifyVarName(ensPerts, 'VV', 'CC')
335 else if ( bhi%momentumControlVar(1) == 'QR' .and. bhi%momentumControlVar(2) == 'DD' ) then
336 call gvt_transform(ensPerts,'UVtoVortDiv')
337 call ens_modifyVarName(ensPerts, 'UU', 'QR')
338 call ens_modifyVarName(ensPerts, 'VV', 'DD')
339 end if
340
341 if (writeEnsPert) then
342 call ens_writeEnsemble(ensPerts, './', 'CTRLVAR', 'CTRLVAR', 'E', &
343 containsFullField_opt = ensContainsFullField)
344 end if
345
346 !
347 !- 7. Compute and remove the ensemble mean; compute the stdDev
348 !
349 call ens_computeMean(ensPerts)
350 call ens_removeMean (ensPerts)
351 call ens_computeStdDev(ensPerts)
352
353 !
354 !- 8. Setup the localization
355 !
356
357 !- 8.1 Setup horizontal localization
358 if ( hLocalize_wind < 0.d0 .and. hLocalize_mass < 0.d0 .and. &
359 hLocalize_humidity < 0.d0 .and. hLocalize_other < 0.d0) then
360 write(*,*)
361 write(*,*) 'csl_setup: NO horizontal correlation localization will be performed'
362 horizLoc=.false.
363 else
364 write(*,*)
365 write(*,*) 'csl_setup: horizontal correlation localization WILL BE performed'
366 horizLoc=.true.
367 end if
368
369 !- 8.2 Setup vertical localization
370 if ( vLocalize_wind < 0.d0 .and. vLocalize_mass < 0.d0 .and. &
371 vLocalize_humidity < 0.d0 .and. vLocalize_other < 0.d0 ) then
372 write(*,*)
373 write(*,*) 'csl_setup: NO vertical correlation localization will be performed'
374 vertLoc=.false.
375 else
376 write(*,*)
377 write(*,*) 'csl_setup: vertical correlation localization WILL BE performed'
378 vertLoc=.true.
379 end if
380
381 !
382 !- 9. Setup pressure profile for vertical localization
383 !
384 if (vco_bhi%vgridPresent) then
385 SurfacePressure = 101000.D0
386 call czp_fetch1DLevels(vco_bhi, SurfacePressure, &
387 profM_opt=pressureProfile_M, profT_opt=pressureProfile_T)
388
389 write(*,*)
390 write(*,*) 'Pressure profile...'
391 do k = 1, vco_bhi%nlev_M
392 write(*,*) k, pressureProfile_M(k) / 100.d0, ' hPa'
393 end do
394
395 write(*,*)
396 write(*,*) 'Pressure profile...'
397 do k = 1, vco_bhi%nlev_T
398 write(*,*) k, pressureProfile_T(k) / 100.d0, ' hPa'
399 end do
400
401 end if
402
403 !
404 !- 10. Setup the scaling
405 !
406 if ( all(scaleFactor(:) == 1.d0) ) then
407 write(*,*)
408 write(*,*) 'csl_setup: NO scaling of the StdDev will be performed'
409 stdDevScaling=.false.
410
411 else
412 write(*,*)
413 write(*,*) 'csl_setup: scaling of the StdDev WILL BE performed'
414 stdDevScaling=.true.
415
416 allocate(scaleFactor_M(vco_bhi%nlev_M))
417 allocate(scaleFactor_T(vco_bhi%nlev_T))
418 do levIndex = 1, vco_bhi%nlev_T
419 if (scaleFactor(levIndex) > 0.0d0) then
420 scaleFactor_T(levIndex) = sqrt(scaleFactor(levIndex))
421 else
422 scaleFactor_T(levIndex) = 0.0d0
423 end if
424 end do
425 scaleFactor_M(1:vco_bhi%nlev_M) = scaleFactor_T(1:vco_bhi%nlev_M)
426 scaleFactor_SF = scaleFactor_T(vco_bhi%nlev_T)
427
428 end if
429
430 !
431 !- 11. Ending
432 !
433 initialized = .true.
434
435 write(*,*)
436 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
437 write(*,*) 'csl_setup: Done!'
438
439 end subroutine csl_setup
440
441 !--------------------------------------------------------------------------
442 ! csl_computeBhi
443 !--------------------------------------------------------------------------
444 subroutine csl_computeBhi
445 !
446 !:Purpose: To compute an homogeneous and isotopic B matrix
447 !
448 implicit none
449
450 ! Locals:
451 real(8), allocatable :: SpVertCorrel(:,:,:)
452 real(8), allocatable :: TotVertCorrel(:,:)
453 real(8), allocatable :: NormB(:,:,:)
454 real(8), allocatable :: NormBsqrt(:,:,:)
455 real(8), allocatable :: PowerSpectrum(:,:)
456 real(8), allocatable :: HorizScale(:)
457 character(len=4), pointer :: varNamesList(:)
458 type(struct_gbi) :: gbi_horizontalMean
459 type(struct_gsv) :: statevector_stdDev
460 type(struct_gsv) :: statevector_mean, statevector_stdDevGridPoint
461 integer :: ier
462
463 write(*,*)
464 write(*,*) 'csl_computeBhi: Starting...'
465 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
466
467 allocate(SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc))
468 allocate(TotVertCorrel(bhi%nVarLev,bhi%nVarLev))
469 allocate(PowerSpectrum(bhi%nVarLev,0:nTrunc))
470 allocate(NormB(bhi%nVarLev,bhi%nVarLev,0:nTrunc))
471 allocate(NormBsqrt(bhi%nVarLev,bhi%nVarLev,0:nTrunc))
472 allocate(HorizScale(bhi%nVarLev))
473
474 !
475 !- 1. Calculate the gridded binned Std. Dev. to be used in the analysis step
476 !
477 nullify(varNamesList)
478 call ens_varNamesList(varNamesList,ensPerts)
479 call gsv_allocate(statevector_stdDev, ens_getNumStep(ensPerts), &
480 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
481 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
482 dataKind_opt=8 )
483
484 call gbi_setup(gbi_horizontalMean, 'HorizontalMean', statevector_stdDev, hco_ens)
485
486 call gbi_stdDev(gbi_horizontalMean, ensPerts, & ! IN
487 statevector_stdDev) ! OUT
488
489 !
490 !- 2. Normalization of the ensemble perturbations
491 !
492 if ( NormByStdDev ) then
493 call ens_normalize(ensPerts)
494 end if
495
496 !
497 !- 3. Covariance statistics in Spectral Space
498 !
499
500 !- 3.1 Vertical correlations and Power Spectra
501 call calcSpectralStats(ensPerts, & ! IN
502 SpVertCorrel, PowerSpectrum, & ! OUT
503 NormB) ! OUT
504
505 if (mmpi_myid == 0) then
506 !- 3.2 Calculate the horiontal correlation lenght scales
507 call calcHorizScale(HorizScale, & ! OUT
508 NormB) ! IN
509
510 !- 3.3 Calculate the total vertical correlation matrix
511 call calcTotVertCorrel(TotVertCorrel, & ! OUT
512 SpVertCorrel, PowerSpectrum) ! IN
513
514 !- 3.4 Set cross-correlations
515 call setSpVertCorrel(NormB) ! INOUT
516
517 !- 3.5 Calculate the square-root of the correlation-based B matrix
518 call calcBsqrt(NormBsqrt, & ! OUT
519 NormB ) ! IN
520 end if
521
522 !- 3.6 Apply scaling
523 if (stdDevScaling) then
524 call scaleStdDev(statevector_stdDev) ! INOUT
525 end if
526
527 call rpn_comm_barrier("GRID",ier)
528
529 !
530 !- 4. Writing statistics to files
531 !
532
533 !- 4.1 Statistics needed by VAR in analysis mode
534 call writeVarStats(NormBsqrt, statevector_stdDev) ! IN
535
536 !- 4.2 Diagnostics fields
537 call ens_copyEnsMean(ensPerts, statevector_mean)
538 call ens_copyEnsStdDev(ensPerts, statevector_stdDevGridPoint)
539
540 call writeDiagStats(NormB, SpVertCorrel, TotVertCorrel, statevector_mean, & ! IN
541 statevector_stdDevGridPoint, PowerSpectrum, HorizScale) ! IN
542
543 !
544 !- 5. Ending
545 !
546 call ens_deallocate(ensPerts)
547 call gsv_deallocate(statevector_mean)
548 call gsv_deallocate(statevector_stdDev)
549 call gsv_deallocate(statevector_stdDevGridPoint)
550
551 deallocate(SpVertCorrel)
552 deallocate(TotVertCorrel)
553 deallocate(PowerSpectrum)
554 deallocate(NormB)
555 deallocate(NormBsqrt)
556 deallocate(HorizScale)
557
558 write(*,*)
559 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
560 write(*,*) 'csl_computeBhi: Done!'
561
562 end subroutine csl_computeBhi
563
564 !--------------------------------------------------------------------------
565 ! csl_toolbox
566 !--------------------------------------------------------------------------
567 subroutine csl_toolbox
568 !
569 !:Purpose: High-level control of various diagnostic tools
570 !
571 implicit none
572
573 ! Locals:
574 real(8), allocatable :: SpVertCorrel(:,:,:)
575 real(8), allocatable :: NormB(:,:,:)
576 real(8), allocatable :: PowerSpectrum(:,:)
577 integer :: nulnam, ier, fnom, fclos
578 type(struct_gsv) :: statevector_stdDev
579 type(struct_gsv) :: statevector_template
580 character(len=4), pointer :: varNamesList(:)
581 character(len=60) :: tool
582
583 NAMELIST /NAMTOOLBOX/tool
584
585 write(*,*)
586 write(*,*) 'csl_toolbox: Starting...'
587 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
588
589 !
590 !- 1. Tool selection
591 !
592 nulnam = 0
593 ier = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
594 read(nulnam,nml=NAMTOOLBOX)
595 write(*,nml=NAMTOOLBOX)
596 ier = fclos(nulnam)
597
598 select case(trim(tool))
599
600 case ('VERTCORREL_GRIDPOINT')
601 write(*,*)
602 write(*,*) 'Computing vertical correlation in grid point space'
603
604 call ens_normalize(ensPerts)
605
606 call calcVertCorrel(ensPerts)
607
608 case ('HORIZCORREL_FUNCTION')
609 write(*,*)
610 write(*,*) 'Computing horizontal correlation functions'
611
612 allocate(SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc))
613 allocate(PowerSpectrum(bhi%nVarLev,0:nTrunc))
614 allocate(NormB(bhi%nVarLev,bhi%nVarLev,0:nTrunc))
615
616 call ens_normalize(ensPerts)
617
618 call calcSpectralStats(ensPerts, & ! IN
619 SpVertCorrel, PowerSpectrum, & ! OUT
620 NormB) ! OUT
621
622 if (mmpi_myid == 0) then
623 call horizCorrelFunction(NormB) ! IN
624 end if
625
626 deallocate(NormB)
627 deallocate(PowerSpectrum)
628 deallocate(SpVertCorrel)
629
630 case ('STDDEV')
631 write(*,*)
632 write(*,*) 'Computing Standard-Deviations'
633
634 nullify(varNamesList)
635 call ens_varNamesList(varNamesList,ensPerts)
636 call gsv_allocate(statevector_stdDev, ens_getNumStep(ensPerts), &
637 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
638 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
639 dataKind_opt=8 )
640
641 call ens_copyEnsStdDev(ensPerts, statevector_stdDev)
642
643 call gio_writeToFile(statevector_stdDev, './stddev.fst', 'STDDEV_GRIDP', &
644 typvar_opt = 'E', numBits_opt = 32)
645
646 call gsv_deallocate(statevector_stdDev)
647
648 case ('HVCORREL_LOCAL')
649 write(*,*)
650 write(*,*) 'Computing Local Correlation'
651 call ens_normalize(ensPerts)
652 call calcLocalCorrelations(ensPerts)
653
654 case ('LOCALIZATIONRADII')
655 write(6,*)
656 write(6,*) 'Estimating the optimal covariance localization radii'
657
658 call ens_copyEnsStdDev(ensPerts, statevector_template)
659
660 call bmd_setup(statevector_template, hco_ens, nEns, pressureProfile_M, pressureProfile_T, 1)
661
662 call bmd_localizationRadii(ensPerts, waveBandIndex_opt=1) ! IN
663
664 case default
665 call utl_abort('csl_toolbox: Unknown TOOL '// trim(tool))
666 end select
667
668 !
669 !- 2. Write the estimated pressure profiles
670 !
671 if (mmpi_myid == 0 .and. vco_bhi%vgridPresent) then
672 call writePressureProfiles
673 end if
674
675 !
676 !- 3. Ending
677 !
678 call ens_deallocate(ensPerts)
679
680 write(*,*)
681 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
682 write(*,*) 'csl_toolbox: Done!'
683
684 end subroutine csl_toolbox
685
686 !--------------------------------------------------------------------------
687 ! calcSpectralStats
688 !--------------------------------------------------------------------------
689 subroutine calcSpectralStats(ensPerts,SpVertCorrel,PowerSpectrum, &
690 NormB)
691 !
692 !:Purpose: To compute background-error covariances in spectral space
693 ! from an ensemble of gridded data
694 !
695 implicit none
696
697 ! Arguments:
698 type(struct_ens), intent(inout) :: ensPerts
699 real(8), intent(out) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
700 real(8), intent(out) :: PowerSpectrum(bhi%nVarLev,0:nTrunc)
701 real(8), intent(out) :: NormB(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
702
703 ! Locals:
704 real(8), allocatable :: NormPowerSpectrum(:,:)
705 real(8), allocatable :: SpectralStateVar(:,:,:)
706 real(8), allocatable :: SpVertCorrel_local(:,:,:)
707 real(8), allocatable :: GridState(:,:,:)
708 real(8), allocatable :: SumWeight(:)
709 real(8), allocatable :: SumWeight_local(:)
710 type(struct_lst) :: lst_bhi ! Spectral transform Parameters
711 real(4), pointer :: ptr4d_r4(:,:,:,:)
712 real(8) :: weight
713 integer :: k1, k2, ens, e, ila, p, k, totwvnb
714 integer :: myLonBeg, myLonEnd, myLatBeg, myLatEnd
715 integer :: nSize, ier
716 character(len=24) :: kind
717
718 write(*,*)
719 write(*,*) 'CalcSpectralStats: Starting...'
720 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
721
722 !
723 !- 1. Setup the (LAM) spectral transform
724 !
725
726 call lst_setup(lst_bhi, & ! OUT
727 hco_bhi%ni, hco_bhi%nj, & ! IN
728 hco_bhi%dlon, nTrunc, & ! IN
729 'LatLonMN', maxlevels_opt=bhi%nVarLev) ! IN
730
731 !
732 !- 2. Calculate the Vertical Covariances in Spectral Space
733 !
734 call ens_getLatLonBounds(ensPerts, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
735
736 allocate( SpVertCorrel_local(bhi%nVarLev, bhi%nVarLev, 0:nTrunc) )
737 SpVertCorrel_local(:,:,:) = 0.d0
738
739 allocate( SpectralStateVar(lst_bhi%nla,lst_bhi%nphase,bhi%nVarLev) )
740 allocate( GridState(myLonBeg:myLonEnd, myLatBeg:myLatEnd, bhi%nVarLev) )
741
742 allocate(SumWeight_local(0:nTrunc))
743 SumWeight_local(:) = 0.d0
744
745 do ens = 1, nEns
746
747 write(6,*) 'ens = ', ens
748 flush(6)
749
750 !- 2.1 Extract fields from ensPerturbations
751 do k = 1, bhi%nVarLev
752 ptr4d_r4 => ens_getOneLev_r4(ensPerts,k)
753 GridState(:,:,k) = real(ptr4d_r4(ens,1,:,:),8)
754 end do
755
756 !- 2.2 Grid Point Space -> Spectral Space
757 kind = 'GridPointToSpectral'
758 call lst_VarTransform(lst_bhi, & ! IN
759 SpectralStateVar, & ! OUT
760 GridState, & ! IN
761 kind, bhi%nVarLev) ! IN
762
763 !- 2.3 Compute the covariances
764 !$OMP PARALLEL DO PRIVATE(totwvnb,weight,e,ila,p,k2,k1)
765 do totwvnb = 0, nTrunc
766 do e = 1, lst_bhi%nePerK(totwvnb)
767 ila = lst_bhi%ilaFromEK(e,totwvnb)
768 do p = 1, lst_bhi%nphase
769 if (trim(SpectralWeights) == 'lst') then
770 weight = lst_bhi%Weight(ila)
771 SumWeight_local(totwvnb) = SumWeight_local(totwvnb) + weight
772 else
773 weight = 2.0d0
774 end if
775 do k2 = 1, bhi%nVarLev
776 do k1 = 1, bhi%nVarLev
777 SpVertCorrel_local(k1,k2,totwvnb) = SpVertCorrel_local(k1,k2,totwvnb) &
778 + weight * (SpectralStateVar(ila,p,k1) * SpectralStateVar(ila,p,k2))
779 end do
780 end do
781 end do
782 end do
783 end do
784 !$OMP END PARALLEL DO
785
786 end do ! Loop in Ensemble
787
788 deallocate(SpectralStateVar)
789 deallocate(GridState)
790
791 ! Gather the all the info in processor 0
792 SpVertCorrel(:,:,:) = 0.d0
793 nSize = bhi%nVarLev * bhi%nVarLev * (nTrunc + 1)
794 call rpn_comm_reduce(SpVertCorrel_local,SpVertCorrel,nsize,"mpi_double_precision","mpi_sum",0,"GRID",ier)
795
796 allocate(SumWeight(0:nTrunc))
797 SumWeight(:) = 0.d0
798 nSize = nTrunc + 1
799 call rpn_comm_reduce(SumWeight_local, SumWeight, nsize,"mpi_double_precision","mpi_sum",0,"GRID",ier)
800
801 deallocate(SumWeight_local)
802 deallocate(SpVertCorrel_local)
803
804 if (mmpi_myid == 0) then
805
806 !- 2.4 Compute the weighted COVARIANCES for each total wavenumber
807 do totwvnb = 0, nTrunc
808 if (trim(SpectralWeights) == 'legacy') then
809 if (totwvnb /= 0 ) then
810 SumWeight(totwvnb) = 2.d0 * real(lst_bhi%nphase * lst_bhi%nePerKglobal(totwvnb),8) - 2.d0
811 else
812 SumWeight(totwvnb) = 1.d0
813 end if
814 SumWeight(totwvnb) = SumWeight(totwvnb)*nEns
815 end if
816
817 if ( SumWeight(totwvnb) /= 0.d0 ) then
818 SpVertCorrel(:,:,totwvnb) = SpVertCorrel(:,:,totwvnb) / SumWeight(totwvnb)
819 else
820 SpVertCorrel(:,:,totwvnb) = 0.d0
821 end if
822
823 end do
824
825 deallocate(SumWeight)
826
827 !- 2.5 Extract the power spectrum (the variances on the diagonal elements)
828 do k = 1, bhi%nVarLev
829 PowerSpectrum(k,:) = SpVertCorrel(k,k,:)
830 end do
831
832 !
833 !- 3. Calculate the Vertical Correlations in Spectral Space
834 !
835 !$OMP PARALLEL DO PRIVATE (totwvnb,k2,k1)
836 do totwvnb = 0, nTrunc
837 do k2 = 1, bhi%nVarLev
838 do k1 = 1, bhi%nVarLev
839 if ( PowerSpectrum(k1,totwvnb) /= 0.d0 .and. &
840 PowerSpectrum(k2,totwvnb) /= 0.d0 ) then
841 SpVertCorrel(k1,k2,totwvnb) = SpVertCorrel(k1,k2,totwvnb) / &
842 sqrt( PowerSpectrum(k1,totwvnb) * PowerSpectrum(k2,totwvnb) )
843 else
844 SpVertCorrel(k1,k2,totwvnb) = 0.d0
845 end if
846 end do
847 end do
848 end do
849 !$OMP END PARALLEL DO
850
851 ! Apply vertical localization (if wanted)
852 if (vertLoc) then
853 call applyVertLoc(SpVertCorrel) ! INOUT
854 end if
855
856 !
857 !- 4. Normalize the power spectrum (i.e. build normalised spectral densities of the variance)
858 !
859 allocate(NormPowerSpectrum(bhi%nVarLev,0:nTrunc))
860
861 call NormalizePowerSpectrum(PowerSpectrum, & ! IN
862 NormPowerSpectrum) ! OUT
863
864 ! Apply horizontal localization (if wanted)
865 if (horizLoc) then
866 call applyHorizLoc(NormPowerSpectrum) ! INOUT
867 end if
868
869 !
870 !- 5. Normalize the spectral vertical correlation matrix to ensure correlations in horizontal
871 !
872
873 !$OMP PARALLEL DO PRIVATE (totwvnb,k2,k1)
874 do totwvnb = 0, nTrunc
875 do k2 = 1, bhi%nVarLev
876 do k1 = 1, bhi%nVarLev
877 NormB(k1,k2,totwvnb) = SpVertCorrel(k1,k2,totwvnb) * &
878 sqrt( NormPowerSpectrum(k1,totwvnb) * NormPowerSpectrum(k2,totwvnb) )
879 end do
880 end do
881 end do
882 !$OMP END PARALLEL DO
883
884 deallocate(NormPowerSpectrum)
885
886 end if ! mmpi_myid == 0
887
888 write(*,*)
889 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
890 write(*,*) 'CalcSpectralStats: Done!'
891
892 end subroutine calcSpectralStats
893
894 !--------------------------------------------------------------------------
895 ! normalizePowerSpectrum
896 !--------------------------------------------------------------------------
897 subroutine normalizePowerSpectrum(PowerSpectrum, NormPowerSpectrum)
898 !
899 !:Purpose: To convert spectral variances into spectral correlations
900 !
901 implicit none
902
903 ! Arguments:
904 real(8), intent(in) :: PowerSpectrum(bhi%nVarLev,0:nTrunc)
905 real(8), intent(out) :: NormPowerSpectrum(bhi%nVarLev,0:nTrunc)
906
907 ! Locals:
908 real(8), allocatable :: SpectralStateVar(:,:,:)
909 real(8), allocatable :: GridState(:,:,:)
910 real(8) :: sum
911 integer :: e, ila, p, k, totwvnb
912 character(len=24) :: kind
913 type(struct_lst) :: lst_norm ! Spectral transform Parameters
914
915 write(*,*)
916 write(*,*) 'NormalizePowerSpectrum: Starting...'
917 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
918
919 !
920 !- 1. Setup the (LAM) spectral transform
921 !
922 call lst_setup(lst_norm, & ! OUT
923 hco_bhi%ni, hco_bhi%nj, & ! IN
924 hco_bhi%dlon, nTrunc, & ! IN
925 'NoMpi')
926
927 !
928 !- 1. Normalize the power spectrum (i.e. build normalised spectral densities of the variance)
929 !
930
931 !- 1.1 Part 1
932
933 !$OMP PARALLEL DO PRIVATE (totwvnb,k,sum)
934 do k = 1, bhi%nVarLev
935 sum = 0.0d0
936 do totwvnb = 0, nTrunc
937 sum = sum + real(totwvnb,8) * PowerSpectrum(k,totwvnb)
938 end do
939 do totwvnb = 0, nTrunc
940 if ( sum /= 0.0d0 ) then
941 NormPowerSpectrum(k,totwvnb) = PowerSpectrum(k,totwvnb) / sum
942 else
943 NormPowerSpectrum(k,totwvnb) = 0.d0
944 end if
945 end do
946 end do
947 !$OMP END PARALLEL DO
948
949 !- 1.2 Part 2
950 allocate( SpectralStateVar(lst_norm%nla,lst_norm%nphase,bhi%nVarLev) )
951 allocate( GridState(hco_bhi%ni, hco_bhi%nj, bhi%nVarLev) )
952
953 !- 1.2.1 Spectral transform of a delta function (at the center of the domain)
954 GridState(:,:,:) = 0.d0
955 GridState(hco_bhi%ni/2,hco_bhi%nj/2,:) = 1.d0
956
957 kind = 'GridPointToSpectral'
958 call lst_VarTransform(lst_norm, & ! IN
959 SpectralStateVar, & ! OUT
960 GridState, & ! IN
961 kind, bhi%nVarLev) ! IN
962
963 !- 1.2.2 Apply the horizontal correlation function
964 !$OMP PARALLEL DO PRIVATE (totwvnb,e,ila,p,k)
965 do totwvnb = 0, nTrunc
966 do e = 1, lst_norm%nePerK(totwvnb)
967 ila = lst_norm%ilaFromEK(e,totwvnb)
968 do p = 1, lst_norm%nphase
969 do k = 1, bhi%nVarLev
970 SpectralStateVar(ila,p,k) = SpectralStateVar(ila,p,k) * NormPowerSpectrum(k,totwvnb) * &
971 lst_norm%NormFactor(ila,p) * lst_norm%NormFactorAd(ila,p)
972 end do
973 end do
974 end do
975 end do
976 !$OMP END PARALLEL DO
977
978 !- 1.2.3 Move back to physical space
979 kind = 'SpectralToGridPoint'
980 call lst_VarTransform(lst_norm, & ! IN
981 SpectralStateVar, & ! IN
982 GridState, & ! OUT
983 kind, bhi%nVarLev) ! IN
984
985 !- 1.2.4 Normalize to 1
986 do k = 1, bhi%nVarLev
987 if ( GridState(hco_bhi%ni/2,hco_bhi%nj/2,k) < 0.d0 ) then
988 write(*,*) 'NormalizePowerSpectrum: Problem in normalization ', k, GridState(hco_bhi%ni/2,hco_bhi%nj/2,k)
989 call utl_abort('aborting in NormalizePowerSpectrum')
990 end if
991
992 if ( GridState(hco_bhi%ni/2,hco_bhi%nj/2,k) /= 0.d0 ) then
993 write(*,*) 'Normalization factor = ', k, GridState(hco_bhi%ni/2,hco_bhi%nj/2,k), 1.d0 / GridState(hco_bhi%ni/2,hco_bhi%nj/2,k)
994 NormPowerSpectrum(k,:) = NormPowerSpectrum(k,:) / GridState(hco_bhi%ni/2,hco_bhi%nj/2,k)
995 else
996 write(*,*) 'Setting NormPowerSpectrum to zero = ', k
997 NormPowerSpectrum(k,:) = 0.d0
998 end if
999 end do
1000
1001 deallocate(SpectralStateVar)
1002 deallocate(GridState)
1003
1004 write(*,*)
1005 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1006 write(*,*) 'NormalizePowerSpectrum: Done!'
1007
1008 end subroutine normalizePowerSpectrum
1009
1010 !--------------------------------------------------------------------------
1011 ! calcHorizScale
1012 !--------------------------------------------------------------------------
1013 subroutine calcHorizScale(HorizScale,SpCovariance)
1014 !
1015 !:Purpose: To compute horizontal lenght scales based on the power spectra
1016 !
1017 implicit none
1018
1019 ! Arguments:
1020 real(8), intent(out) :: HorizScale(bhi%nVarLev)
1021 real(8), intent(in) :: SpCovariance(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1022
1023 ! Locals:
1024 real(8) :: a, b, beta, dx, dist
1025 integer :: totwvnb, k, var
1026
1027 write(*,*)
1028 write(*,*) 'CalcHorizScale: Starting...'
1029
1030 !
1031 !- Computing distance-related variables
1032 !
1033
1034 ! Grid spacing in meters
1035 dx = hco_bhi%dlon * ec_ra
1036 write(*,*)
1037 write(*,*) 'grid spacing (m) =', dx
1038
1039 dist = max(hco_bhi%ni, hco_bhi %nj) * dx
1040 beta = (dist/(2.d0*MPC_PI_R8))**2
1041
1042 !
1043 !- Estimate horizontal correlation scales based on the power spectra
1044 !
1045 do k = 1, bhi%nVarLev
1046 a = 0.d0
1047 b = 0.d0
1048 do totwvnb = 0, nTrunc
1049 a = a + SpCovariance(k,k,totwvnb) * totwvnb
1050 b = b + SpCovariance(k,k,totwvnb) * totwvnb**3
1051 end do
1052 if (b <= 0.d0) then
1053 HorizScale(k) = 0.d0
1054 else
1055 HorizScale(k) = sqrt(2.d0*a*beta/b)
1056 end if
1057 end do
1058
1059 do var = 1, bhi%nControlVariable
1060 write(*,*)
1061 write(*,*) bhi%controlVariable(var)%nomvar(cv_bhi)
1062 do k = bhi%controlVariable(var)%varLevIndexStart, bhi%controlVariable(var)%varLevIndexEnd
1063 write(*,'(i3,2X,f9.2,2X,a2)') k, HorizScale(k)/1000.d0, 'km'
1064 end do
1065 end do
1066
1067 write(*,*)
1068 write(*,*) 'calcHorizScale: Done!'
1069
1070 end subroutine calcHorizScale
1071
1072 !--------------------------------------------------------------------------
1073 ! calcTotVertCorrel
1074 !--------------------------------------------------------------------------
1075 subroutine calcTotVertCorrel(TotVertCorrel, SpVertCorrel, PowerSpectrum)
1076 !
1077 !:Purpose: To compute the total vertical correlations (i.e. gridpoint equivalent)
1078 !
1079 implicit none
1080
1081 ! Arguments:
1082 real(8), intent(out) :: TotVertCorrel(bhi%nVarLev,bhi%nVarLev)
1083 real(8), intent(in) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1084 real(8), intent(in) :: PowerSpectrum(bhi%nVarLev,0:nTrunc)
1085
1086 ! Locals:
1087 real(8), allocatable :: TotVertCov(:,:)
1088 integer :: k1, k2, totwvnb
1089
1090 write(*,*)
1091 write(*,*) 'CalcTotVertCorrel: Starting...'
1092 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1093
1094 !
1095 !- 1. Calculate the total Normalized Covariance Matrix
1096 !
1097 allocate(TotVertCov(bhi%nVarLev,bhi%nVarLev))
1098 TotVertCov(:,:) = 0.d0
1099
1100 do k2 = 1, bhi%nVarLev
1101 do k1 = 1, bhi%nVarLev
1102 do totwvnb = 0, nTrunc
1103 TotVertCov(k1,k2) = TotVertCov(k1,k2) + &
1104 real(totwvnb,8) * SpVertCorrel(k1,k2,totwvnb) * &
1105 sqrt(PowerSpectrum(k1,totwvnb) * PowerSpectrum(k2,totwvnb))
1106 end do
1107 end do
1108 end do
1109
1110 !
1111 !- 2. Transform into correlations
1112 !
1113 do k2 = 1, bhi%nVarLev
1114 do k1 = 1, bhi%nVarLev
1115 if ( TotVertCov(k1,k1) /= 0.d0 .and. &
1116 TotVertCov(k2,k2) /= 0.d0 ) then
1117 TotVertCorrel(k1,k2) = TotVertCov(k1,k2) / &
1118 ( sqrt(TotVertCov(k1,k1)) * sqrt(TotVertCov(k2,k2)) )
1119 else
1120 TotVertCorrel(k1,k2) = 0.d0
1121 end if
1122 end do
1123 end do
1124
1125 deallocate(TotVertCov)
1126
1127 write(*,*)
1128 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1129 write(*,*) 'CalcTotVertCorrel: Done!'
1130
1131 end subroutine calcTotVertCorrel
1132
1133 !--------------------------------------------------------------------------
1134 ! calcBsqrt
1135 !--------------------------------------------------------------------------
1136 subroutine calcBsqrt(Bsqrt,B)
1137 !
1138 !:Purpose: To compute the sqare-root of B
1139 !
1140 implicit none
1141
1142 ! Arguments:
1143 real(8), intent(out) :: Bsqrt(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1144 real(8), intent(in) :: B (bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1145
1146 ! Locals:
1147 integer :: totwvnb
1148
1149 !
1150 !- Calculate B^0.5 for each total wave number
1151 !
1152 Bsqrt(:,:,:) = B(:,:,:)
1153
1154 do totwvnb = 0, nTrunc
1155 call utl_matSqrt(Bsqrt(:,:,totwvnb),bhi%nVarLev,1.0d0,.true.)
1156 end do
1157
1158 end subroutine calcBsqrt
1159
1160 !--------------------------------------------------------------------------
1161 ! setSpVertCorrel
1162 !--------------------------------------------------------------------------
1163 subroutine setSpVertCorrel(SpVertCorrel)
1164 !
1165 !:Purpose: To discard some user-defined cross-correlations
1166 !
1167 implicit none
1168
1169 ! Arguments:
1170 real(8), intent(inout) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1171
1172 ! Locals:
1173 real(8), allocatable :: KeepOrDiscard(:,:)
1174 integer :: totwvnb, var1, var2, k1, k2
1175
1176 write(*,*)
1177 write(*,*) 'SetSpVertCorrel: Starting...'
1178
1179 !
1180 !- Determine which bloc of the correlation matrix to keep/discard
1181 !
1182 allocate(KeepOrDiscard(bhi%nControlVariable,bhi%nControlVariable))
1183
1184 !- Calculate upper half
1185 do var2 = 1, bhi%nControlVariable
1186 do var1 = 1, var2
1187 if (var1 == var2) then
1188 KeepOrDiscard(var1,var2) = 1.d0 ! Keep the Auto-Correlations
1189 elseif( any(correlatedVariables == bhi%controlVariable(var1)%nomvar(cv_bhi)) .and. &
1190 any(correlatedVariables == bhi%controlVariable(var2)%nomvar(cv_bhi)) ) then
1191 KeepOrDiscard(var1,var2) = 1.d0 ! Keep these Cross-Correlations
1192 else
1193 KeepOrDiscard(var1,var2) = 0.d0 ! Discard these Cross-Correlations
1194 end if
1195 write(*,*) var1, var2, bhi%controlVariable(var1)%nomvar(cv_bhi), bhi%controlVariable(var2)%nomvar(cv_bhi), KeepOrDiscard(var1,var2)
1196 end do
1197 end do
1198
1199 ! Symmetrize
1200 do var1 = 2, bhi%nControlVariable
1201 do var2 = 1, var1-1
1202 KeepOrDiscard(var1,var2) = KeepOrDiscard(var2,var1)
1203 end do
1204 end do
1205
1206 !
1207 !- Modify the Vertical Correlation Matrix
1208 !
1209 do totwvnb = 0, nTrunc
1210
1211 do var2 = 1, bhi%nControlVariable
1212 do var1 = 1, bhi%nControlVariable
1213
1214 do k2 = bhi%controlVariable(var2)%varLevIndexStart, bhi%controlVariable(var2)%varLevIndexEnd
1215 do k1 = bhi%controlVariable(var1)%varLevIndexStart, bhi%controlVariable(var1)%varLevIndexEnd
1216 SpVertCorrel(k1,k2,totwvnb) = KeepOrDiscard(var1,var2) * &
1217 SpVertCorrel(k1,k2,totwvnb)
1218 end do
1219 end do
1220
1221 end do
1222 end do
1223
1224 end do ! total wave number
1225
1226 deallocate(KeepOrDiscard)
1227
1228 write(*,*)
1229 write(*,*) 'SetSpVertCorrel: Done!'
1230
1231 end subroutine setSpVertCorrel
1232
1233 !--------------------------------------------------------------------------
1234 ! calcVertCorrel
1235 !--------------------------------------------------------------------------
1236 subroutine calcVertCorrel(ensPerts)
1237 !
1238 !:Purpose: To compute vertical correlations from an ensemble of gridded data
1239 !
1240 implicit none
1241
1242 ! Arguments:
1243 type(struct_ens), intent(inout) :: ensPerts
1244
1245 ! Locals:
1246 real(8), allocatable :: vertCorrel(:,:)
1247 real(4), pointer :: ptr4d_k1_r4(:,:,:,:)
1248 real(4), pointer :: ptr4d_k2_r4(:,:,:,:)
1249 real(8), allocatable :: vertCorrel_local(:,:)
1250 integer :: lonIndex, latIndex, k1, k2, memberIndex
1251 integer :: myLonBeg, myLonEnd, myLatBeg, myLatEnd, nSize, ier
1252 integer :: fstouv, fnom, fstfrm, fclos, iunstats
1253
1254 write(*,*)
1255 write(*,*) 'CalcVertCorrel: Starting...'
1256 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1257
1258 allocate(vertCorrel(bhi%nVarLev,bhi%nVarLev))
1259 vertCorrel(:,:) = 0.d0
1260
1261 allocate(vertCorrel_local(bhi%nVarLev,bhi%nVarLev))
1262 vertCorrel_local(:,:) = 0.d0
1263
1264 !
1265 !- Calculate the Vertical Correlation in GridPoint Space
1266 ! ... we assume that the ensemble grid point mean was removed and that
1267 ! the ensemble values were divided by the grid point std dev.
1268
1269 call ens_getLatLonBounds(ensPerts, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
1270
1271 do memberIndex = 1, nEns
1272
1273 !$OMP PARALLEL DO PRIVATE (k1,k2,ptr4d_k1_r4,ptr4d_k2_r4,latIndex,lonIndex)
1274 do k2 = 1, bhi%nVarLev
1275 do k1 = 1, bhi%nVarLev
1276
1277 ptr4d_k1_r4 => ens_getOneLev_r4(ensPerts,k1)
1278 ptr4d_k2_r4 => ens_getOneLev_r4(ensPerts,k2)
1279
1280 do latIndex = myLatBeg, myLatEnd
1281 do lonIndex = myLonBeg, myLonEnd
1282
1283 if (lonIndex > hco_ens%ni .or. latIndex > hco_ens%nj) cycle ! do not use data in the extension zone
1284
1285 VertCorrel_local(k1,k2) = VertCorrel_local(k1,k2) &
1286 + real(ptr4d_k1_r4(memberIndex,1,lonIndex,latIndex),8) &
1287 * real(ptr4d_k2_r4(memberIndex,1,lonIndex,latIndex),8)
1288 end do
1289 end do
1290
1291 end do
1292 end do
1293 !$OMP END PARALLEL DO
1294
1295 end do ! Loop in Ensemble
1296
1297 !- Communication
1298 nSize = bhi%nVarLev * bhi%nVarLev
1299 call rpn_comm_reduce(vertCorrel_local,vertCorrel,nsize,"mpi_double_precision","mpi_sum",0,"GRID",ier)
1300
1301 deallocate(vertCorrel_local)
1302
1303 !- Conversion to correlation
1304 if (mmpi_myid == 0) then
1305 vertCorrel(:,:) = vertCorrel(:,:) / real((nEns-1)*hco_ens%nj*hco_ens%ni,8)
1306 end if
1307
1308 !- Output
1309 if (mmpi_myid == 0) then
1310 iunstats = 0
1311 ier = fnom(iunstats,'./vertCorrel.fst','RND',0)
1312 ier = fstouv(iunstats,'RND')
1313 call WriteTotVertCorrel(VertCorrel,iunstats,'ZT','GPVCOR_CORE') ! IN
1314 ier = fstfrm(iunstats)
1315 ier = fclos (iunstats)
1316 end if
1317
1318 deallocate(VertCorrel)
1319
1320 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1321 write(*,*) 'CalcVertCorrel: Done!'
1322
1323 end subroutine calcVertCorrel
1324
1325 !--------------------------------------------------------------------------
1326 ! horizCorrelFunction
1327 !--------------------------------------------------------------------------
1328 subroutine horizCorrelFunction(NormB)
1329 !
1330 !:Purpose: To compute and write the horizontal correlation function of
1331 ! every variables and levels in the correlation formulation of
1332 ! the B matrix (i.e. C matrix)
1333 !
1334 implicit none
1335
1336 ! Arguments:
1337 real(8), intent(in) :: NormB(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1338
1339 ! Locals:
1340 real(8), allocatable :: SpectralStateVar(:,:,:)
1341 real(8), allocatable :: GridState(:,:,:)
1342 type(struct_gsv) :: statevector
1343 real(8), pointer :: ptr3d_r8(:,:,:)
1344 character(len=4), pointer :: varNamesList(:)
1345 integer :: e, ila, p, k, totwvnb
1346 type(struct_lst) :: lst_cor ! Spectral transform Parameters
1347 character(len=24) :: kind
1348
1349 nullify(varNamesList)
1350 call ens_varNamesList(varNamesList,ensPerts)
1351 call gsv_allocate(statevector, ens_getNumStep(ensPerts), &
1352 hco_bhi, ens_getVco(ensPerts), varNames_opt=varNamesList, &
1353 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.false., &
1354 dataKind_opt=8 )
1355
1356 !
1357 !- 1. Setup the (LAM) spectral transform
1358 !
1359 call lst_setup(lst_cor, & ! OUT
1360 hco_bhi%ni, hco_bhi%nj, & ! IN
1361 hco_bhi%dlon, nTrunc, & ! IN
1362 'NoMpi')
1363
1364 allocate( SpectralStateVar(lst_cor%nla,lst_cor%nphase,bhi%nVarLev) )
1365 allocate( GridState(hco_bhi%ni, hco_bhi%nj, bhi%nVarLev) )
1366
1367 !- 3.2.1 Spectral transform of a delta function (at the center of the domain)
1368 GridState(:,:,:) = 0.d0
1369 GridState(hco_bhi%ni/2,hco_bhi%nj/2,:) = 1.d0
1370
1371 kind = 'GridPointToSpectral'
1372 call lst_VarTransform(lst_cor, & ! IN
1373 SpectralStateVar, & ! OUT
1374 GridState, & ! IN
1375 kind, bhi%nVarLev) ! IN
1376
1377 !- 3.2.2 Apply the horizontal correlation function
1378 !$OMP PARALLEL DO PRIVATE (totwvnb,e,ila,p,k)
1379 do totwvnb = 0, nTrunc
1380 do e = 1, lst_cor%nePerK(totwvnb)
1381 ila = lst_cor%ilaFromEK(e,totwvnb)
1382 do p = 1, lst_cor%nphase
1383 do k = 1, bhi%nVarLev
1384 SpectralStateVar(ila,p,k) = SpectralStateVar(ila,p,k) * NormB(k,k,totwvnb) * &
1385 lst_cor%NormFactor(ila,p) * lst_cor%NormFactorAd(ila,p)
1386 end do
1387 end do
1388 end do
1389 end do
1390 !$OMP END PARALLEL DO
1391
1392 !- 3.2.3 Move back to physical space
1393 kind = 'SpectralToGridPoint'
1394 call lst_VarTransform(lst_cor, & ! IN
1395 SpectralStateVar, & ! IN
1396 GridState, & ! OUT
1397 kind, bhi%nVarLev) ! IN
1398
1399 !
1400 !- 4. Write to file
1401 !
1402 call gsv_getField(statevector,ptr3d_r8)
1403 ptr3d_r8(:,:,:) = GridState(:,:,:)
1404 call gio_writeToFile(statevector, './horizCorrel.fst', 'CORRELFUNC')
1405
1406 deallocate(SpectralStateVar)
1407 deallocate(GridState)
1408 call gsv_deallocate(statevector)
1409
1410 end subroutine horizCorrelFunction
1411
1412 !--------------------------------------------------------------------------
1413 ! applyHorizLoc
1414 !--------------------------------------------------------------------------
1415 subroutine applyHorizLoc(NormPowerSpectrum)
1416 !
1417 !:Purpose: To apply horizontal localization to the spectral correlations
1418 !
1419 implicit none
1420
1421 ! Arguments:
1422 real(8), intent(inout) :: NormPowerSpectrum(bhi%nVarLev,0:nTrunc)
1423
1424 ! Locals:
1425 real(8), allocatable :: SpectralStateVar(:,:,:)
1426 real(8), allocatable :: GridState(:,:,:)
1427 real(8), allocatable :: GridStateLoc(:,:,:)
1428 real(8), allocatable :: PowerSpectrum(:,:)
1429 real(8), allocatable :: SumWeight(:)
1430 real(8), allocatable :: local_length(:)
1431 type(struct_lst) :: lst_hloc ! Spectral transform Parameters
1432 integer :: totwvnb, var, k, e, ila, p
1433 real(8) :: hlocalize
1434 character(len=24) :: kind
1435
1436 write(*,*)
1437 write(*,*) 'applyHorizLoc: Starting...'
1438
1439 !
1440 !- 1. Setup the (LAM) spectral transform
1441 !
1442 call lst_setup(lst_hloc, & ! OUT
1443 hco_bhi%ni, hco_bhi%nj, & ! IN
1444 hco_bhi%dlon, nTrunc, & ! IN
1445 'NoMpi')
1446
1447 !
1448 !- 1. Get the original gridpoint horizontal correlations
1449 !
1450
1451 allocate( SpectralStateVar(lst_hloc%nla,lst_hloc%nphase,bhi%nVarLev) )
1452 allocate( GridState(hco_bhi%ni, hco_bhi%nj, bhi%nVarLev) )
1453
1454 !- 1.1 Spectral transform of a delta function (at the lower-left of the domain)
1455 GridState(:,:,:) = 0.d0
1456 GridState(1,1,:) = 1.d0
1457
1458 kind = 'GridPointToSpectral'
1459 call lst_VarTransform(lst_hloc, & ! IN
1460 SpectralStateVar, & ! OUT
1461 GridState, & ! IN
1462 kind, bhi%nVarLev) ! IN
1463
1464 !- 1.2 Apply the horizontal correlation function
1465 !$OMP PARALLEL DO PRIVATE (totwvnb,e,ila,p,k)
1466 do totwvnb = 0, nTrunc
1467 do e = 1, lst_hloc%nePerK(totwvnb)
1468 ila = lst_hloc%ilaFromEK(e,totwvnb)
1469 do p = 1, lst_hloc%nphase
1470 do k = 1, bhi%nVarLev
1471 SpectralStateVar(ila,p,k) = SpectralStateVar(ila,p,k) * NormPowerSpectrum(k,totwvnb) * &
1472 lst_hloc%NormFactor(ila,p) * lst_hloc%NormFactorAd(ila,p)
1473 end do
1474 end do
1475 end do
1476 end do
1477 !$OMP END PARALLEL DO
1478
1479 !- 1.3 Move back to physical space
1480 kind = 'SpectralToGridPoint'
1481 call lst_VarTransform(lst_hloc, & ! IN
1482 SpectralStateVar, & ! IN
1483 GridState, & ! OUT
1484 kind, bhi%nVarLev) ! IN
1485
1486 !
1487 !- 2. Create and apply the localization function in gridpoint space
1488 !
1489 allocate(local_length(bhi%nVarLev))
1490 allocate(GridStateLoc(hco_bhi%ni, hco_bhi%nj, bhi%nVarLev) )
1491
1492 do var = 1, bhi%nControlVariable
1493 if ( bhi%controlVariable(var)%nomvar(cv_bhi) == 'CC' .or. &
1494 bhi%controlVariable(var)%nomvar(cv_bhi) == 'PP' .or. &
1495 bhi%controlVariable(var)%nomvar(cv_bhi) == 'QR' .or. &
1496 bhi%controlVariable(var)%nomvar(cv_bhi) == 'DD' ) then
1497 hLocalize = hLocalize_wind
1498 else if ( bhi%controlVariable(var)%nomvar(cv_bhi) == 'TT' .or. &
1499 bhi%controlVariable(var)%nomvar(cv_bhi) == 'TG' .or. &
1500 bhi%controlVariable(var)%nomvar(cv_bhi) == 'P0' ) then
1501 hLocalize = hLocalize_mass
1502 else if ( bhi%controlVariable(var)%nomvar(cv_bhi) == 'LQ' ) then
1503 hLocalize = hLocalize_humidity
1504 else
1505 hLocalize = hLocalize_other
1506 end if
1507
1508 do k = bhi%controlVariable(var)%varLevIndexStart, bhi%controlVariable(var)%varLevIndexEnd
1509 local_length(k) = hLocalize
1510 end do
1511
1512 end do
1513
1514 call lfn_setup('FifthOrder') ! IN
1515 call lfn_CreateBiPerFunction( GridStateLoc, & ! OUT
1516 local_length, hco_bhi%dlon, & ! IN
1517 hco_bhi%ni, hco_bhi%nj, bhi%nVarLev) ! IN
1518
1519 GridState(:,:,:) = GridStateLoc(:,:,:) * GridState(:,:,:)
1520
1521 deallocate(GridStateLoc)
1522 deallocate(local_length)
1523
1524 !
1525 !- 3. Create the localized normalized power spectrum
1526 !
1527 allocate(PowerSpectrum(bhi%nVarLev,0:nTrunc))
1528
1529 !- 3.1 Transform to spectral space
1530 kind = 'GridPointToSpectral'
1531 call lst_VarTransform(lst_hloc, & ! IN
1532 SpectralStateVar, & ! OUT
1533 GridState, & ! IN
1534 kind, bhi%nVarLev) ! IN
1535
1536 !- 3.2 Compute band mean
1537 allocate(SumWeight(0:nTrunc))
1538 SumWeight(:) = 0.d0
1539
1540 PowerSpectrum(:,:) = 0.d0
1541 do totwvnb = 0, nTrunc
1542 do e = 1, lst_hloc%nePerK(totwvnb)
1543 ila = lst_hloc%ilaFromEK(e,totwvnb)
1544 do p = 1, lst_hloc%nphase
1545 SumWeight(totwvnb) = SumWeight(totwvnb) + lst_hloc%Weight(ila)
1546 do k = 1, bhi%nVarLev
1547 PowerSpectrum(k,totwvnb) = PowerSpectrum(k,totwvnb) + &
1548 lst_hloc%Weight(ila) * abs(SpectralStateVar(ila,p,k))
1549 end do
1550 end do
1551 end do
1552 end do
1553
1554 do totwvnb = 0, nTrunc
1555 if (SumWeight(totwvnb) /= 0.d0) then
1556 PowerSpectrum(:,totwvnb) = PowerSpectrum(:,totwvnb) / SumWeight(totwvnb)
1557 else
1558 PowerSpectrum(:,totwvnb) = 0.d0
1559 endif
1560 end do
1561
1562 deallocate(SumWeight)
1563
1564 !- 3.3 Normalize
1565 call normalizePowerSpectrum(PowerSpectrum, & ! IN
1566 NormPowerSpectrum) ! OUT
1567
1568 deallocate(SpectralStateVar)
1569 deallocate(GridState)
1570 deallocate(PowerSpectrum)
1571
1572 write(*,*)
1573 write(*,*) 'applyHorizLoc: Done!'
1574
1575 end subroutine applyHorizLoc
1576
1577 !--------------------------------------------------------------------------
1578 ! applyVertLoc
1579 !--------------------------------------------------------------------------
1580 subroutine applyVertLoc(SpVertCorrel)
1581 !
1582 !:Purpose: To apply vertical localization to the spectral correlations
1583 !
1584 implicit none
1585
1586 ! Arguments:
1587 real(8), intent(inout) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1588
1589 ! Locals:
1590 integer :: totwvnb, var1, var2, k1, k2, lev1, lev2
1591 real(8) :: dist, fact, vLocalize, pres1, pres2, vLocalize1, vLocalize2
1592
1593 write(*,*)
1594 write(*,*) 'applyVertLoc: Starting...'
1595
1596 !
1597 !- 1. Select the localization function
1598 !
1599 call lfn_setup('FifthOrder') ! IN
1600
1601 !
1602 !- 2. Apply localization to the spectral vertical correlations
1603 !
1604
1605 !- 2.1 Loop on control variables
1606 do var2 = 1, bhi%nControlVariable
1607 do var1 = 1, bhi%nControlVariable
1608
1609 !- 2.2 Set the vertical length scale
1610
1611 !- 2.2.1 Select the length scale for control variable 1
1612 if ( bhi%controlVariable(var1)%nomvar(cv_bhi) == 'CC' .or. &
1613 bhi%controlVariable(var1)%nomvar(cv_bhi) == 'PP' .or. &
1614 bhi%controlVariable(var1)%nomvar(cv_bhi) == 'QR' .or. &
1615 bhi%controlVariable(var1)%nomvar(cv_bhi) == 'DD') then
1616 vLocalize1 = vLocalize_wind
1617 else if ( bhi%controlVariable(var1)%nomvar(cv_bhi) == 'TT' .or. &
1618 bhi%controlVariable(var1)%nomvar(cv_bhi) == 'TG' .or. &
1619 bhi%controlVariable(var1)%nomvar(cv_bhi) == 'P0' ) then
1620 vLocalize1 = vLocalize_mass
1621 else if ( bhi%controlVariable(var1)%nomvar(cv_bhi) == 'LQ' ) then
1622 vLocalize1 = vLocalize_humidity
1623 else
1624 vLocalize1 = vLocalize_other
1625 end if
1626
1627 !- 2.2.2 Select the length scale for control variable 2
1628 if ( bhi%controlVariable(var2)%nomvar(cv_bhi) == 'CC' .or. &
1629 bhi%controlVariable(var2)%nomvar(cv_bhi) == 'PP' .or. &
1630 bhi%controlVariable(var2)%nomvar(cv_bhi) == 'QR' .or. &
1631 bhi%controlVariable(var2)%nomvar(cv_bhi) == 'DD') then
1632 vLocalize2 = vLocalize_wind
1633 else if ( bhi%controlVariable(var2)%nomvar(cv_bhi) == 'TT' .or. &
1634 bhi%controlVariable(var2)%nomvar(cv_bhi) == 'TG' .or. &
1635 bhi%controlVariable(var2)%nomvar(cv_bhi) == 'P0' ) then
1636 vLocalize2 = vLocalize_mass
1637 else if ( bhi%controlVariable(var2)%nomvar(cv_bhi) == 'LQ' ) then
1638 vLocalize2 = vLocalize_humidity
1639 else
1640 vLocalize2 = vLocalize_other
1641 end if
1642
1643 !- 2.2.3 Length scale to be use for var1-var2 correlation
1644 vLocalize = (vLocalize1+vLocalize2)/2.d0
1645
1646 !- 2.3 Loop on vertical levels
1647 do k2 = bhi%controlVariable(var2)%varLevIndexStart, bhi%controlVariable(var2)%varLevIndexEnd
1648 do k1 = bhi%controlVariable(var1)%varLevIndexStart, bhi%controlVariable(var1)%varLevIndexEnd
1649
1650 !- 2.4 Set the pressure values
1651
1652 !- 2.4.1 Pressure for control-variable-1 level
1653 if (bhi%controlVariable(var1)%nlev /= 1) then ! variable 3D
1654 lev1 = k1 - bhi%controlVariable(var1)%varLevIndexStart + 1
1655 if (bhi%controlVariable(var1)%GridType == 'TH') then
1656 pres1 = pressureProfile_T(lev1)
1657 else
1658 pres1 = pressureProfile_M(lev1)
1659 end if
1660 else
1661 pres1 = pressureProfile_M(vco_bhi%nlev_M) ! variable 2D
1662 end if
1663
1664 !- 2.4.2 Pressure for control-variable-2 level
1665 if (bhi%controlVariable(var2)%nlev /= 1) then ! variable 3D
1666 lev2 = k2 - bhi%controlVariable(var2)%varLevIndexStart + 1
1667 if (bhi%controlVariable(var2)%GridType == 'TH') then
1668 pres2 = pressureProfile_T(lev2)
1669 else
1670 pres2 = pressureProfile_M(lev2)
1671 end if
1672 else
1673 pres2 = pressureProfile_M(vco_bhi%nlev_M) ! variable 2D
1674 end if
1675
1676 !- 2.5 Compute the localization factor
1677 dist = abs(log(pres2) - log(pres1))
1678 fact = lfn_response(dist,vLocalize)
1679
1680 !- 2.6 Localize each total wavenumber (not scale-dependent!)
1681 do totwvnb = 0, nTrunc
1682 SpVertCorrel(k1,k2,totwvnb) = fact * SpVertCorrel(k1,k2,totwvnb)
1683 end do
1684
1685 end do
1686 end do
1687
1688 end do
1689 end do
1690
1691 deallocate(pressureProfile_M)
1692 deallocate(pressureProfile_T)
1693
1694 write(*,*)
1695 write(*,*) 'applyVertLoc: Done!'
1696
1697 end subroutine applyVertLoc
1698
1699 !--------------------------------------------------------------------------
1700 ! scaleStdDev
1701 !--------------------------------------------------------------------------
1702 subroutine scaleStdDev(statevector_stdDev)
1703 !
1704 !:Purpose: To scale the gridpoint background-error standard deviations
1705 !
1706 implicit none
1707
1708 ! Arguments:
1709 type(struct_gsv), intent(inout) :: statevector_stdDev
1710
1711 ! Locals:
1712 real(8), pointer :: ptr3d_r8(:,:,:)
1713 real(8) :: multFactor
1714 integer :: nVarLev, varLevIndex, levIndex
1715 character(len=4) :: varName
1716
1717 write(*,*)
1718 write(*,*) 'scaleStdDev: Starting...'
1719
1720 nVarLev = gsv_getNumK(statevector_stdDev)
1721
1722 call gsv_getField(statevector_stdDev,ptr3d_r8)
1723
1724 do varLevIndex = 1, nVarLev
1725 varName = gsv_getVarNameFromK(statevector_stdDev,varLevIndex)
1726 levIndex = gsv_getLevFromK(statevector_stdDev,varLevIndex)
1727
1728 if ( vnl_varLevelFromVarname(varName) == 'MM' ) then
1729 multFactor = scaleFactor_M(levIndex)
1730 else if ( vnl_varLevelFromVarname(varName) == 'TH' ) then
1731 multFactor = scaleFactor_T(levIndex)
1732 else ! SF
1733 multFactor = scaleFactor_SF
1734 end if
1735
1736 ptr3d_r8(:,:,varLevIndex) = multFactor * ptr3d_r8(:,:,varLevIndex)
1737 end do
1738
1739 write(*,*)
1740 write(*,*) 'scaleStdDev: Done...'
1741
1742 end subroutine scaleStdDev
1743
1744 !--------------------------------------------------------------------------
1745 ! writeVarStats
1746 !--------------------------------------------------------------------------
1747 subroutine writeVarStats(Bsqrt,statevector_stdDev)
1748 !
1749 !:Purpose: To write data needed for VAR applications, i.e C^1/2 and stdDev
1750 !
1751 implicit none
1752
1753 ! Arguments:
1754 real(8), intent(in) :: Bsqrt(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1755 type(struct_gsv), intent(in) :: statevector_stdDev
1756
1757 ! Locals:
1758 integer :: ier, fstouv, fnom, fstfrm, fclos
1759 integer :: iunstats
1760 character(len=24) :: fileName = './bgcov.fst'
1761
1762 write(*,*)
1763 write(*,*) 'Writing covariance statistics for VAR'
1764
1765 !
1766 !- 1. Add the gridded std dev (and the Tic-Tac-Toc)
1767 !
1768 call gio_writeToFile(statevector_stdDev, trim(fileName), 'STDDEV', &
1769 typvar_opt = 'E', numBits_opt = 32)
1770
1771 !
1772 !- 2. Add C^1/2
1773 !
1774 if (mmpi_myid == 0) then
1775 !- Opening Output file
1776 iunstats = 0
1777 ier = fnom(iunstats,trim(fileName),'RND',0)
1778 ier = fstouv(iunstats,'RND')
1779
1780 !- Add Control Variable Info
1781 call WriteControlVarInfo(iunstats)
1782
1783 !- Bsqrt
1784 call WriteSpVertCorrel(Bsqrt,iunstats,'ZN','B_SQUAREROOT') ! IN
1785
1786 !- Closing output file
1787 ier = fstfrm(iunstats)
1788 ier = fclos (iunstats)
1789 end if
1790
1791 end subroutine writeVarStats
1792
1793 !--------------------------------------------------------------------------
1794 ! writeDiagStats
1795 !--------------------------------------------------------------------------
1796 subroutine writeDiagStats(NormB,SpVertCorrel,TotVertCorrel,statevector_mean, &
1797 statevector_stdDevGridPoint,PowerSpectrum,HorizScale)
1798 !
1799 !:Purpose: To write other relevant data computed during the
1800 ! calculation of Bhi that are not needed for VAR applications
1801 !
1802 implicit none
1803
1804 ! Arguments:
1805 real(8), intent(in) :: NormB(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1806 real(8), intent(in) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1807 real(8), intent(in) :: TotVertCorrel(bhi%nVarLev,bhi%nVarLev)
1808 type(struct_gsv), intent(in) :: statevector_mean
1809 type(struct_gsv), intent(in) :: statevector_stdDevGridPoint
1810 real(8), intent(in) :: PowerSpectrum(bhi%nVarLev,0:nTrunc)
1811 real(8), intent(in) :: HorizScale(bhi%nVarLev)
1812
1813 ! Locals:
1814 integer :: ier, fstouv, fnom, fstfrm, fclos
1815 integer :: iunstats
1816 character(len=24) :: fileName = './bgcov_diag.fst'
1817
1818 write(*,*)
1819 write(*,*) 'Writing Diagnostics'
1820
1821 !
1822 !- 1. Add the gridded mean and std dev (and the Tic-Tac-Toc)
1823 !
1824 call gio_writeToFile(statevector_mean, trim(fileName), 'ENSMEAN', &
1825 typvar_opt = 'E', numBits_opt = 32)
1826 call gio_writeToFile(statevector_stdDevGridPoint, trim(fileName), 'STDDEV_GRIDP', &
1827 typvar_opt = 'E', numBits_opt = 32)
1828
1829 !
1830 !- 2. Add stats in spectral space
1831 !
1832 if (mmpi_myid == 0) then
1833 !- Opening Output file
1834 iunstats = 0
1835 ier = fnom(iunstats,trim(fileName),'RND',0)
1836 ier = fstouv(iunstats,'RND')
1837
1838 !- Spectral Vertical Correlations
1839 call writeSpVertCorrel(SpVertCorrel,iunstats,'ZZ','SPVERTCORREL') ! IN
1840
1841 !- Total Vertical Correlations
1842 call writeTotVertCorrel(TotVertCorrel,iunstats,'ZT','TTVERTCORREL') ! IN
1843
1844 !- Normalized Vertical Correlations
1845 call writeSpVertCorrel(NormB,iunstats,'ZN','NRVERTCORREL') ! IN
1846
1847 !- Power Spectrum
1848 call writePowerSpectrum(PowerSpectrum,iunstats,'POWERSPECT',cv_bhi) ! IN
1849
1850 !- Horizontal Correlation Length scale
1851 call writeHorizScale(HorizScale,iunstats,'HORIZSCALE',cv_bhi) ! IN
1852
1853 !- Closing output file
1854 ier = fstfrm(iunstats)
1855 ier = fclos (iunstats)
1856 end if
1857
1858 end subroutine writeDiagStats
1859
1860 !--------------------------------------------------------------------------
1861 ! writeSpVertCorrel
1862 !--------------------------------------------------------------------------
1863 subroutine writeSpVertCorrel(SpVertCorrel,iun,nomvar_in,etiket_in)
1864 !
1865 !:Purpose: To write vertical correlations in spectral space
1866 !
1867 implicit none
1868
1869 ! Arguments:
1870 real(8), intent(in) :: SpVertCorrel(bhi%nVarLev,bhi%nVarLev,0:nTrunc)
1871 integer, intent(in) :: iun
1872 character(len=*), intent(in) :: nomvar_in
1873 character(len=*), intent(in) :: etiket_in
1874
1875 ! Locals:
1876 real(4), allocatable :: work2d(:,:)
1877 real(4) :: work
1878 integer :: ier, fstecr, totwvnb
1879 integer :: dateo, npak, ni, nj, nk
1880 integer :: ip1, ip2, ip3, deet, npas, datyp
1881 integer :: ig1 ,ig2 ,ig3 ,ig4
1882 character(len=1) :: grtyp
1883 character(len=4) :: nomvar
1884 character(len=2) :: typvar
1885 character(len=12) :: etiket
1886
1887 allocate(work2d(bhi%nVarLev, bhi%nVarLev))
1888
1889 !- Loop over Total Wavenumbers
1890 do totwvnb = 0, nTrunc
1891
1892 npak = -32
1893 dateo = 0
1894 deet = 0
1895 npas = 0
1896 ni = bhi%nVarLev
1897 nj = bhi%nVarLev
1898 nk = 1
1899 ip1 = 0
1900 ip2 = totwvnb
1901 ip3 = nEns
1902 typvar = 'XX'
1903 nomvar = nomvar_in
1904 etiket = etiket_in
1905 grtyp = 'X'
1906 ig1 = 0
1907 ig2 = 0
1908 ig3 = 0
1909 ig4 = 0
1910 datyp = 5
1911
1912 !- Extract from full Matrix
1913 work2d(:,:) = real(SpVertCorrel(:,:,totwvnb),4)
1914
1915 !- Writing
1916 ier = fstecr(work2d, work, npak, iun, dateo, deet, npas, ni, nj, &
1917 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
1918 ig1, ig2, ig3, ig4, datyp, .true.)
1919
1920 end do ! Total Wavenumbers
1921
1922 deallocate(work2d)
1923
1924 end subroutine writeSpVertCorrel
1925
1926 !--------------------------------------------------------------------------
1927 ! writeTotVertCorrel
1928 !--------------------------------------------------------------------------
1929 subroutine writeTotVertCorrel(TotVertCorrel,iun,nomvar_in,etiket_in)
1930 !
1931 !:Purpose: To write the total vertical correlations
1932 ! (i.e. gridpoint equivalent)
1933 !
1934 implicit none
1935
1936 ! Arguments:
1937 real(8), intent(in) :: TotVertCorrel(bhi%nVarLev,bhi%nVarLev)
1938 integer, intent(in) :: iun
1939 character(len=*), intent(in) :: nomvar_in
1940 character(len=*), intent(in) :: etiket_in
1941
1942 ! Locals:
1943 real(4), allocatable :: workecr(:,:)
1944 real(4) :: work
1945 integer :: ier, fstecr
1946 integer :: dateo, npak, ni, nj, nk
1947 integer :: ip1, ip2, ip3, deet, npas, datyp
1948 integer :: ig1 ,ig2 ,ig3 ,ig4
1949 character(len=1) :: grtyp
1950 character(len=4) :: nomvar
1951 character(len=2) :: typvar
1952 character(len=12) :: etiket
1953
1954 allocate(workecr(bhi%nVarLev, bhi%nVarLev))
1955
1956 npak = -32
1957 dateo = 0
1958 deet = 0
1959 npas = 0
1960 ni = bhi%nVarLev
1961 nj = bhi%nVarLev
1962 nk = 1
1963 ip1 = 0
1964 ip2 = 0
1965 ip3 = nEns
1966 typvar = 'XX'
1967 nomvar = nomvar_in
1968 etiket = etiket_in
1969 grtyp = 'X'
1970 ig1 = 0
1971 ig2 = 0
1972 ig3 = 0
1973 ig4 = 0
1974 datyp = 5
1975
1976 !- Covert to real 4
1977 workecr(:,:) = real(TotVertCorrel(:,:),4)
1978
1979 !- Writing
1980 ier = fstecr(workecr, work, npak, iun, dateo, deet, npas, ni, nj, &
1981 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
1982 ig1, ig2, ig3, ig4, datyp, .true.)
1983
1984 deallocate(workecr)
1985
1986 end subroutine writeTotVertCorrel
1987
1988 !--------------------------------------------------------------------------
1989 ! writePowerSpectrum
1990 !--------------------------------------------------------------------------
1991 subroutine writePowerSpectrum(PowerSpectrum,iun,etiket_in,cv_type)
1992 !
1993 !:Purpose: To write the power spectrum
1994 !
1995 implicit none
1996
1997 ! Arguments:
1998 real(8), intent(in) :: PowerSpectrum(bhi%nVarLev,0:nTrunc)
1999 integer, intent(in) :: iun
2000 integer, intent(in) :: cv_type
2001 character(len=*), intent(in) :: Etiket_in
2002
2003 ! Locals:
2004 real(4), allocatable :: workecr(:,:)
2005 real(4) :: work
2006 integer :: ier, fstecr
2007 integer :: var, k, kgdim
2008 integer :: dateo, npak, ni, nj, nk
2009 integer :: ip1, ip2, ip3, deet, npas, datyp
2010 integer :: ig1 ,ig2 ,ig3 ,ig4
2011 character(len=1) :: grtyp
2012 character(len=4) :: nomvar
2013 character(len=2) :: typvar
2014 character(len=12) :: etiket
2015
2016 allocate(workecr(nTrunc+1, 1))
2017
2018 !- Loop over Control Variables
2019 do var = 1, bhi%nControlVariable
2020
2021 !- Loop over vertical Levels
2022 do k = 1, bhi%controlVariable(var)%nlev
2023
2024 npak = -32
2025 dateo = 0
2026 deet = 0
2027 npas = 0
2028 ni = nTrunc + 1
2029 nj = 1
2030 nk = 1
2031 ip1 = bhi%controlVariable(var)%ip1(k)
2032 ip2 = 0
2033 ip3 = 0
2034 typvar = 'E'
2035 nomvar = trim(bhi%controlVariable(var)%nomvar(cv_type))
2036 etiket = trim(Etiket_in)
2037 grtyp = 'X'
2038 ig1 = 0
2039 ig2 = 0
2040 ig3 = 0
2041 ig4 = 0
2042 datyp = 1
2043
2044 !- Extract
2045 kgdim = bhi%controlVariable(var)%varLevIndexStart + k - 1
2046 workecr(:,1) = real(PowerSpectrum(kgdim,:),4)
2047
2048 !- Writing
2049 ier = fstecr(workecr, work, npak, iun, dateo, deet, npas, ni, nj, &
2050 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
2051 ig1, ig2, ig3, ig4, datyp, .true.)
2052
2053 end do ! Vertical Levels
2054
2055 end do ! Variables
2056
2057 deallocate(workecr)
2058
2059 end subroutine writePowerSpectrum
2060
2061 !--------------------------------------------------------------------------
2062 ! writeHorizScale
2063 !--------------------------------------------------------------------------
2064 subroutine writeHorizScale(HorizScale,iun,etiket_in,cv_type)
2065 !
2066 !:Purpose: To write the horizontal lenght scales
2067 !
2068 implicit none
2069
2070 ! Arguments:
2071 real(8), intent(in) :: HorizScale(bhi%nVarLev)
2072 integer, intent(in) :: iun
2073 integer, intent(in) :: cv_type
2074 character(len=*), intent(in) :: Etiket_in
2075
2076 ! Locals:
2077 real(4), allocatable :: workecr(:,:,:)
2078 real(4) :: work
2079 integer :: ier, fstecr, var
2080 integer :: dateo, npak, ni, nj, nk
2081 integer :: ip1, ip2, ip3, deet, npas, datyp
2082 integer :: ig1 ,ig2 ,ig3 ,ig4
2083 character(len=1) :: grtyp
2084 character(len=4) :: nomvar
2085 character(len=2) :: typvar
2086 character(len=12) :: etiket
2087
2088 !- Loop over Control Variables
2089 do var = 1, bhi%nControlVariable
2090
2091 allocate(workecr(1,1,bhi%controlVariable(var)%nlev))
2092
2093 npak = -32
2094 dateo = 0
2095 deet = 0
2096 npas = 0
2097 ni = 1
2098 nj = 1
2099 nk = bhi%controlVariable(var)%nlev
2100 ip1 = 0
2101 ip2 = 0
2102 ip3 = 0
2103 typvar = 'E'
2104 nomvar = trim(bhi%controlVariable(var)%nomvar(cv_type))
2105 etiket = trim(Etiket_in)
2106 grtyp = 'X'
2107 ig1 = 0
2108 ig2 = 0
2109 ig3 = 0
2110 ig4 = 0
2111 datyp = 1
2112
2113 !- Extract
2114 workecr(1,1,:) = real(HorizScale(bhi%controlVariable(var)%varLevIndexStart:bhi%controlVariable(var)%varLevIndexEnd),4)
2115
2116 !- Writing
2117 ier = fstecr(workecr, work, npak, iun, dateo, deet, npas, ni, nj, &
2118 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
2119 ig1, ig2, ig3, ig4, datyp, .true.)
2120
2121 deallocate(workecr)
2122
2123 end do ! Variables
2124
2125 end subroutine writeHorizScale
2126
2127 !--------------------------------------------------------------------------
2128 ! writeControlVarInfo
2129 !--------------------------------------------------------------------------
2130 subroutine writeControlVarInfo(iun)
2131 !
2132 !:Purpose: To write the control variable related info
2133 !
2134 implicit none
2135
2136 ! Arguments:
2137 integer, intent(in) :: iun
2138
2139 ! Locals:
2140 integer :: ier, fstecr, fstecr_s
2141 real(8) :: work
2142 integer :: npak, var, dateo, ni, nj
2143 integer :: ip1,ip2,ip3,deet,npas,datyp,ig1,ig2,ig3,ig4
2144 character(len=1) :: grtyp
2145 character(len=2) :: typvar
2146 character(len=4) :: nomvar
2147 character(len=4) :: ControlModelVarnameList(bhi%nControlVariable)
2148 character(len=4) :: ControlBhiVarnameList (bhi%nControlVariable)
2149 character(len=2) :: ControlVarGridTypeList (bhi%nControlVariable)
2150 integer :: ControlVarNlevList (bhi%nControlVariable)
2151
2152 !
2153 !- 1. Gathering the info
2154 !
2155
2156 !- Loop over Control Variables
2157 do var = 1, bhi%nControlVariable
2158 ControlModelVarnameList(var) = trim(bhi%controlVariable(var)%nomvar(cv_model))
2159 ControlBhiVarnameList(var) = trim(bhi%controlVariable(var)%nomvar(cv_bhi))
2160 ControlVarNlevList(var) = bhi%controlVariable(var)%nlev
2161 ControlVarGridTypeList(var) = bhi%controlVariable(var)%GridType
2162 end do
2163
2164 write(*,*)
2165 write(*,*) 'ControlModelVarnameList = ',ControlModelVarnameList(:)
2166 write(*,*) 'ControlBhiVarnameList = ',ControlBhiVarnameList(:)
2167 write(*,*) 'ControlVarNlevList = ',ControlVarNlevList(:)
2168 write(*,*) 'ControlVarGridTypeList = ',ControlVarGridTypeList(:)
2169
2170 !
2171 !- 2. Writing the list of control variables and number of vertical levels
2172 !
2173 npak = -32
2174 dateo = 0
2175 deet = 0
2176 ip1 = 0
2177 ip2 = 0
2178 ip3 = 0
2179 npas = 0
2180 grtyp = 'X'
2181 typvar = 'X'
2182 ig1 = 0
2183 ig2 = 0
2184 ig3 = 0
2185 ig4 = 0
2186
2187 nomvar = 'CVN'
2188 ni = 4 ! 4 Characters
2189 nj = bhi%nControlVariable
2190 datyp = 7 ! Character
2191
2192 ier = fstecr_s(ControlModelVarnameList, work, npak, &
2193 iun, dateo, deet, npas, ni, nj, 1, ip1, &
2194 ip2, ip3, typvar, nomvar, 'MODEL', grtyp, ig1, &
2195 ig2, ig3, ig4, datyp, .true.)
2196
2197 ier = fstecr_s(ControlBhiVarnameList, work, npak, &
2198 iun, dateo, deet, npas, ni, nj, 1, ip1, &
2199 ip2, ip3, typvar, nomvar, 'B_HI', grtyp, ig1, &
2200 ig2, ig3, ig4, datyp, .true.)
2201
2202 nomvar = 'CVL'
2203 ni = 2 ! 2 Characters
2204 ier = fstecr_s(ControlVarGridTypeList, work, npak, &
2205 iun, dateo, deet, npas, ni, nj, 1, ip1, &
2206 ip2, ip3, typvar, nomvar, 'LEVTYPE', grtyp, ig1, &
2207 ig2, ig3, ig4, datyp, .true.)
2208
2209 datyp = 2 ! Integer
2210 ni = bhi%nControlVariable
2211 nj = 1
2212 ier = fstecr(ControlVarNlevList, work, npak, &
2213 iun, dateo, deet, npas, ni, nj, 1, ip1, &
2214 ip2, ip3, typvar, nomvar, 'NLEV', grtyp, ig1, &
2215 ig2, ig3, ig4, datyp, .true.)
2216
2217 end subroutine writeControlVarInfo
2218
2219 !--------------------------------------------------------------------------
2220 ! writePressureProfiles
2221 !--------------------------------------------------------------------------
2222 subroutine writePressureProfiles
2223 !
2224 !:Purpose: To write the MM and TH pressure profiles used for vertical localization
2225 !
2226 implicit none
2227
2228 ! Locals:
2229 character(len=128) :: outfilename
2230 integer :: jk
2231
2232 outfilename = "./pressureProfile_M.txt"
2233 open (unit=99,file=outfilename,action="write",status="new")
2234 do jk = 1, vco_bhi%nlev_M
2235 write(99,'(I3,2X,F6.1)') jk, pressureProfile_M(jk)/100.d0
2236 end do
2237 close(unit=99)
2238
2239 outfilename = "./pressureProfile_T.txt"
2240 open (unit=99,file=outfilename,action="write",status="new")
2241 do jk = 1, vco_bhi%nlev_T
2242 write(99,'(I3,2X,F6.1)') jk, pressureProfile_T(jk)/100.d0
2243 end do
2244 close(unit=99)
2245
2246 write(6,*) 'finished writing pressure profiles...'
2247 call flush(6)
2248
2249 end subroutine writePressureProfiles
2250
2251 !--------------------------------------------------------------------------
2252 ! calcLocalCorrelations
2253 !--------------------------------------------------------------------------
2254 subroutine calcLocalCorrelations(ensPerts)
2255 !
2256 !:Purpose: To compute the local horizontal correlation for some 'reference' grid points
2257 !
2258 implicit none
2259
2260 ! Arguments:
2261 type(struct_ens), intent(inout) :: ensPerts
2262
2263 ! Locals:
2264 type(struct_gsv) :: statevector_locHorizCor
2265 type(struct_gsv) :: statevector_oneMember
2266 type(struct_gsv) :: statevector_oneMemberTiles
2267 real(8), pointer :: ptr3d_r8(:,:,:)
2268 real(8), pointer :: ptr3d_r8_oneMember(:,:,:)
2269 real(8) :: dnEns
2270 integer :: i, j, k, ens
2271 integer :: blocklength_x, blocklength_y, blockpadding, nirefpoint, njrefpoint
2272 integer :: iref_id, jref_id, iref, jref
2273 integer :: imin, imax, jmin, jmax
2274 character(len=4), pointer :: varNamesList(:)
2275 integer :: ier, fclos, fnom, nulnam
2276
2277 NAMELIST /NAMHVCORREL_LOCAL/nirefpoint, njrefpoint, blockpadding
2278
2279 !
2280 ! To compute the local horizontal correlation for some 'reference' grid point
2281 ! ... we assume that the ensemble grid point mean was removed and that
2282 ! the ensemble values were divided by the grid point std dev.
2283 !
2284
2285 nirefpoint = 4 ! Number of reference grid point in x
2286 njrefpoint = 2 ! Number of reference grid point in y
2287 blockpadding = 4 ! Number of grid point padding between blocks (to set correlation to 0 between each block)
2288
2289 nulnam = 0
2290 ier = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
2291 read(nulnam,nml=NAMHVCORREL_LOCAL)
2292 write(*,nml=NAMHVCORREL_LOCAL)
2293 ier = fclos(nulnam)
2294
2295 blocklength_x = hco_ens%ni / nirefpoint ! Horizontal correlation will be compute blocklength x blocklength gridpoint
2296 ! around each reference point
2297 blocklength_y = hco_ens%nj / njrefpoint ! Horizontal correlation will be compute blocklength x blocklength gridpoint
2298 ! around each reference point
2299
2300 nullify(varNamesList)
2301 call ens_varNamesList(varNamesList,ensPerts)
2302
2303 call gsv_allocate(statevector_locHorizCor, ens_getNumStep(ensPerts), &
2304 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
2305 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
2306 mpi_distribution_opt='VarsLevs', dataKind_opt=8 )
2307
2308 call gsv_allocate(statevector_oneMemberTiles, ens_getNumStep(ensPerts), &
2309 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
2310 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
2311 mpi_distribution_opt='Tiles', dataKind_opt=8 )
2312
2313 call gsv_allocate(statevector_oneMember, ens_getNumStep(ensPerts), &
2314 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
2315 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
2316 mpi_distribution_opt='VarsLevs', dataKind_opt=8 )
2317
2318 call gsv_zero(statevector_locHorizCor)
2319
2320 dnEns = 1.0d0/dble(nEns-1)
2321
2322 call gsv_getField(statevector_locHorizCor,ptr3d_r8)
2323
2324 do ens = 1, nEns
2325 call ens_copyMember(ensPerts, statevector_oneMemberTiles, ens)
2326 call gsv_transposeTilesToVarsLevs(statevector_oneMemberTiles, statevector_oneMember)
2327 call gsv_getField(statevector_oneMember,ptr3d_r8_oneMember)
2328
2329 do k = statevector_locHorizCor%mykBeg, statevector_locHorizCor%mykEnd
2330 do jref_id = 1, njrefpoint
2331 do iref_id = 1, nirefpoint
2332 iref = (2*iref_id-1)*blocklength_x/2
2333 jref = (2*jref_id-1)*blocklength_y/2
2334 jmin = max(jref-(blocklength_y-blockpadding)/2,1)
2335 jmax = min(jref+(blocklength_y-blockpadding)/2,hco_ens%nj)
2336 imin = max(iref-(blocklength_x-blockpadding)/2,1)
2337 imax = min(iref+(blocklength_x-blockpadding)/2,hco_ens%ni)
2338 do j = jmin, jmax
2339 do i = imin, imax
2340 ptr3d_r8(i,j,k) = ptr3d_r8(i,j,k) + &
2341 ptr3d_r8_oneMember(i,j,k)*ptr3d_r8_oneMember(iref,jref,k)
2342 end do
2343 end do
2344 end do
2345 end do
2346 end do
2347
2348 end do
2349
2350 call gsv_scale(statevector_locHorizCor,dnEns)
2351
2352 write(*,*) 'finished computing the local horizontal correlations...'
2353
2354 !
2355 !- 4. Write to file
2356 !
2357 call gio_writeToFile(statevector_locHorizCor, './horizCorrelLocal.fst', 'HCORREL_LOC', &
2358 typvar_opt = 'E', numBits_opt = 32)
2359
2360 call gsv_deallocate(statevector_locHorizCor)
2361 call gsv_deallocate(statevector_oneMember)
2362 call gsv_deallocate(statevector_oneMemberTiles)
2363
2364 end subroutine calcLocalCorrelations
2365
2366end module calcStatsLam_mod
