1module bMatrix1DVar_mod
2 ! MODULE bMatrix1DVar_mod (prefix='bmat1D' category='2. B and R matrices')
3 !
4 !:Purpose: contains all 1Dvar B matrices.
5 !
6 use mathPhysConstants_mod
7 use columnData_mod
8 use columnVariableTransforms_mod
9 use controlVector_mod
10 use gridStatevector_mod
11 use gridstatevectorFileIO_mod
12 use horizontalCoord_mod
13 use midasMpi_mod
14 use obsSpaceData_mod
15 use timeCoord_mod
16 use utilities_mod
17 use verticalCoord_mod
18 use tovsNL_mod
19 use var1D_mod
20 use filenames_mod
21 use localizationFunction_mod
22 use varNameList_mod
23 use ensembleStateVector_mod
24 use stateToColumn_mod
25 use calcHeightAndPressure_mod
26 implicit none
27 save
28 private
29
30 ! public procedures
31 public :: bmat1D_bsetup
32 public :: bmat1D_sqrtB, bmat1D_sqrtBT
33 public :: bmat1D_finalize, bmat1D_get1DVarIncrement
34
35 ! public variables
36 public :: bmat1D_includeAnlVar, bmat1D_numIncludeAnlVar
37
38 integer :: bmat1D_numIncludeAnlVar
39 character(len=4), allocatable :: bmat1D_includeAnlVar(:)
40
41 type(struct_hco), pointer :: hco_yGrid
42 logical :: initialized = .false.
43 integer :: nkgdim
44 integer :: cvDim_mpilocal
45
46 real(8), allocatable :: bSqrtLand(:,:,:), bSqrtSea(:,:,:)
47 real(8), allocatable :: bSqrtEns(:,:,:)
48 real(4), allocatable :: latLand(:), lonLand(:), latSea(:), lonSea(:)
49 integer :: nLonLatPosLand, nLonLatPosSea
50 integer, external :: get_max_rss
51 integer, parameter :: numMasterBmat = 2
52 character(len=4), parameter :: masterBmatTypeList (numMasterBmat) = (/ 'HI', 'ENS' /)
53 character(len=8), parameter :: masterBmatLabelList(numMasterBmat) = (/'B_HI', 'B_ENS' /)
54 logical, parameter :: masterbmatIs3dList (numMasterBmat) = (/.true., .true. /)
55 integer :: numBmat
56 integer, parameter :: numBmatMax = 10
57 character(len=4) :: bmatTypeList (numBmatMax)
58 character(len=9) :: bmatLabelList (numBmatMax)
59 logical :: bmatIs3dList (numBmatMax)
60 logical :: bmatActive (numBmatMax)
61 type(struct_columnData), allocatable :: ensColumns(:)
62 type(struct_columnData) :: meanColumn
63 type(struct_ens) :: ensembles
64
65 ! Namelist variables
66 integer :: nEns ! ensemble size
67 real(8) :: vlocalize ! vertical localization length scale
68 character(len=4) :: includeAnlVar(vnl_numvarmax) ! list of variable names to include in B matrix
69 integer :: numIncludeAnlVar ! MUST NOT BE INCLUDED IN NAMELIST!
70 real(8) :: scaleFactorHI(vco_maxNumLevels) ! scaling factors for HI variances
71 real(8) :: scaleFactorHIHumidity(vco_maxNumLevels) ! scaling factors for HI humidity variances
72 real(8) :: scaleFactorEns(vco_maxNumLevels) ! scaling factors for Ens variances
73 real(8) :: scaleFactorEnsHumidity(vco_maxNumLevels) ! scaling factors for Ens humidity variances
74 logical :: dumpBmatrixTofile ! flag to control output of B matrices to Bmatrix.bin binary file
75 logical :: doAveraging ! flag to control output the average instead of the invidual B matrices
76 real(8) :: latMin ! minimum latitude of the Bmatrix latitude-longitude output box
77 real(8) :: latMax ! maximum latitude of the Bmatrix latitude-longitude output box
78 real(8) :: lonMin ! minimum longitude of the Bmatrix latitude-longitude output box
79 real(8) :: lonMax ! maximum longitude of the Bmatrix latitude-longitude output box
80 NAMELIST /NAMBMAT1D/ scaleFactorHI, scaleFactorHIHumidity, scaleFactorENs, scaleFactorEnsHumidity, nEns, &
81 vLocalize, includeAnlVar, numIncludeAnlVar, dumpBmatrixTofile, latMin, latMax, lonMin, lonMax, &
82 doAveraging
83
84contains
85
86 !--------------------------------------------------------------------------
87 ! bmat1D_bsetup
88 !--------------------------------------------------------------------------
89 subroutine bmat1D_bsetup(vco_in, hco_in, obsSpaceData)
90 !
91 !:Purpose: To initialize the 1Dvar analysis Background term.
92 !
93 implicit none
94
95 ! Arguments:
96 type(struct_vco), pointer, intent(in) :: vco_in
97 type(struct_hco), pointer, intent(in) :: hco_in
98 type (struct_obs), intent(inout) :: obsSpaceData
99
100 ! Locals:
101 integer :: cvdim
102 integer :: masterBmatIndex, bmatIndex
103 logical :: active
104 integer :: nulnam, ierr
105 integer, external :: fnom, fclos
106 integer :: varIndex
107
108 call utl_tmg_start(50, '--Bmatrix')
109
110 ! default values for namelist variables
111 scaleFactorHI(:) = 0.d0
112 scaleFactorHIHumidity(:) = 1.d0
113 scaleFactorEns(:) = 0.d0
114 scaleFactorEnsHumidity(:) = 1.d0
115 nEns = -1
116 vLocalize = -1.d0
117 includeAnlVar(:)= ''
118 numIncludeAnlVar = MPC_missingValue_INT
119 dumpBmatrixTofile = .false.
120 doAveraging = .false.
121 latMin = -1000.d0
122 latMax = 1000.d0
123 lonMin = -1000.d0
124 lonMax = 1000.d0
125
126 nulnam = 0
127 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
128 read(nulnam, nml=nambmat1D, iostat=ierr)
129 if ( ierr /= 0 ) call utl_abort( 'bmat1D_bsetup: Error reading namelist' )
130 if ( mmpi_myid == 0 ) write( *, nml = nambmat1D )
131 ierr = fclos( nulnam )
132 if (numIncludeAnlVar /= MPC_missingValue_INT) then
133 call utl_abort('bmat1D_bsetup: check NAMBMAT1D namelist section: numIncludeAnlVar should be removed')
134 end if
135 numIncludeAnlVar = 0
136 do varIndex = 1, vnl_numvarmax
137 if (trim(includeAnlVar(varIndex)) == '') exit
138 numIncludeAnlVar = numIncludeAnlVar + 1
139 end do
140 bmat1D_numIncludeAnlVar = numIncludeAnlVar
141 allocate( bmat1D_includeAnlVar(bmat1D_numIncludeAnlVar) )
142 bmat1D_includeAnlVar(1:bmat1D_numIncludeAnlVar) = includeAnlVar(1:numIncludeAnlVar)
143
144 !
145 !- 1. Setup the B matrices
146 !
147 do masterBmatIndex = 1, numMasterBmat
148
149 select case( trim(masterBmatTypeList(masterBmatIndex)) )
150 case ('HI')
151 !- 1.1 Time-Mean Homogeneous and Isotropic...
152 write(*,*) 'bmat1D_bsetup: Setting up the modular GLOBAL HI 1D covariances...'
153 call utl_tmg_start(51, '----B_HI_Setup')
154 call bmat1D_SetupBHi(vco_in, obsSpaceData, cvdim)
155 call utl_tmg_stop(51)
156 write(*,*) ' bmat1D_bsetup: cvdim= ', cvdim
157 case ('ENS')
158 !- 1.2 ensemble based
159 write(*,*) 'bmat1D_bsetup: Setting up the ensemble based 1D matrix.'
160 call utl_tmg_start(54, '----B_ENS_Setup')
161 call bmat1D_SetupBEns(vco_in, hco_in, obsSpaceData, cvdim)
162 call utl_tmg_stop(54)
163 write(*,*) ' bmat1D_bsetup: cvdim= ', cvdim
164 case default
165 call utl_abort('bmat1D_bSetup: requested bmatrix type does not exist ' // trim(masterBmatTypeList(masterBmatIndex)))
166 end select
167
168 !- 1.2 Append the info to the B matrix info arrays and setup the proper control sub-vectors
169 numBmat = numBmat + 1
170 bmatLabelList (numBmat) = trim(masterbmatLabelList(masterBmatIndex))
171 bmatTypeList (numBmat) = masterBmatTypeList(masterBmatIndex)
172 bmatIs3dList (numBmat) = masterbmatIs3dList(masterBmatIndex)
173 call cvm_setupSubVector(bmatLabelList(numBmat), bmatTypeList(numBmat), cvdim)
174 end do
175
176 !
177 !- 2. Print a summary and set the active B matrices array
178 !
179 write(*,*)
180 write(*,*) 'bmat1D_bsetup SUMMARY, number of B matrices found = ', numBmat
181 do bmatIndex = 1, numBmat
182 write(*,*) ' B matrix #', bmatIndex
183 active = cvm_subVectorExists(bmatLabelList(bmatIndex))
184 if (active) then
185 write(*,*) ' ACTIVE'
186 else
187 write(*,*) ' NOT USED'
188 end if
189 write(*,*) ' -> label = ', bmatLabelList (bmatIndex)
190 write(*,*) ' -> type = ', bmatTypeList (bmatIndex)
191 if (active) then
192 write(*,*) ' -> is 3D = ', bmatIs3dList (bmatIndex)
193 end if
194 bmatActive(bmatIndex) = active
195 end do
196
197 call utl_tmg_stop(50)
198
199 end subroutine bmat1D_bsetup
200
201 !--------------------------------------------------------------------------
202 ! bmat1D_setupBHi
203 !--------------------------------------------------------------------------
204 subroutine bmat1D_setupBHi(vco_in, obsSpaceData, cvDim_out)
205 !
206 !:Purpose: to setup bmat1D module
207 !
208 implicit none
209
210 ! Arguments:
211 type(struct_vco), pointer, intent(in) :: vco_in
212 type (struct_obs) , intent(in) :: obsSpaceData
213 integer , intent(out) :: cvDim_out
214
215 ! Locals:
216 integer :: levelIndex, ierr
217 integer, external :: fnom, fclos
218 integer :: Vcode_anl
219 logical :: fileExists
220 integer :: nulbgst=0
221 type(struct_vco), pointer :: vco_file => null()
222 type(struct_vco), target :: vco_1Dvar
223 type(struct_vco), pointer :: vco_anl
224 character(len=18) :: oneDBmatLand = './Bmatrix_land.bin'
225 character(len=17) :: oneDBmatSea = './Bmatrix_sea.bin'
226 character(len=4), allocatable :: includeAnlVarHi(:)
227 integer :: extractDate, locationIndex, varIndex, numIncludeAnlVarHi
228 integer :: shiftLevel, varLevIndex1, varLevIndex2
229 integer :: varLevIndexBmat
230 logical, save :: firstCall=.true.
231 real(8), allocatable :: bMatrix(:,:), multFactor(:)
232
233 if (.not. (gsv_varExist(varName='TT') .and. &
234 gsv_varExist(varName='UU') .and. &
235 gsv_varExist(varName='VV') .and. &
236 (gsv_varExist(varName='HU').or.gsv_varExist(varName='LQ')) .and. &
237 gsv_varExist(varName='P0')) ) then
238 call utl_abort('bmat1D_setupBHi: Some or all weather fields are missing. If it is desired to deactivate&
239 & the weather assimilation, then all entries of the array SCALEFACTORHI in the namelist NAMVAR1D&
240 & should be set to zero.')
241 end if
242
243 if (.not. gsv_varExist(varName='TG')) then
244 write(*,*) 'bmat1D_setupBHi: WARNING: The TG field is missing. This must be present when assimilating'
245 write(*,*) 'radiance observations.'
246 end if
247
248 if (firstCall) then
249 call var1D_setup(obsSpaceData)
250 firstCall = .false.
251 end if
252
253 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBHi: Starting'
254 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
255
256 do levelIndex = 1, vco_maxNumLevels
257 if( scaleFactorHI(levelIndex) > 0.0d0 ) then
258 scaleFactorHI(levelIndex) = sqrt( scaleFactorHI(levelIndex))
259 else
260 scaleFactorHI(levelIndex) = 0.0d0
261 end if
262 end do
263
264 do levelIndex = 1, vco_maxNumLevels
265 if(scaleFactorHIHumidity(levelIndex) > 0.0d0) then
266 scaleFactorHIHumidity(levelIndex) = sqrt(scaleFactorHIHumidity(levelIndex))
267 else
268 scaleFactorHIHumidity(levelIndex) = 0.0d0
269 end if
270 end do
271
272 if ( sum(scaleFactorHI(1:vco_maxNumLevels)) == 0.0d0 ) then
273 if ( mmpi_myid == 0 ) write(*,*) 'bmat1D_setupBHi: scaleFactorHI=0, skipping rest of setup'
274 cvDim_out = 0
275 return
276 end if
277
278 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBHi: Read 1DVar background statistics'
279 inquire(file=trim(oneDBmatLand), exist=fileExists)
280 if ( fileExists ) then
281 ierr = fnom(nulbgst, trim(oneDBmatLand), 'FTN+SEQ+UNF+OLD+R/O', 0)
282 else
283 call utl_abort('bmat1D_setupBHi: No 1DVar BACKGROUND STAT FILE ' // trim(oneDBmatLand))
284 end if
285
286 read(nulbgst) extractDate, vco_1Dvar%nLev_T, vco_1Dvar%nLev_M, vco_1Dvar%Vcode, &
287 vco_1Dvar%ip1_sfc, vco_1Dvar%ip1_T_2m, vco_1Dvar%ip1_M_10m, numIncludeAnlVarHi, nkgdim, nLonLatPosLand
288 allocate( vco_1Dvar%ip1_T(vco_1Dvar%nLev_T), vco_1Dvar%ip1_M(vco_1Dvar%nLev_M) )
289 if (numIncludeAnlVarHi /= bmat1D_numIncludeAnlVar) then
290 write(*,*) 'numIncludeAnlVarHi, bmat1D_numIncludeAnlVar= ', numIncludeAnlVarHi, bmat1D_numIncludeAnlVar
291 call utl_abort('bmat1D_setupBHi: incompatible number of 1DVar analyzed variables in ' // trim(oneDBmatLand))
292 end if
293
294 allocate (multFactor(nkgdim))
295 if(vco_1Dvar%Vcode == 5002) then
296 shiftLevel = 1
297 else
298 shiftLevel = 0
299 end if
300 varLevIndexBmat = 0
301 do varIndex = 1, bmat1D_numIncludeAnlVar
302 select case(bmat1D_includeAnlVar(varIndex))
303 case('TT')
304 do levelIndex = 1, vco_1Dvar%nLev_T
305 varLevIndexBmat = varLevIndexBmat + 1
306 multFactor(varLevIndexBmat) = scaleFactorHI(levelIndex)
307 end do
308 case('HU')
309 do levelIndex = 1, vco_1Dvar%nLev_T
310 varLevIndexBmat = varLevIndexBmat + 1
311 multFactor(varLevIndexBmat)= scaleFactorHIHumidity(levelIndex) * scaleFactorHI(levelIndex)
312 end do
313 case('UU','VV')
314 do levelIndex = 1, vco_1Dvar%nLev_M
315 varLevIndexBmat = varLevIndexBmat + 1
316 multFactor(varLevIndexBmat) = scaleFactorHI(levelIndex+shiftLevel)
317 end do
318 case('P0','TG')
319 varLevIndexBmat = varLevIndexBmat + 1
320 multFactor(varLevIndexBmat) = scaleFactorHI(max(vco_1Dvar%nLev_T,vco_1Dvar%nLev_M))
321 case default
322 call utl_abort('bmat1D_setupBHi: unsupported variable ' // bmat1D_includeAnlVar(varIndex))
323 end select
324 end do
325
326 allocate( includeAnlVarHi(bmat1D_numIncludeAnlVar) )
327 allocate( bMatrix(nkgdim,nkgdim) )
328 allocate( latLand(nLonLatPosLand), lonLand(nLonLatPosLand))
329 allocate( bSqrtLand(nLonLatPosLand, nkgdim, nkgdim) )
330 read(nulbgst) vco_1Dvar%ip1_T(:), vco_1Dvar%ip1_M(:), includeAnlVarHi(:)
331 if (any(includeAnlVarHi /= bmat1D_includeAnlVar)) then
332 do varIndex = 1, bmat1D_numIncludeAnlVar
333 write(*,*) varIndex, includeAnlVarHi(varIndex), bmat1D_includeAnlVar(varIndex)
334 end do
335 call utl_abort('bmat1D_setupBHi: incompatible 1DVar analyzed variable list in ' // trim(oneDBmatLand))
336 end if
337 deallocate(includeAnlVarHi)
338
339 do locationIndex = 1, nLonLatPosLand
340 read(nulbgst) latLand(locationIndex), lonLand(locationIndex), bMatrix(:,:)
341 !application of the scaling factor
342 do varLevIndex2 = 1, nkgdim
343 do varLevIndex1 = 1, nkgdim
344 bSqrtLand(locationIndex, varLevIndex1, varLevIndex2) = &
345 bMatrix(varLevIndex1, varLevIndex2) * multFactor(varLevIndex1) * multFactor(varLevIndex2)
346 end do
347 end do
348 call utl_matsqrt(bSqrtLand(locationIndex, :, :), nkgdim, 1.d0, printInformation_opt=.false. )
349 end do
350 ierr = fclos(nulbgst)
351
352 inquire(file=trim(oneDBmatSea), exist=fileExists)
353 if ( fileExists ) then
354 ierr = fnom(nulbgst, trim(oneDBmatSea), 'FTN+SEQ+UNF+OLD+R/O', 0)
355 else
356 call utl_abort('bmat1D_setupBHi: No 1DVar BACKGROUND STAT FILE ' // trim(oneDBmatSea))
357 end if
358 read(nulbgst) extractDate, vco_1Dvar%nLev_T, vco_1Dvar%nLev_M, vco_1Dvar%Vcode, &
359 vco_1Dvar%ip1_sfc, vco_1Dvar%ip1_T_2m, vco_1Dvar%ip1_M_10m, numIncludeAnlVarHi, nkgdim, nLonLatPosSea
360 if (numIncludeAnlVarHi /= bmat1D_numIncludeAnlVar) then
361 write(*,*) 'numIncludeAnlVarHi, bmat1D_numIncludeAnlVar= ', numIncludeAnlVarHi, bmat1D_numIncludeAnlVar
362 call utl_abort('bmat1D_setupBHi: incompatible number of 1DVar analyzed variables in ' // trim(oneDBmatSea))
363 end if
364 allocate( bSqrtSea(nLonLatPosSea, nkgdim, nkgdim) )
365 allocate( latSea(nLonLatPosSea), lonSea(nLonLatPosSea))
366 allocate( includeAnlVarHi(bmat1D_numIncludeAnlVar) )
367 read(nulbgst) vco_1Dvar%ip1_T(:), vco_1Dvar%ip1_M(:), includeAnlVarHi(:)
368 if (any(includeAnlVarHi /= bmat1D_includeAnlVar)) then
369 do varIndex = 1, bmat1D_numIncludeAnlVar
370 write(*,*) varIndex, includeAnlVarHi(varIndex), bmat1D_includeAnlVar(varIndex)
371 end do
372 call utl_abort('bmat1D_setupBHi: incompatible 1DVar analyzed variable list in ' // trim(oneDBmatSea))
373 end if
374 deallocate(includeAnlVarHi)
375 do locationIndex = 1, nLonLatPosSea
376 read(nulbgst) latSea(locationIndex), lonSea(locationIndex), bMatrix(:,:)
377 !application of the scaling factor
378 do varLevIndex2 = 1, nkgdim
379 do varLevIndex1 = 1, nkgdim
380 bSqrtSea(locationIndex, varLevIndex1, varLevIndex2) = &
381 bMatrix(varLevIndex1, varLevIndex2) * multFactor(varLevIndex1) * multFactor(varLevIndex2)
382 end do
383 end do
384 call utl_matsqrt(bSqrtSea(locationIndex,:,:), nkgdim, 1.d0, printInformation_opt=.false. )
385 end do
386 ierr = fclos(nulbgst)
387
388 deallocate(bMatrix, multFactor)
389
390 vco_1Dvar%initialized = .true.
391 vco_1Dvar%vGridPresent = .false.
392 vco_file => vco_1Dvar
393 vco_anl => vco_in
394 if (.not. vco_equal(vco_anl,vco_file)) then
395 call utl_abort('bmat1D_setupBHi: vco from analysisgrid and cov file do not match')
396 end if
397 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBHi: nLev_M, nLev_T=', vco_1Dvar%nLev_M, vco_1Dvar%nLev_T
398 Vcode_anl = vco_anl%vCode
399 if(Vcode_anl /= 5002 .and. Vcode_anl /= 5005) then
400 write(*,*) 'Vcode_anl = ',Vcode_anl
401 call utl_abort('bmat1D_setupBHi: unknown vertical coordinate type!')
402 end if
403
404 cvDim_out = nkgdim * var1D_validHeaderCount
405 cvDim_mpilocal = cvDim_out
406 initialized = .true.
407
408 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBHi: Exiting'
409 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
410
411 end subroutine bmat1D_setupBHi
412
413 !--------------------------------------------------------------------------
414 ! bmat1D_setupBEns
415 !--------------------------------------------------------------------------
416 subroutine bmat1D_setupBEns(vco_in, hco_in, obsSpaceData, cvDim_out)
417 !
418 !:Purpose: to setup bmat1D module
419 !
420 implicit none
421
422 ! Arguments
423 type(struct_vco), pointer, intent(in) :: vco_in
424 type(struct_hco), pointer, intent(in) :: hco_in
425 type (struct_obs) , intent(inout) :: obsSpaceData
426 integer , intent(out) :: cvDim_out
427
428 ! Locals:
429 character(len=256) :: ensPathName = 'ensemble'
430 character(len=256) :: ensFileName
431 type(struct_vco), pointer :: vco_file => null()
432 type(struct_vco), pointer :: vco_ens => null()
433 type(struct_hco), pointer :: hco_ens => null()
434 type(struct_gsv) :: stateVector, stateVectorMean
435 integer, allocatable :: dateStampList(:)
436 character(len=12) :: hInterpolationDegree='LINEAR' ! select degree of horizontal interpolation (if needed)
437 integer :: memberIndex, columnIndex, headerIndex, varIndex, levIndex
438 integer :: varLevIndex1, varLevIndex2
439 integer :: varLevIndexBmat, varLevIndexCol
440 integer :: numStep, levIndexColumn
441 real(8), allocatable :: scaleFactor_M(:), scaleFactor_T(:)
442 real(8) :: scaleFactor_SF, zr
443 logical :: useAnlLevelsOnly, EnsTopMatchesAnlTop
444 real(8), pointer :: pressureProfileFile_M(:), pressureProfileInc_M(:)
445 real(8) :: pSurfRef
446 integer :: nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd
447 integer :: ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd
448 integer :: nLevEns_M, nLevEns_T
449 integer :: nLevInc_M, nLevInc_T
450 real(8) :: logP1, logP2
451 real(8), pointer :: currentProfile(:), meanProfile(:)
452 real(8), allocatable :: lineVector(:,:), meanPressureProfile(:), multFactor(:)
453 integer, allocatable :: varLevColFromVarLevBmat(:)
454 character(len=4), allocatable :: varNameFromVarLevIndexBmat(:)
455 character(len=2) :: varLevel
456
457 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: Starting'
458 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
459
460 if (nEns <= 0) then
461 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: no Ensemble members, skipping rest of setup'
462 cvdim_out = 0
463 return
464 end if
465
466 !- 1.1 Number of time step bins
467 numStep = tim_nstepobsinc
468 if (numStep /= 1 .and. numStep /= 3.and. numStep /= 5 .and. numStep /= 7) then
469 call utl_abort('bmat1D_setupBEns: Invalid value for numStep (choose 1 or 3 or 5 or 7)!')
470 end if
471 allocate(dateStampList(numStep))
472 call tim_getstamplist(dateStampList,numStep,tim_getDatestamp())
473
474 hco_ens => hco_in
475
476 !- 1.2 Horizontal grid
477 ni = hco_ens%ni
478 nj = hco_ens%nj
479 if (hco_ens%global) then
480 if (mmpi_myid == 0) write(*,*)
481 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: GLOBAL mode activated'
482 else
483 if (mmpi_myid == 0) write(*,*)
484 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: LAM mode activated'
485 end if
486
487 !- 1.3 Vertical levels
488 call fln_ensfileName(ensFileName, ensPathName, memberIndex_opt=1)
489 call vco_SetupFromFile(vco_file, ensFileName)
490
491 !- Do we need to read all the vertical levels from the ensemble?
492 useAnlLevelsOnly = vco_subsetOrNot(vco_in, vco_file)
493 if (useAnlLevelsOnly) then
494 write(*,*)
495 write(*,*) 'bmat1D_setupBEns: only the analysis levels will be read in the ensemble '
496 vco_ens => vco_in ! the ensemble target grid is the analysis grid
497 call vco_deallocate(vco_file)
498 vco_file => vco_in ! only the analysis levels will be read in the ensemble
499 EnsTopMatchesAnlTop = .true.
500 else
501 write(*,*)
502 write(*,*) 'bmat1D_setupBEns: all the vertical levels will be read in the ensemble '
503 if ( vco_in%nLev_M > 0 .and. vco_in%vgridPresent ) then
504 pSurfRef = 101000.D0
505 call czp_fetch1DLevels(vco_in, pSurfRef, profM_opt=pressureProfileInc_M)
506 call czp_fetch1DLevels(vco_in, pSurfRef, profM_opt=pressureProfileFile_M)
507
508 EnsTopMatchesAnlTop = abs( log(pressureProfileFile_M(1)) - log(pressureProfileInc_M(1)) ) < 0.1d0
509 write(*,*) 'bmat1D_setupBEns: EnsTopMatchesAnlTop, presEns, presInc = ', &
510 EnsTopMatchesAnlTop, pressureProfileFile_M(1), pressureProfileInc_M(1)
511 deallocate(pressureProfileFile_M)
512 deallocate(pressureProfileInc_M)
513 else
514 ! not sure what this mean when no MM levels
515 write(*,*) 'bmat1D_setupBEns: nLev_M = ', vco_in%nLev_M
516 write(*,*) 'bmat1D_setupBEns: vgridPresent = ', vco_in%vgridPresent
517 EnsTopMatchesAnlTop = .true.
518 end if
519
520 if ( EnsTopMatchesAnlTop ) then
521 if ( mmpi_myid == 0 ) write(*,*) 'bmat1D_setupBEns: top level of ensemble member and analysis grid match'
522 vco_ens => vco_in ! IMPORTANT: top levels DO match, therefore safe
523 ! to force members to be on analysis vertical levels
524 else
525 if ( mmpi_myid == 0 ) write(*,*) 'bmat1D_setupBEns: top level of ensemble member and analysis grid are different, therefore'
526 if ( mmpi_myid == 0 ) write(*,*) ' assume member is already be on correct levels - NO CHECKING IS DONE'
527 vco_ens => vco_file ! IMPORTANT: top levels do not match, therefore must
528 ! assume file is already on correct vertical levels
529 end if
530 end if
531
532 if (vco_in%Vcode /= vco_ens%Vcode) then
533 write(*,*) 'bmat1D_setupBEns: vco_in%Vcode = ', vco_in%Vcode, ', vco_ens%Vcode = ', vco_ens%Vcode
534 call utl_abort('bmat1D_setupBEns: vertical levels of ensemble not compatible with analysis grid')
535 end if
536 nLevEns_M = vco_ens%nLev_M
537 nLevEns_T = vco_ens%nLev_T
538 nLevInc_M = vco_in%nLev_M
539 nLevInc_T = vco_in%nLev_T
540
541 if (vco_in%Vcode == 5002) then
542 if ( (nLevEns_T /= (nLevEns_M+1)) .and. (nLevEns_T /= 1 .or. nLevEns_M /= 1) ) then
543 write(*,*) 'bmat1D_setubBEns: nLevEns_T, nLevEns_M = ',nLevEns_T,nLevEns_M
544 call utl_abort('bmat1D_setubBEns: Vcode=5002, nLevEns_T must equal nLevEns_M+1!')
545 end if
546 else if (vco_in%Vcode == 5005) then
547 if ( nLevEns_T /= nLevEns_M .and. &
548 nLevEns_T /= 0 .and. &
549 nLevEns_M /= 0 ) then
550 write(*,*) 'bmat1D_setubBEns: nLevEns_T, nLevEns_M = ',nLevEns_T,nLevEns_M
551 call utl_abort('bmat1D_setubBEns: Vcode=5005, nLevEns_T must equal nLevEns_M!')
552 end if
553 else if (vco_in%Vcode == 0) then
554 if ( nLevEns_T /= 0 .and. nLevEns_M /= 0 ) then
555 write(*,*) 'bmat1D_setubBEns: nLevEns_T, nLevEns_M = ',nLevEns_T, nLevEns_M
556 call utl_abort('bmat1D_setubBEns: surface-only case (Vcode=0), nLevEns_T and nLevEns_M must equal 0!')
557 end if
558 else
559 write(*,*) 'bmat1D_setubBEns: vco_in%Vcode = ',vco_in%Vcode
560 call utl_abort('bmat1D_setubBEns: unknown vertical coordinate type!')
561 end if
562
563 if (nLevEns_M > nLevInc_M) then
564 call utl_abort('bmat1D_setubBEns: ensemble has more levels than increment - not allowed!')
565 end if
566
567 if (nLevEns_M < nLevInc_M) then
568 if (mmpi_myid == 0) write(*,*) 'bmat1D_setubBEns: ensemble has less levels than increment'
569 if (mmpi_myid == 0) write(*,*) ' some levels near top will have zero increment'
570 end if
571
572 !- 1.4 Bmatrix Weight
573 if (vco_in%Vcode == 5002 .or. vco_in%Vcode == 5005) then
574 allocate(scaleFactor_M(nLevEns_M))
575 allocate(scaleFactor_T(nLevEns_T))
576 do levIndex = 1, nLevEns_T
577 if (scaleFactorEns(levIndex) > 0.0d0) then
578 scaleFactorEns(levIndex) = sqrt(scaleFactorEns(levIndex))
579 else
580 scaleFactorEns(levIndex) = 0.0d0
581 end if
582 end do
583 scaleFactor_T(1:nLevEns_T) = scaleFactorEns(1:nLevEns_T)
584 if (vco_in%Vcode == 5002) then
585 scaleFactor_M(1:nLevEns_M) = scaleFactorEns(2:(nLevEns_M+1))
586 else
587 scaleFactor_M(1:nLevEns_M) = scaleFactorEns(1:nLevEns_M)
588 end if
589
590 do levIndex = 1, nLevEns_T
591 if (scaleFactorEnsHumidity(levIndex) > 0.0d0) then
592 scaleFactorEnsHumidity(levIndex) = sqrt(scaleFactorEnsHumidity(levIndex))
593 else
594 scaleFactorEnsHumidity(levIndex) = 0.0d0
595 end if
596 end do
597
598 scaleFactor_SF = scaleFactor_T(nLevEns_T)
599
600 else ! vco_in%Vcode == 0
601 if (scaleFactorEns(1) > 0.0d0) then
602 scaleFactor_SF = sqrt(scaleFactorEns(1))
603 else
604 call utl_abort('bmat1D_setubBEns: with vCode == 0, the scale factor should never be equal to 0')
605 end if
606 end if
607
608 !- 1.5 Domain Partionning
609 call mmpi_setup_latbands(nj, latPerPE, latPerPEmax, myLatBeg, myLatEnd)
610 call mmpi_setup_lonbands(ni, lonPerPE, lonPerPEmax, myLonBeg, myLonEnd)
611
612 !- 1.6 Localization
613 if ( vLocalize <= 0.0d0 .and. (nLevInc_M > 1 .or. nLevInc_T > 1) ) then
614 call utl_abort('bmat1D_setubBEns: Invalid VERTICAL localization length scale')
615 end if
616
617 call ens_allocate(ensembles, &
618 nEns, numStep, &
619 hco_ens, &
620 vco_ens, dateStampList, &
621 hco_core_opt = hco_in, &
622 varNames_opt = bmat1D_includeAnlVar(1:bmat1D_numIncludeAnlVar), &
623 hInterpolateDegree_opt = hInterpolationDegree)
624 write(*,*) 'Read ensemble members'
625 call ens_readEnsemble(ensembles, ensPathName, biPeriodic=.false., &
626 varNames_opt = bmat1D_includeAnlVar(1:bmat1D_numIncludeAnlVar))
627
628 allocate(ensColumns(nEns))
629 call gsv_allocate(stateVector, numstep, hco_ens, vco_ens, &
630 dateStamp_opt=tim_getDateStamp(), &
631 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
632 dataKind_opt=4, allocHeightSfc_opt=.true.)
633
634 call gsv_allocate(stateVectorMean, numstep, hco_ens, vco_ens, &
635 dateStamp_opt=tim_getDateStamp(), &
636 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
637 dataKind_opt=4, allocHeightSfc_opt=.true.)
638 call gsv_zero(stateVectorMean)
639 do memberIndex = 1, nEns
640 write(*,*) 'Copy member ', memberIndex
641 call gsv_zero(stateVector)
642 call ens_copyMember(ensembles, stateVector, memberIndex)
643 write(*,*) 'interpolate member ', memberIndex
644 call col_setVco(ensColumns(memberIndex), vco_ens)
645 call col_allocate(ensColumns(memberIndex), obs_numheader(obsSpaceData), &
646 mpiLocal_opt=.true., setToZero_opt=.true.)
647 call s2c_nl(stateVector, obsSpaceData, ensColumns(memberIndex), hco_in, &
648 timeInterpType='NEAREST', dealloc_opt=.false. )
649 call gsv_add(statevector, statevectorMean, scaleFactor_opt=(1.d0/nEns))
650 end do
651
652 call col_setVco(meanColumn, vco_ens)
653 call col_allocate(meanColumn, obs_numheader(obsSpaceData), &
654 mpiLocal_opt=.true., setToZero_opt=.true.)
655 call s2c_nl(stateVectorMean, obsSpaceData, meanColumn, hco_in, &
656 timeInterpType='NEAREST' )
657
658 call gsv_deallocate(stateVector)
659 call gsv_deallocate(stateVectorMean)
660
661 nkgdim = 0
662 do varIndex = 1, bmat1D_numIncludeAnlVar
663 currentProfile => col_getColumn(meanColumn, var1D_validHeaderIndex(1), varName_opt=bmat1D_includeAnlVar(varIndex))
664 nkgdim = nkgdim + size(currentProfile)
665 end do
666 write(*,*) 'bmat1D_setupBEns: nkgdim', nkgdim
667 cvDim_out = nkgdim * var1D_validHeaderCount
668
669 currentProfile => col_getColumn(meanColumn, var1D_validHeaderIndex(1))
670 allocate (varLevColFromVarLevBmat(nkgdim))
671 allocate (varNameFromVarLevIndexBmat(nkgdim))
672 allocate (multFactor(nkgdim))
673 varLevIndexBmat = 0
674 do varIndex = 1, bmat1D_numIncludeAnlVar
675 levIndex = 0
676 do varLevIndexCol = 1, size(currentProfile)
677 if ( trim( col_getVarNameFromK(meanColumn,varLevIndexCol) ) == trim( bmat1D_includeAnlVar(varIndex) ) ) then
678 varLevIndexBmat = varLevIndexBmat + 1
679 levIndex = levIndex + 1
680 varLevColFromVarLevBmat(varLevIndexBmat) = varLevIndexCol
681 varNameFromVarLevIndexBmat(varLevIndexBmat) = trim( bmat1D_includeAnlVar(varIndex) )
682 varLevel = vnl_varLevelFromVarname(varNameFromVarLevIndexBmat(varLevIndexBmat))
683 if ( varLevel == 'MM' ) then ! Momentum
684 multFactor(varLevIndexBmat) = scaleFactor_M(levIndex)
685 else if ( varLevel == 'TH' ) then ! Thermo
686 multFactor(varLevIndexBmat) = scaleFactor_T(levIndex)
687 else ! SF
688 multFactor(varLevIndexBmat) = scaleFactor_SF
689 end if
690 if (varNameFromVarLevIndexBmat(varLevIndexBmat) == 'HU') then
691 multFactor(varLevIndexBmat) = multFactor(varLevIndexBmat) * scaleFactorEnsHumidity(levIndex)
692 end if
693 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: bmat1D_includeAnlVar ', bmat1D_includeAnlVar(varIndex), varLevIndexBmat, levIndex
694 end if
695 end do
696 end do
697
698 call ens_deallocate(ensembles)
699 allocate(bSqrtEns(var1D_validHeaderCount,nkgdim,nkgdim))
700 bSqrtEns(:,:,:) = 0.d0
701 allocate(lineVector(1,nkgdim))
702
703 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,memberIndex,meanProfile,currentProfile,lineVector)
704 do columnIndex = 1, var1D_validHeaderCount
705 headerIndex = var1D_validHeaderIndex(columnIndex)
706 meanProfile => col_getColumn(meanColumn, headerIndex)
707 do memberIndex = 1, nEns
708 currentProfile => col_getColumn(ensColumns(memberIndex), headerIndex)
709 lineVector(1,:) = currentProfile(varLevColFromVarLevBmat(:)) - meanProfile(varLevColFromVarLevBmat(:))
710 lineVector(1,:) = lineVector(1,:) * multFactor(:)
711 bSqrtEns(columnIndex,:,:) = bSqrtEns(columnIndex,:,:) + &
712 matmul(transpose(lineVector),lineVector)
713 end do
714 end do
715 !$OMP END PARALLEL DO
716
717 do varLevIndexBmat =1, nkgdim
718 write(*,*) 'bmat1D_setupBEns: variance = ', varLevIndexBmat, &
719 sum(bSqrtEns(:,varLevIndexBmat,varLevIndexBmat))/real((nEns-1)*var1D_validHeaderCount)
720 end do
721
722 deallocate(lineVector)
723 deallocate(multFactor)
724 allocate(meanPressureProfile(nkgdim))
725 call lfn_Setup('FifthOrder')
726
727 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,meanPressureProfile,varLevIndexBmat,varLevIndexCol,varLevIndex1,varLevIndex2,varLevel,levIndexColumn,logP1,logP2,zr)
728 do columnIndex = 1, var1D_validHeaderCount
729 headerIndex = var1D_validHeaderIndex(columnIndex)
730 bSqrtEns(columnIndex,:,:) = bSqrtEns(columnIndex,:,:) / (nEns - 1)
731 if (vLocalize > 0.0d0) then
732 do varLevIndexBmat = 1, nkgdim
733 varLevIndexCol = varLevColFromVarLevBmat(varLevIndexBmat)
734 varLevel = vnl_varLevelFromVarname(varNameFromVarLevIndexBmat(varLevIndexBmat))
735 if (varLevel=='SF') then
736 meanPressureProfile(varLevIndexBmat) = col_getElem(meanColumn, 1, headerIndex, 'P0')
737 else
738 levIndexColumn = col_getLevIndexFromVarLevIndex(meanColumn, varLevIndexCol)
739 meanPressureProfile(varLevIndexBmat) = col_getPressure(meanColumn, levIndexColumn, headerIndex, varLevel)
740 end if
741 end do
742 do varLevIndex1 = 1, nkgdim
743 logP1 = log(meanPressureProfile(varLevIndex1))
744 do varLevIndex2 = 1, nkgdim
745 !- do Schurr product with 5'th order function
746 logP2 = log(meanPressureProfile(varLevIndex2))
747 zr = abs(logP2 - logP1)
748 bSqrtEns(columnIndex, varLevIndex2, varLevIndex1) = &
749 bSqrtEns(columnIndex, varLevIndex2, varLevIndex1) * lfn_response(zr, vLocalize)
750 end do
751 end do
752 end if
753 end do
754 !$OMP END PARALLEL DO
755
756 do varLevIndexBmat =1, nkgdim
757 write(*,*) 'bmat1D_setupBEns: variance (after localization) = ', varLevIndexBmat, &
758 sum(bSqrtEns(:,varLevIndexBmat,varLevIndexBmat))/real(var1D_validHeaderCount)
759 end do
760
761 if (dumpBmatrixTofile) then
762 call dumpBmatrices(vco_in, obsSpaceData, dateStampList, bSqrtEns)
763 end if
764
765 write(*,*) 'bmat1D_setupBEns: computing matrix sqrt'
766 do columnIndex = 1, var1D_validHeaderCount
767 call utl_matsqrt(bSqrtEns(columnIndex, :, :), nkgdim, 1.d0, printInformation_opt=.false. )
768 end do
769
770 deallocate(varLevColFromVarLevBmat)
771 deallocate(varNameFromVarLevIndexBmat)
772 deallocate(meanPressureProfile)
773 call col_deallocate(meanColumn)
774 do memberIndex = 1, nEns
775 call col_deallocate(ensColumns(memberIndex))
776 end do
777 cvDim_mpilocal = cvDim_out
778 initialized = .true.
779
780 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
781 if (mmpi_myid == 0) write(*,*) 'bmat1D_setupBEns: Exiting'
782
783 end subroutine bmat1D_setupBEns
784
785 !--------------------------------------------------------------------------
786 ! dumpBmatrices
787 !--------------------------------------------------------------------------
788 subroutine dumpBmatrices(vco_in, obsSpaceData, dateStampList, bEns)
789 !
790 !:Purpose: output B matrices or their average to binary file Bmatrix.bin
791 !
792 implicit none
793
794 ! Arguments:
795 type(struct_vco), pointer, intent(in) :: vco_in ! Analysis vertical grid descriptor
796 type (struct_obs) , intent(inout) :: obsSpaceData ! ObsSpaceData structure
797 integer , intent(in) :: dateStamplist(:) ! List of datestamps
798 real(8) , intent(in) :: bEns(:,:,:) ! B matrices (first dimension is location index)
799
800 ! Locals:
801 integer :: nulmat, ierr
802 integer :: yyyymmdd, hhmm, countDumped, countDumpedMax, countDumpedMpiGlobal
803 integer :: globalDumpedIndex, countDumpedOut
804 integer :: columnIndex, headerIndex
805 integer :: taskIndex, dumpedIndex
806 integer :: tag, status
807 integer :: numstep, nkgdim
808 integer, external :: fnom, fclos, newdate
809 integer :: obsOffset(0:mmpi_nprocs-1)
810 integer :: countDumpedAllTasks(mmpi_nprocs)
811 integer, allocatable :: listColumnDumped(:)
812 real(8) :: latitude, longitude
813 real(8), allocatable :: tempoBmatrix(:,:)
814 real(8), allocatable :: outLats(:),outLons(:), outBmatrix(:,:,:)
815
816 if (mmpi_myid == 0) write(*,*) 'dumpBmatrices: Starting'
817 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
818
819 countDumped = 0
820 do columnIndex = 1, var1D_validHeaderCount
821 headerIndex = var1D_validHeaderIndex(columnIndex)
822 latitude = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex) * MPC_DEGREES_PER_RADIAN_R8
823 longitude = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex) * MPC_DEGREES_PER_RADIAN_R8
824 if (latitude >= latMin .and. &
825 latitude <= latMax .and. &
826 longitude >= lonMin .and. &
827 longitude <= lonMax) then
828 countDumped = countDumped + 1
829 end if
830 end do
831 call rpn_comm_barrier('GRID', ierr)
832
833 call rpn_comm_allgather(countDumped, 1, 'MPI_INTEGER', countDumpedAllTasks, 1,'MPI_INTEGER', 'GRID', ierr)
834 countDumpedMpiGlobal = sum( countDumpedAllTasks(:) )
835 countDumpedMax = maxval( countDumpedAllTasks(:) )
836 obsOffset(0) = 0
837 do taskIndex = 1, mmpi_nprocs - 1
838 obsOffset(taskIndex) = obsOffset(taskIndex - 1) + countDumpedAllTasks(taskIndex)
839 end do
840 write(*,*) 'dumpBmatrices: obsOffset: ', obsOffset(:)
841
842 if (countDumpedMax > 0) then
843 allocate(listColumnDumped(countDumpedMax))
844 listColumnDumped(:) = -1
845 countDumped = 0
846 do columnIndex = 1, var1D_validHeaderCount
847 headerIndex = var1D_validHeaderIndex(columnIndex)
848 latitude = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex) * MPC_DEGREES_PER_RADIAN_R8
849 longitude = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex) * MPC_DEGREES_PER_RADIAN_R8
850 if (latitude >= latMin .and. &
851 latitude <= latMax .and. &
852 longitude >= lonMin .and. &
853 longitude <= lonMax) then
854 countDumped = countDumped + 1
855 listColumnDumped(countDumped) = columnIndex
856 end if
857 end do
858 nkgdim = size(bEns, dim=2)
859 if (mmpi_myId == 0) then
860 nulmat = 0
861 ierr = fnom(nulmat, './Bmatrix.bin', 'FTN+SEQ+UNF', 0)
862 numstep = size(dateStampList)
863 ierr = newdate(dateStampList(1 + numstep / 2), yyyymmdd, hhmm, -3)
864 countDumpedOut = countDumped
865 if (doAveraging) countDumpedOut = 1
866 write(nulmat) yyyymmdd * 100 + nint(hhmm/100.), vco_in%nlev_T, vco_in%nlev_M, vco_in%Vcode, &
867 vco_in%ip1_sfc, vco_in%ip1_T_2m, vco_in%ip1_M_10m, bmat1D_numIncludeAnlVar, nkgdim, countDumpedOut
868 write(nulmat) vco_in%ip1_T(:), vco_in%ip1_M(:), bmat1D_includeAnlVar(1:bmat1D_numIncludeAnlVar)
869 allocate(outLats(countDumpedMpiGlobal), outLons(countDumpedMpiGlobal), OutBmatrix(countDumpedMpiGlobal, nkgdim, nkgdim))
870 allocate(tempoBmatrix(nkgdim, nkgdim))
871 end if
872 end if
873 do dumpedIndex = 1, countDumpedMax
874 columnIndex = listColumnDumped(dumpedIndex)
875 if (columnIndex > 0) then
876 headerIndex = var1D_validHeaderIndex(columnIndex)
877 latitude = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex)
878 longitude = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex)
879 else
880 latitude = MPC_missingValue_R8
881 longitude = MPC_missingValue_R8
882 end if
883 if (mmpi_myId == 0) then
884 if (latitude /= MPC_missingValue_R8 .and. longitude /= MPC_missingValue_R8) then
885 outBmatrix(dumpedIndex, :, :) = bEns(columnIndex, :, :)
886 outLats(dumpedIndex) = latitude
887 outLons(dumpedIndex) = longitude
888 end if
889 do taskIndex = 1, mmpi_nprocs - 1
890 tag = 3 * taskIndex
891 call rpn_comm_recv(latitude, 1, 'mpi_real8', taskIndex, tag, 'GRID', status, ierr)
892 call rpn_comm_recv(longitude, 1, 'mpi_real8', taskIndex, tag+1, 'GRID', status, ierr)
893 if (latitude /= MPC_missingValue_R8 .and. longitude /= MPC_missingValue_R8) then
894 call rpn_comm_recv(tempoBmatrix(:,:), nkgdim*nkgdim, 'mpi_real8', taskIndex, tag + 2, 'GRID', status, ierr)
895 globalDumpedIndex = dumpedIndex + obsOffset(taskIndex)
896 outBmatrix(globalDumpedIndex, :, :) = tempoBmatrix(:, :)
897 outLats(globalDumpedIndex) = latitude
898 outLons(globalDumpedIndex) = longitude
899 end if
900 end do
901 else
902 tag = 3 * mmpi_myID
903 call rpn_comm_send(latitude, 1, 'mpi_real8', 0, tag, 'GRID', ierr)
904 call rpn_comm_send(longitude, 1, 'mpi_real8', 0, tag + 1, 'GRID', ierr)
905 if (columnIndex > 0) then
906 call rpn_comm_send(bEns(columnIndex, :, :), nkgdim*nkgdim, 'mpi_real8', 0, tag + 2, 'GRID', ierr)
907 end if
908 end if
909 end do
910
911 if (mmpi_myId == 0) then
912 if (doAveraging) then
913 write(nulmat) real(sum(outLats(:))/countDumpedMpiGlobal, kind=4), real(sum(outLons(:))/countDumpedMpiGlobal, kind=4), &
914 sum(outBmatrix(:, :, :),dim=1)/countDumpedMpiGlobal
915 else
916 do dumpedIndex = 1, countDumpedMpiGlobal
917 write(nulmat) real(outLats(dumpedIndex), kind=4), real(outLons(dumpedIndex), kind=4), outBmatrix(dumpedIndex, :, :)
918 end do
919 end if
920 ierr = fclos(nulmat)
921 deallocate(outLats)
922 deallocate(outLons)
923 deallocate(outBmatrix)
924 deallocate(tempoBmatrix)
925 end if
926 if (allocated(listColumnDumped)) deallocate(listColumnDumped)
927 end subroutine dumpBmatrices
928
929 !--------------------------------------------------------------------------
930 ! bmat1D_bSqrtHi
931 !--------------------------------------------------------------------------
932 subroutine bmat1D_bSqrtHi(controlVector_in, column, obsSpaceData)
933 !
934 !:Purpose: HI component of B square root in 1DVar mode
935 !
936 implicit none
937
938 ! Arguments:
939 real(8), intent(in) :: controlVector_in(cvDim_mpilocal)
940 type(struct_columnData), intent(inout) :: column
941 type(struct_obs), intent(in) :: obsSpaceData
942
943 ! Locals:
944 integer :: headerIndex, latitudeBandIndex(1), varIndex, columnIndex
945 real(8), pointer :: currentColumn(:)
946 real(8), allocatable :: oneDProfile(:)
947 real(8) :: latitude
948 integer :: surfaceType, offset
949
950 if (mmpi_myid == 0) write(*,*) 'bmat1D_bsqrtHi: starting'
951 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
952
953 if (.not. initialized) then
954 if (mmpi_myid == 0) write(*,*) 'bmat1D_bsqrtHi: 1Dvar B matrix not initialized'
955 return
956 end if
957 allocate(oneDProfile(nkgdim))
958
959 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,oneDProfile,offset,varIndex,currentColumn,latitude,surfaceType,latitudeBandIndex)
960 do columnIndex = 1, var1D_validHeaderCount
961 headerIndex = var1D_validHeaderIndex(columnIndex)
962 latitude = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex) !radian
963 surfaceType = tvs_ChangedStypValue(obsSpaceData, headerIndex)
964 if (surfaceType == 1) then !Sea
965 latitudeBandIndex = minloc( abs( latitude - latSea(:)) )
966 oneDProfile(:) = matmul(bSqrtSea(latitudeBandIndex(1), :, :), controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim))
967 else ! Land or Sea Ice
968 latitudeBandIndex = minloc( abs( latitude - latLand(:)) )
969 oneDProfile(:) = matmul(bSqrtLand(latitudeBandIndex(1), :, :), controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim))
970 end if
971 offset = 0
972 do varIndex = 1, bmat1D_numIncludeAnlVar
973 currentColumn => col_getColumn(column, headerIndex, varName_opt=bmat1D_includeAnlVar(varIndex))
974 currentColumn(:) = currentColumn(:) + oneDProfile(offset+1:offset+size(currentColumn))
975 offset = offset + size(currentColumn)
976 end do
977 if (offset /= nkgdim) then
978 write(*,*) 'bmat1D_bsqrtHi: offset, nkgdim', offset, nkgdim
979 call utl_abort('bmat1D_bSqrtHi: inconsistency between Bmatrix and statevector size')
980 end if
981 end do
982 !$OMP END PARALLEL DO
983
984 deallocate(oneDProfile)
985 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
986 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtHi: done'
987
988 end subroutine bmat1D_bSqrtHi
989
990 !--------------------------------------------------------------------------
991 ! bmat1D_bSqrtHiAd
992 !--------------------------------------------------------------------------
993 subroutine bmat1D_bSqrtHiAd(controlVector_in, column, obsSpaceData)
994 !
995 !:Purpose: HI component of B square root adjoint in 1DVar mode
996 !
997 implicit none
998
999 ! Arguments:
1000 real(8), intent(inout) :: controlVector_in(cvDim_mpilocal)
1001 type(struct_columnData), intent(inout) :: column
1002 type (struct_obs), intent(in) :: obsSpaceData
1003
1004 ! Locals:
1005 integer :: headerIndex, latitudeBandIndex(1), varIndex, columnIndex
1006 real(8), pointer :: currentColumn(:)
1007 real(8), allocatable :: oneDProfile(:)
1008 real(8) :: latitude
1009 integer :: surfaceType, offset
1010
1011 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtHiAd: starting'
1012 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
1013 if (.not. initialized) then
1014 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtHiAd: 1dvar Bmatrix not initialized'
1015 return
1016 end if
1017 allocate(oneDProfile(nkgdim))
1018
1019 controlVector_in(:) = 0.d0
1020
1021 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,oneDProfile,offset,varIndex,currentColumn,latitude,surfaceType,latitudeBandIndex)
1022 do columnIndex = 1, var1D_validHeaderCount
1023 headerIndex = var1D_validHeaderIndex(columnIndex)
1024 offset = 0
1025 do varIndex = 1, bmat1D_numIncludeAnlVar
1026 currentColumn => col_getColumn(column, headerIndex, varName_opt=bmat1D_includeAnlVar(varIndex))
1027 oneDProfile(offset+1:offset+size(currentColumn)) = currentColumn(:)
1028 offset = offset + size(currentColumn)
1029 end do
1030 if (offset /= nkgdim) then
1031 write(*,*) 'bmat1D_bSqrtHiAd: offset, nkgdim', offset, nkgdim
1032 call utl_abort('bmat1D_bSqrtHiAd: inconsistency between Bmatrix and statevector size')
1033 end if
1034 latitude = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex) !radians
1035 surfaceType = tvs_ChangedStypValue(obsSpaceData, headerIndex)
1036 if (surfaceType == 1) then !Sea
1037 latitudeBandIndex = minloc( abs( latitude - latSea(:)) )
1038 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) = &
1039 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) + &
1040 matmul(bSqrtSea(latitudeBandIndex(1),:,:), oneDProfile)
1041 else ! Land or Sea Ice
1042 latitudeBandIndex = minloc( abs( latitude - latLand(:)) )
1043 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) = &
1044 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) + &
1045 matmul(bSqrtLand(latitudeBandIndex(1),:,:), oneDProfile)
1046 end if
1047 end do
1048 !$OMP END PARALLEL DO
1049
1050 deallocate(oneDProfile)
1051 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
1052 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtHiAd: done'
1053
1054 end subroutine bmat1D_bSqrtHiAd
1055
1056 !--------------------------------------------------------------------------
1057 ! bmat1D_bSqrtEns
1058 !--------------------------------------------------------------------------
1059 subroutine bmat1D_bSqrtEns(controlVector_in, column)
1060 !
1061 !:Purpose: Ensemble component of B square root in 1DVar mode
1062 !
1063 implicit none
1064
1065 ! Arguments:
1066 real(8), intent(in) :: controlVector_in(cvDim_mpilocal)
1067 type(struct_columnData), intent(inout) :: column
1068
1069 ! Locals:
1070 integer :: headerIndex, varIndex, columnIndex
1071 real(8), pointer :: currentColumn(:)
1072 real(8), allocatable :: oneDProfile(:)
1073 integer :: offset
1074
1075 allocate(oneDProfile(nkgdim))
1076
1077 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,oneDProfile,offset,varIndex,currentColumn)
1078 do columnIndex = 1, var1D_validHeaderCount
1079 headerIndex = var1D_validHeaderIndex(columnIndex)
1080 oneDProfile(:) = matmul(bSqrtEns(columnIndex, :, :), controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim))
1081 offset = 0
1082 do varIndex = 1, bmat1D_numIncludeAnlVar
1083 currentColumn => col_getColumn(column, headerIndex, varName_opt=bmat1D_includeAnlVar(varIndex))
1084 currentColumn(:) = currentColumn(:) + oneDProfile(offset+1:offset+size(currentColumn))
1085 offset = offset + size(currentColumn)
1086 end do
1087 if (offset /= nkgdim) then
1088 write(*,*) 'bmat1D_bsqrtEns: offset, nkgdim', offset, nkgdim
1089 call utl_abort('bmat1D_bSqrtEns: inconsistency between Bmatrix and statevector size')
1090 end if
1091 end do
1092 !$OMP END PARALLEL DO
1093
1094 deallocate(oneDProfile)
1095
1096 end subroutine bmat1D_bSqrtEns
1097
1098 !--------------------------------------------------------------------------
1099 ! bmat1D_bSqrtEnsAd
1100 !--------------------------------------------------------------------------
1101 subroutine bmat1D_bSqrtEnsAd(controlVector_in, column)
1102 !
1103 !:Purpose: Ensemble component of B square root in 1DVar mode
1104 !
1105 implicit none
1106
1107 ! Arguments:
1108 real(8), intent(inout) :: controlVector_in(cvDim_mpilocal)
1109 type(struct_columnData), intent(inout) :: column
1110
1111 ! Locals:
1112 integer :: headerIndex, varIndex, columnIndex
1113 real(8), pointer :: currentColumn(:)
1114 real(8), allocatable :: oneDProfile(:)
1115 integer :: offset
1116
1117 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtEnsAd: starting'
1118 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
1119 if (.not. initialized) then
1120 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtEnsAd: 1dvar Bmatrix not initialized'
1121 return
1122 end if
1123 allocate(oneDProfile(nkgdim))
1124
1125 controlVector_in(:) = 0.d0
1126
1127 !$OMP PARALLEL DO PRIVATE (columnIndex,headerIndex,oneDProfile,offset,varIndex,currentColumn)
1128 do columnIndex = 1, var1D_validHeaderCount
1129 headerIndex = var1D_validHeaderIndex(columnIndex)
1130 offset = 0
1131 do varIndex = 1, bmat1D_numIncludeAnlVar
1132 currentColumn => col_getColumn(column, headerIndex, varName_opt=bmat1D_includeAnlVar(varIndex))
1133 oneDProfile(offset+1:offset+size(currentColumn)) = currentColumn(:)
1134 offset = offset + size(currentColumn)
1135 end do
1136 if (offset /= nkgdim) then
1137 write(*,*) 'bmat1D_bSqrtHiAd: offset, nkgdim', offset, nkgdim
1138 call utl_abort('bmat1D_bSqrtEnsAd: inconsistency between Bmatrix and statevector size')
1139 end if
1140 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) = &
1141 controlVector_in(1+(columnIndex-1)*nkgdim:columnIndex*nkgdim) + &
1142 matmul(bSqrtEns(columnIndex,:,:), oneDProfile)
1143 end do
1144 !$OMP END PARALLEL DO
1145
1146 deallocate(oneDProfile)
1147 if (mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
1148 if (mmpi_myid == 0) write(*,*) 'bmat1D_bSqrtEnsAd: done'
1149
1150 end subroutine bmat1D_bSqrtEnsAd
1151
1152 !--------------------------------------------------------------------------
1153 ! bmat1D_sqrtB
1154 !--------------------------------------------------------------------------
1155 subroutine bmat1D_sqrtB(controlVector, cvdim, column, obsSpaceData)
1156 !
1157 !:Purpose: To transform model state from control-vector space to grid-point
1158 ! space.
1159 !
1160 implicit none
1161
1162 ! Arguments:
1163 integer, intent(in) :: cvdim
1164 real(8), intent(in) :: controlVector(cvdim)
1165 type(struct_columnData), intent(inout) :: column
1166 type(struct_obs), intent(in) :: obsSpaceData
1167
1168 ! Locals:
1169 integer :: bmatIndex
1170 real(8), pointer :: subVector(:)
1171
1172 !
1173 !- 1. Compute the analysis increment
1174 !
1175 call col_zero(column)
1176
1177 bmat_loop: do bmatIndex = 1, numBmat
1178 if ( .not. bmatActive(bmatIndex) ) cycle bmat_loop
1179 subVector => cvm_getSubVector( controlVector, bmatLabelList(bmatIndex) )
1180
1181 call utl_tmg_start(50, '--Bmatrix')
1182 select case( trim(bmatTypeList(bmatIndex)) )
1183 case ('HI')
1184 !- 1.1 Time-Mean Homogeneous and Isotropic...
1185 call utl_tmg_start(52, '----B_HI_TL')
1186 call bmat1D_bsqrtHi(subVector, & ! IN
1187 column, & ! OUT
1188 obsspacedata ) ! IN
1189 call utl_tmg_stop(52)
1190 case ('ENS')
1191 !- 1.2 Ensemble based
1192 call utl_tmg_start(57, '----B_ENS_TL')
1193 call bmat1D_bsqrtEns(subVector, & ! IN
1194 column) ! OUT
1195 call utl_tmg_stop(57)
1196 case default
1197 call utl_abort( 'bmat1D_sqrtB: requested bmatrix type does not exist ' // trim(bmatTypeList(bmatIndex)) )
1198 end select
1199 call utl_tmg_stop(50)
1200
1201 end do bmat_loop
1202
1203 end subroutine bmat1D_sqrtB
1204
1205 !--------------------------------------------------------------------------
1206 ! bmat1D_sqrtBT
1207 !--------------------------------------------------------------------------
1208 subroutine bmat1D_sqrtBT(controlVector, cvdim, column, obsSpaceData)
1209 !
1210 !:Purpose: To transform model state from grid-point space to
1211 ! error-covariance space.
1212 !
1213 implicit none
1214
1215 ! Arguments:
1216 integer, intent(in) :: cvdim
1217 real(8), intent(in) :: controlVector(cvdim)
1218 type(struct_columnData), intent(inout) :: column
1219 type(struct_obs), intent(in) :: obsSpaceData
1220
1221 ! Locals:
1222 integer :: bmatIndex
1223 real(8), pointer :: subVector(:)
1224
1225 ! Process components in opposite order as forward calculation
1226 bmat_loop: do bmatIndex = numBmat, 1, -1
1227 if ( .not. bmatActive(bmatIndex) ) cycle bmat_loop
1228 subVector => cvm_getSubVector( controlVector, bmatLabelList(bmatIndex) )
1229
1230 call utl_tmg_start(50, '--Bmatrix')
1231 select case( trim(bmatTypeList(bmatIndex)) )
1232 case ('HI')
1233 !- Time-Mean Homogeneous and Isotropic...
1234 call utl_tmg_start(53, '----B_HI_AD')
1235 call bmat1D_bsqrtHiAd(subvector, & ! IN
1236 column, & ! OUT
1237 obSSpaceData ) ! IN
1238 call utl_tmg_stop(53)
1239 case ('ENS')
1240 !- Ensemble based
1241 call utl_tmg_start(61, '----B_ENS_AD')
1242 call bmat1D_bsqrtEnsAd(subvector, & ! IN
1243 column ) ! OUT
1244 call utl_tmg_stop(61)
1245 case default
1246 call utl_abort( 'bmat1D_sqrtBT: requested bmatrix type does not exist ' // trim(bmatTypeList(bmatIndex)) )
1247 end select
1248 call utl_tmg_stop(50)
1249
1250 end do bmat_loop
1251
1252 end subroutine bmat1D_sqrtBT
1253
1254 !--------------------------------------------------------------------------
1255 ! bmat1D_get1DVarIncrement
1256 !--------------------------------------------------------------------------
1257 subroutine bmat1D_get1DVarIncrement(incr_cv, column, columnTrlOnAnlIncLev, &
1258 obsSpaceData, nvadim_mpilocal)
1259 !
1260 !:Purpose: to compute 1Dvar increment from control vector
1261 !
1262 implicit none
1263
1264 ! Arguments:
1265 real(8), intent(in) :: incr_cv(:)
1266 type(struct_columnData), intent(inout) :: column
1267 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1268 type(struct_obs), intent(in) :: obsSpaceData
1269 integer, intent(in) :: nvadim_mpilocal
1270
1271 ! compute increment from control vector (multiply by B^1/2)
1272 call bmat1D_sqrtB(incr_cv, nvadim_mpilocal, column, obsSpaceData)
1273 call cvt_transform(column, 'ZandP_tl', columnTrlOnAnlIncLev)
1274
1275 end subroutine bmat1D_get1DVarIncrement
1276
1277 !--------------------------------------------------------------------------
1278 ! bmat1D_Finalize
1279 !--------------------------------------------------------------------------
1280 subroutine bmat1D_Finalize()
1281 !
1282 !:Purpose: to deallocate memory used by internal module structures
1283 !
1284 implicit none
1285
1286 if (initialized) then
1287 if (allocated(bSqrtLand)) then
1288 deallocate( bSqrtLand )
1289 deallocate( bSqrtSea )
1290 deallocate( latLand, lonLand, latSea, lonSea )
1291 end if
1292 if (allocated(bSqrtEns)) deallocate( bSqrtEns )
1293 call var1D_finalize()
1294 end if
1295
1296 end subroutine bmat1D_Finalize
1297
1298end module bMatrix1DVar_mod
