1MODULE bMatrixHI_mod
2 ! MODULE bMatrixHI_mod (prefix='bhi' category='2. B and R matrices')
3 !
4 !:Purpose: Performs transformation from control vector to analysis increment
5 ! (and adjoint transformation) using the background-error covariance
6 ! matrix based on homogeneous and isotropic correlations. This is
7 ! the Global version. A separate module exists for limited-area applications.
8 !
9 use midasMpi_mod
10 use earthConstants_mod
11 use gridStateVector_mod
12 use globalSpectralTransform_mod
13 use horizontalCoord_mod
14 use verticalCoord_mod
15 use timeCoord_mod
16 use varNameList_mod
17 use utilities_mod
18 use gridVariableTransforms_mod
19 use interpolation_mod
20 use calcHeightAndPressure_mod
21 implicit none
22 save
23 private
24
25 ! public procedures
26 public :: bhi_Setup,bhi_BSqrt,bhi_BSqrtAd,bhi_Finalize,bhi_expandToMPIglobal,bhi_expandToMPIglobal_r4,bhi_reduceToMPIlocal,bhi_reduceToMPIlocal_r4
27 public :: bhi_getScaleFactor,bhi_truncateCV
28
29
30 logical :: initialized = .false.
31 integer :: nj_l,ni_l
32 integer :: AnalGridID ! EZscintID
33 integer :: nlev_M,nlev_T,nlev_T_even,nkgdim,nkgdim2,nkgdimSqrt
34 integer :: nla_mpiglobal,nla_mpilocal
35 integer :: cvDim_mpilocal,cvDim_mpiglobal
36 integer :: gstID, gstID2
37 integer :: nlev_bdl
38 type(struct_vco),pointer :: vco_anl
39
40 real(8),allocatable :: tantheta(:,:)
41 real(8),allocatable :: PtoT(:,:,:)
42
43 real(8),pointer :: rgsig(:,:)
44 real(8),pointer :: rgsiguu(:,:),rgsigvv(:,:),rgsigtt(:,:),rgsigtb(:,:),rgsigq(:,:)
45 real(8),pointer :: rgsigps(:),rgsigpsb(:)
46 real(8),allocatable :: tgstdbg(:,:)
47
48 real(8),allocatable :: corns(:,:,:)
49 real(8),allocatable :: rstddev(:,:)
50
51 ! namelist variables:
52 integer :: ntrunc ! spectral trunction
53 real(8) :: scaleFactor(vco_maxNumLevels) ! scale factor applied to variances (all variables)
54 real(8) :: scaleFactorLQ(vco_maxNumLevels) ! scale factor applied to humidity
55 real(8) :: scaleFactorCC(vco_maxNumLevels) ! scale factor applied to velocity potential
56 logical :: scaleTG ! scale factor applied to skin temperature
57 logical :: TweakTG ! adjust skin temp variance based on land-sea mask and sea ice
58 logical :: ReadWrite_sqrt ! choose to read or write the sqrt of correlations
59 logical :: squareSqrt ! choose to use the 'square' formulation of corr matrix (not used)
60 character(len=4) :: stddevMode ! can be 'GD2D' or 'SP2D'
61 integer :: numModeZero ! number of eigenmodes to set to zero
62
63 ! constants
64 real(8) :: rcscltg(1)=100000.d0
65 real(8) :: rlimsuptg=3.0d0
66 logical :: llimtg=.true.
67 integer :: nulbgst=0
68 integer :: nLevPtoT
69 real(8) :: rvlocbalt = 6.0d0
70 real(8) :: rvlocpsi = 6.0d0
71 real(8) :: rvlocchi = 6.0d0
72 real(8) :: rvlocpsitt = 6.0d0
73 real(8) :: rvlocunbalt = 4.0d0
74 real(8) :: rvloclq = 4.0d0
75 real(8) :: rlimlv_bdl = 85000.0d0
76
77 ! this should come from state vector object
78 integer :: numvar3d
79 integer :: numvar2d
80 integer :: nspositVO
81 integer :: nspositDI
82 integer :: nspositTT
83 integer :: nspositQ
84 integer :: nspositPS
85 integer :: nspositTG
86
87 real(8), pointer :: pressureProfile_M(:),pressureProfile_T(:)
88
89 integer :: mymBeg,mymEnd,mymSkip,mymCount
90 integer :: mynBeg,mynEnd,mynSkip,mynCount
91 integer :: maxMyNla
92 integer :: myLatBeg,myLatEnd
93 integer :: myLonBeg,myLonEnd
94 integer, pointer :: ilaList_mpiglobal(:)
95 integer, pointer :: ilaList_mpilocal(:)
96
97 integer,external :: get_max_rss
98
99
100CONTAINS
101
102 SUBROUTINE BHI_setup(hco_in,vco_in,CVDIM_OUT, mode_opt)
103 implicit none
104
105 ! Arguments:
106 type(struct_hco), pointer, intent(in) :: hco_in
107 type(struct_vco), pointer, intent(in) :: vco_in
108 integer , intent(out) :: cvDim_out
109 character(len=*), optional, intent(in) :: mode_opt
110
111 ! Locals:
112 character(len=15) :: bhi_mode
113 integer :: jlev, nulnam, ierr, fnom, fclos, fstouv, fstfrm
114 integer :: jm, jn, latPerPE, lonPerPE, latPerPEmax, lonPerPEmax, Vcode_anl
115 logical :: llfound, lExists
116 real(8) :: zps
117 type(struct_vco),pointer :: vco_file => null()
118 character(len=8) :: bFileName = './bgcov'
119
120 NAMELIST /NAMBHI/ntrunc,scaleFactor,scaleFactorLQ,scaleFactorCC,scaleTG,numModeZero,squareSqrt,TweakTG,ReadWrite_sqrt,stddevMode
121
122 if(mmpi_myid == 0) write(*,*) 'bhi_setup: starting'
123 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
124
125 if ( present(mode_opt) ) then
126 if ( trim(mode_opt) == 'Analysis' .or. trim(mode_opt) == 'BackgroundCheck') then
127 bhi_mode = trim(mode_opt)
128 if(mmpi_myid == 0) write(*,*)
129 if(mmpi_myid == 0) write(*,*) 'bmatrixHI: Mode activated = ', trim(bhi_mode)
130 else
131 write(*,*)
132 write(*,*) 'mode = ', trim(mode_opt)
133 call utl_abort('bmatrixHI: unknown mode')
134 end if
135 else
136 bhi_mode = 'Analysis'
137 if(mmpi_myid == 0) write(*,*)
138 if(mmpi_myid == 0) write(*,*) 'bmatrixHI: Analysis mode activated (by default)'
139 end if
140
141 ! default values for namelist variables
142 ntrunc = 108
143 scaleFactor(:) = 0.0d0
144 scaleFactorLQ(:) = 1.0d0
145 scaleFactorCC(:) = 1.0d0
146 scaleTG = .true.
147 numModeZero = 0
148 squareSqrt = .false.
149 TweakTG = .false.
150 ReadWrite_sqrt = .false.
151 stddevMode = 'SP2D'
152
153 nulnam = 0
154 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
155 read(nulnam,nml=nambhi,iostat=ierr)
156 if ( ierr /= 0 ) call utl_abort( 'bhi_setup: Error reading namelist' )
157 if ( mmpi_myid == 0 ) write( *, nml = nambhi )
158 ierr = fclos( nulnam )
159
160 do jlev = 1, vco_maxNumLevels
161 if( scaleFactor(jlev) > 0.0d0 ) then
162 scaleFactor(jlev) = sqrt(scaleFactor(jlev))
163 else
164 scaleFactor(jlev) = 0.0d0
165 endif
166 enddo
167
168 if ( sum(scaleFactor(1:vco_maxNumLevels)) == 0.0d0 ) then
169 if ( mmpi_myid == 0 ) write(*,*) 'bmatrixHI: scaleFactor=0, skipping rest of setup'
170 cvdim_out = 0
171 return
172 end if
173
174 vco_anl => vco_in
175 nLev_M = vco_anl%nlev_M
176 nLev_T = vco_anl%nlev_T
177 ! need an even number of levels for spectral transform (gstID2)
178 if(mod(nLev_T,2).ne.0) then
179 nLev_T_even = nLev_T+1
180 else
181 nLev_T_even = nLev_T
182 endif
183 if(mmpi_myid == 0) write(*,*) 'BHI_setup: nLev_M, nLev_T, nLev_T_even=',nLev_M, nLev_T, nLev_T_even
184
185 ! check if analysisgrid and covariance file have the same vertical levels
186 call vco_SetupFromFile( vco_file, & ! OUT
187 bFileName ) ! IN
188 if (.not. vco_equal(vco_anl,vco_file)) then
189 call utl_abort('bmatrixHI: vco from analysisgrid and cov file do not match')
190 end if
191
192 Vcode_anl = vco_anl%Vcode
193 if(Vcode_anl .ne. 5002 .and. Vcode_anl .ne. 5005) then
194 write(*,*) 'Vcode_anl = ',Vcode_anl
195 call utl_abort('bmatrixHI: unknown vertical coordinate type!')
196 endif
197
198 if (.not. (gsv_varExist(varName='TT') .and. &
199 gsv_varExist(varName='UU') .and. &
200 gsv_varExist(varName='VV') .and. &
201 (gsv_varExist(varName='HU').or.gsv_varExist(varName='LQ')) .and. &
202 gsv_varExist(varName='P0')) ) then
203 call utl_abort('bmatrixHI: Some or all weather fields are missing. If it is desired to deactivate the weather assimilation, then all entries of the array SCALEFACTOR in the namelist NAMBHI should be set to zero.')
204 end if
205 if (.not. gsv_varExist(varName='TG')) then
206 write(*,*) 'bmatrixHI: WARNING: The TG field is missing. This must be present when assimilating'
207 write(*,*) ' radiance observations'
208 end if
209
210 do jlev = 1, max(nLev_M,nLev_T)
211 if(scaleFactorLQ(jlev).gt.0.0d0) then
212 scaleFactorLQ(jlev) = sqrt(scaleFactorLQ(jlev))
213 else
214 scaleFactorLQ(jlev) = 0.0d0
215 endif
216 enddo
217
218 do jlev = 1, max(nLev_M,nLev_T)
219 if(scaleFactorCC(jlev).gt.0.0d0) then
220 scaleFactorCC(jlev) = sqrt(scaleFactorCC(jlev))
221 else
222 scaleFactorCC(jlev) = 0.0d0
223 endif
224 enddo
225
226 if ( trim(bhi_mode) == 'BackgroundCheck' ) then
227 cvDim_out = 9999 ! Dummy value > 0 to indicate to the background check (s/r ose_compute_HBHT_ensemble)
228 ! that Bhi is used
229 return
230 end if
231
232 numvar3d = 4
233 numvar2d = 2
234
235 nLevPtot = nLev_M-1 ! ignore streamfunction at hyb=1, since highly correlated with next level
236 nspositVO = 1
237 nspositDI = 1*nLev_M+1
238 nspositTT = 2*nLev_M+1
239 nspositQ = 2*nLev_M+1*nLev_T+1
240 nspositPS = 2*nLev_M+2*nLev_T+1
241 nspositTG = 2*nLev_M+2*nLev_T+2
242 nkgdim = nLev_M*2 + nLev_T*2 + numvar2d
243 nkgdim2 = nkgdim + nLev_T
244 if(squareSqrt) then
245 nkgdimSqrt = nkgdim2
246 else
247 nkgdimSqrt = nkgdim
248 endif
249 nla_mpiglobal = (ntrunc+1)*(ntrunc+2)/2
250
251 ni_l = hco_in%ni
252 nj_l = hco_in%nj
253 AnalGridID = hco_in%EZscintID
254
255 gstID = gst_setup(ni_l,nj_l,ntrunc,nkgdim)
256 gstID2 = gst_setup(ni_l,nj_l,ntrunc,nlev_T_even)
257 if(mmpi_myid == 0) write(*,*) 'BHI:returned value of gstID =',gstID
258 if(mmpi_myid == 0) write(*,*) 'BHI:returned value of gstID2=',gstID2
259
260 call mmpi_setup_latbands(nj_l, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
261 call mmpi_setup_lonbands(ni_l, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
262
263 call mmpi_setup_m(ntrunc,mymBeg,mymEnd,mymSkip,mymCount)
264 call mmpi_setup_n(ntrunc,mynBeg,mynEnd,mynSkip,mynCount)
265
266 call gst_ilaList_mpiglobal(ilaList_mpiglobal,nla_mpilocal,maxMyNla,gstID,mymBeg,mymEnd,mymSkip,mynBeg,mynEnd,mynSkip)
267 call gst_ilaList_mpilocal(ilaList_mpilocal,gstID,mymBeg,mymEnd,mymSkip,mynBeg,mynEnd,mynSkip)
268
269 ! compute mpilocal control vector size
270 cvDim_mpilocal = 0
271 do jm = mymBeg, mymEnd, mymSkip
272 do jn = mynBeg, mynEnd, mynSkip
273 if(jm.le.jn) then
274 if(jm == 0) then
275 ! only real component for jm=0
276 cvDim_mpilocal = cvDim_mpilocal + 1*nkgdimSqrt
277 else
278 ! both real and imaginary components for jm>0
279 cvDim_mpilocal = cvDim_mpilocal + 2*nkgdimSqrt
280 endif
281 endif
282 enddo
283 enddo
284 cvDim_out = cvDim_mpilocal
285
286 ! also compute mpiglobal control vector dimension
287 call rpn_comm_allreduce(cvDim_mpilocal,cvDim_mpiglobal,1,"mpi_integer","mpi_sum","GRID",ierr)
288
289 allocate(PtoT(nlev_T+1,nlev_M,nj_l))
290 allocate(tantheta(nlev_M,nj_l))
291 allocate(rgsig(nj_l,nkgdim))
292 allocate(tgstdbg(ni_l,nj_l))
293 rgsiguu => rgsig(1:nj_l,nspositVO:nspositVO+nlev_M-1)
294 rgsigvv => rgsig(1:nj_l,nspositDI:nspositDI+nlev_M-1)
295 rgsigtt => rgsig(1:nj_l,nspositTT:nspositTT+nlev_T-1)
296 rgsigq => rgsig(1:nj_l,nspositQ :nspositQ +nlev_T-1)
297 rgsigps => rgsig(1:nj_l,nspositPS)
298 allocate(rgsigtb(nj_l,nlev_T))
299 allocate(rgsigpsb(nj_l))
300 allocate(corns(nkgdim2,nkgdim2,0:ntrunc))
301 allocate(rstddev(nkgdim2,0:ntrunc))
302
303 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
304
305 zps = 101000.D0
306 call czp_fetch1DLevels(vco_anl, zps, &
307 profM_opt=pressureProfile_M, profT_opt=pressureProfile_T)
308
309 llfound = .false.
310 nlev_bdl = 0
311 do jlev = 1, nlev_M
312 if(.not.llfound .and. (pressureProfile_M(jlev) .ge. rlimlv_bdl )) then
313 nlev_bdl = jlev
314 llfound = .true.
315 endif
316 enddo
317
318 inquire(file=bFileName,exist=lExists)
319 IF ( lexists )then
320 ierr = fnom(nulbgst,bFileName,'RND+OLD+R/O',0)
321 if ( ierr == 0 ) then
322 ierr = fstouv(nulbgst,'RND+OLD')
323 else
324 call utl_abort('BHI_setup:NO BACKGROUND STAT FILE!!')
325 endif
326 endif
327
328 call BHI_rdspPtoT
329
330 call BHI_readcorns2
331 if(mmpi_myid == 0) write(*,*) 'Memory Used (after readcorns2): ',get_max_rss()/1024,'Mb'
332
333 call BHI_sutg
334
335 if(stddevmode == 'GD2D') then
336 call BHI_rdstd
337 elseif(stddevMode == 'SP2D') then
338 call BHI_rdspstd_newfmt
339 else
340 call utl_abort('BHI_setup: unknown stddevMode')
341 endif
342
343 call BHI_scalestd
344
345 call BHI_sucorns2
346 if(mmpi_myid == 0) write(*,*) 'Memory Used (after sucorns2): ',get_max_rss()/1024,'Mb'
347
348 ierr = fstfrm(nulbgst)
349 ierr = fclos(nulbgst)
350
351 if(mmpi_myid == 0) write(*,*) 'END OF BHI_SETUP'
352
353 initialized = .true.
354
355 END SUBROUTINE BHI_setup
356
357
358 subroutine bhi_getScaleFactor(scaleFactor_out)
359 implicit none
360
361 ! Arguments:
362 real(8), intent(out) :: scaleFactor_out(:)
363
364 ! Locals:
365 integer :: jlev
366
367 do jlev = 1, max(nLev_M,nLev_T)
368 scaleFactor_out(jlev) = scaleFactor(jlev)
369 enddo
370
371 end subroutine bhi_getScaleFactor
372
373
374 SUBROUTINE BHI_scalestd
375 implicit none
376
377 ! Locals:
378 integer :: jlev, jlon, jlat, shift_level, Vcode_anl
379
380 Vcode_anl = vco_anl%Vcode
381 if(Vcode_anl == 5002) then
382 shift_level = 1
383 else
384 shift_level = 0
385 endif
386
387 do jlev = 1, nlev_M
388 do jlat = 1, nj_l
389 rgsiguu(jlat,jlev) = scaleFactor(jlev+shift_level)*rgsiguu(jlat,jlev)
390 rgsigvv(jlat,jlev) = scaleFactorCC(jlev+shift_level)*scaleFactor(jlev+shift_level)*rgsigvv(jlat,jlev)
391 enddo
392 enddo
393 do jlev = 1, nlev_T
394 do jlat = 1, nj_l
395 rgsigtt(jlat,jlev) = scaleFactor(jlev)*rgsigtt(jlat,jlev)
396 rgsigq(jlat,jlev) = scaleFactorLQ(jlev)*scaleFactor(jlev)*rgsigq(jlat,jlev)
397 rgsigtb(jlat,jlev) = scaleFactor(jlev)*rgsigtb(jlat,jlev)
398 enddo
399 enddo
400 do jlat = 1, nj_l
401 rgsigpsb(jlat) = scaleFactor(max(nLev_M,nLev_T))*rgsigpsb(jlat)
402 rgsigps(jlat) = scaleFactor(max(nLev_M,nLev_T))*rgsigps(jlat)
403 enddo
404 ! User has the option to not scale down the STDDEV of TG (because underestimated in Benkf)
405 if(scaleTG) then
406 do jlat = 1, nj_l
407 do jlon = 1, ni_l
408 tgstdbg(jlon,jlat) = scaleFactor(max(nLev_M,nLev_T))*tgstdbg(jlon,jlat)
409 enddo
410 enddo
411 endif
412
413 END SUBROUTINE BHI_scalestd
414
415
416 SUBROUTINE BHI_SUCORNS2
417 implicit none
418
419 ! Locals:
420 real(8) :: eigenval(nkgdim2), eigenvec(nkgdim2,nkgdim2), result(nkgdim2,nkgdim2)
421 real(8) :: eigenvalsqrt(nkgdim2), eigenvec2(nkgdim2,nkgdim2), eigenvalsqrt2(nkgdim2)
422 integer :: jlat,jn,jk1,jk2,jk3
423 integer :: ilwork,info,klatPtoT
424 integer :: iulcorvert, ikey, nsize
425 real(8) :: zwork(2*4*nkgdim2)
426 real(8) :: ztt(nlev_T,nlev_T,(ntrunc+1)),ztpsi(nlev_T,nlev_M,(ntrunc+1))
427 real(8) :: ztlen,zcorr,zr,zpres1,zpres2
428 real(8) :: zfact,zfact2,zcoriolis,zpsips(nLevPtoT)
429 real(8) :: zpsi(nlev_M,nlev_M),zfacttb(nj_l,nlev_T),zfactpsb(nj_l)
430 real(8) :: corvert(nkgdim2,nkgdim2)
431 real(8),allocatable :: corns_temp(:,:,:)
432 logical :: lldebug, lfound_sqrt
433 ! standard file variables
434 integer :: ini,inj,ink, inpas, inbits, idatyp, ideet
435 integer :: ip1,ip2,ip3,ig1,ig2,ig3,ig4,iswa,ilength,idltf
436 integer :: iubc,iextr1,iextr2,iextr3,ierr
437 integer :: idateo
438 character(len=2) :: cltypvar
439 character(len=1) :: clgrtyp
440 character(len=4) :: clnomvar
441 character(len=12) :: cletiket
442 integer :: fstprm,fstinf
443 integer :: fnom,fstouv,fstfrm,fclos
444
445 lldebug = .false.
446
447 iulcorvert = 0
448 if(mmpi_myid==0) then
449 ierr = fnom(iulcorvert,'corvert_modular.fst','RND',0)
450 ierr = fstouv(iulcorvert,'RND')
451 endif
452
453 klatPtoT = 1
454 zfactpsb(:) = 0.0d0
455 zfacttb(:,:) = 0.0d0
456
457 if(lldebug) then
458 do jk1 = 1, nlev_T
459 do jk2 = 1, nlevPtoT
460 write(622,*) jk1,jk2,klatPtoT,PtoT(jk1,jk2,klatPtoT)
461 enddo
462 enddo
463 endif
464
465 ! explicitly compute the balanced temperature and temperature-psi correlations
466
467 do jn = 0, ntrunc
468
469 ztpsi(:,:,jn+1) = 0.0d0
470 ztt(:,:,jn+1) = 0.0d0
471 do jk1 = 1, nlevPtoT
472 do jk2 = 1, nlev_T
473 do jk3 = 1, nlevPtoT
474 ztpsi(jk2,jk1,jn+1) = ztpsi(jk2,jk1,jn+1)+PtoT(jk2,jk3,klatPtoT)*corns(jk3,jk1,jn)
475 enddo
476 enddo
477 enddo
478 if(nlevPtoT.lt.nlev_M) then
479 do jk1 = (nlevPtoT+1), nlev_M
480 do jk2 = 1, nlev_T
481 ztpsi(jk2,jk1,jn+1) = ztpsi(jk2,nlevPtoT,jn+1)
482 enddo
483 enddo
484 endif
485 do jk1 = 1, nlev_T
486 do jk2 = 1, nlev_T
487 do jk3 = 1, nlevPtoT
488 ztt(jk2,jk1,jn+1) = ztt(jk2,jk1,jn+1)+ztpsi(jk2,jk3,jn+1)*PtoT(jk1,jk3,klatPtoT)
489 enddo
490 enddo
491 enddo
492 enddo
493
494 if(lldebug) then
495 write(620,*) ztt
496 write(621,*) ztpsi
497 endif
498
499 ! fill in blocks for balance temperature
500
501 do jn = 0, ntrunc
502 do jk1 = 1, nlev_T
503 do jk2 = 1, nlev_T
504 corns(nkgdim+jk2,nkgdim+jk1,jn) = ztt(jk2,jk1,jn+1)
505 enddo
506 enddo
507 do jk1 = 1, nlev_M
508 do jk2 = 1, nlev_T
509 corns( jk1,nkgdim+jk2,jn) = ztpsi(jk2,jk1,jn+1)
510 corns(nkgdim+jk2, jk1,jn) = ztpsi(jk2,jk1,jn+1)
511 enddo
512 enddo
513 enddo
514
515 ! Save un-localized PSI correlations
516 do jk2 = 1, nlev_M
517 do jk1 = 1, nlev_M
518 zpsi(jk1,jk2) = 0.0d0
519 do jn = 0, ntrunc
520 zpsi(jk1,jk2) = zpsi(jk1,jk2)+((2*jn+1)*corns(jk1,jk2,jn))
521 enddo
522 enddo
523 enddo
524
525 ! Apply vertical localization to corrns
526
527 ! unbalanced temperature
528 ztlen = rvlocunbalt
529 if(ztlen.gt.0.0d0) then
530 ! calculate 5'th order function (from Gaspari and Cohn)
531 do jk1 = 1, nlev_T
532 zpres1 = log(pressureProfile_T(jk1))
533 do jk2 = 1, nlev_T
534 zpres2 = log(pressureProfile_T(jk2))
535 zr = abs(zpres2 - zpres1)
536 zcorr = gasparicohn(ztlen,zr)
537 do jn = 0, ntrunc
538 corns(jk1+2*nlev_M,jk2+2*nlev_M,jn) = &
539 corns(jk1+2*nlev_M,jk2+2*nlev_M,jn)*zcorr
540 enddo
541 enddo
542 enddo
543 endif
544
545 ! balanced temperature
546 ztlen = rvlocbalt
547 if(ztlen.gt.0.0d0) then
548 ! calculate 5'th order function (from Gaspari and Cohn)
549 do jk1 = 1, nlev_T
550 zpres1 = log(pressureProfile_T(jk1))
551 do jk2 = 1, nlev_T
552 zpres2 = log(pressureProfile_T(jk2))
553 zr = abs(zpres2 - zpres1)
554 zcorr = gasparicohn(ztlen,zr)
555 do jn = 0, ntrunc
556 corns(jk1+nkgdim,jk2+nkgdim,jn) = &
557 corns(jk1+nkgdim,jk2+nkgdim,jn)*zcorr
558 enddo
559 enddo
560 enddo
561 endif
562
563 ! streamfunction
564 ztlen = rvlocpsi ! specify length scale (in units of ln(Pressure))
565 if(ztlen.gt.0.0d0) then
566 ! calculate 5'th order function (from Gaspari and Cohn)
567 do jk1 = 1, nlev_M
568 zpres1 = log(pressureProfile_M(jk1))
569 do jk2 = 1, nlev_M
570 zpres2 = log(pressureProfile_M(jk2))
571 zr = abs(zpres2 - zpres1)
572 zcorr = gasparicohn(ztlen,zr)
573 do jn = 0, ntrunc
574 corns(jk1,jk2,jn) = corns(jk1,jk2,jn)*zcorr
575 enddo
576 enddo
577 enddo
578 endif
579
580 ! temp-psi cross-correlations
581 ztlen = rvlocpsitt ! specify length scale (in units of ln(Pressure))
582 if(ztlen.gt.0.0d0) then
583 ! calculate 5'th order function (from Gaspari and Cohn)
584 do jk1 = 1, nlev_M
585 zpres1 = log(pressureProfile_M(jk1))
586 do jk2 = 1, nlev_T
587 zpres2 = log(pressureProfile_T(jk2))
588 zr = abs(zpres2 - zpres1)
589 zcorr = gasparicohn(ztlen,zr)
590 do jn = 0, ntrunc
591 corns(jk1,jk2+nkgdim,jn) = corns(jk1,jk2+nkgdim,jn)*zcorr
592 corns(jk2+nkgdim,jk1,jn) = corns(jk2+nkgdim,jk1,jn)*zcorr
593 enddo
594 enddo
595 enddo
596 endif
597
598 ! velocity potential (unbalanced)
599 ztlen = rvlocchi ! specify length scale (in units of ln(Pressure))
600 if(ztlen.gt.0.0d0) then
601 ! calculate 5'th order function (from Gaspari and Cohn)
602 do jk1 = 1, nlev_M
603 zpres1 = log(pressureProfile_M(jk1))
604 do jk2 = 1, nlev_M
605 zpres2 = log(pressureProfile_M(jk2))
606 zr = abs(zpres2 - zpres1)
607 zcorr = gasparicohn(ztlen,zr)
608 do jn = 0, ntrunc
609 corns(jk1+nlev_M,jk2+nlev_M,jn) = corns(jk1+nlev_M,jk2+nlev_M,jn)*zcorr
610 enddo
611 enddo
612 enddo
613 endif
614
615 ! cross-correlation t'-ps'
616 if(.true.) then
617 ztlen = rvlocunbalt ! specify length scale (in units of ln(Pressure))
618 if(ztlen.gt.0.0d0) then
619 ! calculate 5'th order function (from Gaspari and Cohn)
620 zpres1 = log(pressureProfile_T(nlev_T))
621 do jk2 = 1, nlev_T
622 zpres2 = log(pressureProfile_T(jk2))
623 zr = abs(zpres2 - zpres1)
624 zcorr = gasparicohn(ztlen,zr)
625 do jn = 0, ntrunc
626 corns(1+2*nlev_M+2*nlev_T,jk2+2*nlev_M,jn) = &
627 corns(1+2*nlev_M+2*nlev_T,jk2+2*nlev_M,jn)*zcorr
628 corns(jk2+2*nlev_M,1+2*nlev_M+2*nlev_T,jn) = &
629 corns(jk2+2*nlev_M,1+2*nlev_M+2*nlev_T,jn)*zcorr
630 enddo
631 enddo
632 endif
633 endif
634
635 ! humidity
636 ztlen = rvloclq ! specify length scale (in units of ln(Pressure))
637 if(ztlen.gt.0.0d0) then
638 ! calculate 5'th order function (from Gaspari and Cohn)
639 do jk1 = 1, nlev_T
640 zpres1 = log(pressureProfile_T(jk1))
641 do jk2 = 1, nlev_T
642 zpres2 = log(pressureProfile_T(jk2))
643 zr = abs(zpres2 - zpres1)
644 zcorr = gasparicohn(ztlen,zr)
645 do jn = 0, ntrunc
646 corns(jk1+2*nlev_M+nlev_T,jk2+2*nlev_M+nlev_T,jn) = &
647 corns(jk1+2*nlev_M+nlev_T,jk2+2*nlev_M+nlev_T,jn)*zcorr
648 enddo
649 enddo
650 enddo
651 endif
652
653 ! compute total vertical correlations (including for balanced temperature)
654 if(.true.) then
655 do jk2 = 1, nkgdim2
656 do jk1 = 1, nkgdim2
657 corvert(jk1,jk2) = 0.0d0
658 do jn = 0, ntrunc
659 corvert(jk1,jk2) = corvert(jk1,jk2)+((2*jn+1)*corns(jk1,jk2,jn))
660 enddo
661 enddo
662 enddo
663
664 if(lldebug) then
665 write(701,*) corvert
666 write(702,*) zpsi
667 endif
668
669 if(mmpi_myid == 0) then
670 ikey = fstinf(NULBGST,ini,inj,ink,-1,'CORRNS',-1,0,-1,' ','ZZ')
671 ierr = fstprm(ikey,idateo,ideet,inpas,ini,inj,ink, inbits &
672 ,idatyp,ip1,ip2,ip3,cltypvar,clnomvar,cletiket,clgrtyp &
673 ,ig1,ig2,ig3,ig4,iswa,ilength,idltf,iubc,iextr1,iextr2 &
674 ,iextr3)
675
676 ini = nkgdim2
677 inj = nkgdim2
678 ink = 1
679 ip1 = 0
680 ip2 = ntrunc
681 ip3 = 0
682 clnomvar = 'ZV'
683 cletiket = 'CORVERT'
684 idatyp = 5
685
686 ierr = utl_fstecr(corvert, -inbits, iulcorvert, idateo &
687 , ideet,inpas, ini, inj, ink, ip1, ip2, ip3, cltypvar, &
688 clnomvar,cletiket,clgrtyp,ig1, ig2, ig3, ig4, idatyp, &
689 .true.)
690
691 endif
692
693 ! Modify RGSIGTB to obtain correct sigma_Tb
694 do jk1 = 1, nlev_T
695 zfact = corvert(jk1+nkgdim,jk1+nkgdim)
696 do jlat = 1, nj_l
697 zcoriolis = abs(2.d0*ec_Omega*gst_getrmu(jlat,gstID))
698 if(zfact.gt.0.0d0.and.zcoriolis.ne.0.0d0) then
699 zfact2 = 1.0d0/(zfact*zcoriolis*zcoriolis)
700 else
701 zfact2 = 0.0d0
702 endif
703 zfacttb(jlat,jk1) = zfacttb(jlat,jk1)+zfact2
704 enddo
705 enddo
706
707 ! Modify RGSIGPSB to obtain correct sigma_PSb
708 do jlat = 1, nj_l
709 do jk2 = 1, nlevPtoT
710 zpsips(jk2) = 0.0d0
711 do jk1 = 1, nlevPtoT
712 zpsips(jk2) = zpsips(jk2)+PtoT(nlev_T+1,jk1,klatPtoT)*zpsi(jk1,jk2)
713 enddo
714 enddo
715 zfact = 0.0d0
716 do jk1 = 1, nlevPtoT
717 zfact = zfact+PtoT(nlev_T+1,jk1,klatPtoT)*zpsips(jk1)
718 enddo
719 zcoriolis = abs(2.d0*ec_Omega*gst_getrmu(jlat,gstID))
720 if(zfact.gt.0.0d0.and.zcoriolis.ne.0.0d0) then
721 zfact2 = 1.0d0/(zfact*zcoriolis*zcoriolis)
722 else
723 zfact2 = 0.0d0
724 endif
725 zfactpsb(jlat) = zfactpsb(jlat)+zfact2
726 enddo
727 endif
728
729 ! Modify RGSIGTB and RGSIGPSB to obtain correct sigma_Tb and sigma_Psb
730 do jlat = 1, nj_l
731 if(zfactpsb(jlat).gt.0.0d0) then
732 rgsigpsb(jlat) = rgsigpsb(jlat)*sqrt(zfactpsb(jlat))
733 else
734 rgsigpsb(jlat) = 0.0d0
735 endif
736 do jk1 = 1, nlev_T
737 if(zfacttb(jlat,jk1).gt.0.0d0) then
738 rgsigtb(jlat,jk1) = rgsigtb(jlat,jk1)*sqrt(zfacttb(jlat,jk1))
739 else
740 rgsigtb(jlat,jk1) = 0.0d0
741 endif
742 enddo
743 enddo
744
745 ! compute square-root of corns for each total wavenumber
746 allocate(corns_temp(nkgdim2,nkgdim2,0:ntrunc))
747 corns_temp(:,:,:)=0.0d0
748 do jn = mmpi_myid, ntrunc, mmpi_nprocs
749
750 do jk1 = 1, nkgdim2
751 do jk2 = 1, nkgdim2
752 eigenvec(jk2,jk1) = corns(jk2,jk1,jn)
753 enddo
754 enddo
755
756 ! CALCULATE EIGENVALUES AND EIGENVECTORS.
757 ilwork = 4*nkgdim2*2
758 call dsyev('V','U',nkgdim2,eigenvec,nkgdim2,eigenval,zwork,ilwork,info)
759 if(info.ne.0) then
760 write(*,*) 'bhi_sucorns2: non-zero value of info =',info,' returned by dsyev for wavenumber ',jn
761 call utl_abort('BHI_SUCORNS')
762 endif
763
764 ! set selected number of eigenmodes to zero
765 if(numModeZero.gt.0) then
766 write(*,*) 'bhi_sucorns2: setting ',numModeZero,' eigenvalues to zero for wavenumber n=',jn
767 write(*,*) 'bhi_sucorns2: original eigenvalues=',eigenval(:)
768 do jk1 = 1, numModeZero
769 eigenval(jk1) = 0.0d0
770 enddo
771 write(*,*) 'bhi_sucorns2: modified eigenvalues=',eigenval(:)
772 endif
773
774 do jk1 = 1, nkgdim2
775 if(eigenval(jk1).lt.1.0d-15) then
776 eigenvalsqrt(jk1) = 0.0d0
777 else
778 eigenvalsqrt(jk1) = sqrt(eigenval(jk1))
779 endif
780 enddo
781
782 ! Reverse the order of E-Values if old formulation (for compatibility)
783 if(.not.squareSqrt) then
784 do jk1 = 1, nkgdim2
785 eigenvalsqrt2(jk1) = eigenvalsqrt(nkgdim2-jk1+1)
786 do jk2 = 1, nkgdim2
787 eigenvec2(jk2,jk1) = eigenvec(jk2,nkgdim2-jk1+1)
788 enddo
789 enddo
790 eigenvalsqrt(:) = eigenvalsqrt2(:)
791 eigenvec(:,:) = eigenvec2(:,:)
792 endif
793
794 ! compute E * lambda^1/2
795 result(:,:) = 0.0d0
796 do jk1 = 1, nkgdimSqrt
797 do jk2 = 1, nkgdim2
798 result(jk2,jk1) = eigenvec(jk2,jk1)*eigenvalsqrt(jk1)
799 enddo
800 enddo
801
802 ! compute (E * lambda^1/2) * E^T if new formulation
803 if(squareSqrt) then
804 do jk1 = 1, nkgdim2
805 do jk2 = 1, nkgdim2
806 do jk3 = 1, nkgdim2
807 corns_temp(jk2,jk1,jn) = corns_temp(jk2,jk1,jn) + result(jk2,jk3)*eigenvec(jk1,jk3)
808 enddo
809 enddo
810 enddo
811 else
812 corns_temp(:,:,jn) = result(:,:)
813 endif
814
815 !if(jn == 30) then
816 ! write(200,*) corns(:,:,jn)
817 ! write(201,*) corns_temp(:,:,jn)
818 ! write(202,*) eigenval(:)
819 ! write(203,*) eigenvec(:,:)
820 ! call flush(200)
821 ! call flush(201)
822 ! call flush(202)
823 ! call flush(203)
824 !endif
825
826 enddo ! jn
827
828 nsize = nkgdim2*nkgdim2*(ntrunc+1)
829 call rpn_comm_allreduce(corns_temp,corns,nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
830 deallocate(corns_temp)
831
832 if(mmpi_myid==0) then
833 ierr = fstfrm(iulcorvert)
834 ierr = fclos(iulcorvert)
835 endif
836
837 if(ReadWrite_sqrt) then
838 ! if desired, read precomputed sqrt of corns which overwrites computed matrix
839 call readcorns_sqrt(lfound_sqrt)
840 if(.not.lfound_sqrt) then
841 ! if precomputed sqrt does not exist in stats, then write it out to a separate file
842 call writecorns_sqrt
843 endif
844 endif
845
846 END SUBROUTINE BHI_SUCORNS2
847
848
849 SUBROUTINE WRITECORNS_SQRT
850 implicit none
851
852 ! Locals:
853 integer :: jn, nulcorns_sqrt, ierr
854 ! standard file variables
855 integer :: ip1,ip2,ip3
856 integer :: idateo, ipak, idatyp
857 integer :: fnom, fstouv, fstfrm, fclos
858
859 write(*,*) 'WRITECORNS_SQRT: CORNS_SQRT will be written to file corns_sqrt.fst for NTRUNC =', ntrunc
860
861 if(mmpi_myid==0) then
862 nulcorns_sqrt = 0
863 ierr = fnom(nulcorns_sqrt,'corns_sqrt.fst','RND',0)
864 ierr = fstouv(nulcorns_sqrt,'RND')
865
866 ipak = -32
867 idatyp = 5
868 ip1 = 0
869 ip3 = ntrunc
870 idateo = 0
871
872 do jn = 0, ntrunc
873 ip2 = jn
874 ierr = utl_fstecr(corns(1,1,jn),ipak,nulcorns_sqrt,idateo,0,0,nkgdim2,nkgdim2,1, &
875 ip1,ip2,ip3,'X','ZZ','CORNS_SQRT','X',0,0,0,0,idatyp,.true.)
876 enddo
877
878 ierr = fstfrm(nulcorns_sqrt)
879 ierr = fclos(nulcorns_sqrt)
880 endif
881
882 END SUBROUTINE WRITECORNS_SQRT
883
884
885 SUBROUTINE READCORNS_SQRT(lfound_sqrt)
886 implicit none
887
888 ! Arguments:
889 logical, intent(out) :: lfound_sqrt
890
891 ! Locals:
892 integer :: jn, icornskey
893 integer :: jcol,jrow
894 real(8), allocatable :: zcornssrc(:,:)
895 ! standard file variables
896 integer :: ini,inj,ink
897 integer :: ip1,ip2,ip3
898 integer :: idateo
899 character(len=2) :: cltypvar
900 character(len=4) :: clnomvar
901 character(len=12) :: cletiket
902
903 allocate(zcornssrc(nkgdim2,nkgdim2))
904
905 write(*,*) 'READCORNS_SQRT: looking for CORNS_SQRT with NTRUNC =', ntrunc
906 do jn = 0, ntrunc
907
908 ! Looking for FST record parameters..
909 idateo = -1
910 cletiket = 'CORNS_SQRT'
911 ip1 = -1
912 ip2 = jn
913 ip3 = ntrunc
914 cltypvar = 'X'
915 clnomvar = 'ZZ'
916 icornskey = utl_fstlir(ZCORNSSRC,nulbgst,INI,INJ,INK,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
917
918 if( jn == 0 ) then
919 if(icornskey .lt.0 ) then
920 write(*,*) 'READCORNS_SQRT: CORNS_SQRT not found in stats file, use computed sqrt'
921 lfound_sqrt = .false.
922 return
923 else
924 write(*,*) 'READCORNS_SQRT: CORNS_SQRT found in stats file, will use it instead of computed sqrt'
925 lfound_sqrt = .true.
926 endif
927 endif
928
929 if (ini .ne. nkgdim2 .or. inj .ne. nkgdim2) then
930 call utl_abort('READCORNS2: BG stat levels inconsitencies')
931 endif
932
933 do jcol = 1, nkgdim2
934 do jrow = 1, nkgdim2
935 corns(jrow,jcol,jn) = zcornssrc(jrow,jcol)
936 enddo
937 enddo
938
939 enddo
940
941 deallocate(zcornssrc)
942
943 END SUBROUTINE READCORNS_SQRT
944
945
946 FUNCTION GASPARICOHN(ztlen,zr)
947 implicit none
948
949 ! Arguments:
950 real(8), intent(in) :: ztlen
951 real(8), intent(in) :: zr
952 ! Result:
953 real(8) :: gasparicohn
954
955 ! Locals:
956 real(8) :: zlc
957
958 zlc = ztlen/2.0d0
959 if(zr.le.zlc) then
960 gasparicohn = -0.250d0*(zr/zlc)**5+0.5d0*(zr/zlc)**4 &
961 +0.625d0*(zr/zlc)**3-(5.0d0/3.0d0)*(zr/zlc)**2+1.0d0
962 elseif(zr.le.(2.0d0*zlc)) then
963 gasparicohn = (1.0d0/12.0d0)*(zr/zlc)**5-0.5d0*(zr/zlc)**4 &
964 +0.625d0*(zr/zlc)**3+(5.0d0/3.0d0)*(zr/zlc)**2 &
965 -5.0d0*(zr/zlc)+4.0d0-(2.0d0/3.0d0)*(zlc/zr)
966 else
967 gasparicohn = 0.0d0
968 endif
969 if(gasparicohn.lt.0.0d0) gasparicohn = 0.0d0
970
971 END FUNCTION GASPARICOHN
972
973
974 SUBROUTINE BHI_CALCCORR(zgd,pcscl,klev)
975 implicit none
976
977 ! Arguments:
978 integer, intent(in) :: klev
979 real(8), intent(out) :: zgd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,klev)
980 real(8), intent(in) :: pcscl(klev)
981
982 ! Locals:
983 integer :: jlev, jlat, jlon
984 real(8) :: zr, dlfac, dltemp, dln, dlcsurn, dlc, dlcorr
985 ! parameters that define the correlation function
986 integer :: ntoar = 3
987 real(8) :: dlalpha = 0.2d0
988 integer :: kcorrtyp = 1
989
990 dlfac = 1.d0/(1.d0+dlalpha)
991 dln = 1.d0*real(ntoar,8)
992 dltemp = (3.d0*(1.d0 + dlalpha))/(1.d0 + dlalpha/(dln*dln))
993 dltemp = dsqrt(dltemp)
994
995 if (kcorrtyp == 1) then
996 ! Gaussian correlation
997 do jlev = 1, klev
998 dlc = 1.d0/dble(pcscl(jlev))
999 dlc = 0.5d0*dlc*dlc
1000 do jlat = myLatBeg, myLatEnd
1001 zr = ec_ra * acos(gst_getRmu(jlat,gstID))
1002 dlcorr = dexp(-(zr**2)*dlc)
1003 do jlon = myLonBeg, myLonEnd
1004 zgd(jlon,jlat,jlev) = dlcorr
1005 enddo
1006 enddo
1007 enddo
1008 elseif (kcorrtyp == 2) then
1009 ! Autoregressive (SOAR) correlation
1010 do jlev = 1, klev
1011 dlc = dltemp/dble(pcscl(jlev))
1012 dlcsurn = dlc/dln
1013 do jlat = myLatBeg, myLatEnd
1014 zr = ec_ra * acos(gst_getRmu(jlat,gstID))
1015 dlcorr = (1.d0 + dlc*zr + zr*dlc*zr*dlc/3.d0)*dexp(-zr*dlc) &
1016 + dlalpha*(1.d0 + dlcsurn*zr + zr*dlcsurn*zr*dlcsurn/3.d0)*dexp(-zr*dlcsurn)
1017 dlcorr = dlcorr*dlfac
1018 do jlon = myLonBeg, myLonEnd
1019 zgd(jlon,jlat,jlev) = dlcorr
1020 enddo
1021 enddo
1022 enddo
1023 else
1024 call utl_abort('CALCCORR- Undefined correlation type')
1025 endif
1026
1027 END SUBROUTINE BHI_calcCorr
1028
1029
1030 SUBROUTINE BHI_SUTG
1031 implicit none
1032
1033 ! Locals:
1034 logical :: llpb
1035 integer :: ikey, jlat, jlon, jla, ezgprm, ezqkdef
1036 integer :: jn, jm, ila_mpilocal, ila_mpiglobal, inlev, itggid, inmxlev, iset, nsize
1037 integer :: ezdefset
1038 integer :: ip1style,ip1kind
1039 integer :: koutmpg
1040 real(8), allocatable :: dltg(:,:), tgstdbg_tmp(:,:)
1041 real(8) :: cortgg(nla_mpiglobal,2)
1042 real(8) :: rcscltg_vec(nlev_T_even)
1043 real(8) :: zabs, zpole, dlfac
1044 real(8) :: zsp_mpilocal(nla_mpilocal,2,nlev_T_even)
1045 real(8) :: zgd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlev_T_even)
1046 real(8) :: zsp_mpiglobal(nla_mpiglobal,2,1)
1047 real(8),allocatable :: my_zsp_mpiglobal(:,:,:)
1048 real(4), allocatable :: TrialLandSeaMask(:,:), TrialSeaIceMask(:,:)
1049 real(4), allocatable :: AnalLandSeaMask(:,:), AnalSeaIceMask(:,:)
1050 ! standard file variables
1051 integer :: ini,inj,ink, idatyp
1052 integer :: ip1,ip2,ip3
1053 integer :: ierr,ntrials
1054 integer :: idateo
1055 integer :: fstprm,fstinf,iultg,fnom,fclos,fstouv,fstfrm
1056 integer :: ip1s(1), nulbgsts(1)
1057 integer :: TrlmNumberWanted
1058 integer :: fstlir, key, nultrl, ni_trial, nj_trial
1059 integer :: deet, npas, nbits, datyp
1060 integer :: ig1,ig2,ig3,ig4,swa, lng, dltf, ubc
1061 integer :: extra1, extra2, extra3
1062 integer :: TrialGridID
1063 character(len=2) :: cltypvar
1064 character(len=1) :: clgrtyp
1065 character(len=4) :: clnomvar
1066 character(len=12) :: cletiket
1067 character(len=2) :: flnum
1068 character(len=128) :: trialfile
1069 logical :: trialExists
1070
1071 !
1072 !- 1. Reading and processing TG standard deviations
1073 !
1074
1075 !- 1.1 Read the Std. Dev. field from ./bgcov file
1076 clnomvar = 'TG'
1077 idateo = -1
1078 inmxlev = 1
1079 ntrials = 1
1080 nulbgsts(1) = nulbgst
1081
1082 call utl_getfldprm(IP1S, IP2, IP3, INLEV, CLETIKET, CLTYPVAR, ITGGID, &
1083 clnomvar, idateo, inmxlev, nulbgsts, ip1style, ip1kind, &
1084 ntrials, koutmpg)
1085 ip1 = ip1s(1)
1086
1087 ierr = ezgprm(itggid,CLGRTYP,INI,INJ,IG1,IG2,IG3,IG4)
1088 allocate(dltg(ini,inj))
1089
1090 ikey = utl_fstlir(dltg,koutmpg,ini,inj,ink,idateo,cletiket,ip1, &
1091 ip2, ip3, cltypvar, clnomvar)
1092
1093 !- 1.2 Rearrange the data according to the analysis grid (if necessary)
1094 if(clgrtyp == 'G' .and. ni_l == ini .and. nj_l == inj .and. ig1 == 0 &
1095 .and. ig2 ==0 .and. ig3 == 0 .and.ig4 == 0) then
1096 !- 1.2.1 The std. dev. on the analysis grid
1097 do jlat = 1, nj_l
1098 do jlon = 1,ni_l
1099 tgstdbg(jlon,jlat) = dltg(jlon,jlat)
1100 enddo
1101 enddo
1102
1103 elseif(clgrtyp == 'G' .and. ni_l == ini .and. nj_l == inj .and. ig1 == &
1104 0 .and. ig2 ==1 .and. ig3 == 0 .and.ig4 == 0) then
1105 !- 1.2.2 flipped Gaussian grid no longer supported
1106 call utl_abort('bhi_sutg: The flipped Gaussian grid is no longer supported!')
1107
1108 else
1109 !- 1.2.3 The std. dev. are NOT on the analysis grid. Interpolation is needed
1110 iset = ezdefset(AnalGridID,itggid)
1111 if ( TweakTG ) then
1112 ierr = int_hInterpScalar(tgstdbg,dltg,interpDegree='NEAREST')
1113 else
1114 ierr = int_hInterpScalar(tgstdbg,dltg,interpDegree='CUBIC')
1115 end if
1116
1117 end if
1118
1119 !- 1.3 Tweaking the Std Dev.
1120
1121 !- 1.3.1 If specified at the top of the module, do not accept TG errors of more than value specified above
1122 if ( llimtg ) then
1123 if ( mmpi_myid == 0 ) write(*,*)
1124 if ( mmpi_myid == 0 ) write(*,*) 'Capping TG Std. Dev. using a max value (K) = ', rlimsuptg
1125 where ( tgstdbg > rlimsuptg) tgstdbg = rlimsuptg
1126 end if
1127
1128 !- 1.3.2 Take into account the Land-Sea mask and the Sea-Ice mask of the day
1129 if ( TweakTG ) then
1130
1131 if ( mmpi_myid == 0 ) write(*,*)
1132 if ( mmpi_myid == 0 ) write(*,*) 'Adjusting TG Std Dev based on LandSea and SeaIce masks'
1133
1134 !- Read MG and GL in the middle of the assimilation time window
1135 if ( tim_nStepObs == 1 ) then
1136 TrlmNumberWanted = 1
1137 else
1138 TrlmNumberWanted = nint( (tim_nStepObs + 1.d0) / 2.d0)
1139 end if
1140
1141 write(flnum,'(I2.2)') TrlmNumberWanted
1142 trialfile='./trlm_'//trim(flnum)
1143 inquire(file=trim(trialfile),exist=trialExists)
1144
1145 if ( .not. trialExists ) then
1146 if ( mmpi_myid == 0 ) write(*,*)
1147 if ( mmpi_myid == 0 ) write(*,*) 'Trial file not found = ', trialfile
1148 if ( mmpi_myid == 0 ) write(*,*) 'Look for an ensemble of trial files '
1149
1150 trialfile='./trlm_'//trim(flnum)//'_0001'
1151 inquire(file=trim(trialfile),exist=trialExists)
1152 if ( .not. trialExists ) then
1153 if ( mmpi_myid == 0 ) write(*,*) 'Ensemble trial file not found = ', trialfile
1154 call utl_abort('BMatrixHI : DID NOT FIND A TRIAL FIELD FILE')
1155 end if
1156 end if
1157
1158 nultrl = 0
1159 ierr = fnom(nultrl,trim(trialfile),'RND+OLD+R/O',0)
1160 ierr = fstouv(nultrl,'RND+OLD')
1161
1162 !- Determine grid size and EZSCINT ID
1163 idateo = -1
1164 cletiket = ' '
1165 ip1 = -1
1166 ip2 = -1
1167 ip3 = -1
1168 cltypvar = ' '
1169 clnomvar = 'MG'
1170
1171 key = fstinf( nultrl, & ! IN
1172 ni_trial, nj_trial, ink, & ! OUT
1173 idateo, cletiket, ip1, ip2, ip3, cltypvar, clnomvar ) ! IN
1174
1175 if (key < 0) then
1176 write(*,*)
1177 write(*,*) 'bMatrixHI: Unable to find trial field = ',clnomvar
1178 call utl_abort('BMatrixHI')
1179 end if
1180
1181 ierr = fstprm( key, & ! IN
1182 idateo, deet, npas, ni_trial, nj_trial, ink, nbits, & ! OUT
1183 datyp, ip1, ip2, ip3, cltypvar, clnomvar, cletiket, & ! OUT
1184 clgrtyp, ig1, ig2, ig3, & ! OUT
1185 ig4, swa, lng, dltf, ubc, extra1, extra2, extra3 ) ! OUT
1186
1187 allocate(TrialLandSeaMask(ni_trial, nj_trial))
1188 allocate(TrialSeaIceMask(ni_trial, nj_trial))
1189
1190 idateo = -1
1191 cletiket = ' '
1192 ip1 = -1
1193 ip2 = -1
1194 ip3 = -1
1195 cltypvar = ' '
1196 clnomvar = 'MG'
1197 ierr = fstlir(TrialLandSeaMask, nultrl, ini, inj, ink, &
1198 idateo ,cletiket, ip1, ip2, ip3, cltypvar, clnomvar)
1199 if ( ierr < 0 ) then
1200 write(*,*)
1201 write(*,*) 'bMatrixHI: Unable to read trial field = ',clnomvar
1202 call utl_abort('BMatrixHI : fstlir failed')
1203 end if
1204
1205 if (ini /= ni_trial .or. inj /= nj_trial) then
1206 write(*,*)
1207 write(*,*) 'bMatrixHI: Invalid dimensions for ...'
1208 write(*,*) 'nomvar =', trim(clnomvar)
1209 write(*,*) 'etiket =', trim(cletiket)
1210 write(*,*) 'ip1 =', ip1
1211 write(*,*) 'Found ni,nj =', ini, inj
1212 write(*,*) 'Should be =', ni_trial, nj_trial
1213 call utl_abort('bMatrixHI')
1214 end if
1215
1216 clnomvar = 'GL'
1217 ierr = fstlir(TrialSeaIceMask, nultrl, ini, inj, ink, &
1218 idateo ,cletiket, ip1, ip2, ip3, cltypvar, clnomvar)
1219 if ( ierr < 0 ) then
1220 write(*,*)
1221 write(*,*) 'bMatrixHI: Unable to read trial field = ',clnomvar
1222 call utl_abort('BMatrixHI : fstlir failed')
1223 end if
1224
1225 if (ini /= ni_trial .or. inj /= nj_trial) then
1226 write(*,*)
1227 write(*,*) 'bMatrixHI: Invalid dimensions for ...'
1228 write(*,*) 'nomvar =', trim(clnomvar)
1229 write(*,*) 'etiket =', trim(cletiket)
1230 write(*,*) 'ip1 =', ip1
1231 write(*,*) 'Found ni,nj =', ini, inj
1232 write(*,*) 'Should be =', ni_trial, nj_trial
1233 call utl_abort('bMatrixHI')
1234 end if
1235
1236 TrialGridID = ezqkdef( ni_trial, nj_trial, clgrtyp, ig1, ig2, ig3, ig4, nultrl ) ! IN
1237
1238 ierr = fstfrm(nultrl)
1239 ierr = fclos(nultrl)
1240
1241 !- Interpolate to the Analysis Grid
1242 allocate(AnalLandSeaMask(ni_l, nj_l))
1243 allocate(AnalSeaIceMask(ni_l, nj_l))
1244
1245 ierr = ezdefset(AnalGridID , TrialGridID ) ! IN, IN
1246
1247 ! Nearest-neighbor interpolation
1248 ierr = int_hInterpScalar(AnalLandSeaMask, TrialLandSeaMask, interpDegree='NEAREST') ! OUT, IN
1249 ierr = int_hInterpScalar(AnalSeaIceMask , TrialSeaIceMask, interpDegree='NEAREST') ! OUT, IN
1250
1251 deallocate(TrialLandSeaMask)
1252 deallocate(TrialSeaIceMask)
1253
1254 !- Modify the input/regridded TG Std. Dev.
1255 do jlat = 1, nj_l
1256 do jlon = 1,ni_l
1257
1258 if ( AnalLandSeaMask(jlon,jlat) > 0.1 ) then
1259 ! We take this as a land point.
1260 ! Force std. dev. to capping value
1261 tgstdbg(jlon,jlat) = rlimsuptg
1262 else if ( AnalSeaIceMask(jlon,jlat) > 0.2 ) then
1263 ! We have significant sea ice on this sea point.
1264 ! Force std. dev. to capping value
1265 tgstdbg(jlon,jlat) = rlimsuptg
1266 else
1267 ! We have an open water point. Make sure that the std. dev. is realistic
1268 ! 1.55 is slightly above the max value over the ocean in the legacy TG Std. Dev. field (in 2014)
1269 if ( tgstdbg(jlon,jlat) > 1.55d0 ) then
1270 tgstdbg(jlon,jlat) = 0.9d0
1271 end if
1272 end if
1273
1274 end do
1275 end do
1276
1277 !- Write the modified Std. Dev.
1278 if ( mmpi_myid == 0 ) then
1279 allocate(tgstdbg_tmp(ni_l,nj_l))
1280 do jlat = 1, nj_l
1281 do jlon = 1,ni_l
1282 tgstdbg_tmp(jlon,jlat) = tgstdbg(jlon,jlat)
1283 end do
1284 end do
1285
1286 iultg = 0
1287 ierr = fnom(iultg,'tg_stddev_of_the_day.fst','RND',0)
1288 ierr = fstouv(iultg,'RND')
1289
1290 ini = ni_l
1291 inj = nj_l
1292 ink = 1
1293 ip1 = 0
1294 ip2 = 0
1295 ip3 = 0
1296 idateo = 0
1297 cltypvar = 'E'
1298 clnomvar = 'TG'
1299 cletiket = 'TWEAK_STDDEV'
1300 clgrtyp = 'G'
1301 ig1 = 0
1302 ig2 = 0
1303 ig3 = 0
1304 ig4 = 0
1305 idatyp = 1
1306
1307 ierr = utl_fstecr(tgstdbg_tmp, -32, iultg, idateo, &
1308 0, 0, ini, inj, ink, ip1, ip2, ip3, cltypvar, &
1309 clnomvar,cletiket,clgrtyp,ig1, ig2, ig3, ig4, idatyp, &
1310 .true.)
1311
1312 do jlat = 1, nj_l
1313 do jlon = 1,ni_l
1314 tgstdbg_tmp(jlon,jlat) = AnalLandSeaMask(jlon,jlat)
1315 end do
1316 end do
1317 cltypvar = 'P'
1318 clnomvar = 'MG'
1319 cletiket = 'TRIAL2ANAL'
1320 ierr = utl_fstecr(tgstdbg_tmp, -32, iultg, idateo, &
1321 0, 0, ini, inj, ink, ip1, ip2, ip3, cltypvar, &
1322 clnomvar,cletiket,clgrtyp,ig1, ig2, ig3, ig4, idatyp, &
1323 .true.)
1324
1325 do jlat = 1, nj_l
1326 do jlon = 1,ni_l
1327 tgstdbg_tmp(jlon,jlat) = AnalSeaIceMask(jlon,jlat)
1328 end do
1329 end do
1330 cltypvar = 'P'
1331 clnomvar = 'GL'
1332 cletiket = 'TRIAL2ANAL'
1333 ierr = utl_fstecr(tgstdbg_tmp, -32, iultg, idateo, &
1334 0, 0, ini, inj, ink, ip1, ip2, ip3, cltypvar, &
1335 clnomvar,cletiket,clgrtyp,ig1, ig2, ig3, ig4, idatyp, &
1336 .true.)
1337
1338 ierr = fstfrm(iultg)
1339 ierr = fclos(iultg)
1340
1341 deallocate(tgstdbg_tmp)
1342 end if
1343
1344 deallocate(AnalLandSeaMask)
1345 deallocate(AnalSeaIceMask)
1346
1347 end if ! if ( TweakTG )
1348
1349 !
1350 !- 2. Compute correlations in spectral space (CORNS)
1351 !
1352 zgd(:,:,:) = 0.0d0
1353 zsp_mpilocal(:,:,:) = 0.0d0
1354 allocate(my_zsp_mpiglobal(nla_mpiglobal,2,1))
1355 my_zsp_mpiglobal(:,:,:) = 0.0d0
1356
1357 do jla = 1, nla_mpiglobal
1358 cortgg(jla,1) = 0.0d0
1359 cortgg(jla,2) = 0.0d0
1360 enddo
1361
1362 ! 2.4.2 Compute correlations in physical space
1363 rcscltg_vec(:) = rcscltg(1)
1364 call BHI_calccorr(zgd,rcscltg_vec,nlev_T_even)
1365
1366 ! 2.4.3 Bring back the result in spectral space
1367 call gst_setID(gstID2)
1368 call gst_reespe(zsp_mpilocal,zgd)
1369
1370 ! and make the result mpiglobal
1371 do jm = mymBeg, mymEnd, mymSkip
1372 do jn = mynBeg, mynEnd, mynSkip
1373 if(jm.le.jn) then
1374 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
1375 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
1376 my_zsp_mpiglobal(ila_mpiglobal,:,1) = zsp_mpilocal(ila_mpilocal,:,1)
1377 endif
1378 enddo
1379 enddo
1380 nsize = 2*nla_mpiglobal
1381 call rpn_comm_allreduce(my_zsp_mpiglobal(:,:,1),zsp_mpiglobal(:,:,1),nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
1382 deallocate(my_zsp_mpiglobal)
1383 ! 2.4.4 Check positiveness
1384 llpb = .false.
1385 do jla = 1, ntrunc+1
1386 zabs = abs(zsp_mpiglobal(jla,1,1))
1387 llpb = llpb.or.((zsp_mpiglobal(jla,1,1).lt.0.).and.(zabs.gt.epsilon(zabs)))
1388 enddo
1389 if(llpb) then
1390 call utl_abort(' AUTOCORRELATION NEGATIVES')
1391 endif
1392 do jla = 1, ntrunc+1
1393 zsp_mpiglobal(jla,1,1) = abs(zsp_mpiglobal(jla,1,1))
1394 enddo
1395
1396 zpole = 0.d0
1397 do jla = 1, ntrunc+1
1398 jn = jla-1
1399 zpole = zpole + zsp_mpiglobal(jla,1,1)*sqrt((2.d0*jn+1.d0)/2.d0)
1400 enddo
1401 if(zpole.le.0.d0) then
1402 call utl_abort('POLE VALUE NEGATIVE IN SUTG')
1403 endif
1404 do jla = 1, ntrunc+1
1405 zsp_mpiglobal(jla,1,1) = zsp_mpiglobal(jla,1,1)/zpole
1406 zsp_mpiglobal(jla,2,1) = zsp_mpiglobal(jla,2,1)/zpole
1407 enddo
1408
1409 ! 2.4.5 Correlation
1410 do jm = 0, ntrunc
1411 do jn = jm, ntrunc
1412 jla = gst_getNIND(jm,gstID) + jn - jm
1413 dlfac = 0.5d0/dsqrt((2*jn+1.d0)/2.d0)
1414 cortgg(jla,1) = dlfac * zsp_mpiglobal(jn+1,1,1)
1415 cortgg(jla,2) = dlfac * zsp_mpiglobal(jn+1,1,1)
1416 enddo
1417 enddo
1418
1419 ! 2.5. For zonal modes : set to zero the imaginary part and set the correct factor 1.0 for the real part
1420 do jla = 1, ntrunc + 1
1421 cortgg(jla,1) = 0.5d0*cortgg(jla,1)
1422 cortgg(jla,2) = 0.0d0
1423 enddo
1424
1425 ! 2.6. Result in corns array
1426 do jn = 0, ntrunc
1427 ila_mpiglobal = jn + 1
1428 corns(nspositTG,nspositTG,jn) = 2.d0*cortgg(ila_mpiglobal,1)
1429 enddo
1430
1431 deallocate(dltg)
1432 !write(*,*)'DONE in SUTG'
1433
1434 END SUBROUTINE BHI_sutg
1435
1436
1437 SUBROUTINE BHI_convol
1438 implicit none
1439
1440 ! Locals:
1441 real(8) dlfact2,dlc,dsummed
1442 real(8) dtlen,zr,dlfact
1443 integer jn,jlat,jk
1444 real(8) zsp(0:ntrunc,nkgdim),zgr(nj_l,nkgdim)
1445 real(8) dlwti(nj_l),zrmu(nj_l)
1446 real(8) :: RPORVO = 6000.D3
1447 real(8) :: RPORDI = 6000.D3
1448 real(8) :: RPORTT = 3000.D3
1449 real(8) :: RPORQ = 3000.D3
1450 real(8) :: RPORPS = 3000.D3
1451
1452 do jlat = 1, nj_l
1453 dlwti(jlat) = gst_getrwt(jlat,gstID)
1454 zrmu(jlat) = gst_getrmu(jlat,gstID)
1455 end do
1456
1457! 1.2 CONVERT THE CORRELATIONS IN SPECTRAL SPACE INTO SPECTRAL
1458! COEFFICIENTS OF THE CORRELATION FUNCTION AND FUNCTION TO BE
1459! SELF-CONVOLVED
1460 do jn = 0, ntrunc
1461 dlfact = ((2.0d0*jn +1.0d0)/2.0d0)**(0.25d0)
1462 dlfact2= ((2.0d0*jn +1.0d0)/2.0d0)**(0.25d0)
1463 do jk = 1, nkgdim
1464 zsp(jn,jk) = rstddev(jk,jn)*dlfact*dlfact2
1465 enddo
1466 enddo
1467
1468 ! Transform to physical space
1469 call gst_zleginv(gstID,zgr,zsp,nkgdim)
1470
1471 ! Truncate in horizontal extent with Gaussian window
1472 do jk = 1, nkgdim
1473 if (jk.ge.nspositVO.and.jk.lt.nspositVO+nlev_M) then
1474 dtlen = rporvo
1475 elseif (jk.ge.nspositDI.and.jk.lt.nspositDI+nlev_M) then
1476 dtlen = rpordi
1477 elseif (jk.ge.nspositTT.and.jk.lt.nspositTT+nlev_T) then
1478 dtlen = rportt
1479 elseif (jk.ge.nspositQ.and.jk.lt.nspositQ+nlev_T) then
1480 dtlen = rporq
1481 elseif (jk == nspositPS) then
1482 dtlen = rporps
1483 endif
1484
1485 if(dtlen.gt.0.0d0) then
1486 dlc = 1.d0/dble(dtlen)
1487 dlc = 0.5d0*dlc*dlc
1488 do jlat = 1, nj_l
1489 zr = ec_ra * acos(zrmu(jlat))
1490 dlfact = dexp(-(zr**2)*dlc)
1491 zgr(jlat,jk) = dlfact*zgr(jlat,jk)
1492 enddo
1493 endif
1494
1495 !write(*,*) 'zeroing length (km)=',jk,dtlen/1000.0
1496 enddo
1497
1498 ! Transform back to spectral space
1499 call gst_zlegdir(gstID,zgr,zsp,nkgdim)
1500
1501 ! Convert back to correlations
1502 do jk = 1, nkgdim
1503 do jn = 0, ntrunc
1504 zsp(jn,jk) = zsp(jn,jk)*(2.0d0/(2.0d0*jn+1.0d0))**(0.25d0)
1505 enddo
1506 enddo
1507
1508 ! PUT BACK INTO RSTDDEV
1509 do jn = 0, ntrunc
1510 do jk = 1, nkgdim
1511 rstddev(jk,jn) = zsp(jn,jk)
1512 enddo
1513 enddo
1514
1515 ! Re-normalize to ensure correlations
1516 do jk = 1, nkgdim
1517 dsummed = 0.d0
1518 do jn = 0, ntrunc
1519 dsummed = dsummed+ dble(rstddev(jk,jn)**2)*sqrt(((2.d0*jn)+1.d0)/2.d0)
1520 enddo
1521 dsummed = sqrt(dsummed)
1522 do jn = 0, ntrunc
1523 if(dsummed.gt.1.d-30) rstddev(jk,jn) = rstddev(jk,jn)/dsummed
1524 enddo
1525 enddo
1526
1527 ! CONVERT THE SPECTRAL COEFFICIENTS OF THE CORRELATION FUNCTION
1528 ! . BACK INTO CORRELATIONS OF SPECTRAL COMPONENTS
1529 do jn = 0, ntrunc
1530 dlfact = sqrt(0.5d0)*(1.0d0/((2.0d0*jn+1)/2.0d0))**0.25d0
1531 do jk = 1, nkgdim
1532 rstddev(jk,jn) = rstddev(jk,jn)*dlfact
1533 enddo
1534 enddo
1535
1536 END SUBROUTINE BHI_convol
1537
1538
1539 SUBROUTINE BHI_setCrossCorr(kn)
1540 implicit none
1541
1542 ! Arguments:
1543 integer, intent(in) :: kn
1544
1545 ! Locals:
1546 integer :: jblock1, inbrblock, jblock2
1547 integer :: jk1, jk2, nlev_all(numvar3d), levOffset(numvar3d+1)
1548
1549 inbrblock = numvar3d
1550 nlev_all(1) = nLev_M
1551 nlev_all(2) = nLev_M
1552 nlev_all(3) = nLev_T
1553 nlev_all(4) = nLev_T
1554 levOffset(1) = 0
1555 levOffset(2) = 1*nLev_M
1556 levOffset(3) = 2*nLev_M
1557 levOffset(4) = 2*nLev_M+1*nLev_T
1558 levOffset(5) = 2*nLev_M+2*nLev_T
1559
1560 ! Set cross-variable correlations to 0 ...
1561 do jblock1 = 1, inbrblock
1562 do jblock2 = 1, inbrblock
1563 if (jblock1.ne.jblock2) then
1564 do jk2 = 1, nlev_all(jblock2)
1565 do jk1 = 1, nlev_all(jblock1)
1566 corns(jk1 + levOffset(jblock1),jk2 + levOffset(jblock2),kn) = 0.0d0
1567 enddo
1568 enddo
1569 endif
1570 enddo
1571 enddo
1572
1573 ! ... but T'ln(ps') correlations
1574 do jk2 = 1, nkgdim
1575 do jk1 = levOffset(5)+1, levOffset(5)+numvar2d
1576 if ((jk1.ne.nspositPS.or.jk2.lt.nspositTT.or. &
1577 jk2.ge.(nspositTT+nlev_T)).and.(jk1.ne.jk2)) then
1578 corns(jk1,jk2,kn) = 0.0d0
1579 endif
1580 enddo
1581 enddo
1582
1583 do jk2 = levOffset(5)+1, levOffset(5)+numvar2d
1584 do jk1 = 1, nkgdim
1585 if ((jk2.ne.nspositPS.or.jk1.lt.nspositTT.or. &
1586 jk1.ge.(nspositTT+nlev_T)) .and.(jk1.ne.jk2)) then
1587 corns(jk1,jk2,kn) = 0.0d0
1588 endif
1589 enddo
1590 enddo
1591
1592 END SUBROUTINE BHI_setCrossCorr
1593
1594
1595 SUBROUTINE BHI_READCORNS2
1596 implicit none
1597
1598 ! Locals:
1599 integer :: kip1
1600 integer :: jn, istdkey,icornskey
1601 integer :: iksdim,jcol,jrow
1602 real(8), allocatable, dimension(:) :: zstdsrc
1603 real(8), allocatable, dimension(:,:) :: zcornssrc
1604 ! standard file variables
1605 integer :: ini,inj,ink
1606 integer :: ip1,ip2,ip3
1607 integer :: idateo
1608 character(len=2) :: cltypvar
1609 character(len=4) :: clnomvar
1610 character(len=12) :: cletiket
1611
1612 iksdim = 2*nlev_M+2*nlev_T+1 ! assume 4 3d variables and 1 2d variable (TG not included)
1613 allocate(zcornssrc(iksdim,iksdim))
1614 allocate(zstdsrc(iksdim))
1615
1616 kip1 = -1
1617
1618 do jn = 0, ntrunc
1619
1620 ! Looking for FST record parameters..
1621
1622 idateo = -1
1623 cletiket = 'RSTDDEV'
1624 ip1 = kip1
1625 ip2 = jn
1626 ip3 = -1
1627 cltypvar = 'X'
1628 clnomvar = 'SS'
1629
1630 istdkey = utl_fstlir(ZSTDSRC,nulbgst,INI,INJ,INK,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1631
1632 if(istdkey .lt.0 ) then
1633 call utl_abort('READCORNS2: Problem with background stat file')
1634 endif
1635
1636 if (ini .ne. iksdim) then
1637 call utl_abort('READCORNS2: BG stat levels inconsitencies')
1638 endif
1639
1640 ! Looking for FST record parameters..
1641
1642 idateo = -1
1643 cletiket = 'CORRNS'
1644 ip1 = kip1
1645 IP2 = JN
1646 ip3 = -1
1647 cltypvar = 'X'
1648 clnomvar = 'ZZ'
1649 icornskey = utl_fstlir(ZCORNSSRC,nulbgst,INI,INJ,INK,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1650
1651 if(icornskey .lt.0 ) then
1652 call utl_abort('READCORNS2: Problem with background stat file')
1653 endif
1654
1655 if (ini .ne. iksdim .or. inj .ne. iksdim) then
1656 call utl_abort('READCORNS2: BG stat levels inconsitencies')
1657 endif
1658
1659 do jcol = 1, nkgdim2
1660 rstddev(jcol,jn) = 0.0d0
1661 do jrow = 1, nkgdim2
1662 corns(jrow,jcol,jn) = 0.0d0
1663 enddo
1664 enddo
1665
1666 do jcol = 1, iksdim
1667 do jrow = 1, iksdim
1668 corns(jrow,jcol,jn) = zcornssrc(jrow,jcol)
1669 enddo
1670 enddo
1671
1672 ! Set cross-variable correlations to zero except between T' and ln(ps')
1673 call BHI_setcrosscorr(jn)
1674
1675 do jrow = 1, iksdim
1676 rstddev(jrow,jn) = zstdsrc(jrow)
1677 enddo
1678
1679 enddo
1680
1681 ! Apply convolution to RSTDDEV correlations
1682
1683 call BHI_convol
1684
1685 do jn = 0, ntrunc
1686
1687 ! Re-build of correlation matrix: factorization of corns with convoluted RSTDDEV
1688 do jcol = 1, nkgdim
1689 do jrow = 1, nkgdim
1690 corns(jrow,jcol,jn) = rstddev(jrow,jn) * corns(jrow,jcol,jn)* rstddev(jcol,jn)
1691 enddo
1692 enddo
1693
1694 enddo
1695
1696 deallocate(zcornssrc)
1697 deallocate(zstdsrc)
1698
1699 !write(*,*) 'Done in READCORNS2'
1700 END SUBROUTINE BHI_READCORNS2
1701
1702
1703 SUBROUTINE BHI_RDSPSTD
1704 implicit none
1705
1706 ! Locals:
1707 integer, parameter :: inbrvar3d=5
1708 integer, parameter :: inbrvar2d=2
1709 integer :: jvar,jn,inix,injx,inkx
1710 integer :: ikey, jlevo, jlat,firstn,lastn
1711 real(8) :: zsp(0:ntrunc,max(nlev_M,nlev_T)),zspbuf(max(nlev_M,nlev_T))
1712 real(8) :: zgr(nj_l,max(nlev_M,nlev_T))
1713 character(len=4) :: varName3d(inbrvar3d),varName2d(inbrvar2d)
1714 ! standard file variables
1715 integer :: ini,inj,ink
1716 integer :: ip1,ip2,ip3
1717 integer :: idate(100),nlev_MT
1718 character(len=2) :: cltypvar
1719 character(len=4) :: clnomvar
1720 character(len=12) :: cletiket
1721 integer :: fstinf
1722
1723 data varName3d/'PP ','UC ','UT ','LQ ','TB '/
1724 data varName2d/'UP ','PB '/
1725
1726 rgsig(:,:) = 0.0d0
1727 rgsigtb(:,:) = 0.0d0
1728 rgsigpsb(:) = 0.0d0
1729
1730! 2. Reading the data
1731
1732 idate(1) = -1
1733 ip1 = -1
1734 ip2 = -1
1735 ip3 = -1
1736
1737 cletiket = 'SPSTDDEV'
1738 cltypvar = 'X'
1739
1740 do jvar = 1, inbrvar3d
1741 clnomvar = varName3d(jvar)
1742 if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
1743 nlev_MT = nlev_M
1744 else
1745 nlev_MT = nlev_T
1746 endif
1747 firstn = -1
1748 do jn = 0, ntrunc
1749 ip2 = jn
1750 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1751
1752 if(ikey .ge.0 ) then
1753 ikey = utl_fstlir(zspbuf(1:nlev_MT),nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1754 else
1755 if(firstn == -1) firstn = jn
1756 lastn = jn
1757 zspbuf(:) = 0.0d0
1758 endif
1759
1760 if (ini .ne. nlev_MT) then
1761 call utl_abort('RDSPSTD: BG stat levels inconsitencies')
1762 endif
1763
1764 do jlevo = 1, nlev_MT
1765 zsp(jn,jlevo) = zspbuf(jlevo)
1766 enddo
1767 enddo
1768 if(mmpi_myid == 0.and.firstn.ne.-1) then
1769 write(*,*) 'WARNING: CANNOT FIND SPSTD FOR ',clnomvar, &
1770 ' AT N BETWEEN ',firstn,' AND ',lastn,', SETTING TO ZERO!!!'
1771 endif
1772
1773 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
1774
1775 if(clnomvar == 'PP') then
1776 do jlat = 1, nj_l
1777 do jlevo = 1, nlev_M
1778 rgsiguu(jlat,jlevo) = zgr(jlat,jlevo)
1779 enddo
1780 enddo
1781 elseif(clnomvar == 'UC' .or. clnomvar == 'CC') then
1782 do jlat = 1, nj_l
1783 do jlevo = 1, nlev_M
1784 rgsigvv(jlat,jlevo) = zgr(jlat,jlevo)
1785 enddo
1786 enddo
1787 elseif(clnomvar == 'UT') then
1788 do jlat = 1, nj_l
1789 do jlevo = 1, nlev_T
1790 rgsigtt(jlat,jlevo) = zgr(jlat,jlevo)
1791 enddo
1792 enddo
1793 elseif(clnomvar == 'TB') then
1794 do jlat = 1, nj_l
1795 do jlevo = 1, nlev_T
1796 rgsigtb(jlat,jlevo) = zgr(jlat,jlevo)
1797 enddo
1798 enddo
1799 elseif(clnomvar == 'LQ') then
1800 do jlat = 1, nj_l
1801 do jlevo = 1, nlev_T
1802 rgsigq(jlat,jlevo) = zgr(jlat,jlevo)
1803 enddo
1804 enddo
1805 endif
1806
1807 enddo
1808
1809 nlev_MT = 1
1810 do jvar = 1, inbrvar2d
1811 clnomvar = varName2d(jvar)
1812 firstn = -1
1813 do jn = 0, ntrunc
1814 ip2 = jn
1815 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1816
1817 if(ikey .ge.0 ) then
1818 ikey = utl_fstlir(zspbuf,nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1819 else
1820 if(firstn == -1) firstn = jn
1821 lastn = jn
1822 zspbuf(:) = 0.0d0
1823 endif
1824
1825 zsp(jn,1) = zspbuf(1)
1826
1827 enddo
1828 if(mmpi_myid == 0.and.firstn.ne.-1) then
1829 write(*,*) 'WARNING: CANNOT FIND SPSTD FOR ',clnomvar, &
1830 ' AT N BETWEEN ',firstn,' AND ',lastn,', SETTING TO ZERO!!!'
1831
1832 endif
1833
1834 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
1835
1836 if(clnomvar == 'UP') then
1837 do jlat = 1, nj_l
1838 rgsigps(jlat) = zgr(jlat,1)*100.0d0
1839 enddo
1840 endif
1841 if(clnomvar == 'PB') then
1842 do jlat = 1, nj_l
1843 rgsigpsb(jlat) = zgr(jlat,1)*100.0d0
1844 enddo
1845 endif
1846
1847 enddo
1848
1849 END SUBROUTINE BHI_RDSPSTD
1850
1851
1852 SUBROUTINE BHI_RDSTD
1853 implicit none
1854
1855 ! Locals:
1856 integer, parameter :: inbrvar3d=5
1857 integer, parameter :: inbrvar2d=2
1858 integer :: jvar,inix,injx,inkx
1859 integer :: ikey
1860 real(8) :: zgr(nj_l,max(nlev_M,nlev_T))
1861 character(len=4) :: varName3d(inbrvar3d),varName2d(inbrvar2d)
1862 ! standard file variables
1863 integer :: ini,inj,ink
1864 integer :: ip1,ip2,ip3
1865 integer :: idate(100),nlev_MT
1866 character(len=2) :: cltypvar
1867 character(len=4) :: clnomvar
1868 character(len=12) :: cletiket
1869 integer :: fstinf
1870
1871 data varName3d/'PP ','UC ','UT ','LQ ','TB '/
1872 data varName2d/'UP ','PB '/
1873
1874 rgsig(:,:) = 0.0d0
1875 rgsigtb(:,:) = 0.0d0
1876 rgsigpsb(:) = 0.0d0
1877
1878! 2. Reading the data
1879
1880 idate(1) = -1
1881 ip1 = -1
1882 ip2 = -1
1883 ip3 = -1
1884
1885 cletiket = 'STDDEV'
1886 cltypvar = 'E'
1887
1888 do jvar = 1, inbrvar3d
1889 clnomvar = varName3d(jvar)
1890 if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
1891 nlev_MT = nlev_M
1892 else
1893 nlev_MT = nlev_T
1894 endif
1895
1896 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1897
1898 if(ikey .ge.0 ) then
1899 ikey = utl_fstlir(zgr(:,1:nlev_MT),nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1900 else
1901 write(*,*) 'RDSTD: could not read varName=',clnomvar
1902 call utl_abort('RDSTD')
1903 endif
1904
1905 if (ink .ne. nlev_MT) then
1906 write(*,*) 'ink, nlev_MT=', ink, nlev_MT
1907 call utl_abort('RDSPSTD: BG stat levels inconsitencies')
1908 endif
1909
1910 if(clnomvar == 'PP') then
1911 rgsiguu(:,:) = zgr(:,1:nlev_M)
1912 elseif(clnomvar == 'UC' .or. clnomvar == 'CC') then
1913 rgsigvv(:,:) = zgr(:,1:nlev_M)
1914 elseif(clnomvar == 'UT') then
1915 rgsigtt(:,:) = zgr(:,1:nlev_T)
1916 elseif(clnomvar == 'TB') then
1917 rgsigtb(:,:) = zgr(:,1:nlev_T)
1918 elseif(clnomvar == 'LQ') then
1919 rgsigq(:,:) = zgr(:,1:nlev_T)
1920 endif
1921
1922 enddo
1923
1924 nlev_MT = 1
1925 do jvar = 1, inbrvar2d
1926 clnomvar = varName2d(jvar)
1927
1928 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1929
1930 if(ikey .ge.0 ) then
1931 ikey = utl_fstlir(zgr,nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1932 else
1933 write(*,*) 'RDSTD: could not read varName=',clnomvar
1934 call utl_abort('RDSTD')
1935 endif
1936
1937 if(clnomvar == 'UP') then
1938 rgsigps(:) = zgr(:,1)*100.0d0
1939 endif
1940 if(clnomvar == 'PB') then
1941 rgsigpsb(:) = zgr(:,1)*100.0d0
1942 endif
1943
1944 enddo
1945
1946 END SUBROUTINE BHI_RDSTD
1947
1948
1949 SUBROUTINE BHI_RDSPSTD_NEWFMT
1950 implicit none
1951
1952 ! Locals:
1953 integer, parameter :: inbrvar3d=5
1954 integer, parameter :: inbrvar2d=2
1955 integer :: jvar,jn,inix,injx,inkx,ntrunc_file
1956 integer :: ikey,jlevo,jlat
1957 real(8) :: zsp(0:ntrunc,max(nlev_M,nlev_T))
1958 real(8), allocatable :: zspbuf(:)
1959 real(8) :: zgr(nj_l,max(nlev_M,nlev_T))
1960 character(len=4) :: varName3d(inbrvar3d),varName2d(inbrvar2d)
1961 ! standard file variables
1962 integer :: ini,inj,ink
1963 integer :: ip1,ip2,ip3
1964 integer :: idate(100),nlev_MT
1965 character(len=2) :: cltypvar
1966 character(len=4) :: clnomvar
1967 character(len=12) :: cletiket
1968 integer :: fstinf
1969
1970 data varName3d/'PP ','UC ','UT ','LQ ','TB '/
1971 data varName2d/'UP ','PB '/
1972
1973 rgsig(:,:) = 0.0d0
1974 rgsigtb(:,:) = 0.0d0
1975 rgsigpsb(:) = 0.0d0
1976
1977! 2. Reading the data
1978
1979 idate(1) = -1
1980 ip2 = -1
1981 ip3 = -1
1982
1983 cletiket = 'SPSTDDEV'
1984 cltypvar = 'X'
1985
1986 ! check if file is old format
1987 ip1 = -1
1988 clnomvar = varName3d(1)
1989 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1990 write(*,*) 'ini,inj,ink=',inix,injx,inkx
1991 if(inix.gt.1) then
1992 write(*,*) 'BHI_RDSPSTD_NEWFMT: ini>1, SPSTDDEV is in old format, calling BHI_RDSPSTD...'
1993 call bhi_rdspstd
1994 return
1995 endif
1996
1997 !write(*,*) 'Reading 3D variables'
1998 do jvar = 1, inbrvar3d
1999 clnomvar = varName3d(jvar)
2000 if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
2001 nlev_MT = nlev_M
2002 else
2003 nlev_MT = nlev_T
2004 endif
2005 !write(*,*)'Reading ',clnomvar
2006 do jlevo = 1, nlev_MT
2007 if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
2008 ip1 = vco_anl%ip1_M(jlevo)
2009 else
2010 ip1 = vco_anl%ip1_T(jlevo)
2011 endif
2012 ikey = fstinf(nulbgst,inix,ntrunc_file,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2013 ntrunc_file = ntrunc_file-1
2014
2015 allocate(zspbuf(0:ntrunc_file))
2016 if(ikey .ge.0 ) then
2017 ikey = utl_fstlir(zspbuf(0:ntrunc_file),nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2018 else
2019 write(*,*) 'RDSPSTD_NEWFMT: ',jvar,clnomvar,nlev_MT,jlevo,ikey,ntrunc,ntrunc_file
2020 call utl_abort('RDSPSTD_NEWFMT: SPSTDDEV record not found')
2021 endif
2022
2023 zsp(:,jlevo) = 0.0d0
2024 do jn = 0, min(ntrunc,ntrunc_file)
2025 zsp(jn,jlevo) = zspbuf(jn)
2026 enddo
2027 deallocate(zspbuf)
2028
2029 enddo
2030
2031 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
2032
2033 if(clnomvar == 'PP') then
2034 do jlat = 1, nj_l
2035 do jlevo = 1, nlev_M
2036 rgsiguu(jlat,jlevo) = zgr(jlat,jlevo)
2037 enddo
2038 enddo
2039 elseif(clnomvar == 'UC' .or. clnomvar == 'CC') then
2040 do jlat = 1, nj_l
2041 do jlevo = 1, nlev_M
2042 rgsigvv(jlat,jlevo) = zgr(jlat,jlevo)
2043 enddo
2044 enddo
2045 elseif(clnomvar == 'UT') then
2046 do jlat = 1, nj_l
2047 do jlevo = 1, nlev_T
2048 rgsigtt(jlat,jlevo) = zgr(jlat,jlevo)
2049 enddo
2050 enddo
2051 elseif(clnomvar == 'TB') then
2052 do jlat = 1, nj_l
2053 do jlevo = 1, nlev_T
2054 rgsigtb(jlat,jlevo) = zgr(jlat,jlevo)
2055 enddo
2056 enddo
2057 elseif(clnomvar == 'LQ') then
2058 do jlat = 1, nj_l
2059 do jlevo = 1, nlev_T
2060 rgsigq(jlat,jlevo) = zgr(jlat,jlevo)
2061 enddo
2062 enddo
2063 endif
2064
2065 enddo
2066
2067 nlev_MT = 1
2068 do jvar = 1, inbrvar2d
2069 clnomvar = varName2d(jvar)
2070 ip1 = -1
2071 ikey = fstinf(nulbgst,inix,ntrunc_file,inkx,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2072 ntrunc_file = ntrunc_file-1
2073
2074 allocate(zspbuf(0:ntrunc_file))
2075
2076 if(ikey .ge.0 ) then
2077 ikey = utl_fstlir(zspbuf(0:ntrunc_file),nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2078 else
2079 write(*,*) 'RDSPSTD_NEWFMT: ',jvar,clnomvar,nlev_MT,jlevo,ikey,ntrunc,ntrunc_file
2080 call utl_abort('RDSPSTD_NEWFMT: SPSTDDEV record not found')
2081 endif
2082
2083 zsp(:,1) = 0.0d0
2084 do jn = 0, min(ntrunc,ntrunc_file)
2085 zsp(jn,1) = zspbuf(jn)
2086 enddo
2087 deallocate(zspbuf)
2088
2089 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
2090
2091 if(clnomvar == 'UP') then
2092 do jlat = 1, nj_l
2093 rgsigps(jlat) = zgr(jlat,1)*100.0d0
2094 enddo
2095 endif
2096 if(clnomvar == 'PB') then
2097 do jlat = 1, nj_l
2098 rgsigpsb(jlat) = zgr(jlat,1)*100.0d0
2099 enddo
2100 endif
2101
2102 enddo
2103
2104 END SUBROUTINE BHI_RDSPSTD_NEWFMT
2105
2106
2107 SUBROUTINE BHI_RDSPPTOT
2108 IMPLICIT NONE
2109
2110 ! Locals:
2111 integer :: jn, jk1, jk2, ikey, ilen,jlat,inix,injx,inkx
2112 real(8) :: zsptheta(0:ntrunc,nlev_M)
2113 real(8) :: zgrtheta(nj_l,nlev_M)
2114 real(8) :: zPtoTsrc(nlev_T+1,nlev_M)
2115 real(8) :: zspPtoT(0:ntrunc,nlev_T+1,nlev_M)
2116 real(8) :: zgrPtoT(nj_l,nlev_T+1,nlev_M)
2117 real(8) :: ztheta(nlev_M)
2118 ! standard file variables
2119 integer :: ini,inj,ink
2120 integer :: ip1,ip2,ip3
2121 integer :: idateo
2122 character(len=2) :: cltypvar
2123 character(len=4) :: clnomvar
2124 character(len=12) :: cletiket
2125 integer :: fstinf
2126
2127 ip1 = -1
2128 ip3 = -1
2129 idateo = -1
2130 cletiket = 'SP_THETA'
2131 cltypvar = 'X'
2132 clnomvar = 'ZZ'
2133
2134 ! read spectral coefficients for theta
2135
2136 do jn = 0, ntrunc
2137 ip2 = jn
2138 ikey = fstinf(nulbgst,inix,injx,inkx,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2139
2140 if(ikey .ge.0 ) then
2141 ikey = utl_fstlir(ztheta,nulbgst,ini,inj,ink,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2142 else
2143 if(mmpi_myid == 0) write(*,*) 'WARNING: CANNOT FIND THETA FOR ',jn,' SETTING TO ZERO!!!'
2144 ztheta(:) = 0.0d0
2145 endif
2146
2147 do jk1 = 1, nlev_M
2148 zsptheta(jn,jk1) = ztheta(jk1)
2149 enddo
2150
2151 enddo
2152
2153 ! converting theta in physical space
2154
2155 call gst_zleginv(gstID,zgrtheta,zsptheta,nlev_M)
2156
2157 do jlat = 1, nj_l
2158 do jk1 = 1, nlev_M
2159 tantheta(jk1,jlat) = tan(zgrtheta(jlat,jk1))
2160 end do
2161 end do
2162
2163 ip1 = -1
2164 ip2 = -1
2165 ip3 = -1
2166 idateo = -1
2167 cletiket = 'SP_PtoT'
2168 cltypvar = 'X'
2169 clnomvar = 'ZZ'
2170
2171 ! read of spectral coefficients for P to T operator
2172
2173 do jn = 0, ntrunc
2174 ip2 = jn
2175 ikey = fstinf(nulbgst,inix,injx,inkx,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2176
2177 if(ikey .ge.0 ) then
2178 ikey = utl_fstlir(zPtoTsrc,nulbgst,ini,inj,ink,idateo,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
2179 else
2180 if(mmpi_myid == 0) write(*,*) 'WARNING: CANNOT FIND P_to_T FOR ',jn,' SETTING TO ZERO!!!'
2181 zPtoTsrc(:,:) = 0.0d0
2182 endif
2183
2184 do jk2 = 1, nlev_M
2185 do jk1 = 1, nlev_T+1
2186 zspPtoT(jn,jk1,jk2) = zPtoTsrc(jk1,jk2)
2187 enddo
2188 enddo
2189
2190 enddo
2191
2192 ilen = nlev_M*(nlev_T+1)
2193 call gst_zleginv(gstID,zgrPtoT,zspPtoT,ilen)
2194
2195 do jlat = 1, nj_l
2196 do jk2 = 1, nlev_M
2197 do jk1 = 1, nlev_T+1
2198 PtoT(jk1,jk2,jlat) = zgrPtoT(jlat,jk1,jk2)
2199 end do
2200 end do
2201 enddo
2202
2203 END SUBROUTINE BHI_RDSPPTOT
2204
2205
2206 SUBROUTINE BHI_truncateCV(controlVector_inout,ntruncCut)
2207 implicit none
2208 !
2209 !:Purpose: Set to zero all coefficients with total wavenumber > ntruncCut
2210 !
2211
2212 ! Arguments:
2213 real(8), pointer, intent(inout) :: controlVector_inout(:)
2214 integer, intent(in) :: ntruncCut
2215
2216 ! Locals:
2217 integer :: jn, jm, ila_mpiglobal, ila_mpilocal, jlev, jdim
2218
2219 if(.not. initialized) then
2220 if(mmpi_myid == 0) write(*,*) 'bhi_truncateCV: bMatrixHI not initialized'
2221 return
2222 endif
2223
2224 if(ntruncCut.gt.ntrunc) then
2225 write(*,*) ntruncCut, ntrunc
2226 call utl_abort('bhi_truncateCV: ntruncCut is greater than ntrunc!')
2227 endif
2228
2229 jdim = 0
2230 do jlev = 1, nkgdimSqrt
2231 do jm = mymBeg, mymEnd, mymSkip
2232 do jn = mynBeg, mynEnd, mynSkip
2233 if(jm.le.jn) then
2234 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
2235 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
2236 if(jm == 0) then
2237 ! only real component for jm=0
2238 jdim = jdim + 1
2239 if(jn.gt.ntruncCut) controlVector_inout(jdim) = 0.0d0
2240 else
2241 ! both real and imaginary components for jm>0
2242 jdim = jdim + 1
2243 if(jn.gt.ntruncCut) controlVector_inout(jdim) = 0.0d0
2244 jdim = jdim + 1
2245 if(jn.gt.ntruncCut) controlVector_inout(jdim) = 0.0d0
2246 endif
2247 endif
2248 enddo
2249 enddo
2250 enddo
2251
2252 END SUBROUTINE BHI_truncateCV
2253
2254
2255 SUBROUTINE BHI_bSqrt(controlVector_in, statevector, stateVectorRef_opt)
2256 implicit none
2257
2258 ! Arguments:
2259 real(8), intent(in) :: controlVector_in(cvDim_mpilocal)
2260 type(struct_gsv), intent(inout) :: statevector
2261 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
2262
2263 ! Locals:
2264 real(8),allocatable :: gd_out(:,:,:)
2265 real(8) :: hiControlVector(nla_mpilocal,2,nkgdimSqrt)
2266
2267 if(mmpi_myid == 0) write(*,*) 'bhi_bsqrt: starting'
2268 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2269
2270 if(.not. initialized) then
2271 if(mmpi_myid == 0) write(*,*) 'bMatrixHI not initialized'
2272 return
2273 endif
2274
2275 allocate(gd_out(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim))
2276
2277 call bhi_cain(controlVector_in,hiControlVector)
2278 call bhi_spa2gd(hiControlVector,gd_out)
2279
2280 call copyToStatevector(statevector,gd_out)
2281
2282 if (gsv_varExist(statevector,'HU')) then
2283 call gvt_transform( statevector, & ! INOUT
2284 'LQtoHU_tlm', & ! IN
2285 stateVectorRef_opt=stateVectorRef_opt ) ! IN
2286 end if
2287
2288 deallocate(gd_out)
2289
2290 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2291 if(mmpi_myid == 0) write(*,*) 'bhi_bsqrt: done'
2292
2293 END SUBROUTINE BHI_bSqrt
2294
2295
2296 SUBROUTINE BHI_bSqrtAd(statevector, controlVector_out, stateVectorRef_opt)
2297 implicit none
2298
2299 ! Arguments:
2300 real(8), intent(inout) :: controlVector_out(cvDim_mpilocal)
2301 type(struct_gsv), intent(inout) :: statevector
2302 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
2303
2304 ! Locals:
2305 real(8), allocatable :: gd_in(:,:,:)
2306 real(8) :: hiControlVector(nla_mpilocal,2,nkgdimSqrt)
2307
2308 if(.not. initialized) then
2309 if(mmpi_myid == 0) write(*,*) 'bMatrixHI not initialized'
2310 return
2311 endif
2312
2313 if(mmpi_myid == 0) write(*,*) 'bhi_bsqrtad: starting'
2314 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2315
2316 allocate(gd_in(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim))
2317 gd_in(:,:,:) = 0.d0
2318 if (gsv_varExist(statevector,'HU')) then
2319 call gvt_transform( statevector, & ! INOUT
2320 'LQtoHU_ad', & ! IN
2321 stateVectorRef_opt=stateVectorRef_opt ) ! IN
2322 end if
2323
2324 call copyFromStatevector(statevector,gd_in)
2325
2326 call bhi_spa2gdad(gd_in,hiControlVector)
2327
2328 call bhi_cainad(hiControlVector,controlVector_out)
2329
2330 deallocate(gd_in)
2331
2332 if(mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2333 if(mmpi_myid == 0) write(*,*) 'bhi_bsqrtad: done'
2334
2335 END SUBROUTINE BHI_bSqrtAd
2336
2337
2338 SUBROUTINE copyToStatevector(statevector,gd)
2339 implicit none
2340
2341 ! Arguments:
2342 type(struct_gsv), intent(inout) :: statevector
2343 real(8), intent(in) :: gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
2344
2345 ! Locals:
2346 integer :: jlon, jlev, jlev2, jlat, jvar, ilev1, ilev2
2347 real(8), pointer :: field_r8(:,:,:)
2348 real(4), pointer :: field_r4(:,:,:)
2349
2350 do jvar = 1, vnl_numvarmax
2351 if(gsv_varExist(statevector,vnl_varNameList(jvar))) then
2352 if (gsv_getDataKind(statevector) == 8) then
2353 call gsv_getField(statevector,field_r8,vnl_varNameList(jvar))
2354 else
2355 call gsv_getField(statevector,field_r4,vnl_varNameList(jvar))
2356 end if
2357 if(vnl_varNameList(jvar) == 'UU ') then
2358 ilev1 = nspositVO
2359 elseif(vnl_varNameList(jvar) == 'VV ') then
2360 ilev1 = nspositDI
2361 elseif(vnl_varNameList(jvar) == 'TT ') then
2362 ilev1 = nspositTT
2363 elseif(vnl_varNameList(jvar) == 'HU ' .or. vnl_varNameList(jvar) == 'LQ ') then
2364 ilev1 = nspositQ
2365 elseif(vnl_varNameList(jvar) == 'P0 ') then
2366 ilev1 = nspositPS
2367 elseif(vnl_varNameList(jvar) == 'TG ') then
2368 ilev1 = nspositTG
2369 else
2370 ! Cycle (instead of abort) to allow for non-NWP assimilation (e.g. chemical data assimilation)
2371! call utl_abort('bmatrixhi_mod: copyToStatevector: No covariances available for variable:' // vnl_varNameList(jvar))
2372 cycle
2373 endif
2374 ilev2 = ilev1 - 1 + gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
2375 if (gsv_getDataKind(statevector) == 8) then
2376 do jlev = ilev1, ilev2
2377 jlev2 = jlev-ilev1+1
2378 do jlat = myLatBeg, myLatEnd
2379 do jlon = myLonBeg, myLonEnd
2380 field_r8(jlon,jlat,jlev2) = gd(jlon,jlat,jlev)
2381 enddo
2382 enddo
2383 enddo
2384 else
2385 do jlev = ilev1, ilev2
2386 jlev2 = jlev-ilev1+1
2387 do jlat = myLatBeg, myLatEnd
2388 do jlon = myLonBeg, myLonEnd
2389 field_r4(jlon,jlat,jlev2) = gd(jlon,jlat,jlev)
2390 enddo
2391 enddo
2392 enddo
2393 end if
2394 endif
2395 enddo
2396
2397 END SUBROUTINE copyToStatevector
2398
2399
2400 SUBROUTINE copyFromStatevector(statevector,gd)
2401 implicit none
2402
2403 ! Arguments:
2404 type(struct_gsv), intent(inout) :: statevector
2405 real(8), intent(out) :: gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
2406
2407 ! Locals:
2408 integer :: jlon, jlev, jlev2, jlat, jvar, ilev1, ilev2
2409 real(8), pointer :: field_r8(:,:,:)
2410 real(4), pointer :: field_r4(:,:,:)
2411
2412 do jvar = 1, vnl_numvarmax
2413 if(gsv_varExist(statevector,vnl_varNameList(jvar))) then
2414 if (gsv_getDataKind(statevector) == 8) then
2415 call gsv_getField(statevector,field_r8,vnl_varNameList(jvar))
2416 else
2417 call gsv_getField(statevector,field_r4,vnl_varNameList(jvar))
2418 end if
2419 if(vnl_varNameList(jvar) == 'UU ') then
2420 ilev1 = nspositVO
2421 elseif(vnl_varNameList(jvar) == 'VV ') then
2422 ilev1 = nspositDI
2423 elseif(vnl_varNameList(jvar) == 'TT ') then
2424 ilev1 = nspositTT
2425 elseif(vnl_varNameList(jvar) == 'HU ') then
2426 ilev1 = nspositQ
2427 elseif(vnl_varNameList(jvar) == 'P0 ') then
2428 ilev1 = nspositPS
2429 elseif(vnl_varNameList(jvar) == 'TG ') then
2430 ilev1 = nspositTG
2431 else
2432 ! Cycle (instead of abort) to allow for non-NWP assimilation (e.g. chemical data assimilation)
2433 cycle
2434 endif
2435 ilev2 = ilev1 - 1 + gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
2436 if (gsv_getDataKind(statevector) == 8) then
2437 do jlev = ilev1, ilev2
2438 jlev2 = jlev-ilev1+1
2439 do jlat = myLatBeg, myLatEnd
2440 do jlon = myLonBeg, myLonEnd
2441 gd(jlon,jlat,jlev) = field_r8(jlon,jlat,jlev2)
2442 enddo
2443 enddo
2444 enddo
2445 else
2446 do jlev = ilev1, ilev2
2447 jlev2 = jlev-ilev1+1
2448 do jlat = myLatBeg, myLatEnd
2449 do jlon = myLonBeg, myLonEnd
2450 gd(jlon,jlat,jlev) = field_r4(jlon,jlat,jlev2)
2451 enddo
2452 enddo
2453 enddo
2454 end if
2455 endif
2456 enddo
2457
2458 END SUBROUTINE copyFromStatevector
2459
2460
2461 SUBROUTINE BHI_reduceToMPILocal(cv_mpilocal,cv_mpiglobal)
2462 implicit none
2463
2464 ! Arguments:
2465 real(8), intent(out) :: cv_mpilocal(cvDim_mpilocal)
2466 real(8), intent(in) :: cv_mpiglobal(:)
2467
2468 ! Locals:
2469 real(8), allocatable :: cv_allMaxMpilocal(:,:)
2470 integer, allocatable :: cvDim_allMpilocal(:), displs(:)
2471 integer, allocatable :: allnBeg(:),allnEnd(:),allnSkip(:)
2472 integer, allocatable :: allmBeg(:),allmEnd(:),allmSkip(:)
2473 integer :: jproc,cvDim_maxmpilocal,ierr
2474 integer :: jlev,jn,jm,ila_mpiglobal,jdim_mpilocal,jdim_mpiglobal
2475
2476 call rpn_comm_allreduce(cvDim_mpilocal, cvDim_maxmpilocal, &
2477 1,"MPI_INTEGER","MPI_MAX","GRID",ierr)
2478
2479 if(mmpi_myid == 0) then
2480 allocate(cvDim_allMpiLocal(mmpi_nprocs))
2481 else
2482 allocate(cvDim_allMpiLocal(1))
2483 end if
2484
2485 call rpn_comm_gather(cvDim_mpiLocal ,1,"mpi_integer", &
2486 cvDim_allMpiLocal,1,"mpi_integer",0,"GRID",ierr)
2487
2488 if(mmpi_myid == 0) then
2489 allocate(allnBeg(mmpi_nprocs))
2490 allocate(allnEnd(mmpi_nprocs))
2491 allocate(allnSkip(mmpi_nprocs))
2492 allocate(allmBeg(mmpi_nprocs))
2493 allocate(allmEnd(mmpi_nprocs))
2494 allocate(allmSkip(mmpi_nprocs))
2495 else
2496 allocate(allnBeg(1))
2497 allocate(allnEnd(1))
2498 allocate(allnSkip(1))
2499 allocate(allmBeg(1))
2500 allocate(allmEnd(1))
2501 allocate(allmSkip(1))
2502 end if
2503
2504 call rpn_comm_gather(mynBeg ,1,"mpi_integer", &
2505 allnBeg ,1,"mpi_integer",0,"GRID",ierr)
2506 call rpn_comm_gather(mynEnd ,1,"mpi_integer", &
2507 allnEnd ,1,"mpi_integer",0,"GRID",ierr)
2508 call rpn_comm_gather(mynSkip ,1,"mpi_integer", &
2509 allnSkip,1,"mpi_integer",0,"GRID",ierr)
2510
2511 call rpn_comm_gather(mymBeg ,1,"mpi_integer", &
2512 allmBeg ,1,"mpi_integer",0,"GRID",ierr)
2513 call rpn_comm_gather(mymEnd ,1,"mpi_integer", &
2514 allmEnd ,1,"mpi_integer",0,"GRID",ierr)
2515 call rpn_comm_gather(mymSkip ,1,"mpi_integer", &
2516 allmSkip,1,"mpi_integer",0,"GRID",ierr)
2517
2518 ! Prepare to data to be distributed
2519 if (mmpi_myid == 0) then
2520
2521 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
2522
2523 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,jlev,jm,jn,ila_mpiglobal,jdim_mpiglobal)
2524 do jproc = 0, (mmpi_nprocs-1)
2525 cv_allmaxmpilocal(:,jproc+1) = 0.d0
2526 jdim_mpilocal = 0
2527
2528 do jlev = 1, nkgdimSqrt
2529 do jm = allmBeg(jproc+1), allmEnd(jproc+1), allmSkip(jproc+1)
2530 do jn = allnBeg(jproc+1), allnEnd(jproc+1), allnSkip(jproc+1)
2531
2532 if(jm.le.jn) then
2533
2534 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
2535
2536 ! figure out index into global control vector
2537 if(jm == 0) then
2538 ! for jm=0 only real part
2539 jdim_mpiglobal = ila_mpiglobal
2540 else
2541 ! for jm>0 both real and imaginary part
2542 jdim_mpiglobal = 2*ila_mpiglobal-1 - (ntrunc+1)
2543 endif
2544 ! add offset for level
2545 jdim_mpiglobal = jdim_mpiglobal + (jlev-1) * (ntrunc+1)*(ntrunc+1)
2546
2547 ! index into local control vector computer as in cain
2548 if(jm == 0) then
2549 ! only real component for jm=0
2550 jdim_mpilocal = jdim_mpilocal + 1
2551 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
2552 else
2553 ! both real and imaginary components for jm>0
2554 jdim_mpilocal = jdim_mpilocal + 1
2555 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
2556 jdim_mpilocal = jdim_mpilocal + 1
2557 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal+1)
2558 endif
2559
2560 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
2561 write(*,*)
2562 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_mpilocal
2563 write(*,*) ' proc, jlev, jn, jm = ',jproc, jlev, jn, jm
2564 call utl_abort('bhi_reduceToMPILocal')
2565 end if
2566 if (jdim_mpiglobal > cvDim_mpiglobal) then
2567 write(*,*)
2568 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
2569 write(*,*) ' proc, jlev, jn, jm = ',jproc, jlev, jn, jm
2570 call utl_abort('bhi_reduceToMPILocal')
2571 end if
2572
2573 endif
2574 enddo
2575 enddo
2576 enddo
2577
2578 end do ! jproc
2579 !$OMP END PARALLEL DO
2580
2581 else
2582 allocate(cv_allmaxmpilocal(1,1))
2583 end if
2584
2585 !- Distribute
2586 allocate(displs(mmpi_nprocs))
2587 do jproc = 0, (mmpi_nprocs-1)
2588 displs(jproc+1) = jproc*cvDim_maxMpiLocal ! displacement wrt cv_allMaxMpiLocal from which
2589 ! to take the outgoing data to process jproc
2590 end do
2591
2592 call rpn_comm_scatterv(cv_allMaxMpiLocal, cvDim_allMpiLocal, displs, "mpi_double_precision", &
2593 cv_mpiLocal, cvDim_mpiLocal, "mpi_double_precision", &
2594 0, "GRID", ierr)
2595
2596 !- End
2597 deallocate(displs)
2598 deallocate(cv_allMaxMpiLocal)
2599 deallocate(cvDim_allMpiLocal)
2600 deallocate(allnBeg)
2601 deallocate(allnEnd)
2602 deallocate(allnSkip)
2603 deallocate(allmBeg)
2604 deallocate(allmEnd)
2605 deallocate(allmSkip)
2606
2607 END SUBROUTINE BHI_reduceToMPILocal
2608
2609
2610 SUBROUTINE BHI_reduceToMPILocal_r4(cv_mpilocal,cv_mpiglobal)
2611 implicit none
2612
2613 ! Arguments:
2614 real(4), intent(out) :: cv_mpilocal(cvDim_mpilocal)
2615 real(4), intent(in) :: cv_mpiglobal(:)
2616
2617 ! Locals:
2618 real(4), allocatable :: cv_allMaxMpilocal(:,:)
2619 integer, allocatable :: cvDim_allMpilocal(:), displs(:)
2620 integer, allocatable :: allnBeg(:),allnEnd(:),allnSkip(:)
2621 integer, allocatable :: allmBeg(:),allmEnd(:),allmSkip(:)
2622 integer :: jproc,cvDim_maxmpilocal,ierr
2623 integer :: jlev,jn,jm,ila_mpiglobal,jdim_mpilocal,jdim_mpiglobal
2624
2625 call rpn_comm_allreduce(cvDim_mpilocal, cvDim_maxmpilocal, &
2626 1,"MPI_INTEGER","MPI_MAX","GRID",ierr)
2627
2628 if(mmpi_myid == 0) then
2629 allocate(cvDim_allMpiLocal(mmpi_nprocs))
2630 else
2631 allocate(cvDim_allMpiLocal(1))
2632 end if
2633
2634 call rpn_comm_gather(cvDim_mpiLocal ,1,"mpi_integer", &
2635 cvDim_allMpiLocal,1,"mpi_integer",0,"GRID",ierr)
2636
2637 if(mmpi_myid == 0) then
2638 allocate(allnBeg(mmpi_nprocs))
2639 allocate(allnEnd(mmpi_nprocs))
2640 allocate(allnSkip(mmpi_nprocs))
2641 allocate(allmBeg(mmpi_nprocs))
2642 allocate(allmEnd(mmpi_nprocs))
2643 allocate(allmSkip(mmpi_nprocs))
2644 else
2645 allocate(allnBeg(1))
2646 allocate(allnEnd(1))
2647 allocate(allnSkip(1))
2648 allocate(allmBeg(1))
2649 allocate(allmEnd(1))
2650 allocate(allmSkip(1))
2651 end if
2652
2653 call rpn_comm_gather(mynBeg ,1,"mpi_integer", &
2654 allnBeg ,1,"mpi_integer",0,"GRID",ierr)
2655 call rpn_comm_gather(mynEnd ,1,"mpi_integer", &
2656 allnEnd ,1,"mpi_integer",0,"GRID",ierr)
2657 call rpn_comm_gather(mynSkip ,1,"mpi_integer", &
2658 allnSkip,1,"mpi_integer",0,"GRID",ierr)
2659
2660 call rpn_comm_gather(mymBeg ,1,"mpi_integer", &
2661 allmBeg ,1,"mpi_integer",0,"GRID",ierr)
2662 call rpn_comm_gather(mymEnd ,1,"mpi_integer", &
2663 allmEnd ,1,"mpi_integer",0,"GRID",ierr)
2664 call rpn_comm_gather(mymSkip ,1,"mpi_integer", &
2665 allmSkip,1,"mpi_integer",0,"GRID",ierr)
2666
2667 ! Prepare to data to be distributed
2668 if (mmpi_myid == 0) then
2669
2670 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
2671
2672 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,jlev,jm,jn,ila_mpiglobal,jdim_mpiglobal)
2673 do jproc = 0, (mmpi_nprocs-1)
2674 cv_allmaxmpilocal(:,jproc+1) = 0.d0
2675 jdim_mpilocal = 0
2676
2677 do jlev = 1, nkgdimSqrt
2678 do jm = allmBeg(jproc+1), allmEnd(jproc+1), allmSkip(jproc+1)
2679 do jn = allnBeg(jproc+1), allnEnd(jproc+1), allnSkip(jproc+1)
2680
2681 if(jm.le.jn) then
2682
2683 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
2684
2685 ! figure out index into global control vector
2686 if(jm == 0) then
2687 ! for jm=0 only real part
2688 jdim_mpiglobal = ila_mpiglobal
2689 else
2690 ! for jm>0 both real and imaginary part
2691 jdim_mpiglobal = 2*ila_mpiglobal-1 - (ntrunc+1)
2692 endif
2693 ! add offset for level
2694 jdim_mpiglobal = jdim_mpiglobal + (jlev-1) * (ntrunc+1)*(ntrunc+1)
2695
2696 ! index into local control vector computer as in cain
2697 if(jm == 0) then
2698 ! only real component for jm=0
2699 jdim_mpilocal = jdim_mpilocal + 1
2700 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
2701 else
2702 ! both real and imaginary components for jm>0
2703 jdim_mpilocal = jdim_mpilocal + 1
2704 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
2705 jdim_mpilocal = jdim_mpilocal + 1
2706 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal+1)
2707 endif
2708
2709 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
2710 write(*,*)
2711 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_mpilocal
2712 write(*,*) ' proc, jlev, jn, jm = ',jproc, jlev, jn, jm
2713 call utl_abort('bhi_reduceToMPILocal')
2714 end if
2715 if (jdim_mpiglobal > cvDim_mpiglobal) then
2716 write(*,*)
2717 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
2718 write(*,*) ' proc, jlev, jn, jm = ',jproc, jlev, jn, jm
2719 call utl_abort('bhi_reduceToMPILocal')
2720 end if
2721
2722 endif
2723 enddo
2724 enddo
2725 enddo
2726
2727 end do ! jproc
2728 !$OMP END PARALLEL DO
2729
2730 else
2731 allocate(cv_allmaxmpilocal(1,1))
2732 end if
2733
2734 !- Distribute
2735 allocate(displs(mmpi_nprocs))
2736 do jproc = 0, (mmpi_nprocs-1)
2737 displs(jproc+1) = jproc*cvDim_maxMpiLocal ! displacement wrt cv_allMaxMpiLocal from which
2738 ! to take the outgoing data to process jproc
2739 end do
2740
2741 call rpn_comm_scatterv(cv_allMaxMpiLocal, cvDim_allMpiLocal, displs, "mpi_real4", &
2742 cv_mpiLocal, cvDim_mpiLocal, "mpi_real4", &
2743 0, "GRID", ierr)
2744
2745 !- End
2746 deallocate(displs)
2747 deallocate(cv_allMaxMpiLocal)
2748 deallocate(cvDim_allMpiLocal)
2749 deallocate(allnBeg)
2750 deallocate(allnEnd)
2751 deallocate(allnSkip)
2752 deallocate(allmBeg)
2753 deallocate(allmEnd)
2754 deallocate(allmSkip)
2755
2756 END SUBROUTINE BHI_reduceToMPILocal_r4
2757
2758
2759 SUBROUTINE BHI_expandToMPIGlobal(cv_mpilocal,cv_mpiglobal)
2760 implicit none
2761
2762 ! Arguments:
2763 real(8), intent(in) :: cv_mpilocal(cvDim_mpilocal)
2764 real(8), intent(out) :: cv_mpiglobal(:)
2765
2766 ! Locals:
2767 real(8), allocatable :: cv_maxmpilocal(:)
2768 real(8), pointer :: cv_allmaxmpilocal(:,:) => null()
2769 integer, allocatable :: allnBeg(:),allnEnd(:),allnSkip(:)
2770 integer, allocatable :: allmBeg(:),allmEnd(:),allmSkip(:)
2771 integer :: jlev, jn, jm, jproc, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal, ierr, cvDim_maxmpilocal
2772
2773 !
2774 !- 1. Gather all local control vectors onto mpi task 0
2775 !
2776 call rpn_comm_allreduce(cvDim_mpilocal,cvDim_maxmpilocal,1,"mpi_integer","mpi_max","GRID",ierr)
2777
2778 allocate(cv_maxmpilocal(cvDim_maxmpilocal))
2779
2780 nullify(cv_allmaxmpilocal)
2781 if(mmpi_myid == 0) then
2782 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
2783 else
2784 allocate(cv_allmaxmpilocal(1,1))
2785 end if
2786
2787 cv_maxmpilocal(:) = 0.0d0
2788 cv_maxmpilocal(1:cvDim_mpilocal) = cv_mpilocal(1:cvDim_mpilocal)
2789
2790 call rpn_comm_gather(cv_maxmpilocal, cvDim_maxmpilocal, "mpi_double_precision", &
2791 cv_allmaxmpilocal, cvDim_maxmpilocal, "mpi_double_precision", 0, "GRID", ierr )
2792
2793 deallocate(cv_maxmpilocal)
2794
2795 !
2796 !- 2. Reorganize gathered mpilocal control vectors into the mpiglobal control vector
2797 !
2798 if(mmpi_myid == 0) then
2799 allocate(allnBeg(mmpi_nprocs))
2800 allocate(allnEnd(mmpi_nprocs))
2801 allocate(allnSkip(mmpi_nprocs))
2802 allocate(allmBeg(mmpi_nprocs))
2803 allocate(allmEnd(mmpi_nprocs))
2804 allocate(allmSkip(mmpi_nprocs))
2805 else
2806 allocate(allnBeg(1))
2807 allocate(allnEnd(1))
2808 allocate(allnSkip(1))
2809 allocate(allmBeg(1))
2810 allocate(allmEnd(1))
2811 allocate(allmSkip(1))
2812 end if
2813
2814 call rpn_comm_gather(mynBeg ,1,"mpi_integer", &
2815 allnBeg ,1,"mpi_integer",0,"GRID",ierr)
2816 call rpn_comm_gather(mynEnd ,1,"mpi_integer", &
2817 allnEnd ,1,"mpi_integer",0,"GRID",ierr)
2818 call rpn_comm_gather(mynSkip ,1,"mpi_integer", &
2819 allnSkip,1,"mpi_integer",0,"GRID",ierr)
2820
2821 call rpn_comm_gather(mymBeg ,1,"mpi_integer", &
2822 allmBeg ,1,"mpi_integer",0,"GRID",ierr)
2823 call rpn_comm_gather(mymEnd ,1,"mpi_integer", &
2824 allmEnd ,1,"mpi_integer",0,"GRID",ierr)
2825 call rpn_comm_gather(mymSkip ,1,"mpi_integer", &
2826 allmSkip,1,"mpi_integer",0,"GRID",ierr)
2827
2828 if(mmpi_myid == 0) then
2829 cv_mpiglobal(:) = 0.0d0
2830
2831 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,jlev,jm,jn,ila_mpiglobal,jdim_mpiglobal)
2832 do jproc = 0, (mmpi_nprocs-1)
2833 jdim_mpilocal = 0
2834
2835 do jlev = 1, nkgdimSqrt
2836 do jm = allmBeg(jproc+1), allmEnd(jproc+1), allmSkip(jproc+1)
2837 do jn = allnBeg(jproc+1), allnEnd(jproc+1), allnSkip(jproc+1)
2838 if(jm.le.jn) then
2839
2840 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
2841
2842 ! figure out index into global control vector
2843 if(jm == 0) then
2844 ! for jm=0 only real part
2845 jdim_mpiglobal = ila_mpiglobal
2846 else
2847 ! for jm>0 both real and imaginary part
2848 jdim_mpiglobal = 2*ila_mpiglobal-1 - (ntrunc+1)
2849 endif
2850 ! add offset for level
2851 jdim_mpiglobal = jdim_mpiglobal + (jlev-1) * (ntrunc+1)*(ntrunc+1)
2852
2853 ! index into local control vector
2854 if(jm == 0) then
2855 ! only real component for jm=0
2856 jdim_mpilocal = jdim_mpilocal + 1
2857 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2858 else
2859 ! both real and imaginary components for jm>0
2860 jdim_mpilocal = jdim_mpilocal + 1
2861 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2862 jdim_mpilocal = jdim_mpilocal + 1
2863 cv_mpiglobal(jdim_mpiglobal+1) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2864 endif
2865
2866 if(jdim_mpiglobal.gt.cvDim_mpiglobal) &
2867 write(*,*) 'ERROR: jdim,cvDim,mpiglobal=',jdim_mpiglobal,cvDim_mpiglobal,jlev,jn,jm
2868
2869 endif
2870 enddo
2871 enddo
2872 enddo
2873 enddo ! jproc
2874 !$OMP END PARALLEL DO
2875
2876 endif ! myid == 0
2877
2878 deallocate(allnBeg)
2879 deallocate(allnEnd)
2880 deallocate(allnSkip)
2881 deallocate(allmBeg)
2882 deallocate(allmEnd)
2883 deallocate(allmSkip)
2884 deallocate(cv_allmaxmpilocal)
2885
2886 end SUBROUTINE BHI_expandToMPIGlobal
2887
2888
2889 SUBROUTINE BHI_expandToMPIGlobal_r4(cv_mpilocal,cv_mpiglobal)
2890 implicit none
2891
2892 ! Arguments:
2893 real(4), intent(in) :: cv_mpilocal(cvDim_mpilocal)
2894 real(4), intent(out) :: cv_mpiglobal(:)
2895
2896 ! Locals:
2897 real(4), allocatable :: cv_maxmpilocal(:)
2898 real(4), pointer :: cv_allmaxmpilocal(:,:) => null()
2899 integer, allocatable :: allnBeg(:),allnEnd(:),allnSkip(:)
2900 integer, allocatable :: allmBeg(:),allmEnd(:),allmSkip(:)
2901 integer :: jlev, jn, jm, jproc, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal, ierr, cvDim_maxmpilocal
2902
2903 !
2904 !- 1. Gather all local control vectors onto mpi task 0
2905 !
2906 call rpn_comm_allreduce(cvDim_mpilocal,cvDim_maxmpilocal,1,"mpi_integer","mpi_max","GRID",ierr)
2907
2908 allocate(cv_maxmpilocal(cvDim_maxmpilocal))
2909
2910 nullify(cv_allmaxmpilocal)
2911 if(mmpi_myid == 0) then
2912 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
2913 else
2914 allocate(cv_allmaxmpilocal(1,1))
2915 end if
2916
2917 cv_maxmpilocal(:) = 0.0d0
2918 cv_maxmpilocal(1:cvDim_mpilocal) = cv_mpilocal(1:cvDim_mpilocal)
2919
2920 call rpn_comm_gather(cv_maxmpilocal, cvDim_maxmpilocal, "mpi_real4", &
2921 cv_allmaxmpilocal, cvDim_maxmpilocal, "mpi_real4", 0, "GRID", ierr )
2922
2923 deallocate(cv_maxmpilocal)
2924
2925 !
2926 !- 2. Reorganize gathered mpilocal control vectors into the mpiglobal control vector
2927 !
2928 if(mmpi_myid == 0) then
2929 allocate(allnBeg(mmpi_nprocs))
2930 allocate(allnEnd(mmpi_nprocs))
2931 allocate(allnSkip(mmpi_nprocs))
2932 allocate(allmBeg(mmpi_nprocs))
2933 allocate(allmEnd(mmpi_nprocs))
2934 allocate(allmSkip(mmpi_nprocs))
2935 else
2936 allocate(allnBeg(1))
2937 allocate(allnEnd(1))
2938 allocate(allnSkip(1))
2939 allocate(allmBeg(1))
2940 allocate(allmEnd(1))
2941 allocate(allmSkip(1))
2942 end if
2943
2944 call rpn_comm_gather(mynBeg ,1,"mpi_integer", &
2945 allnBeg ,1,"mpi_integer",0,"GRID",ierr)
2946 call rpn_comm_gather(mynEnd ,1,"mpi_integer", &
2947 allnEnd ,1,"mpi_integer",0,"GRID",ierr)
2948 call rpn_comm_gather(mynSkip ,1,"mpi_integer", &
2949 allnSkip,1,"mpi_integer",0,"GRID",ierr)
2950
2951 call rpn_comm_gather(mymBeg ,1,"mpi_integer", &
2952 allmBeg ,1,"mpi_integer",0,"GRID",ierr)
2953 call rpn_comm_gather(mymEnd ,1,"mpi_integer", &
2954 allmEnd ,1,"mpi_integer",0,"GRID",ierr)
2955 call rpn_comm_gather(mymSkip ,1,"mpi_integer", &
2956 allmSkip,1,"mpi_integer",0,"GRID",ierr)
2957
2958 if(mmpi_myid == 0) then
2959 cv_mpiglobal(:) = 0.0d0
2960
2961 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,jlev,jm,jn,ila_mpiglobal,jdim_mpiglobal)
2962 do jproc = 0, (mmpi_nprocs-1)
2963 jdim_mpilocal = 0
2964
2965 do jlev = 1, nkgdimSqrt
2966 do jm = allmBeg(jproc+1), allmEnd(jproc+1), allmSkip(jproc+1)
2967 do jn = allnBeg(jproc+1), allnEnd(jproc+1), allnSkip(jproc+1)
2968 if(jm.le.jn) then
2969
2970 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
2971
2972 ! figure out index into global control vector
2973 if(jm == 0) then
2974 ! for jm=0 only real part
2975 jdim_mpiglobal = ila_mpiglobal
2976 else
2977 ! for jm>0 both real and imaginary part
2978 jdim_mpiglobal = 2*ila_mpiglobal-1 - (ntrunc+1)
2979 endif
2980 ! add offset for level
2981 jdim_mpiglobal = jdim_mpiglobal + (jlev-1) * (ntrunc+1)*(ntrunc+1)
2982
2983 ! index into local control vector
2984 if(jm == 0) then
2985 ! only real component for jm=0
2986 jdim_mpilocal = jdim_mpilocal + 1
2987 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2988 else
2989 ! both real and imaginary components for jm>0
2990 jdim_mpilocal = jdim_mpilocal + 1
2991 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2992 jdim_mpilocal = jdim_mpilocal + 1
2993 cv_mpiglobal(jdim_mpiglobal+1) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
2994 endif
2995
2996 if(jdim_mpiglobal.gt.cvDim_mpiglobal) &
2997 write(*,*) 'ERROR: jdim,cvDim,mpiglobal=',jdim_mpiglobal,cvDim_mpiglobal,jlev,jn,jm
2998
2999 endif
3000 enddo
3001 enddo
3002 enddo
3003 enddo ! jproc
3004 !$OMP END PARALLEL DO
3005
3006 endif ! myid == 0
3007
3008 deallocate(allnBeg)
3009 deallocate(allnEnd)
3010 deallocate(allnSkip)
3011 deallocate(allmBeg)
3012 deallocate(allmEnd)
3013 deallocate(allmSkip)
3014 deallocate(cv_allmaxmpilocal)
3015
3016 end SUBROUTINE BHI_expandToMPIGlobal_r4
3017
3018
3019 SUBROUTINE BHI_cain(controlVector_in,hiControlVector_out)
3020 implicit none
3021
3022 ! Arguments:
3023 real(8), intent(in) :: controlVector_in(cvDim_mpilocal)
3024 real(8), intent(out) :: hiControlVector_out(nla_mpilocal,2,nkgdimSqrt)
3025
3026 ! Locals:
3027 integer :: jdim, jlev, jm, jn, ila_mpilocal, ila_mpiglobal
3028
3029 jdim = 0
3030 hiControlVector_out(:,:,:) = 0.0d0
3031 do jlev = 1, nkgdimSqrt
3032 do jm = mymBeg, mymEnd, mymSkip
3033 do jn = mynBeg, mynEnd, mynSkip
3034 if(jm.le.jn) then
3035 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
3036 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3037 if(jm == 0) then
3038 ! only real component for jm=0
3039 jdim = jdim + 1
3040 hiControlVector_out(ila_mpilocal,1,jlev) = controlVector_in(jdim)
3041 else
3042 ! both real and imaginary components for jm>0
3043 jdim = jdim + 1
3044 hiControlVector_out(ila_mpilocal,1,jlev) = controlVector_in(jdim)
3045 jdim = jdim + 1
3046 hiControlVector_out(ila_mpilocal,2,jlev) = controlVector_in(jdim)
3047 endif
3048 endif
3049 enddo
3050 enddo
3051 enddo
3052
3053 end SUBROUTINE BHI_cain
3054
3055
3056 SUBROUTINE BHI_cainAd(hiControlVector_in,controlVector_out)
3057 IMPLICIT NONE
3058
3059 ! Arguments:
3060 real(8), intent(out) :: controlVector_out(cvDim_mpilocal)
3061 real(8), intent(in) :: hiControlVector_in(nla_mpilocal,2,nkgdimSqrt)
3062
3063 ! Locals:
3064 integer :: jdim, jlev, jm, jn, ila_mpilocal, ila_mpiglobal
3065
3066 jdim = 0
3067 do jlev = 1, nkgdimSqrt
3068 do jm = mymBeg, mymEnd, mymSkip
3069 do jn = mynBeg, mynEnd, mynSkip
3070 if(jm.le.jn) then
3071 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
3072 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3073 if(jm == 0) then
3074 ! only real component for jm=0
3075 jdim = jdim + 1
3076 controlVector_out(jdim) = controlVector_out(jdim) + hiControlVector_in(ila_mpilocal,1,jlev)
3077 else
3078 ! both real and imaginary components for jm>0
3079 jdim = jdim + 1
3080 controlVector_out(jdim) = controlVector_out(jdim) + hiControlVector_in(ila_mpilocal,1,jlev)*2.0d0
3081 jdim = jdim + 1
3082 controlVector_out(jdim) = controlVector_out(jdim) + hiControlVector_in(ila_mpilocal,2,jlev)*2.0d0
3083 endif
3084 endif
3085 enddo
3086 enddo
3087 enddo
3088
3089 END SUBROUTINE BHI_cainAd
3090
3091
3092 SUBROUTINE BHI_SPA2GD(hiControlVector_in,gd_out)
3093 IMPLICIT NONE
3094
3095 ! Arguments:
3096 real(8), intent(in) :: hiControlVector_in(nla_mpilocal,2,nkgdimSqrt)
3097 real(8), intent(out) :: gd_out(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
3098
3099 ! Locals:
3100 real(8) :: sptb(nla_mpilocal,2,nlev_T_even),sp(nla_mpilocal,2,nkgdim)
3101 real(8) :: tb0(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlev_T_even)
3102 integer :: jn,jm,ila_mpilocal,ila_mpiglobal,icount
3103 real(8) :: sq2, zp
3104 real(8) , allocatable :: zsp(:,:,:), zsp2(:,:,:)
3105 integer :: jlev, jlon, jlat, jla_mpilocal, klatPtoT
3106 real(8), pointer :: zgdpsi(:,:,:),zgdchi(:,:,:)
3107 real(8), target :: gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
3108 real(8) :: dla2, dl1sa2, zcoriolis, zpsb(myLonBeg:myLonEnd,myLatBeg:myLatEnd)
3109
3110 klatPtoT = 1
3111
3112 ! maybe not needed:
3113 sp(:,:,:) = 0.0d0
3114 sptb(:,:,:) = 0.0d0
3115
3116 sq2 = sqrt(2.0d0)
3117 allocate(zsp(nkgdimSqrt,2,mymCount))
3118 allocate(zsp2(nkgdim2,2,mymCount))
3119 !$OMP PARALLEL DO PRIVATE(jn,jm,jlev,ila_mpiglobal,ila_mpilocal,zsp2,zsp,icount)
3120 do jn = mynBeg, mynEnd, mynSkip
3121
3122 icount = 0
3123 do jm = mymBeg, mymEnd, mymSkip
3124 if(jm.le.jn) then
3125 icount = icount+1
3126 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
3127 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3128 do jlev = 1, nkgdimSqrt
3129 zsp(jlev,1,icount) = hiControlVector_in(ila_mpilocal,1,jlev)
3130 zsp(jlev,2,icount) = hiControlVector_in(ila_mpilocal,2,jlev)
3131 enddo
3132 endif
3133 enddo
3134
3135 if(icount.gt.0) then
3136
3137 CALL DGEMM('N','N',nkgdim2,2*icount,nkgdimSqrt,1.0d0,corns(1,1,jn),nkgdim2,zsp(1,1,1),nkgdimSqrt,0.0d0,zsp2(1,1,1),nkgdim2)
3138
3139 icount = 0
3140 do jm = mymBeg, mymEnd, mymSkip
3141 if(jm.le.jn) then
3142 icount = icount+1
3143 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
3144 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3145 do jlev = 1, nkgdim
3146 sp(ila_mpilocal,1,jlev) = zsp2(jlev,1,icount)
3147 sp(ila_mpilocal,2,jlev) = zsp2(jlev,2,icount)
3148 enddo
3149 do jlev = 1, nlev_T
3150 sptb(ila_mpilocal,1,jlev) = zsp2(jlev+nkgdim,1,icount)
3151 sptb(ila_mpilocal,2,jlev) = zsp2(jlev+nkgdim,2,icount)
3152 enddo
3153 endif
3154 enddo
3155
3156 endif
3157
3158 ! make adjustments for jm=0
3159 if(mymBeg == 0) then
3160
3161 ila_mpiglobal = gst_getNind(0,gstID) + jn
3162 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3163
3164 do jlev = 1, nkgdim
3165 sp(ila_mpilocal,1,jlev) = sp(ila_mpilocal,1,jlev)*sq2
3166 sp(ila_mpilocal,2,jlev) = 0.0d0
3167 enddo
3168 do jlev = 1, nlev_T
3169 sptb(ila_mpilocal,1,jlev) = sptb(ila_mpilocal,1,jlev)*sq2
3170 sptb(ila_mpilocal,2,jlev) = 0.0d0
3171 enddo
3172
3173 endif
3174
3175 enddo
3176 !$OMP END PARALLEL DO
3177 deallocate(zsp)
3178 deallocate(zsp2)
3179
3180 !$OMP PARALLEL DO PRIVATE(JLAT,JLEV,JLON)
3181 do jlev = 1, nkgdim
3182 do jlat = myLatBeg, myLatEnd
3183 do jlon = myLonBeg, myLonEnd
3184 gd(jlon,jlat,jlev) = 0.0d0
3185 enddo
3186 enddo
3187 enddo
3188 !$OMP END PARALLEL DO
3189
3190 !$OMP PARALLEL DO PRIVATE(JLAT,JLEV,JLON)
3191 do jlev = 1, nlev_T_even
3192 do jlat = myLatBeg, myLatEnd
3193 do jlon = myLonBeg, myLonEnd
3194 tb0(jlon,jlat,jlev) = 0.0d0
3195 enddo
3196 enddo
3197 enddo
3198 !$OMP END PARALLEL DO
3199
3200 call gst_setID(gstID)
3201 call gst_speree(sp,gd)
3202 call gst_setID(gstID2)
3203 call gst_speree(sptb,tb0)
3204
3205 !$OMP PARALLEL DO PRIVATE(jlat,zcoriolis,jlev,jlon,zp)
3206 do jlat = myLatBeg, myLatEnd
3207 zcoriolis = 2.d0*ec_Omega*gst_getRmu(jlat,gstID)
3208 do jlon = myLonBeg, myLonEnd
3209 zpsb(jlon,jlat) = 0.0d0
3210 do jlev = 1, nlevPtoT
3211 zp = zcoriolis*gd(jlon,jlat,nspositVO+jlev-1)
3212 zpsb(jlon,jlat) = zpsb(jlon,jlat) + PtoT(nlev_T+1,jlev,klatPtoT)*zp
3213 enddo
3214 enddo
3215
3216 do jlev = 1, nlev_T
3217 do jlon = myLonBeg, myLonEnd
3218 tb0(jlon,jlat,jlev) = zcoriolis*tb0(jlon,jlat,jlev)
3219 enddo
3220 enddo
3221
3222 do jlev = 1, nkgdim
3223 do jlon = myLonBeg, myLonEnd
3224 if(jlev.ne.nspositTG) then
3225 gd(jlon,jlat,jlev) = gd(jlon,jlat,jlev)*rgsig(jlat,jlev)
3226 else
3227 gd(jlon,jlat,jlev) = gd(jlon,jlat,jlev)*tgstdbg(jlon,jlat)
3228 endif
3229 enddo
3230 enddo
3231
3232 do jlev = 1, nlev_T
3233 do jlon = myLonBeg, myLonEnd
3234 tb0(jlon,jlat,jlev) = tb0(jlon,jlat,jlev)*rgsigtb(jlat,jlev)
3235 gd(jlon,jlat,nspositTT+jlev-1) = gd(jlon,jlat,nspositTT+jlev-1)+tb0(jlon,jlat,jlev)
3236 enddo
3237 enddo
3238 do jlon = myLonBeg, myLonEnd
3239 zpsb(jlon,jlat) = zpsb(jlon,jlat)*rgsigpsb(jlat)
3240 gd(jlon,jlat,nspositPS) = gd(jlon,jlat,nspositPS)+zpsb(jlon,jlat)
3241 enddo
3242 enddo ! jlat
3243 !$OMP END PARALLEL DO
3244
3245 zgdpsi(myLonBeg:,myLatBeg:,1:) => gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nspositVO:(nspositVO+nlev_M-1))
3246 zgdchi(myLonBeg:,myLatBeg:,1:) => gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nspositDI:(nspositDI+nlev_M-1))
3247 !$OMP PARALLEL DO PRIVATE(jlat,jlev,jlon)
3248 do jlev = nlev_bdl, nlev_M
3249 do jlat = myLatBeg, myLatEnd
3250 do jlon = myLonBeg, myLonEnd
3251 zgdchi(jlon,jlat,jlev) = zgdchi(jlon,jlat,jlev) - tantheta(jlev,jlat)*zgdpsi(jlon,jlat,jlev)
3252 enddo
3253 enddo
3254 enddo
3255 !$OMP END PARALLEL DO
3256
3257 sp(:,:,:) = 0.0d0
3258
3259 call gst_setID(gstID)
3260 call gst_reespe(sp,gd)
3261
3262 dla2 = ec_ra * ec_ra
3263 dl1sa2 = 1.d0/dla2
3264 !$OMP PARALLEL DO PRIVATE(JLEV,JLA_MPILOCAL,ILA_MPIGLOBAL)
3265 do jlev = 1, nlev_M
3266 do jla_mpilocal = 1, nla_mpilocal
3267 ila_mpiglobal = ilaList_mpiglobal(jla_mpilocal)
3268 sp(jla_mpilocal,1,nspositVO+jlev-1) = sp(jla_mpilocal,1,nspositVO+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3269 sp(jla_mpilocal,2,nspositVO+jlev-1) = sp(jla_mpilocal,2,nspositVO+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3270 sp(jla_mpilocal,1,nspositDI+jlev-1) = sp(jla_mpilocal,1,nspositDI+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3271 sp(jla_mpilocal,2,nspositDI+jlev-1) = sp(jla_mpilocal,2,nspositDI+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3272 enddo
3273 enddo
3274 !$OMP END PARALLEL DO
3275
3276 !$OMP PARALLEL DO PRIVATE(JLAT,JLEV,JLON)
3277 do jlev = 1, nkgdim
3278 do jlat = myLatBeg, myLatEnd
3279 do jlon = myLonBeg, myLonEnd
3280 gd(jlon,jlat,jlev) = 0.0d0
3281 enddo
3282 enddo
3283 enddo
3284 !$OMP END PARALLEL DO
3285
3286 call gst_setID(gstID)
3287 call gst_spgd(sp,gd,nlev_M)
3288
3289 !$OMP PARALLEL DO PRIVATE(JLAT,JLEV,JLON)
3290 do jlev = 1, nkgdim
3291 do jlat = myLatBeg, myLatEnd
3292 do jlon = myLonBeg, myLonEnd
3293 gd_out(jlon,jlat,jlev) = gd(jlon,jlat,jlev)
3294 enddo
3295 enddo
3296 enddo
3297 !$OMP END PARALLEL DO
3298
3299 END SUBROUTINE BHI_SPA2GD
3300
3301
3302 SUBROUTINE BHI_SPA2GDAD(gd_in,hiControlVector_out)
3303 implicit none
3304
3305 ! Arguments:
3306 real(8), intent(inout) :: hiControlVector_out(nla_mpilocal,2,nkgdimSqrt)
3307 real(8), intent(in) :: gd_in(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
3308
3309 ! Locals:
3310 real(8) :: sptb(nla_mpilocal,2,nlev_T_even)
3311 real(8) :: sp(nla_mpilocal,2,nkgdim)
3312 real(8) :: tb0(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nlev_T_even)
3313 integer :: jn, jm, ila_mpilocal, ila_mpiglobal, icount
3314 real(8) :: sq2, zp
3315 real(8) ,allocatable :: zsp(:,:,:), zsp2(:,:,:)
3316 integer :: jlev, jlon, jlat, jla_mpilocal, klatPtoT
3317 real(8) :: dl1sa2, dla2, zcoriolis, zpsb(myLonBeg:myLonEnd,myLatBeg:myLatEnd)
3318 real(8),pointer :: zgdpsi(:,:,:) ,zgdchi(:,:,:)
3319 real(8), target :: gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nkgdim)
3320
3321 klatPtoT = 1
3322
3323 !$OMP PARALLEL DO PRIVATE(JLAT,JLEV,JLON)
3324 do jlev = 1, nkgdim
3325 do jlat = myLatBeg, myLatEnd
3326 do jlon = myLonBeg, myLonEnd
3327 gd(jlon,jlat,jlev) = gd_in(jlon,jlat,jlev)
3328 enddo
3329 enddo
3330 enddo
3331 !$OMP END PARALLEL DO
3332
3333 call gst_setID(gstID)
3334 call gst_spgda(sp,gd,nlev_M)
3335
3336 dla2 = ec_ra * ec_ra
3337 dl1sa2 = 1.d0/dla2
3338 !$OMP PARALLEL DO PRIVATE(JLEV,JLA_MPILOCAL,ILA_MPIGLOBAL)
3339 do jlev = 1, nlev_M
3340 do jla_mpilocal = 1, nla_mpilocal
3341 ila_mpiglobal = ilaList_mpiglobal(jla_mpilocal)
3342 sp(jla_mpilocal,1,nspositVO+jlev-1) = sp(jla_mpilocal,1,nspositVO+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3343 sp(jla_mpilocal,2,nspositVO+jlev-1) = sp(jla_mpilocal,2,nspositVO+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3344 sp(jla_mpilocal,1,nspositDI+jlev-1) = sp(jla_mpilocal,1,nspositDI+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3345 sp(jla_mpilocal,2,nspositDI+jlev-1) = sp(jla_mpilocal,2,nspositDI+jlev-1)*dl1sa2*gst_getrnnp1(ila_mpiglobal,gstID)
3346 enddo
3347 enddo
3348 !$OMP END PARALLEL DO
3349
3350 call gst_setID(gstID)
3351 call gst_speree(sp,gd)
3352
3353 zgdpsi(myLonBeg:,myLatBeg:,1:) => gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nspositVO:(nspositVO+nlev_M-1))
3354 zgdchi(myLonBeg:,myLatBeg:,1:) => gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nspositDI:(nspositDI+nlev_M-1))
3355 !$OMP PARALLEL DO PRIVATE(jlat,jlev,jlon)
3356 do jlev = nlev_bdl, nlev_M
3357 do jlat = myLatBeg, myLatEnd
3358 do jlon = myLonBeg, myLonEnd
3359 zgdpsi(jlon,jlat,jlev) = zgdpsi(jlon,jlat,jlev) - tantheta(jlev,jlat)*zgdchi(jlon,jlat,jlev)
3360 enddo
3361 enddo
3362 enddo
3363 !$OMP END PARALLEL DO
3364
3365 !$OMP PARALLEL DO PRIVATE(jlat,zcoriolis,jlev,jlon,zp)
3366 do jlat = myLatBeg, myLatEnd
3367 zcoriolis = 2.d0*ec_Omega*gst_getRMU(jlat,gstID)
3368 tb0(:,jlat,:) = 0.0d0
3369 do jlev = 1, nlev_T
3370 do jlon = myLonBeg, myLonEnd
3371 tb0(jlon,jlat,jlev) = gd(jlon,jlat,nspositTT+jlev-1)
3372 tb0(jlon,jlat,jlev) = tb0(jlon,jlat,jlev)*rgsigtb(jlat,jlev)
3373 enddo
3374 enddo
3375 do jlon = myLonBeg, myLonEnd
3376 zpsb(jlon,jlat) = gd(jlon,jlat,nspositPS)
3377 zpsb(jlon,jlat) = zpsb(jlon,jlat)*rgsigpsb(jlat)
3378 enddo
3379
3380 do jlev = 1, nkgdim
3381 do jlon = myLonBeg, myLonEnd
3382 if(jlev.ne.nspositTG) then
3383 gd(jlon,jlat,jlev) = gd(jlon,jlat,jlev)*rgsig(jlat,jlev)
3384 else
3385 gd(jlon,jlat,nspositTG) = gd(jlon,jlat,nspositTG)*tgstdbg(jlon,jlat)
3386 endif
3387 enddo
3388 enddo
3389
3390 do jlev = 1, nlev_T
3391 do jlon = myLonBeg, myLonEnd
3392 tb0(jlon,jlat,jlev) = zcoriolis*tb0(jlon,jlat,jlev)
3393 enddo
3394 enddo
3395
3396 do jlev = 1, nlevPtoT
3397 do jlon = myLonBeg, myLonEnd
3398 zp = PtoT(nlev_T+1,jlev,klatPtoT)*zpsb(jlon,jlat)
3399 gd(jlon,jlat,nspositVO+jlev-1) = zcoriolis*zp+gd(jlon,jlat,nspositVO+jlev-1)
3400 enddo
3401 enddo
3402 enddo
3403 !$OMP END PARALLEL DO
3404
3405 call gst_setID(gstID)
3406 call gst_reespe(sp,gd)
3407 call gst_setID(gstID2)
3408 call gst_reespe(sptb,tb0)
3409
3410 hiControlVector_out(:,:,:) = 0.0d0
3411 sq2 = sqrt(2.0d0)
3412 allocate(zsp(nkgdimSqrt,2,mymCount))
3413 allocate(zsp2(nkgdim2,2,mymCount))
3414 !$OMP PARALLEL DO PRIVATE(JN,JM,JLEV,ILA_MPILOCAL,ILA_MPIGLOBAL,zsp,zsp2,icount)
3415 do jn = mynBeg, mynEnd, mynSkip
3416
3417 icount = 0
3418 do jm = mymBeg, mymEnd, mymSkip
3419 if(jm.le.jn) then
3420 icount = icount+1
3421 ila_mpiglobal = gst_getNind(jm,gstID) + jn - jm
3422 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3423 do jlev = 1, nkgdim
3424 zsp2(jlev,1,icount) = sp(ila_mpilocal,1,jlev)
3425 zsp2(jlev,2,icount) = sp(ila_mpilocal,2,jlev)
3426 enddo
3427 do jlev = 1, nlev_T
3428 zsp2(jlev+nkgdim,1,icount) = sptb(ila_mpilocal,1,jlev)
3429 zsp2(jlev+nkgdim,2,icount) = sptb(ila_mpilocal,2,jlev)
3430 enddo
3431 endif
3432 enddo
3433
3434 if(icount.gt.0) then
3435
3436 CALL DGEMM('T','N',nkgdimSqrt,2*icount,nkgdim2,1.0d0,corns(1,1,jn),nkgdim2,zsp2(1,1,1),nkgdim2,0.0d0,zsp(1,1,1),nkgdimSqrt)
3437
3438 icount = 0
3439 do jm = mymBeg, jn, mymSkip
3440 icount=icount+1
3441 ila_mpiglobal = gst_getNIND(jm,gstID) + jn - jm
3442 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3443 do jlev = 1, nkgdimSqrt
3444 hiControlVector_out(ila_mpilocal,1,jlev) = zsp(jlev,1,icount)
3445 hiControlVector_out(ila_mpilocal,2,jlev) = zsp(jlev,2,icount)
3446 enddo
3447 enddo
3448
3449 endif
3450
3451 ! make adjustments for jm=0
3452 if(mymBeg == 0) then
3453
3454 ila_mpiglobal = gst_getNIND(0,gstID) + jn
3455 ila_mpilocal = ilaList_mpilocal(ila_mpiglobal)
3456
3457 do jlev = 1, nkgdimSqrt
3458 hiControlVector_out(ila_mpilocal,1,jlev) = hiControlVector_out(ila_mpilocal,1,jlev)*sq2
3459 hiControlVector_out(ila_mpilocal,2,jlev) = hiControlVector_out(ila_mpilocal,2,jlev)*sq2
3460 enddo
3461
3462 endif
3463
3464 enddo
3465 !$OMP END PARALLEL DO
3466 deallocate(zsp)
3467 deallocate(zsp2)
3468
3469 END SUBROUTINE BHI_SPA2GDAD
3470
3471
3472 SUBROUTINE BHI_Finalize()
3473 implicit none
3474
3475 if (initialized) then
3476 deallocate(pressureProfile_M)
3477 deallocate(pressureProfile_T)
3478 deallocate(PtoT)
3479 deallocate(tantheta)
3480 deallocate(rgsig)
3481 deallocate(tgstdbg)
3482 deallocate(rgsigtb)
3483 deallocate(rgsigpsb)
3484 deallocate(corns)
3485 deallocate(rstddev)
3486 end if
3487
3488 END SUBROUTINE BHI_Finalize
3489
3490END MODULE bMatrixHI_mod
