1MODULE bMatrixDiff_mod
2 ! MODULE bMatrixDiff_mod (prefix='bdiff' category='2. B and R matrices')
3 !
4 !:Purpose: Performs transformation from control vector to analysis increment
5 ! (and the adjoint transformation) using the background-error
6 ! covariance matrix based on correlations modelled using a
7 ! diffusion operator.
8 !
9 use midasMpi_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use horizontalCoord_mod
13 use verticalCoord_mod
14 use varNameList_mod
15 use physicsFunctions_mod
16 use utilities_mod
17 use diffusion_mod
18 implicit none
19 save
20 private
21
22 ! public procedures
23 public :: bdiff_Setup, bdiff_BSqrt, bdiff_BSqrtAd, bdiff_Finalize
24 public :: bdiff_getScaleFactor, bdiff_reduceToMPILocal
25
26 logical :: initialized = .false.
27 integer :: nj_l, ni_l
28 integer :: cvDim_mpilocal, cvDim_mpiglobal
29
30 integer, allocatable :: diffID(:)
31
32 ! Background-error covariance matrix elements.
33 real(8), allocatable :: stddev(:,:,:)
34
35 integer, parameter :: maxNumVars = 200
36 real(8) :: scaleFactor_sigma(maxNumVars)
37
38 ! read in from the namelist:
39 real(8) :: scaleFactor(maxNumVars) ! scale factor applied to variances
40 character(len=4) :: stddevMode ! can be 'GD2D' or 'HOMO'
41 real(8) :: homogeneous_std(maxNumVars) ! homogeneous standard deviation (when stddevMode is 'HOMO')
42
43 ! Number of incremental variables/fields
44 integer :: numvar2d
45 ! Start position of each field in composite arrays
46 integer, allocatable :: nsposit(:)
47 ! Name list of incremental variables/fields
48 character(len=4), allocatable :: bdiff_varNameList(:)
49
50 integer :: myLatBeg, myLatEnd
51 integer :: myLonBeg, myLonEnd
52 integer :: latPerPE, latPerPEmax, lonPerPE, lonPerPEmax
53
54 integer,external :: get_max_rss
55
56 integer :: nulbgst = 0
57
58CONTAINS
59
60 !--------------------------------------------------------------------------
61 ! bdiff_setup
62 !--------------------------------------------------------------------------
63 subroutine bdiff_setup ( hco_in, vco_in, cvDim_out, mode_opt )
64 !
65 !:Purpose: Setup the diffusion B matrix
66 !
67 implicit none
68
69 ! Arguments:
70 type(struct_hco), pointer, intent(in) :: hco_in
71 type(struct_vco), pointer, intent(in) :: vco_in
72 integer , intent(out) :: cvDim_out
73 character(len=*), optional, intent(in) :: mode_opt
74
75 ! Locals:
76 character(len=15) :: bdiff_mode
77 type(struct_vco), pointer :: vco_anl
78 integer :: nulnam, ierr, fnom, fclos
79 integer :: variableIndex, latIndex, latIndexIgnore
80 real(8) :: maxDistance
81 real(8), allocatable :: distance(:)
82 character(len=*), parameter :: myName = 'bdiff_setup'
83 ! namelist variables
84 real :: corr_len( maxNumVars ) ! Horizontal correlation length scale (km)
85 real :: stab( maxNumVars ) ! Stability criteria (definitely < 0.5)
86 integer :: nsamp(maxNumVars) ! Number of samples in the estimation of the normalization factors by randomization
87 real(8) :: latIgnoreFraction ! Relative zonal grid spacing limit where lats near each numerical pole are ignored
88 logical :: useImplicit(maxNumVars) ! choose to use implicit formulation of diffusion operator (.true.) or explicit version (.false.)
89
90 NAMELIST /NAMBDIFF/ corr_len, stab, nsamp, useImplicit, scaleFactor, stddevMode, homogeneous_std, latIgnoreFraction
91
92 call utl_tmg_start(65,'----B_DIFF_Setup')
93 if(mmpi_myid == 0) write(*,*) myName//': starting'
94 if(mmpi_myid == 0) write(*,*) myName//': Memory Used: ',get_max_rss()/1024,'Mb'
95
96 ! default values for namelist variables
97 corr_len(:) = 10.0
98 stab(:) = 0.2
99 nsamp(:) = 10000
100 useImplicit(:) = .false.
101 scaleFactor(:) = 0.0d0
102 stddevMode = 'GD2D'
103 homogeneous_std(:) = -1.0d0
104 latIgnoreFraction = 1.0d6 ! large value so that nothing is ignored by default
105
106 if ( .not. utl_isNamelistPresent('NAMBDIFF','./flnml') ) then
107 if ( mmpi_myid == 0 ) then
108 write(*,*) 'bdiff_setup: nambdiff is missing in the namelist.'
109 write(*,*) ' The default values will be taken.'
110 end if
111 else
112 nulnam = 0
113 ierr = fnom( nulnam,'./flnml','FTN+SEQ+R/O',0)
114 read( nulnam, nml = nambdiff, iostat = ierr )
115 if ( ierr /= 0 ) call utl_abort( myName//': Error reading namelist')
116 if ( mmpi_myid == 0) write( *, nml = nambdiff )
117 ierr = fclos( nulnam )
118 end if
119
120 if ( sum(scaleFactor(:) ) == 0.0d0 ) then
121 if( mmpi_myid == 0) write(*,*) myName//': scaleFactor=0, skipping rest of setup'
122 cvdim_out = 0
123 call utl_tmg_stop(65)
124 return
125 end if
126
127 if ( present( mode_opt ) ) then
128 if ( trim( mode_opt ) == 'Analysis' .or. trim( mode_opt ) == 'BackgroundCheck' ) then
129 bdiff_mode = trim( mode_opt )
130 if( mmpi_myid == 0 ) write(*,*)
131 if( mmpi_myid == 0 ) write(*,*) myName//': Mode activated = ', trim(bdiff_mode)
132 else
133 write(*,*)
134 write(*,*) myName//'mode = ', trim(mode_opt)
135 call utl_abort( myName//': unknown mode' )
136 end if
137 else
138 bdiff_mode = 'Analysis'
139 if( mmpi_myid == 0 ) write(*,*)
140 if( mmpi_myid == 0 ) write(*,*) myName//': analysis mode activated (by default)'
141 end if
142
143 vco_anl => vco_in
144 if (vco_anl%Vcode /= 5002 .and. vco_anl%Vcode /= 5005 .and. vco_anl%Vcode /= 0 ) then
145 write(*,*) myName//'vco_anl%Vcode = ', vco_anl%Vcode
146 call utl_abort( myName//': unknown vertical coordinate type!')
147 end if
148
149 numvar2d = 0
150
151 allocate( bdiff_varNameList( vnl_numvarmax ) )
152 bdiff_varNameList( : )=''
153 allocate( nsposit( vnl_numvarmax + 1 ) )
154 nsposit(1) = 1
155
156 ! Find the 2D variables (within NAMSTATE namelist)
157
158 if ( gsv_varExist( varName = 'GL ' )) then
159
160 numvar2d = numvar2d + 1
161 nsposit( numvar2d + 1 ) = nsposit( numvar2d ) + 1
162 bdiff_varNameList( numvar2d ) = 'GL '
163
164 end if
165
166 if ( gsv_varExist( varName = 'TM ')) then
167
168 numvar2d = numvar2d + 1
169 nsposit( numvar2d + 1 ) = nsposit( numvar2d ) + 1
170 bdiff_varNameList( numvar2d ) = 'TM '
171
172 end if
173
174 if ( numvar2d == 0) then
175
176 if ( mmpi_myid == 0) then
177 write(*,*) myName//': Bdiff matrix not produced.'
178 write(*,*) myName//': END'
179 end if
180 call utl_tmg_stop(65)
181 cvdim_out = 0
182 return
183
184 else if (mmpi_myid == 0) then
185
186 write(*,*) myName//': number of 2D variables', numvar2d, bdiff_varNameList( 1 : numvar2d )
187
188 end if
189
190 if ( trim(bdiff_mode) == 'BackgroundCheck' ) then
191 cvDim_out = 9999 ! Dummy value > 0 to indicate to the background check (s/r ose_compute_HBHT_ensemble)
192 ! that Diff is used
193 call utl_tmg_stop(65)
194 return
195 end if
196
197 ! Assumes the input 'scalefactor' is a scaling factor of the variances.
198
199 do variableIndex = 1, numvar2d
200 if ( scaleFactor( variableIndex ) > 0.0d0 ) then
201 scaleFactor_sigma( variableIndex ) = sqrt( scaleFactor( variableIndex ) )
202 else
203 scaleFactor_sigma( variableIndex ) = 0.0d0
204 end if
205 end do
206
207 ni_l = hco_in%ni
208 nj_l = hco_in%nj
209
210 ! Compute latIndexIgnore from latIgnoreFraction
211 if (latIgnoreFraction < 1.0) then
212 allocate(distance(nj_l))
213 do latIndex = 1, nj_l
214 distance(latIndex) = &
215 phf_calcDistance(real(hco_in%lat2d_4(ni_l/2,latIndex),8), real(hco_in%lon2d_4((ni_l/2)+1,latIndex),8), &
216 real(hco_in%lat2d_4(ni_l/2,latIndex),8), real(hco_in%lon2d_4((ni_l/2) ,latIndex),8))
217 end do
218 maxDistance = maxval(distance(:))
219 if (mmpi_myid==0) write(*,*) myName//': maxDistance = ', maxDistance
220 latIndexIgnore = 0
221 latLoop: do latIndex = 1, nj_l
222 if (distance(latIndex)/maxDistance > latIgnoreFraction) then
223 if (mmpi_myid==0) then
224 write(*,*) ' latIndex-1, distance, fraction = ', latIndex-1, distance(latIndex-1), distance(latIndex-1)/maxDistance
225 write(*,*) '***latIndex , distance, fraction = ', latIndex , distance(latIndex ), distance(latIndex )/maxDistance
226 write(*,*) ' latIndex+1, distance, fraction = ', latIndex+1, distance(latIndex+1), distance(latIndex+1)/maxDistance
227 end if
228 latIndexIgnore = latIndex
229 exit latLoop
230 end if
231 end do latLoop
232 deallocate(distance)
233 else
234 latIndexIgnore = 0
235 end if
236 if (mmpi_myid==0) write(*,*) myName//': latIndexIgnore = ', latIndexIgnore
237
238 allocate( diffID( numvar2d ) )
239 do variableIndex = 1, numvar2d
240 write(*,*) myName//': setup the diffusion operator for the variable ', bdiff_varNameList( variableIndex )
241 diffID( variableIndex ) = diff_setup ( variableIndex, bdiff_varNameList(1:numvar2d), hco_in, vco_in, corr_len( variableIndex ), &
242 stab( variableIndex ), nsamp( variableIndex ), useImplicit( variableIndex ), &
243 latIndexIgnore )
244 end do
245
246 call mmpi_setup_latbands( nj_l, latPerPE, latPerPEmax, myLatBeg, myLatEnd )
247 call mmpi_setup_lonbands( ni_l, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd )
248
249 ! compute mpilocal control vector size
250 cvDim_mpilocal = lonPerPE * latPerPE * numvar2d
251 cvDim_out = cvDim_mpilocal
252
253 ! also compute mpiglobal control vector dimension
254 call rpn_comm_allreduce( cvDim_mpilocal, cvDim_mpiglobal, 1, "mpi_integer", "mpi_sum", "GRID", ierr )
255
256 allocate(stddev(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d))
257
258 call bdiff_rdstats( hco_in, vco_in )
259
260 if(mmpi_myid == 0) write(*,*) myName//': Memory Used: ',get_max_rss()/1024,'Mb'
261
262 if(mmpi_myid == 0) write(*,*) myName//': END'
263
264 initialized = .true.
265
266 call utl_tmg_stop(65)
267
268 end subroutine bdiff_setup
269
270 !--------------------------------------------------------------------------
271 ! bdiff_getScaleFactor
272 !--------------------------------------------------------------------------
273 subroutine bdiff_getScaleFactor( scaleFactor_out )
274 !
275 !:Purpose: Return the specified scaleFactor.
276 !
277 implicit none
278
279 ! Arguments:
280 real(8), intent(out) :: scaleFactor_out(:)
281
282 ! Locals:
283 integer :: variableIndex
284
285 do variableIndex = 1, numvar2d
286
287 scaleFactor_out( variableIndex ) = scaleFactor( variableIndex )
288
289 end do
290
291 end subroutine bdiff_getScaleFactor
292
293 !--------------------------------------------------------------------------
294 ! bdiff_rdstats
295 !--------------------------------------------------------------------------
296 subroutine bdiff_rdstats( hco_in, vco_in )
297 !
298 !:Purpose: To read background-error stats file.
299 !
300 implicit none
301
302 ! Arguments:
303 type(struct_hco), pointer, intent(in) :: hco_in
304 type(struct_vco), pointer, intent(in) :: vco_in
305
306 ! Locals:
307 integer :: ierr, nmax, fnom, fstouv, fstfrm, fclos
308 integer :: variableIndex
309 logical :: lExists
310 character(len=12) :: bFileName1 = './bgstddev'
311 character(len=8) :: bFileName2 = './bgcov'
312 character(len=*), parameter :: myName = 'bdiff_rdstats'
313
314 write(*,*) myName//': stddevMode is ', stddevMode
315 write(*,*) myName//': Number of 2D variables', numvar2d, bdiff_varNameList( 1 : numvar2d )
316
317 if ( stddevMode == 'GD2D' ) then
318
319 inquire( file = bFileName1, exist = lExists )
320
321 if ( lexists ) then
322 ierr = fnom(nulbgst, bFileName1, 'RND+OLD+R/O', 0)
323 if ( ierr == 0 ) then
324 nmax = fstouv(nulbgst, 'RND+OLD')
325 else
326 call utl_abort( myName//': error opening file '//trim(bFileName1))
327 end if
328 else
329 ! Assume background-error stats in file bgcov.
330 inquire( file = bFileName2, exist = lExists )
331 if ( lexists ) then
332 ierr = fnom( nulbgst, bFileName2, 'RND+OLD+R/O', 0 )
333 if ( ierr == 0 ) then
334 nmax = fstouv(nulbgst, 'RND+OLD')
335 else
336 call utl_abort( myName//': error opening file '//trim(bFileName2) )
337 end if
338 else
339 call utl_abort( myName//': no background error statistics file found!!' )
340 end if
341 end if
342
343 call bdiff_readBGstdField( hco_in, vco_in )
344
345 ierr = fstfrm(nulbgst)
346 ierr = fclos(nulbgst)
347
348 else if ( stddevMode == 'HOMO' ) then
349
350 do variableIndex = 1, numvar2d
351 write(*,*) myName//': stdev = ', homogeneous_std( variableIndex ), ' for variable ', bdiff_varNameList( variableIndex )
352 stddev( :, :, variableIndex) = homogeneous_std( variableIndex )
353 end do
354
355 else
356
357 call utl_abort( myName//': unknown stddevMode: '//trim(stddevMode) )
358
359 end if
360
361 call bdiff_scalestd
362
363 write(*,*) myName//': END '
364
365 end subroutine bdiff_rdstats
366
367 !--------------------------------------------------------------------------
368 ! bdiff_readBGstdField
369 !--------------------------------------------------------------------------
370 subroutine bdiff_readBGstdField( hco_in, vco_in )
371 !
372 !:Purpose: to read 2D background error standard deviation field
373 ! stored on Z, U or G grid and interpolate it to the analysis grid
374 !
375 implicit none
376
377 ! Arguments:
378 type(struct_hco), pointer, intent(in) :: hco_in
379 type(struct_vco), pointer, intent(in) :: vco_in
380
381 ! Locals:
382 integer :: variableIndex, ierr
383 type(struct_gsv) :: statevector
384 real(4), pointer :: field3D_r4_ptr(:,:,:)
385 real(8) :: minStddev, maxStddev
386 character(len=*), parameter :: myName = 'bdiff_readBGstdField'
387
388 write(*,*) myName//': Reading 2D fields from ./bgstddev...'
389 write(*,*) myName//': Number of 2D variables', numvar2d, bdiff_varNameList( 1 : numvar2d )
390
391 call gsv_allocate( statevector, 1, hco_in, vco_in, dateStamp_opt=-1, &
392 dataKind_opt=4, mpi_local_opt=.true., &
393 hInterpolateDegree_opt='LINEAR', &
394 varNames_opt=bdiff_varNameList(1:numvar2d) )
395 call gsv_zero( statevector )
396 call gio_readFromFile(statevector, './bgstddev', 'STDDEV', ' ', unitConversion_opt = .false. )
397
398 do variableIndex = 1, numvar2d
399 call gsv_getField( statevector, field3D_r4_ptr, bdiff_varNameList( variableIndex ) )
400 stddev( :, :, variableIndex ) = dble( field3D_r4_ptr( :, :, 1 ) )
401 if (mmpi_nprocs > 1) then
402 call rpn_comm_allreduce(minval(stddev(:,:,variableIndex)),minStddev,1,'mpi_real8','mpi_min','GRID',ierr)
403 call rpn_comm_allreduce(maxval(stddev(:,:,variableIndex)),maxStddev,1,'mpi_real8','mpi_max','GRID',ierr)
404 else
405 minStddev = minval(stddev(:,:,variableIndex))
406 maxStddev = maxval(stddev(:,:,variableIndex))
407 end if
408 write(*,*) myName//': variable ', bdiff_varNameList( variableIndex ),' min/max: ', &
409 minStddev, maxStddev
410 end do
411
412 call gsv_deallocate( statevector )
413
414 end subroutine bdiff_readBGstdField
415
416 !--------------------------------------------------------------------------
417 ! bdiff_scalestd
418 !--------------------------------------------------------------------------
419 subroutine bdiff_scalestd
420 !
421 !:Purpose: To scale background-error standard-deviation values.
422 !
423 implicit none
424
425 ! Locals:
426 integer :: variableIndex
427 character(len=*), parameter :: myName = 'bdiff_scalestd'
428
429 do variableIndex = 1, numvar2d
430
431 write(*,*) myName//': scaling ', bdiff_varNameList( variableIndex ), ' STD field with the factor ', scaleFactor_sigma( variableIndex )
432
433 stddev( :, :, variableIndex ) = scaleFactor_sigma( variableIndex ) * stddev( : , : , variableIndex )
434
435 end do
436
437 end subroutine bdiff_scalestd
438
439 !--------------------------------------------------------------------------
440 ! bdiff_bSqrt
441 !--------------------------------------------------------------------------
442 subroutine bdiff_bSqrt( controlVector_in, statevector )
443 !
444 !:Purpose: Apply the sqrt of the B matrix
445 !
446 implicit none
447
448 ! Arguments:
449 real(8), intent(in) :: controlVector_in( cvDim_mpilocal )
450 type(struct_gsv), intent(inout) :: statevector
451
452 ! Locals:
453 real(8) :: gd_in( myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
454 real(8) :: gd_out(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
455 integer :: variableIndex
456 character(len=*), parameter :: myName = 'bdiff_bSqrt'
457
458 if( .not. initialized) then
459 if( mmpi_myid == 0 ) write(*,*) myName//': bMatrixDIFF not initialized'
460 return
461 end if
462
463 if(mmpi_myid == 0) write(*,*) myName//': starting'
464
465 call bdiff_cain( controlVector_in, gd_in )
466
467 do variableIndex = 1, numvar2d
468
469 ! Apply square root of the diffusion operator.
470 call diff_Csqrt( diffID( variableIndex ), gd_in( :, :, variableIndex ), gd_out( :, :, variableIndex ) )
471
472 ! Multiply by the diagonal matrix of background error standard deviations.
473 gd_out( :, :, variableIndex ) = gd_out( :, :, variableIndex ) * stddev( :, :, variableIndex )
474
475 end do
476
477 call bdiff_copyToStatevector( statevector, gd_out )
478
479 if(mmpi_myid == 0) write(*,*) myName//': Memory Used: ',get_max_rss()/1024,'Mb'
480 if(mmpi_myid == 0) write(*,*) myName//': done'
481
482 end subroutine bdiff_bSqrt
483
484 !--------------------------------------------------------------------------
485 ! bdiff_bSqrtAd
486 !--------------------------------------------------------------------------
487 subroutine bdiff_bSqrtAd( statevector, controlVector_out )
488 !
489 !:Purpose: Apply the adjoint (i.e. transpose) of the sqrt of the B matrix
490 !
491 implicit none
492
493 ! Arguments:
494 type(struct_gsv), intent(in) :: statevector
495 real(8), intent(out) :: controlVector_out(cvDim_mpilocal)
496
497 ! Locals:
498 real(8) :: gd_in( myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
499 real(8) :: gd_out(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
500 integer :: variableIndex
501 character(len=*), parameter :: myName = 'bdiff_bSqrtAd'
502
503 if ( .not. initialized ) then
504 if ( mmpi_myid == 0 ) write(*,*) myName//': bMatrixDIFF not initialized'
505 return
506 end if
507
508 if(mmpi_myid == 0) write(*,*) myName//': starting'
509
510 call bdiff_copyFromStatevector( statevector, gd_in )
511
512 do variableIndex = 1, numvar2d
513
514 ! Multiply by the diagonal matrix of background-error standard deviations.
515 gd_in( :, :, variableIndex ) = gd_in( :, :, variableIndex ) * stddev( :, :, variableIndex )
516
517 ! Apply the adjoint of the square root of the diffusion operator.
518 call diff_Csqrtadj( diffID( variableIndex ), gd_in( :, :, variableIndex ), gd_out( :, :, variableIndex) )
519
520 end do
521
522 call bdiff_cainad(gd_out, controlVector_out)
523
524 if ( mmpi_myid == 0) write(*,*) myName//': Memory Used: ', get_max_rss()/1024,'Mb'
525 if ( mmpi_myid == 0) write(*,*) myNAme//': done'
526
527 end subroutine bdiff_bSqrtAd
528
529 !--------------------------------------------------------------------------
530 ! bdiff_copyToStatevector
531 !--------------------------------------------------------------------------
532 subroutine bdiff_copyToStatevector ( statevector, gd )
533 !
534 !:Purpose: Copy the working array to a statevector object
535 !
536 implicit none
537
538 ! Arguments:
539 real(8), intent(in) :: gd(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
540 type(struct_gsv), intent(inout) :: statevector
541
542 ! Locals:
543 integer :: jlon, jlev, jlev2, jlat, variableIndex, ilev1, ilev2
544 real(4), pointer :: field_r4(:,:,:)
545 real(8), pointer :: field_r8(:,:,:)
546 character(len=*), parameter :: myName = 'bdiff_copyToStatevector'
547
548 do variableIndex = 1, numvar2d
549
550 ilev1 = nsposit( variableIndex )
551 ilev2 = nsposit( variableIndex + 1 ) - 1
552
553 if ( mmpi_myid == 0) write(*,*) myName//': ',bdiff_varNameList( variableIndex )
554
555 if (gsv_getDataKind(statevector) == 4) then
556 call gsv_getField( statevector, field_r4, bdiff_varNameList( variableIndex ) )
557 do jlev = ilev1, ilev2
558 jlev2 = jlev-ilev1+1
559 do jlat = myLatBeg, myLatEnd
560 do jlon = myLonBeg, myLonEnd
561 field_r4(jlon,jlat,jlev2) = gd(jlon,jlat,jlev)
562 end do
563 end do
564 end do
565 else
566 call gsv_getField( statevector, field_r8, bdiff_varNameList( variableIndex ) )
567 do jlev = ilev1, ilev2
568 jlev2 = jlev-ilev1+1
569 do jlat = myLatBeg, myLatEnd
570 do jlon = myLonBeg, myLonEnd
571 field_r8(jlon,jlat,jlev2) = gd(jlon,jlat,jlev)
572 end do
573 end do
574 end do
575 end if
576
577 end do
578
579 end subroutine bdiff_copyToStatevector
580
581 !--------------------------------------------------------------------------
582 ! bdiff_copyFromStatevector
583 !--------------------------------------------------------------------------
584 subroutine bdiff_copyFromStatevector( statevector, gd )
585 !
586 !:Purpose: Copy the contents of the statevector object to the working array
587 !
588 implicit none
589
590 ! Arguments:
591 type(struct_gsv), intent(in) :: statevector
592 real(8), intent(out) :: gd(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
593
594 ! Locals:
595 integer :: jlon, jlev, jlev2, jlat, variableIndex, ilev1, ilev2
596 real(4), pointer :: field_r4(:,:,:)
597 real(8), pointer :: field_r8(:,:,:)
598 character(len=*), parameter :: myName = 'bdiff_copyFromStatevector'
599
600 do variableIndex = 1, numvar2d
601
602 ilev1 = nsposit( variableIndex )
603 ilev2 = nsposit( variableIndex + 1 ) - 1
604
605 if ( mmpi_myid == 0) write(*,*) myName//': ',bdiff_varNameList( variableIndex )
606 if (gsv_getDataKind(statevector) == 4) then
607 call gsv_getField(statevector, field_r4, bdiff_varNameList( variableIndex ))
608 do jlev = ilev1, ilev2
609 jlev2 = jlev-ilev1+1
610 do jlat = myLatBeg, myLatEnd
611 do jlon = myLonBeg, myLonEnd
612 gd( jlon, jlat, jlev ) = field_r4( jlon, jlat, jlev2 )
613 end do
614 end do
615 end do
616 else
617 call gsv_getField(statevector, field_r8, bdiff_varNameList( variableIndex ))
618 do jlev = ilev1, ilev2
619 jlev2 = jlev-ilev1+1
620 do jlat = myLatBeg, myLatEnd
621 do jlon = myLonBeg, myLonEnd
622 gd( jlon, jlat, jlev ) = field_r8( jlon, jlat, jlev2 )
623 end do
624 end do
625 end do
626 end if
627 end do
628
629 end subroutine bdiff_copyFromStatevector
630
631 !--------------------------------------------------------------------------
632 ! bdiff_reduceToMPILocal
633 !--------------------------------------------------------------------------
634 subroutine bdiff_reduceToMPILocal(cv_mpilocal,cv_mpiglobal)
635 !
636 !:Purpose: Extract the subset of the global control vector needed for local MPI task
637 !
638 implicit none
639
640 ! Arguments:
641 real(8), intent(out) :: cv_mpilocal(:)
642 real(8), intent(in) :: cv_mpiglobal(:)
643
644 ! Locals:
645 integer :: jn, jlat, jlon, jlev, ierr
646 real(8), allocatable :: gd_mpiGlobal(:,:,:)
647
648 allocate(gd_mpiGlobal(ni_l,nj_l,numvar2d))
649 gd_mpiGlobal(:,:,:) = 0.0d0
650
651 jn = 0
652 if (mmpi_myid == 0) then
653 do jlev = 1, numvar2d
654 do jlat = 1, nj_l
655 do jlon = 1, ni_l
656 jn = jn + 1
657 gd_mpiGlobal( jlon, jlat, jlev ) = cv_mpiglobal( jn )
658 end do
659 end do
660 end do
661 end if
662 call rpn_comm_bcast(gd_mpiGlobal, size(gd_mpiGlobal), 'MPI_REAL8', 0, 'GRID', ierr)
663
664 jn = 0
665 do jlev = 1, numvar2d
666 do jlat = myLatBeg, myLatEnd
667 do jlon = myLonBeg, myLonEnd
668 jn = jn + 1
669 cv_mpilocal(jn) = gd_mpiGlobal(jlon,jlat,jlev)
670 end do
671 end do
672 end do
673
674 deallocate(gd_mpiGlobal)
675
676 end subroutine bdiff_reduceToMPILocal
677
678 !--------------------------------------------------------------------------
679 ! bdiff_cain
680 !--------------------------------------------------------------------------
681 subroutine bdiff_cain( controlVector_in, gd_out )
682 !
683 !:Purpose: Transform from control vector to working array
684 !
685 implicit none
686
687 ! Arguments:
688 real(8), intent(in) :: controlVector_in(cvDim_mpilocal)
689 real(8), intent(out) :: gd_out(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
690
691 ! Locals:
692 integer :: jn, jlev, jlon, jlat
693
694 jn = 0
695 do jlev = 1, numvar2d
696 do jlat = myLatBeg, myLatEnd
697 do jlon = myLonBeg, myLonEnd
698 jn = jn + 1
699 gd_out( jlon, jlat, jlev ) = ControlVector_in( jn )
700 end do
701 end do
702 end do
703
704 end subroutine bdiff_cain
705
706 !--------------------------------------------------------------------------
707 ! bdiff_cainAd
708 !--------------------------------------------------------------------------
709 subroutine bdiff_cainAd( gd_in, diffControlVector_out )
710 !
711 !:Purpose: Transform from working array to control vector
712 !
713 implicit none
714
715 ! Arguments:
716 real(8), intent(in) :: gd_in(myLonBeg:myLonEnd, myLatBeg:myLatEnd, numvar2d)
717 real(8), intent(out) :: diffControlVector_out(cvDim_mpilocal)
718
719 ! Locals:
720 integer :: jn, jlev, jlon, jlat
721
722 jn = 0
723 do jlev = 1, numvar2d
724 do jlat = myLatBeg, myLatEnd
725 do jlon = myLonBeg, myLonEnd
726 jn = jn + 1
727 diffControlVector_out(jn) = gd_in(jlon,jlat,jlev)
728 end do
729 end do
730 end do
731
732 end subroutine bdiff_cainAd
733
734 !--------------------------------------------------------------------------
735 ! bdiff_Finalize
736 !--------------------------------------------------------------------------
737 subroutine bdiff_Finalize()
738 !
739 !:Purpose: Deallocate some arrays after we don't need the B matrix anymore.
740 !
741 implicit none
742
743 if ( initialized ) then
744 initialized = .false.
745 deallocate( stddev )
746 deallocate( diffID )
747 deallocate( nsposit )
748 deallocate( bdiff_varNameList )
749 end if
750
751 end subroutine bdiff_Finalize
752
753end module bMatrixDiff_mod
