1module calcStatsGlb_mod
2 ! MODULE calcStatsGlb_mod (prefix='csg' 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 (global
6 ! version).
7 !
8 use codePrecision_mod
9 use midasMpi_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use ensembleStateVector_mod
13 use globalSpectralTransform_mod
14 use mathPhysConstants_mod
15 use horizontalCoord_mod
16 use verticalCoord_mod
17 use varNameList_mod
18 use calcHeightAndPressure_mod
19 use earthConstants_mod
20 use utilities_mod
21 use spectralFilter_mod
22 use menetrierDiag_mod
23 use fileNames_mod
24 use timeCoord_mod
25 use gridVariableTransforms_mod
26 use gridBinning_mod
27 use verticalModes_mod
28 implicit none
29 save
30 private
31
32 ! Public Subroutines
33 public :: csg_setup, csg_computeBhi, csg_toolbox
34
35 type(struct_hco), pointer :: hco_ens ! Ensemble horizontal grid parameters
36 type(struct_vco), pointer :: vco_ens ! Ensemble horizontal grid parameters
37
38 integer :: nens,ni,nj,nLevEns_M,nLevEns_T,nLevPtoT,nkgdimEns,varLevOffset(6),nla
39 integer :: myLatBeg, myLatEnd, myLonBeg, myLonEnd, latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
40 integer :: mymBeg, mymEnd, mymSkip, mymCount, mynBeg, mynEnd, mynSkip, mynCount
41 integer :: myMemberBeg, myMemberEnd, myMemberCount, maxMyMemberCount
42 integer :: nEnsOverDimension
43 integer :: nla_mpilocal, maxMyNla
44 integer, pointer :: ilaList_mpiglobal(:), ilaList_mpilocal(:)
45
46 character(len=256), allocatable :: cflensin(:)
47 integer :: gstID_nkgdimEns, gstID_nLevEns_M, gstID_nLevEns_T_P1
48 integer, allocatable :: nip1_M(:),nip1_T(:)
49 real(8), pointer :: pressureProfile_M(:), pressureProfile_T(:)
50 integer,external :: get_max_rss
51
52 integer,parameter :: nvar3d=4, nvar2d=1, nvar=nvar3d+nvar2d
53 character(len=4) :: nomvar3d(nvar3d,3), nomvar2d(nvar2d,3), nomvar(nvar,3)
54 integer, parameter :: modelSpace = 1
55 integer, parameter :: cvSpace = 2
56 integer, parameter :: cvUnbalSpace = 3
57
58 real(8) :: gridSpacingInKm
59
60 ! For wave band decomposition
61 integer, parameter :: maxNumLocalLength = 20
62 integer :: nWaveBand
63 integer :: nVertWaveBand
64
65 logical :: initialized = .false.
66
67 ! Namelist variables
68 integer :: ntrunc
69 integer :: waveBandPeaks(maxNumLocalLength) ! For horizontal wave band decomposition
70 integer :: vertWaveBandPeaks(maxNumLocalLength) ! For vertical wave band decomposition
71
72 contains
73
74 !--------------------------------------------------------------------------
75 ! CSG_SETUP
76 !--------------------------------------------------------------------------
77 subroutine csg_setup(nens_in, hco_in, vco_in)
78 !
79 !:Purpose: Main setup routine for this module
80 !
81 implicit none
82
83 ! Arguments:
84 integer, intent(in) :: nens_in
85 type(struct_vco), pointer, intent(in) :: vco_in
86 type(struct_hco), pointer, intent(in) :: hco_in
87
88 ! Locals:
89 integer :: nulnam, ierr, waveBandIndex, memberIndex
90 integer :: fclos, fnom
91 real(8) :: zps
92
93 NAMELIST /NAMCALCSTATS_GLB/ntrunc,waveBandPeaks,vertWaveBandPeaks
94
95 write(*,*)
96 write(*,*) 'csg_setup: Starting...'
97
98 nens=nens_in
99 allocate(cflensin(nens))
100 do memberIndex = 1, nEns
101 call fln_ensfileName(cflensin(memberIndex), 'ensemble', memberIndex_opt=memberIndex)
102 end do
103
104 if ( mmpi_myid == 0 ) then
105 call mpc_printConstants(6)
106 call pre_printPrecisions
107 end if
108
109 ! parameters from namelist (date in filename should come directly from sequencer?)
110 ntrunc=108
111 waveBandPeaks(:) = -1.0d0
112 vertWaveBandPeaks(:) = -1.0d0
113
114 nulnam=0
115 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
116 read(nulnam,nml=NAMCALCSTATS_GLB)
117 if (ierr /= 0) call utl_abort('csg_setup: Error reading namelist NAMCALCSTATS_GLB')
118 if (mmpi_myid == 0) write(*,nml=NAMCALCSTATS_GLB)
119 ierr=fclos(nulnam)
120
121 !- Setup horizontal grid
122 hco_ens => hco_in
123 ni=hco_in%ni
124 nj=hco_in%nj
125
126 gridSpacingInKm = ec_ra * hco_in%dlon / 1000.d0
127 write(*,*) 'Grid Spacing in Km = ', gridSpacingInKm
128
129 ! setup mpi local grid parameters
130 call mmpi_setup_latbands(nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
131 call mmpi_setup_lonbands(ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
132
133 !- Setup vertical levels
134 vco_ens => vco_in
135 nLevEns_M=vco_in%nlev_M
136 nLevEns_T=vco_in%nlev_T
137 nLevPtot=nLevEns_M-1 ! ignore streamfunction at hyb=1, since highly correlated with next level
138 varLevOffset(1) = 0
139 varLevOffset(2) = 1*nLevEns_M
140 varLevOffset(3) = 2*nLevEns_M
141 varLevOffset(4) = 2*nLevEns_M+1*nLevEns_T
142 varLevOffset(5) = 2*nLevEns_M+2*nLevEns_T
143 nkgdimEns=nLevEns_M*2+nLevEns_T*2+1 ! NO TG !!!
144 nla=(ntrunc+1)*(ntrunc+2)/2
145
146 !- Setup the global spectral transform
147 gstID_nkgdimEns = gst_setup(ni,nj,ntrunc,nkgdimEns)
148 gstID_nLevEns_M = gst_setup(ni,nj,ntrunc,nLevEns_M)
149 gstID_nLevEns_T_P1 = gst_setup(ni,nj,ntrunc,nLevEns_T+1)
150
151 ! setup mpi local spectral vector parameters
152 call mmpi_setup_m(ntrunc, mymBeg, mymEnd, mymSkip, mymCount)
153 call mmpi_setup_n(ntrunc, mynBeg, mynEnd, mynSkip, mynCount)
154 call gst_ilaList_mpiglobal(ilaList_mpiglobal, nla_mpilocal, maxMyNla, &
155 gstID_nkgdimEns, mymBeg, mymEnd, mymSkip, mynBeg, mynEnd, mynSkip)
156 call gst_ilaList_mpilocal(ilaList_mpilocal, &
157 gstID_nkgdimEns, mymBeg, mymEnd, mymSkip, mynBeg, mynEnd, mynSkip)
158
159 ! setup ensemble members mpi partinionning (when working with struct_ens)
160 call mmpi_setup_levels(nEns,myMemberBeg,myMemberEnd,myMemberCount)
161 call rpn_comm_allreduce(myMemberCount, maxMyMemberCount, &
162 1,"MPI_INTEGER","MPI_MAX","GRID",ierr)
163 nEnsOverDimension = mmpi_npex * maxMyMemberCount
164
165 !- Setup ip1s
166 allocate(nip1_M(nLevEns_M))
167 nip1_M(:)=vco_in%ip1_M(:)
168 allocate(nip1_T(nLevEns_T))
169 nip1_T(:)=vco_in%ip1_T(:)
170
171 !- Estimate the pressure profile for each vertical grid
172 zps = 101000.D0
173 call czp_fetch1DLevels(vco_in, zps, &
174 profM_opt=pressureProfile_M, profT_opt=pressureProfile_T)
175
176 !
177 !- Setup variable names
178 !
179 nomvar3d(1,modelSpace)='UU'
180 nomvar3d(2,modelSpace)='VV'
181 nomvar3d(3,modelSpace)='TT'
182 nomvar3d(4,modelSpace)='LQ'
183 nomvar2d(1,modelSpace)='P0'
184
185 nomvar3d(1,cvSpace)='PP'
186 nomvar3d(2,cvSpace)='CC'
187 nomvar3d(3,cvSpace)='TT'
188 nomvar3d(4,cvSpace)='LQ'
189 nomvar2d(1,cvSpace)='P0'
190
191 nomvar3d(1,cvUnbalSpace)='PP'
192 nomvar3d(2,cvUnbalSpace)='UC'
193 nomvar3d(3,cvUnbalSpace)='UT'
194 nomvar3d(4,cvUnbalSpace)='LQ'
195 nomvar2d(1,cvUnbalSpace)='UP'
196
197 nomvar(1:4,:) = nomvar3d(1:4,:)
198 nomvar(5,:) = nomvar2d(1,:)
199
200 !
201 !- Horizontal wave band decomposition option
202 !
203 nWaveBand = count(waveBandPeaks >= 0)
204 if ( nWaveBand < 1 ) then
205 nWaveBand = 1
206 else if (nWaveBand == 1) then
207 write(*,*) 'You have specified only ONE waveBandPeaks'
208 call utl_abort('calbmatrix_glb')
209 else
210 write(*,*)
211 write(*,*) 'Horizontal waveBand decomposition is ACTIVATED'
212 end if
213
214 ! Make sure that the wavenumbers are in the correct (decreasing) order
215 do waveBandIndex = 1, nWaveBand-1
216 if ( waveBandPeaks(waveBandIndex)-waveBandPeaks(waveBandIndex+1) <= 0 ) then
217 write(*,*) 'csg_setup: waveBandPeaks are not in decreasing wavenumber order'
218 call utl_abort('calbmatrix_glb')
219 end if
220 end do
221
222 ! Make sure the truncation is compatible with the waveBandPeaks
223 if ( ntrunc < nj-1 .and. nWaveBand > 1 ) then
224 write(*,*) 'csg_setup: The truncation is not compatible with wave band decomposition'
225 write(*,*) ' ntrunc should = ', nj-1
226 call utl_abort('calbmatrix_glb')
227 end if
228
229 !
230 !- Vertical wave band decomposition option
231 !
232 nVertWaveBand = count(vertWaveBandPeaks >= 0)
233 if ( nVertWaveBand < 1 ) then
234 nVertWaveBand = 1
235 else if (nVertWaveBand == 1) then
236 write(*,*) 'You have specified only ONE waveBandPeaks'
237 call utl_abort('calbmatrix_glb')
238 else
239 write(*,*)
240 write(*,*) 'Vertical waveBand decomposition is ACTIVATED'
241 end if
242
243 ! Make sure that the wavenumbers are in the correct (decreasing) order
244 do waveBandIndex = 1, nVertWaveBand-1
245 if ( vertWaveBandPeaks(waveBandIndex)-vertWaveBandPeaks(waveBandIndex+1) <= 0 ) then
246 write(*,*) 'csg_setup: vertWaveBandPeaks are not in decreasing mode order'
247 call utl_abort('calbmatrix_glb')
248 end if
249 end do
250
251 !
252 !- Ending
253 !
254 initialized = .true.
255
256 write(*,*)
257 write(*,*) 'csg_setup: Done!'
258
259 end subroutine csg_setup
260
261 !--------------------------------------------------------------------------
262 ! csg_computeBhi
263 !--------------------------------------------------------------------------
264 subroutine csg_computeBhi()
265 !
266 !:Purpose: Master routine for Bhi computation in global mode
267 !
268 implicit none
269
270 ! Locals:
271 integer :: ierr, nulnam
272 integer :: fclos, fnom
273 character(len=12) :: formulation ! Bhi formulation
274
275 NAMELIST /NAMCOMPUTEBHI/formulation
276
277 ! parameters from namelist
278 formulation='legacy'
279
280 nulnam=0
281 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
282 read(nulnam,nml=NAMCOMPUTEBHI)
283 if (ierr /= 0) call utl_abort('csg_computeBhi: Error reading namelist NAMCOMPUTEBHI')
284 if (mmpi_myid == 0) write(*,nml=NAMCOMPUTEBHI)
285 ierr=fclos(nulnam)
286
287 select case(trim(formulation))
288 case ('legacy')
289 call csg_computeBhiLegacy
290 case ('latbands')
291 call csg_computeBhiLatBands
292 case default
293 write(*,*)
294 write(*,*) 'csg_computeBhi: Unknown value of FORMULATION for Bhi computation: ',formulation
295 write(*,*) 'Please select legacy or latbands'
296 call utl_abort('csg_computeBhi')
297 end select
298
299 end subroutine csg_computeBhi
300
301 !--------------------------------------------------------------------------
302 ! csg_computeBhiLegacy
303 !--------------------------------------------------------------------------
304 subroutine csg_computeBhiLegacy
305 !
306 !:Purpose: Computation of Bhi using the legacy formulation
307 !
308 implicit none
309
310 ! Locals:
311 real(4), pointer :: ensPerturbations(:,:,:,:), ens_ptr(:,:,:,:)
312 real(4), pointer :: ensBalPerturbations(:,:,:,:)
313 real(8), pointer :: stddev3d(:,:,:), stddev3dBal(:,:,:), stddev3dUnbal(:,:,:)
314 real(8), pointer :: stddev3d_ptr(:,:,:)
315 real(8), pointer :: stddevZonAvg(:,:), stddevZonAvgBal(:,:), stddevZonAvgUnbal(:,:)
316 real(8), allocatable :: PtoT(:,:,:),theta1(:,:),theta2(:,:)
317 real(8), allocatable :: corns(:,:,:),rstddev(:,:)
318 integer :: variableType
319
320 allocate(ensPerturbations(myLonBeg:myLonEnd, myLatBeg:myLatEnd, nkgdimEns, nens))
321 allocate(ensBalPerturbations(myLonBeg:myLonEnd, myLatBeg:myLatEnd, nLevEns_T+1, nens))
322 allocate(theta1(nlevEns_M,nj))
323 allocate(stddev3d(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdimEns))
324 allocate(stddevZonAvg(nkgdimEns,nj))
325 allocate(PtoT(nlevEns_T+1,nlevEns_M,nj))
326 allocate(theta2(nlevEns_M,nj))
327 allocate(stddev3dBal(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLevEns_T+1))
328 allocate(stddev3dUnbal(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdimEns))
329 allocate(stddevZonAvgBal(nLevEns_T+1,nj))
330 allocate(stddevZonAvgUnbal(nkgdimEns,nj))
331 allocate(corns(nkgdimEns,nkgdimEns,0:ntrunc))
332 allocate(rstddev(nkgdimEns,0:ntrunc))
333
334 write(*,*) 'Initializing ensemble arrays to claim memory'
335 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
336 ensPerturbations(:,:,:,:) = 0.0
337 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
338 ensBalPerturbations(:,:,:,:) = 0.0
339 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
340
341 call readEnsemble(ensPerturbations)
342 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
343
344 call removeMean(ensPerturbations)
345
346 call uv_to_psichi(ensPerturbations)
347
348 call calcStddev3d(ensPerturbations,stddev3d,nkgdimens)
349
350 call calcZonAvg(stddevZonAvg,stddev3d,nkgdimens)
351
352 call calcTheta(ensPerturbations,theta1) ! theta1 is put in glbcov and used for analysis!
353 if (mmpi_myid == 0) write(301,*) theta1
354
355 call removeBalancedChi(ensPerturbations,theta1)
356
357 call normalize3d(ensPerturbations,stddev3d)
358
359 call calcPtoT(ensPerturbations,PtoT)
360 if (mmpi_myid == 0) write(303,*) PTOT(:,:,1)
361 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
362
363! call calcTheta(ensPerturbations,theta2) ! theta2 is used previously for computing unbalanced Chi!
364! if (mmpi_myid == 0) write(302,*) theta2
365
366 call removeBalancedT_Ps(ensPerturbations,ensBalPerturbations,PtoT)
367 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
368
369! call removeBalancedChi(ensPerturbations,theta2)
370
371 call multiply3d(ensPerturbations,stddev3d,nkgdimens)
372
373 ens_ptr(myLonBeg:,myLatBeg:,1:,1:) => ensBalPerturbations(:,:,1:nLevEns_T,:)
374 stddev3d_ptr(myLonBeg:,myLatBeg:,1:) => stddev3d(:,:,(2*nLevEns_M+1):(2*nLevEns_M+nLevEns_T))
375 call multiply3d(ens_ptr, stddev3d_ptr, nLevEns_T)
376
377 ens_ptr(myLonBeg:,myLatBeg:,1:,1:) => ensBalPerturbations(:,:,(nLevEns_T+1):(nLevEns_T+1),:)
378 stddev3d_ptr(myLonBeg:,myLatBeg:,1:) => stddev3d(:,:,(2*nLevEns_M+2*nLevEns_T+1):(2*nLevEns_M+2*nLevEns_T+1))
379 call multiply3d(ens_ptr, stddev3d_ptr, 1)
380
381 call spectralFilter(ensPerturbations,nkgdimens)
382
383 call spectralFilter(ensBalPerturbations,nLevEns_T+1)
384
385 call calcStddev3d(ensPerturbations,stddev3dUnbal,nkgdimens)
386
387 call calcStddev3d(ensBalPerturbations,stddev3dBal,nLevEns_T+1)
388
389 call calcZonAvg(stddevZonAvgUnbal,stddev3dUnbal,nkgdimens)
390
391 call calcZonAvg(stddevZonAvgBal,stddev3dBal,nLevEns_T+1)
392
393 call normalize3d(ensPerturbations,stddev3dUnbal)
394
395 call removeGlobalMean(ensPerturbations)
396
397 call calcCorrelations(ensPerturbations,corns,rstddev)
398
399 variableType = cvUnbalSpace
400 call writeStats(corns,rstddev,ptot,theta1)
401
402 call writeStddev(stddevZonAvg,stddev3d,stddevZonAvgUnbal,stddev3dUnbal)
403
404 call writeStddevBal(stddevZonAvgBal,stddev3dBal)
405
406 call writeSpStats(ptot,theta1)
407 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
408
409 if (mmpi_myid == 0) then
410 write(200,*) stddevZonAvg(1:nlevEns_M,:)
411 write(201,*) stddevZonAvg((1+1*nlevEns_M):(2*nlevEns_M),:)
412 write(202,*) stddevZonAvg((1+2*nlevEns_M):(3*nlevEns_T),:)
413 write(203,*) stddevZonAvg((1+2*nlevEns_M+1*nlevEns_T):(2*nlevEns_M+2*nlevEns_T),:)
414 write(204,*) stddevZonAvg((1+2*nlevEns_M+2*nlevEns_T),:)/1.0d2
415
416 write(400,*) stddevZonAvgUnbal(1:nlevEns_M,:)
417 write(401,*) stddevZonAvgUnbal((1+1*nlevEns_M):(2*nlevEns_M),:)
418 write(402,*) stddevZonAvgUnbal((1+2*nlevEns_M):(3*nlevEns_T),:)
419 write(403,*) stddevZonAvgUnbal((1+2*nlevEns_M+1*nlevEns_T):(2*nlevEns_M+2*nlevEns_T),:)
420 write(404,*) stddevZonAvgUnbal((1+2*nlevEns_M+2*nlevEns_T),:)/1.0d2
421 end if
422
423 end subroutine csg_computeBhiLegacy
424
425 !--------------------------------------------------------------------------
426 ! csg_computeBhiLatBands
427 !--------------------------------------------------------------------------
428 subroutine csg_computeBhiLatBands
429 !
430 !:Purpose: Computation of Bhi on a set of latitude bands
431 !
432 implicit none
433
434 ! Locals:
435 integer :: variableType, latIndex, jlatband, lat1, lat2, lat3
436 real(4), pointer :: ensPerturbations(:,:,:,:)
437 real(8), pointer :: stddev3d(:,:,:)
438 real(8), pointer :: stddevZonAvg(:,:)
439 real(8), allocatable :: corns(:,:,:),rstddev(:,:)
440 real(8) :: latMask(nj)
441
442 allocate(ensPerturbations(ni,nj,nkgdimEns,nens))
443 allocate(stddev3d(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdimEns))
444 allocate(stddevZonAvg(nkgdimEns,nj))
445 allocate(corns(nkgdimEns,nkgdimEns,0:ntrunc))
446 allocate(rstddev(nkgdimEns,0:ntrunc))
447
448 call readEnsemble(ensPerturbations)
449 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
450
451 call removeMean(ensPerturbations)
452
453 call uv_to_psichi(ensPerturbations)
454
455 call calcStddev3d(ensPerturbations,stddev3d,nkgdimens)
456
457 call calcZonAvg(stddevZonAvg,stddev3d,nkgdimens)
458
459 call normalize3d(ensPerturbations,stddev3d)
460
461 call removeGlobalMean(ensPerturbations)
462
463 do jlatband = 1, 3
464 write(*,*) 'csg_computeBhiLatBands: selected LATBAND = ',jlatband
465 lat1=nj/4
466 lat2=nj/2
467 lat3=3*nj/4
468 write(*,*) 'lat1,2,3=',lat1,lat2,lat3
469 if(jlatband==1) then
470 ! Southern extratropics
471 latMask(1:lat1) = 1.0d0
472 do latIndex = lat1, lat2
473 !latMask(latIndex) = sqrt(dble((lat2-latIndex)*4)/dble(nj))
474 latMask(latIndex) = sqrt(0.5d0*(1.0d0+cos(dble((latIndex-lat1)*4)*MPC_PI_R8/dble(nj))))
475 end do
476 latMask(lat2:nj) = 0.0d0
477 else if(jlatband==2) then
478 ! Tropics
479 !latMask(1:lat1) = 0.0d0
480 !do latIndex = lat1, lat2
481 ! latMask(latIndex) = sqrt(dble((latIndex-lat1)*4)/dble(nj))
482 !end do
483 !do latIndex = lat2,lat3
484 ! latMask(latIndex) = sqrt(dble((lat3-latIndex)*4)/dble(nj))
485 !end do
486 !latMask(lat3:nj) = 0.0d0
487
488 ! NOTE: use much broader band for tropics to avoid shortening horizontal correlations
489 ! ok, since the masks do not have to sum to one for calculation, but they do
490 ! when used in bmatrixhi_mod
491 ! Tropics
492 do latIndex = 1, lat1
493 !latMask(latIndex) = sqrt(dble((latIndex-1)*4)/dble(nj))
494 latMask(latIndex) = sqrt(0.5d0*(1.0d0+cos(dble((lat1-latIndex)*4)*MPC_PI_R8/dble(nj))))
495 end do
496 latMask(lat1:lat3) = 1.0d0
497 do latIndex = lat3,nj
498 !latMask(latIndex) = sqrt(dble((nj-latIndex)*4)/dble(nj))
499 latMask(latIndex) = sqrt(0.5d0*(1.0d0+cos(dble((latIndex-lat3)*4)*MPC_PI_R8/dble(nj))))
500 end do
501 else if(jlatband==3) then
502 ! Northern extratropics
503 latMask(1:lat2) = 0.0d0
504 do latIndex = lat2, lat3
505 !latMask(latIndex) = sqrt(dble((latIndex-lat2)*4)/dble(nj))
506 latMask(latIndex) = sqrt(0.5d0*(1.0d0+cos(dble((lat3-latIndex)*4)*MPC_PI_R8/dble(nj))))
507 end do
508 latMask(lat3:nj) = 1.0d0
509 end if
510 write(*,*) 'latMask = ',latMask(:)
511 call calcCorrelations(ensPerturbations,corns,rstddev,latMask_opt=latMask)
512
513 variableType = cvSpace
514 call writeStats(corns,rstddev,latBand_opt=jlatBand)
515 end do
516
517 call writeStddev(stddevZonAvg,stddev3d)
518
519 if (mmpi_myid == 0) then
520 write(200,*) stddevZonAvg(1:nlevEns_M,:)
521 write(201,*) stddevZonAvg((1+1*nlevEns_M):(2*nlevEns_M),:)
522 write(202,*) stddevZonAvg((1+2*nlevEns_M):(3*nlevEns_T),:)
523 write(203,*) stddevZonAvg((1+2*nlevEns_M+1*nlevEns_T):(2*nlevEns_M+2*nlevEns_T),:)
524 write(204,*) stddevZonAvg((1+2*nlevEns_M+2*nlevEns_T),:)/1.0d2
525 end if
526
527 end subroutine csg_computeBhiLatBands
528
529 !--------------------------------------------------------------------------
530 ! CSG_TOOLBOOX
531 !--------------------------------------------------------------------------
532 subroutine csg_toolbox
533 !
534 !:Purpose: High-level routine to do a variety of diagnostic operations
535 ! on the ensemble of error samples
536 !
537 implicit none
538
539 ! NOTE: The diagnostic computed here are in model variable space
540 ! (no variable transform!!!)
541
542 ! Locals:
543 integer :: waveBandIndex
544 integer :: nulnam, ierr, fclos, fnom, numStep
545 real(8), allocatable :: corns(:,:,:), rstddev(:,:), powerSpec(:,:)
546 integer, allocatable :: dateStampList(:)
547 type(struct_ens), target :: ensPerts, ensPertsFilt
548 type(struct_ens), pointer :: ensPerts_ptr
549 type(struct_gsv) :: statevector_template
550 type(struct_gbi) :: gbi_zonalMean
551 type(struct_gbi) :: gbi_globalMean
552 type(struct_vms) :: vModes
553 integer :: variableType
554 logical :: ensContainsFullField
555 logical :: makeBiPeriodic
556 logical :: doSpectralFilter
557 real(8) :: vertModesLengthScale(2)
558 real(8) :: lengthScaleTop, lengthScaleBot
559 character(len=60) :: tool
560 character(len=2) :: wbnum
561 character(len=2) :: ctrlVarHumidity
562
563 NAMELIST /NAMTOOLBOX/tool, ensContainsFullField, doSpectralFilter, vertModesLengthScale, &
564 ctrlVarHumidity
565
566 write(*,*)
567 write(*,*) 'csg_toolbox'
568 write(*,*)
569
570 !
571 !- Set options
572 !
573 variableType = modelSpace ! hardwired
574 ensContainsFullField = .true. ! default value
575 doSpectralFilter = .true.
576 vertModesLengthScale(1) = 6.d0
577 vertModesLengthScale(2) = -1.d0
578 ctrlVarHumidity = 'HU'
579
580 nulnam = 0
581 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
582 read(nulnam,nml=NAMTOOLBOX)
583 if (ierr /= 0) call utl_abort('csg_toolbox: Error reading namelist NAMTOOLBOX')
584 if (mmpi_myid == 0) write(*,nml=NAMTOOLBOX)
585 ierr = fclos(nulnam)
586
587 lengthScaleTop = vertModesLengthScale(1)
588 if (vertModesLengthScale(2) == -1.d0) then
589 lengthScaleBot = lengthScaleTop
590 else
591 lengthScaleBot = vertModesLengthScale(2)
592 end if
593
594 !
595 !- Read ensemble
596 !
597 numStep = 1
598 allocate(dateStampList(numStep))
599 dateStampList(:) = -1
600 call ens_allocate(ensPerts, nEns, numStep, hco_ens, vco_ens, dateStampList)
601
602 if (nWaveBand > 1 .or. nVertWaveBand > 1) then
603 call ens_allocate(ensPertsFilt, nEns, numStep, hco_ens, vco_ens, dateStampList)
604 ensPerts_ptr => ensPertsFilt
605 else
606 ensPerts_ptr => ensPerts
607 end if
608
609 makeBiPeriodic = .false.
610 call ens_readEnsemble(ensPerts, './ensemble', makeBiPeriodic, &
611 containsFullField_opt=ensContainsFullField)
612
613 if ( ctrlVarHumidity == 'LQ' .and. ensContainsFullField ) then
614 call gvt_transform(ensPerts,'HUtoLQ')
615 call ens_modifyVarName(ensPerts, 'HU', 'LQ')
616 if (nWaveBand > 1 .or. nVertWaveBand > 1) then
617 call ens_modifyVarName(ensPerts_ptr, 'HU', 'LQ')
618 end if
619 end if
620
621 !
622 !- Compute and remove the ensemble mean; compute the stdDev
623 !
624 call ens_computeMean(ensPerts)
625 call ens_removeMean (ensPerts)
626
627 !
628 !- Tool selection
629 !
630 select case(trim(tool))
631 case ('HVCORREL_HI')
632 write(*,*)
633 write(*,*) 'Computing Homogeneous and Isotropic Correlation'
634
635 if (mmpi_nprocs > 1) then
636 call utl_abort('csg_toolbox: this tool is not yet MPI capable') ! only due to horizCorrelFunction
637 end if
638
639 call spectralFilter2(ensPerts, & ! IN
640 ensPerts_ptr, & ! OUT
641 waveBandIndex_opt=1) ! IN
642 call ens_computeStdDev(ensPerts_ptr)
643 call ens_normalize(ensPerts_ptr)
644 call ens_removeGlobalMean(ensPerts_ptr)
645
646 allocate(corns(nkgdimEns,nkgdimEns,0:ntrunc))
647 allocate(rstddev(nkgdimEns,0:ntrunc))
648
649 call calcCorrelations2(ensPerts_ptr, & ! IN
650 corns, & ! OUT (vertical correlation in spectral space)
651 rstddev) ! OUT ( sqrt(normalized power spectrum) )
652
653 call writeStats(corns,rstddev,waveBandIndex_opt=1) ! IN
654 call calcHorizScale(rstddev,variableType,waveBandIndex_opt=1) ! IN
655 call horizCorrelFunction(rstddev,variableType,waveBandIndex_opt=1) ! IN
656
657 deallocate(rstddev)
658 deallocate(corns)
659
660 case ('HVCORREL_LOCAL')
661 write(*,*)
662 write(*,*) 'Computing Local Correlation'
663
664 call ens_removeGlobalMean(ensPerts)
665 call spectralFilter2(ensPerts, & ! IN
666 ensPerts_ptr, & ! OUT
667 waveBandIndex_opt=1) ! IN
668 call ens_computeStdDev(ensPerts_ptr)
669 call ens_normalize(ensPerts_ptr)
670 call calcLocalCorrelations(ensPerts_ptr) ! IN
671
672 case ('VCORRMATRIX')
673 write(*,*)
674 write(*,*) 'Computing Local Correlation'
675
676 if (doSpectralFilter) then
677 call ens_removeGlobalMean(ensPerts)
678 call spectralFilter2(ensPerts, & ! IN
679 ensPerts_ptr, & ! OUT
680 waveBandIndex_opt=1) ! IN
681 end if
682 call ens_computeStdDev(ensPerts_ptr)
683 call ens_normalize(ensPerts_ptr)
684 call calcLocalVertCorrMatrix(ensPerts_ptr) ! IN
685
686 case ('LOCALIZATIONRADII')
687 write(*,*)
688 write(*,*) 'Estimating the optimal covariance localization radii'
689
690 if (nWaveBand > 1 .and. nVertWaveBand > 1) then
691 call utl_abort('csg_toolbox: cannot do horizontal AND vertical scale-decomposition')
692 end if
693
694 call ens_removeGlobalMean(ensPerts)
695
696 do waveBandIndex = 1, max(nWaveBand,nVertWaveBand)
697
698 if (nVertWaveBand > 1) then
699 if (waveBandIndex == 1) then
700 call vms_computeModesFromFunction(vco_ens, lengthScaleTop, lengthScaleBot, & ! IN
701 vModes) ! OUT
702 end if
703 call vertModesFilter(ensPerts, & ! IN
704 ensPerts_ptr, & ! OUT
705 vModes, waveBandIndex) ! IN
706 else
707 call spectralFilter2(ensPerts, & ! IN
708 ensPerts_ptr, & ! OUT
709 waveBandIndex_opt=waveBandIndex) ! IN
710 end if
711
712 call ens_computeStdDev(ensPerts_ptr)
713
714 if (waveBandIndex == 1) then
715 call ens_copyEnsStdDev(ensPerts_ptr, statevector_template) ! IN
716 call bmd_setup(statevector_template, hco_ens, nEns, pressureProfile_M, & ! IN
717 pressureProfile_T, max(nWaveBand,nVertWaveBand)) ! IN
718 end if
719
720 call bmd_localizationRadii(ensPerts_ptr, waveBandIndex_opt=waveBandIndex) ! IN
721
722 end do
723
724 case ('STDDEV')
725 write(*,*)
726 write(*,*) 'Computing standard deviations using PSI/CHI for winds'
727
728 ! u,v to psi,chi
729 call gvt_setup(hco_ens, hco_ens, vco_ens)
730 call gvt_transform(ensPerts_ptr, 'UVtoPsiChi')
731 if (.not. ens_varExist(ensPerts_ptr,'PP') .and. &
732 .not. ens_varExist(ensPerts_ptr,'CC') ) then
733 call ens_modifyVarName(ensPerts_ptr, 'UU', 'PP')
734 call ens_modifyVarName(ensPerts_ptr, 'VV', 'CC')
735 end if
736
737 ! Compute the grid point std dev
738 call ens_computeStdDev(ensPerts_ptr)
739 call ens_copyEnsStdDev(ensPerts_ptr, statevector_template) ! IN
740 call gio_writeToFile(statevector_template, './stddev.fst', 'STDDEV_GRIDP', &
741 typvar_opt = 'E', numBits_opt = 32)
742
743 ! Compute the zonal std dev
744 call gbi_setup(gbi_zonalMean, 'YrowBand', statevector_template, hco_ens)
745 call gbi_stdDev(gbi_zonalMean, ensPerts_ptr, & ! IN
746 statevector_template) ! OUT
747 call gio_writeToFile(statevector_template, './stddev.fst', 'STDDEV_ZONAL', &
748 typvar_opt = 'E', numBits_opt = 32)
749
750 case ('POWERSPEC')
751 write(*,*)
752 write(*,*) 'Computing power spectra'
753
754 call calcPowerSpec(ensPerts_ptr, & ! IN
755 powerSpec) ! OUT
756 call writePowerSpec(powerSpec, modelSpace) ! IN
757
758 case ('VERTMODES_SPEC')
759 write(*,*)
760 write(*,*) 'Computing vertical modes spectra'
761
762 if (nWaveBand > 1 .or. nVertWaveBand > 1) then
763 call utl_abort('csg_toolbox: waveband decomposition cannot be use when TOOLBOX=VERTMODES_SPEC')
764 end if
765
766 call vms_computeModesFromFunction(vco_ens, lengthScaleTop, lengthScaleBot, & ! IN
767 vModes) ! OUT
768 call vms_writeModes(vModes)
769
770 call calcVertModesSpec(ensPerts_ptr,vModes) ! IN
771
772 case ('VERTMODES_WAVEBAND')
773 write(*,*)
774 write(*,*) 'Computing vertical-scale decomposed ensemble perturbations'
775
776 if (nVertWaveBand == 1) then
777 call utl_abort('csg_toolbox: please specify a vertical waveband decomposition when TOOLBOX=VERTMODES_WAVEBAND')
778 end if
779
780 call vms_computeModesFromFunction(vco_ens, lengthScaleTop, lengthScaleBot, & ! IN
781 vModes) ! OUT
782 call vms_writeModes(vModes)
783
784 do waveBandIndex = 1, nVertWaveBand
785
786 call vertModesFilter(ensPerts, & ! IN
787 ensPerts_ptr, & ! OUT
788 vModes, waveBandIndex) ! IN
789
790 write(wbnum,'(I2.2)') waveBandIndex
791
792 ! Compute the grid point std dev
793 call ens_computeStdDev(ensPerts_ptr)
794 call ens_copyEnsStdDev(ensPerts_ptr, & ! IN
795 statevector_template) ! OUT
796 call gio_writeToFile(statevector_template, './stddev_'//trim(wbnum)//'.fst', 'STD_GRIDP_'//trim(wbnum), &
797 typvar_opt = 'E', numBits_opt = 32)
798
799 ! Compute the global std dev
800 if (waveBandIndex == 1) then
801 call gbi_setup(gbi_globalMean, 'HorizontalMean', statevector_template, hco_ens)
802 end if
803 call gbi_stdDev(gbi_globalMean, ensPerts_ptr, & ! IN
804 statevector_template) ! OUT
805 call gio_writeToFile(statevector_template, './stddev_'//trim(wbnum)//'.fst', 'STD_GLB_'//trim(wbnum), &
806 typvar_opt = 'E', numBits_opt = 32)
807
808 ! Write the scale-decomposed ensemble perturbations for member=1
809 call ens_copyMember(ensPerts_ptr, stateVector_template, 1)
810 call gio_writeToFile(statevector_template, './ensPert_0001_'//trim(wbnum)//'.fst', 'ENSPERT_'//trim(wbnum), &
811 numBits_opt = 32)
812
813 end do
814
815 ! Write the full ensemble perturbations for member=1
816 call ens_copyMember(ensPerts, stateVector_template, 1)
817 call gio_writeToFile(statevector_template, './ensPert_0001.fst', 'ENSPERT', &
818 numBits_opt = 32)
819
820 ! Compute the grid point std dev for the full perturbations
821 call ens_computeStdDev(ensPerts)
822 call ens_copyEnsStdDev(ensPerts, & ! IN
823 statevector_template) ! OUT
824 call gio_writeToFile(statevector_template, './stddev.fst', 'STD_GRIDP', &
825 typvar_opt = 'E', numBits_opt = 32)
826
827 ! Compute the global std dev for the full perturbations
828 call gbi_stdDev(gbi_globalMean, ensPerts, & ! IN
829 statevector_template) ! OUT
830 call gio_writeToFile(statevector_template, './stddev.fst', 'STD_GLB', &
831 typvar_opt = 'E', numBits_opt = 32)
832
833 case default
834 write(*,*)
835 write(*,*) 'Unknown TOOL in csg_toolbox : ', trim(tool)
836 call utl_abort('csg_toolbox')
837 end select
838
839 !
840 !- Write the estimated pressure profiles
841 !
842 if (vco_ens%vgridPresent) then
843 call writePressureProfiles
844 end if
845
846 !
847 !- Ending
848 !
849 call ens_deallocate(ensPerts)
850
851 write(*,*)
852 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
853 write(*,*) 'csl_toolbox: Done!'
854
855 end subroutine csg_toolbox
856
857 !--------------------------------------------------------------------------
858 ! WRITESPSTATS
859 !--------------------------------------------------------------------------
860 subroutine writeSpStats(ptot,theta)
861 !
862 !:Purpose: Write the spectral representation of PtoT and THETA to files
863 !
864 implicit none
865
866 ! Arguments:
867 real(8), intent(in) :: PtoT(:,:,:)
868 real(8), intent(in) :: theta(:,:)
869
870 ! Locals:
871 integer jn,ierr,ipak,latIndex,levIndex1,levIndex2,nlev
872 integer fstouv,fnom,fstfrm,fclos
873 integer ip1,ip3,kni,knj,idatyp,idateo
874 integer :: nulstats
875 real(8) :: bufz(nLevEns_M),bufyz(nj,nLevEns_M),zsp(0:ntrunc,nLevEns_M)
876 real(8) :: bufptot(nj,(nLevEns_T+1)*nLevEns_M),spptot(0:ntrunc,(nLevEns_T+1)*nLevEns_M)
877 real(8) :: zspptot(nLevEns_T+1,nLevEns_M)
878
879 if (mmpi_myid /= 0) return
880
881 nulstats=0
882 ierr = fnom (nulstats,'./stats_sp.fst','RND',0)
883 ierr = fstouv(nulstats,'RND')
884
885 ipak = -32
886 idatyp = 5
887 ip1 = 0
888 ip3 = nens
889 idateo = 0
890
891 ! write out SP_THETA
892
893 do latIndex = 1, nj
894 do levIndex1 = 1, nLevEns_M
895 bufyz(latIndex,levIndex1) = theta(levIndex1,latIndex)
896 end do
897 end do
898
899 call gst_zlegdir(gstID_nkgdimEns,bufyz,zsp,nLevEns_M)
900
901 do jn = 0, ntrunc
902 do levIndex1=1, nLevEns_M
903 bufz(levIndex1) = zsp(jn,levIndex1)
904 end do
905
906 ierr = utl_fstecr(bufz,ipak,nulstats,idateo,0,0,nlevEns_M,1,1, &
907 ip1,jn,ip3,'X','ZZ','SP_THETA','X',0,0,0,0,idatyp,.true.)
908
909 end do
910
911 ! write out SP_PTOT
912
913 do latIndex = 1, nj
914 do levIndex1 = 1, (nLevEns_T+1)
915 do levIndex2 = 1, nLevEns_M
916 bufptot(latIndex,(levIndex2-1)*(nLevEns_T+1)+levIndex1) = PtoT(levIndex1,levIndex2,latIndex)
917 end do
918 end do
919 end do
920
921 nlev=(nLevEns_T+1)*nLevEns_M
922 call gst_zlegdir(gstID_nkgdimEns,bufptot,spptot,nLev)
923
924 do jn = 0, ntrunc
925 do levIndex1 = 1, (nLevEns_T+1)
926 do levIndex2 = 1, nLevEns_M
927 zspptot(levIndex1,levIndex2) = spptot(jn,(levIndex2-1)*(nLevEns_T+1)+levIndex1)
928 end do
929 end do
930
931 kni=nLevEns_T+1
932 knj=nLevEns_M
933 ierr = utl_fstecr(zspptot,ipak,nulstats,idateo,0,0,kni,knj,1, &
934 ip1,jn,ip3,'X','ZZ','SP_PTOT ','X',0,0,0,0,idatyp,.true.)
935 end do
936
937 ierr = fstfrm(nulstats)
938 ierr = fclos (nulstats)
939
940 write(*,*) 'finished writing statistics...'
941
942 end subroutine writeSpStats
943
944 !--------------------------------------------------------------------------
945 ! REMOVEBALANCEDCHI
946 !--------------------------------------------------------------------------
947 subroutine removeBalancedChi(ensPerturbations,theta)
948 !
949 !:Purpose: Subtract the balanced components of velocity potential
950 ! from the full variable
951 !
952 implicit none
953
954 ! Arguments:
955 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
956 real(8), intent(in) :: theta(:,:)
957
958 ! Locals:
959 real(4), pointer :: psi_ptr(:,:,:), chi_ptr(:,:,:)
960 integer :: ensIndex,latIndex,levIndex,lonIndex
961
962 do ensIndex = 1,nens
963 psi_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,1:nlevEns_M,ensIndex)
964 chi_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,(nlevEns_M+1):(2*nlevEns_M),ensIndex)
965
966 do levIndex = 1, nLevEns_M
967 do latIndex = myLatBeg, myLatEnd
968 do lonIndex = myLonBeg, myLonEnd
969 chi_ptr(lonIndex,latIndex,levIndex) = chi_ptr(lonIndex,latIndex,levIndex) + &
970 tan(theta(levIndex,latIndex))*psi_ptr(lonIndex,latIndex,levIndex)
971 end do
972 end do
973 end do
974
975 end do
976
977 write(*,*) 'finished removing balanced chi...'
978
979 end subroutine removeBalancedChi
980
981 !--------------------------------------------------------------------------
982 ! REMOVEBALANCEDT_PS
983 !--------------------------------------------------------------------------
984 subroutine removeBalancedT_Ps(ensPerturbations,ensBalPerturbations,PtoT)
985 !
986 !:Purpose: Subtract the balanced components of temperature and surface
987 ! pressure from the full variables
988 !
989 implicit none
990
991 ! Arguments:
992 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
993 real(4), pointer, intent(inout) :: ensBalPerturbations(:,:,:,:)
994 real(8), intent(in) :: PtoT(:,:,:)
995
996 ! Locals:
997 real(4),pointer :: tt_ptr(:,:,:), ps_ptr(:,:,:), ttb_ptr(:,:,:), psb_ptr(:,:,:)
998 real(8) :: spectralState(nla_mpilocal,2,nLevEns_M), spBalancedP(nla_mpilocal,2,nlevEns_M)
999 real(8) :: balancedP(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlevEns_M)
1000 real(8) :: psi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLevEns_M)
1001 integer :: ensIndex, latIndex, lonIndex, jk1, jk2
1002
1003 do ensIndex=1,nens
1004
1005 write(*,*) 'removing balanced T and Ps for member ', ensIndex
1006
1007 psi(:,:,:) = ensPerturbations(:,:,1:nlevEns_M,ensIndex)
1008 call gst_setID(gstID_nLevEns_M)
1009 call gst_reespe(spectralState,psi)
1010 call calcBalancedP(spectralState,spBalancedP)
1011 call gst_speree(spBalancedP,balancedP)
1012
1013 tt_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,(1+2*nLevEns_M):(2*nLevEns_M+1*nLevEns_T),ensIndex)
1014 ps_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,(1+2*nLevEns_M+2*nLevEns_T):(1+2*nLevEns_M+2*nLevEns_T),ensIndex)
1015 ttb_ptr(myLonBeg:,myLatBeg:,1:) => ensBalPerturbations(:,:,1:nLevEns_T,ensIndex)
1016 psb_ptr(myLonBeg:,myLatBeg:,1:) => ensBalPerturbations(:,:,(1+nLevEns_T):(1+nLevEns_T),ensIndex)
1017
1018 ttb_ptr(:,:,:)=0.0d0
1019 psb_ptr(:,:,:)=0.0d0
1020
1021 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,jk1,jk2)
1022 do latIndex = myLatBeg, myLatEnd
1023 do lonIndex = myLonBeg, myLonEnd
1024 do jk1 = 1, nLevEns_T
1025 do jk2 = 1, nlevptot
1026 ttb_ptr(lonIndex,latIndex,jk1) = ttb_ptr(lonIndex,latIndex,jk1) + PtoT(jk1,jk2,latIndex)*balancedP(lonIndex,latIndex,jk2)
1027 end do
1028 tt_ptr(lonIndex,latIndex,jk1) = tt_ptr(lonIndex,latIndex,jk1) - ttb_ptr(lonIndex,latIndex,jk1)
1029 end do
1030 do jk2 = 1, nlevptot
1031 psb_ptr(lonIndex,latIndex,1) = psb_ptr(lonIndex,latIndex,1) + PtoT(nLevEns_T+1,jk2,latIndex)*balancedP(lonIndex,latIndex,jk2)
1032 end do
1033 ps_ptr(lonIndex,latIndex,1) = ps_ptr(lonIndex,latIndex,1) - psb_ptr(lonIndex,latIndex,1)
1034 end do
1035 end do
1036 !$OMP END PARALLEL DO
1037
1038 end do
1039
1040 write(*,*) 'finished removing balanced T and Ps...'
1041
1042 end subroutine removeBalancedT_Ps
1043
1044 !--------------------------------------------------------------------------
1045 ! CALCCORRELATIONS
1046 !--------------------------------------------------------------------------
1047 subroutine calcCorrelations(ensPerturbations,corns,rstddev,latMask_opt)
1048 !
1049 !:Purpose: Calculate the homogeneous and isotropic correlations in spectral space
1050 !
1051 implicit none
1052
1053 ! Arguments:
1054 real(4), pointer, intent(in) :: ensPerturbations(:,:,:,:)
1055 real(8), intent(out) :: corns(nkgdimEns,nkgdimEns,0:ntrunc)
1056 real(8), intent(out) :: rstddev(nkgdimEns,0:ntrunc)
1057 real(8), optional, intent(in) :: latMask_opt(:)
1058
1059 ! Locals:
1060 real(8) :: spectralState(nla_mpilocal,2,nkgdimEns)
1061 real(8) :: corns_mpiglobal(nkgdimEns,nkgdimEns,0:ntrunc)
1062 real(8) :: gridState(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdimEns)
1063 real(8) :: dfact,dfact2,dsummed
1064 integer :: ensIndex,ila_mpilocal,ila_mpiglobal,jn,jm,jk1,jk2,latIndex,nsize,ierr
1065
1066 call utl_tmg_start(120,'--CalcStats_Corr')
1067
1068 corns(:,:,:) = 0.0d0
1069 do ensIndex = 1, nens
1070
1071 write(*,*) 'calcCorrelations: processing member ',ensIndex
1072
1073 gridState(:,:,:) = ensPerturbations(:,:,:,ensIndex)
1074 if(present(latMask_opt)) then
1075 do latIndex = myLatBeg, myLatEnd
1076 gridState(:,latIndex,:) = latMask_opt(latIndex)*gridState(:,latIndex,:)
1077 end do
1078 end if
1079 call gst_setID(gstID_nkgdimEns)
1080 call gst_reespe(spectralState,gridState)
1081
1082 !$OMP PARALLEL DO PRIVATE (jn,jm,dfact,ila_mpilocal,ila_mpiglobal,jk1,jk2)
1083 do jn = mynBeg, mynEnd, mynSkip
1084 do jm = mymBeg, mymEnd, mymSkip
1085 if(jm.le.jn) then
1086 dfact = 2.0d0
1087 if (jm.eq.0) dfact = 1.0d0
1088 ila_mpiglobal = gst_getNind(jm,gstID_nkgdimEns) + jn - jm
1089 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1090 do jk1 = 1, nkgdimEns
1091 do jk2 = 1, nkgdimEns
1092 corns(jk1,jk2,jn) = corns(jk1,jk2,jn) + &
1093 dfact*( spectralState(ila_mpilocal,1,jk1)*spectralState(ila_mpilocal,1,jk2) + &
1094 spectralState(ila_mpilocal,2,jk1)*spectralState(ila_mpilocal,2,jk2) )
1095 end do
1096 end do
1097 end if
1098 end do
1099 end do
1100 !$OMP END PARALLEL DO
1101
1102 end do
1103
1104 ! communicate between all tasks
1105 nsize = nkgdimEns*nkgdimEns*(1+ntrunc)
1106 call rpn_comm_allreduce(corns,corns_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
1107 corns(:,:,:) = corns_mpiglobal(:,:,:)
1108
1109 !$OMP PARALLEL DO PRIVATE (jn,jk1)
1110 do jn = 0, ntrunc
1111 do jk1 = 1, nkgdimEns
1112 if(abs(corns(jk1,jk1,jn)).gt.0.0d0) then
1113 rstddev(jk1,jn) = dsqrt(abs(corns(jk1,jk1,jn)))
1114 else
1115 rstddev(jk1,jn) = 0.0d0
1116 end if
1117 end do
1118 end do
1119 !$OMP END PARALLEL DO
1120
1121 !$OMP PARALLEL DO PRIVATE (jn,jk1,jk2)
1122 do jn = 0, ntrunc
1123 do jk1 = 1, nkgdimEns
1124 do jk2 = 1, nkgdimEns
1125 if(rstddev(jk1,jn).ne.0..and.rstddev(jk2,jn).ne.0.) then
1126 corns(jk1,jk2,jn) = corns(jk1,jk2,jn)/(rstddev(jk1,jn)*rstddev(jk2,jn))
1127 else
1128 corns(jk1,jk2,jn) = 0.0d0
1129 end if
1130 end do
1131 end do
1132 end do
1133 !$OMP END PARALLEL DO
1134
1135 dfact2 = 1.0d0/sqrt(dble(nens-1))
1136 do jn = 0, ntrunc
1137 dfact = 1.0d0/sqrt(2.0d0*dble(jn) + 1.0d0)
1138 do jk1 = 1, nkgdimEns
1139 rstddev(jk1,jn) = rstddev(jk1,jn)*dfact2*dfact
1140 end do
1141 end do
1142
1143 ! Normalize to ensure correlations in horizontal and Multiply by sqrt(0.5) to make valid for m.ne.0
1144 !$OMP PARALLEL DO PRIVATE (jk1,jn,dsummed)
1145 do jk1 = 1, nkgdimEns
1146 dsummed=0.0d0
1147 do jn = 0, ntrunc
1148 dsummed=dsummed + (rstddev(jk1,jn)**2)*((2.0d0*dble(jn))+1.0d0)/2.0d0
1149 end do
1150 do jn = 0, ntrunc
1151 if(dsummed.gt.0.0d0) rstddev(jk1,jn)=rstddev(jk1,jn)*sqrt(0.5d0/dsummed)
1152 end do
1153 end do
1154 !$OMP END PARALLEL DO
1155
1156 call utl_tmg_stop(120)
1157 write(*,*) 'finished computing correlations...'
1158
1159 end subroutine calcCorrelations
1160
1161 !--------------------------------------------------------------------------
1162 ! CALCCORRELATIONS2
1163 !--------------------------------------------------------------------------
1164 subroutine calcCorrelations2(ensPerts,corns,rstddev,latMask_opt)
1165 !
1166 !:Purpose: Calculate the homogeneous and isotropic correlations in spectral space
1167 !
1168 implicit none
1169
1170 ! Arguments:
1171 type(struct_ens), intent(inout) :: ensPerts
1172 real(8), intent(out) :: corns(nkgdimEns,nkgdimEns,0:ntrunc)
1173 real(8), intent(out) :: rstddev(nkgdimEns,0:ntrunc)
1174 real(8), optional, intent(in) :: latMask_opt(:)
1175
1176 ! Locals:
1177 real(4), pointer :: ptr4d_r4(:,:,:,:)
1178 real(8) :: corns_mpiglobal(nkgdimEns,nkgdimEns,0:ntrunc)
1179 real(8) :: spectralState(nla_mpilocal,2,nkgdimEns)
1180 real(8) :: gridState(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdimEns)
1181 real(8) :: dfact, dfact2, dsummed
1182 integer :: ensIndex, ila_mpilocal, ila_mpiglobal, jn, jm, jk1, jk2
1183 integer :: levIndex, latIndex, nsize, ierr
1184
1185 call utl_tmg_start(120,'--CalcStats_Corr')
1186
1187 corns(:,:,:) = 0.0d0
1188 do ensIndex = 1, nens
1189
1190 write(*,*) 'calcCorrelations: processing member ',ensIndex
1191
1192 !- 2.1 Extract fields from ensPerturbations
1193 do levIndex = 1, ens_getNumK(ensPerts)
1194 ptr4d_r4 => ens_getOneLev_r4(ensPerts,levIndex)
1195 gridState(:,:,levIndex) = real(ptr4d_r4(ensIndex,1,:,:),8)
1196 end do
1197 if ( present(latMask_opt) ) then
1198 do latIndex = myLatBeg, myLatEnd
1199 gridState(:,latIndex,:) = latMask_opt(latIndex)*gridState(:,latIndex,:)
1200 end do
1201 end if
1202 call gst_setID(gstID_nkgdimEns)
1203 call gst_reespe(spectralState,gridState)
1204
1205 !$OMP PARALLEL DO PRIVATE (jn,jm,dfact,ila_mpilocal,ila_mpiglobal,jk1,jk2)
1206 do jn = mynBeg, mynEnd, mynSkip
1207 do jm = mymBeg, mymEnd, mymSkip
1208 if(jm.le.jn) then
1209 dfact = 2.0d0
1210 if (jm.eq.0) dfact = 1.0d0
1211 ila_mpiglobal = gst_getNind(jm,gstID_nkgdimEns) + jn - jm
1212 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1213 do jk1 = 1, nkgdimEns
1214 do jk2 = 1, nkgdimEns
1215 corns(jk1,jk2,jn) = corns(jk1,jk2,jn) + &
1216 dfact*( spectralState(ila_mpilocal,1,jk1)*spectralState(ila_mpilocal,1,jk2) + &
1217 spectralState(ila_mpilocal,2,jk1)*spectralState(ila_mpilocal,2,jk2) )
1218 end do
1219 end do
1220 end if
1221 end do
1222 end do
1223 !$OMP END PARALLEL DO
1224
1225 end do
1226
1227 ! communicate between all tasks
1228 nsize = nkgdimEns*nkgdimEns*(1+ntrunc)
1229 call rpn_comm_allreduce(corns,corns_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
1230 corns(:,:,:) = corns_mpiglobal(:,:,:)
1231
1232 !$OMP PARALLEL DO PRIVATE (jn,jk1)
1233 do jn = 0, ntrunc
1234 do jk1 = 1, nkgdimEns
1235 if(abs(corns(jk1,jk1,jn)).gt.0.0d0) then
1236 rstddev(jk1,jn) = dsqrt(abs(corns(jk1,jk1,jn)))
1237 else
1238 rstddev(jk1,jn) = 0.0d0
1239 end if
1240 end do
1241 end do
1242 !$OMP END PARALLEL DO
1243
1244 !$OMP PARALLEL DO PRIVATE (jn,jk1,jk2)
1245 do jn = 0, ntrunc
1246 do jk1 = 1, nkgdimEns
1247 do jk2 = 1, nkgdimEns
1248 if(rstddev(jk1,jn).ne.0..and.rstddev(jk2,jn).ne.0.) then
1249 corns(jk1,jk2,jn) = corns(jk1,jk2,jn)/(rstddev(jk1,jn)*rstddev(jk2,jn))
1250 else
1251 corns(jk1,jk2,jn) = 0.0d0
1252 end if
1253 end do
1254 end do
1255 end do
1256 !$OMP END PARALLEL DO
1257
1258 dfact2 = 1.0d0/sqrt(dble(nens-1))
1259 do jn = 0, ntrunc
1260 dfact = 1.0d0/sqrt(2.0d0*dble(jn) + 1.0d0)
1261 do jk1 = 1, nkgdimEns
1262 rstddev(jk1,jn) = rstddev(jk1,jn)*dfact2*dfact
1263 end do
1264 end do
1265
1266 ! Normalize to ensure correlations in horizontal and Multiply by sqrt(0.5) to make valid for m.ne.0
1267 !$OMP PARALLEL DO PRIVATE (jk1,jn,dsummed)
1268 do jk1 = 1, nkgdimEns
1269 dsummed=0.0d0
1270 do jn = 0, ntrunc
1271 dsummed=dsummed + (rstddev(jk1,jn)**2)*((2.0d0*dble(jn))+1.0d0)/2.0d0
1272 end do
1273 do jn = 0, ntrunc
1274 if(dsummed.gt.0.0d0) rstddev(jk1,jn)=rstddev(jk1,jn)*sqrt(0.5d0/dsummed)
1275 end do
1276 end do
1277 !$OMP END PARALLEL DO
1278
1279 call utl_tmg_stop(120)
1280 write(*,*) 'finished computing correlations...'
1281
1282 end subroutine calcCorrelations2
1283
1284 !--------------------------------------------------------------------------
1285 ! CALCPOWERSPEC
1286 !--------------------------------------------------------------------------
1287 subroutine calcPowerSpec(ensPerts,powerSpec)
1288 !
1289 !:Purpose: Calculate the horizontal power spectrum of the ensemble
1290 ! of error samples
1291 !
1292 implicit none
1293
1294 ! Arguments:
1295 type(struct_ens), intent(inout) :: ensPerts
1296 real(8), allocatable, intent(out) :: powerSpec(:,:)
1297
1298 ! Locals:
1299 real(8), allocatable :: ensPertSP(:,:,:)
1300 real(8), allocatable :: ensPertGD(:,:,:)
1301 real(4), pointer :: ptr4d_r4(:,:,:,:)
1302 real(8) :: dfact, dfact2
1303 integer :: gstPowerSpecID
1304 integer :: memberIndex, levIndex, latIndex, lonIndex
1305 integer :: jn, jm, ila_mpilocal, ila_mpiglobal
1306
1307 allocate(powerSpec (ens_getNumK(ensPerts),0:ntrunc))
1308
1309 !
1310 !- Spectral decomposition and spectral coefficient summation
1311 !
1312 gstPowerSpecID = gst_setup(ni,nj,nTrunc,nEnsOverDimension)
1313
1314 allocate(ensPertSP(nla_mpilocal,2,nEnsOverDimension))
1315 allocate(ensPertGD(nEnsOverDimension,myLonBeg:myLonEnd,myLatBeg:myLatEnd))
1316
1317 powerSpec(:,:)=0.0d0
1318
1319 do levIndex = 1, ens_getNumK(ensPerts) ! Loop on variables and vertical levels
1320
1321 ptr4d_r4 => ens_getOneLev_r4(ensPerts,levIndex)
1322
1323 do latIndex = myLatBeg, myLatEnd
1324 ensPertGD(:,:,latIndex) = 0.0d0
1325 end do
1326
1327 !$OMP PARALLEL DO PRIVATE (memberIndex,latIndex,lonIndex)
1328 do latIndex = myLatBeg, myLatEnd
1329 do lonIndex = myLonBeg, myLonEnd
1330 do memberIndex = 1, nEns
1331 ensPertGD(memberIndex,lonIndex,latIndex) = dble(ptr4d_r4(memberIndex,1,lonIndex,latIndex))
1332 end do
1333 end do
1334 end do
1335 !$OMP END PARALLEL DO
1336
1337 !- GridPoint space -> Spectral Space
1338 call gst_setID(gstPowerSpecID) ! IN
1339 call gst_reespe_kij(ensPertSP, & ! OUT
1340 ensPertGD) ! IN
1341
1342 !$OMP PARALLEL DO PRIVATE (memberIndex,jn,jm,dfact,ila_mpilocal,ila_mpiglobal)
1343 do jn = mynBeg, mynEnd, mynSkip
1344 do jm = mymBeg, mymEnd, mymSkip
1345 if (jm .le. jn) then
1346 dfact = 2.0d0
1347 if (jm.eq.0) dfact = 1.0d0
1348 ila_mpiglobal = gst_getNind(jm,gstPowerSpecID) + jn - jm
1349 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1350 do memberIndex = 1, nEns
1351 powerSpec(levIndex,jn) = powerSpec(levIndex,jn) + &
1352 dfact*( ensPertSP(ila_mpilocal,1,memberIndex)**2 + ensPertSP(ila_mpilocal,2,memberIndex)**2 )
1353 end do
1354 end if
1355 end do
1356 end do
1357 !$OMP END PARALLEL DO
1358
1359 end do
1360
1361 !
1362 !- Communicate between all tasks
1363 !
1364 call mmpi_allreduce_sumR8_2d(powerSpec, "GRID")
1365
1366 !
1367 !- Apply the appropriate scaling
1368 !
1369 dfact2 = 1.0d0/sqrt(dble(nens-1))
1370 do jn = 0, ntrunc
1371 dfact = 1.0d0/sqrt(2.0d0*dble(jn) + 1.0d0)
1372 do levIndex = 1, ens_getNumK(ensPerts)
1373 powerSpec(levIndex,jn) = powerSpec(levIndex,jn)*dfact2*dfact
1374 end do
1375 end do
1376
1377 write(*,*) 'finished computing power spectrum...'
1378
1379 end subroutine calcPowerSpec
1380
1381 !--------------------------------------------------------------------------
1382 ! WRITESTATS
1383 !--------------------------------------------------------------------------
1384 subroutine writeStats(corns, rstddev, ptot_opt, theta_opt, waveBandIndex_opt, latBand_opt)
1385 !
1386 !:Purpose: Write several components of the BHI matrix to a file
1387 !
1388 implicit none
1389
1390 ! Arguments:
1391 real(8), intent(in) :: corns(nkgdimEns,nkgdimEns,0:ntrunc)
1392 real(8), intent(in) :: rstddev(nkgdimEns,0:ntrunc)
1393 real(8), optional, intent(in) :: PtoT_opt(:,:,:)
1394 real(8), optional, intent(in) :: theta_opt(:,:)
1395 integer, optional, intent(in) :: waveBandIndex_opt
1396 integer, optional, intent(in) :: latBand_opt
1397
1398 ! Locals:
1399 real(8) :: prcor(nkgdimEns,nkgdimEns)
1400 integer :: jn,ierr,ipak,jk,jl
1401 integer :: fstouv,fnom,fstfrm,fclos
1402 integer :: ip1,ip2,ip3,idatyp,idateo
1403 integer :: nulstats
1404 character(len=128) :: outfilename
1405 character(len=2) :: wbnum
1406
1407 if (mmpi_myid /= 0) return
1408
1409 nulstats=0
1410 if ( nWaveBand == 1 ) then
1411 outfilename='./bgcov.fst'
1412 else
1413 if (.not. present(waveBandIndex_opt)) then
1414 write(*,*) 'writeStats: No waveBandIndex was supplied!!!'
1415 call utl_abort('calbmatrix_glb')
1416 end if
1417 write(wbnum,'(I2.2)') waveBandIndex_opt
1418 outfilename='./bgcov_'//trim(wbnum)//'.fst'
1419 end if
1420 ierr = fnom (nulstats,trim(outfilename),'RND',0)
1421 ierr = fstouv(nulstats,'RND')
1422
1423 ipak = -32
1424 idatyp = 5
1425 ip1 = 0
1426 ip2 = 0
1427 ip3 = nens
1428 idateo = 0
1429
1430 if(present(latBand_opt)) ip1 = latBand_opt
1431
1432 if (present(ptot_opt)) then
1433 ierr = utl_fstecr(ptot_opt,ipak,nulstats,idateo,0,0,nlevEns_T+1,nlevEns_M,nj, &
1434 ip1,ip2,ip3,'X','ZZ','P_to_T ','X',0,0,0,0,idatyp,.true.)
1435 end if
1436 if (present(theta_opt)) then
1437 ierr = utl_fstecr(theta_opt,ipak,nulstats,idateo,0,0,nlevEns_M,nj,1, &
1438 ip1,ip2,ip3,'X','ZZ','THETA ','X',0,0,0,0,idatyp,.true.)
1439 end if
1440
1441 do jn = 0, ntrunc
1442 ip2 = jn
1443 ierr = utl_fstecr(corns(:,:,jn),ipak,nulstats,idateo,0,0,nkgdimEns,nkgdimEns,1, &
1444 ip1,ip2,ip3,'X','ZZ','CORRNS ','X',0,0,0,0,idatyp,.true.)
1445 end do
1446
1447 do jn = 0, ntrunc
1448 ip2 = jn
1449 ierr = utl_fstecr(rstddev(:,jn),ipak,nulstats,idateo,0,0,nkgdimEns,1,1, &
1450 ip1,ip2,ip3,'X','SS','RSTDDEV ','X',0,0,0,0,idatyp,.true.)
1451 end do
1452
1453
1454 ! Computing the total vertical correlation matrix
1455 do jk = 1, nkgdimEns
1456 do jl = 1, nkgdimEns
1457 prcor(jk,jl) = 0.0d0
1458 do jn = 0, ntrunc
1459 prcor(jk,jl) = prcor(jk,jl) + ((2*jn+1)*rstddev(jk,jn)*rstddev(jl,jn)*corns(jk,jl,jn))
1460 end do
1461 end do
1462 end do
1463
1464 do jk = 1, nkgdimEns
1465 do jl = 1, nkgdimEns
1466 if(prcor(jk,jk)*prcor(jl,jl) .gt. 0.0d0) then
1467 prcor(jk,jl) = prcor(jk,jl) / (sqrt(prcor(jk,jk)*prcor(jl,jl)))
1468 else
1469 prcor(jk,jl) = 0.0d0
1470 end if
1471 end do
1472 end do
1473 ip2 =0
1474 ierr = utl_fstecr(prcor(:,:),ipak,nulstats,idateo,0,0,nkgdimEns,nkgdimEns,1, &
1475 ip1,ip2,ip3,'X','ZV','CORVERT ','X',0,0,0,0,idatyp,.true.)
1476
1477 ierr = fstfrm(nulstats)
1478 ierr = fclos (nulstats)
1479
1480 write(*,*) 'finished writing statistics...'
1481
1482 end subroutine writeStats
1483
1484 !--------------------------------------------------------------------------
1485 ! WRITEPRESSUREPROFILES
1486 !--------------------------------------------------------------------------
1487 subroutine writePressureProfiles
1488 !
1489 !:Purpose: Write the profiles of pressure to ascii files
1490 !
1491 implicit none
1492
1493 ! Locals:
1494 character(len=128) :: outfilename
1495 integer :: jk
1496
1497 if (mmpi_myid /= 0) return
1498
1499 outfilename = "./pressureProfile_M.txt"
1500 open (unit=99,file=outfilename,action="write",status="new")
1501 do jk = 1, nLevEns_M
1502 write(99,'(I3,2X,F7.2)') jk, pressureProfile_M(jk)/100.d0
1503 end do
1504 close(unit=99)
1505
1506 outfilename = "./pressureProfile_T.txt"
1507 open (unit=99,file=outfilename,action="write",status="new")
1508 do jk = 1, nLevEns_T
1509 write(99,'(I3,2X,F7.2)') jk, pressureProfile_T(jk)/100.d0
1510 end do
1511 close(unit=99)
1512
1513 write(*,*) 'finished writing pressure profiles...'
1514
1515 end subroutine writePressureProfiles
1516
1517 !--------------------------------------------------------------------------
1518 ! WRITESTDDEV
1519 !--------------------------------------------------------------------------
1520 subroutine writeStddev(stddevZonAvg,stddev3d,stddevZonAvgUnbal_opt,stddev3dUnbal_opt)
1521 !
1522 !:Purpose: Write the stddev to a file
1523 !
1524 implicit none
1525
1526 ! Arguments:
1527 real(8), pointer, intent(in) :: stddevZonAvg(:,:)
1528 real(8), pointer, intent(in) :: stddev3d(:,:,:)
1529 real(8), pointer, optional, intent(in) :: stddevZonAvgUnbal_opt(:,:)
1530 real(8), pointer, optional, intent(in) :: stddev3dUnbal_opt(:,:,:)
1531
1532 ! Locals:
1533 type(struct_gsv) :: stateVector
1534 real(8) :: dfact, zbufyz(nj,max(nLevEns_M,nLevens_T)), zbufy(nj)
1535 integer :: latIndex, levIndex, ierr, varIndex, varIndexStddev, nLevEns, numVarToWrite
1536 integer :: ip1, ip2, ip3, idatyp, idateo, numBits
1537 integer :: nulstats
1538 real(8), pointer :: field(:,:,:)
1539 integer, allocatable :: dateStampList(:)
1540 character(len=4) :: nomVarToWrite(1:20)
1541 integer :: fstouv, fnom, fstfrm, fclos
1542
1543 numBits = 32
1544 idatyp = 5
1545 ip1 = 0
1546 ip2 = 0
1547 ip3 = nens
1548 idateo = 0
1549
1550 ! figure out full list of variables to write
1551 numVarToWrite = size(nomvar,1)
1552 nomVarToWrite(1:numVarToWrite) = nomvar(:,cvSpace)
1553 if (present(stddevZonAvgUnbal_opt) .and. present(stddev3dUnbal_opt)) then
1554 do varIndex = 1, nvar
1555 if( all(nomvar(:,cvSpace) /= nomvar(varIndex,cvUnbalSpace)) ) then
1556 numVarToWrite = numVarToWrite + 1
1557 nomVarToWrite(numVarToWrite) = nomvar(varIndex,cvUnbalSpace)
1558 end if
1559 end do
1560 end if
1561
1562 ! first, write stddev3d
1563
1564 ! allocate stateVector used for writing stddev3d
1565 allocate(dateStampList(1))
1566 dateStampList(:) = 0
1567 call gsv_allocate( stateVector, 1, hco_ens, vco_ens, &
1568 mpi_local_opt=.true., dateStampList_opt=dateStampList, &
1569 varNames_opt=nomVarToWrite(1:numVarToWrite) )
1570 do varIndex = 1, numVarToWrite
1571 nLevEns = gsv_getNumLevFromVarName(stateVector,nomVarToWrite(varIndex))
1572 call gsv_getField(stateVector,field,nomVarToWrite(varIndex))
1573 if ( any(nomVarToWrite(varIndex) == nomvar(:,cvSpace)) ) then
1574 field(:,:,:) = stddev3d(:, :, (varLevOffset(varIndex)+1):(varLevOffset(varIndex)+nLevEns) )
1575 else if(present(stddevZonAvgUnbal_opt) .and. present(stddev3dUnbal_opt)) then
1576 varIndexStddev = findloc(nomvar(:,cvUnbalSpace), value=nomVarToWrite(varIndex), dim=1)
1577 field(:,:,:) = stddev3dUnbal_opt(:, :, (varLevOffset(varIndexStddev)+1):(varLevOffset(varIndexStddev)+nLevEns) )
1578 end if
1579 end do
1580 call gio_writeToFile(stateVector, './stddev.fst', etiket_in='STDDEV3D', ip3_opt=ip3, &
1581 typvar_opt='E', numBits_opt=numBits, containsFullField_opt=.false.)
1582
1583 ! second, write stddev (zonal average)
1584
1585 if (mmpi_myid == 0) then
1586 nulstats=0
1587 ierr = fnom (nulstats,'./stddev.fst','RND',0)
1588 ierr = fstouv(nulstats,'RND')
1589
1590 ! do 3d variables
1591 do varIndex=1,nvar3d
1592 nLevEns = gsv_getNumLevFromVarName(stateVector,nomvar(varIndex,cvSpace))
1593
1594 dfact=1.0d0
1595 do levIndex = 1, nlevEns
1596 do latIndex = 1, nj
1597 zbufyz(latIndex,levIndex)=dfact*stddevZonAvg(varLevOffset(varIndex)+levIndex,latIndex)
1598 end do
1599 end do
1600 ierr = utl_fstecr(zbufyz(:,1:nLevEns),-numBits,nulstats,idateo,0,0,1,nj,nlevEns,ip1,ip2,ip3, &
1601 'E',nomvar3d(varIndex,cvSpace),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1602
1603 if (present(stddevZonAvgUnbal_opt) .and. present(stddev3dUnbal_opt)) then
1604 if(nomvar3d(varIndex,cvSpace).ne.nomvar3d(varIndex,cvUnbalSpace)) then
1605 dfact=1.0d0
1606 do levIndex = 1, nlevEns
1607 do latIndex = 1, nj
1608 zbufyz(latIndex,levIndex)=dfact*stddevZonAvgUnbal_opt(varLevOffset(varIndex)+levIndex,latIndex)
1609 end do
1610 end do
1611 ierr = utl_fstecr(zbufyz(:,1:nLevEns),-numBits,nulstats,idateo,0,0,1,nj,nlevEns,ip1,ip2,ip3, &
1612 'E',nomvar3d(varIndex,cvUnbalSpace),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1613 end if
1614 end if
1615
1616 end do
1617
1618 ! now do 2D variables
1619 do varIndex=1,nvar2d
1620 if(nomvar2d(varIndex,cvSpace).eq.'P0') then
1621 dfact=1.0d0/1.0d2
1622 else
1623 dfact=1.0d0
1624 end if
1625
1626 do latIndex = 1, nj
1627 zbufy(latIndex)=dfact*stddevZonAvg(varLevOffset(nvar3d+1)+varIndex,latIndex)
1628 end do
1629 ierr = utl_fstecr(zbufy,-numBits,nulstats,idateo,0,0,1,nj,1,ip1,ip2,ip3, &
1630 'E',nomvar2d(varIndex,cvSpace),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1631
1632 if (present(stddevZonAvgUnbal_opt) .and. present(stddev3dUnbal_opt)) then
1633 if(nomvar2d(varIndex,cvSpace).ne.nomvar2d(varIndex,cvUnbalSpace)) then
1634 if(nomvar2d(varIndex,cvUnbalSpace).eq.'UP') then
1635 dfact=1.0d0/1.0d2
1636 else
1637 dfact=1.0d0
1638 end if
1639
1640 do latIndex = 1, nj
1641 zbufy(latIndex)=dfact*stddevZonAvgUnbal_opt(varLevOffset(nvar3d+1)+varIndex,latIndex)
1642 end do
1643 ierr = utl_fstecr(zbufy,-numBits,nulstats,idateo,0,0,1,nj,1,ip1,ip2,ip3, &
1644 'E',nomvar2d(varIndex,cvUnbalSpace),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1645 end if
1646 end if
1647
1648 end do
1649
1650 ierr = fstfrm(nulstats)
1651 ierr = fclos (nulstats)
1652 end if
1653
1654 call gsv_deallocate(stateVector)
1655
1656 write(*,*) 'finished writing stddev...'
1657
1658 end subroutine writeStddev
1659
1660 !--------------------------------------------------------------------------
1661 ! WRITESTDDEVBAL
1662 !--------------------------------------------------------------------------
1663 subroutine writeStddevBal(stddevZonAvgBal,stddev3dBal)
1664 !
1665 !:Purpose: Write the stddev of the balanced variables to a file
1666 !
1667 implicit none
1668
1669 ! Arguments:
1670 real(8), intent(in) :: stddevZonAvgBal(:,:)
1671 real(8), intent(in) :: stddev3dBal(:,:,:)
1672
1673 ! Locals:
1674 type(struct_gsv) :: stateVector
1675 real(8) :: dfact, zbufyz(nj,max(nLevEns_M,nLevens_T)), zbufy(nj)
1676 integer :: latIndex, levIndex, ierr, varIndex, nLevEns
1677 integer :: fstouv, fnom, fstfrm, fclos
1678 integer :: ip1, ip2, ip3, idatyp, idateo, numBits
1679 integer :: nulstats
1680 real(8), pointer :: field(:,:,:)
1681 integer, allocatable :: dateStampList(:)
1682 integer, parameter :: nvar3d=1, nvar2d=1, nvar=nvar3d+nvar2d
1683 character(len=4) :: nomVarToWrite(nvar)
1684 character(len=4) :: nomvar3dBal(nvar3d), nomvar2dBal(nvar2d)
1685 integer :: varLevOffsetBal(nvar)
1686
1687 nomvar3dBal(1)='TB'
1688 nomvar2dBal(1)='PB'
1689 varLevOffsetBal(1) = 0
1690 varLevOffsetBal(2) = 1*nLevEns_T
1691
1692 numBits = 32
1693 idatyp = 5
1694 ip1 = 0
1695 ip2 = 0
1696 ip3 = nens
1697 idateo = 0
1698
1699 nomVarToWrite(1:nvar3d) = nomvar3dBal(:)
1700 nomVarToWrite((1+nvar3d):nvar) = nomvar2dBal(:)
1701
1702 ! first, write stddev3d
1703
1704 ! allocate stateVector used for writing stddev3d
1705 allocate(dateStampList(1))
1706 dateStampList(:) = 0
1707 call gsv_allocate( stateVector, 1, hco_ens, vco_ens, &
1708 mpi_local_opt=.true., dateStampList_opt=dateStampList, &
1709 varNames_opt=nomVarToWrite(1:nvar) )
1710 do varIndex = 1, nvar
1711 nLevEns = gsv_getNumLevFromVarName(stateVector,nomVarToWrite(varIndex))
1712 call gsv_getField(stateVector,field,nomVarToWrite(varIndex))
1713 field(:,:,:) = stddev3dBal(:, :, (varLevOffsetBal(varIndex)+1):(varLevOffsetBal(varIndex)+nLevEns) )
1714 end do
1715 call gio_writeToFile(stateVector, './stddev_balanced.fst', etiket_in='STDDEV3D', ip3_opt=ip3, &
1716 typvar_opt='E', numBits_opt=numBits, containsFullField_opt=.false.)
1717
1718 ! second, write stddev (zonal average)
1719
1720 if (mmpi_myid == 0) then
1721 nulstats=0
1722 ierr = fnom (nulstats,'./stddev_balanced.fst','RND',0)
1723 ierr = fstouv(nulstats,'RND')
1724
1725 ! do 3d variables
1726 do varIndex=1,nvar3d
1727 nLevEns = gsv_getNumLevFromVarName(stateVector,nomVar3dBal(varIndex))
1728 !nip1_l(1:nLevEns_T)=nip1_T(1:nLevEns_T)
1729 dfact=1.0d0
1730 do levIndex = 1, nlevEns
1731 do latIndex = 1, nj
1732 zbufyz(latIndex,levIndex)=dfact*stddevZonAvgBal(varLevOffsetBal(varIndex)+levIndex,latIndex)
1733 end do
1734 end do
1735 ierr = utl_fstecr(zbufyz(:,1:nLevEns),-numBits,nulstats,idateo,0,0,1,nj,nlevEns,ip1,ip2,ip3, &
1736 'E',nomvar3dBal(varIndex),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1737 end do
1738
1739 ! now do 2D variables
1740 do varIndex=1,nvar2d
1741 dfact=1.0d0/1.0d2
1742 do latIndex = 1, nj
1743 zbufy(latIndex)=dfact*stddevZonAvgBal(varLevOffsetBal(nvar3d+1)+varIndex,latIndex)
1744 end do
1745 ierr = utl_fstecr(zbufy,-numBits,nulstats,idateo,0,0,1,nj,1,ip1,ip2,ip3, &
1746 'E',nomvar2dBal(varIndex),'STDDEV ','X',0,0,0,0,idatyp,.true.)
1747 end do
1748
1749 ierr = fstfrm(nulstats)
1750 ierr = fclos (nulstats)
1751 end if
1752
1753 call gsv_deallocate(stateVector)
1754
1755 write(*,*) 'finished writing stddev...'
1756
1757 end subroutine writeStddevBal
1758
1759 !--------------------------------------------------------------------------
1760 ! SPECTRALFILTER
1761 !--------------------------------------------------------------------------
1762 subroutine spectralFilter(ensPerturbations,nlev,waveBandIndex_opt)
1763 !
1764 !:Purpose: Apply a spectral filter
1765 !
1766 implicit none
1767
1768 ! Arguments:
1769 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
1770 integer, intent(in) :: nlev
1771 integer, optional, intent(in) :: waveBandIndex_opt
1772
1773 ! Locals:
1774 real(8) :: spectralState(nla_mpilocal,2,nlev)
1775 real(8) :: member(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlev)
1776 integer :: ensIndex, jk, jn, jm, ila_mpilocal, ila_mpiglobal
1777 real(8), allocatable :: ResponseFunction(:)
1778 real(8) :: waveLength
1779 character(len=128) :: outfilename
1780 character(len=2) :: wbnum
1781
1782 if ( nWaveBand /= 1 ) then
1783 write(*,*) 'Bandpass filtering step'
1784 if (.not. present(waveBandIndex_opt)) then
1785 write(*,*) 'Error: No waveBandIndex was supplied!!!'
1786 call utl_abort('calbmatrix_glb')
1787 end if
1788 allocate(ResponseFunction(0:ntrunc))
1789 write(wbnum,'(I2.2)') waveBandIndex_opt
1790 outfilename = "./ResponseFunction_"//wbnum//".txt"
1791 open (unit=99,file=outfilename,action="write",status="new")
1792 do jn = 0, ntrunc
1793 ResponseFunction(jn) = spf_filterResponseFunction(dble(jn),waveBandIndex_opt, waveBandPeaks, nWaveBand)
1794 if ( jn /= 0) then
1795 waveLength=2*MPC_PI_R8*ec_ra/dble(jn)
1796 else
1797 waveLength=0.d0
1798 end if
1799 write(* ,'(I4,2X,F7.1,2X,F5.3)') jn, waveLength/1000.d0, ResponseFunction(jn)
1800 write(99,'(I4,2X,F7.1,2X,F5.3)') jn, waveLength/1000.d0, ResponseFunction(jn)
1801 end do
1802 close(unit=99)
1803 end if
1804
1805 do ensIndex=1,nens
1806 member(:,:,:)=dble(ensPerturbations(:,:,:,ensIndex))
1807 if(nlev.eq.nkgdimEns) then
1808 call gst_setID(gstID_nkgdimEns)
1809 else if(nlev.eq.nLevEns_T+1) then
1810 call gst_setID(gstID_nLevEns_T_P1)
1811 else
1812 write(*,*) 'spectralFilter: nlev = ',nlev
1813 call utl_abort('spectralFilter: spectral transform not initialized for this number of levels')
1814 end if
1815 call gst_reespe(spectralState,member)
1816 if ( nWaveBand /= 1 ) then
1817 !$OMP PARALLEL DO PRIVATE (jk,jn,jm,ila_mpilocal,ila_mpiglobal)
1818 do jk = 1, nlev
1819 do jn = mynBeg, mynEnd, mynSkip
1820 do jm = mymBeg, mymEnd, mymSkip
1821 if(jm.le.jn) then
1822 ila_mpiglobal = gst_getNind(jm,gstID_nkgdimEns) + jn - jm
1823 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1824 spectralState(ila_mpilocal,1,jk) = spectralState(ila_mpilocal,1,jk) * ResponseFunction(jn)
1825 spectralState(ila_mpilocal,2,jk) = spectralState(ila_mpilocal,2,jk) * ResponseFunction(jn)
1826 end if
1827 end do
1828 end do
1829 end do
1830 !$OMP END PARALLEL DO
1831 end if
1832 call gst_speree(spectralState,member)
1833 ensPerturbations(:,:,:,ensIndex)=sngl(member(:,:,:))
1834 end do
1835
1836 if ( nWaveBand /= 1 ) then
1837 deallocate(ResponseFunction)
1838 end if
1839
1840 write(*,*) 'finished applying spectral filter...'
1841
1842 end subroutine spectralFilter
1843
1844 !--------------------------------------------------------------------------
1845 ! SPECTRALFILTER2
1846 !--------------------------------------------------------------------------
1847 subroutine spectralFilter2(ensPerts_in,ensPerts_out,waveBandIndex_opt)
1848 !
1849 !:Purpose: Apply a spectral filter
1850 !
1851 implicit none
1852
1853 ! Arguments:
1854 type(struct_ens), intent(inout) :: ensPerts_in
1855 type(struct_ens), intent(inout) :: ensPerts_out
1856 integer, optional, intent(in) :: waveBandIndex_opt
1857
1858 ! Locals:
1859 real(8), allocatable :: ensPertSP(:,:,:)
1860 real(8), allocatable :: ensPertGD(:,:,:)
1861 real(4), pointer :: ptr4d_r4(:,:,:,:)
1862 integer :: memberIndex, levIndex, latIndex, lonIndex
1863 integer :: jn, jm, ila_mpilocal, ila_mpiglobal
1864 integer :: gstFilterID
1865 real(8), allocatable :: ResponseFunction(:)
1866 real(8) :: waveLength
1867 character(len=128) :: outfilename
1868 character(len=2) :: wbnum
1869
1870 !
1871 !- 1. Pre-compute the response function (if needed)
1872 !
1873 if ( nWaveBand /= 1 ) then
1874 write(*,*) 'Bandpass filtering step'
1875 if (.not. present(waveBandIndex_opt)) then
1876 write(*,*) 'Error: No waveBandIndex was supplied!!!'
1877 call utl_abort('calbmatrix_glb: spectralFilter2')
1878 end if
1879 allocate(ResponseFunction(0:ntrunc))
1880 write(wbnum,'(I2.2)') waveBandIndex_opt
1881 outfilename = "./ResponseFunction_"//wbnum//".txt"
1882 if (mmpi_myid == 0) then
1883 open (unit=99,file=outfilename,action="write",status="new")
1884 end if
1885 do jn = 0, ntrunc
1886 ResponseFunction(jn) = spf_filterResponseFunction(dble(jn),waveBandIndex_opt, waveBandPeaks, nWaveBand)
1887 if ( jn /= 0) then
1888 waveLength=2*MPC_PI_R8*ec_ra/dble(jn)
1889 else
1890 waveLength=0.d0
1891 end if
1892 write(* ,'(I4,2X,F7.1,2X,F5.3)') jn, waveLength/1000.d0, ResponseFunction(jn)
1893 if (mmpi_myid == 0) then
1894 write(99,'(I4,2X,F7.1,2X,F5.3)') jn, waveLength/1000.d0, ResponseFunction(jn)
1895 end if
1896 end do
1897 if (mmpi_myid == 0) then
1898 close(unit=99)
1899 end if
1900 end if
1901
1902 !
1903 !- 2. Apply Filter
1904 !
1905 gstFilterID = gst_setup(ni,nj,nTrunc,nEnsOverDimension)
1906
1907 allocate(ensPertSP(nla_mpilocal,2,nEnsOverDimension))
1908 allocate(ensPertGD(nEnsOverDimension,myLonBeg:myLonEnd,myLatBeg:myLatEnd))
1909 ensPertSP(:,:,:) = 0.0d0
1910
1911 do levIndex = 1, ens_getNumK(ensPerts_in) ! Loop on variables and vertical levels
1912
1913 ptr4d_r4 => ens_getOneLev_r4(ensPerts_in,levIndex)
1914
1915 do latIndex = myLatBeg, myLatEnd
1916 ensPertGD(:,:,latIndex) = 0.0d0
1917 end do
1918
1919 !$OMP PARALLEL DO PRIVATE (memberIndex,latIndex,lonIndex)
1920 do latIndex = myLatBeg, myLatEnd
1921 do lonIndex = myLonBeg, myLonEnd
1922 do memberIndex = 1, nEns
1923 ensPertGD(memberIndex,lonIndex,latIndex) = dble(ptr4d_r4(memberIndex,1,lonIndex,latIndex))
1924 end do
1925 end do
1926 end do
1927 !$OMP END PARALLEL DO
1928
1929 !- GridPoint space -> Spectral Space
1930 call gst_setID(gstFilterID) ! IN
1931 call gst_reespe_kij(ensPertSP, & ! OUT
1932 ensPertGD) ! IN
1933
1934 !- Filtering
1935 if ( nWaveBand /= 1 ) then
1936 !$OMP PARALLEL DO PRIVATE (memberIndex,jn,jm,ila_mpilocal,ila_mpiglobal)
1937 do memberIndex = 1, nEns
1938 do jn = mynBeg, mynEnd, mynSkip
1939 do jm = mymBeg, mymEnd, mymSkip
1940 if (jm <= jn) then
1941 ila_mpiglobal = gst_getNind(jm,gstFilterID) + jn - jm
1942 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1943 ensPertSP(ila_mpilocal,1,memberIndex) = ensPertSP(ila_mpilocal,1,memberIndex) * ResponseFunction(jn)
1944 ensPertSP(ila_mpilocal,2,memberIndex) = ensPertSP(ila_mpilocal,2,memberIndex) * ResponseFunction(jn)
1945 end if
1946 end do
1947 end do
1948 end do
1949 !$OMP END PARALLEL DO
1950 end if
1951
1952 ! Spectral Space -> GridPoint space
1953 call gst_setID(gstFilterID) ! IN
1954 call gst_speree_kij(ensPertSP, & ! IN
1955 ensPertGD) ! OUT
1956
1957 ptr4d_r4 => ens_getOneLev_r4(ensPerts_out,levIndex)
1958
1959 !$OMP PARALLEL DO PRIVATE (memberIndex,latIndex,lonIndex)
1960 do latIndex = myLatBeg, myLatEnd
1961 do lonIndex = myLonBeg, myLonEnd
1962 do memberIndex = 1, nEns
1963 ptr4d_r4(memberIndex,1,lonIndex,latIndex) = sngl(ensPertGD(memberIndex,lonIndex,latIndex))
1964 end do
1965 end do
1966 end do
1967 !$OMP END PARALLEL DO
1968
1969 end do
1970
1971 if ( nWaveBand /= 1 ) then
1972 deallocate(ResponseFunction)
1973 end if
1974 deallocate(ensPertGD)
1975 deallocate(ensPertSP)
1976
1977 write(*,*) 'finished applying spectral filter...'
1978
1979 end subroutine spectralFilter2
1980
1981 !--------------------------------------------------------------------------
1982 ! CALCTHETA
1983 !--------------------------------------------------------------------------
1984 subroutine calcTheta(ensPerturbations,theta)
1985 !
1986 !:Purpose: Calculate the Theta turning angle according to Ekman balance
1987 !
1988 implicit none
1989
1990 ! Arguments:
1991 real(4), pointer, intent(in) :: ensPerturbations(:,:,:,:)
1992 real(8), intent(out) :: theta(:,:)
1993
1994 ! Locals:
1995 real(8) :: zchipsi(nLevEns_M,nj), zpsipsi(nLevEns_M,nj)
1996 real(8) :: zchipsi_mpiglobal(nLevEns_M,nj), zpsipsi_mpiglobal(nLevEns_M,nj)
1997 real(4), pointer :: psi_ptr(:,:,:), chi_ptr(:,:,:)
1998 integer :: latIndex,lonIndex,levIndex,ensIndex, ierr, nsize
1999
2000 theta(:,:) = 0.0d0
2001 zchipsi(:,:) = 0.0d0
2002 zpsipsi(:,:) = 0.0d0
2003 zchipsi_mpiglobal(:,:) = 0.0d0
2004 zpsipsi_mpiglobal(:,:) = 0.0d0
2005
2006 do ensIndex = 1,nens
2007 psi_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,1:nlevEns_M,ensIndex)
2008 chi_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,(nlevEns_M+1):(2*nlevEns_M),ensIndex)
2009
2010 ! update zchipsi and zpsipsi covariances
2011 do latIndex = myLatBeg, myLatEnd
2012 do lonIndex = myLonBeg, myLonEnd
2013 do levIndex = 1, nLevEns_M
2014 zpsipsi(levIndex,latIndex) = zpsipsi(levIndex,latIndex) + &
2015 psi_ptr(lonIndex,latIndex,levIndex) * psi_ptr(lonIndex,latIndex,levIndex)
2016 zchipsi(levIndex,latIndex) = zchipsi(levIndex,latIndex) + &
2017 chi_ptr(lonIndex,latIndex,levIndex) * psi_ptr(lonIndex,latIndex,levIndex)
2018 end do
2019 end do
2020 end do
2021 end do
2022
2023 nsize = nLevEns_M*nj
2024 call rpn_comm_allreduce(zchipsi,zchipsi_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2025 call rpn_comm_allreduce(zpsipsi,zpsipsi_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2026
2027 ! calculate THETA
2028 do latIndex = 1, nj
2029 do levIndex = 1, nLevEns_M
2030 theta(levIndex,latIndex) = atan(-zchipsi_mpiglobal(levIndex,latIndex) / &
2031 zpsipsi_mpiglobal(levIndex,latIndex))
2032 end do
2033 end do
2034
2035 write(*,*) 'finished computing theta...'
2036
2037 end subroutine calcTheta
2038
2039 !--------------------------------------------------------------------------
2040 ! CALCPTOT
2041 !--------------------------------------------------------------------------
2042 subroutine calcPtoT(ensPerturbations,PtoT)
2043 !
2044 !:Purpose: Calculate the "P" to Temperature transform matrix using
2045 ! a regression analysis of the "P" and temperature samples
2046 !
2047 implicit none
2048
2049 ! Arguments:
2050 real(4), pointer, intent(in) :: ensPerturbations(:,:,:,:)
2051 real(8), intent(out) :: PtoT(:,:,:)
2052
2053 ! Locals:
2054 real(8) :: spectralState(nla_mpilocal,2,nLevEns_M), spBalancedP(nla_mpilocal,2,nlevEns_M)
2055 real(8) :: balancedP(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlevEns_M)
2056 real(8) :: psi(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLevEns_M)
2057 real(4), pointer :: tt_ptr(:,:,:),ps_ptr(:,:)
2058 INTEGER :: ensIndex, JK1, JK2, nsize
2059 INTEGER :: IERR, JK, latIndex, lonIndex, JB, JPNLATBND
2060 PARAMETER (JPNLATBND = 3)
2061 REAL(8) :: ZFACT,zlat(nj)
2062 REAL(8) :: ZFACTTOT
2063 REAL(8) :: ZM1(NLEVENS_T+1,NLEVENS_M,JPNLATBND), ZM2(NLEVPTOT,NLEVPTOT,JPNLATBND)
2064 REAL(8) :: ZM1_mpiglobal(NLEVENS_T+1,NLEVENS_M,JPNLATBND), ZM2_mpiglobal(NLEVPTOT,NLEVPTOT,JPNLATBND)
2065 REAL(8) :: ZPTOTBND(NLEVENS_T+1,NLEVENS_M)
2066 REAL(8) :: ZM2INV(NLEVPTOT,NLEVPTOT,JPNLATBND)
2067 REAL(8) :: DLA2, DL1SA2
2068 REAL(8) :: DLLATMIN(JPNLATBND), DLLATMAX(JPNLATBND)
2069 real(8) :: zeigwrk(4*nlevPtoT),zeigen(nlevPtoT,nlevPtoT),zeigenv(nlevPtoT)
2070 real(8) :: zeigenvi(nlevPtoT)
2071 integer :: iwork,info
2072
2073 DATA DLLATMIN / -60.0D0, -30.0D0, 30.0D0 /
2074 DATA DLLATMAX / -30.0D0, 30.0D0, 60.0D0 /
2075
2076 DLA2 = ec_ra * ec_ra
2077 DL1SA2 = 1.D0/DLA2
2078
2079 ! 1. Initialize P_to_T, ZM1, ZM2
2080
2081 ZFACTTOT = 0.0D0
2082 DO latIndex = 1, nj
2083 ZFACTTOT = ZFACTTOT + cos(GST_GETRLATI(latIndex))
2084 END DO
2085 ZFACTTOT = NJ/ZFACTTOT
2086
2087 PtoT(:,:,:) = 0.0d0
2088 ZM1(:,:,:) = 0.0d0
2089 ZM2(:,:,:) = 0.0d0
2090 ZPTOTBND(:,:) = 0.0d0
2091 ZM1_mpiglobal(:,:,:) = 0.0d0
2092 ZM2_mpiglobal(:,:,:) = 0.0d0
2093
2094 do ensIndex = 1, nens
2095
2096 write(*,*) 'calcPtoT: processing member ',ensIndex
2097
2098 psi(:,:,:) = ensPerturbations(:,:,1:nlevEns_M,ensIndex)
2099 call gst_setID(gstID_nLevEns_M)
2100 call gst_reespe(spectralState,psi)
2101 call calcBalancedP(spectralState,spBalancedP)
2102 call gst_speree(spBalancedP,balancedP)
2103
2104 tt_ptr(myLonBeg:,myLatBeg:,1:) => ensPerturbations(:,:,(2*nLevEns_M+1):(2*nLevEns_M+nLevEns_T),ensIndex)
2105 ps_ptr(myLonBeg:,myLatBeg:) => ensPerturbations(:,:,2*nLevEns_M+2*nLevEns_T+1,ensIndex)
2106
2107 do latIndex = myLatBeg, myLatEnd
2108 zlat(latIndex)=GST_GETRLATI(latIndex)
2109 end do
2110
2111 !$OMP PARALLEL DO PRIVATE (JK1,JB,latIndex,ZFACT,lonIndex,JK2)
2112 DO JK1 = 1, (nLevEns_T+1)
2113 DO JB=1,JPNLATBND
2114 DO latIndex = myLatBeg, myLatEnd
2115 if ((ZLAT(latIndex) .gt. 2.D0*MPC_PI_R8*DLLATMIN(JB)/360.D0) .and. &
2116 (ZLAT(latIndex) .le. 2.D0*MPC_PI_R8*DLLATMAX(JB)/360.D0)) then
2117 ZFACT = cos(ZLAT(latIndex))*ZFACTTOT
2118 DO lonIndex = myLonBeg, myLonEnd
2119
2120 ! update ZM1 = sum_over_t_x_y[vec(T lnPs) vec(P_b)^T]
2121 DO JK2 = 1, nLevEns_M
2122 IF(JK1.LE.nLevEns_T) THEN
2123 zm1(jk1,jk2,jb) = zm1(jk1,jk2,jb) + zfact * tt_ptr(lonIndex,latIndex,jk1) * balancedP(lonIndex,latIndex,jk2)
2124 ELSE
2125 zm1(jk1,jk2,jb) = zm1(jk1,jk2,jb) + zfact * ps_ptr(lonIndex,latIndex) * balancedP(lonIndex,latIndex,jk2)
2126 END IF
2127 END DO
2128
2129 ! update ZM2 = sum_over_t_x_y[vec(P_b) vec(P_b)^T]
2130 IF(JK1.LE.NLEVPTOT) THEN
2131 DO JK2 = 1, NLEVPTOT
2132 zm2(jk1,jk2,jb) = zm2(jk1,jk2,jb) + zfact * balancedP(lonIndex,latIndex,jk1) * balancedP(lonIndex,latIndex,jk2)
2133 END DO
2134 END IF
2135
2136 END DO
2137 end if
2138 END DO
2139 END DO ! Loop on JPNLATBND
2140 END DO ! Loop on JK1
2141 !$OMP END PARALLEL DO
2142
2143 END DO
2144
2145 ! communicate matrices to have global result on all tasks
2146 nsize = (NLEVENS_T+1)*NLEVENS_M*JPNLATBND
2147 call rpn_comm_allreduce(zm1,zm1_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2148 nsize = NLEVPTOT*NLEVPTOT*JPNLATBND
2149 call rpn_comm_allreduce(zm2,zm2_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2150
2151 ! SET ZM1_MPIGLOBAL, ZM2_MPIGLOBAL EQUAL FOR ALL THREE REGIONS
2152 DO JK1 = 1, NLEVPTOT
2153 DO JK2 = 1, NLEVPTOT
2154 ZM2_MPIGLOBAL(JK1,JK2,1)=ZM2_MPIGLOBAL(JK1,JK2,1)+ZM2_MPIGLOBAL(JK1,JK2,3)
2155 ZM2_MPIGLOBAL(JK1,JK2,2)=ZM2_MPIGLOBAL(JK1,JK2,1)
2156 ZM2_MPIGLOBAL(JK1,JK2,3)=ZM2_MPIGLOBAL(JK1,JK2,1)
2157 END DO
2158 END DO
2159 DO JK1 = 1, (nLevEns_T+1)
2160 DO JK2 = 1, NLEVPTOT
2161 ZM1_MPIGLOBAL(JK1,JK2,1)=ZM1_MPIGLOBAL(JK1,JK2,1)+ZM1_MPIGLOBAL(JK1,JK2,3)
2162 ZM1_MPIGLOBAL(JK1,JK2,2)=ZM1_MPIGLOBAL(JK1,JK2,1)
2163 ZM1_MPIGLOBAL(JK1,JK2,3)=ZM1_MPIGLOBAL(JK1,JK2,1)
2164 END DO
2165 END DO
2166
2167 DO JK1=1,NLEVPTOT
2168 DO JK2=1,NLEVPTOT
2169 ZEIGEN(JK1,JK2)=ZM2_MPIGLOBAL(JK1,JK2,1)
2170 END DO
2171 END DO
2172 IWORK=4*NLEVPTOT
2173 CALL DSYEV('V','U',NLEVPTOT,ZEIGEN,NLEVPTOT,ZEIGENV,ZEIGWRK,IWORK,INFO)
2174
2175 write(*,*) 'calcPtot: info=',info
2176 write(*,*) 'calcPtot: eigen values=',zeigenv(:)
2177
2178 do JK1=1,NLEVPTOT
2179 if (ZEIGENV(JK1).gt.0.0d0) then
2180 ZEIGENVI(JK1)=1.0d0/ZEIGENV(JK1)
2181 else
2182 ZEIGENVI(JK1)=0.0d0
2183 end if
2184 end do
2185
2186 DO JK1=1,NLEVPTOT
2187 DO JK2=1,NLEVPTOT
2188 ZM2INV(JK1,JK2,1)=0.0d0
2189 DO JK=1,NLEVPTOT
2190 ZM2INV(JK1,JK2,1)=ZM2INV(JK1,JK2,1)+ZEIGEN(JK1,JK)*ZEIGENVI(JK)*ZEIGEN(JK2,JK)
2191 END DO
2192 END DO
2193 END DO
2194
2195
2196 ! Calculate A = ZM1_MPIGLOBAL*inv(ZM2)
2197 DO JK1 = 1, (nLevEns_T+1)
2198 DO JK2 = 1, NLEVPTOT
2199 DO JK = 1, NLEVPTOT
2200 ZPTOTBND(JK1,JK2) = ZPTOTBND(JK1,JK2) + ZM1_MPIGLOBAL(JK1,JK,1) * ZM2INV(JK,JK2,1)
2201 END DO
2202 END DO
2203 END DO
2204
2205 DO JK1 = 1, nLevEns_T+1
2206 DO JK2 = 1, NLEVPTOT
2207 DO latIndex = 1, nj
2208 PTOT(JK1,JK2,latIndex) = ZPTOTBND(JK1,JK2)
2209 END DO
2210 END DO
2211 END DO
2212
2213 write(*,*) 'finished computing PtoT...'
2214
2215 end subroutine calcPtoT
2216
2217 !--------------------------------------------------------------------------
2218 ! REMOVEGLOBALMEAN
2219 !--------------------------------------------------------------------------
2220 subroutine removeGlobalMean(ensPerturbations)
2221 !
2222 !:Purpose: Calculate and subtract the horizontal mean value from the ensemble
2223 ! of samples for each member, vertical level and variable
2224 !
2225 implicit none
2226
2227 ! Arguments:
2228 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
2229
2230 ! Locals:
2231 integer :: lonIndex,latIndex,levIndex,ensIndex,ierr
2232 real(8) :: dmean, dmean_mpiglobal
2233
2234 do ensIndex = 1, nens
2235 do levIndex = 1, nkgdimEns
2236 dmean = 0.0d0
2237 do latIndex = myLatBeg, myLatEnd
2238 do lonIndex = myLonBeg, myLonEnd
2239 dmean = dmean + ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)
2240 end do
2241 end do
2242 call rpn_comm_allreduce(dmean, dmean_mpiglobal,1,"mpi_double_precision","mpi_sum","GRID",ierr)
2243 dmean_mpiglobal = dmean_mpiglobal/(dble(ni)*dble(nj))
2244 do latIndex = myLatBeg, myLatEnd
2245 do lonIndex = myLonBeg, myLonEnd
2246 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex) = &
2247 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex) - dmean_mpiglobal
2248 end do
2249 end do
2250 end do
2251 end do
2252
2253 write(*,*) 'finished removing global mean...'
2254
2255 end subroutine removeGlobalMean
2256
2257 !--------------------------------------------------------------------------
2258 ! CALCZONAVG
2259 !--------------------------------------------------------------------------
2260 subroutine calcZonAvg(fieldsZonAvg_mpiglobal,fields3D,nlev)
2261 !
2262 !:Purpose: Calculate the zonal average of the supplied 3D fields
2263 !
2264 implicit none
2265
2266 ! Arguments:
2267 real(8), pointer, intent(inout) :: fieldsZonAvg_mpiglobal(:,:)
2268 real(8), pointer, intent(in) :: fields3D(:,:,:)
2269 integer, intent(in) :: nlev
2270
2271 ! Locals:
2272 integer :: lonIndex, latIndex, levIndex, ierr, nsize
2273 real(8) :: dfact
2274 real(8), allocatable :: fieldsZonAvg(:,:)
2275
2276 allocate(fieldsZonAvg(nlev,nj))
2277 fieldsZonAvg(:,:)=0.0d0
2278
2279 dfact=1.0d0/dble(ni)
2280 !$OMP PARALLEL DO PRIVATE (levIndex,latIndex,lonIndex)
2281 do levIndex = 1, nlev
2282 do latIndex = myLatBeg, myLatEnd
2283 do lonIndex = myLonBeg, myLonEnd
2284 fieldsZonAvg(levIndex,latIndex) = fieldsZonAvg(levIndex,latIndex) + &
2285 dfact*fields3D(lonIndex,latIndex,levIndex)
2286 end do
2287 end do
2288 end do
2289 !$OMP END PARALLEL DO
2290
2291 ! combine info from all mpi tasks
2292 nsize = nlev*nj
2293 call rpn_comm_allreduce(fieldsZonAvg,fieldsZonAvg_mpiglobal,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2294
2295 deallocate(fieldsZonAvg)
2296
2297 write(*,*) 'finished computing the zonal average...'
2298
2299 end subroutine calcZonAvg
2300
2301 !--------------------------------------------------------------------------
2302 ! CALCSTDDEV3D
2303 !--------------------------------------------------------------------------
2304 subroutine calcStddev3d(ensPerturbations,stddev3d,nlev)
2305 !
2306 !:Purpose: Calculate the 3d stddev field
2307 !
2308 implicit none
2309
2310 ! Arguments:
2311 real(8), pointer, intent(inout) :: stddev3d(:,:,:)
2312 real(4), pointer, intent(in) :: ensPerturbations(:,:,:,:)
2313 integer, intent(in) :: nlev
2314
2315 ! Locals:
2316 integer :: lonIndex,latIndex,levIndex,ensIndex
2317 real(8) :: dnens
2318
2319 write(*,*) 'started computing the stddev...'
2320 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2321
2322 stddev3d(:,:,:) = 0.0d0
2323 dnens = 1.0d0/dble(nens-1)
2324 !$OMP PARALLEL DO PRIVATE (levIndex,ensIndex,latIndex,lonIndex)
2325 do levIndex = 1, nlev
2326 do ensIndex = 1, nens
2327 do latIndex = myLatBeg, myLatEnd
2328 do lonIndex = myLonBeg, myLonEnd
2329 stddev3d(lonIndex,latIndex,levIndex)=stddev3d(lonIndex,latIndex,levIndex)+ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)**2
2330 end do
2331 end do
2332 end do
2333 do latIndex = myLatBeg, myLatEnd
2334 do lonIndex = myLonBeg, myLonEnd
2335 if(stddev3d(lonIndex,latIndex,levIndex).gt.0.0d0) then
2336 stddev3d(lonIndex,latIndex,levIndex)=sqrt(stddev3d(lonIndex,latIndex,levIndex)*dnens)
2337 else
2338 stddev3d(lonIndex,latIndex,levIndex)=0.0d0
2339 end if
2340 end do
2341 end do
2342 end do
2343 !$OMP END PARALLEL DO
2344
2345 write(*,*) 'finished computing the stddev...'
2346
2347 end subroutine calcStddev3d
2348
2349 !--------------------------------------------------------------------------
2350 ! CALCBALANCEDP
2351 !--------------------------------------------------------------------------
2352 subroutine calcBalancedP(sppsi,spgz)
2353 !
2354 !:Purpose: Calculate the balanced "P" variable using geostrophy
2355 !
2356 implicit none
2357
2358 ! Arguments:
2359 real(8), intent(in) :: sppsi(:,:,:)
2360 real(8), intent(out) :: spgz(:,:,:)
2361
2362 ! Locals:
2363 real(8) :: spvor_mpiglobal(nla,2,nlevEns_M)
2364 real(8) :: spvor_mpiglobal2(nla,2,nlevEns_M)
2365 real(8) :: spgz_mpiglobal(nla,2,nlevEns_M)
2366 integer :: ia, ib, ji, jm, levIndex, jla_mpilocal, ila_mpiglobal
2367 integer :: ierr, nsize
2368 real(8) :: zn, zm, zenm, zenmp1, zcon, dl1sa2
2369
2370 ! convert PSI to vorticity
2371 dl1sa2 = 1.0d0/(ec_ra*ec_ra)
2372 spvor_mpiglobal(:,:,:) = 0.0d0
2373 spvor_mpiglobal2(:,:,:) = 0.0d0
2374 do levIndex = 1, nlevEns_M
2375 do jla_mpilocal = 1, nla_mpilocal
2376 ila_mpiglobal = ilaList_mpiglobal(jla_mpilocal)
2377 spvor_mpiglobal(ila_mpiglobal,1,levIndex) = sppsi(jla_mpilocal,1,levIndex)*dl1sa2*gst_getRnnp1(ila_mpiglobal)
2378 spvor_mpiglobal(ila_mpiglobal,2,levIndex) = sppsi(jla_mpilocal,2,levIndex)*dl1sa2*gst_getRnnp1(ila_mpiglobal)
2379 end do
2380 ! ensure input field is zero for spectral component (0,0)
2381 spvor_mpiglobal(1,1,levIndex) = 0.0D0
2382 spvor_mpiglobal(1,2,levIndex) = 0.0D0
2383 end do
2384
2385 nsize = nla*2*nlevEns_M
2386 call rpn_comm_allreduce(spvor_mpiglobal,spvor_mpiglobal2,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
2387
2388 ! initialize output field to zero
2389 spgz_mpiglobal(:,:,:)=0.0d0
2390
2391 ! loop over levels and zonal wavenumbers
2392 ! n.b.: at the tip of the triangle, no contributions
2393
2394 zcon = -2.D0*ec_ROmega*ec_ra**2
2395 do levIndex = 1, nlevEns_M
2396
2397 ! the base address ia will point to the spherical harmonic
2398 ! coefficient (m,m), in the input field
2399 ia = 1
2400 do jm = 0, ntrunc-1
2401 ib = ia + ntrunc - jm
2402 zm = dble(jm)
2403
2404 ! at the base, contributions from n+1 coeff only
2405 zn = zm
2406 zenmp1 = sqrt ( ((zn+1)**2-zm**2)/(4.D0*(zn+1)**2-1.D0) )
2407 spgz_mpiglobal(ia,1,levIndex)=zcon*spvor_mpiglobal2(ia+1,1,levIndex)*zenmp1/((zn+1.0D0)**2)
2408 spgz_mpiglobal(ia,2,levIndex)=zcon*spvor_mpiglobal2(ia+1,2,levIndex)*zenmp1/((zn+1.0D0)**2)
2409
2410 zn = zn+1
2411 do ji = ia+1, ib-1
2412 zenm = sqrt ( (zn**2-zm**2)/(4.D0*zn**2-1.D0) )
2413 zenmp1 = sqrt ( ((zn+1)**2-zm**2)/(4.D0*(zn+1)**2-1.D0) )
2414 spgz_mpiglobal(ji,1,levIndex)=spvor_mpiglobal2(ji-1,1,levIndex)*zenm/(zn**2)
2415 spgz_mpiglobal(ji,2,levIndex)=spvor_mpiglobal2(ji-1,2,levIndex)*zenm/(zn**2)
2416 spgz_mpiglobal(ji,1,levIndex)=zcon*(spgz_mpiglobal(ji,1,levIndex)+spvor_mpiglobal2(ji+1,1,levIndex)*zenmp1/((zn+1.0D0)**2))
2417 spgz_mpiglobal(ji,2,levIndex)=zcon*(spgz_mpiglobal(ji,2,levIndex)+spvor_mpiglobal2(ji+1,2,levIndex)*zenmp1/((zn+1.0D0)**2))
2418 zn = zn + 1.0D0
2419 end do
2420
2421 ! at the top, contributions from n-1 coeff only
2422 zenm = sqrt ( (zn**2-zm**2)/(4.D0*zn**2-1.D0) )
2423 spgz_mpiglobal(ib,1,levIndex) = zcon*spvor_mpiglobal2(ib-1,1,levIndex)*zenm/(zn**2)
2424 spgz_mpiglobal(ib,2,levIndex) = zcon*spvor_mpiglobal2(ib-1,2,levIndex)*zenm/(zn**2)
2425 ia = ib + 1
2426 end do
2427 end do
2428
2429 ! ensure correct value for mass spectral-coefficient for m=n=0
2430 do levIndex = 1, nlevens_M
2431 spgz_mpiglobal(1,1,levIndex) = 0.0D0
2432 spgz_mpiglobal(1,2,levIndex) = 0.0D0
2433 end do
2434
2435 ! copy to mpilocal output array
2436 do levIndex = 1, nlevEns_M
2437 do jla_mpilocal = 1, nla_mpilocal
2438 ila_mpiglobal = ilaList_mpiglobal(jla_mpilocal)
2439 spgz(jla_mpilocal,1,levIndex) = spgz_mpiglobal(ila_mpiglobal,1,levIndex)
2440 spgz(jla_mpilocal,2,levIndex) = spgz_mpiglobal(ila_mpiglobal,2,levIndex)
2441 end do
2442 end do
2443
2444 end subroutine calcBalancedP
2445
2446 !--------------------------------------------------------------------------
2447 ! NORMALIZED3D
2448 !--------------------------------------------------------------------------
2449 subroutine normalize3d(ensPerturbations,stddev3d)
2450 !
2451 !:Purpose: Divide the ensemble perturbations by the supplied 3d stddev field
2452 !
2453 implicit none
2454
2455 ! Arguments:
2456 real(8), pointer, intent(in) :: stddev3d(:,:,:)
2457 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
2458
2459 ! Locals:
2460 integer :: lonIndex,latIndex,levIndex,ensIndex
2461 real(8) :: dfact
2462
2463 !$OMP PARALLEL DO PRIVATE (levIndex,ensIndex,latIndex,lonIndex,DFACT)
2464 do levIndex = 1, nkgdimEns
2465 do latIndex = myLatBeg, myLatEnd
2466 do lonIndex = myLonBeg, myLonEnd
2467 if(stddev3d(lonIndex,latIndex,levIndex).gt.0.0d0) then
2468 dfact=1.0d0/stddev3d(lonIndex,latIndex,levIndex)
2469 else
2470 dfact=0.0d0
2471 end if
2472 do ensIndex = 1, nens
2473 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)=ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)*dfact
2474 end do
2475 end do
2476 end do
2477 end do
2478 !$OMP END PARALLEL DO
2479
2480 write(*,*) 'finished normalizing by stddev3D...'
2481
2482 end subroutine normalize3d
2483
2484 !--------------------------------------------------------------------------
2485 ! MULTIPLY3D
2486 !--------------------------------------------------------------------------
2487 subroutine multiply3d(ensPerturbations,stddev3d,nlev)
2488 !
2489 !:Purpose: Multiply the ensemble perturbations by the supplied 3d stddev field
2490 !
2491 implicit none
2492
2493 ! Arguments:
2494 real(8), pointer, intent(in) :: stddev3d(:,:,:)
2495 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
2496 integer, intent(in) :: nlev
2497
2498 ! Locals:
2499 integer :: lonIndex,latIndex,levIndex,ensIndex
2500
2501 !$OMP PARALLEL DO PRIVATE (levIndex,ensIndex,latIndex,lonIndex)
2502 do ensIndex = 1, nens
2503 do levIndex = 1, nlev
2504 do latIndex = myLatBeg, myLatEnd
2505 do lonIndex = myLonBeg, myLonEnd
2506 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)=ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)*stddev3d(lonIndex,latIndex,levIndex)
2507 end do
2508 end do
2509 end do
2510 end do
2511 !$OMP END PARALLEL DO
2512
2513 write(*,*) 'finished multiplying by stddev3D...'
2514
2515 end subroutine multiply3d
2516
2517 !--------------------------------------------------------------------------
2518 ! READENSEMBLE
2519 !--------------------------------------------------------------------------
2520 subroutine readEnsemble(ensPerturbations)
2521 !
2522 !:Purpose: Read the ensemble of error samples from files
2523 !
2524 implicit none
2525
2526 ! Arguments:
2527 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
2528
2529 ! Locals:
2530 integer :: lonIndex, latIndex, levIndex, ensIndex, numStep
2531 integer, allocatable :: dateStampList(:)
2532 real(4), pointer :: field_r4(:,:,:)
2533 logical :: makeBiPeriodic
2534 type(struct_ens) :: ensPerts
2535 type(struct_gsv) :: stateVector
2536
2537 write(*,*) 'Before reading the ensemble:'
2538 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2539
2540 numStep = 1
2541 allocate(dateStampList(numStep))
2542 dateStampList(:) = -1
2543 call ens_allocate(ensPerts, nEns, numStep, hco_ens, vco_ens, dateStampList)
2544
2545 makeBiPeriodic = .false.
2546 call ens_readEnsemble(ensPerts, './ensemble', makeBiPeriodic, &
2547 containsFullField_opt=.false.)
2548
2549 call gsv_allocate(stateVector, 1, hco_ens, vco_ens, dateStampList_opt=dateStampList, &
2550 mpi_local_opt=.true., dataKind_opt=4)
2551
2552 do ensIndex = 1, nEns
2553 write(*,*) 'readEnsemble: copying over member ', ensIndex
2554 call ens_copyMember(ensPerts, stateVector, ensIndex)
2555
2556 call gsv_getField(stateVector,field_r4,'UU')
2557 do levIndex = 1, nLevEns_M
2558 do latIndex = myLatBeg, myLatEnd
2559 do lonIndex = myLonBeg, myLonEnd
2560 ensPerturbations(lonIndex,latIndex,levIndex+varLevOffset(1),ensIndex)= field_r4(lonIndex,latIndex,levIndex)
2561 end do
2562 end do
2563 end do
2564
2565 call gsv_getField(stateVector,field_r4,'VV')
2566 do levIndex = 1, nLevEns_M
2567 do latIndex = myLatBeg, myLatEnd
2568 do lonIndex = myLonBeg, myLonEnd
2569 ensPerturbations(lonIndex,latIndex,levIndex+varLevOffset(2),ensIndex)= field_r4(lonIndex,latIndex,levIndex)
2570 end do
2571 end do
2572 end do
2573
2574 call gsv_getField(stateVector,field_r4,'TT')
2575 do levIndex = 1, nLevEns_T
2576 do latIndex = myLatBeg, myLatEnd
2577 do lonIndex = myLonBeg, myLonEnd
2578 ensPerturbations(lonIndex,latIndex,levIndex+varLevOffset(3),ensIndex)= field_r4(lonIndex,latIndex,levIndex)
2579 end do
2580 end do
2581 end do
2582
2583 call gsv_getField(stateVector,field_r4,'LQ')
2584 do levIndex=1,nLevEns_T
2585 do latIndex = myLatBeg, myLatEnd
2586 do lonIndex = myLonBeg, myLonEnd
2587 ensPerturbations(lonIndex,latIndex,levIndex+varLevOffset(4),ensIndex) = field_r4(lonIndex,latIndex,levIndex)
2588 end do
2589 end do
2590 end do
2591
2592 call gsv_getField(stateVector,field_r4,'P0')
2593 do latIndex = myLatBeg, myLatEnd
2594 do lonIndex = myLonBeg, myLonEnd
2595 ensPerturbations(lonIndex,latIndex,1+varLevOffset(5),ensIndex)= field_r4(lonIndex,latIndex,1)
2596 end do
2597 end do
2598
2599 end do
2600
2601 call gsv_deallocate(stateVector)
2602 call ens_deallocate(ensPerts)
2603
2604 write(*,*) 'After reading the ensemble:'
2605 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2606
2607 write(*,*) 'finished reading ensemble members...'
2608
2609 end subroutine readEnsemble
2610
2611 !--------------------------------------------------------------------------
2612 ! UV_TO_PSICHI
2613 !--------------------------------------------------------------------------
2614 subroutine uv_to_psichi(ensPerturbations)
2615 !
2616 !:Purpose: Transform wind components to Psi and Chi
2617 !
2618 implicit none
2619
2620 ! Arguments:
2621 real(4), intent(inout) :: ensPerturbations(:,:,:,:)
2622
2623 ! Locals:
2624 integer :: ensIndex, levIndex, jla_mpilocal, ila_mpiglobal
2625 real(8) :: dla2
2626 real(8) :: spectralState(nla_mpilocal,2,nkgdimEns)
2627 real(8) :: member(myLonBeg:myLonEnd,myLatBeg:myLatend,nkgdimens)
2628
2629 ! Convert from U/V to PSI/CHI and spectrally filter all fields
2630 call utl_tmg_start(121,'--CalcStats_UVtoPsiChi')
2631 dla2 = ec_ra * ec_ra
2632 do ensIndex=1,nens
2633 write(*,*) ' doing u/v -> psi/chi and spectral filter for member ', ensIndex
2634 member(:,:,:)=dble(ensPerturbations(:,:,:,ensIndex))
2635 call gst_setID(gstID_nkgdimEns)
2636 call gst_gdsp(spectralState,member,nlevEns_M)
2637 do levIndex = 1, nlevEns_M
2638 do jla_mpilocal = 1, nla_mpilocal
2639 ila_mpiglobal = ilaList_mpiglobal(jla_mpilocal)
2640 spectralState(jla_mpilocal,1,levIndex) = spectralState(jla_mpilocal,1,levIndex) * dla2*gst_getR1snp1(ila_mpiglobal)
2641 spectralState(jla_mpilocal,2,levIndex) = spectralState(jla_mpilocal,2,levIndex) * dla2*gst_getR1snp1(ila_mpiglobal)
2642 spectralState(jla_mpilocal,1,levIndex+nlevEns_M) = spectralState(jla_mpilocal,1,levIndex+nlevEns_M) * dla2*gst_getR1snp1(ila_mpiglobal)
2643 spectralState(jla_mpilocal,2,levIndex+nlevEns_M) = spectralState(jla_mpilocal,2,levIndex+nlevEns_M) * dla2*gst_getR1snp1(ila_mpiglobal)
2644 end do
2645 end do
2646 call gst_speree(spectralState,member)
2647 ensPerturbations(:,:,:,ensIndex)=sngl(member(:,:,:))
2648 end do
2649
2650 call utl_tmg_stop(121)
2651 write(*,*) 'finished doing u/v -> psi/chi and spectral filter...'
2652
2653 end subroutine uv_to_psichi
2654
2655 !--------------------------------------------------------------------------
2656 ! REMOVEMEAN
2657 !--------------------------------------------------------------------------
2658 subroutine removeMean(ensPerturbations)
2659 !
2660 !:Purpose: Compute and subtract the ensemble mean
2661 !
2662 implicit none
2663
2664 ! Arguments:
2665 real(4), pointer, intent(inout) :: ensPerturbations(:,:,:,:)
2666
2667 ! Locals:
2668 integer :: ensIndex, levIndex, latIndex, lonIndex
2669 real(8) :: dnens, gd2d(myLonBeg:myLonEnd,myLatBeg:myLatEnd)
2670
2671 ! remove mean and divide by sqrt(2*(NENS-1)) - extra 2 is needed?
2672 dnens=1.0d0/dble(nens)
2673 !$OMP PARALLEL DO PRIVATE (levIndex,gd2d,ensIndex,latIndex,lonIndex)
2674 do levIndex = 1, nkgdimEns
2675 gd2d(:,:)=0.0d0
2676 do ensIndex = 1, nens
2677 do latIndex = myLatBeg, myLatEnd
2678 do lonIndex = myLonBeg, myLonEnd
2679 gd2d(lonIndex,latIndex)=gd2d(lonIndex,latIndex)+ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)
2680 end do
2681 end do
2682 end do
2683 do latIndex = myLatBeg, myLatEnd
2684 do lonIndex = myLonBeg, myLonEnd
2685 gd2d(lonIndex,latIndex)=gd2d(lonIndex,latIndex)*dnens
2686 end do
2687 end do
2688 do ensIndex = 1, nens
2689 do latIndex = myLatBeg, myLatEnd
2690 do lonIndex = myLonBeg, myLonEnd
2691 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)= &
2692 ensPerturbations(lonIndex,latIndex,levIndex,ensIndex)-gd2d(lonIndex,latIndex)
2693 end do
2694 end do
2695 end do
2696 end do
2697 !$OMP END PARALLEL DO
2698
2699 write(*,*) 'finished removing the ensemble mean...'
2700
2701 end subroutine removeMean
2702
2703 !--------------------------------------------------------------------------
2704 ! HORIZCORRELFUNCTION
2705 !--------------------------------------------------------------------------
2706 subroutine horizCorrelFunction(rstddev,variableType, waveBandIndex_opt)
2707 !
2708 !:Purpose: Compute homogeneous-isotropic horizontal correlation
2709 ! function from spectral variances
2710 !
2711 implicit none
2712
2713 ! Arguments:
2714 real(8), intent(in) :: rstddev(nkgdimEns,0:ntrunc)
2715 integer, intent(in) :: variableType
2716 integer,optional, intent(in) :: waveBandIndex_opt
2717
2718 ! Locals:
2719 real(8) :: spectralState(nla,2,nkgdimEns)
2720 real(8) :: gridState(ni,nj,nkgdimEns)
2721 integer :: ji, jk, jn, jm, ila, iref, jref
2722 integer :: nLevEns, nLevStart, nLevEnd, varIndex, iStart, iEnd
2723 character(len=128) :: outfilename
2724 character(len=2) :: wbnum
2725
2726 write(*,*)
2727 write(*,*) 'Computing horizontal correlation functions'
2728
2729 !
2730 !- 1. Spectral transform of a delta function (at the center of the domain)
2731 !
2732
2733 !- 1.1 Create the delta function
2734 iref = ni/2
2735 jref = nj/2
2736
2737 GridState(:,:,:) = 0.d0
2738 GridState(iref,jref,:) = 1.d0
2739
2740 !- 1.2 Adjoint of the identity (change of norm)
2741 GridState(iref,jref,:) = GridState(iref,jref,:) * real(ni,8) / gst_getRWT(jref)
2742
2743 !- 1.3 Move to spectral space
2744 call gst_setID(gstID_nkgdimEns)
2745 call gst_reespe(spectralState, & ! OUT
2746 GridState) ! IN
2747
2748 !
2749 !- 2. Apply the horizontal correlation function
2750 !
2751 !$OMP PARALLEL DO PRIVATE (jk,jn,jm,ila)
2752 do jk = 1, nkgdimEns
2753 do jn = 0, ntrunc
2754 do jm = 0, jn
2755 ila = gst_getNind(jm) + jn - jm
2756 if (jm.eq.0) then
2757 spectralState(ila,1,jk) = spectralState(ila,1,jk) * rstddev(jk,jn)**2 * 2.0d0
2758 spectralState(ila,2,jk) = 0.d0
2759 else
2760 spectralState(ila,1,jk) = spectralState(ila,1,jk) * rstddev(jk,jn)**2 * 2.0d0
2761 spectralState(ila,2,jk) = spectralState(ila,2,jk) * rstddev(jk,jn)**2 * 2.0d0
2762 end if
2763 end do
2764 end do
2765 end do
2766 !$OMP END PARALLEL DO
2767
2768 !
2769 !- 3. Move back to grid point space
2770 !
2771 call gst_setID(gstID_nkgdimEns)
2772 call gst_speree(spectralState, & ! IN
2773 GridState) ! OUT
2774
2775 !
2776 !- 4. Write to file
2777 !
2778 if (mmpi_myid == 0) then
2779 if ( nWaveBand /= 1 ) then
2780 if (.not. present(waveBandIndex_opt)) then
2781 write(*,*) 'horizCorrelFunction: No waveBandIndex was supplied!!!'
2782 call utl_abort('calbmatrix_glb')
2783 end if
2784 write(wbnum,'(I2.2)') waveBandIndex_opt
2785 end if
2786
2787 !- 4.1 2D correlation function in fst format
2788 if ( nWaveBand == 1 ) then
2789 outfilename = "./horizCorrel.fst"
2790 else
2791 outfilename = "./horizCorrel_"//wbnum//".fst"
2792 end if
2793 call write3d(GridState,outfilename,'HORIZCORFUNC',variableType)
2794
2795 !- 4.2 1D correlation function in txt format (for plotting purposes)
2796 iStart=iref
2797 iEnd=3*ni/4 ! About 10 000 km away from the center of the domain
2798
2799 do varIndex = 1, nvar3d
2800 if ( nWaveBand == 1 ) then
2801 outfilename = "./horizCorrel_"//trim(nomvar3d(varIndex,variableType))//".txt"
2802 else
2803 outfilename = "./horizCorrel_"//trim(nomvar3d(varIndex,variableType))//"_"//wbnum//".txt"
2804 end if
2805 open (unit=99,file=outfilename,action="write",status="new")
2806
2807 if(vnl_varLevelFromVarName(nomvar3d(varIndex,variableType)).eq.'MM') then
2808 nLevEns = nLevEns_M
2809 else
2810 nLevEns = nLevEns_T
2811 end if
2812 nLevStart = varLevOffset(varIndex)+ 1
2813 nLevEnd = varLevOffset(varIndex)+ nLevEns
2814
2815 do ji=iStart,iEnd
2816 do jk = nLevStart,nLevEnd
2817 if ( jk == nLevStart ) then
2818 write(99,'(I7,2X,F7.1,2X,F6.4,$)') ji-iStart, (ji-iStart)*gridSpacingInKm, GridState(ji,jref,jk)
2819 else if ( jk == nLevEnd ) then
2820 write(99,'(2X,F6.4)') GridState(ji,jref,jk) ! Saut de ligne
2821 else
2822 write(99,'(2X,F6.4,$)') GridState(ji,jref,jk)
2823 end if
2824 end do
2825 end do
2826 close(unit=99)
2827 end do
2828
2829 do varIndex = 1, nvar2d
2830 if ( nWaveBand == 1 ) then
2831 outfilename = "./horizCorrel_"//trim(nomvar2d(varIndex,variableType))//".txt"
2832 else
2833 outfilename = "./horizCorrel_"//trim(nomvar2d(varIndex,variableType))//"_"//wbnum//".txt"
2834 end if
2835 open (unit=99,file=outfilename,action="write",status="new")
2836 do ji=iStart,iEnd
2837 write(99,'(I7,2X,F7.1,2X,F6.4)') ji-iStart, (ji-iStart)*gridSpacingInKm, GridState(ji,jref,varLevOffset(nvar3d+1)+varIndex)
2838 end do
2839 close(unit=99)
2840 end do
2841
2842 end if
2843
2844 end subroutine horizCorrelFunction
2845
2846 !--------------------------------------------------------------------------
2847 ! WRITE3D
2848 !--------------------------------------------------------------------------
2849 subroutine write3d(gridpoint3d,filename,etiket_in,variableType)
2850 !
2851 !:Purpose: Write the 3D stddev fields for all variables
2852 !
2853 implicit none
2854
2855 ! Arguments:
2856 real(8), intent(in) :: gridpoint3d(:,:,:)
2857 character(len=*), intent(in) :: filename
2858 character(len=*), intent(in) :: etiket_in
2859 integer, intent(in) :: variableType
2860
2861 ! Locals:
2862 real(8) :: dfact,zbuf(ni,nj)
2863 integer latIndex,lonIndex,levIndex,ierr,varIndex,nLevEns
2864 integer fstouv,fnom,fstfrm,fclos
2865 integer ip1,ip2,ip3,idatyp,idateo,ipak,nip1_l(max(nLevEns_M,nLevens_T))
2866 integer :: nulstats
2867 character(len=12) :: etiket
2868
2869 etiket=trim(etiket_in)
2870
2871 nulstats=0
2872 ierr = fnom (nulstats,filename,'RND',0)
2873 ierr = fstouv(nulstats,'RND')
2874
2875 ipak = -32
2876 idatyp = 5
2877 ip1 = 0
2878 ip2 = 0
2879 ip3 = nens
2880 idateo = 0
2881
2882 ! do 3d variables
2883 do varIndex = 1, nvar3d
2884
2885 if(vnl_varLevelFromVarName(nomvar3d(varIndex,variableType)).eq.'MM') then
2886 nLevEns = nLevEns_M
2887 nip1_l(1:nLevEns_M)=nip1_M(1:nLevEns_M)
2888 else
2889 nLevEns = nLevEns_T
2890 nip1_l(1:nLevEns_T)=nip1_T(1:nLevEns_T)
2891 end if
2892 dfact=1.0d0
2893
2894 do levIndex = 1, nlevEns
2895 do latIndex = myLatBeg, myLatEnd
2896 do lonIndex = myLonBeg, myLonEnd
2897 zbuf(lonIndex,latIndex) = dfact*gridpoint3d(lonIndex,latIndex,varLevOffset(varIndex)+levIndex)
2898 end do
2899 end do
2900 ierr = utl_fstecr(zbuf,ipak,nulstats,idateo,0,0,ni,nj,1,nip1_l(levIndex),ip2,ip3, &
2901 'E',nomvar3d(varIndex,variableType),etiket,'G',0,0,0,0,idatyp,.true.)
2902 end do
2903
2904 end do
2905
2906 ! now do 2D variables
2907 do varIndex = 1, nvar2d
2908 !if(nomvar2d(varIndex,variableType).eq.'P0') then
2909 ! dfact=1.0d0/1.0d2
2910 !else
2911 dfact = 1.0d0
2912 !end if
2913
2914 do latIndex = myLatBeg, myLatEnd
2915 do lonIndex = myLonBeg, myLonEnd
2916 zbuf(lonIndex,latIndex) = dfact*gridpoint3d(lonIndex,latIndex,varLevOffset(nvar3d+1)+varIndex)
2917 end do
2918 end do
2919 ierr = utl_fstecr(zbuf,ipak,nulstats,idateo,0,0,ni,nj,1,0,ip2,ip3, &
2920 'E',nomvar2d(varIndex,variableType),etiket,'G',0,0,0,0,idatyp,.true.)
2921
2922 end do
2923
2924 ierr = fstfrm(nulstats)
2925 ierr = fclos (nulstats)
2926
2927 write(*,*) 'finished writing 3d array...'
2928
2929 end subroutine write3d
2930
2931 !--------------------------------------------------------------------------
2932 ! CALCHORIZSCALE
2933 !--------------------------------------------------------------------------
2934 subroutine CalcHorizScale(rstddev,variableType,waveBandIndex_opt)
2935 !
2936 !:Purpose: Calculate the horizontal correlation length scale
2937 !
2938 implicit none
2939
2940 ! Based on subroutine corrlength.ftn in the "old" var code
2941
2942 ! Arguments:
2943 real(8), intent(in) :: rstddev(nkgdimEns,0:ntrunc)
2944 integer, intent(in) :: variableType
2945 integer, optional, intent(in) :: waveBandIndex_opt
2946
2947 ! Locals:
2948 real(8) :: HorizScale(nkgdimEns)
2949 real(8), pointer :: PressureProfile(:)
2950 integer :: jk, jn, nLevEns, varIndex
2951 real(8) :: rjn, fact, temp, a, b
2952 character(len=128) :: outfilename
2953 character(len=2) :: wbnum
2954
2955 if (mmpi_myid /= 0) return
2956
2957 write(*,*)
2958 write(*,*) 'Computing horizontal correlation lengthscales'
2959
2960 !
2961 !- 1. Compute the horizontal correlation lengthscales
2962 !
2963 do jk = 1, nkgdimEns
2964 a = 0.d0
2965 b = 0.d0
2966 do jn = 0, ntrunc
2967 rjn = dble(jn)
2968 fact = (2.0d0*rjn + 1.0d0)/2.0d0
2969 temp = (rstddev(jk,jn)**2) * fact
2970 a = a + temp
2971 if (jn /= 0) then
2972 b = b - temp*rjn*(rjn+1.d0)
2973 end if
2974 end do
2975 if ( a > 0.d0 .and. b /= 0.d0 ) then
2976 HorizScale(jk) = ec_ra * sqrt(-2.0d0*a/b)
2977 else
2978 HorizScale(jk) = 0.d0
2979 end if
2980 end do
2981
2982 !
2983 !- 2. Write the results
2984 !
2985 if ( nWaveBand /= 1 ) then
2986 if (.not. present(waveBandIndex_opt)) then
2987 write(*,*) 'CalcHorizScale: No waveBandIndex was supplied!!!'
2988 call utl_abort('calbmatrix_glb')
2989 end if
2990 write(wbnum,'(I2.2)') waveBandIndex_opt
2991 end if
2992
2993 do varIndex = 1, nvar3d
2994 write(*,*)
2995 write(*,*) nomvar3d(varIndex,variableType)
2996
2997 if ( nWaveBand == 1 ) then
2998 outfilename = "./horizScale_"//trim(nomvar3d(varIndex,variableType))//".txt"
2999 else
3000 outfilename = "./horizScale_"//trim(nomvar3d(varIndex,variableType))//"_"//wbnum//".txt"
3001 end if
3002 open (unit=99,file=outfilename,action="write",status="new")
3003
3004 if(vnl_varLevelFromVarName(nomvar3d(varIndex,variableType)).eq.'MM') then
3005 nLevEns = nLevEns_M
3006 PressureProfile => pressureProfile_M
3007 else
3008 nLevEns = nLevEns_T
3009 PressureProfile => pressureProfile_T
3010 end if
3011 do jk=1,nlevEns
3012 write(* ,'(I3,2X,F6.1,2X,F6.1)') jk, PressureProfile(jk)/100.d0, HorizScale(varLevOffset(varIndex)+jk)/1000.d0
3013 write(99,'(I3,2X,F6.1,2X,F6.1)') jk, PressureProfile(jk)/100.d0, HorizScale(varLevOffset(varIndex)+jk)/1000.d0
3014 end do
3015
3016 close(unit=99)
3017 end do
3018
3019 do varIndex = 1, nvar2d
3020 write(*,*)
3021 write(*,*) nomvar2d(varIndex,variableType)
3022
3023 if ( nWaveBand == 1 ) then
3024 outfilename = "./horizScale_"//trim(nomvar2d(varIndex,variableType))//".txt"
3025 else
3026 outfilename = "./horizScale_"//trim(nomvar2d(varIndex,variableType))//"_"//wbnum//".txt"
3027 end if
3028 open (unit=99,file=outfilename,action="write",status="new")
3029
3030 write(* ,'(I3,2X,F6.1,2X,F6.1)') 1, 1010.0, HorizScale(varLevOffset(nvar3d+1)+varIndex)/1000.d0
3031 write(99,'(I3,2X,F6.1,2X,F6.1)') 1, 1010.0, HorizScale(varLevOffset(nvar3d+1)+varIndex)/1000.d0
3032
3033 close(unit=99)
3034 end do
3035
3036 end subroutine CalcHorizScale
3037
3038 !--------------------------------------------------------------------------
3039 ! WRITEPOWERSPEC
3040 !--------------------------------------------------------------------------
3041 subroutine writePowerSpec(powerSpec,variableType)
3042 !
3043 !:Purpose: Write the computed power spectrum to an ascii file
3044 !
3045 implicit none
3046
3047 ! Arguments:
3048 real(8), intent(in) :: powerSpec(nkgdimEns,0:ntrunc)
3049 integer, intent(in) :: variableType
3050
3051 ! Locals:
3052 integer :: jk, nLevEns, nLevStart, nLevEnd, varIndex, jn
3053 real(8) :: waveLength
3054 character(len=128) :: outfilename
3055
3056 if (mmpi_myid /= 0) return
3057
3058 !- Write to txt files
3059 do varIndex = 1, nvar3d
3060
3061 outfilename = "./PowerSpec_"//trim(nomvar3d(varIndex,variableType))//".txt"
3062 open (unit=99,file=outfilename,action="write",status="new")
3063
3064 if(vnl_varLevelFromVarName(nomvar3d(varIndex,variableType)).eq.'MM') then
3065 nLevEns = nLevEns_M
3066 else
3067 nLevEns = nLevEns_T
3068 end if
3069 nLevStart = varLevOffset(varIndex)+ 1
3070 nLevEnd = varLevOffset(varIndex)+ nLevEns
3071
3072 do jn = 0, ntrunc
3073 if ( jn /= 0) then
3074 waveLength=2*MPC_PI_R8*ec_ra/dble(jn)
3075 else
3076 waveLength=0.d0
3077 end if
3078 do jk = nLevStart,nLevEnd
3079 if ( jk == nLevStart ) then
3080 write(99,'(I4,2X,F7.1,2X,e10.3,$)') jn, waveLength/1000.d0, sngl(powerSpec(jk,jn))
3081 else if ( jk == nLevEnd ) then
3082 write(99,'(2X,e10.3)') sngl(powerSpec(jk,jn)) ! Saut de ligne
3083 else
3084 write(99,'(2X,e10.3,$)') sngl(powerSpec(jk,jn))
3085 end if
3086 end do
3087 end do
3088 close(unit=99)
3089 end do
3090
3091 do varIndex = 1, nvar2d
3092
3093 outfilename = "./PowerSpec_"//trim(nomvar2d(varIndex,variableType))//".txt"
3094 open (unit=99,file=outfilename,action="write",status="new")
3095 do jn = 0, ntrunc
3096 if ( jn /= 0) then
3097 waveLength=2*MPC_PI_R8*ec_ra/dble(jn)
3098 else
3099 waveLength=0.d0
3100 end if
3101 write(99,'(I4,2X,F7.1,2X,e10.3)') jn, waveLength/1000.d0, sngl(powerSpec(varLevOffset(nvar3d+1)+varIndex,jn))
3102 end do
3103 close(unit=99)
3104 end do
3105
3106 end subroutine writePowerSpec
3107
3108 !--------------------------------------------------------------------------
3109 ! calcLocalCorrelations (identical to the routine with the same name in calcstatslam)
3110 !--------------------------------------------------------------------------
3111 subroutine calcLocalCorrelations(ensPerts)
3112 !
3113 !:Purpose: Compute local horizontal correlations
3114 !
3115 implicit none
3116
3117 ! Arguments:
3118 type(struct_ens), intent(in) :: ensPerts
3119
3120 ! Locals:
3121 type(struct_gsv) :: statevector_locHorizCor
3122 type(struct_gsv) :: statevector_oneMember
3123 type(struct_gsv) :: statevector_oneMemberTiles
3124 real(8), pointer :: ptr3d_r8(:,:,:)
3125 real(8), pointer :: ptr3d_r8_oneMember(:,:,:)
3126 real(8) :: dnEns
3127 integer :: i, j, k, ens
3128 integer :: blocklength_x, blocklength_y
3129 integer :: iref_id, jref_id, iref, jref
3130 integer :: imin, imax, jmin, jmax
3131 character(len=4), pointer :: varNamesList(:)
3132 integer :: ierr, fclos, fnom, nulnam
3133
3134 ! Namelist variables
3135 integer :: blockpadding
3136 integer :: nirefpoint
3137 integer :: njrefpoint
3138
3139 NAMELIST /NAMHVCORREL_LOCAL/nirefpoint, njrefpoint, blockpadding
3140
3141 !
3142 ! To compute the local horizontal correlation for some 'reference' grid point
3143 ! ... we assume that the ensemble grid point mean was removed and that
3144 ! the ensemble values were divided by the grid point std dev.
3145 !
3146
3147 nirefpoint = 4 ! Number of reference grid point in x
3148 njrefpoint = 2 ! Number of reference grid point in y
3149 blockpadding = 4 ! Number of grid point padding between blocks (to set correlation to 0 between each block)
3150
3151 nulnam = 0
3152 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
3153 read(nulnam,nml=NAMHVCORREL_LOCAL)
3154 if (ierr /= 0) call utl_abort('calcLocalCorrelations: Error reading namelist NAMHVCORREL_LOCAL')
3155 if (mmpi_myid == 0) write(*,nml=NAMHVCORREL_LOCAL)
3156 ierr = fclos(nulnam)
3157
3158 blocklength_x = hco_ens%ni / nirefpoint ! Horizontal correlation will be compute blocklength x blocklength gridpoint
3159 ! around each reference point
3160 blocklength_y = hco_ens%nj / njrefpoint ! Horizontal correlation will be compute blocklength x blocklength gridpoint
3161 ! around each reference point
3162
3163 nullify(varNamesList)
3164 call ens_varNamesList(varNamesList,ensPerts)
3165
3166 call gsv_allocate(statevector_locHorizCor, ens_getNumStep(ensPerts), &
3167 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3168 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
3169 mpi_distribution_opt='VarsLevs', dataKind_opt=8 )
3170
3171 call gsv_allocate(statevector_oneMemberTiles, ens_getNumStep(ensPerts), &
3172 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3173 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
3174 mpi_distribution_opt='Tiles', dataKind_opt=8 )
3175
3176 call gsv_allocate(statevector_oneMember, ens_getNumStep(ensPerts), &
3177 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3178 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
3179 mpi_distribution_opt='VarsLevs', dataKind_opt=8 )
3180
3181 call gsv_zero(statevector_locHorizCor)
3182
3183 dnEns = 1.0d0/dble(nEns-1)
3184
3185 call gsv_getField(statevector_locHorizCor,ptr3d_r8)
3186
3187 do ens = 1, nEns
3188 call ens_copyMember(ensPerts, statevector_oneMemberTiles, ens)
3189 call gsv_transposeTilesToVarsLevs(statevector_oneMemberTiles, statevector_oneMember)
3190 call gsv_getField(statevector_oneMember,ptr3d_r8_oneMember)
3191
3192 do k = statevector_locHorizCor%mykBeg, statevector_locHorizCor%mykEnd
3193 do jref_id = 1, njrefpoint
3194 do iref_id = 1, nirefpoint
3195 iref = (2*iref_id-1)*blocklength_x/2
3196 jref = (2*jref_id-1)*blocklength_y/2
3197 jmin = max(jref-(blocklength_y-blockpadding)/2,1)
3198 jmax = min(jref+(blocklength_y-blockpadding)/2,hco_ens%nj)
3199 imin = max(iref-(blocklength_x-blockpadding)/2,1)
3200 imax = min(iref+(blocklength_x-blockpadding)/2,hco_ens%ni)
3201 do j = jmin, jmax
3202 do i = imin, imax
3203 ptr3d_r8(i,j,k) = ptr3d_r8(i,j,k) + &
3204 ptr3d_r8_oneMember(i,j,k)*ptr3d_r8_oneMember(iref,jref,k)
3205 end do
3206 end do
3207 end do
3208 end do
3209 end do
3210
3211 end do
3212
3213 call gsv_scale(statevector_locHorizCor,dnEns)
3214
3215 write(*,*) 'finished computing the local horizontal correlations...'
3216
3217 !
3218 !- 4. Write to file
3219 !
3220 call gio_writeToFile(statevector_locHorizCor, './horizCorrelLocal.fst', 'HCORREL_LOC', &
3221 typvar_opt = 'E', numBits_opt = 32)
3222
3223 call gsv_deallocate(statevector_locHorizCor)
3224 call gsv_deallocate(statevector_oneMember)
3225 call gsv_deallocate(statevector_oneMemberTiles)
3226
3227 end subroutine calcLocalCorrelations
3228
3229 !--------------------------------------------------------------------------
3230 ! calcLocalVertCorrMatrix
3231 !--------------------------------------------------------------------------
3232 subroutine calcLocalVertCorrMatrix(ensPerts)
3233 !
3234 !:Purpose: Compute all vertical and between-variable local correlations
3235 !
3236 implicit none
3237
3238 ! Arguments:
3239 type(struct_ens), intent(in) :: ensPerts
3240
3241 ! Locals:
3242 type(struct_gsv) :: statevector_vertCorr
3243 type(struct_gsv) :: statevector_oneMember
3244 real(8), pointer :: ptr3d_r8(:,:,:)
3245 real(8), pointer :: ptr3d_r8_oneMember(:,:,:)
3246 real(8) :: dnEns
3247 integer :: lonIndex, latIndex, varLevIndex1, varLevIndex2, memberIndex
3248 integer :: levIndex1, dateStamp
3249 character(len=4), pointer :: varNamesList(:)
3250 character(len=12) :: etiket
3251 character(len=4) :: varName
3252 character(len=3) :: levIndexStr
3253 character(len=256) :: ensFileName, outFileName
3254
3255 ! Set the dateStamp using the first ensemble member
3256 call fln_ensfileName(ensFileName, ens_getPathName(ensPerts), memberIndex_opt=1)
3257 dateStamp = tim_getDatestampFromFile(ensFileName)
3258
3259 ! Get list of variable names in the ensemble
3260 nullify(varNamesList)
3261 call ens_varNamesList(varNamesList,ensPerts)
3262
3263 call gsv_allocate(statevector_vertCorr, ens_getNumStep(ensPerts), &
3264 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3265 datestamp_opt=dateStamp, mpi_local_opt=.true., &
3266 mpi_distribution_opt='Tiles', dataKind_opt=8 )
3267
3268 call gsv_allocate(statevector_oneMember, ens_getNumStep(ensPerts), &
3269 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3270 datestamp_opt=dateStamp, mpi_local_opt=.true., &
3271 mpi_distribution_opt='Tiles', dataKind_opt=8 )
3272
3273 dnEns = 1.0d0/dble(nEns-1)
3274
3275 call gsv_getField(statevector_vertCorr,ptr3d_r8)
3276
3277 ! Loop over all vertical levels and variables
3278 varLev1: do varLevIndex1 = statevector_vertCorr%mykBeg, statevector_vertCorr%mykEnd
3279
3280 varName = ens_getVarNameFromK(ensPerts,varLevIndex1)
3281 levIndex1 = ens_getLevFromK(ensPerts,varLevIndex1)
3282
3283 ! Compute vertical correlations relative to varLev1
3284 call gsv_zero(statevector_vertCorr)
3285 member: do memberIndex = 1, nEns
3286 call ens_copyMember(ensPerts, statevector_oneMember, memberIndex)
3287 call gsv_getField(statevector_oneMember,ptr3d_r8_oneMember)
3288 varLev2: do varLevIndex2 = statevector_vertCorr%mykBeg, statevector_vertCorr%mykEnd
3289 do latIndex = statevector_vertCorr%myLatBeg, statevector_vertCorr%myLatEnd
3290 do lonIndex = statevector_vertCorr%myLonBeg, statevector_vertCorr%myLonEnd
3291 ptr3d_r8(lonIndex,latIndex,varLevIndex2) = &
3292 ptr3d_r8(lonIndex,latIndex,varLevIndex2) + &
3293 ptr3d_r8_oneMember(lonIndex,latIndex,varLevIndex1)* &
3294 ptr3d_r8_oneMember(lonIndex,latIndex,varLevIndex2)
3295 end do
3296 end do
3297 end do varLev2
3298 end do member
3299 call gsv_scale(statevector_vertCorr,dnEns)
3300
3301 ! Write to file the correlation matrix 'row' for this value of varLev1
3302 write(levIndexStr,'(i3.3)') levIndex1
3303 etiket = 'VCOR_' // trim(varName) // levIndexStr
3304 outFileName = './vertCorr_' // trim(varName) // levIndexStr // '.fst'
3305 call gio_writeToFile(statevector_vertCorr, &
3306 trim(outFileName), etiket_in = etiket, &
3307 typvar_opt = 'E', numBits_opt = 32)
3308 write(*,*) 'calcLocalVertCorrMatrix: finished variable/level =', varName, levIndex1
3309 end do varLev1
3310
3311 write(*,*) 'calcLocalVertCorrMatrix: finished computing the local vertical correlation matrix...'
3312
3313 call gsv_deallocate(statevector_vertCorr)
3314 call gsv_deallocate(statevector_oneMember)
3315
3316 end subroutine calcLocalVertCorrMatrix
3317
3318 !--------------------------------------------------------------------------
3319 ! calcVertModesSpec
3320 !--------------------------------------------------------------------------
3321 subroutine calcVertModesSpec(ensPerts, vModes)
3322 !
3323 !:Purpose: Compute the amplitude of the ensemble perturbations after their
3324 ! projection onto the vertical modes. This is the vertical equivalent
3325 ! of power spectra in the horizontal.
3326 !
3327 implicit none
3328
3329 ! Arguments:
3330 type(struct_ens), intent(inout) :: ensPerts
3331 type(struct_vms), intent(in) :: vModes
3332
3333 ! Locals:
3334 type(struct_gsv) :: gridStateVector_oneMember
3335 real(8), allocatable :: vertModesState3d(:,:,:)
3336 real(8), pointer :: gridState3d(:,:,:)
3337 real(8), allocatable :: powerSpec(:,:)
3338 real(8), allocatable :: latWeight(:) ! Weight given to grid point in the statistic computation
3339 real(8) :: sumWeight
3340 character(len=4), pointer :: varNamesList(:)
3341 character(len=128) :: outfilename
3342 integer :: numVar, nLev, varIndex, memberIndex
3343 integer :: nMode, modeIndex, latIndex, lonIndex
3344
3345 !
3346 !- Setup
3347 !
3348 nullify(varNamesList)
3349 call ens_varNamesList(varNamesList,ensPerts)
3350 numVar = size(varNamesList)
3351
3352 allocate(powerSpec(numVar,max(nLevEns_M,nLevEns_T)))
3353 powerSpec(:,:)=0.0d0
3354
3355 call gsv_allocate(gridStateVector_oneMember, ens_getNumStep(ensPerts), &
3356 ens_getHco(ensPerts), ens_getVco(ensPerts), varNames_opt=varNamesList, &
3357 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
3358 mpi_distribution_opt='Tiles', dataKind_opt=8 )
3359
3360 allocate(latWeight(hco_ens%nj))
3361 do latIndex = 1, hco_ens%nj
3362 latWeight(latIndex) = cos(hco_ens%lat(latIndex))
3363 if (mmpi_myid == 0) then
3364 write(*,*) latIndex, hco_ens%lat(latIndex), cos(hco_ens%lat(latIndex))
3365 end if
3366 end do
3367
3368 sumWeight = 0.d0
3369 do latIndex = myLatBeg, myLatEnd
3370 do lonIndex = myLonBeg, myLonEnd
3371 sumWeight = sumWeight + latWeight(latIndex)
3372 end do
3373 end do
3374
3375 !
3376 !- Vertical modes decomposition and coefficient summation
3377 !
3378 do memberIndex = 1, nEns
3379 write(*,*) 'calcVertModesSpec, member index = ', memberIndex
3380 call ens_copyMember(ensPerts, gridStateVector_oneMember, memberIndex)
3381
3382 do varIndex = 1, numVar
3383 write(*,*) 'calcVertModesSpec, varName = ', varIndex, varNamesList(varIndex)
3384 if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))).eq.'MM') then
3385 nLev = nLevEns_M
3386 nMode = nLev
3387 else if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))).eq.'TH') then
3388 nLev = nLevEns_T
3389 nMode = nLev
3390 else
3391 if (memberIndex == 1) then
3392 write(*,*) ' Skipping this variable not on momentum or thermodynamic levels'
3393 end if
3394 cycle
3395 end if
3396
3397 nullify(gridState3d)
3398 call gsv_getField(gridStateVector_oneMember,gridState3d,varName_opt=varNamesList(varIndex))
3399
3400 if (allocated(vertModesState3d)) deallocate(vertModesState3d)
3401 allocate(vertModesState3d(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nMode))
3402
3403 call vms_transform(vModes, vertModesState3d, gridState3d, &
3404 'GridPointToVertModes', myLonBeg, myLonEnd, &
3405 myLatBeg, myLatEnd, nLev, varNamesList(varIndex))
3406
3407 !$OMP PARALLEL DO PRIVATE (modeIndex,latIndex,lonIndex)
3408 do modeIndex = 1, nMode
3409 do latIndex = myLatBeg, myLatEnd
3410 do lonIndex = myLonBeg, myLonEnd
3411 powerSpec(varIndex,modeIndex) = powerSpec(varIndex,modeIndex) &
3412 + vertModesState3d(lonIndex,latIndex,modeIndex)**2 &
3413 * latWeight(latIndex)
3414 end do
3415 end do
3416 end do
3417 !$OMP END PARALLEL DO
3418
3419 end do
3420
3421 end do
3422
3423 deallocate(vertModesState3d)
3424 call gsv_deallocate(gridStateVector_oneMember)
3425
3426 !
3427 !- Communicate between all tasks
3428 !
3429 call mmpi_allreduce_sumR8_2d(powerSpec, "GRID")
3430 call mmpi_allreduce_sumreal8scalar(sumWeight, "GRID")
3431
3432 !
3433 !- Apply the appropriate scaling
3434 !
3435 powerSpec(:,:) = powerSpec(:,:)/(dble(nEns-1)*sumWeight)
3436
3437 !
3438 !- Write to file
3439 !
3440 if (mmpi_myid == 0) then
3441 do varIndex = 1, numVar
3442
3443 if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))) /= 'MM' .and. &
3444 vnl_varLevelFromVarName(trim(varNamesList(varIndex))) /= 'TH') cycle
3445
3446 outfilename = "./vertModesSpec_"//trim(varNamesList(varIndex))//".txt"
3447 write(*,*) 'calcVertModesSpec, writing ', trim(outfilename)
3448 open (unit=99,file=outfilename,action="write",status="new")
3449
3450 if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))) == 'MM') then
3451 nMode = nLevEns_M
3452 else
3453 nMode = nLevEns_T
3454 end if
3455
3456 do modeIndex = 1, nMode
3457 write(99,'(I4,2X,E10.4)') modeIndex, powerSpec(varIndex,modeIndex)
3458 end do
3459 close(unit=99)
3460 end do
3461 end if
3462
3463 end subroutine calcVertModesSpec
3464
3465 !--------------------------------------------------------------------------
3466 ! vertModesFilter
3467 !--------------------------------------------------------------------------
3468 subroutine vertModesFilter(ensPerts_in, ensPerts_out, vModes, waveBandIndex)
3469 !
3470 !:Purpose: Filter the ensemble perturbations to isolate a given vertical
3471 ! waveband of vertical modes.
3472 !
3473 implicit none
3474
3475 ! Arguments:
3476 type(struct_ens), intent(inout) :: ensPerts_in
3477 type(struct_ens), intent(inout) :: ensPerts_out
3478 type(struct_vms), intent(in) :: vModes
3479 integer, intent(in) :: waveBandIndex
3480
3481 ! Locals:
3482 type(struct_gsv) :: gridStateVector_oneMember
3483 real(8), allocatable :: vertModesState3d(:,:,:)
3484 real(8), pointer :: gridState3d(:,:,:)
3485 real(8), allocatable :: ResponseFunction(:)
3486 character(len=128) :: outfilename
3487 character(len=2) :: wbnum
3488 character(len=4), pointer :: varNamesList(:)
3489 integer :: nMode, nModeMax, modeIndex, latIndex, lonIndex
3490 integer :: numVar, nLev, varIndex, memberIndex
3491
3492 !
3493 !- 1. Pre-compute the response function
3494 !
3495 nModeMax=max(nLevEns_M,nLevEns_T)
3496 allocate(ResponseFunction(nModeMax))
3497 write(wbnum,'(I2.2)') waveBandIndex
3498 outfilename = "./vertResponseFunction_"//wbnum//".txt"
3499 if (mmpi_myid == 0) then
3500 open (unit=99,file=outfilename,action="write",status="new")
3501 end if
3502 do modeIndex = 1, nModeMax
3503 ResponseFunction(modeIndex) = spf_filterResponseFunction(dble(modeIndex), waveBandIndex, &
3504 vertWaveBandPeaks, nVertWaveBand)
3505 write(* ,'(I4,2X,F5.3)') modeIndex, ResponseFunction(modeIndex)
3506 if (mmpi_myid == 0) then
3507 write(99,'(I4,2X,F5.3)') modeIndex, ResponseFunction(modeIndex)
3508 end if
3509 end do
3510 if (mmpi_myid == 0) then
3511 close(unit=99)
3512 end if
3513
3514 !
3515 !- 2. Apply Filter
3516 !
3517 nullify(varNamesList)
3518 call ens_varNamesList(varNamesList,ensPerts_in)
3519 numVar = size(varNamesList)
3520
3521 call gsv_allocate(gridStateVector_oneMember, ens_getNumStep(ensPerts_in), &
3522 ens_getHco(ensPerts_in), ens_getVco(ensPerts_in), varNames_opt=varNamesList, &
3523 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
3524 mpi_distribution_opt='Tiles', dataKind_opt=8 )
3525
3526 do memberIndex = 1, nEns
3527 call ens_copyMember(ensPerts_in, gridStateVector_oneMember, memberIndex)
3528
3529 do varIndex = 1, numVar
3530
3531 nullify(gridState3d)
3532 call gsv_getField(gridStateVector_oneMember,gridState3d,varName_opt=varNamesList(varIndex))
3533
3534 if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))).eq.'MM') then
3535 nLev = nLevEns_M
3536 nMode = nLev
3537 else if (vnl_varLevelFromVarName(trim(varNamesList(varIndex))).eq.'TH') then
3538 nLev = nLevEns_T
3539 nMode = nLev
3540 else
3541 gridState3d(:,:,:) = 0.d0
3542 cycle
3543 end if
3544
3545 if (allocated(vertModesState3d)) deallocate(vertModesState3d)
3546 allocate(vertModesState3d(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nMode))
3547
3548 !- GridPoint space -> Vertical modes space
3549 call vms_transform(vModes, vertModesState3d, gridState3d, &
3550 'GridPointToVertModes', myLonBeg, myLonEnd, &
3551 myLatBeg, myLatEnd, nLev, varNamesList(varIndex))
3552
3553 !- Filtering
3554 !$OMP PARALLEL DO PRIVATE (modeIndex,latIndex,lonIndex)
3555 do modeIndex = 1, nMode
3556 do latIndex = myLatBeg, myLatEnd
3557 do lonIndex = myLonBeg, myLonEnd
3558 vertModesState3d(lonIndex,latIndex,modeIndex) = ResponseFunction(modeIndex) * &
3559 vertModesState3d(lonIndex,latIndex,modeIndex)
3560 end do
3561 end do
3562 end do
3563 !$OMP END PARALLEL DO
3564
3565 !- Vertical modes space -> GridPoint space
3566 call vms_transform(vModes, vertModesState3d, gridState3d, &
3567 'VertModesToGridPoint', myLonBeg, myLonEnd, &
3568 myLatBeg, myLatEnd, nLev, varNamesList(varIndex))
3569
3570 end do
3571
3572 call ens_insertMember(ensPerts_out, gridStateVector_oneMember, memberIndex)
3573
3574 end do
3575
3576 deallocate(vertModesState3d)
3577 call gsv_deallocate(gridStateVector_oneMember)
3578 deallocate(ResponseFunction)
3579
3580 end subroutine vertModesFilter
3581
3582end module calcStatsGlb_mod
