1module menetrierDiag_mod
2 ! MODULE menetrierDiag_mod (prefix='bmd' category='1. High-level functionality')
3 !
4 !:Purpose: To compute optimal localization radii according to the theory
5 ! developed by Benjamin Menetrier (Meteo-France) and reported
6 ! in Menetrier, Michel, Montmerle and Berre, 2015, Parts 1 and 2.
7 !
8 use earthConstants_mod
9 use utilities_mod
10 use localizationFunction_mod
11 use varNameList_mod
12 use horizontalCoord_mod
13 use gridStatevector_mod
14 use ensembleStatevector_mod
15 use timeCoord_mod
16 use midasMpi_mod
17 implicit none
18 save
19 private
20
21 ! Public Subroutines
22 public :: bmd_setup, bmd_localizationRadii
23
24 real(8), pointer :: pressureProfile_M(:), pressureProfile_T(:)
25
26 type(struct_hco), pointer :: hco_ens ! Ensemble horizontal grid parameters
27
28 integer :: nens, ni, nj, nLevEns_M, nLevEns_T, nkgdimEns
29 integer :: nvar3d, nvar2d, nWaveBand
30
31 integer, allocatable :: varLevOffset(:)
32
33 character(len=4), allocatable :: nomvar3d(:), nomvar2d(:)
34
35 logical :: global
36
37 logical :: initialized = .false.
38
39 ! Namelist variables:
40 integer :: strideForHLoc
41 integer :: strideForVloc
42 integer :: horizPadding
43 logical :: hLoc
44 logical :: vLoc
45
46contains
47
48 !--------------------------------------------------------------------------
49 ! bmd_setup
50 !--------------------------------------------------------------------------
51 subroutine bmd_setup(statevector_template, hco_core_in, nens_in, pressureProfile_M_in, &
52 pressureProfile_T_in, nWaveBand_in)
53 implicit none
54
55 ! Arguments:
56 type(struct_gsv), intent(in) :: statevector_template
57 type(struct_hco), pointer, intent(in) :: hco_core_in
58 integer, intent(in) :: nens_in
59 integer, intent(in) :: nWaveBand_in
60 real(8), pointer, intent(inout) :: pressureProfile_M_in(:)
61 real(8), pointer, intent(inout) :: pressureProfile_T_in(:)
62
63 ! Locals:
64 integer :: nulnam, ierr, fclos, fnom
65 integer :: nVar, varNameIndex, var2dIndex, var3dIndex
66 character(len=4), pointer :: varNamesList(:)
67
68 NAMELIST /NAMLOCALIZATIONRADII/strideForHLoc,strideForVloc,hLoc,vLoc,horizPadding
69
70 !
71 !- 1. Input parameters
72 !
73 hco_ens => hco_core_in
74 nens = nens_in
75 ni = hco_ens%ni
76 nj = hco_ens%nj
77 nLevEns_M = gsv_getNumLev(statevector_template,'MM') !nLevEns_M_in
78 nLevEns_T = gsv_getNumLev(statevector_template,'TH') !nLevEns_T_in
79 nkgdimEns = statevector_template%nk
80 pressureProfile_M => pressureProfile_M_in
81 pressureProfile_T => pressureProfile_T_in
82 nWaveBand = nWaveBand_in
83 global = hco_ens%global
84
85 nullify(varNamesList)
86 call gsv_varNamesList(varNamesList, statevector_template)
87
88 nVar = size(varNamesList)
89 allocate(varLevOffset(nVar))
90 nVar3d = 0
91 nVar2d = 0
92 do varNameIndex = 1, size(varNamesList)
93 varLevOffset(varNameIndex) = gsv_getOffsetFromVarName(statevector_template,varNamesList(varNameIndex))
94 if (vnl_varLevelFromVarname(varNamesList(varNameIndex)) == 'SF' ) then
95 nVar2d = nVar2d + 1
96 else
97 nVar3d = nVar3d + 1
98 end if
99 end do
100
101 allocate(nomvar3d(nvar3d))
102 allocate(nomvar2d(nvar2d))
103
104 var2dIndex = 0
105 var3dIndex = 0
106 do varNameIndex = 1, size(varNamesList)
107 if (vnl_varLevelFromVarname(varNamesList(varNameIndex)) == 'SF' ) then
108 var2dIndex = var2dIndex + 1
109 nomvar2d(var2dIndex) = varNamesList(varNameIndex)
110 else
111 var3dIndex = var3dIndex + 1
112 nomvar3d(var3dIndex) = varNamesList(varNameIndex)
113 end if
114 end do
115
116 !
117 !- 2. Read Namelist
118 !
119 hLoc = .true.
120 vLoc = .true.
121 strideForHLoc = 100 ! Horizontal correlations will be computed every "stride" gridpoint in x and y
122 strideForVLoc = 50 ! Vertical correlations will be computed every "stride" gridpoint in x and y
123 horizPadding = 0 ! Number of grid point to discard along the horizontal edges (only for LAM)
124
125 nulnam = 0
126 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
127 read(nulnam,nml=NAMLOCALIZATIONRADII)
128 write(*,nml=NAMLOCALIZATIONRADII)
129 ierr = fclos(nulnam)
130
131 !
132 !- 3. Ending
133 !
134 initialized = .true.
135
136 end subroutine bmd_setup
137
138 !--------------------------------------------------------------------------
139 ! bmd_localizationRadii
140 !--------------------------------------------------------------------------
141 subroutine bmd_localizationRadii(ensPerts,waveBandIndex_opt)
142 implicit none
143
144 ! Arguments:
145 type(struct_ens), intent(inout) :: ensPerts
146 integer,optional, intent(in) :: waveBandIndex_opt
147
148 !
149 !- Estimate the horionzontal and vertical localization radii
150 !
151 if ( hLoc ) then
152 call calcHorizLocalizationRadii(ensPerts, strideForHLoc, & ! IN
153 waveBandIndex_opt=waveBandIndex_opt) ! IN
154 end if
155 if ( vLoc ) then
156 call calcVertLocalizationRadii(ensPerts, strideForVloc, & ! IN
157 waveBandIndex_opt=waveBandIndex_opt) ! IN
158 end if
159
160 end subroutine bmd_localizationRadii
161
162 !--------------------------------------------------------------------------
163 ! CALCHORIZLOCALIZATIONRADII
164 !--------------------------------------------------------------------------
165 subroutine calcHorizLocalizationRadii(ensPerts,stride,waveBandIndex_opt)
166 implicit none
167
168 ! Arguments:
169 type(struct_ens), intent(inout) :: ensPerts
170 integer, intent(in) :: stride
171 integer, optional, intent(in) :: waveBandIndex_opt
172
173 ! Locals:
174 type(struct_gsv) :: statevector_ensStdDev
175 type(struct_gsv) :: statevector_ensStdDev_tiles
176 type(struct_gsv) :: statevector_oneMemberTiles
177 type(struct_gsv) :: statevector_oneMember(nens)
178 real(8), pointer :: ensStdDev(:,:,:)
179 real(4), pointer :: ptr3d_r4(:,:,:)
180 real(4), allocatable :: ensPert_local(:,:,:)
181 real(8), allocatable :: meanCorrel(:,:)
182 real(8), allocatable :: meanCorrelSquare(:,:)
183 real(8), allocatable :: meanVarianceProduct(:,:)
184 real(8), allocatable :: meanCovarianceSquare(:,:)
185 real(8), allocatable :: meanFourthMoment(:,:)
186 real(8), allocatable :: meanCorrel_local(:,:)
187 real(8), allocatable :: meanCorrelSquare_local(:,:)
188 real(8), allocatable :: meanVarianceProduct_local(:,:)
189 real(8), allocatable :: meanCovarianceSquare_local(:,:)
190 real(8), allocatable :: meanFourthMoment_local(:,:)
191 real(8), allocatable :: localizationFunctions(:,:,:) ! Eq. 19-21 in MMMB 2015 Part 2
192 real(8), allocatable :: localizationRadii(:,:)
193 real(8), allocatable :: distanceBinThresholds(:) ! Maximum distance for each distance-bin
194 real(8), allocatable :: distanceBinMean(:) ! Mean distance for each distance-bin
195 real(8), allocatable :: distanceBinWeight(:) ! Weight given to each bin in the curve fitting step
196 real(8), allocatable :: gridPointWeight(:,:,:) ! Weight given to grid point in the statistic computation
197 real(8), allocatable :: sumWeight(:,:) ! Sample size for each distance-bin
198 real(8), allocatable :: sumWeight_local(:,:)
199 real(8), pointer :: PressureProfile(:)
200 logical, allocatable :: gridPointAlreadyUsed(:,:)
201 real(8) :: dnens, correlation, covariance, fourthMoment, distance, maxDistance, weight
202 real(8) :: t1, t2, t3, rmse
203 integer :: i, j, k, f, ens, bin, numbins, numFunctions, nSize
204 integer :: iref, jref, ier
205 integer :: nLevEns, jvar, mykBeg, mykEnd
206 character(len=128) :: outfilename
207 character(len=2) :: wbnum
208 character(len=4), pointer :: varNamesList(:)
209
210 !
211 !- 1. Setup
212 !
213 call ens_copyEnsStdDev(ensPerts, statevector_ensStdDev_tiles)
214
215 nullify(varNamesList)
216 call ens_varNamesList(varNamesList,ensPerts)
217 call gsv_allocate(statevector_ensStdDev, ens_getNumStep(ensPerts), &
218 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
219 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
220 mpi_distribution_opt='VarsLevs', dataKind_opt=8 )
221 call gsv_transposeTilesToVarsLevs(statevector_ensStdDev_tiles, statevector_ensStdDev)
222 call gsv_getField(statevector_ensStdDev,ensStdDev)
223
224 call gsv_deallocate(statevector_ensStdDev_tiles)
225
226 mykBeg = statevector_ensStdDev%mykBeg
227 mykEnd = statevector_ensStdDev%mykEnd
228
229 numFunctions = 3
230
231 if ( global ) then
232 numBins = ni/4 ! 1/4 of the Earth
233 if ( horizPadding /= 0 ) then
234 write(*,*)
235 write(*,*) 'WARNING: horizPadding /= 0. Force it to zero since we are in global mode.'
236 horizPadding = 0
237 end if
238 else
239 numBins = ni-(2*horizPadding)
240 if ( horizPadding == 0 ) then
241 write(*,*)
242 write(*,*) 'WARNING: horizPadding == 0. The rim and blending zone WILL HAVE an impact on the results.'
243 end if
244 end if
245
246 allocate(meanCorrelSquare(numBins,nkgdimEns))
247 allocate(meanCorrel(numBins,nkgdimEns))
248 allocate(meanVarianceProduct(numBins,nkgdimEns))
249 allocate(meanCovarianceSquare(numBins,nkgdimEns))
250 allocate(meanFourthMoment(numBins,nkgdimEns))
251 allocate(sumWeight(numBins,nkgdimEns))
252
253 allocate(meanCorrelSquare_local(numBins,nkgdimEns))
254 allocate(meanCorrel_local(numBins,nkgdimEns))
255 allocate(meanVarianceProduct_local(numBins,nkgdimEns))
256 allocate(meanCovarianceSquare_local(numBins,nkgdimEns))
257 allocate(meanFourthMoment_local(numBins,nkgdimEns))
258 allocate(sumWeight_local(numBins,nkgdimEns))
259
260 allocate(gridPointWeight(ni,nj,nkgdimEns))
261
262 allocate(distanceBinThresholds(numBins))
263
264 ! Assign the (upper) threshold to each separation-distance-bin
265 write(*,*)
266 write(*,*) 'Separation distance bin'
267 do bin = 1, numbins
268 distanceBinThresholds(bin) = calcDistance(hco_ens%lat(nj/2),hco_ens%lon(1+horizPadding),hco_ens%lat(nj/2),hco_ens%lon(bin+horizPadding))
269 write(*,*) bin, hco_ens%lat(nj/2),hco_ens%lon(1+horizPadding),hco_ens%lon(bin+horizPadding), distanceBinThresholds(bin)/1000.d0
270 end do
271
272 maxDistance=distanceBinThresholds(numBins)
273
274 dnens = 1.0d0/dble(nens-1)
275
276 ! Grid point Weight
277 write(*,*)
278 write(*,*) 'Grid point Weight'
279 do j = 1, nj
280 gridPointWeight(:,j,:) = cos(hco_ens%lat(j))
281 write(*,*) j, hco_ens%lat(j), cos(hco_ens%lat(j))
282 end do
283
284 !
285 !- 2. Estimation of localization functions
286 !
287
288 !- 2.1 Computation of various statistics for different separation distances
289 meanCorrelSquare_local(:,:) = 0.d0
290 meanCorrel_local(:,:) = 0.d0
291 meanCovarianceSquare_local(:,:) = 0.d0
292 meanVarianceProduct_local(:,:) = 0.d0
293 meanFourthMoment_local(:,:) = 0.d0
294 sumWeight_local(:,:) = 0.d0
295
296 allocate(ensPert_local(nens,ni,nj))
297 allocate(gridPointAlreadyUsed(ni,nj))
298
299 call gsv_allocate(statevector_oneMemberTiles, ens_getNumStep(ensPerts), &
300 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
301 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
302 mpi_distribution_opt='Tiles', dataKind_opt=4 )
303
304 do ens = 1, nens
305 call gsv_allocate(statevector_oneMember(ens), ens_getNumStep(ensPerts), &
306 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
307 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
308 mpi_distribution_opt='VarsLevs', dataKind_opt=4 )
309 call ens_copyMember(ensPerts, statevector_oneMemberTiles, ens)
310 call gsv_transposeTilesToVarsLevs(statevector_oneMemberTiles, statevector_oneMember(ens))
311 end do
312
313 call gsv_deallocate(statevector_oneMemberTiles)
314
315 !$OMP PARALLEL DO PRIVATE (ptr3d_r4,k,iref,jref,j,i,ens,correlation,covariance,fourthMoment,distance,bin,weight,ensPert_local,gridPointAlreadyUsed)
316 do k = mykBeg, mykEnd
317 write(*,*) 'Computing distance-bin statistics for ensemble level: ', k
318
319 !- Select data needed to speed up the process (ensemble member index must come first in ensPert_Local
320 ! because ensemble member is the last loop index below)
321 do ens = 1, nens
322 call gsv_getField(statevector_oneMember(ens),ptr3d_r4)
323 do j = 1, nj
324 do i = 1, ni
325 ensPert_local(ens,i,j) = ptr3d_r4(i,j,k)
326 end do
327 end do
328 end do
329
330 gridPointAlreadyUsed(:,:) = .false.
331
332 do jref = nint(stride/2.0)+horizPadding, nj-horizPadding, stride ! Pick every stride point to save cost.
333 do iref = nint(stride/2.0)+horizPadding, ni-horizPadding, stride ! Pick every stride point to save cost.
334
335 if (k == 1) then
336 write(*,*) 'grid point info', iref, jref
337 end if
338
339 do j = 1+horizPadding, nj-horizPadding
340 do i = 1+horizPadding, ni-horizPadding
341
342 if ( gridPointAlreadyUsed(i,j) ) cycle ! prevent using the same pair of points more than once
343
344 distance=calcDistance(hco_ens%lat(jref),hco_ens%lon(iref),hco_ens%lat(j),hco_ens%lon(i))
345 if (distance <= maxDistance .and. gridPointWeight(i,j,k) > 0.d0) then
346 covariance = 0.d0
347 fourthMoment = 0.d0
348 do ens = 1, nens
349 covariance = covariance + &
350 ensPert_local(ens,i,j) * ensPert_local(ens,iref,jref)
351 fourthMoment = fourthMoment + &
352 (ensPert_local(ens,i,j) * ensPert_local(ens,iref,jref))**2
353 end do
354 covariance = covariance * dnens
355 fourthMoment = fourthMoment / real(nens,8)
356 if ( ensStdDev(i,j,k) > 0.d0 .and. ensStdDev(iref,jref,k) > 0.d0 ) then
357 correlation = covariance / (ensStdDev(i,j,k)*ensStdDev(iref,jref,k))
358 else
359 correlation = 0.d0
360 end if
361
362 bin=findBinIndex(distance,distanceBinThresholds,numBins)
363
364 weight = sqrt(gridPointWeight(iref,jref,k)) * sqrt(gridPointWeight(i,j,k))
365
366 sumWeight_local(bin,k) = sumWeight_local(bin,k) + weight
367
368 meanCorrel_local(bin,k) = meanCorrel_local(bin,k) + &
369 correlation * weight
370 meanCorrelSquare_local(bin,k) = meanCorrelSquare_local(bin,k) + &
371 correlation**2 * weight
372 meanCovarianceSquare_local(bin,k) = meanCovarianceSquare_local(bin,k) + &
373 covariance**2 * weight
374 meanVarianceProduct_local(bin,k) = meanVarianceProduct_local(bin,k) + &
375 ensStdDev(i,j,k)**2 * ensStdDev(iref,jref,k)**2 * weight
376 meanFourthMoment_local(bin,k) = meanFourthMoment_local(bin,k) + &
377 fourthMoment * weight
378 end if
379 end do
380 end do
381
382 ! From now on, omit the current reference point
383 gridPointAlreadyUsed(iref,jref) = .true.
384
385 end do ! iref
386 end do ! jref
387
388 end do ! nkgdimEns
389 !$OMP END PARALLEL DO
390
391 deallocate(ensPert_local)
392 deallocate(gridPointAlreadyUsed)
393 do ens = 1, nens
394 call gsv_deallocate(statevector_oneMember(ens))
395 end do
396
397 !- 2.2 Gather the all the info in processor 0
398 nSize = nkgdimEns * numbins
399 call rpn_comm_reduce(sumWeight_local ,sumWeight ,nSize, &
400 "mpi_double_precision","mpi_sum",0,"GRID",ier)
401 call rpn_comm_reduce(meanCorrel_local ,meanCorrel ,nSize, &
402 "mpi_double_precision","mpi_sum",0,"GRID",ier)
403 call rpn_comm_reduce(meanCorrelSquare_local ,meanCorrelSquare ,nSize, &
404 "mpi_double_precision","mpi_sum",0,"GRID",ier)
405 call rpn_comm_reduce(meanVarianceProduct_local ,meanVarianceProduct ,nSize, &
406 "mpi_double_precision","mpi_sum",0,"GRID",ier)
407 call rpn_comm_reduce(meanFourthMoment_local ,meanFourthMoment ,nSize, &
408 "mpi_double_precision","mpi_sum",0,"GRID",ier)
409 call rpn_comm_reduce(meanCovarianceSquare_local,meanCovarianceSquare,nSize, &
410 "mpi_double_precision","mpi_sum",0,"GRID",ier)
411
412 deallocate(sumWeight_local)
413 deallocate(meanCorrel_local)
414 deallocate(meanCorrelSquare_local)
415 deallocate(meanVarianceProduct_local)
416 deallocate(meanFourthMoment_local)
417 deallocate(meanCovarianceSquare_local)
418
419 if (mmpi_myid == 0) then
420 !- 2.3 Computation of the localization functions
421 allocate(localizationFunctions(numFunctions,numBins,nkgdimEns))
422
423 t1=dble((nens-1)**2)/dble(nens*(nens-3))
424 t2=dble(nens)/dble((nens-2)*(nens-3))
425 t3=dble(nens-1)/dble(nens*(nens-2)*(nens-3))
426 !$OMP PARALLEL DO PRIVATE (k,bin)
427 do k = 1, nkgdimEns
428 do bin = 1, numbins
429
430 meanCorrel(bin,k) = meanCorrel(bin,k) / sumWeight(bin,k)
431 meanCorrelSquare(bin,k) = meanCorrelSquare(bin,k) / sumWeight(bin,k)
432 meanCovarianceSquare(bin,k) = meanCovarianceSquare(bin,k) / sumWeight(bin,k)
433 meanVarianceProduct(bin,k) = meanVarianceProduct(bin,k) / sumWeight(bin,k)
434 meanFourthMoment(bin,k) = meanFourthMoment(bin,k) / sumWeight(bin,k)
435
436 if ( meanCovarianceSquare(bin,k) /= 0.d0 ) then
437 ! Form 1: General formulation (Eq. 19 in MMMB 2015 Part 2)
438 localizationFunctions(1,bin,k) = t1 - t2*meanFourthMoment(bin,k)/meanCovarianceSquare(bin,k) + &
439 t3*meanVarianceProduct(bin,k)/meanCovarianceSquare(bin,k)
440 ! Form 2: Gaussian sample distribution (Eq. 20 in MMMB 2015 Part 2)
441 localizationFunctions(2,bin,k) = dble(nens-1)/dble((nens+1)*(nens-2)) * &
442 (dble(nens-1)-meanVarianceProduct(bin,k)/meanCovarianceSquare(bin,k))
443 else
444 write(*,*) ' !!! Warning !!! meanCovarianceSquare = 0 in bin, level = ', bin, k
445 localizationFunctions(1,bin,k) = 0.d0
446 localizationFunctions(2,bin,k) = 0.d0
447 end if
448 ! Form 3: Gaussian sample distribution and correlation-based formulation (Eq. 21 in MMMB 2015 Part 2)
449 if ( meanCorrelSquare(bin,k) /= 0.d0 ) then
450 localizationFunctions(3,bin,k) = dble(nens-1)/dble((nens+1)*(nens-2)) * &
451 (dble(nens-1)-1.d0/meanCorrelSquare(bin,k))
452 else
453 write(*,*) ' !!! Warning !!! meanCorrelSquare = 0 in bin, level = ', bin, k
454 localizationFunctions(3,bin,k) = 0.d0
455 end if
456 end do
457 end do
458 !$OMP END PARALLEL DO
459
460 !
461 !- 3. Estimation of localization radii (curve fitting)
462 !
463 allocate(localizationRadii(numFunctions,nkgdimEns))
464 allocate(distanceBinMean(numBins))
465 allocate(distanceBinWeight(numBins))
466
467 call lfn_setup('FifthOrder') ! IN
468
469 localizationRadii(:,:) = 2000.d0*1000.d0 ! First Guess (meter)
470 distanceBinWeight(:) = 1.d0 ! Even weight
471 do bin = 1, numBins
472 if (bin == 1) then
473 distanceBinMean(bin) = distanceBinThresholds(bin) ! = 0 m
474 else
475 distanceBinMean(bin) = 0.5d0*(distanceBinThresholds(bin)+distanceBinThresholds(bin-1))
476 end if
477 end do
478
479 do f = 1, numFunctions
480 write(*,*)
481 write(*,*) 'Localization Function : ', f
482 do k = 1, nkgdimEns
483 write(*,*) ' ----- EnsLev : ', k
484 call lfn_lengthscale( localizationRadii(f,k), & ! INOUT
485 rmse, & ! OUT
486 localizationFunctions(f,:,k), & ! IN
487 distanceBinMean, & ! IN
488 distanceBinWeight, & ! IN
489 numbins ) ! IN
490 end do
491 end do
492
493 !
494 !- 4. Write to file
495 !
496 if ( nWaveBand /= 1 ) then
497 if (.not. present(waveBandIndex_opt)) then
498 write(*,*) 'calcLocalizationRadii: No waveBandIndex was supplied!!!'
499 call utl_abort('calbmatrix_glb')
500 end if
501 write(wbnum,'(I2.2)') waveBandIndex_opt
502 end if
503
504 !- 4.1 Localization functions in txt format (for plotting purposes)
505 do jvar = 1, nvar3d
506 if ( nWaveBand == 1 ) then
507 outfilename = "./horizLocalizationFunctions_"//trim(nomvar3d(jvar))//".txt"
508 else
509 outfilename = "./horizLocalizationFunctions_"//trim(nomvar3d(jvar))//"_"//wbnum//".txt"
510 end if
511 open (unit=99,file=outfilename,action="write",status="new")
512
513 if(vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
514 nLevEns = nLevEns_M
515 else
516 nLevEns = nLevEns_T
517 endif
518 do k=1,nlevEns
519 do bin = 1, numbins
520 write(99,'(I3,2X,I3,2X,F7.1,2X,I7,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,E10.3,2X,E10.3,2X,E10.3)') k, bin, &
521 distanceBinThresholds(bin)/1000.d0, nint(sumWeight(bin,varLevOffset(jvar)+k)), &
522 meanCorrel(bin,varLevOffset(jvar)+k), &
523 localizationFunctions(3,bin,varLevOffset(jvar)+k), &
524 localizationFunctions(2,bin,varLevOffset(jvar)+k), &
525 localizationFunctions(1,bin,varLevOffset(jvar)+k), &
526 meanCorrelSquare(bin,varLevOffset(jvar)+k), &
527 meanCovarianceSquare(bin,varLevOffset(jvar)+k), &
528 meanVarianceProduct(bin,varLevOffset(jvar)+k), &
529 meanFourthMoment(bin,varLevOffset(jvar)+k)
530 end do
531 end do
532 close(unit=99)
533 end do
534
535 do jvar = 1, nvar2d
536 k = varLevOffset(nvar3d+1)+jvar
537 if ( nWaveBand == 1 ) then
538 outfilename = "./horizLocalizationFunctions_"//trim(nomvar2d(jvar))//".txt"
539 else
540 outfilename = "./horizLocalizationFunctions_"//trim(nomvar2d(jvar))//"_"//wbnum//".txt"
541 end if
542 open (unit=99,file=outfilename,action="write",status="new")
543 do bin = 1, numbins
544 write(99,'(I3,2X,I3,2X,F7.1,2X,I7,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,E10.3,2X,E10.3,2X,E10.3)') 1, bin, &
545 distanceBinThresholds(bin)/1000.d0, nint(sumWeight(bin,k)), &
546 meanCorrel(bin,k), &
547 localizationFunctions(3, bin,k), &
548 localizationFunctions(2, bin,k), &
549 localizationFunctions(1, bin,k), &
550 meanCorrelSquare(bin,k), &
551 meanCovarianceSquare(bin,k), &
552 meanVarianceProduct(bin,k), &
553 meanFourthMoment(bin,k)
554 end do
555 close(unit=99)
556 end do
557
558 !- 4.2 Localization radii in txt format (for plotting purposes)
559 do jvar = 1, nvar3d
560 if ( nWaveBand == 1 ) then
561 outfilename = "./horizLocalizationRadii_"//trim(nomvar3d(jvar))//".txt"
562 else
563 outfilename = "./horizLocalizationRadii_"//trim(nomvar3d(jvar))//"_"//wbnum//".txt"
564 end if
565 open (unit=99,file=outfilename,action="write",status="new")
566
567 if(vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
568 nLevEns = nLevEns_M
569 PressureProfile => pressureProfile_M
570 else
571 nLevEns = nLevEns_T
572 PressureProfile => pressureProfile_T
573 endif
574 do k=1,nlevEns
575 write(99,'(I3,2X,F7.2,2X,F7.1,2X,F7.1,2X,F7.1)') k, PressureProfile(k)/100.d0, &
576 min(localizationRadii(3,varLevOffset(jvar)+k)/1000.d0,99999.9d0), &
577 min(localizationRadii(2,varLevOffset(jvar)+k)/1000.d0,99999.9d0), &
578 min(localizationRadii(1,varLevOffset(jvar)+k)/1000.d0,99999.9d0)
579 end do
580 close(unit=99)
581 end do
582
583 do jvar = 1, nvar2d
584 k = varLevOffset(nvar3d+1)+jvar
585 if ( nWaveBand == 1 ) then
586 outfilename = "./horizLocalizationRadii_"//trim(nomvar2d(jvar))//".txt"
587 else
588 outfilename = "./horizLocalizationRadii_"//trim(nomvar2d(jvar))//"_"//wbnum//".txt"
589 end if
590 open (unit=99,file=outfilename,action="write",status="new")
591 write(99,'(I3,2X,F7.2,2X,F7.1,2X,F7.1,2X,F7.1)') 1, 1010.0, &
592 min(localizationRadii(3,k)/1000.d0,99999.9d0), &
593 min(localizationRadii(2,k)/1000.d0,99999.9d0), &
594 min(localizationRadii(1,k)/1000.d0,99999.9d0)
595 close(unit=99)
596 end do
597
598 deallocate(localizationFunctions)
599 deallocate(localizationRadii)
600 deallocate(distanceBinMean)
601 deallocate(distanceBinWeight)
602
603 end if ! mmpi_myid == 0
604
605 deallocate(meanCorrelSquare)
606 deallocate(meanCorrel)
607 deallocate(meanCovarianceSquare)
608 deallocate(meanVarianceProduct)
609 deallocate(meanFourthMoment)
610 deallocate(distanceBinThresholds)
611 deallocate(sumWeight)
612 deallocate(gridPointWeight)
613
614 write(*,*)
615 write(*,*) 'finished estimating the horizontal localization radii...'
616
617 end subroutine calcHorizLocalizationRadii
618
619 !--------------------------------------------------------------------------
620 ! FINDBININDEX
621 !--------------------------------------------------------------------------
622 function findBinIndex(distance,distanceBinThresholds,numBins) result(binIndex)
623 implicit none
624
625 ! Arguments:
626 integer, intent(in) :: numBins
627 real(8), intent(in) :: distance
628 real(8), intent(in) :: distanceBinThresholds(numBins)
629 ! Result:
630 integer :: binIndex
631
632 ! Locals:
633 integer :: bin
634
635 binIndex = -1
636 do bin = 1, numbins
637 if ( distance <= distanceBinThresholds(bin) ) then
638 binIndex = bin
639 exit
640 end if
641 end do
642
643 if (binIndex == -1) then
644 write(*,*) 'findBinIndex: No match found! ABORTING'
645 call utl_abort('findBinIndex')
646 end if
647
648 end function findBinIndex
649
650 !--------------------------------------------------------------------------
651 ! DISTANCE
652 !--------------------------------------------------------------------------
653 function calcDistance(lat2, lon2, lat1, lon1) result(distanceInM)
654 !:Purpose: To compute the distance between two points on Earth: (lat1,lon1)
655 ! and (lat2,lon2). Calcul utilisant la Formule d'Haversine
656 ! Reference: R.W. Sinnott,'Virtues of Haversine',Sky and Telescope,
657 ! vol.68, no.2, 1984, p.159)
658 implicit none
659
660 ! Arguments:
661 real(8), intent(in) :: lat1
662 real(8), intent(in) :: lon1
663 real(8), intent(in) :: lat2
664 real(8), intent(in) :: lon2
665 ! Result:
666 real(8) :: distanceInM
667
668 ! Locals:
669 real(8) :: dlat, dlon, a, c
670
671 dlon = lon2 - lon1
672 dlat = lat2 - lat1
673
674 a = (sin(dlat/2.d0))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2.d0))**2
675 c = 2.d0 * atan2(sqrt(a),sqrt(1.d0-a))
676 distanceInM = ec_ra * c
677
678 end function calcDistance
679
680 !--------------------------------------------------------------------------
681 ! CALCVERTLOCALIZATIONRADII
682 !--------------------------------------------------------------------------
683 subroutine calcVertLocalizationRadii(ensPerts,stride,waveBandIndex_opt)
684 implicit none
685
686 ! Arguments:
687 type(struct_ens), intent(inout) :: ensPerts
688 integer, intent(in) :: stride
689 integer,optional, intent(in) :: waveBandIndex_opt
690
691 ! Locals:
692 type(struct_gsv) :: statevector_ensStdDev
693 real(8), pointer :: ensStdDev(:,:,:)
694 real(4), pointer :: ptr4d_r4(:,:,:,:)
695 real(4), allocatable :: ensPert_local(:,:)
696 real(8), allocatable :: meanCorrel(:,:)
697 real(8), allocatable :: meanCorrelSquare(:,:)
698 real(8), allocatable :: meanVarianceProduct(:,:)
699 real(8), allocatable :: meanCovarianceSquare(:,:)
700 real(8), allocatable :: meanFourthMoment(:,:)
701 real(8), allocatable :: meanCorrel_local(:,:)
702 real(8), allocatable :: meanCorrelSquare_local(:,:)
703 real(8), allocatable :: meanVarianceProduct_local(:,:)
704 real(8), allocatable :: meanCovarianceSquare_local(:,:)
705 real(8), allocatable :: meanFourthMoment_local(:,:)
706 real(8), allocatable :: localizationFunctions(:,:,:) ! Eq. 19-21 in MMMB 2015 Part 2
707 real(8), allocatable :: localizationRadii(:,:)
708 real(8), allocatable, target :: distanceBinInLnP_T(:,:) ! Distance between each pair of thermo vertical levels in ln(Pressure)
709 real(8), allocatable, target :: distanceBinInLnP_M(:,:) ! Distance between each pair of momentum vertical levels in ln(Pressure)
710 real(8), allocatable :: distanceBinWeight(:) ! Weight given to each bin in the curve fitting step
711 real(8), allocatable :: gridPointWeight(:,:) ! Weight given to grid point in the statistic computation
712 real(8), allocatable :: sumWeight(:,:) ! Sample size for each distance-bin
713 real(8), allocatable :: sumWeight_local(:,:)
714 real(8), pointer :: PressureProfile(:)
715 real(8), pointer :: distanceBinInLnP(:,:)
716 real(8) :: dnens, correlation, covariance, fourthMoment, weight
717 real(8) :: t1, t2, t3, rmse
718 integer :: k, k2, kens, f, ens, bin, numbins, numFunctions
719 integer :: iref, jref
720 integer :: nLevEns, nLevStart, nLevEnd, jvar
721 integer :: myLonBeg, myLonEnd, myLatBeg, myLatEnd, nSize, ier
722 character(len=128) :: outfilename
723 character(len=2) :: wbnum
724
725 !
726 !- 1. Setup
727 !
728 call ens_copyEnsStdDev(ensPerts, statevector_ensStdDev)
729 call gsv_getField(statevector_ensStdDev,ensStdDev)
730
731 numFunctions = 3
732
733 numBins=max(nLevEns_M,nLevEns_T)
734
735 call ens_getLatLonBounds(ensPerts, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
736
737 allocate(meanCorrelSquare(numBins,nkgdimEns))
738 allocate(meanCorrel(numBins,nkgdimEns))
739 allocate(meanVarianceProduct(numBins,nkgdimEns))
740 allocate(meanCovarianceSquare(numBins,nkgdimEns))
741 allocate(meanFourthMoment(numBins,nkgdimEns))
742 allocate(sumWeight(numBins,nkgdimEns))
743
744 allocate(meanCorrelSquare_local(numBins,nkgdimEns))
745 allocate(meanCorrel_local(numBins,nkgdimEns))
746 allocate(meanVarianceProduct_local(numBins,nkgdimEns))
747 allocate(meanCovarianceSquare_local(numBins,nkgdimEns))
748 allocate(meanFourthMoment_local(numBins,nkgdimEns))
749 allocate(sumWeight_local(numBins,nkgdimEns))
750
751 allocate(gridPointWeight(ni,nj))
752
753 allocate(distanceBinInLnP_M(nLevEns_M,nLevEns_M))
754 allocate(distanceBinInLnP_T(nLevEns_T,nLevEns_T))
755
756 ! Compute the separation-distance in ln(p) between vertical levels
757 do k = 1, nLevEns_M
758 do k2 = 1, nLevEns_M
759 distanceBinInLnP_M(k,k2) = abs(log(pressureProfile_M(k))-log(pressureProfile_M(k2)))
760 end do
761 end do
762 do k = 1, nLevEns_T
763 do k2 = 1, nLevEns_T
764 distanceBinInLnP_T(k,k2) = abs(log(pressureProfile_T(k))-log(pressureProfile_T(k2)))
765 end do
766 end do
767
768 dnens = 1.0d0/dble(nens-1)
769
770 ! Grid point Weight
771 do jref = 1, nj
772 gridPointWeight(:,jref) = cos(hco_ens%lat(jref))
773 end do
774
775 !
776 !- 2. Estimation of localization functions
777 !
778
779 !- 2.1 Computation of various statistics for different separation distances
780 meanCorrelSquare_local(:,:) = 0.d0
781 meanCorrel_local(:,:) = 0.d0
782 meanCovarianceSquare_local(:,:) = 0.d0
783 meanVarianceProduct_local(:,:) = 0.d0
784 meanFourthMoment_local(:,:) = 0.d0
785 sumWeight_local(:,:) = 0.d0
786
787 allocate(ensPert_local(nens,nkgdimEns))
788
789 do jref = nint(stride/2.0)+horizPadding, nj-horizPadding, stride ! Pick every stride point to save cost.
790 do iref = nint(stride/2.0)+horizPadding, ni-horizPadding, stride ! Pick every stride point to save cost.
791
792 if (iref < myLonBeg .or. iref > myLonEnd .or. &
793 jref < myLatBeg .or. jref > myLatEnd ) cycle
794
795 !- Select data needed to speed up the process (ensemble member index must come first in ensPert_Local
796 ! because ensemble member is the last loop index below)
797 do ens = 1, nens
798 do k = 1, nkgdimEns
799 ptr4d_r4 => ens_getOneLev_r4(ensPerts,k)
800 ensPert_local(ens,k) = ptr4d_r4(ens,1,iref,jref)
801 end do
802 end do
803
804 ! Loop on all 3D variables
805 do jvar = 1, nvar3d
806
807 if (vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
808 nLevEns = nLevEns_M
809 else
810 nLevEns = nLevEns_T
811 endif
812 nLevStart = varLevOffset(jvar)+ 1
813 nLevEnd = varLevOffset(jvar)+ nLevEns
814
815 !$OMP PARALLEL DO PRIVATE (k,k2,ens,correlation,covariance,fourthMoment,bin,weight)
816 do k = nLevStart, nLevEnd
817 do k2 = nLevStart, nLevEnd
818
819 covariance = 0.d0
820 fourthMoment = 0.d0
821 do ens = 1, nens
822 covariance = covariance + &
823 ensPert_local(ens,k2) * ensPert_local(ens,k)
824 fourthMoment = fourthMoment + &
825 (ensPert_local(ens,k2) * ensPert_local(ens,k))**2
826 end do
827 covariance = covariance * dnens
828 fourthMoment = fourthMoment / real(nens,8)
829 if ( ensStdDev(iref,jref,k2) > 0.d0 .and. ensStdDev(iref,jref,k) > 0.d0 ) then
830 correlation = covariance / (ensStdDev(iref,jref,k2)*ensStdDev(iref,jref,k))
831 else
832 correlation = 0.d0
833 end if
834
835 bin=k2-nLevStart+1
836
837 weight = gridPointWeight(iref,jref)
838
839 sumWeight_local(bin,k) = sumWeight_local(bin,k) + weight
840
841 meanCorrel_local(bin,k) = meanCorrel_local(bin,k) + &
842 correlation * weight
843 meanCorrelSquare_local(bin,k) = meanCorrelSquare_local(bin,k) + &
844 correlation**2 * weight
845 meanCovarianceSquare_local(bin,k) = meanCovarianceSquare_local(bin,k) + &
846 covariance**2 * weight
847 meanVarianceProduct_local(bin,k) = meanVarianceProduct_local(bin,k) + &
848 ensStdDev(iref,jref,k2)**2 * ensStdDev(iref,jref,k)**2 * weight
849 meanFourthMoment_local(bin,k) = meanFourthMoment_local(bin,k) + &
850 fourthMoment * weight
851 end do
852 end do
853 !$OMP END PARALLEL DO
854 end do ! var3D
855
856 end do ! iref
857 end do ! jref
858
859 deallocate(ensPert_local)
860
861 !- 2.2 Gather the all the info in processor 0
862 nSize = nkgdimEns * numbins
863 call rpn_comm_reduce(sumWeight_local ,sumWeight ,nSize, &
864 "mpi_double_precision","mpi_sum",0,"GRID",ier)
865 call rpn_comm_reduce(meanCorrel_local ,meanCorrel ,nSize, &
866 "mpi_double_precision","mpi_sum",0,"GRID",ier)
867 call rpn_comm_reduce(meanCorrelSquare_local ,meanCorrelSquare ,nSize, &
868 "mpi_double_precision","mpi_sum",0,"GRID",ier)
869 call rpn_comm_reduce(meanVarianceProduct_local ,meanVarianceProduct ,nSize, &
870 "mpi_double_precision","mpi_sum",0,"GRID",ier)
871 call rpn_comm_reduce(meanFourthMoment_local ,meanFourthMoment ,nSize, &
872 "mpi_double_precision","mpi_sum",0,"GRID",ier)
873 call rpn_comm_reduce(meanCovarianceSquare_local,meanCovarianceSquare,nSize, &
874 "mpi_double_precision","mpi_sum",0,"GRID",ier)
875
876 deallocate(sumWeight_local)
877 deallocate(meanCorrel_local)
878 deallocate(meanCorrelSquare_local)
879 deallocate(meanVarianceProduct_local)
880 deallocate(meanFourthMoment_local)
881 deallocate(meanCovarianceSquare_local)
882
883 if (mmpi_myid == 0) then
884 !- 2.3 Computation of the localization functions
885 allocate(localizationFunctions(numFunctions,numBins,nkgdimEns))
886
887 t1=dble((nens-1)**2)/dble(nens*(nens-3))
888 t2=dble(nens)/dble((nens-2)*(nens-3))
889 t3=dble(nens-1)/dble(nens*(nens-2)*(nens-3))
890
891 do jvar = 1, nvar3d
892 if (vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
893 nLevEns = nLevEns_M
894 else
895 nLevEns = nLevEns_T
896 endif
897 nLevStart = varLevOffset(jvar)+ 1
898 nLevEnd = varLevOffset(jvar)+ nLevEns
899 !$OMP PARALLEL DO PRIVATE (k,bin)
900 do k = nLevStart, nLevEnd
901 do bin = 1, nLevEns
902 meanCorrel(bin,k) = meanCorrel(bin,k) / sumWeight(bin,k)
903 meanCorrelSquare(bin,k) = meanCorrelSquare(bin,k) / sumWeight(bin,k)
904 meanCovarianceSquare(bin,k) = meanCovarianceSquare(bin,k) / sumWeight(bin,k)
905 meanVarianceProduct(bin,k) = meanVarianceProduct(bin,k) / sumWeight(bin,k)
906 meanFourthMoment(bin,k) = meanFourthMoment(bin,k) / sumWeight(bin,k)
907
908 if ( meanCovarianceSquare(bin,k) /= 0.d0 ) then
909 ! Form 1: General formulation (Eq. 19 in MMMB 2015 Part 2)
910 localizationFunctions(1,bin,k) = t1 - t2*meanFourthMoment(bin,k)/meanCovarianceSquare(bin,k) + &
911 t3*meanVarianceProduct(bin,k)/meanCovarianceSquare(bin,k)
912 ! Form 2: Gaussian sample distribution (Eq. 20 in MMMB 2015 Part 2)
913 localizationFunctions(2,bin,k) = dble(nens-1)/dble((nens+1)*(nens-2)) * &
914 (dble(nens-1)-meanVarianceProduct(bin,k)/meanCovarianceSquare(bin,k))
915 else
916 write(*,*) ' !!! Warning !!! meanCovarianceSquare = 0 in bin, level = ', bin, k
917 localizationFunctions(1,bin,k) = 0.d0
918 localizationFunctions(2,bin,k) = 0.d0
919 end if
920 ! Form 3: Gaussian sample distribution and correlation-based formulation (Eq. 21 in MMMB 2015 Part 2)
921 if ( meanCorrelSquare(bin,k) /= 0.d0 ) then
922 localizationFunctions(3,bin,k) = dble(nens-1)/dble((nens+1)*(nens-2)) * &
923 (dble(nens-1)-1.d0/meanCorrelSquare(bin,k))
924 else
925 write(*,*) ' !!! Warning !!! meanCorrelSquare = 0 in bin, level = ', bin, k
926 localizationFunctions(3,bin,k) = 0.d0
927 end if
928 end do
929 end do
930 !$OMP END PARALLEL DO
931 end do
932
933 !
934 !- 3. Estimation of localization radii (curve fitting)
935 !
936 allocate(localizationRadii(numFunctions,nkgdimEns))
937 allocate(distanceBinWeight(numBins))
938
939 call lfn_setup('FifthOrder') ! IN
940
941 localizationRadii(:,:) = 2.d0 ! First Guess (in ln(p) distance)
942 distanceBinWeight(:) = 1.d0 ! Even weight
943
944 do jvar = 1, nvar3d
945 if (vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
946 nLevEns = nLevEns_M
947 distanceBinInLnP => distanceBinInLnP_M
948 else
949 nLevEns = nLevEns_T
950 distanceBinInLnP => distanceBinInLnP_T
951 endif
952 nLevStart = varLevOffset(jvar)+ 1
953 nLevEnd = varLevOffset(jvar)+ nLevEns
954
955 write(*,*)
956 write(*,*) nomvar3d(jvar)
957
958 do f = 1, numFunctions
959
960 write(*,*)
961 write(*,*) 'Localization Function : ', f
962 do k = nLevStart, nLevEnd
963 kens = k-nLevStart+1
964 write(*,*) ' ----- EnsLev : ', k
965 call lfn_lengthscale( localizationRadii(f,k), & ! INOUT
966 rmse, & ! OUT
967 localizationFunctions(f,1:nLevEns,k), & ! IN
968 distanceBinInLnP(kens,:), & ! IN
969 distanceBinWeight(1:nLevEns), & ! IN
970 nLevEns ) ! IN
971 end do
972 end do
973 end do
974
975 !
976 !- 4. Write to file
977 !
978 if ( nWaveBand /= 1 ) then
979 if (.not. present(waveBandIndex_opt)) then
980 write(*,*) 'calcLocalizationRadii: No waveBandIndex was supplied!!!'
981 call utl_abort('calbmatrix_glb')
982 end if
983 write(wbnum,'(I2.2)') waveBandIndex_opt
984 end if
985
986 !- 4.1 Localization functions in txt format (for plotting purposes)
987 do jvar = 1, nvar3d
988 if ( nWaveBand == 1 ) then
989 outfilename = "./vertLocalizationFunctions_"//trim(nomvar3d(jvar))//".txt"
990 else
991 outfilename = "./vertLocalizationFunctions_"//trim(nomvar3d(jvar))//"_"//wbnum//".txt"
992 end if
993 open (unit=99,file=outfilename,action="write",status="new")
994
995 if(vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
996 nLevEns = nLevEns_M
997 PressureProfile => pressureProfile_M
998 distanceBinInLnP => distanceBinInLnP_M
999 else
1000 nLevEns = nLevEns_T
1001 PressureProfile => pressureProfile_T
1002 distanceBinInLnP => distanceBinInLnP_T
1003 endif
1004
1005 do k=1,nlevEns
1006 do bin = 1, nlevEns
1007 write(99,'(I3,2X,I3,2X,F7.2,2X,F7.3,2X,I7,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,F6.4,2X,E10.3,2X,E10.3,2X,E10.3)') k, bin, &
1008 PressureProfile(bin)/100.d0, distanceBinInLnP(k,bin), nint(sumWeight(bin,varLevOffset(jvar)+k)), &
1009 meanCorrel(bin,varLevOffset(jvar)+k), &
1010 localizationFunctions(3,bin,varLevOffset(jvar)+k), &
1011 localizationFunctions(2,bin,varLevOffset(jvar)+k), &
1012 localizationFunctions(1,bin,varLevOffset(jvar)+k), &
1013 meanCorrelSquare(bin,varLevOffset(jvar)+k), &
1014 meanCovarianceSquare(bin,varLevOffset(jvar)+k), &
1015 meanVarianceProduct(bin,varLevOffset(jvar)+k), &
1016 meanFourthMoment(bin,varLevOffset(jvar)+k)
1017 end do
1018 end do
1019 close(unit=99)
1020 end do
1021
1022 !- 4.2 Localization radii in txt format (for plotting purposes)
1023 do jvar = 1, nvar3d
1024 if ( nWaveBand == 1 ) then
1025 outfilename = "./vertLocalizationRadii_"//trim(nomvar3d(jvar))//".txt"
1026 else
1027 outfilename = "./vertLocalizationRadii_"//trim(nomvar3d(jvar))//"_"//wbnum//".txt"
1028 end if
1029 open (unit=99,file=outfilename,action="write",status="new")
1030
1031 if(vnl_varLevelFromVarName(nomvar3d(jvar)).eq.'MM') then
1032 nLevEns = nLevEns_M
1033 PressureProfile => pressureProfile_M
1034 else
1035 nLevEns = nLevEns_T
1036 PressureProfile => pressureProfile_T
1037 endif
1038 do k=1,nlevEns
1039 write(99,'(I3,2X,F7.2,2X,F7.2,2X,F7.2,2X,F7.2)') k, PressureProfile(k)/100.d0, &
1040 localizationRadii(3,varLevOffset(jvar)+k), &
1041 localizationRadii(2,varLevOffset(jvar)+k), &
1042 localizationRadii(1,varLevOffset(jvar)+k)
1043 end do
1044 close(unit=99)
1045 end do
1046
1047 deallocate(localizationFunctions)
1048 deallocate(localizationRadii)
1049 deallocate(distanceBinWeight)
1050
1051 end if ! mmpi_myid == 0
1052
1053 deallocate(meanCorrelSquare)
1054 deallocate(meanCorrel)
1055 deallocate(meanCovarianceSquare)
1056 deallocate(meanVarianceProduct)
1057 deallocate(meanFourthMoment)
1058 deallocate(distanceBinInLnP_M)
1059 deallocate(distanceBinInLnP_T)
1060 deallocate(sumWeight)
1061 deallocate(gridPointWeight)
1062
1063 write(*,*) 'finished estimating the vertical localization radii...'
1064
1065 end subroutine calcVertLocalizationRadii
1066
1067end module menetrierDiag_mod
