1module lamBmatrixHI_mod
2 ! MODULE lamBmatrixHI_mod (prefix='lbhi' category='2. B and R matrices')
3 !
4 !:Purpose: Performs transformation from control vector to analysis increment
5 ! using the homogeneous and isotropic background error covariance
6 ! matrix.
7 !
8 use midasMpi_mod
9 use horizontalCoord_mod
10 use verticalCoord_mod
11 use LamSpectralTransform_mod
12 use gridStateVector_mod
13 use lamAnalysisGridTransforms_mod
14 use utilities_mod
15 use gridVariableTransforms_mod
16 use varNameList_mod
17 use interpolation_mod
18 implicit none
19 save
20 private
21
22 ! public procedures
23 public :: lbhi_Setup, lbhi_bSqrt, lbhi_bSqrtAdj, lbhi_Finalize
24 public :: lbhi_expandToMPIglobal, lbhi_expandToMPIglobal_r4, lbhi_reduceToMPIlocal, lbhi_reduceToMPIlocal_r4
25
26 integer, parameter :: cv_model = 1
27 integer, parameter :: cv_bhi = 2
28 type :: lbhi_cv
29 character(len=4) :: NomVar(2)
30 character(len=2) :: GridType ! TT=Thermo, MM=Momentum, NS=Non-staggered
31 integer :: nlev
32 integer :: kDimStart
33 integer :: kDimEnd
34 real(8), allocatable :: GpStdDev(:,:,:)
35 integer, allocatable :: ip1(:)
36 end type lbhi_cv
37
38 integer,parameter :: nMaxControlVar = 10
39 type(lbhi_cv) :: ControlVariable(nMaxControlVar)
40
41 integer :: UWindID = -1
42 integer :: VWindID = -1
43
44 character(len=12) :: WindTransform
45
46 type(struct_hco), pointer :: hco_bhi ! Analysis horizontal grid parameters
47 type(struct_hco), pointer :: hco_bstats ! Grid-point std dev horizontal grid parameters
48
49 type(struct_lst) :: lst_bhi ! Spectral transform Parameters for B
50 type(struct_lst) :: lst_wind ! Spectral transform Parameters for Vort/Div -> Psi/Chi
51
52 real(8), allocatable :: bsqrt (:,:,:) ! B^1/2
53
54 integer :: nControlVariable
55 integer :: trunc
56 integer :: nksdim
57 integer :: nkgdim
58 integer :: cvDim
59 integer :: cvDim_mpiglobal
60
61 integer :: nlev_M
62 integer :: nlev_T
63
64 logical :: regrid
65 logical :: initialized = .false.
66
67 integer :: LatPerPE, LatPerPEmax, myLatBeg, myLatEnd
68 integer :: LonPerPE, LonPerPEmax, myLonBeg, myLonEnd
69
70 real(8) :: scaleFactor(vco_maxNumLevels)
71
72 character(len=8), parameter :: BStatsFilename = './bgcov'
73
74contains
75
76 !--------------------------------------------------------------------------
77 ! lbhi_SETUP
78 !--------------------------------------------------------------------------
79 subroutine lbhi_Setup( hco_anl_in, hco_core_in, vco_anl_in, cvDim_out )
80 implicit none
81
82 ! Arguments:
83 type(struct_hco), pointer, intent(in) :: hco_anl_in
84 type(struct_hco), pointer, intent(in) :: hco_core_in
85 type(struct_vco), pointer, intent(in) :: vco_anl_in
86 integer, intent(out) :: cvDim_out
87
88 ! Locals:
89 integer :: var
90 integer :: iu_bstats = 0
91 integer :: iu_flnml = 0
92 integer :: ier, fnom, fstouv, fstfrm, fclos, levIndex, nLev
93 logical :: FileExist
94 type(struct_vco), pointer :: vco_file => null()
95 integer :: ntrunc
96
97 NAMELIST /NAMBHI/ntrunc,scaleFactor
98
99 write(*,*)
100 write(*,*) 'lbhi_Setup: Starting...'
101
102 !
103 !- 0. Read namelist options
104 !
105 ntrunc = 75 ! default values
106 scaleFactor(:) = 0.0d0 ! default values
107
108 ier = fnom(iu_flnml,'./flnml','FTN+SEQ+R/O',0)
109 write(*,*)
110 write(*,*) 'lbhi_setup: Reading namelist, ier = ',ier
111 read(iu_flnml,nml=nambhi)
112 write(*,nml=nambhi)
113 ier = fclos(iu_flnml)
114
115 nLev = max(max(vco_anl_in%nlev_M,vco_anl_in%nlev_T),1)
116
117 do levIndex = 1, nLev
118 if ( scaleFactor(levIndex) > 0.0d0 ) then
119 scaleFactor(levIndex) = sqrt(scaleFactor(levIndex))
120 else
121 scaleFactor(levIndex) = 0.0d0
122 end if
123 end do
124
125 write(*,*) ' sum(scaleFactor) : ',sum(scaleFactor(1:nLev))
126 if ( sum(scaleFactor(1:nLev)) == 0.0d0 ) then
127 write(*,*) 'lambmatrixHI: scaleFactor=0, skipping rest of setup'
128 cvDim_out = 0
129 return
130 end if
131
132 !- Setup the LAM analysis grid metrics
133 call lgt_SetupFromHCO( hco_anl_in, hco_core_in ) ! IN
134
135 trunc = ntrunc
136 write(*,*)
137 write(*,*) 'Spectral TRUNCATION = ', trunc
138
139 !
140 !- 1. Open the background stats file
141 !
142 inquire(file=trim(BStatsFilename), exist=FileExist)
143
144 if ( FileExist ) then
145 ier = fnom(iu_bstats,trim(BStatsFilename),'RND+OLD+R/O',0)
146 if ( ier == 0 ) then
147 write(*,*)
148 write(*,*) 'Background Stats File :', trim(BStatsFilename)
149 write(*,*) 'opened as unit file ',iu_bstats
150 ier = fstouv(iu_bstats,'RND+OLD')
151 else
152 write(*,*)
153 write(*,*) 'lbhi_Setup: Error in opening the background stats file'
154 write(*,*) trim(BStatsFilename)
155 call utl_abort('lbhi_Setup')
156 end if
157 else
158 write(*,*)
159 write(*,*) 'lbhi_Setup: The background stats file DOES NOT EXIST'
160 write(*,*) trim(BStatsFilename)
161 call utl_abort('lbhi_Setup')
162 end if
163
164 ! Check if analysisgrid and covariance file have the same vertical levels
165 call vco_SetupFromFile( vco_file, & ! OUT
166 BStatsFilename ) ! IN
167 if (.not. vco_equal(vco_anl_in,vco_file)) then
168 call utl_abort('lamBmatrixHI: vco from analysisgrid and cov file do not match')
169 end if
170
171 !
172 !- 2. Set some variables
173 !
174 hco_bhi => hco_anl_in
175 nlev_M = vco_anl_in%nlev_M
176 nlev_T = vco_anl_in%nlev_T
177
178 !- Read variables and vertical grid info from the background stats file
179 call lbhi_GetControlVariableInfo( iu_bstats ) ! IN
180 call lbhi_GetHorizGridInfo()
181
182 nkgdim = 0
183 do var = 1, nControlVariable
184 allocate( ControlVariable(var)%GpStdDev (1:hco_bhi%ni, 1:hco_bhi%nj, 1:ControlVariable(var)%nlev) )
185 allocate( ControlVariable(var)%ip1 (1:ControlVariable(var)%nlev) )
186
187 if (ControlVariable(var)%GridType == 'TH') then
188 ControlVariable(var)%ip1(:) = vco_anl_in%ip1_T(:)
189 elseif (ControlVariable(var)%GridType == 'MM') then
190 ControlVariable(var)%ip1(:) = vco_anl_in%ip1_M(:)
191 else
192 ControlVariable(var)%ip1(:) = 0
193 end if
194
195 ControlVariable(var)%kDimStart = nkgdim + 1
196 nkgdim = nkgdim + ControlVariable(var)%nlev
197 ControlVariable(var)%kDimEnd = nkgdim
198 end do
199
200 nksdim = nkgdim ! + nlev
201
202 allocate( bsqrt (1:nksdim, 1:nksdim ,0:trunc) )
203
204 !- 2.2 Initialized the LAM spectral transform
205 call mmpi_setup_lonbands(hco_bhi%ni, & ! IN
206 lonPerPE, lonPerPEmax, myLonBeg, myLonEnd ) ! OUT
207
208 call mmpi_setup_latbands(hco_bhi%nj, & ! IN
209 latPerPE, latPerPEmax, myLatBeg, myLatEnd ) ! OUT
210
211 call lst_Setup(lst_bhi, & ! OUT
212 hco_bhi%ni, hco_bhi%nj, & ! IN
213 hco_bhi%dlon, trunc, & ! IN
214 'LatLonMN', maxlevels_opt=nksdim) ! IN
215
216 cvDim = nkgdim * lst_bhi%nla * lst_bhi%nphase
217 cvDim_out = cvDim
218
219 ! also compute mpiglobal control vector dimension
220 call rpn_comm_allreduce(cvDim,cvDim_mpiglobal,1,"mpi_integer","mpi_sum","GRID",ier)
221
222 !- 2.3 Initialized the Wind spectral transform
223 if ( trim(WindTransform) == 'VortDiv' ) then
224 call lst_Setup(lst_wind, & ! OUT
225 hco_bhi%ni, hco_bhi%nj, & ! IN
226 hco_bhi%dlon, trunc, & ! IN
227 'LatLonMN', maxlevels_opt=nlev_M) ! IN
228 end if
229
230 !
231 !- 3. Read info from the background error statistics file
232 !
233 call lbhi_ReadStats( iu_bstats ) ! IN
234
235 !
236 !- 4. Close the background stats files
237 !
238 ier = fstfrm(iu_bstats)
239 ier = fclos (iu_bstats)
240
241 !
242 !- 6. Ending
243 !
244 initialized = .true.
245
246 write(*,*)
247 write(*,*) 'lbhi_Setup: Done!'
248
249 end subroutine lbhi_Setup
250
251 !--------------------------------------------------------------------------
252 ! lbhi_GetControlVariableInfo
253 !--------------------------------------------------------------------------
254 subroutine lbhi_GetControlVariableInfo( iu_bstats )
255 implicit none
256
257 ! Arguments:
258 integer, intent(in) :: iu_bstats
259
260 ! Locals:
261 integer :: key, fstinf, fstlir, fstlir_s
262 integer :: ni, nj, nlev
263 integer :: dateo, nk
264 integer :: ip1, ip2, ip3
265 integer :: var
266 character(len=4 ) :: nomvar
267 character(len=2 ) :: typvar
268 character(len=12) :: etiket
269 character(len=4), allocatable :: ControlModelVarnameList(:)
270 character(len=4), allocatable :: ControlBhiVarnameList (:)
271 character(len=2), allocatable :: ControlVarGridTypeList (:)
272 integer, allocatable :: ControlVarNlevList (:)
273
274 !
275 !- 1. How Many Control Variables do we have?
276 !
277 dateo = -1
278 etiket = 'NLEV'
279 ip1 = -1
280 ip2 = -1
281 ip3 = -1
282 typvar = ' '
283 nomvar = 'CVL'
284
285 key = fstinf( iu_bstats, & ! IN
286 ni, nj, nk, & ! OUT
287 dateo, etiket, ip1, ip2, ip3, typvar, nomvar )! IN
288
289 if (key < 0) then
290 write(*,*)
291 write(*,*) 'lbhi_GetControlVariableInfo: Unable to find variable =',nomvar
292 call utl_abort('lbhi_GetControlVariableInfo')
293 end if
294
295 nControlVariable = ni
296 write(*,*)
297 write(*,*) 'Number of Control Variables found = ', nControlVariable
298
299 allocate(ControlModelVarnameList(nControlVariable))
300 allocate(ControlBhiVarnameList (nControlVariable))
301 allocate(ControlVarGridTypeList (nControlVariable))
302 allocate(ControlVarNlevList (nControlVariable))
303
304 !
305 !- 2. Read the info from the input file
306 !
307 nomvar = 'CVN'
308
309 etiket = 'MODEL'
310 key = fstlir_s(ControlModelVarnameList, & ! OUT
311 iu_bstats, & ! IN
312 ni, nj, nlev, & ! OUT
313 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
314 if (key < 0) then
315 write(*,*)
316 write(*,*) 'lbhi_GetControlVariableInfo: Cannot find variable ', nomvar
317 call utl_abort('lbhi_GetControlVariableInfo')
318 end if
319
320 etiket = 'B_HI'
321 key = fstlir_s(ControlBhiVarnameList, & ! OUT
322 iu_bstats, & ! IN
323 ni, nj, nlev, & ! OUT
324 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
325 if (key < 0) then
326 write(*,*)
327 write(*,*) 'lbhi_GetControlVariableInfo: Cannot find variable ', nomvar
328 call utl_abort('lbhi_GetControlVariableInfo')
329 end if
330
331
332 nomvar = 'CVL'
333
334 etiket = 'NLEV'
335 key = fstlir (ControlVarNlevList, & ! OUT
336 iu_bstats, & ! IN
337 ni, nj, nlev, & ! OUT
338 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
339 if (key < 0) then
340 write(*,*)
341 write(*,*) 'lbhi_GetControlVariableInfo: Cannot find variable ', nomvar
342 call utl_abort('lbhi_GetControlVariableInfo')
343 end if
344
345 etiket = 'LEVTYPE'
346 key = fstlir_s(ControlVarGridTypeList, & ! OUT
347 iu_bstats, & ! IN
348 ni, nj, nlev, & ! OUT
349 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
350 if (key < 0) then
351 write(*,*)
352 write(*,*) 'lbhi_GetControlVariableInfo: Cannot find variable ', nomvar
353 call utl_abort('lbhi_GetControlVariableInfo')
354 end if
355
356 !
357 !- 3. Introduce the info in the ControlVariable structure
358 !
359 do var = 1, nControlVariable
360 ControlVariable(var)%nomvar(cv_model)= trim(ControlModelVarnameList(var))
361 ControlVariable(var)%nomvar(cv_bhi) = trim(ControlBhiVarnameList(var))
362 ControlVariable(var)%nlev = ControlVarNlevList(var)
363 ControlVariable(var)%GridType = trim(ControlVarGridTypeList(var))
364
365 if (trim(ControlVariable(var)%nomvar(cv_model)) == 'UU' ) then
366 UWindID = var
367 else if (trim(ControlVariable(var)%nomvar(cv_model)) == 'VV' ) then
368 VWindID = var
369 end if
370
371 if ( trim(ControlVariable(var)%nomvar(cv_model)) == 'LQ' ) then
372 ControlVariable(var)%nomvar(cv_model) = 'HU' ! PATCH because gridStateVector uses HU
373 end if
374
375 write(*,*)
376 write(*,*) 'nomvar(cv_model) = ', ControlVariable(var)%nomvar(cv_model)
377 write(*,*) 'nomvar(cv_bhi) = ', ControlVariable(var)%nomvar(cv_bhi)
378 write(*,*) 'nlev = ', ControlVariable(var)%nlev
379 write(*,*) 'gridtype = ', ControlVariable(var)%GridType
380 end do
381
382 deallocate(ControlModelVarnameList)
383 deallocate(ControlBhiVarnameList)
384 deallocate(ControlVarGridTypeList)
385 deallocate(ControlVarNlevList)
386
387 !
388 !- 4. Error traps
389 !
390
391 !- Make sure that all the control variables are present in GridStateVector
392 do var = 1, nControlVariable
393 if ( .not. gsv_varExist(varName=ControlVariable(var)%nomvar(cv_model)) ) then
394 write(*,*)
395 write(*,*) 'lbhi_GetControlVariableInfo: The following variable is MISSING in GridStateVector'
396 write(*,*) trim(ControlVariable(var)%nomvar(cv_model))
397 call utl_abort('lbhi_GetControlVariableInfo')
398 end if
399 end do
400
401 !
402 !- 5. Set the type of momentum control variables
403 !
404 if (UWindID /= -1 .and. VWindID /= -1) then
405 if ( trim(ControlVariable(UWindID)%nomvar(cv_model)) == 'UU' .and. &
406 trim(ControlVariable(UWindID)%nomvar(cv_bhi)) == 'PP' .and. &
407 trim(ControlVariable(VWindID)%nomvar(cv_model)) == 'VV' .and. &
408 trim(ControlVariable(VWindID)%nomvar(cv_bhi)) == 'CC' ) then
409 write(*,*)
410 write(*,*) 'Momentum Control Variables = Psi-Chi'
411 WindTransform = 'PsiChi'
412 else if ( trim(ControlVariable(UWindID)%nomvar(cv_model)) == 'UU' .and. &
413 trim(ControlVariable(UWindID)%nomvar(cv_bhi)) == 'QR' .and. &
414 trim(ControlVariable(VWindID)%nomvar(cv_model)) == 'VV' .and. &
415 trim(ControlVariable(VWindID)%nomvar(cv_bhi)) == 'DD' ) then
416 write(*,*)
417 write(*,*) 'Momentum Control Variables = Vort-Div'
418 WindTransform = 'VortDiv'
419 else if ( trim(ControlVariable(UWindID)%nomvar(cv_model)) == 'UU' .and. &
420 trim(ControlVariable(UWindID)%nomvar(cv_bhi)) == 'UU' .and. &
421 trim(ControlVariable(VWindID)%nomvar(cv_model)) == 'VV' .and. &
422 trim(ControlVariable(VWindID)%nomvar(cv_bhi)) == 'VV' ) then
423 write(*,*)
424 write(*,*) 'Momentum Control Variables = U-V '
425 WindTransform = 'UV'
426 else
427 call utl_abort('lbhi_GetControlVariableInfo: Unkown Wind Tranform')
428 end if
429 end if
430
431 end subroutine lbhi_GetControlVariableInfo
432
433 !--------------------------------------------------------------------------
434 ! lbhi_GetHorizGridInfo
435 !--------------------------------------------------------------------------
436 subroutine lbhi_GetHorizGridInfo()
437 implicit none
438
439 ! Locals:
440 character(len=4) :: varName
441
442 !
443 !- 1. Get horizontal grid parameters
444 !
445 varName = ControlVariable(1)%nomvar(cv_bhi)
446
447 call hco_setupFromFile(hco_bstats, BStatsFilename, ' ', 'bstats', & ! IN
448 varName_opt=varName) ! IN
449
450 !- 1.3 Regridding needed ?
451 if ( hco_equal(hco_bstats,hco_bhi) ) then
452 Regrid = .false.
453 write(*,*)
454 write(*,*) 'lbhi_GetHorizGridInfo: No Horizontal regridding needed'
455 else
456 Regrid = .true.
457 write(*,*)
458 write(*,*) 'lbhi_GetHorizGridInfo: Horizontal regridding is needed'
459 end if
460
461 end subroutine lbhi_GetHorizGridInfo
462
463 !--------------------------------------------------------------------------
464 ! lbhi_READSTATS
465 !--------------------------------------------------------------------------
466 subroutine lbhi_ReadStats( iu_bstats )
467 implicit none
468
469 ! Arguments:
470 integer, intent(in) :: iu_bstats
471
472 !
473 !- 1. Read the background error statistics
474 !
475
476 !- 1.1 Verical correlations of control variables in spectral space
477 call lbhi_ReadBSqrt( iu_bstats ) ! IN
478
479 !- 1.2 Mass - Rotational wind statistical linear balance operator
480
481 ! JFC: Pas encore code
482 !if ( usePtoT ) then
483 ! call lbhi_ReadPtoT( iu_bstats, PtoT_Type )
484 !end if
485
486 !- 1.3 Read grid-point standard deviations of control variables
487 call lbhi_ReadGridPointStdDev(iu_bstats ) ! IN
488
489 end subroutine lbhi_ReadStats
490
491 !--------------------------------------------------------------------------
492 ! lbhi_ReadBSqrt
493 !--------------------------------------------------------------------------
494 subroutine lbhi_ReadBSqrt( iu_bstats )
495 implicit none
496
497 ! Arguments:
498 integer, intent(in) :: iu_bstats
499
500 ! Locals:
501 real(8), allocatable :: bsqrt2d (:,:)
502 integer :: key, fstinf, fstinl, totwvnb, infon
503 integer, parameter :: nmax=2000
504 integer :: liste(nmax)
505 integer :: ip1, ip2, ip3
506 integer :: ni_t, nj_t, nlev_t, dateo
507 character(len=4 ) :: nomvar
508 character(len=2 ) :: typvar
509 character(len=12) :: etiket
510
511 dateo = -1
512 etiket = 'B_SQUAREROOT'
513 ip1 = -1
514 ip3 = -1
515 typvar = ' '
516 nomvar = 'ZN'
517
518 !
519 !- 1. Find the truncation in the stats file
520 !
521 ip2 = -1
522
523 key = fstinl(iu_bstats, & ! IN
524 ni_t, nj_t, nlev_t, & ! OUT
525 dateo, etiket, ip1, ip2, ip3, typvar, nomvar, & ! IN
526 liste, infon, & ! OUT
527 nmax ) ! IN
528
529 if (key >= 0) then
530 !- 1.2 Ensure spectral trunctation are the same
531 if ( infon - 1 /= trunc ) then
532 write(*,*)
533 write(*,*) 'lbhi_ReadBSqrt: Truncation here and on stats file different'
534 write(*,*) 'VAR truncation = ', trunc
535 write(*,*) 'Stats file truncation = ', infon-1
536 call utl_abort('lbhi_ReadBSqrt')
537 end if
538 else
539 write(*,*)
540 write(*,*) 'lbhi_ReadBSqrt: Cannot find B square-root ', nomvar
541 call utl_abort('lbhi_ReadBSqrt')
542 end if
543
544 !
545 !- 2. Read B^0.5
546 !
547 allocate( bsqrt2d (1:nksdim, 1:nksdim) )
548
549 do totwvnb = 0, trunc
550
551 ip2 = totwvnb
552
553 !- 2.1 Check if field exists and its dimensions
554 key = fstinf( iu_bstats, & ! IN
555 ni_t, nj_t, nlev_t, & ! OUT
556 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
557
558 if (key >= 0) then
559 !- 2.2 Ensure that the number of vertical levels are compatible
560 if ( ni_t /= nksdim .or. nj_t /= nksdim ) then
561 write(*,*)
562 write(*,*) 'lbhi_ReadBSqrt: BG stat levels inconsitencies'
563 write(*,*) 'for BSQRT: ni_t, nj_t, nksdim =', ni_t, nj_t, nksdim
564 call utl_abort('lbhi_ReadBSqrt')
565 endif
566
567 !- 2.3 Reading
568 key = utl_fstlir( bsqrt2d, & ! OUT
569 iu_bstats, & ! IN
570 ni_t, nj_t, nlev_t, & ! OUT
571 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
572 else
573 write(*,*)
574 write(*,*) 'lbhi_ReadBSqrt: Cannot find BSQRT for totwvnb = ', totwvnb
575 call utl_abort('lbhi_ReadBSqrt')
576 end if
577
578 !- 2.4 Transfer to a 3D array
579 bsqrt(:,:,totwvnb) = bsqrt2d(:,:)
580
581 end do
582
583 deallocate( bsqrt2d )
584
585 end subroutine lbhi_ReadBSqrt
586
587 !--------------------------------------------------------------------------
588 ! lbhi_ReadGridPointStdDev
589 !--------------------------------------------------------------------------
590 subroutine lbhi_ReadGridPointStdDev(iu_bstats)
591 implicit none
592
593 ! Locals:
594 integer, intent(in) :: iu_bstats
595
596 ! Locals:
597 real(8), allocatable :: StdDev2D(:,:)
598 real(8), allocatable :: StdDev2D_Regrid(:,:)
599 integer :: ezdefset, ier
600 integer :: ni_t, nj_t, nlev_t, var, k
601 integer :: dateo, ip1,ip2,ip3
602 character(len=4 ) :: nomvar
603 character(len=2 ) :: typvar
604 character(len=12) :: etiket
605 real(8) :: UnitConv
606
607 !
608 !- 1. Read grid point standard deviations
609 !
610 allocate( StdDev2D(1:hco_bstats%ni,1:hco_bstats%nj) )
611 if (Regrid) then
612 allocate( StdDev2D_Regrid(1:hco_bhi%ni, 1:hco_bhi%nj) )
613 ier = ezdefset(hco_bhi%EZscintID, hco_bstats%EZscintID) ! IN
614 end if
615
616 !- 1.1 Loop over Control Variables
617 do var = 1, nControlVariable
618
619 !- 1.2 Loop over vertical Levels
620 do k = 1, ControlVariable(var)%nlev
621 dateo = -1
622 ip1 = ControlVariable(var)%ip1(k)
623 ip2 = -1
624 ip3 = -1
625 typvar = ' '
626 nomvar = trim(ControlVariable(var)%nomvar(cv_bhi))
627 etiket = 'STDDEV'
628 if ( trim(nomvar) == 'P0') then
629 UnitConv = 100.0d0 ! hPa -> Pa
630 else
631 UnitConv = 1.0d0
632 end if
633
634 !- 1.2.1 Reading
635 ier = utl_fstlir( StdDev2D, & ! OUT
636 iu_bstats, & ! IN
637 ni_t, nj_t, nlev_t, & ! OUT
638 dateo, etiket, ip1, ip2, ip3, typvar,nomvar) ! IN
639
640 if (ier < 0) then
641 write(*,*)
642 write(*,*) 'lbhi_ReadGridPointStdDev: Cannot find Std Deviations'
643 write(*,*) 'nomvar =', trim(ControlVariable(var)%nomvar(cv_bhi))
644 write(*,*) 'etiket =', trim(etiket)
645 write(*,*) 'ip1 =', ControlVariable(var)%ip1(k)
646 call utl_abort('lbhi_ReadGridPointStdDev')
647 end if
648
649 if (ni_t /= hco_bstats%ni .or. nj_t /= hco_bstats%nj) then
650 write(*,*)
651 write(*,*) 'lbhi_ReadGridPointStdDev: Invalid dimensions for ...'
652 write(*,*) 'nomvar =', trim(ControlVariable(var)%nomvar(cv_bhi))
653 write(*,*) 'etiket =', trim(etiket)
654 write(*,*) 'ip1 =', ControlVariable(var)%ip1(k)
655 write(*,*) 'Found ni,nj =', ni_t, nj_t
656 write(*,*) 'Should be =', hco_bstats%ni, hco_bstats%nj
657 call utl_abort('lbhi_ReadGridPointStdDev')
658 end if
659
660 !- 1.2.2 Regrid (if necessary) and transfer to 3D array
661 if ( .not. Regrid) then
662 ControlVariable(var)%GpStdDev(:,:,k) = StdDev2D(:,:)
663 else
664 ! Note: EZSCINT setup was done above
665 ier = int_hInterpScalar(StdDev2D_Regrid, StdDev2D, interpDegree='LINEAR')
666 ControlVariable(var)%GpStdDev(:,:,k) = StdDev2D_Regrid(:,:)
667 end if
668
669 !- 1.3 Scaling
670 ControlVariable(var)%GpStdDev(:,:,k) = ControlVariable(var)%GpStdDev(:,:,k) * &
671 UnitConv * scaleFactor(k)
672
673 end do
674
675 end do
676
677 deallocate( StdDev2D )
678 if (Regrid) then
679 deallocate( StdDev2D_Regrid )
680 end if
681
682 end subroutine lbhi_ReadGridPointStdDev
683
684 !--------------------------------------------------------------------------
685 ! lbhi_bSqrt
686 !--------------------------------------------------------------------------
687 subroutine lbhi_bSqrt(controlVector_in, statevector, stateVectorRef_opt)
688 implicit none
689
690 ! Arguments:
691 real(8), intent(in) :: controlVector_in(cvDim)
692 type(struct_gsv), intent(inout) :: statevector
693 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
694
695 ! Locals:
696 real(8), allocatable :: gd_out(:,:,:)
697 real(8), allocatable :: hiControlVector(:,:,:)
698
699 if ( .not. initialized ) then
700 if(mmpi_myid == 0) write(*,*) 'lbhi_bSqrt: LAM_bMatrixHI not initialized'
701 return
702 endif
703
704 write(*,*)
705 write(*,*) 'lbhi_bSqrt: Starting ...'
706
707 !
708 !- 1. Extract data from the 1D controlVector array
709 !
710 allocate( hiControlVector(lst_bhi%nla, lst_bhi%nphase, nksdim) )
711
712 call lbhi_cain( controlVector_in, & ! IN
713 hiControlVector ) ! OUT
714
715 !
716 !- 2. Move from control variables space to model variables space
717 !
718 allocate( gd_out (myLonBeg:myLonEnd, myLatBeg:myLatEnd, nksdim) )
719
720 call lbhi_cv2gd( hiControlVector, & ! IN
721 gd_out ) ! OUT
722
723 deallocate(hiControlVector)
724
725 !
726 !- 3. Transfer results to statevector structure
727 !
728 call StatevectorInterface( statevector, & ! INOUT
729 gd_out, & ! IN
730 'ToStateVector' ) ! IN
731 deallocate(gd_out)
732
733 !
734 !- 4. Convert LQ_inc to HU_inc
735 !
736 if ( gsv_varExist(varName='HU') ) then
737 call gvt_transform( statevector, & ! INOUT
738 'LQtoHU_tlm', & ! IN
739 stateVectorRef_opt=stateVectorRef_opt ) ! IN
740 end if
741
742 write(*,*)
743 write(*,*) 'lbhi_bSqrt: Done'
744
745 end subroutine lbhi_bSqrt
746
747 !--------------------------------------------------------------------------
748 ! lbhi_bSqrtAdj
749 !--------------------------------------------------------------------------
750 subroutine lbhi_bSqrtAdj(statevector, controlVector_out, stateVectorRef_opt)
751 implicit none
752
753 ! Arguments:
754 real(8), intent(out) :: controlVector_out(cvDim)
755 type(struct_gsv), intent(inout) :: statevector
756 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
757
758 ! Locals:
759 real(8), allocatable :: gd_in(:,:,:)
760 real(8), allocatable :: hiControlVector(:,:,:)
761
762 if ( .not. initialized ) then
763 if(mmpi_myid == 0) write(*,*) 'lbhi_bSqrtAdj: LAM_bMatrixHI not initialized'
764 return
765 endif
766
767 write(*,*)
768 write(*,*) 'lbhi_bSqrtAdj: Starting ...'
769
770 !
771 !- 4. Convert LQ_inc to HU_inc
772 !
773 if ( gsv_varExist(varName='HU') ) then
774 call gvt_transform( statevector, & ! INOUT
775 'LQtoHU_ad', & ! IN
776 stateVectorRef_opt=stateVectorRef_opt ) ! IN
777 end if
778
779 !
780 !- 3. Extract data from the StateVector
781 !
782 allocate( gd_in(myLonBeg:myLonEnd, myLatBeg:myLatEnd, nksdim) )
783
784 call StatevectorInterface ( statevector, & ! IN
785 gd_in, & ! OUT
786 'FromStateVector' ) ! IN
787
788 !
789 !- 2. Move from model variables space to control variables space
790 !
791 allocate( hiControlVector(lst_bhi%nla, lst_bhi%nphase, nksdim) )
792 hiControlVector(:,:,:) = 0.d0
793
794 call lbhi_cv2gdAdj( hiControlVector, & ! OUT
795 gd_in ) ! IN
796
797 !
798 !- 1. Put data into the 1D controlVector array
799 !
800 controlVector_out(:) = 0.d0
801 call lbhi_cainAdj(controlVector_out, hiControlVector)
802
803 deallocate(gd_in)
804 deallocate(hiControlVector)
805
806 write(*,*)
807 write(*,*) 'lbhi_bSqrtAdj: Done'
808
809 end subroutine lbhi_bSqrtAdj
810
811 !--------------------------------------------------------------------------
812 ! lbhi_cv2gd
813 !--------------------------------------------------------------------------
814 subroutine lbhi_cv2gd(hiControlVector_in, gd_out)
815 implicit none
816
817 ! Arguments:
818 real(8), intent(inout) :: hiControlVector_in(lst_bhi%nla, lst_bhi%nphase, nksdim)
819 real(8), intent(out) :: gd_out(myLonBeg:myLonEnd ,myLatBeg:myLatEnd ,1:nksdim)
820
821 ! Locals:
822 real(8), allocatable :: uphy(:,:,:)
823 real(8), allocatable :: vphy(:,:,:)
824 real(8), allocatable :: psi(:,:,:)
825 real(8), allocatable :: chi(:,:,:)
826 integer :: kstart, kend, var
827 character(len=19) :: kind
828
829 !
830 !- 1. B^1/2 * xi (in spectral space)
831 !
832 call lbhi_bSqrtXi(hiControlVector_in) ! INOUT
833
834 !
835 !- 2. Spectral Space -> Grid Point Space
836 !
837 kind = 'SpectralToGridPoint'
838 call lst_VarTransform(lst_bhi, & ! IN
839 hiControlVector_in, & ! IN
840 gd_out, & ! OUT
841 kind, nksdim ) ! IN
842
843 !
844 !- 3. Multiply by the grid point standard deviations
845 !
846 !$OMP PARALLEL DO PRIVATE(var,kstart,kend)
847 do var = 1, nControlVariable
848 kstart = ControlVariable(var)%kDimStart
849 kend = ControlVariable(var)%kDimEnd
850 gd_out(:,:,kstart:kend) = gd_out(:,:,kstart:kend) * ControlVariable(var)%GpStdDev(myLonBeg:myLonEnd,myLatBeg:myLatEnd,:)
851 end do
852 !$OMP END PARALLEL DO
853
854 !
855 !- 4. Momentum variables transform
856 !
857 if ( UWindID /= -1 .and. VWindID /= -1) then
858 if ( ControlVariable(UWindID)%nlev /= nlev_M .or. &
859 ControlVariable(VWindID)%nlev /= nlev_M ) then
860 call utl_abort('lbhi_cv2gd: Error in Wind related parameters')
861 end if
862
863 if ( trim(WindTransform) /= 'UV') then
864
865 allocate(uphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
866 allocate(vphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
867 allocate(psi (myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
868 allocate(chi (myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
869
870 !- 4.1 Vort / Div -> Psi / Chi
871 if ( trim(WindTransform) == 'VortDiv') then
872
873 ! 4.1.1 Putting vorticity in psi array and divergence in chi array
874 psi(:,:,:) = gd_out(:,:,ControlVariable(UWindID)%kDimStart:ControlVariable(UWindID)%kDimEnd)
875 chi(:,:,:) = gd_out(:,:,ControlVariable(VWindID)%kDimStart:ControlVariable(VWindID)%kDimEnd)
876
877 ! 4.1.2 Vort -> Psi
878 call lst_Laplacian(lst_wind, & ! IN
879 Psi, & ! INOUT
880 'Inverse', nlev_M) ! IN
881
882 ! 4.1.3 Div -> Chi
883 call lst_Laplacian(lst_wind, & ! IN
884 Chi, & ! INOUT
885 'Inverse', nlev_M) ! IN
886
887 end if
888
889 !- 4.2 Psi / Chi -> U-wind / V-wind
890 if ( trim(WindTransform) == 'PsiChi') then
891 psi(:,:,:) = gd_out(:,:,ControlVariable(UWindID)%kDimStart:ControlVariable(UWindID)%kDimEnd)
892 chi(:,:,:) = gd_out(:,:,ControlVariable(VWindID)%kDimStart:ControlVariable(VWindID)%kDimEnd)
893 end if
894
895 ! 4.2.1 Do Transform
896 call lgt_PsiChiToUV(psi, chi, & ! IN
897 uphy, vphy, & ! OUT
898 nlev_M) ! IN
899
900 ! 4.2.2 Insert results in gd_out and deallocate memories
901 gd_out(:,:,1 : nlev_M) = uphy(:,:,:)
902 gd_out(:,:,nlev_M+1:2*nlev_M) = vphy(:,:,:)
903
904 deallocate(chi)
905 deallocate(psi)
906 deallocate(vphy)
907 deallocate(uphy)
908
909 end if
910
911 end if
912
913 end subroutine lbhi_cv2gd
914
915 !--------------------------------------------------------------------------
916 ! lbhi_cv2gdAdj
917 !--------------------------------------------------------------------------
918 subroutine lbhi_cv2gdAdj(hiControlVector_out, gd_in)
919 implicit none
920
921 ! Arguments:
922 real(8), intent(out) :: hiControlVector_out(lst_bhi%nla, lst_bhi%nphase, nksdim)
923 real(8), intent(inout) :: gd_in(myLonBeg:myLonEnd, myLatBeg:myLatEnd ,1:nksdim)
924
925 ! Locals:
926 real(8), allocatable :: uphy(:,:,:)
927 real(8), allocatable :: vphy(:,:,:)
928 real(8), allocatable :: psi(:,:,:)
929 real(8), allocatable :: chi(:,:,:)
930 integer :: kstart, kend, var
931 character(len=19) :: kind
932
933 !
934 !- 4. Momentum variables transform
935 !
936 if ( UWindID /= -1 .and. VWindID /= -1) then
937 if ( ControlVariable(UWindID)%nlev /= nlev_M .or. &
938 ControlVariable(VWindID)%nlev /= nlev_M ) then
939 call utl_abort('lbhi_cv2gdadj: Error in Wind related parameters')
940 end if
941
942 if ( trim(WindTransform) /= 'UV' ) then
943
944 allocate(uphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
945 allocate(vphy(myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
946 allocate(psi (myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
947 allocate(chi (myLonBeg:myLonEnd,myLatBeg:myLatEnd,1:nlev_M))
948
949 !- 4.2 U-wind / V-wind -> Psi / Chi
950
951 ! 4.2.2 Extract winds
952 uphy(:,:,:) = gd_in(:,:,1 : nlev_M)
953 vphy(:,:,:) = gd_in(:,:,nlev_M+1:2*nlev_M)
954
955 ! 4.2.1 Do Transform
956 call lgt_PsiChiToUVAdj( psi, chi, & ! OUT
957 uphy, vphy, & ! IN
958 nlev_M) ! IN
959
960 !- 4.1 Psi / Chi -> Vort / Div (Laplacian and inverse Laplacian are auto-adjoint operators)
961 if ( trim(WindTransform) == 'VortDiv' ) then
962
963 ! 4.1.2 Psi -> Vort
964 call lst_Laplacian(lst_wind, & ! IN
965 Psi, & ! INOUT
966 'Inverse', nlev_M) ! IN
967
968 ! 4.1.3 Chi -> Vort
969 call lst_Laplacian(lst_wind, & ! IN
970 Chi, & ! INOUT
971 'Inverse', nlev_M) ! IN
972
973 end if
974
975 !- Insert results in gd and deallocate moemories
976 gd_in(:,:,ControlVariable(UWindID)%kDimStart:ControlVariable(UWindID)%kDimEnd) = psi(:,:,:)
977 gd_in(:,:,ControlVariable(VWindID)%kDimStart:ControlVariable(VWindID)%kDimEnd) = chi(:,:,:)
978
979 deallocate(chi)
980 deallocate(psi)
981 deallocate(vphy)
982 deallocate(uphy)
983
984 end if
985
986 end if
987
988 !
989 !- 3. Multiply by the grid point standard deviations
990 !
991
992 !$OMP PARALLEL DO PRIVATE(var,kstart,kend)
993 do var = 1, nControlVariable
994 kstart = ControlVariable(var)%kDimStart
995 kend = ControlVariable(var)%kDimEnd
996 gd_in(:,:,kstart:kend) = gd_in(:,:,kstart:kend) * ControlVariable(var)%GpStdDev(myLonBeg:myLonEnd,myLatBeg:myLatEnd,:)
997 end do
998 !$OMP END PARALLEL DO
999
1000 !
1001 !- 2. Grid Point Space -> Spectral Space
1002 !
1003 kind = 'GridPointToSpectral'
1004 call lst_VarTransform(lst_bhi, & ! IN
1005 hiControlVector_out, & ! OUT
1006 gd_in, & ! IN
1007 kind, nksdim ) ! IN
1008
1009 !
1010 !- 1. B^1/2 * xi (in spectral space)
1011 !
1012 call lbhi_bSqrtXi( hiControlVector_out ) ! INOUT
1013
1014 end subroutine lbhi_cv2gdAdj
1015
1016 !--------------------------------------------------------------------------
1017 ! lbhi_bSqrtXi
1018 !--------------------------------------------------------------------------
1019 subroutine lbhi_bSqrtXi(hiControlVector_in)
1020 implicit none
1021
1022 ! Arguments:
1023 real(8), intent(inout) :: hiControlVector_in(lst_bhi%nla, lst_bhi%nphase, nksdim)
1024
1025 ! Locals:
1026 real(8), allocatable :: sp_in (:,:,:)
1027 real(8), allocatable :: sp_out(:,:,:)
1028 integer :: totwvnb, e, k, ila
1029 integer :: m, n, lda, ldb, ldc
1030
1031 !
1032 !- 1. B^1/2 * xi (in spectral space)
1033 !
1034 do totwvnb = 0, trunc
1035
1036 if ( lst_bhi%nePerK(totwvnb) == 0 ) then
1037 cycle
1038 end if
1039
1040 allocate( sp_in (nksdim,lst_bhi%nphase,lst_bhi%nePerK(totwvnb)) )
1041 allocate( sp_out(nksdim,lst_bhi%nphase,lst_bhi%nePerK(totwvnb)) )
1042
1043 !- 1.1 Select spectral elements associated with the total wavenumber
1044 !$OMP PARALLEL DO PRIVATE(e,ila,k)
1045 do e = 1, lst_bhi%nePerK(totwvnb)
1046 ila = lst_bhi%ilaFromEK(e,totwvnb)
1047 do k = 1, nksdim
1048 sp_in(k,1:lst_bhi%nphase,e) = hiControlVector_in(ila,1:lst_bhi%nphase,k)
1049 end do
1050 end do
1051 !$OMP END PARALLEL DO
1052
1053 !- 1.2 Compute bsqrt * sp_in using DGEMM
1054
1055 ! For documentation on dgemm, see: http://www.netlib.org/blas/dgemm.f
1056 ! Matrix A = BSQRT(:,:,totwvnb)
1057 ! Matrix B = SP_IN
1058 ! Matrix C = SP_OUT
1059 m = nksdim
1060 n = lst_bhi%nphase * lst_bhi%nePerK(totwvnb)
1061 k = nksdim
1062 lda = nksdim
1063 ldb = nksdim
1064 ldc = nksdim
1065
1066 call dgemm( 'N', 'N', m, n, k, 1.d0, & ! IN
1067 bsqrt(:,:,totwvnb), lda, sp_in, ldb, 0.d0, & ! IN
1068 sp_out, & ! OUT
1069 ldc ) ! IN
1070
1071 !- 1.3 Replace sp values with output matrix
1072 !$OMP PARALLEL DO PRIVATE(e,ila,k)
1073 do e = 1, lst_bhi%nePerK(totwvnb)
1074 ila = lst_bhi%ilaFromEK(e,totwvnb)
1075 do k = 1, nksdim
1076 hiControlVector_in(ila,1:lst_bhi%nphase,k) = sp_out(k,1:lst_bhi%nphase,e)
1077 end do
1078 end do
1079 !$OMP END PARALLEL DO
1080
1081 deallocate(sp_in)
1082 deallocate(sp_out)
1083
1084 end do ! Total Wavenumber
1085
1086 end subroutine lbhi_bSqrtXi
1087
1088 !--------------------------------------------------------------------------
1089 ! lbhi_cain
1090 !--------------------------------------------------------------------------
1091 SUBROUTINE lbhi_cain(controlVector_in, hiControlVector_out)
1092 implicit none
1093
1094 ! Arguments:
1095 real(8), intent(in) :: controlVector_in(cvDim)
1096 real(8), intent(out) :: hiControlVector_out(lst_bhi%nla,lst_bhi%nphase,nksdim)
1097
1098 ! Locals:
1099 integer :: dim, k, ila, p
1100
1101 dim = 0
1102 hiControlVector_out(:,:,:) = 0.0d0
1103 do k = 1, nksdim
1104 do ila = 1, lst_bhi%nla
1105 do p = 1, lst_bhi%nphase
1106 dim = dim + 1
1107 hiControlVector_out(ila,p,k) = controlVector_in(dim) * lst_bhi%NormFactor(ila,p)
1108 end do
1109 end do
1110 end do
1111
1112 end SUBROUTINE lbhi_cain
1113
1114 !--------------------------------------------------------------------------
1115 ! lbhi_cainAdj
1116 !--------------------------------------------------------------------------
1117 SUBROUTINE lbhi_cainAdj(controlVector_out, hiControlVector_in)
1118 implicit none
1119
1120 ! Arguments:
1121 real(8), intent(out) :: controlVector_out(cvDim)
1122 real(8), intent(in ) :: hiControlVector_in(lst_bhi%nla,lst_bhi%nphase,nksdim)
1123
1124 ! Locals:
1125 integer :: dim, k, ila, p
1126
1127 dim = 0
1128 do k = 1, nksdim
1129 do ila = 1, lst_bhi%nla
1130 do p = 1, lst_bhi%nphase
1131 dim = dim + 1
1132 controlVector_out(dim) = controlVector_out(dim) + &
1133 hiControlVector_in(ila,p,k) * lst_bhi%NormFactorAd(ila,p)
1134 end do
1135 end do
1136 end do
1137
1138 end SUBROUTINE lbhi_cainAdj
1139
1140 !--------------------------------------------------------------------------
1141 ! StatevectorInterface
1142 !--------------------------------------------------------------------------
1143 subroutine StatevectorInterface(statevector, gd, Direction)
1144 implicit none
1145
1146 ! Arguments:
1147 type(struct_gsv), intent(inout) :: statevector
1148 real(8), intent(inout) :: gd(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nksdim)
1149 character(len=*), intent(in) :: Direction
1150
1151 ! Locals:
1152 integer :: var
1153 integer :: kgdStart, kgdEnd, i, j, k, kgd, nlev
1154 real(4), pointer :: field_r4(:,:,:)
1155 real(8), pointer :: field_r8(:,:,:)
1156 character(len=4 ) :: varname
1157 logical :: ToStateVector
1158
1159 select case ( trim(Direction) )
1160 case ('ToStateVector')
1161 ToStateVector = .true.
1162 case ('FromStateVector')
1163 ToStateVector = .false.
1164 case default
1165 write(*,*)
1166 write(*,*) 'StatevectorInterface: Unknown Direction ', trim(Direction)
1167 call utl_abort('StatevectorInterface')
1168 end select
1169
1170 do var = 1, nControlVariable
1171
1172 varname = ControlVariable(var)%nomvar(cv_model)
1173
1174 if (.not. gsv_varExist(statevector,varname) ) then
1175 write(*,*)
1176 write(*,*) 'StatevectorInterface: The following variable is MISSING in GridStateVector'
1177 write(*,*) varname
1178 call utl_abort('StatevectorInterface')
1179 end if
1180
1181 if (gsv_getDataKind(statevector) == 4) then
1182 call gsv_getField(statevector,field_r4,varname)
1183 else
1184 call gsv_getField(statevector,field_r8,varname)
1185 end if
1186
1187 kgdStart = ControlVariable(var)%kDimStart
1188 kgdEnd = ControlVariable(var)%kDimEnd
1189
1190 nlev = gsv_getNumLev(statevector,vnl_varLevelFromVarname(varname))
1191 if ( kgdEnd - kgdStart + 1 /= nlev ) then
1192 write(*,*)
1193 write(*,*) 'StatevectorInterface: Number of vertical level mismatch'
1194 write(*,*) kgdEnd - kgdStart + 1, nlev
1195 call utl_abort('StatevectorInterface')
1196 end if
1197
1198 if ( ToStateVector ) then
1199 if (gsv_getDataKind(statevector) == 4) then
1200 !$OMP PARALLEL DO PRIVATE(j,kgd,k,i)
1201 do kgd = kgdStart, kgdEnd
1202 k = kgd - kgdStart + 1
1203 do j = myLatBeg, myLatEnd
1204 do i = myLonBeg, myLonEnd
1205 field_r4(i,j,k) = gd(i,j,kgd)
1206 end do
1207 end do
1208 end do
1209 !$OMP END PARALLEL DO
1210 else
1211 !$OMP PARALLEL DO PRIVATE(j,kgd,k,i)
1212 do kgd = kgdStart, kgdEnd
1213 k = kgd - kgdStart + 1
1214 do j = myLatBeg, myLatEnd
1215 do i = myLonBeg, myLonEnd
1216 field_r8(i,j,k) = gd(i,j,kgd)
1217 end do
1218 end do
1219 end do
1220 !$OMP END PARALLEL DO
1221 end if
1222 else
1223 if (gsv_getDataKind(statevector) == 4) then
1224 !$OMP PARALLEL DO PRIVATE(j,kgd,k,i)
1225 do kgd = kgdStart, kgdEnd
1226 k = kgd - kgdStart + 1
1227 do j = myLatBeg, myLatEnd
1228 do i = myLonBeg, myLonEnd
1229 gd(i,j,kgd) = field_r4(i,j,k)
1230 end do
1231 end do
1232 end do
1233 !$OMP END PARALLEL DO
1234 else
1235 !$OMP PARALLEL DO PRIVATE(j,kgd,k,i)
1236 do kgd = kgdStart, kgdEnd
1237 k = kgd - kgdStart + 1
1238 do j = myLatBeg, myLatEnd
1239 do i = myLonBeg, myLonEnd
1240 gd(i,j,kgd) = field_r8(i,j,k)
1241 end do
1242 end do
1243 end do
1244 !$OMP END PARALLEL DO
1245 end if
1246 end if
1247 end do
1248
1249 end subroutine StatevectorInterface
1250
1251 !--------------------------------------------------------------------------
1252 ! lbhi_reduceToMPILocal
1253 !--------------------------------------------------------------------------
1254 SUBROUTINE lbhi_reduceToMPILocal(cv_mpilocal,cv_mpiglobal)
1255 implicit none
1256
1257 ! Arguments:
1258 real(8), intent(out) :: cv_mpilocal(cvDim)
1259 real(8), intent(in) :: cv_mpiglobal(:)
1260
1261 ! Locals:
1262 real(8), allocatable :: cv_allmaxmpilocal(:,:)
1263 integer, allocatable :: cvDim_allMpilocal(:), displs(:)
1264 integer, allocatable :: ilaGlobal(:), allnlaLocal(:)
1265 integer, allocatable :: allilaGlobal(:,:)
1266 integer :: k, ila, p, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal
1267 integer :: ier, nlaMax, cvDim_maxmpilocal, jproc
1268
1269 call rpn_comm_allreduce(cvDim, cvDim_maxmpilocal, &
1270 1,"MPI_INTEGER","MPI_MAX","GRID",ier)
1271
1272 allocate(cvDim_allMpiLocal(mmpi_nprocs))
1273
1274 call rpn_comm_allgather(cvDim ,1,"mpi_integer", &
1275 cvDim_allMpiLocal,1,"mpi_integer","GRID",ier)
1276
1277 call rpn_comm_allreduce(lst_bhi%nla,nlaMax,1,"mpi_integer","mpi_max","GRID",ier)
1278
1279 if (mmpi_myid == 0) then
1280 allocate(allnlaLocal(mmpi_nprocs))
1281 allocate(allilaGlobal(nlaMax,mmpi_nprocs))
1282 else
1283 allocate(allnlaLocal(1))
1284 allocate(allilaGlobal(1,1))
1285 end if
1286
1287 allocate(ilaGlobal(nlaMax))
1288 ilaGlobal(:) = -1
1289 ilaGlobal(1:lst_bhi%nla) = lst_bhi%ilaGlobal(:)
1290
1291 call rpn_comm_gather(lst_bhi%nla, 1, "mpi_integer", &
1292 allnlaLocal, 1, "mpi_integer", 0, "GRID", ier)
1293 call rpn_comm_gather(ilaGlobal , nlaMax, "mpi_integer", &
1294 allilaGlobal, nlaMax, "mpi_integer",0 ,"GRID", ier)
1295
1296 deallocate(ilaGlobal)
1297
1298 ! assign part of mpiglobal vector from current mpi process
1299 if (mmpi_myid == 0) then
1300
1301 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
1302
1303 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,k,ila,p,ila_mpiglobal,jdim_mpiglobal)
1304 do jproc = 0, (mmpi_nprocs-1)
1305 cv_allmaxmpilocal(:,jproc+1) = 0.d0
1306
1307 do k = 1, nksdim
1308 do ila = 1, allnlaLocal(jproc+1)
1309 do p = 1, lst_bhi%nphase
1310
1311 jdim_mpilocal = ( (k-1) * allnlaLocal(jproc+1) * lst_bhi%nphase ) + &
1312 ( (ila-1) * lst_bhi%nphase ) + p
1313
1314 ila_mpiglobal = allilaGlobal(ila,jproc+1)
1315 if ( ila_mpiglobal <= 0 ) then
1316 write(*,*) 'lbhi_reduceToMPILocal: invalid ila_mpiglobal index ', ila_mpiglobal
1317 call utl_abort('lbhi_reduceToMPILocal')
1318 end if
1319
1320 jdim_mpiglobal = ( (k-1) * lst_bhi%nlaGlobal * lst_bhi%nphase ) + &
1321 ( (ila_mpiglobal-1) * lst_bhi%nphase ) + p
1322
1323 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
1324 write(*,*)
1325 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_allMpiLocal(jproc+1)
1326 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1327 call utl_abort('lbhi_reduceToMPILocal')
1328 end if
1329 if (jdim_mpiglobal > cvDim_mpiglobal) then
1330 write(*,*)
1331 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
1332 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1333 call utl_abort('lbhi_reduceToMPILocal')
1334 end if
1335
1336 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
1337
1338 end do
1339 end do
1340 end do
1341 end do
1342 !$OMP END PARALLEL DO
1343
1344 else
1345 allocate(cv_allmaxmpilocal(1,1))
1346 end if
1347
1348 deallocate(allnlaLocal)
1349 deallocate(allilaGlobal)
1350
1351 !- Distribute
1352 allocate(displs(mmpi_nprocs))
1353 !$OMP PARALLEL DO PRIVATE(jproc)
1354 do jproc = 0, (mmpi_nprocs-1)
1355 displs(jproc+1) = jproc*cvDim_maxMpiLocal ! displacement wrt cv_allMaxMpiLocal from which
1356 ! to take the outgoing data to process jproc
1357 end do
1358 !$OMP END PARALLEL DO
1359
1360 call rpn_comm_scatterv(cv_allMaxMpiLocal, cvDim_allMpiLocal, displs, "mpi_double_precision", &
1361 cv_mpiLocal , cvDim , "mpi_double_precision", &
1362 0, "GRID", ier)
1363
1364 deallocate(displs)
1365 deallocate(cv_allMaxMpiLocal)
1366 deallocate(cvDim_allMpiLocal)
1367
1368
1369 END SUBROUTINE lbhi_reduceToMPILocal
1370
1371 !--------------------------------------------------------------------------
1372 ! lbhi_reduceToMPILocal_r4
1373 !--------------------------------------------------------------------------
1374 SUBROUTINE lbhi_reduceToMPILocal_r4(cv_mpilocal,cv_mpiglobal)
1375 implicit none
1376
1377 ! Arguments:
1378 real(4), intent(out) :: cv_mpilocal(cvDim)
1379 real(4), intent(in) :: cv_mpiglobal(:)
1380
1381 ! Locals:
1382 real(4), allocatable :: cv_allmaxmpilocal(:,:)
1383 integer, allocatable :: cvDim_allMpilocal(:), displs(:)
1384 integer, allocatable :: ilaGlobal(:), allnlaLocal(:)
1385 integer, allocatable :: allilaGlobal(:,:)
1386 integer :: k, ila, p, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal
1387 integer :: ier, nlaMax, cvDim_maxmpilocal, jproc
1388
1389 call rpn_comm_allreduce(cvDim, cvDim_maxmpilocal, &
1390 1,"MPI_INTEGER","MPI_MAX","GRID",ier)
1391
1392 allocate(cvDim_allMpiLocal(mmpi_nprocs))
1393
1394 call rpn_comm_allgather(cvDim ,1,"mpi_integer", &
1395 cvDim_allMpiLocal,1,"mpi_integer","GRID",ier)
1396
1397 call rpn_comm_allreduce(lst_bhi%nla,nlaMax,1,"mpi_integer","mpi_max","GRID",ier)
1398
1399 if (mmpi_myid == 0) then
1400 allocate(allnlaLocal(mmpi_nprocs))
1401 allocate(allilaGlobal(nlaMax,mmpi_nprocs))
1402 else
1403 allocate(allnlaLocal(1))
1404 allocate(allilaGlobal(1,1))
1405 end if
1406
1407 allocate(ilaGlobal(nlaMax))
1408 ilaGlobal(:) = -1
1409 ilaGlobal(1:lst_bhi%nla) = lst_bhi%ilaGlobal(:)
1410
1411 call rpn_comm_gather(lst_bhi%nla, 1, "mpi_integer", &
1412 allnlaLocal, 1, "mpi_integer", 0, "GRID", ier)
1413 call rpn_comm_gather(ilaGlobal , nlaMax, "mpi_integer", &
1414 allilaGlobal, nlaMax, "mpi_integer",0 ,"GRID", ier)
1415
1416 deallocate(ilaGlobal)
1417
1418 ! assign part of mpiglobal vector from current mpi process
1419 if (mmpi_myid == 0) then
1420
1421 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
1422
1423 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,k,ila,p,ila_mpiglobal,jdim_mpiglobal)
1424 do jproc = 0, (mmpi_nprocs-1)
1425 cv_allmaxmpilocal(:,jproc+1) = 0.d0
1426
1427 do k = 1, nksdim
1428 do ila = 1, allnlaLocal(jproc+1)
1429 do p = 1, lst_bhi%nphase
1430
1431 jdim_mpilocal = ( (k-1) * allnlaLocal(jproc+1) * lst_bhi%nphase ) + &
1432 ( (ila-1) * lst_bhi%nphase ) + p
1433
1434 ila_mpiglobal = allilaGlobal(ila,jproc+1)
1435 if ( ila_mpiglobal <= 0 ) then
1436 write(*,*) 'lbhi_reduceToMPILocal: invalid ila_mpiglobal index ', ila_mpiglobal
1437 call utl_abort('lbhi_reduceToMPILocal')
1438 end if
1439
1440 jdim_mpiglobal = ( (k-1) * lst_bhi%nlaGlobal * lst_bhi%nphase ) + &
1441 ( (ila_mpiglobal-1) * lst_bhi%nphase ) + p
1442
1443 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
1444 write(*,*)
1445 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_allMpiLocal(jproc+1)
1446 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1447 call utl_abort('lbhi_reduceToMPILocal')
1448 end if
1449 if (jdim_mpiglobal > cvDim_mpiglobal) then
1450 write(*,*)
1451 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
1452 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1453 call utl_abort('lbhi_reduceToMPILocal')
1454 end if
1455
1456 cv_allmaxmpilocal(jdim_mpilocal,jproc+1) = cv_mpiglobal(jdim_mpiglobal)
1457
1458 end do
1459 end do
1460 end do
1461 end do
1462 !$OMP END PARALLEL DO
1463
1464 else
1465 allocate(cv_allmaxmpilocal(1,1))
1466 end if
1467
1468 deallocate(allnlaLocal)
1469 deallocate(allilaGlobal)
1470
1471 !- Distribute
1472 allocate(displs(mmpi_nprocs))
1473 !$OMP PARALLEL DO PRIVATE(jproc)
1474 do jproc = 0, (mmpi_nprocs-1)
1475 displs(jproc+1) = jproc*cvDim_maxMpiLocal ! displacement wrt cv_allMaxMpiLocal from which
1476 ! to take the outgoing data to process jproc
1477 end do
1478 !$OMP END PARALLEL DO
1479
1480 call rpn_comm_scatterv(cv_allMaxMpiLocal, cvDim_allMpiLocal, displs, "mpi_real4", &
1481 cv_mpiLocal , cvDim , "mpi_real4", &
1482 0, "GRID", ier)
1483
1484 deallocate(displs)
1485 deallocate(cv_allMaxMpiLocal)
1486 deallocate(cvDim_allMpiLocal)
1487
1488
1489 END SUBROUTINE lbhi_reduceToMPILocal_r4
1490
1491 !--------------------------------------------------------------------------
1492 ! lbhi_expandToMPIGlobal
1493 !--------------------------------------------------------------------------
1494 SUBROUTINE lbhi_expandToMPIGlobal(cv_mpilocal,cv_mpiglobal)
1495 implicit none
1496
1497 ! Arguments:
1498 real(8), intent(in) :: cv_mpilocal(cvDim)
1499 real(8), intent(out) :: cv_mpiglobal(:)
1500
1501 ! Locals:
1502 real(8), allocatable :: cv_maxmpilocal(:)
1503 real(8), pointer :: cv_allmaxmpilocal(:,:) => null()
1504 integer, allocatable :: cvDim_allMpilocal(:)
1505 integer, allocatable :: ilaGlobal(:), allnlaLocal(:)
1506 integer, allocatable :: allilaGlobal(:,:)
1507 integer :: k, ila, p, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal
1508 integer :: ier, nlaMax, cvDim_maxmpilocal, jproc
1509
1510 !
1511 !- 1. Gather all local control vectors onto mpi task 0
1512 !
1513 allocate(cvDim_allMpiLocal(mmpi_nprocs))
1514 call rpn_comm_allgather(cvDim ,1,"mpi_integer", &
1515 cvDim_allMpiLocal,1,"mpi_integer","GRID",ier)
1516
1517 call rpn_comm_allreduce(cvDim,cvDim_maxmpilocal,1,"mpi_integer","mpi_max","GRID",ier)
1518
1519 allocate(cv_maxmpilocal(cvDim_maxmpilocal))
1520
1521 cv_maxmpilocal(:) = 0.0d0
1522 cv_maxmpilocal(1:cvDim) = cv_mpilocal(1:cvDim)
1523
1524 nullify(cv_allmaxmpilocal)
1525 if (mmpi_myid == 0) then
1526 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
1527 else
1528 allocate(cv_allmaxmpilocal(1,1))
1529 end if
1530
1531 call rpn_comm_gather(cv_maxmpilocal, cvDim_maxmpilocal, "mpi_double_precision", &
1532 cv_allmaxmpilocal, cvDim_maxmpilocal, "mpi_double_precision", 0, "GRID", ier )
1533
1534 deallocate(cv_maxmpilocal)
1535
1536 !
1537 !- 2. Reorganize gathered mpilocal control vectors into the mpiglobal control vector
1538 !
1539
1540 call rpn_comm_allreduce(lst_bhi%nla,nlaMax,1,"mpi_integer","mpi_max","GRID",ier)
1541
1542 if (mmpi_myid == 0) then
1543 allocate(allnlaLocal(mmpi_nprocs))
1544 allocate(allilaGlobal(nlaMax,mmpi_nprocs))
1545 else
1546 allocate(allnlaLocal(1))
1547 allocate(allilaGlobal(1,1))
1548 end if
1549
1550 allocate(ilaGlobal(nlaMax))
1551 ilaGlobal(:) = -1
1552 ilaGlobal(1:lst_bhi%nla) = lst_bhi%ilaGlobal(:)
1553
1554 call rpn_comm_gather(lst_bhi%nla, 1, "mpi_integer", &
1555 allnlaLocal, 1, "mpi_integer", 0, "GRID", ier)
1556 call rpn_comm_gather(ilaGlobal , nlaMax, "mpi_integer", &
1557 allilaGlobal, nlaMax, "mpi_integer",0 ,"GRID", ier)
1558
1559 deallocate(ilaGlobal)
1560
1561 if (mmpi_myid == 0) then
1562 cv_mpiglobal(:) = 0.0d0
1563
1564 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,k,ila,p,ila_mpiglobal,jdim_mpiglobal)
1565 do jproc = 0, (mmpi_nprocs-1)
1566 do k = 1, nksdim
1567 do ila = 1, allnlaLocal(jproc+1)
1568 do p = 1, lst_bhi%nphase
1569
1570 jdim_mpilocal = ( (k-1) * allnlaLocal(jproc+1) * lst_bhi%nphase ) + &
1571 ( (ila-1) * lst_bhi%nphase ) + p
1572
1573 ila_mpiglobal = allilaGlobal(ila,jproc+1)
1574 if ( ila_mpiglobal <= 0 ) then
1575 write(*,*) 'lbhi_expandToMPIGlobal: invalid ila_mpiglobal index ', ila_mpiglobal
1576 call utl_abort('lbhi_expandToMPIGlobal')
1577 end if
1578
1579 jdim_mpiglobal = ( (k-1) * lst_bhi%nlaGlobal * lst_bhi%nphase ) + &
1580 ( (ila_mpiglobal-1) * lst_bhi%nphase ) + p
1581
1582 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
1583 write(*,*)
1584 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_allMpiLocal(jproc+1)
1585 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1586 call utl_abort('lbhi_expandToMPIGlobal')
1587 end if
1588 if (jdim_mpiglobal > cvDim_mpiglobal) then
1589 write(*,*)
1590 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
1591 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1592 call utl_abort('lbhi_expandToMPIGlobal')
1593 end if
1594
1595 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
1596
1597 end do
1598 end do
1599 end do
1600 end do
1601 !$OMP END PARALLEL DO
1602
1603 end if
1604
1605 deallocate(allnlaLocal)
1606 deallocate(allilaGlobal)
1607 deallocate(cv_allmaxmpilocal)
1608 deallocate(cvDim_allMpiLocal)
1609
1610 end SUBROUTINE lbhi_expandToMPIGlobal
1611
1612 !--------------------------------------------------------------------------
1613 ! lbhi_expandToMPIGlobal_r4
1614 !--------------------------------------------------------------------------
1615 SUBROUTINE lbhi_expandToMPIGlobal_r4(cv_mpilocal,cv_mpiglobal)
1616 implicit none
1617
1618 ! Arguments:
1619 real(4), intent(in) :: cv_mpilocal(cvDim)
1620 real(4), intent(out) :: cv_mpiglobal(:)
1621
1622 ! Locals:
1623 real(4), allocatable :: cv_maxmpilocal(:)
1624 real(4), pointer :: cv_allmaxmpilocal(:,:) => null()
1625 integer, allocatable :: cvDim_allMpilocal(:)
1626 integer, allocatable :: ilaGlobal(:), allnlaLocal(:)
1627 integer, allocatable :: allilaGlobal(:,:)
1628 integer :: k, ila, p, ila_mpiglobal, jdim_mpilocal, jdim_mpiglobal
1629 integer :: ier, nlaMax, cvDim_maxmpilocal, jproc
1630
1631 !
1632 !- 1. Gather all local control vectors onto mpi task 0
1633 !
1634 allocate(cvDim_allMpiLocal(mmpi_nprocs))
1635 call rpn_comm_allgather(cvDim ,1,"mpi_integer", &
1636 cvDim_allMpiLocal,1,"mpi_integer","GRID",ier)
1637
1638 call rpn_comm_allreduce(cvDim,cvDim_maxmpilocal,1,"mpi_integer","mpi_max","GRID",ier)
1639
1640 allocate(cv_maxmpilocal(cvDim_maxmpilocal))
1641
1642 cv_maxmpilocal(:) = 0.0d0
1643 cv_maxmpilocal(1:cvDim) = cv_mpilocal(1:cvDim)
1644
1645 nullify(cv_allmaxmpilocal)
1646 if (mmpi_myid == 0) then
1647 allocate(cv_allmaxmpilocal(cvDim_maxmpilocal,mmpi_nprocs))
1648 else
1649 allocate(cv_allmaxmpilocal(1,1))
1650 end if
1651
1652 call rpn_comm_gather(cv_maxmpilocal, cvDim_maxmpilocal, "mpi_real4", &
1653 cv_allmaxmpilocal, cvDim_maxmpilocal, "mpi_real4", 0, "GRID", ier )
1654
1655 deallocate(cv_maxmpilocal)
1656
1657 !
1658 !- 2. Reorganize gathered mpilocal control vectors into the mpiglobal control vector
1659 !
1660
1661 call rpn_comm_allreduce(lst_bhi%nla,nlaMax,1,"mpi_integer","mpi_max","GRID",ier)
1662
1663 if (mmpi_myid == 0) then
1664 allocate(allnlaLocal(mmpi_nprocs))
1665 allocate(allilaGlobal(nlaMax,mmpi_nprocs))
1666 else
1667 allocate(allnlaLocal(1))
1668 allocate(allilaGlobal(1,1))
1669 end if
1670
1671 allocate(ilaGlobal(nlaMax))
1672 ilaGlobal(:) = -1
1673 ilaGlobal(1:lst_bhi%nla) = lst_bhi%ilaGlobal(:)
1674
1675 call rpn_comm_gather(lst_bhi%nla, 1, "mpi_integer", &
1676 allnlaLocal, 1, "mpi_integer", 0, "GRID", ier)
1677 call rpn_comm_gather(ilaGlobal , nlaMax, "mpi_integer", &
1678 allilaGlobal, nlaMax, "mpi_integer",0 ,"GRID", ier)
1679
1680 deallocate(ilaGlobal)
1681
1682 if (mmpi_myid == 0) then
1683 cv_mpiglobal(:) = 0.0d0
1684
1685 !$OMP PARALLEL DO PRIVATE(jproc,jdim_mpilocal,k,ila,p,ila_mpiglobal,jdim_mpiglobal)
1686 do jproc = 0, (mmpi_nprocs-1)
1687 do k = 1, nksdim
1688 do ila = 1, allnlaLocal(jproc+1)
1689 do p = 1, lst_bhi%nphase
1690
1691 jdim_mpilocal = ( (k-1) * allnlaLocal(jproc+1) * lst_bhi%nphase ) + &
1692 ( (ila-1) * lst_bhi%nphase ) + p
1693
1694 ila_mpiglobal = allilaGlobal(ila,jproc+1)
1695 if ( ila_mpiglobal <= 0 ) then
1696 write(*,*) 'lbhi_expandToMPIGlobal: invalid ila_mpiglobal index ', ila_mpiglobal
1697 call utl_abort('lbhi_expandToMPIGlobal')
1698 end if
1699
1700 jdim_mpiglobal = ( (k-1) * lst_bhi%nlaGlobal * lst_bhi%nphase ) + &
1701 ( (ila_mpiglobal-1) * lst_bhi%nphase ) + p
1702
1703 if (jdim_mpilocal > cvDim_allMpiLocal(jproc+1)) then
1704 write(*,*)
1705 write(*,*) 'ERROR: jdim_mpilocal > cvDim_allMpiLocal(jproc+1)', jdim_mpilocal, cvDim_allMpiLocal(jproc+1)
1706 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1707 call utl_abort('lbhi_expandToMPIGlobal')
1708 end if
1709 if (jdim_mpiglobal > cvDim_mpiglobal) then
1710 write(*,*)
1711 write(*,*) 'ERROR: jdim_mpiglobal > cvDim_mpiglobal', jdim_mpiglobal, cvDim_mpiglobal
1712 write(*,*) ' proc, k, ila, p = ',jproc,k,ila,p
1713 call utl_abort('lbhi_expandToMPIGlobal')
1714 end if
1715
1716 cv_mpiglobal(jdim_mpiglobal) = cv_allmaxmpilocal(jdim_mpilocal,jproc+1)
1717
1718 end do
1719 end do
1720 end do
1721 end do
1722 !$OMP END PARALLEL DO
1723
1724 end if
1725
1726 deallocate(allnlaLocal)
1727 deallocate(allilaGlobal)
1728 deallocate(cv_allmaxmpilocal)
1729 deallocate(cvDim_allMpiLocal)
1730
1731 end SUBROUTINE lbhi_expandToMPIGlobal_r4
1732
1733 !--------------------------------------------------------------------------
1734 ! lbhi_Finalize
1735 !--------------------------------------------------------------------------
1736 subroutine lbhi_Finalize
1737 implicit none
1738
1739 ! Locals:
1740 integer :: var
1741
1742 if (initialized) then
1743 deallocate(bsqrt)
1744 do var = 1, nControlVariable
1745 deallocate(ControlVariable(var)%GpStdDev)
1746 deallocate(ControlVariable(var)%ip1 )
1747 end do
1748 end if
1749
1750 end subroutine lbhi_Finalize
1751
1752end module lamBmatrixHI_mod
