1module bCovarSetupChem_mod
2 ! MODULE bCovarSetupChem_mod (prefix='bcsc' category='6. High-level data objects')
3 !
4 !:Purpose: Contains routines for the reading and preparation of
5 ! background-error covariance elements. Correlation matrices
6 ! based on horizontally homogeneous/isotropic correlations.
7 !
8 !:Comments:
9 !
10 ! 1. Covariances uncoupled from weather variable.
11 !
12 ! 2. Handles univariate and multivariate covariances.
13 ! See routines bcsc_readcorns2 and bcsc_sucorns2.
14 !
15 ! 3. For multiple univariate variables (or univarite blocks of one to
16 ! multiple variables), one can alternatively have multiple sets of
17 ! covariance matrices within this module instead of a single covariance
18 ! matrix setup (similarly to what was done for corvert*).
19 !
20
21 ! Public Subroutines:
22 !
23 ! bcsc_setupCH: Must be called first. Sets of background covariance
24 ! matrix (and balance operators if any are eventually
25 ! added)
26 ! bcsc_getCovarCH: Provides covariances and related
27 ! bcsc_getScaleFactor : Provides std dev scaling factor
28 ! elements for background check and obs operators.
29 ! bcsc_finalize Deallocate internal module arrays.
30 ! bcsc_StatsExistForVarName: Determine is Covar available for specified variable
31 ! bcsc_getBgStddev: Interpolation background error std dev to obs location
32 ! bcsc_retrieveBgStddev: Retrieve previously saved background stddev
33 ! profiles in bgStddev from the header index.
34 ! bcsc_addBgStddev: Add background stddev profiles (and inverse) to
35 ! bgStddev which can be retrieved later using a header index.
36
37 use midasMpi_mod
38 use MathPhysConstants_mod
39 use earthConstants_mod
40 use obsSubSpaceData_mod
41 use gridStateVector_mod
42 use globalSpectralTransform_mod
43 use horizontalCoord_mod
44 use verticalCoord_mod
45 use varNameList_mod
46 use utilities_mod
47 use calcHeightAndPressure_mod
48
49 implicit none
50 save
51 private
52
53 ! public procedures
54 ! ------------------
55
56 public :: bcsc_setupCH,bcsc_finalize
57 public :: bcsc_getCovarCH
58 public :: bcsc_getScaleFactor
59 public :: bcsc_StatsExistForVarName
60 public :: bcsc_retrieveBgStddev,bcsc_addBgStddev,bcsc_getBgStddev
61
62 ! public types
63 ! ------------
64
65 public :: struct_bcsc_bgStats
66
67 ! module shared variables
68 ! -----------------------
69
70 integer :: nlev_M,nlev_T
71 integer :: gstID
72 real(8), allocatable :: rlatr(:),rlongr(:)
73 type(struct_vco),pointer :: vco_anl
74
75 character(len=15) :: bcsc_mode
76
77 integer, external :: get_max_rss
78 integer :: nulbgst=0
79
80 ! Background error covariance matrix elements.
81 ! One could add an additional dimension to corns
82 ! for separate block-univariate correlation matrices.
83 ! This would also permit merging of bmatrixhi_mod and bmatrixchem_mod
84 ! into one module.
85
86 real(8),allocatable :: stddev(:,:,:)
87 real(8),allocatable :: rstddev(:,:)
88
89 ! Parameters of the NAMBCHM namelist
90 integer :: ntrunc ! spectral truncation
91 real(8) :: rpor(vnl_numvarmax) ! horizontal localization distance (Gaussian)
92 real(8) :: rvloc(vnl_numvarmax) ! vertical localization distance (GC)
93 real(8) :: scaleFactor(vnl_numvarmax,vco_maxNumLevels) ! variable and level dependent scale factor applied to variances
94 integer :: numModeZero ! number of eigenmodes to set to zero
95 logical :: ReadWrite_sqrt ! choose to read and/or write sqrt of correlations
96 logical :: getPhysSpaceStats ! choose to calculate/save physical space cov (stddev, corverti)
97 logical :: getPhysSpaceHCorrel ! calculate correlation lengths from spectral cov (needed for some CH obs operator settings)
98 logical :: WritePhysSpaceStats ! choose to output physical space stats in 'bCovarSetupChem_out.fst'
99 character(len=4) :: stddevMode ! can be 'GD2D', 'GD3D' or 'SP2D'
100 character(len=4) :: IncludeAnlVarKindCH(vnl_numvarmax) ! list of CH variable names to consider
101 character(len=4) :: CrossCornsVarKindCH(vnl_numvarmax) ! not sure what this is for...
102 character(len=20) :: TransformVarKindCH ! name of variable transform to apply to chemistry variables
103
104 ! Square root of scaleFactor
105 real(8) :: scaleFactor_stddev(vnl_numvarmax,vco_maxNumLevels)
106
107 real(8), parameter :: zps = 101000.D0 ! Reference surface pressure
108
109 ! module structures
110 ! -----------------
111
112 type :: struct_bcsc_bgStats
113
114 ! Structure storing background error stats
115 ! and related elements
116
117 ! logical indicating if stats are available
118 logical :: initialized=.false.
119
120 integer :: numvar3d ! number of 3D fields
121 integer :: numvar2d ! number of 2D fields
122 integer :: ni ! number of latitudes
123 integer :: nj ! number of longitudes
124 integer :: nlev ! number of vertical levels
125
126 ! Total number of elements over all verticallevels and
127 ! variables (varNameList)
128 integer :: nkgdim
129
130 integer :: ntrunc ! spectral dimension
131 character(len=4), allocatable :: varNameList(:) ! list of variable names
132
133 integer, allocatable :: nsposit(:) ! start positions of fields
134 real(8), allocatable :: lat(:) ! grid lat in radians
135 real(8), allocatable :: lon(:) ! grid lon in radians
136 real(8), allocatable :: vlev(:) ! vertical levels
137 real(8), allocatable :: corns(:,:,:) ! spectral space correlation matrix
138
139 ! Phys. Space vertical correlation matrix
140 real(8), allocatable :: corvert(:,:,:)
141
142 real(8), allocatable :: corverti(:,:,:) ! Inverse of 'corvert'
143
144 real(8), allocatable :: stddev(:,:,:) ! error std dev
145 real(8), allocatable :: hcorrlen(:,:) ! horizontal correlation lengths
146
147 end type struct_bcsc_bgStats
148
149 ! Assigned type variables
150 ! -----------------------
151
152 type(struct_bcsc_bgStats) :: bgStats
153 type(struct_oss_obsdata) :: bgStddev ! Arrays for background error
154 ! std dev in obs space
155
156 !*************************************************************************
157
158 contains
159
160 !--------------------------------------------------------------------------
161 ! bcsc_setupCH
162 !--------------------------------------------------------------------------
163 subroutine bcsc_setupCH(hco_in,vco_in,covarNeeded,mode)
164 !
165 !:Purpose: Set up for constituents static background error covariances.
166 !
167 implicit none
168
169 ! Arguments:
170 type(struct_hco), pointer, intent(in) :: hco_in
171 type(struct_vco), pointer, intent(in) :: vco_in
172 logical, intent(out) :: covarNeeded
173 character(len=*), intent(in) :: mode ! 'Analysis' or 'BackgroundCheck'
174
175 ! Locals:
176 integer :: nulnam, ierr, fnom, fclos
177 integer :: varIndex,nChmVars,varIndex2
178 character(len=4) :: BchmVars(vnl_numvarmax)
179 real(8), pointer :: pressureProfile_T(:)
180
181 NAMELIST /NAMBCHM/ntrunc,rpor,rvloc,scaleFactor,numModeZero,ReadWrite_sqrt, &
182 stddevMode,IncludeAnlVarKindCH,getPhysSpaceHCorrel, &
183 CrossCornsVarKindCH,WritePhysSpaceStats, &
184 TransformVarKindCH,getPhysSpaceStats
185
186 write(*,*) 'Started bcsc_setupCH'
187
188 ! First check if there are any CH fields
189
190 covarNeeded = .true.
191 varIndex2=0
192 do varIndex = 1, vnl_numvarmax
193 if (gsv_varExist(varName = vnl_varNameList(varIndex))) then
194 if (vnl_varKindFromVarname(vnl_varNameList(varIndex)) == 'CH') then
195 varIndex2 = 1
196 exit
197 end if
198 end if
199 end do
200 if (varIndex2 == 0) then
201 ! Assume there is no need for Bchm
202 covarNeeded = .false.
203 return
204 end if
205
206 bgStats%numvar3d = 0
207 bgStats%numvar2d = 0
208
209 allocate(bgStats%varNameList(vnl_numvarmax))
210 bgStats%varNameList(:) = ''
211 allocate(bgStats%nsposit(vnl_numvarmax+1))
212 bgStats%nsposit(1) = 1
213
214 if ( trim(mode) == 'Analysis' .or. trim(mode) == 'BackgroundCheck') then
215 bcsc_mode = trim(mode)
216 if(mmpi_myid == 0) write(*,*)
217 if(mmpi_myid == 0) write(*,*) 'bcsc_setupCH: Mode activated = ', &
218 trim(bcsc_mode)
219 else
220 write(*,*)
221 write(*,*) 'mode = ', trim(mode)
222 call utl_abort('bcsc_setupCH: unknown mode')
223 end if
224
225 ! Initialization of namelist NAMBCHM parameters
226
227 ntrunc=108
228 rpor(:)=3000.D3
229 rvloc(:)=4.0D0
230 scaleFactor(:,:) = 0.0d0
231 numModeZero = 0
232 ReadWrite_sqrt = .false.
233 WritePhysSpaceStats = .false.
234 getPhysSpaceHCorrel = .false.
235 getPhysSpaceStats = .false.
236 stddevMode = 'GD3D'
237 IncludeAnlVarKindCH(:) = ''
238 CrossCornsVarKindCH(:) = ''
239 TransformVarKindCH = ''
240
241 ! Read namelist input
242
243 nulnam = 0
244 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
245 read(nulnam,nml=NAMBCHM,iostat=ierr)
246 if(ierr /= 0) call utl_abort('bcsc_setupCH: Error reading namelist')
247 if(mmpi_myid == 0) write(*,nml=NAMBCHM)
248 ierr = fclos(nulnam)
249
250 ! Set BchmVars
251 nChmVars=0
252 BChmVars(:)=''
253 if (trim(IncludeAnlVarKindCH(1)) == '') then
254 do varIndex = 1, vnl_numvarmax
255 if (.not. gsv_varExist(varName = vnl_varNameList(varIndex))) cycle
256 if (vnl_varKindFromVarname(vnl_varNameList(varIndex)) /= 'CH') cycle
257 nChmVars = nChmVars+1
258 BchmVars(nChmVars) = trim(vnl_varNameList(varIndex))
259 end do
260 else
261 do varIndex = 1, vnl_numvarmax
262 if (.not. gsv_varExist(varName = vnl_varNameList(varIndex))) cycle
263 if (vnl_varKindFromVarname(vnl_varNameList(varIndex)) /= 'CH') cycle
264 do varIndex2 = 1, vnl_numvarmax
265 if (trim(IncludeAnlVarKindCH(varIndex2)) == '') exit
266 if (trim(vnl_varNameList(varIndex)) == &
267 trim(IncludeAnlVarKindCH(varIndex2))) then
268 nChmVars = nChmVars+1
269 BchmVars(nChmVars)= trim(vnl_varNameList(varIndex))
270 exit
271 end if
272 end do
273 end do
274 end if
275
276 if (nChmVars == 0) then
277 if(mmpi_myid == 0) then
278 write(*,*) 'Size of BchmVars is zero. B matrix for CH family ', &
279 'not produced.'
280 write(*,*) 'No chemical assimilation to be performed.'
281 write(*,*) 'Completed bcsc_setupCH'
282 end if
283 covarNeeded = .false.
284 return
285 end if
286
287 ! Set vertical dimensions
288
289 vco_anl => vco_in
290 nLev_M = vco_anl%nlev_M
291 nLev_T = vco_anl%nlev_T
292
293 ! Find the 3D variables (within NAMSTATE namelist)
294
295 do varIndex = 1, vnl_numvarmax3D
296 if (gsv_varExist(varName=vnl_varNameList3D(varIndex)) .and. &
297 any(trim(vnl_varNameList3D(varIndex))==BchmVars(1:nChmVars)) ) then
298
299 if (vnl_varKindFromVarname(vnl_varNameList3D(varIndex)) /= 'CH') cycle
300
301 bgStats%numvar3d = bgStats%numvar3d + 1
302 bgStats%nsposit(bgStats%numvar3d+1)= &
303 bgStats%nsposit(bgStats%numvar3d)+nLev_T
304 bgStats%varNameList(bgStats%numvar3d)= &
305 vnl_varNameList3D(varIndex)
306 end if
307 end do
308
309 ! Find the 2D variables (within NAMSTATE namelist)
310
311 do varIndex = 1, vnl_numvarmax2D
312 if (gsv_varExist(varName=vnl_varNameList2D(varIndex)) .and. &
313 any(trim(vnl_varNameList2D(varIndex)) == BchmVars(1:nChmVars)) ) then
314
315 if (vnl_varKindFromVarname(vnl_varNameList2D(varIndex)) /= 'CH') cycle
316 bgStats%numvar2d = bgStats%numvar2d + 1
317 bgStats%nsposit(bgStats%numvar3d+bgStats%numvar2d+1) = &
318 bgStats%nsposit(bgStats%numvar3d+bgStats%numvar2d)+1
319 bgStats%varNameList(bgStats%numvar2d) = vnl_varNameList2D(varIndex)
320 end if
321 end do
322
323 if (bgStats%numvar3d+bgStats%numvar2d == 0) then
324 if (mmpi_myid == 0) then
325 write(*,*) 'B matrix for CH family not produced.'
326 write(*,*) 'No chemical assimilation to be performed.'
327 write(*,*) 'Completed bcsc_setupCH'
328 end if
329 covarNeeded = .false.
330 return
331 else if (mmpi_myid == 0) then
332 if (bgStats%numvar3d > 0) then
333 write(*,*) 'bcsc_setupCH: Number of 3D variables', &
334 bgStats%numvar3d,bgStats%varNameList(1:bgStats%numvar3d)
335 end if
336 if (bgStats%numvar2d > 0) then
337 write(*,*) 'bcsc_setupCH: Number of 2D variables', &
338 bgStats%numvar2d,bgStats%varNameList(bgStats%numvar3d+1: &
339 bgStats%numvar3d+bgStats%numvar2d)
340 end if
341 end if
342
343 bgStats%nkgdim = &
344 bgStats%nsposit(bgStats%numvar3d+bgStats%numvar2d+1)-1
345
346 ! Scalefactors must be > 0 until ensembles for constituents can be used.
347
348 if (bgStats%numvar3d > 0) then
349 if (any(scaleFactor(1:bgStats%numvar3d,1:nLev_T) <= 0.0D0)) then
350 write(*,*) 'Scalefactors: ',scaleFactor(1:bgStats%numvar3d,1:nLev_T)
351 call utl_abort('bcsc_setupCH: Scalefactors values must be > 0 for now.')
352 end if
353 end if
354 if (bgStats%numvar2d > 0) then
355 if (any(scaleFactor(1:bgStats%numvar2d,1) <= 0.0D0)) then
356 write(*,*) 'Scalefactors: ',scaleFactor(1:bgStats%numvar3d,1:nLev_T)
357 call utl_abort('bcsc_setupCH: Scalefactors values must be > 0 for now.')
358 end if
359 end if
360
361 ! Set scalefactor_stddev
362
363 scaleFactor_stddev(1:bgStats%numvar3d+bgStats%numvar2d, &
364 1:max(nLev_M,nLev_T))=0.0d0
365
366 where (scaleFactor(1:bgStats%numvar3d+bgStats%numvar2d, &
367 1:max(nLev_M,nLev_T)) > 0.0d0)
368
369 scaleFactor_stddev(1:bgStats%numvar3d+bgStats%numvar2d, &
370 1:max(nLev_M,nLev_T)) = sqrt(scaleFactor(1:bgStats%numvar3d &
371 + bgStats%numvar2d,1:max(nLev_M,nLev_T)))
372
373 end where
374
375 bgStats%ni = hco_in%ni
376 bgStats%nj = hco_in%nj
377 bgStats%ntrunc = ntrunc
378 gstID = gst_setup(bgStats%ni,bgStats%nj,bgStats%ntrunc, &
379 bgStats%nkgdim)
380
381 if (allocated(bgStats%lat)) deallocate(bgStats%lat)
382 if (allocated(bgStats%lon)) deallocate(bgStats%lon)
383 allocate(bgStats%lat(bgStats%nj))
384 allocate(bgStats%lon(bgStats%ni+1))
385
386 bgStats%lat(1:bgStats%nj) = hco_in%lat(1:bgStats%nj)
387 bgStats%lon(1:bgStats%ni) = hco_in%lon(1:bgStats%ni)
388 bgStats%lon(bgStats%ni+1) = 2*MPC_PI_R8
389
390 ! Assign sizes and transfer some fields
391
392 if ( bgStats%numvar3d > 0 ) then
393 bgStats%nlev = nlev_T
394 else if (bgStats%numvar2d > 0 ) then
395 bgStats%nlev = 1
396 else
397 bgStats%nlev = 0
398 end if
399
400 allocate(stddev(bgStats%ni+1, bgStats%nj, bgStats%nkgdim))
401 allocate(rstddev(bgStats%nkgdim, 0:bgStats%ntrunc))
402
403 allocate(bgStats%corns(bgStats%nkgdim,bgStats%nkgdim, &
404 0:bgStats%ntrunc))
405 allocate(bgStats%stddev(bgStats%ni+1,bgStats%nj, &
406 bgStats%nkgdim ))
407 if ( bgStats%nlev > 1 ) then
408 allocate(bgStats%vlev(bgStats%nlev))
409 if (getPhysSpaceStats .or. trim(bcsc_mode) == 'BackgroundCheck') then
410 allocate(bgStats%corvert(bgStats%nlev,bgStats%nlev, &
411 bgStats%numvar3d))
412 end if
413 if (getPhysSpaceStats) then
414 allocate(bgStats%corverti(bgStats%nlev,bgStats%nlev,&
415 bgStats%numvar3d))
416 end if
417 end if
418
419 ! Get vertical levels
420
421 if (bgStats%nlev > 1) then
422 call czp_fetch1DLevels(vco_anl, zps, profT_opt=pressureProfile_T)
423 bgStats%vlev(1:bgStats%nlev) = pressureProfile_T(1:bgStats%nlev)
424 else if (bgStats%nlev == 1) then
425 bgStats%vlev(1)=zps
426 end if
427
428 ! Read covar stats, scale standard deviations, and apply localization
429 ! to vertical correlation matrices in horizontal spectral space
430
431 call bcsc_rdstats
432 bgStats%stddev(:,:,:)=stddev(:,:,:)
433 if (allocated(stddev)) deallocate(stddev)
434
435 ! Generate or read correlation matrix square roots
436
437 call bcsc_sucorns2
438
439 if(mmpi_myid == 0) write(*,*) 'Completed bcsc_setupCH'
440
441 bgStats%initialized = .true.
442 if (allocated(rstddev)) deallocate(rstddev)
443
444 end subroutine bcsc_setupCH
445
446 !--------------------------------------------------------------------------
447 ! bcsc_getScaleFactor
448 !--------------------------------------------------------------------------
449 subroutine bcsc_getScaleFactor(scaleFactorOut)
450 !
451 !:Purpose: To set scaling factors for background error std. dev.
452 !
453 implicit none
454
455 ! Arguments:
456 real(8), intent(out) :: scaleFactorOut(:,:) ! Error std. dev. scale factor
457
458 ! Locals:
459 integer :: levelIndex,varIndex
460
461 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
462 do levelIndex = 1, bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
463 scaleFactorOut(varIndex,levelIndex) = scaleFactor_stddev(varIndex,levelIndex)
464 end do
465 end do
466
467 end subroutine bcsc_getScaleFactor
468
469 !--------------------------------------------------------------------------
470 ! bcsc_rdstats
471 !--------------------------------------------------------------------------
472 subroutine bcsc_rdstats
473 !
474 !:Purpose: To read chemical constituents background stats file.
475 !
476 implicit none
477
478 ! Locals:
479 integer :: ierr, fnom, fstouv, fstfrm, fclos
480 logical :: lExists
481 character(len=12) :: bFileName = './bgchemcov'
482 type(struct_vco),pointer :: vco_file => null()
483
484 inquire(file=bFileName,exist=lExists)
485 if ( lexists ) then
486 ierr = fnom(nulbgst,bFileName,'RND+OLD+R/O',0)
487 if ( ierr == 0 ) then
488 ierr = fstouv(nulbgst,'RND+OLD')
489 else
490 call utl_abort('bcsc_rdstats: Problem in opening the background ' // &
491 'chemical constituent stat file')
492 end if
493 else
494 call utl_abort('bcsc_rdstats: Background chemical constituent stat ' // &
495 'file is missing')
496 end if
497
498 ! check if analysisgrid and covariance file have the same vertical levels
499 call vco_SetupFromFile( vco_file, & ! OUT
500 bFileName ) ! IN
501 if (.not. vco_equal(vco_anl,vco_file)) then
502 call utl_abort('bcsc_rstats: vco from analysisgrid and chem cov file ' // &
503 'do not match')
504 end if
505
506 ! Read spectral space correlations
507
508 call bcsc_readcorns2
509
510 ! Read error standard deviations
511
512 call bcsc_rdstddev
513
514 ! Scale error standard deviations (and save to file if requested)
515
516 call bcsc_scalestd
517
518 ierr = fstfrm(nulbgst)
519 ierr = fclos(nulbgst)
520
521 end subroutine bcsc_rdstats
522
523 !--------------------------------------------------------------------------
524 ! bcsc_scalestd
525 !--------------------------------------------------------------------------
526 subroutine bcsc_scalestd
527 !
528 !:Purpose: To scale error standard-deviation values.
529 !
530 implicit none
531
532 ! Locals:
533 integer :: lonIndex, latIndex, varIndex, levelIndex, nlev, nulsig
534 integer :: ierr, fnom, fstouv, fstfrm, fclos
535
536 if (WritePhysSpaceStats .and. mmpi_myid == 0) then
537 nulsig = 0
538 ierr = fnom(nulsig,'bCovarSetupChem_out.fst','STD+RND',0)
539 ierr = fstouv(nulsig,'RND')
540 ierr = utl_fstecr(bgStats%vlev,-32,nulsig,0,0,0,1,1,bgStats%nlev, &
541 0,0,0,'X','PX','Pressure','X',0,0,0,0,5,.true.)
542 ierr = utl_fstecr(bgStats%lat,-32,nulsig,0,0,0,1,bgStats%nj,1, &
543 0,0,0,'X','^^','latitude','X',0,0,0,0,5,.true.)
544 ierr = utl_fstecr(bgStats%lon,-32,nulsig,0,0,0,bgStats%ni+1,1,1, &
545 0,0,0,'X','>>','longitude','X',0,0,0,0,5,.true.)
546 end if
547
548 do varIndex = 1,bgStats%numvar3d+bgStats%numvar2d
549 nlev=bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
550 do lonIndex = 1, bgStats%ni+1
551 do latIndex = 1, bgStats%nj
552 stddev(lonIndex,latIndex,bgStats%nsposit(varIndex): &
553 bgStats%nsposit(varIndex+1)-1) = &
554 scaleFactor_stddev(varIndex,1:nlev)* &
555 stddev(lonIndex,latIndex,bgStats%nsposit(varIndex): &
556 bgStats%nsposit(varIndex+1)-1)
557 end do
558 end do
559
560 if (WritePhysSpaceStats .and. mmpi_myid == 0) then
561 do levelIndex=1,nlev
562 ierr = utl_fstecr(stddev(1:bgStats%ni+1,1:bgStats%nj, &
563 bgStats%nsposit(varIndex)-1+levelIndex),-32,nulsig, &
564 0,0,0,bgStats%ni+1,bgStats%nj,1,levelIndex,0,nlev, &
565 'X',bgStats%varNameList(varIndex),'STDDEV','X',0,0,0,0, &
566 5,.true.)
567 end do
568 end if
569
570 end do
571
572 if (WritePhysSpaceStats .and. mmpi_myid == 0) then
573 ierr = fstfrm(nulsig)
574 ierr = fclos(nulsig)
575 end if
576
577 end subroutine bcsc_scalestd
578
579 !--------------------------------------------------------------------------
580 ! bcsc_readCorns2
581 !--------------------------------------------------------------------------
582 subroutine bcsc_readCorns2
583 !
584 !:Purpose: To read correlation information and to form the correlation
585 ! matrix.
586 !
587 !:Notes: Can read distinct block diagonal matrices with or without
588 ! cross-correlations.
589 !
590
591 ! Based on bhi_readcorns2.
592 implicit none
593
594 ! Locals:
595 integer :: jn,ierr,varIndex,varIndex2
596 integer :: jcol,jrow,jstart,jnum,jstart2,jnum2
597 real(8), allocatable, dimension(:) :: zstdsrc
598 real(8), allocatable, dimension(:,:) :: zcornssrc
599
600 ! Standard file variables
601 integer :: ini,inj,ink
602 integer :: ip1,ip2,ip3,idateo
603 character(len=2) :: cltypvar
604 character(len=4) :: clnomvar
605 character(len=4), allocatable :: clnomvarCrosscorns(:)
606 character(len=12) :: cletiket
607 integer :: fstinf
608
609 rstddev(:,:) = 0.0d0
610 bgStats%corns(:,:,:) = 0.0d0
611 if (any(CrossCornsVarKindCH(:) /= '')) then
612 allocate(clnomvarCrosscorns(bgStats%numvar3d+bgStats%numvar2d))
613 clnomvarCrosscorns(:)=''
614 end if
615
616 ! Read auto-correlations matrices
617
618 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
619
620 clnomvar = bgStats%varNameList(varIndex)
621 jstart = bgStats%nsposit(varIndex)
622 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
623 allocate(zcornssrc(jnum,jnum))
624 allocate(zstdsrc(jnum))
625
626 do jn = 0, bgStats%ntrunc
627
628 ! Looking for FST record parameters.
629
630 idateo = -1
631 cletiket = 'RSTDDEV'
632 ip1 = -1
633 ip2 = jn
634 ip3 = -1
635 cltypvar = 'X'
636
637 ierr = utl_fstlir(ZSTDSRC,nulbgst,INI,INJ,INK,idateo,cletiket,ip1, &
638 ip2,ip3,cltypvar,clnomvar)
639
640 if (ierr < 0 .and. ip2 < 10 .and. all(CrossCornsVarKindCH(:) == '')) then
641 write(*,*) 'bcsc_readcorns2: RSTDDEV ',ip2,jnum,clnomvar
642 call utl_abort('bcsc_readcorns2: Problem with constituents ' // &
643 'background stat file')
644 else if (ierr < 0 .and. ip2 == 0 .and. any(CrossCornsVarKindCH(:) /= '')) then
645 write(*,*) 'bcsc_readcorns2: Assumes content from cross-corrns ' // &
646 'input for ',clnomvar
647 clnomvarCrosscorns(varIndex)=clnomvar
648 exit
649 end if
650 if (ini /= jnum) then
651 call utl_abort('bcsc_readcorns2: Constituents ' // &
652 'background stat levels inconsistencies')
653 end if
654
655 ! Looking for FST record parameters..
656
657 if (ierr >= 0) then
658 idateo = -1
659 cletiket = 'CORRNS'
660 ip1 = -1
661 ip2 = jn
662 ip3 = -1
663 cltypvar = 'X'
664 ierr = utl_fstlir(ZCORNSSRC,nulbgst,INI,INJ,INK,idateo,cletiket, &
665 ip1,ip2,ip3,cltypvar,clnomvar)
666
667 if (ierr < 0) then
668 write(*,*) 'bcsc_readcorns2: CORRNS ',ip2,jnum,clnomvar
669 call utl_abort('bcsc_readcorns2: Problem with constituents ' // &
670 'background stat file')
671 end if
672 if (ini /= jnum .and. inj /= jnum) then
673 call utl_abort('bcsc_readcorns2: Constituents BG stat levels ' // &
674 'inconsistencies')
675 end if
676 else
677 write(*,*) 'WARNING from bcsc_readcorns2: Setting RSDTDDEV to ' // &
678 '1.D-15 for NOMVAR and JN: ',clnomvar,' ',jn
679 zstdsrc(1:jnum) = 1.D-15
680 zcornssrc(1:jnum,1:jnum) = 0.0D0
681 do jrow = 1, jnum
682 zcornssrc(jrow,jrow) = 1.0D0
683 end do
684 end if
685
686 rstddev(jstart:jstart+jnum-1,jn) = zstdsrc(1:jnum)
687 bgStats%corns(jstart:jstart+jnum-1,jstart:jstart+jnum-1,jn)= &
688 zcornssrc(1:jnum,1:jnum)
689
690 end do
691
692 deallocate(zcornssrc)
693 deallocate(zstdsrc)
694
695 end do
696
697 ! Read cross-correlation matrices
698
699 if (any(CrossCornsVarKindCH(:) /= '')) then
700
701 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
702
703 if (all(CrossCornsVarKindCH(:) /= bgStats%varNameList(varIndex))) cycle
704
705 clnomvar = bgStats%varNameList(varIndex)
706 jstart = bgStats%nsposit(varIndex)
707 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
708
709 if (clnomvarCrosscorns(varIndex) == clnomvar) then
710 clnomvarCrosscorns(varIndex)=''
711 end if
712
713 do varIndex2 = 1, bgStats%numvar3d+bgStats%numvar2d
714
715 if (varIndex == varIndex2) cycle
716
717 cletiket='CORRNS '//bgStats%varNameList(varIndex2)
718 ierr = fstinf(nulbgst,INI,INJ,INK,-1,cletiket,-1,-1,-1,'X',clnomvar)
719 if (ierr < 0 ) cycle
720
721 jstart2 = bgStats%nsposit(varIndex2)
722 jnum2 = bgStats%nsposit(varIndex2+1)-bgStats%nsposit(varIndex2)
723 allocate(zcornssrc(jnum,jnum2))
724
725 if (clnomvarCrosscorns(varIndex2) == bgStats%varNameList(varIndex2)) then
726 clnomvarCrosscorns(varIndex2)=''
727 end if
728
729 do jn = 0, bgStats%ntrunc
730
731 ierr = utl_fstlir(ZCORNSSRC,nulbgst,INI,INJ,INK,-1,cletiket,-1, &
732 jn,-1,'X',clnomvar)
733
734 if (ierr < 0) then
735 if (jn < 10) then
736 write(*,*) 'bcsc_readcorns2: CORRNS ',ip2,jnum,clnomvar, &
737 bgStats%varNameList(varIndex2)
738 call utl_abort('bcsc_readcorns2: Problem with constituents ' // &
739 'background stat file')
740 else
741 exit
742 end if
743 end if
744 if (ini /= jnum .and. inj /= jnum2) then
745 call utl_abort('bcsc_readcorns2: Constituents BG2 stat ' // &
746 'levels inconsistencies')
747 end if
748
749 bgStats%corns(jstart:jstart+jnum-1,jstart2:jstart2+jnum2-1,jn)= &
750 zcornssrc(1:jnum,1:jnum2)
751 bgStats%corns(jstart2:jstart2+jnum2-1,jstart:jstart+jnum-1,jn)= &
752 transpose(zcornssrc(1:jnum,1:jnum2))
753
754 end do
755 deallocate(zcornssrc)
756 end do
757 end do
758 end if
759
760 if (any(CrossCornsVarKindCH(:) /= '')) then
761 if (any(clnomvarCrosscorns(:) /= '')) then
762 write(*,*) 'bcsc_readcorns2: Missing matrix for ', &
763 clnomvarCrosscorns(1:bgStats%numvar3d+bgStats%numvar2d)
764 call utl_abort('bcsc_readcorns2: Missing correlations matrix')
765 end if
766 deallocate(clnomvarCrosscorns)
767 end if
768
769 ! Apply convolution to RSTDDEV correlation
770
771 call bcsc_convol
772
773 do jn = 0, bgStats%ntrunc
774
775 ! Re-build correlation matrix: factorization of corns with convoluted RSTDDEV
776 do jcol = 1, bgStats%nkgdim
777 bgStats%corns(1:bgStats%nkgdim,jcol,jn) = &
778 rstddev(1:bgStats%nkgdim,jn) * &
779 bgStats%corns(1:bgStats%nkgdim,jcol,jn) * rstddev(jcol,jn)
780 end do
781
782 end do
783
784 end subroutine bcsc_readCorns2
785
786 !--------------------------------------------------------------------------
787 ! bcsc_convol
788 !--------------------------------------------------------------------------
789 subroutine bcsc_convol
790 implicit none
791
792 ! Locals:
793 real(8) :: dlfact2,dlc,dsummed
794 real(8) :: dtlen,zr,dlfact
795 integer :: jn,latIndex,jk,varIndex,levelIndex,nsize,ierr
796 real(8) :: zleg(0:bgStats%ntrunc,bgStats%nj)
797 real(8) :: zsp(0:bgStats%ntrunc,bgStats%nkgdim)
798 real(8) :: zgr(bgStats%nj,bgStats%nkgdim)
799 real(8) :: zrmu(bgStats%nj)
800 integer :: nlev_MT,ini,inj,ink,nulcorns
801 real(8), allocatable :: wtemp(:,:,:)
802 real(8), allocatable :: hcorrel(:,:,:),hdist(:)
803 logical :: lfound
804 integer :: fnom, fstouv, fstfrm, fclos
805
806 do latIndex = 1, bgStats%nj
807 zrmu(latIndex) = gst_getrmu(latIndex,gstID)
808 end do
809
810 ! CONVERT THE CORRELATIONS IN SPECTRAL SPACE INTO SPECTRAL
811 ! COEFFICIENTS OF THE CORRELATION FUNCTION AND FUNCTION TO BE
812 ! SELF-CONVOLVED
813
814 do jn = 0, bgStats%ntrunc
815 dlfact = ((2.0d0*jn+1.0d0)/2.0d0)**(0.25d0)
816 dlfact2 = ((2.0d0*jn +1.0d0)/2.0d0)**(0.25d0)
817 do jk = 1, bgStats%nkgdim
818 zsp(jn,jk) = rstddev(jk,jn)*dlfact*dlfact2
819 end do
820 end do
821
822 ! Transform to physical space
823 call gst_zleginv(gstID,zgr,zsp,bgStats%nkgdim)
824
825 ! Truncate in horizontal extent with Gaussian window
826
827 varIndex=1
828 do jk = 1, bgStats%nkgdim
829 if (jk == bgStats%nsposit(varIndex)) then
830 dtlen = rpor(varIndex)
831 varIndex=varIndex+1
832 end if
833 if (dtlen > 0.0d0) then
834 dlc = 1.d0/dble(dtlen)
835 dlc = 0.5d0*dlc*dlc
836 do latIndex = 1, bgStats%nj
837 zr = ec_ra * acos(zrmu(latIndex))
838 dlfact = dexp(-(zr**2)*dlc)
839 zgr(latIndex,jk) = dlfact*zgr(latIndex,jk)
840 end do
841 end if
842
843 !write(*,*) 'zeroing length (km)=',jk,dtlen/1000.0
844 end do
845
846 ! Transform back to spectral space
847 call gst_zlegdir(gstID,zgr,zsp,bgStats%nkgdim)
848
849 ! Convert back to correlations
850 do jk = 1, bgStats%nkgdim
851 do jn = 0, bgStats%ntrunc
852 zsp(jn,jk) = zsp(jn,jk)*(2.0d0/(2.0d0*jn+1.0d0))**(0.25d0)
853 end do
854 end do
855
856 ! PUT BACK INTO RSTDDEV
857 do jn = 0, bgStats%ntrunc
858 do jk = 1, bgStats%nkgdim
859 rstddev(jk,jn) = zsp(jn,jk)
860 end do
861 end do
862
863 ! Re-normalize to ensure correlations
864 do jk = 1, bgStats%nkgdim
865 dsummed = 0.d0
866 do jn = 0, bgStats%ntrunc
867 dsummed = dsummed+ dble(rstddev(jk,jn)**2)*sqrt(((2.d0*jn)+1.d0)/2.d0)
868 end do
869 dsummed = sqrt(dsummed)
870 do jn = 0, bgStats%ntrunc
871 if(dsummed > 1.d-30) rstddev(jk,jn) = rstddev(jk,jn)/dsummed
872 end do
873 end do
874
875 ! CONVERT THE SPECTRAL COEFFICIENTS OF THE CORRELATION FUNCTION
876 ! BACK INTO CORRELATIONS OF SPECTRAL COMPONENTS
877 do jn = 0, bgStats%ntrunc
878 dlfact = sqrt(0.5d0)*(1.0d0/((2.0d0*jn+1)/2.0d0))**0.25d0
879 do jk = 1, bgStats%nkgdim
880 rstddev(jk,jn) = rstddev(jk,jn)*dlfact
881 end do
882 end do
883
884 if ( .not.getPhysSpaceHCorrel .or. .not.getPhysSpaceStats ) return
885
886 ! Compute resultant physical space horizontal correlations and
887 ! 1/e correlation length from correlation array if not available
888
889 if (allocated(hcorrel)) deallocate(hcorrel)
890 allocate(hcorrel(bgStats%nj,bgStats%nlev, &
891 bgStats%numvar3d+bgStats%numvar2d))
892 if (allocated(wtemp)) deallocate(wtemp)
893 allocate(wtemp(0:bgStats%ntrunc,bgStats%nj,1))
894 if (allocated(bgStats%hcorrlen)) deallocate(bgStats%hcorrlen)
895 allocate(bgStats%hcorrlen(bgStats%nlev, &
896 bgStats%numvar3d+bgStats%numvar2d))
897 if (allocated(hdist)) deallocate(hdist)
898 allocate(hdist(bgStats%nj))
899
900 lfound=.true.
901 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
902 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
903 do levelIndex = 1, nlev_MT
904 ierr = utl_fstlir(hcorrel(:,levelIndex,varIndex),nulbgst,INI,INJ, &
905 INK,-1,'HCORREL',levelIndex,-1,-1,'X', &
906 bgStats%varNameList(varIndex))
907 if (ierr < 0 ) then
908 lfound=.false.
909 exit
910 end if
911 end do
912 ierr = utl_fstlir(bgStats%hcorrlen(1:nlev_MT,varIndex),nulbgst,INI, &
913 INJ,INK,-1,'HCORRLEN',-1,-1,-1,'X',bgStats%varNameList(varIndex))
914 if (ierr < 0 ) then
915 lfound=.false.
916 exit
917 end if
918 end do
919
920 if (lfound) return
921
922 do latIndex = 1, bgStats%nj
923 hdist(latIndex)=ec_ra*acos(zrmu(latIndex))
924 end do
925
926 zleg(:,:)=0.0d0
927 wtemp(:,:,:)=0.0d0
928 hcorrel(:,:,:)=0.0d0
929 bgStats%hcorrlen(:,:)=0.0
930
931 do latIndex = mmpi_myid+1, bgStats%nJ, mmpi_nprocs
932 do jn = 0, bgStats%ntrunc
933 wtemp(jn,latIndex,1) = gst_getzleg(jn,latIndex,gstID)
934 end do
935 end do
936
937 nsize=bgStats%nJ*(bgStats%ntrunc+1)
938 call rpn_comm_allreduce(wtemp(0:bgStats%ntrunc,1:bgStats%nJ,1), &
939 zleg(0:bgStats%ntrunc,1:bgStats%nJ),nsize,"mpi_double_precision", &
940 "mpi_sum","GRID",ierr)
941
942 deallocate(wtemp)
943 allocate(wtemp(bgStats%nj, bgStats%nlev, &
944 bgStats%numvar3d+bgStats%numvar2d))
945 wtemp(:,:,:)=0.0
946
947 varIndex = 1
948 levelIndex = 1
949 do jk = 1, bgStats%nkgdim
950 if (jk == bgStats%nsposit(varIndex+1)) then
951 varIndex = varIndex+1
952 levelIndex = 1
953 end if
954
955 do latIndex = mmpi_myid+1, bgStats%nj, mmpi_nprocs
956 do jn = 0, bgStats%ntrunc
957 wtemp(latIndex,levelIndex,varIndex) = wtemp(latIndex,levelIndex, &
958 varIndex)+rstddev(jk,jn)*rstddev(jk,jn)* &
959 sqrt(2.0)*sqrt(2.0*jn+1.0)*zleg(jn,latIndex)
960 end do
961 end do
962 levelIndex = levelIndex+1
963 end do
964
965 nsize=bgStats%nj*bgStats%nkgdim
966 call rpn_comm_allreduce(wtemp,hcorrel,nsize,"mpi_double_precision", &
967 "mpi_sum","GRID",ierr)
968 deallocate(wtemp)
969
970 varIndex = 1
971 levelIndex = 1
972 do jk = 1, bgStats%nkgdim
973 if (jk == bgStats%nsposit(varIndex+1)) then
974 varIndex = varIndex+1
975 levelIndex = 1
976 end if
977 do latIndex=bgStats%nj-1,2,-1
978 if (hcorrel(latIndex,levelIndex,varIndex) <= 0.368) then ! 1/e ~ 0.368
979 bgStats%hcorrlen(levelIndex,varIndex) = (hdist(latIndex)* &
980 (hcorrel(latIndex+1,levelIndex,varIndex)-0.368) &
981 + hdist(latIndex+1)*(0.368-hcorrel(latIndex, &
982 levelIndex,varIndex))) &
983 /(hcorrel(latIndex+1,levelIndex,varIndex)- &
984 hcorrel(latIndex,levelIndex,varIndex))
985 exit
986 end if
987 end do
988 levelIndex = levelIndex+1
989 end do
990
991
992 if ( mmpi_myid == 0 ) then
993
994 if (WritePhysSpaceStats) then
995 nulcorns = 0
996 ierr = fnom(nulcorns,'bCovarSetupChem_out.fst','STD+RND',0)
997 ierr = fstouv(nulcorns,'RND')
998
999 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
1000 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1001 do levelIndex = 1, nlev_MT
1002 ierr = utl_fstecr(hcorrel(:,levelIndex,varIndex),-32,nulcorns,0,0,0, &
1003 1,bgStats%nj,1,levelIndex,0,nlev_MT,'X', &
1004 bgStats%varNameList(varIndex), &
1005 'HCORREL','X',0,0,0,0,5,.true.)
1006 end do
1007 ierr = utl_fstecr(bgStats%hcorrlen(1:nlev_MT,varIndex),-32,nulcorns, &
1008 0,0,0,1,1,nlev_MT,0,0,0,'X',bgStats%varNameList(varIndex), &
1009 'HCORRLEN','X',0,0,0,0,5,.true.)
1010 ierr = utl_fstecr(hdist(1:bgStats%nj),-32,nulcorns,0,0,0,1, &
1011 bgStats%nj,1,0,0,0,'X',bgStats%varNameList(varIndex), &
1012 'HDIST','X',0,0,0,0,5,.true.)
1013 end do
1014
1015 ierr = fstfrm(nulcorns)
1016 ierr = fclos(nulcorns)
1017
1018 end if
1019
1020 write(*,*)
1021 write(*,*) 'bcsc_convol: Horizontal correlations'
1022 write(*,*)
1023 write(*,*) 'Separation distances (km)'
1024 write(*,*) bgStats%nj-bgStats%nj*4/5+1
1025 write(*,*) hdist(bgStats%nj*4/5:bgStats%nj)/1000.00
1026 write(*,*)
1027 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
1028 if (varIndex <= bgStats%numvar3d) then
1029 write(*,*) bgStats%varNameList(varIndex), bgStats%nlev
1030 do levelIndex = 1, bgStats%nlev
1031 write(*,'(i3,f10.2,3000f6.2)') levelIndex, &
1032 bgStats%hcorrlen(levelIndex,varIndex)/1000.00, &
1033 hcorrel(bgStats%nj*4/5:bgStats%nj,levelIndex,varIndex)
1034 end do
1035 else
1036 write(*,*) bgStats%varNameList(varIndex), 1
1037 write(*,'(i3,f10.2,3000f6.2)') 1, &
1038 bgStats%hcorrlen(1,varIndex)/1000.00, &
1039 hcorrel(bgStats%nj*4/5:bgStats%nj,1,varIndex)
1040 end if
1041 write(*,*)
1042 end do
1043 end if
1044
1045 if (allocated(hcorrel)) deallocate(hcorrel)
1046 if (allocated(hdist)) deallocate(hdist)
1047
1048 end subroutine bcsc_convol
1049
1050 !--------------------------------------------------------------------------
1051 ! bcsc_rdstddev
1052 !--------------------------------------------------------------------------
1053 subroutine bcsc_rdstddev
1054 !
1055 !:Purpose: To read stddev and to set as 3D fields.
1056 !
1057 implicit none
1058
1059 ! Locals:
1060 integer :: ikey
1061 real(8) :: rcoord(10000)
1062 ! standard file variables
1063 integer :: ini,inj,ink
1064 integer :: ip1,ip2,ip3
1065 integer :: idate(100)
1066
1067 ! Reading the data
1068
1069 idate(1) = -1
1070 ip2 = -1
1071 ip3 = -1
1072
1073 ! Get latitudes and longitudes if available
1074
1075 ip1=-1
1076 ikey = utl_fstlir(rcoord,nulbgst,ini,inj,ink, &
1077 idate(1),' ',ip1,ip2,ip3,' ','^^')
1078
1079 if (ikey >= 0) then
1080 if (allocated(rlatr)) deallocate(rlatr)
1081 allocate(rlatr(inj))
1082 rlatr(1:inj) = rcoord(1:inj)
1083 rlatr(1:inj) = rlatr(1:inj)*MPC_RADIANS_PER_DEGREE_R8
1084 else
1085 ! Assume same as bgStats%lat
1086 if (allocated(rlatr)) deallocate(rlatr)
1087 allocate(rlatr(bgStats%nj))
1088 inj = bgStats%nj
1089 rlatr(1:inj) = bgStats%lat(1:inj)
1090 end if
1091
1092 ikey = utl_fstlir(rcoord,nulbgst,ini,inj,ink, &
1093 idate(1),' ',ip1,ip2,ip3,' ','>>')
1094
1095 if (ikey >= 0) then
1096 if (allocated(rlongr)) deallocate(rlongr)
1097 allocate(rlongr(ini+1))
1098 rlongr(1:ini) = rcoord(1:ini)
1099 rlongr(1:ini) = rlongr(1:ini)*MPC_RADIANS_PER_DEGREE_R8
1100 else if (stddevMode /= 'SP2D') then
1101 ! Assume same as bgStats%lon
1102 if (allocated(rlongr)) deallocate(rlongr)
1103 allocate(rlongr(bgStats%ni+1))
1104 ini = bgStats%ni
1105 rlongr(1:ini) = bgStats%lon(1:ini)
1106 end if
1107 rlongr(ini+1) = 360.*MPC_RADIANS_PER_DEGREE_R8
1108
1109 ! Read specified input type for error std. dev.
1110
1111 if(stddevMode == 'GD3D') then
1112 call bcsc_rdstd3D
1113 elseif(stddevMode == 'GD2D') then
1114 call bcsc_rdstd
1115 elseif(stddevMode == 'SP2D') then
1116 call bcsc_rdspstd
1117 else
1118 call utl_abort('bcsc_rdstddev: unknown stddevMode')
1119 end if
1120
1121 end subroutine bcsc_rdstddev
1122
1123 !--------------------------------------------------------------------------
1124 ! bcsc_rdspstd
1125 !--------------------------------------------------------------------------
1126 subroutine bcsc_rdspstd
1127 implicit none
1128
1129 ! Locals:
1130 integer :: varIndex,jn,inix,injx,inkx
1131 integer :: ikey, levelIndexo, firstn,lastn
1132 real(8) :: zsp(0:bgStats%ntrunc,max(nlev_M,nlev_T))
1133 real(8) :: zspbuf(max(nlev_M,nlev_T))
1134 real(8) :: zgr(bgStats%nj,max(nlev_M,nlev_T))
1135 ! standard file variables
1136 integer :: ini,inj,ink
1137 integer :: ip1,ip2,ip3
1138 integer :: idate(100),nlev_MT
1139 character(len=2) :: cltypvar
1140 character(len=4) :: clnomvar
1141 character(len=12) :: cletiket
1142 integer :: fstinf
1143
1144 stddev(:,:,:) = 0.0d0
1145
1146 ! Reading the Legendre poly representation of the 2D background
1147 ! error std. dev. field
1148
1149 idate(1) = -1
1150 ip1 = -1
1151 ip2 = -1
1152 ip3 = -1
1153
1154 cletiket = 'SPSTDDEV'
1155 cltypvar = 'X'
1156
1157 do varIndex = 1,bgStats%numvar3d+bgStats%numvar2d
1158 clnomvar = bgStats%varNameList(varIndex)
1159 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1160
1161 firstn = -1
1162 do jn = 0, bgStats%ntrunc
1163 ip2 = jn
1164 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3, &
1165 cltypvar,clnomvar)
1166
1167 if(ikey >= 0 ) then
1168 ikey = utl_fstlir(zspbuf(1:nlev_MT),nulbgst,ini,inj,ink,idate(1), &
1169 cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1170 else
1171 if(firstn == -1) firstn = jn
1172 lastn = jn
1173 zspbuf(:) = 0.0d0
1174 end if
1175
1176 if (ini /= nlev_MT) then
1177 if (ini == nlev_MT-1) then
1178 zspbuf(nlev_MT) = zspbuf(ini)
1179 write(*,*) 'WARNING in CHM_RDSPSTD: ini and nlev_MT not ', &
1180 'same size - ',ini,nlev_MT
1181 else
1182 write(*,*) 'JN, INI, nlev_MT, NOMVAR: ',jn,ini,nlev_MT,' ',clnomvar
1183 call utl_abort('CHM_RDSPSTD: Constituents background stats ' // &
1184 'levels inconsistency')
1185 end if
1186 end if
1187
1188 do levelIndexo = 1, nlev_MT
1189 zsp(jn,levelIndexo) = zspbuf(levelIndexo)
1190 end do
1191 end do
1192
1193 if (mmpi_myid == 0 .and. firstn /= -1) then
1194 write(*,*) 'WARNING: CANNOT FIND SPSTD FOR ',clnomvar, &
1195 ' AT N BETWEEN ',firstn,' AND ',lastn,', SETTING TO ZERO!!!'
1196 end if
1197
1198 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
1199
1200 do jn = 1, bgStats%ni+1
1201 stddev(jn,1:bgStats%nj, &
1202 bgStats%nsposit(varIndex):bgStats%nsposit(varIndex+1)-1) = &
1203 zgr(1:bgStats%nj,1:nlev_MT)
1204 end do
1205
1206 end do
1207
1208 end subroutine bcsc_rdspstd
1209
1210 !--------------------------------------------------------------------------
1211 ! bcsc_rdspstd_newfmt
1212 !--------------------------------------------------------------------------
1213 subroutine bcsc_rdspstd_newfmt
1214 implicit none
1215
1216 ! Locals:
1217 integer :: varIndex,jn,inix,injx,inkx,ntrunc_file
1218 integer :: ikey,levelIndexo
1219 real(8) :: zsp(0:bgStats%ntrunc,max(nlev_M,nlev_T))
1220 real(8), allocatable :: zspbuf(:)
1221 real(8) :: zgr(bgStats%nj,max(nlev_M,nlev_T))
1222 ! standard file variables
1223 integer :: ini,inj,ink
1224 integer :: ip1,ip2,ip3
1225 integer :: idate(100),nlev_MT
1226 character(len=2) :: cltypvar
1227 character(len=4) :: clnomvar
1228 character(len=12) :: cletiket
1229 integer :: fstinf
1230
1231 stddev(:,:,:) = 0.0d0
1232
1233 ! Reading the data
1234
1235 idate(1) = -1
1236 ip2 = -1
1237 ip3 = -1
1238
1239 cletiket = 'SPSTDDEV'
1240 cltypvar = 'X'
1241
1242 ! Check if file is old format
1243
1244 ip1 = -1
1245 clnomvar = bgStats%varNameList(1)
1246 ikey = fstinf(nulbgst,inix,injx,inkx,idate(1),cletiket,ip1,ip2,ip3, &
1247 cltypvar,clnomvar)
1248 write(*,*) 'ini,inj,ink=',inix,injx,inkx
1249 if(inix > 1) then
1250 write(*,*) 'bcsc_rdspstd_newfmt: ini>1, SPSTDDEV is in old format, ', &
1251 ' calling bcsc_RDSPSTD...'
1252 call bcsc_rdspstd
1253 return
1254 end if
1255
1256 ! write(*,*) 'Reading for 3D and 2D variables'
1257
1258 do varIndex = 1,bgStats%numvar3d+bgStats%numvar2d
1259 clnomvar = bgStats%varNameList(varIndex)
1260 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1261
1262 !write(*,*) 'Reading ',clnomvar
1263
1264 do levelIndexo = 1, nlev_MT
1265 if (nlev_MT == 1) then
1266 ip1 = -1
1267 else if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
1268 ip1 = vco_anl%ip1_M(levelIndexo)
1269 else
1270 ip1 = vco_anl%ip1_T(levelIndexo)
1271 end if
1272
1273 ikey = fstinf(nulbgst,inix,ntrunc_file,inkx,idate(1),cletiket, &
1274 ip1,ip2,ip3,cltypvar,clnomvar)
1275 ntrunc_file = ntrunc_file-1
1276
1277 allocate(zspbuf(0:ntrunc_file))
1278
1279 if(ikey >= 0 ) then
1280 ikey = utl_fstlir(zspbuf(0:ntrunc_file),nulbgst,ini,inj,ink, &
1281 idate(1),cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
1282 else
1283 write(*,*) 'bcsc_rdspstd_newfmt: ',varIndex,clnomvar,nlev_MT, &
1284 levelIndexo,ikey,bgStats%ntrunc,ntrunc_file
1285 call utl_abort('bcsc_rdspstd_newfmt: SPSTDDEV record not found')
1286 end if
1287
1288 zsp(:,levelIndexo) = 0.0d0
1289 do jn = 0, min(bgStats%ntrunc,ntrunc_file)
1290 zsp(jn,levelIndexo) = zspbuf(jn)
1291 end do
1292 deallocate(zspbuf)
1293 end do
1294
1295 call gst_zleginv(gstID,zgr(:,1:nlev_MT),zsp(:,1:nlev_MT),nlev_MT)
1296
1297 do jn = 1, bgStats%ni+1
1298 stddev(jn,1:bgStats%nj,bgStats%nsposit(varIndex): &
1299 bgStats%nsposit(varIndex+1)-1) = zgr(1:bgStats%nj,1:nlev_MT)
1300 end do
1301
1302 end do
1303
1304 end subroutine bcsc_rdspstd_newfmt
1305
1306 !--------------------------------------------------------------------------
1307 ! bcsc_rdstd
1308 !--------------------------------------------------------------------------
1309 subroutine bcsc_rdstd
1310 !
1311 !:Purpose: To read 2D stddev and to store as 3D
1312 !
1313 implicit none
1314
1315 ! Locals:
1316 integer :: varIndex,in
1317 integer :: ikey
1318 real(8), allocatable :: stddev3d(:,:,:)
1319 ! standard file variables
1320 integer :: ini,inj,ink
1321 integer :: ip1,ip2,ip3
1322 integer :: idate(100),nlev_MT
1323 character(len=2) :: cltypvar
1324 character(len=4) :: clnomvar
1325 character(len=12) :: cletiket
1326 integer :: fstinf
1327 real(8), allocatable :: vlev(:),vlevout(:)
1328
1329 stddev(:,:,:) = 0.0d0
1330
1331 ! Reading the data
1332
1333 idate(1) = -1
1334 ip1 = -1
1335 ip2 = -1
1336 ip3 = -1
1337
1338 cletiket = 'STDDEV'
1339 cltypvar=' '
1340
1341 ! Reading for 3D and 2D variables
1342
1343 do varIndex = 1,bgStats%numvar3d+bgStats%numvar2d
1344 clnomvar = bgStats%varNameList(varIndex)
1345 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1346
1347 ikey = fstinf(nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3, &
1348 cltypvar,clnomvar)
1349 if (ikey < 0 .or. ini > 1 .or. ink /= nlev_MT) then
1350 write(*,*) 'bcsc_rdstd: ',varIndex,clnomvar,ikey,ini,ink,nlev_MT
1351 call utl_abort(': bcsc_rdstd record not found or incorrect')
1352 end if
1353
1354 allocate(stddev3d(1,inj,ink))
1355 stddev3d(1,:,:) = 0.0D0
1356 allocate(vlev(ink),vlevout(nlev_MT))
1357 vlev(:) = 1.0D0
1358 vlevout(:) = 1.0D0
1359
1360 ikey = utl_fstlir(stddev3d(1,:,:), nulbgst, ini, inj, ink, &
1361 idate(1), cletiket, ip1, ip2, ip3, cltypvar, clnomvar)
1362
1363 if (ikey < 0) then
1364 write(*,*) 'bcsc_rdstd: ',varIndex,clnomvar,nlev_MT,ikey
1365 call utl_abort(': bcsc_rdstd record not found')
1366 end if
1367
1368 ! Extend to 3D
1369 if (inj == bgStats%nj) then
1370 do in = 1, bgStats%ni+1
1371 stddev(in,:,bgStats%nsposit(varIndex): &
1372 bgStats%nsposit(varIndex+1)-1) = stddev3d(1,:,:)
1373 end do
1374 else
1375 ! Interpolate in lat
1376 call gsv_field3d_hbilin(stddev3d, 1, inj, ink, rlongr, rlatr, vlev, &
1377 stddev(:,:,bgStats%nsposit(varIndex): &
1378 bgStats%nsposit(varIndex+1)-1), bgStats%ni+1, bgStats%ni, &
1379 nlev_MT, bgStats%lon, bgStats%lat, vlevout)
1380 end if
1381 deallocate(stddev3d, vlev, vlevout)
1382
1383 end do
1384
1385 end subroutine bcsc_rdstd
1386
1387 !--------------------------------------------------------------------------
1388 ! bcsc_rdstd3d
1389 !--------------------------------------------------------------------------
1390 subroutine bcsc_rdstd3d
1391 !
1392 !:Purpose: To read 3D stddev.
1393 !
1394 ! Originally based on bcsc_rdspstd_newfmt
1395 !
1396 implicit none
1397
1398 ! Locals:
1399 integer :: varIndex
1400 integer :: ikey,levelIndexo
1401 real(8), allocatable :: stddev3d(:,:,:)
1402 ! standard file variables
1403 integer :: ini,inj,ink
1404 integer :: ip1,ip2,ip3
1405 integer :: idate(100),nlev_MT
1406 character(len=2) :: cltypvar
1407 character(len=4) :: clnomvar
1408 character(len=12) :: cletiket
1409 integer :: fstinf
1410 real(8) :: vlev(1),vlevout(1)
1411
1412 stddev(:,:,:) = 0.0d0
1413 vlev(:) = 1.0D0
1414 vlevout(:) = 1.0D0
1415
1416 ! Reading the data
1417
1418 idate(1) = -1
1419 ip2 = -1
1420 ip3 = -1
1421
1422 cletiket = 'STDDEV3D'
1423 cltypvar=' '
1424
1425 ! Reading for 3D and 2D variables
1426
1427 do varIndex = 1,bgStats%numvar3d+bgStats%numvar2d
1428 clnomvar = bgStats%varNameList(varIndex)
1429 nlev_MT = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1430
1431 !write(*,*) 'Reading ',clnomvar
1432
1433 do levelIndexo = 1, nlev_MT
1434 if (nlev_MT == 1) then
1435 ip1 = -1
1436 else if(vnl_varLevelFromVarName(clnomvar) == 'MM') then
1437 ip1 = vco_anl%ip1_M(levelIndexo)
1438 else
1439 ip1 = vco_anl%ip1_T(levelIndexo)
1440 end if
1441
1442 ikey = fstinf(nulbgst,ini,inj,ink,idate(1),cletiket,ip1,ip2,ip3, &
1443 cltypvar,clnomvar)
1444 if (ikey < 0) then
1445 write(*,*) 'bcsc_rdstd3d: ',varIndex,clnomvar,ikey,levelIndexo
1446 call utl_abort(': bcsc_RDSTD record not foun0d')
1447 end if
1448
1449 allocate(stddev3d(ini+1, inj, 1))
1450 stddev3d(:,:,1) = 0.0D0
1451
1452 ikey = utl_fstlir(stddev3d(1:ini,:,1), nulbgst, ini, inj, ink, &
1453 idate(1), cletiket, ip1, ip2, ip3, cltypvar, clnomvar)
1454 if (ikey < 0) then
1455 write(*,*) 'bcsc_rdstd3d: ',varIndex,clnomvar,nlev_MT,levelIndexo, &
1456 ikey,ip1
1457 call utl_abort(': bcsc_RDSTD3D record not found')
1458 end if
1459
1460 ! Move to stddev
1461 if (inj == bgStats%nj .and. ini == bgStats%ni) then
1462 stddev(1:bgStats%ni,:,bgStats%nsposit(varIndex)+(levelIndexo-1)) = &
1463 stddev3d(1:bgStats%ni,:,1)
1464 stddev(bgStats%ni+1,:,bgStats%nsposit(varIndex)+(levelIndexo-1)) = &
1465 stddev3d(1,:,1)
1466 else
1467 ! Interpolate in lat and long
1468 stddev3d(ini+1,:,1) = stddev3d(1,:,1)
1469 call gsv_field3d_hbilin(stddev3d, ini+1, inj, 1, rlongr, rlatr, vlev, &
1470 stddev(:,:,bgStats%nsposit(varIndex)+(levelIndexo-1)), &
1471 bgStats%ni+1, bgStats%nj, 1, &
1472 bgStats%lon, bgStats%lat, vlevout)
1473 end if
1474 deallocate(stddev3d)
1475
1476 end do
1477 end do
1478
1479 end subroutine bcsc_rdstd3d
1480
1481 !--------------------------------------------------------------------------
1482 ! bcsc_sucorns2
1483 !--------------------------------------------------------------------------
1484 subroutine bcsc_sucorns2
1485 implicit none
1486
1487 ! Locals:
1488 real(8) :: eigenval(bgStats%nkgdim)
1489 real(8) :: eigenvalsqrt(bgStats%nkgdim)
1490 real(8), allocatable :: eigenvec(:,:),result(:,:)
1491 integer :: jn,jk1,jk2,varIndex,ierr
1492 integer :: ilwork,info,jnum,jstart,nsize
1493 real(8) :: zwork(2*4*bgStats%nkgdim)
1494 real(8) :: ztlen,zcorr,zr,zpres1,zpres2,eigenvalmax
1495 real(8), allocatable :: corns_temp(:,:,:)
1496 logical, allocatable :: lfound_sqrt(:)
1497
1498 ! Apply vertical localization to correlations of 3D fields.
1499 ! Currently assumes no-cross correlations for variables (block diagonal matrix)
1500
1501 if ( bgStats%numvar3d > 0 ) then
1502 do varIndex = 1, bgStats%numvar3d
1503 ztlen = rvloc(varIndex) ! specify length scale (in units of ln(Pressure))
1504
1505 jstart = bgStats%nsposit(varIndex)
1506 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1507
1508 if(ztlen > 0.0d0) then
1509 ! calculate 5'th order function (from Gaspari and Cohn)
1510 do jk1 = 1, jnum
1511 zpres1 = log(bgStats%vlev(jk1))
1512 do jk2 = 1, jnum
1513 zpres2 = log(bgStats%vlev(jk2))
1514 zr = abs(zpres2 - zpres1)
1515 zcorr = gasparicohn(ztlen,zr)
1516 bgStats%corns(jstart-1+jk1,jstart-1+jk2,0:bgStats%ntrunc) = &
1517 bgStats%corns(jstart-1+jk1,jstart-1+jk2, &
1518 0:bgStats%ntrunc)*zcorr
1519 end do
1520 end do
1521 end if
1522 end do
1523 end if
1524
1525 ! Compute total vertical correlations and its inverse
1526 ! (currently for each block matrix).
1527
1528 if (getPhysSpaceStats .or. trim(bcsc_mode) == 'BackgroundCheck') then
1529 call bcsc_corvertSetup
1530 end if
1531
1532 if (trim(bcsc_mode) == 'BackgroundCheck') return
1533
1534 allocate(lfound_sqrt(bgStats%numvar3d+bgStats%numvar2d))
1535 if (ReadWrite_sqrt) then
1536 ! if desired, read precomputed sqrt of corns
1537 call readcorns(lfound_sqrt,bgStats%ntrunc,'CORNS_SQRT')
1538 else
1539 lfound_sqrt(:)=.false.
1540 end if
1541
1542 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
1543
1544 if (any(CrossCornsVarKindCH(:) /= '')) then
1545 jstart=1
1546 jnum=bgStats%nkgdim
1547 else
1548 jstart = bgStats%nsposit(varIndex)
1549 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1550 end if
1551
1552 if (.not.lfound_sqrt(varIndex)) then
1553
1554 ! compute square-root of corns for each total wavenumber
1555
1556 allocate(corns_temp(jnum,jnum,0:bgStats%ntrunc))
1557 allocate(eigenvec(jnum,jnum),result(jnum,jnum))
1558 corns_temp(:,:,:) = 0.0d0
1559 do jn = mmpi_myid, bgStats%ntrunc, mmpi_nprocs
1560
1561 eigenvec(1:jnum,1:jnum) = bgStats%corns(jstart:jstart-1+jnum, &
1562 jstart:jstart-1+jnum,jn)
1563
1564 ! CALCULATE EIGENVALUES AND EIGENVECTORS.
1565 ilwork = 4*jnum*2
1566 call dsyev('V','U',jnum,eigenvec,jnum,eigenval,zwork,ilwork,info)
1567 if(info /= 0) then
1568 write(*,*) 'bcsc_sucorns2: non-zero value of info =',info, &
1569 ' returned by dsyev for wavenumber ',jn,varIndex
1570 call utl_abort('bcsc_sucorns2')
1571 end if
1572
1573 ! set selected number of eigenmodes to zero
1574 if(numModeZero > 0) then
1575 do jk1 = 1, numModeZero
1576 eigenval(jk1) = 0.0d0
1577 end do
1578 end if
1579
1580 eigenvalmax = maxval(eigenval(1:jnum))
1581 do jk1 = 1, jnum
1582 ! if(eigenval(jk1) < 1.0d-15) then
1583 if(eigenval(jk1) < 1.0d-8*eigenvalmax) then
1584 eigenvalsqrt(jk1) = 0.0d0
1585 else
1586 eigenvalsqrt(jk1) = sqrt(eigenval(jk1))
1587 end if
1588 end do
1589
1590 ! E * lambda^1/2
1591 do jk1 = 1, jnum
1592 result(1:jnum,jk1) = eigenvec(1:jnum,jk1)*eigenvalsqrt(jk1)
1593 end do
1594
1595 ! (E * lambda^1/2) * E^T
1596 do jk1 = 1, jnum
1597 do jk2 = 1, jnum
1598 corns_temp(1:jnum,jk1,jn) = corns_temp(1:jnum,jk1,jn) &
1599 + result(1:jnum,jk2)*eigenvec(jk1,jk2)
1600 end do
1601 end do
1602
1603 end do ! jn
1604
1605 nsize = jnum*jnum*(bgStats%ntrunc+1)
1606 call rpn_comm_allreduce(corns_temp, &
1607 bgStats%corns(jstart:jstart+jnum-1,jstart:jstart+jnum-1,:), &
1608 nsize,"mpi_double_precision","mpi_sum","GRID",ierr)
1609 deallocate(corns_temp,eigenvec,result)
1610
1611 end if
1612
1613 if (any(CrossCornsVarKindCH(:) /= '')) exit
1614
1615 end do
1616
1617 deallocate(lfound_sqrt)
1618 if(ReadWrite_sqrt) then
1619 ! Write computed sqrt to a separate file.
1620 call writecorns(bgStats%ntrunc,'CORNS_SQRT',-1)
1621 end if
1622
1623 end subroutine bcsc_sucorns2
1624
1625 !--------------------------------------------------------------------------
1626 ! bcsc_corvertSetup
1627 !--------------------------------------------------------------------------
1628 subroutine bcsc_corvertSetup
1629 !
1630 !:Purpose: To compute total vertical correlations (bcsc_corvert) and its
1631 ! inverse (bgStats%corverti; currently for each block matrix).
1632 !
1633 !
1634 !:Note: Currently assumes no (or neglects) cross-correlations
1635 !
1636 implicit none
1637
1638 ! Locals:
1639 real(8) :: eigenval(bgStats%nkgdim)
1640 real(8), allocatable :: eigenvec(:,:),result(:,:)
1641 integer :: jn,jk1,jk2,varIndex,numvartot,jstart
1642 integer :: ilwork,info,jnum,nvlev,ierr
1643 real(8) :: zwork(2*4*bgStats%nkgdim)
1644 real(8) :: eigenvalmax
1645 integer iulcorvert
1646 ! Standard file variables
1647 character(len=4) :: clnomvar
1648 character(len=12) :: cletiket
1649 integer :: fnom,fstouv,fstfrm,fclos
1650 logical, allocatable :: lfound(:)
1651
1652 ! Compute total vertical correlations and its inverse
1653 ! (currently for each block matrix).
1654
1655 if (mmpi_myid == 0 .and. WritePhysSpaceStats) then
1656 iulcorvert = 0
1657 ierr = fnom(iulcorvert,'bCovarSetupChem_out.fst','STD+RND',0)
1658 ierr = fstouv(iulcorvert,'RND')
1659 end if
1660
1661 nvlev=-1
1662 numvartot=bgStats%numvar3d
1663 do varIndex = 1, numvartot
1664
1665 jstart = bgStats%nsposit(varIndex)
1666 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1667
1668 bgStats%corvert(1:jnum,1:jnum,varIndex) = 0.0d0
1669 do jn = 0, bgStats%ntrunc
1670 bgStats%corvert(1:jnum,1:jnum,varIndex) = &
1671 bgStats%corvert(1:jnum,1:jnum,varIndex)+ ((2*jn+1)* &
1672 bgStats%corns(bgStats%nsposit(varIndex): &
1673 bgStats%nsposit(varIndex+1)-1,&
1674 bgStats%nsposit(varIndex):bgStats%nsposit(varIndex+1)-1,jn))
1675 end do
1676
1677 ! Inverse (and possible vertical interpolation to model levels)
1678 ! not needed if not in analysis mode
1679
1680 if (trim(bcsc_mode) == 'BackgroundCheck') cycle
1681
1682 if (varIndex == 1) then
1683 allocate(lfound(numvartot))
1684 call readcorns(lfound,0,'CORVERTI')
1685 end if
1686
1687 if (.not.lfound(varIndex)) then
1688
1689 write(*,*) 'bcsc_corvertSetup: Calculating CORVERTI ', &
1690 'for varIndex =',varIndex
1691
1692 allocate(eigenvec(jnum,jnum),result(jnum,jnum))
1693 eigenvec(1:jnum,1:jnum)=bgStats%corvert(1:jnum,1:jnum,varIndex)
1694
1695 ! CALCULATE EIGENVALUES AND EIGENVECTORS.
1696 ilwork = 4*jnum*2
1697 call dsyev('V','U',jnum,eigenvec,jnum,eigenval,zwork,ilwork,info)
1698 if (info /= 0) then
1699 write(*,*) 'bcsc_corvertSetup: non-zero value of info =',info, &
1700 ' returned by dsyev for wavenumber ',jn
1701 call utl_abort('bcsc_corvertSetup')
1702 end if
1703
1704 ! Set selected number of eigenmodes to zero
1705 if (numModeZero > 0) then
1706 do jk1 = 1, numModeZero
1707 eigenval(jk1) = 0.0d0
1708 end do
1709 write(*,*) 'bcsc_corvertSetup: modified eigenvalues=',eigenval(:)
1710 end if
1711
1712 ! E * lambda^{-1}
1713 eigenvalmax=maxval(eigenval(1:jnum))
1714 do jk1 = 1, jnum
1715 !! if (eigenval(jk1) > 1.0d-8*eigenvalmax) then
1716 ! Threshold increased below to reduce numerical error effects
1717 ! relevant when using corverti in applications
1718 ! (Threshold factor should not be smaller than 1.0D-4.)
1719 ! As consequence C*C^-1 will differ from the identity matrix
1720 ! while being nearly diagonal.
1721 if (eigenval(jk1) > 1.0D-4*eigenvalmax) then
1722 result(1:jnum,jk1) = eigenvec(1:jnum,jk1)/eigenval(jk1)
1723 else
1724 result(1:jnum,jk1) = 0.0D0
1725 end if
1726 end do
1727
1728 ! E * lambda^{-1} * E^T
1729 bgStats%corverti(1:jnum,1:jnum,varIndex)=0.0D0
1730 do jk1 = 1, jnum
1731 do jk2 = 1, jnum
1732 bgStats%corverti(1:jnum,jk1,varIndex) = &
1733 bgStats%corverti(1:jnum,jk1,varIndex) + result(1:jnum,jk2)* &
1734 eigenvec(jk1,jk2)
1735 end do
1736 end do
1737
1738 ! zr=maxval(abs(bgStats%corverti(1:jnum,1:jnum,varIndex)))
1739 ! where (abs(bgStats%corverti(1:jnum,1:jnum,varIndex)) < 1.E-5*zr) &
1740 ! bgStats%corverti(1:jnum,1:jnum,varIndex)=0.0D0
1741
1742 ! Check inverse (for output when mmpi_myid is 0)
1743 result(1:jnum,1:jnum)=0.0D0
1744 do jk1 = 1, jnum
1745 do jk2 = 1, jnum
1746 result(1:jnum,jk1) = result(1:jnum,jk1) + &
1747 bgStats%corvert(1:jnum,jk2,varIndex)* &
1748 bgStats%corverti(jk2,jk1,varIndex)
1749 end do
1750 end do
1751
1752 cletiket = 'C*C^-1'
1753 clnomvar = bgStats%varNameList(varIndex)
1754
1755 if (mmpi_myid == 0 .and. writePhysSpaceStats) then
1756 ierr = utl_fstecr(result(1:jnum,1:jnum),-32,iulcorvert,0,0,0,jnum,jnum, &
1757 1,0,0,bgStats%ntrunc,'X',clnomvar,cletiket,'X',0,0,0,0,5,.true.)
1758 end if
1759
1760 deallocate(eigenvec,result)
1761 end if
1762
1763 end do
1764
1765 if (mmpi_myid /= 0 .or. .not.WritePhysSpaceStats) return
1766
1767 ierr = fstfrm(iulcorvert)
1768 ierr = fclos(iulcorvert)
1769
1770 jnum=nvlev
1771 call writecorns(0,'CORVERT',nvlev)
1772 if ( allocated(lfound) ) then
1773 if (any(.not.lfound(1:numvartot))) call writecorns(0,'CORVERTI',-1)
1774 deallocate(lfound)
1775 end if
1776
1777 end subroutine bcsc_corvertSetup
1778
1779 !--------------------------------------------------------------------------
1780 ! writecorns
1781 !--------------------------------------------------------------------------
1782 subroutine writecorns(nmat,cletiket,nlev)
1783
1784 implicit none
1785
1786 ! Arguments:
1787 character(len=*), intent(in) :: cletiket
1788 integer, intent(in) :: nmat
1789 integer, intent(in) :: nlev
1790
1791 ! Locals:
1792 integer :: jn, nulcorns,ierr,varIndex,jstart,jnum,numvartot
1793 ! standard file variables
1794 integer :: ip1,ip2,ip3
1795 integer :: idateo, ipak, idatyp
1796 character(len=4) :: clnomvar
1797 integer :: fnom, fstouv, fstfrm, fclos
1798
1799 if(mmpi_myid==0) then
1800
1801 write(*,*) 'WRITECORNS: ', trim(cletiket), ' being written to file ', &
1802 'bCovaarSetupChem_out.fst for number of matrices - 1 =', nmat
1803
1804 nulcorns = 0
1805 ierr = fnom(nulcorns,'bCovarSetupChem_out.fst','STD+RND',0)
1806 ierr = fstouv(nulcorns,'RND')
1807
1808 ipak = -32
1809 idatyp = 5
1810 ip1 = 0
1811 ip2 = 0
1812 ip3 = bgStats%ntrunc
1813 idateo = 0
1814
1815 if (nmat == 0) then
1816 numvartot=bgStats%numvar3d
1817 else
1818 numvartot=bgStats%numvar3d+bgStats%numvar2d
1819 end if
1820
1821 do varIndex = 1, numvartot
1822
1823 if (any(CrossCornsVarKindCH(:) /= '') .and. nmat > 0) then
1824 jstart=1
1825 jnum=bgStats%nkgdim
1826 clnomvar = 'ZZ'
1827 else
1828 jstart = bgStats%nsposit(varIndex)
1829 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1830 if (nlev > 0) jnum=nlev
1831 clnomvar = bgStats%varNameList(varIndex)
1832 end if
1833
1834 if (nmat > 0 ) then
1835 do jn = 0, nmat
1836 ip2 = jn
1837 ierr = utl_fstecr(bgStats%corns(jstart:jstart+jnum-1, &
1838 jstart:jstart+jnum-1,jn),ipak,nulcorns,idateo,0,0,jnum, &
1839 jnum,1,ip1,ip2,ip3,'X',clnomvar,cletiket,'X',0,0,0,0, &
1840 idatyp,.true.)
1841 end do
1842 else
1843 if (trim(cletiket) == 'CORVERT') then
1844 ierr = utl_fstecr(bgStats%corvert(1:jnum,1:jnum,varIndex), &
1845 ipak,nulcorns,idateo,0,0,jnum,jnum,1, &
1846 ip1,ip2,ip3,'X',clnomvar,cletiket,'X',0,0,0,0,idatyp,.true.)
1847 else
1848 ierr = utl_fstecr(bgStats%corverti(1:jnum,1:jnum,varIndex), &
1849 ipak,nulcorns,idateo,0,0,jnum,jnum,1, &
1850 ip1,ip2,ip3,'X',clnomvar,cletiket,'X',0,0,0,0,idatyp,.true.)
1851 end if
1852 end if
1853
1854 if (any(CrossCornsVarKindCH(:) /= '') .and. nmat > 0) exit
1855
1856 end do
1857
1858 ierr = fstfrm(nulcorns)
1859 ierr = fclos(nulcorns)
1860 end if
1861
1862 end subroutine writecorns
1863
1864 !--------------------------------------------------------------------------
1865 ! readcorns
1866 !--------------------------------------------------------------------------
1867 subroutine readcorns(lfound,nmat,cletiket)
1868
1869 implicit none
1870
1871 ! Arguments:
1872 logical, intent(out) :: lfound(:)
1873 integer, intent(in) :: nmat
1874 character(len=*), intent(in) :: cletiket
1875
1876 ! Locals:
1877 integer :: jn, icornskey,varIndex,jnum,jstart,numvartot
1878 real(8), allocatable :: zcornssrc(:,:)
1879 ! standard file variables
1880 integer :: ini,inj,ink
1881 integer :: ip1,ip2,ip3
1882 integer :: idateo
1883 character(len=2) :: cltypvar
1884 character(len=4) :: clnomvar
1885
1886 if (nmat == 0) then
1887 numvartot=bgStats%numvar3d
1888 if (trim(cletiket) == 'CORVERT') then
1889 bgStats%corvert(:,:,:)=0.0d0
1890 else if (trim(cletiket) == 'CORVERTI') then
1891 bgStats%corverti(:,:,:)=0.0d0
1892 end if
1893 else if (nmat == bgStats%ntrunc) then
1894 numvartot=bgStats%numvar3d+bgStats%numvar2d
1895 end if
1896
1897 write(*,*) 'readcorns: ', trim(cletiket), &
1898 ' being searched for number of matrices -1 =',nmat
1899
1900 idateo = -1
1901 ip1 = -1
1902 ip3 = bgStats%ntrunc
1903 cltypvar = 'X'
1904
1905 lfound(:)=.false.
1906 VARCYCLE: do varIndex = 1, numvartot
1907
1908 if (any(CrossCornsVarKindCH(:) /= '') .and. nmat > 0) then
1909 jstart=1
1910 jnum=bgStats%nkgdim
1911 clnomvar = 'ZZ'
1912 else
1913 jstart = bgStats%nsposit(varIndex)
1914 jnum = bgStats%nsposit(varIndex+1)-bgStats%nsposit(varIndex)
1915 clnomvar = bgStats%varNameList(varIndex)
1916 end if
1917 allocate(zcornssrc(jnum,jnum))
1918
1919 do jn = 0, nmat
1920 ip2 = jn
1921
1922 ! Looking for FST record parameters..
1923
1924 icornskey = utl_fstlir(ZCORNSSRC,nulbgst,INI,INJ,INK,idateo,cletiket, &
1925 ip1,ip2,ip3,cltypvar,clnomvar)
1926
1927 if(icornskey < 0 ) then
1928 write(*,*) 'readcorns: matrix not found in stats file for variable ', &
1929 clnomvar
1930 deallocate(zcornssrc)
1931 cycle VARCYCLE
1932 end if
1933
1934 if (ini /= jnum .or. inj /= jnum) then
1935 call utl_abort('readcorns: BG stat levels inconsitencies')
1936 end if
1937
1938 if (nmat > 0) then
1939 if (jn == 0) then
1940 bgStats%corns(jstart:jstart+jnum-1,jstart:jstart+jnum-1,:)=0.0d0
1941 end if
1942 bgStats%corns(jstart:jstart+jnum-1,jstart:jstart+jnum-1,jn)= &
1943 zcornssrc(1:jnum,1:jnum)
1944 else if (trim(cletiket) == 'CORVERT') then
1945 bgStats%corvert(1:jnum,1:jnum,varIndex)=zcornssrc(1:jnum,1:jnum)
1946 else
1947 bgStats%corverti(1:jnum,1:jnum,varIndex)=zcornssrc(1:jnum,1:jnum)
1948 end if
1949 end do
1950 lfound(varIndex)=.true.
1951
1952 deallocate(zcornssrc)
1953
1954 if (any(CrossCornsVarKindCH(:) /= '') .and. nmat > 0) exit
1955
1956 end do VARCYCLE
1957
1958 end subroutine readcorns
1959
1960 !--------------------------------------------------------------------------
1961 ! gaspariCohn
1962 !--------------------------------------------------------------------------
1963 function gaspariCohn(ztlen,zr)
1964 implicit none
1965
1966 ! Arguments:
1967 real(8), intent(in) :: ztlen,zr
1968 ! Result:
1969 real(8) :: gasparicohn
1970
1971 ! Locals:
1972 real(8) :: zlc
1973
1974 zlc = ztlen/2.0d0
1975 if(zr <= zlc) then
1976 gasparicohn = -0.250d0*(zr/zlc)**5+0.5d0*(zr/zlc)**4 &
1977 +0.625d0*(zr/zlc)**3-(5.0d0/3.0d0)*(zr/zlc)**2+1.0d0
1978 elseif(zr <= (2.0d0*zlc)) then
1979 gasparicohn = (1.0d0/12.0d0)*(zr/zlc)**5-0.5d0*(zr/zlc)**4 &
1980 +0.625d0*(zr/zlc)**3+(5.0d0/3.0d0)*(zr/zlc)**2 &
1981 -5.0d0*(zr/zlc)+4.0d0-(2.0d0/3.0d0)*(zlc/zr)
1982 else
1983 gasparicohn = 0.0d0
1984 end if
1985 if(gasparicohn < 0.0d0) gasparicohn = 0.0d0
1986
1987 end function gaspariCohn
1988
1989 !--------------------------------------------------------------------------
1990 ! bcsc_finalize
1991 !--------------------------------------------------------------------------
1992 subroutine bcsc_finalize()
1993 implicit none
1994
1995 if(bgStats%initialized) then
1996 if (allocated(bgStats%nsposit)) deallocate(bgStats%nsposit)
1997 if (allocated(bgStats%vlev)) deallocate(bgStats%vlev)
1998 if (allocated(bgStats%lat)) deallocate(bgStats%lat)
1999 if (allocated(bgStats%lon)) deallocate(bgStats%lon)
2000 if (allocated(bgStats%stddev)) deallocate(bgStats%stddev)
2001 if (allocated(bgStats%corns)) deallocate(bgStats%corns)
2002 if (allocated(bgStats%hcorrlen)) deallocate(bgStats%hcorrlen)
2003 if (allocated(bgStats%corvert)) deallocate(bgStats%corvert)
2004 if (allocated(bgStats%corverti)) deallocate(bgStats%corverti)
2005 end if
2006
2007 end subroutine bcsc_finalize
2008
2009 !--------------------------------------------------------------------------
2010 ! bcsc_getCovarCH
2011 !--------------------------------------------------------------------------
2012 subroutine bcsc_getCovarCH(bgStatsOut,transformVarKind_opt)
2013 !
2014 !:Purpose: Pass on covariances in bgStats
2015 !
2016 implicit none
2017
2018 ! Arguments:
2019 type(struct_bcsc_bgStats), intent(out) :: bgStatsOut ! Structure with covariance elements
2020 character(len=*), optional, intent(out) :: TransformVarKind_opt ! Name of variable transform to apply to chemistry variables
2021
2022 if (present(TransformVarKind_opt)) TransformVarKind_opt=TransformVarKindCH
2023
2024 if(.not.bgStats%initialized) then
2025 bgStatsOut%initialized=.false.
2026 call utl_abort('bcsc_getCovarCH: Covariances not set up.')
2027 end if
2028
2029 bgStatsOut=bgStats
2030
2031 end subroutine bcsc_getCovarCH
2032
2033 !--------------------------------------------------------------------------
2034 ! bcsc_resetCorvert
2035 !--------------------------------------------------------------------------
2036 subroutine bcsc_resetCorvert(nlev_T,vlev_T)
2037 !
2038 !:Purpose: Vertically interpolate error correlation matrix fields to
2039 ! generate approximate matrices/vectors in trial field vertical
2040 ! levels via interpolation. No need to make matrix positive
2041 ! definite for this approximation.
2042 !
2043 implicit none
2044
2045 ! Arguments:
2046 integer, intent(in) :: nlev_T ! Number of vertical levels for trial fields
2047 real(8), pointer, intent(in) :: vlev_T(:) ! Trial field vertical levels
2048
2049 ! Locals:
2050 integer :: ilev1,ilev2,j,d1,d2
2051 real(8) :: dz
2052 real(8), allocatable :: wtemp(:,:,:)
2053
2054 d2=0
2055 d1=nlev_T-bgStats%nlev
2056 if (d1 < 0) then
2057 d1=0
2058 d2=d1
2059 end if
2060
2061 allocate(wtemp(nlev_T, nlev_T,bgStats%numvar3d))
2062 wtemp(:,:,:)=0.0d0
2063
2064 ! Apply interpolation
2065
2066 do ilev1 = 1, nlev_T
2067 if (bgStats%vlev(1) >= vlev_T(ilev1)) then
2068 wtemp(ilev1,1:bgStats%nlev+d2,1:bgStats%numvar3d)= &
2069 bgStats%corvert(1,1:bgStats%nlev+d2,1:bgStats%numvar3d)
2070 wtemp(1:bgStats%nlev+d2,ilev1,1:bgStats%numvar3d)= &
2071 bgStats%corvert(1:bgStats%nlev+d2,1,1:bgStats%numvar3d)
2072 else if (bgStats%vlev(bgStats%nlev) <= vlev_T(ilev1)) then
2073 wtemp(ilev1,1+ilev1-bgStats%nlev-d2:ilev1,1:bgStats%numvar3d)= &
2074 bgStats%corvert(bgStats%nlev,1-d2:bgStats%nlev, &
2075 1:bgStats%numvar3d)
2076 wtemp(1+ilev1-bgStats%nlev-d2:ilev1,ilev1,1:bgStats%numvar3d)= &
2077 bgStats%corvert(1-d2:bgStats%nlev,bgStats%nlev, &
2078 1:bgStats%numvar3d)
2079 else
2080 do ilev2=1,bgStats%nlev-1
2081 if (vlev_T(ilev1) >= bgStats%vlev(ilev2) .and. &
2082 vlev_T(ilev1) < bgStats%vlev(ilev2+1)) exit
2083 end do
2084
2085 dz=log(vlev_T(ilev1)/bgStats%vlev(ilev2)) &
2086 /log(bgStats%vlev(ilev2+1)/bgStats%vlev(ilev2))
2087
2088 do j=1-max(ilev1-d1,1),nlev_T-min(d1+ilev1,nlev_T)
2089 wtemp(ilev1,ilev1+j,1:bgStats%numvar3d)= &
2090 (bgStats%corvert(ilev2,ilev2+j,1:bgStats%numvar3d) &
2091 *(1.0-dz) + bgStats%corvert(ilev2+1,ilev2+1+j, &
2092 1:bgStats%numvar3d)*dz)
2093 wtemp(ilev1+j,ilev1,1:bgStats%numvar3d)= &
2094 wtemp(ilev1,ilev1+j,1:bgStats%numvar3d)
2095 end do
2096 end if
2097 end do
2098
2099 bgStats%corvert(1:nlev_T,1:nlev_T,:) = wtemp(:,:,:)
2100
2101 deallocate(wtemp)
2102
2103 end subroutine bcsc_resetCorvert
2104
2105 !--------------------------------------------------------------------------
2106 ! bcsc_resetHcorrlen
2107 !--------------------------------------------------------------------------
2108 subroutine bcsc_resetHcorrlen(nlev_T,vlev_T)
2109 !
2110 !:Purpose: To interpolate horizontal correlation length
2111 ! to trial field vertical levels.
2112 !
2113 implicit none
2114
2115 ! Arguments:
2116 integer, intent(in) :: nlev_T ! Number of target vertical levels
2117 real(8), intent(in) :: vlev_T(:) ! Target vertical levels
2118
2119 ! Locals:
2120 real(8) :: hcorrlen(nlev_T,bgStats%numvar3d) ! correlation lengths
2121 integer :: nlev,ilev1,ilev2
2122 real(8) :: dz
2123
2124 nlev = bgStats%nlev
2125
2126 ! Apply interpolation
2127
2128 do ilev1=1, nlev_T
2129 if (bgStats%vlev(1) >= vlev_T(ilev1)) then
2130 hcorrlen(ilev1,:)=bgStats%hcorrlen(1,1:bgStats%numvar3d)
2131 else if (bgStats%vlev(nlev) <= vlev_T(ilev1)) then
2132 hcorrlen(ilev1,:)=bgStats%hcorrlen(nlev,1:bgStats%numvar3d)
2133 else
2134 do ilev2=1,nlev-1
2135 if (vlev_T(ilev1) >= bgStats%vlev(ilev2) .and. &
2136 vlev_T(ilev1) < bgStats%vlev(ilev2+1)) exit
2137 end do
2138 dz=log(vlev_T(ilev1)/ bgStats%vlev(ilev2)) &
2139 /log( bgStats%vlev(ilev2+1)/ bgStats%vlev(ilev2))
2140 hcorrlen(ilev1,:)=bgStats%hcorrlen(ilev2,1:bgStats%numvar3d)* &
2141 (1.0-dz)+bgStats%hcorrlen(ilev2+1, &
2142 1:bgStats%numvar3d)*dz
2143 end if
2144 end do
2145
2146 bgStats%hcorrlen(1:nlev_T,1:bgStats%numvar3d) = hcorrlen(:,:)
2147
2148 end subroutine bcsc_resetHcorrlen
2149
2150 !--------------------------------------------------------------------------
2151 ! bcsc_StatsExistForVarName
2152 !--------------------------------------------------------------------------
2153 logical function bcsc_StatsExistForVarName(varName)
2154 !
2155 !:Purpose: To check whether covariances have been made available for the
2156 ! specified variable
2157 !
2158 implicit none
2159
2160 ! Arguments:
2161 character(len=*), intent(in) :: varName
2162
2163 if (allocated(bgStats%varNameList)) then
2164 if (any(bgStats%varNameList(:) == varName)) then
2165 bcsc_StatsExistForVarName = .true.
2166 else
2167 bcsc_StatsExistForVarName = .false.
2168 end if
2169 else
2170 bcsc_StatsExistForVarName = .false.
2171 end if
2172
2173 end function bcsc_StatsExistForVarName
2174
2175 !--------------------------------------------------------------------------
2176 ! bcsc_getBgStddev
2177 !--------------------------------------------------------------------------
2178 subroutine bcsc_getBgStddev(varName,maxsize,xlat,xlong,stddevOut,vlev_opt)
2179 !
2180 !:Purpose: To interpolate error std. dev. to obs location.
2181 !
2182 implicit none
2183
2184 ! Arguments:
2185 character(len=*), intent(in) :: varName ! Variable name
2186 integer, intent(in) :: maxsize ! Max array size
2187 real(8), intent(in) :: xlat ! Target latitude
2188 real(8), intent(in) :: xlong ! Target longitude
2189 real(8), intent(out) :: stddevOut(:) ! Error std. dev.
2190 real(8), optional, intent(in) :: vlev_opt(:) ! Target vertical levels
2191
2192 ! Locals:
2193 integer :: varIndex,latIndex,lonIndex,levIndex,nlev,startPosition
2194 real(8) :: work(maxsize,2),zc1,zc2,rlat1,rlat2,rlong1,rlong2,zd1,zd2
2195 real(8) :: stddev_max
2196 integer :: ilev1,ilev2
2197 real(8) :: dz
2198 integer, parameter :: itype=0
2199
2200 if (.not.bgStats%initialized) return
2201
2202 ! Determine location and size of background error std. dev.
2203
2204 do varIndex = 1, bgStats%numvar3d+bgStats%numvar2d
2205 if (trim(varName) == trim(bgStats%varNameList(varIndex))) then
2206 if (varIndex <= bgStats%numvar3d ) then
2207 nlev = bgStats%nlev
2208 else
2209 nlev = 1
2210 end if
2211 startPosition = bgStats%nsposit(varIndex)
2212 exit
2213 end if
2214 end do
2215 if (varIndex > bgStats%numvar3d+bgStats%numvar2d) &
2216 call utl_abort('bcsc_getbgStddev: Variable not found')
2217
2218 if (.not.present(vlev_opt) .and. nlev /= maxsize ) then
2219 write(*,*) 'nlev, maxsize: ',nlev,maxsize
2220 call utl_abort('bcsc_getbgStddev: Inconsistent size')
2221 end if
2222
2223 ! Determine reference longitude index
2224
2225 lonIndex = 2
2226 do while (xlong > bgStats%lon(lonIndex) .and. lonIndex < bgStats%ni+1)
2227 lonIndex = lonIndex+1
2228 end do
2229
2230 ! Set interpolation weights
2231
2232 rlong2 = bgStats%lon(lonIndex)
2233 rlong1 = bgStats%lon(lonIndex-1)
2234
2235 zd2 = (xlong-rlong1)/(rlong2-rlong1)
2236 zd1 = 1.0-zd2
2237
2238 ! Determine reference latitude index
2239
2240 latIndex = 2
2241 do while (xlat > bgStats%lat(latIndex) .and. latIndex < bgStats%nj)
2242 latIndex = latIndex+1
2243 end do
2244
2245 ! Set interpolation weights
2246
2247 rlat2 = bgStats%lat(latIndex)
2248 rlat1 = bgStats%lat(latIndex-1)
2249
2250 zc2 = (xlat-rlat1)/(rlat2-rlat1)
2251 zc1 = 1.0-zc2
2252
2253 if (xlat <= rlat1) then
2254 zc1 = 1.0
2255 zc2 = 0.0
2256 else if (xlat >= rlat2) then
2257 zc1 = 0.0
2258 zc2 = 1.0
2259 end if
2260
2261 ! Apply interpolation
2262
2263 if (itype == 0) then
2264
2265 ! Interpolation of variances and take square root
2266
2267 work(1:nlev,1) = zc1*bgStats%stddev(lonIndex-1,latIndex-1, &
2268 startPosition:startPosition-1+nlev)**2 + &
2269 zc2*bgStats%stddev(lonIndex-1,latIndex, &
2270 startPosition:startPosition-1+nlev)**2
2271
2272 work(1:nlev,2) = zc1*bgStats%stddev(lonIndex,latIndex-1, &
2273 startPosition:startPosition-1+nlev)**2 + &
2274 zc2*bgStats%stddev(lonIndex,latIndex, &
2275 startPosition:startPosition-1+nlev)**2
2276
2277 work(1:nlev,1) = zd1*work(1:nlev,1) + zd2*work(1:nlev,2)
2278
2279 if (nlev /= maxsize) then
2280 do ilev1=1, maxsize
2281 if (bgStats%vlev(1) >= vlev_opt(ilev1)) then
2282 stddevOut(ilev1)=sqrt(work(1,1))
2283 else if (bgStats%vlev(nlev) <= vlev_opt(ilev1)) then
2284 stddevOut(ilev1)=sqrt(work(nlev,1))
2285 else
2286 do ilev2=1,nlev-1
2287 if (vlev_opt(ilev1) >= bgStats%vlev(ilev2) .and. &
2288 vlev_opt(ilev1) < bgStats%vlev(ilev2+1)) exit
2289 end do
2290 dz=log(vlev_opt(ilev1)/ bgStats%vlev(ilev2)) &
2291 /log( bgStats%vlev(ilev2+1)/ bgStats%vlev(ilev2))
2292 stddevOut(ilev1)=sqrt(work(ilev2,1)*(1.0-dz)+work(ilev2+1,1)*dz)
2293 end if
2294 end do
2295 else
2296 stddevOut(1:nlev) = sqrt(work(1:nlev,1))
2297 end if
2298
2299 else
2300
2301 ! Interpolation of std. dev. (to reduce execution time)
2302
2303 work(1:nlev,1) = zc1*bgStats%stddev(lonIndex-1,latIndex-1, &
2304 startPosition:startPosition-1+nlev)**2 + &
2305 zc2*bgStats%stddev(lonIndex-1,latIndex, &
2306 startPosition:startPosition-1+nlev)**2
2307
2308 work(1:nlev,2) = zc1*bgStats%stddev(lonIndex,latIndex-1, &
2309 startPosition:startPosition-1+nlev)**2 + &
2310 zc2*bgStats%stddev(lonIndex,latIndex, &
2311 startPosition:startPosition-1+nlev)**2
2312
2313 work(1:nlev,1) = zd1*work(1:nlev,1) + zd2*work(1:nlev,2)
2314
2315 if (nlev /= maxsize) then
2316 do ilev1=1, maxsize
2317 if (bgStats%vlev(1) >= vlev_opt(ilev1)) then
2318 stddevOut(ilev1)=work(1,1)
2319 else if (bgStats%vlev(nlev) <= vlev_opt(ilev1)) then
2320 stddevOut(ilev1)=work(nlev,1)
2321 else
2322 do ilev2=1,nlev-1
2323 if (vlev_opt(ilev1) >= bgStats%vlev(ilev2) .and. &
2324 vlev_opt(ilev1) < bgStats%vlev(ilev2+1)) exit
2325 end do
2326 dz=log(vlev_opt(ilev1)/bgStats%vlev(ilev2)) &
2327 /log(bgStats%vlev(ilev2+1)/bgStats%vlev(ilev2))
2328 stddevOut(ilev1)=work(ilev2,1)*(1.0-dz)+work(ilev2+1,1)*dz
2329 end if
2330 end do
2331 else
2332 stddevOut(1:nlev) = work(1:nlev,1)
2333 end if
2334
2335 end if
2336
2337 stddev_max = maxval(bgStats%stddev(lonIndex-1:lonIndex, &
2338 latIndex-1:latIndex, &
2339 startPosition:startPosition-1+nlev))
2340 do levIndex = 1, nlev
2341 if (stddevOut(levIndex) < 0.0 .or. stddevOut(levIndex) > 1.1*stddev_max) then
2342 write(*,*) 'bcsc_getbgStddev: Interpolated std dev incorrect:'
2343 write(*,*) 'bcsc_getbgStddev: zc1,zc2,zd1,zd2 = ',zc1,zc2,zd1,zd2
2344 write(*,*) 'bcsc_getbgStddev: stddevOut,stddev_max = ', &
2345 stddevOut(levIndex),stddev_max
2346 write(*,*) 'bcsc_getbgStddev: lonIndex,latIndex,j,xlong,xlat = ', &
2347 lonIndex,latIndex,levIndex,xlong,xlat
2348 write(*,*) 'bcsc_getbgStddev: rlong2,rlong1,rlat1,rlat2 = ', &
2349 rlong2,rlong1,rlat1,rlat2
2350 call utl_abort('bcsc_getbgStddev')
2351 end if
2352 end do
2353
2354 end subroutine bcsc_getBgStddev
2355
2356 !--------------------------------------------------------------------------
2357 ! bcsc_addBgStddev
2358 !--------------------------------------------------------------------------
2359 subroutine bcsc_addBgStddev(headerIndex,stddevIn,obsdataMaxsize)
2360 !
2361 !:Purpose: To add background stddev profiles (and inverse) to bgStddev
2362 ! which can be retrieved later using a header index.
2363 !
2364 implicit none
2365
2366 ! Arguments:
2367 integer, intent(in) :: headerIndex
2368 real(8), intent(in) :: stddevIn(:,:)
2369 integer, intent(in) :: obsdataMaxsize
2370
2371 if (.not.associated(bgStddev%data2d)) then
2372 call oss_obsdata_alloc(bgStddev, obsdataMaxsize, dim1= &
2373 size(stddevIn,dim=1), dim2_opt=size(stddevIn,dim=2))
2374 bgStddev%nrep = 0
2375 end if
2376
2377 ! In this case nrep will count the number of filled reps in the data arrays
2378 bgStddev%nrep = bgStddev%nrep+1
2379
2380 if (bgStddev%nrep > obsdataMaxsize) then
2381 call utl_abort('bcsc_addbgStddev: Reached max size of array ' // &
2382 trim(utl_str(obsdataMaxsize)) )
2383 end if
2384
2385 ! Use the header number as the unique code for this obs data
2386 write(bgStddev%code(bgStddev%nrep),'(I22)') headerIndex
2387
2388 bgStddev%data2d(:,1,bgStddev%nrep) = stddevIn(:,1)
2389
2390 where (stddevIn(:,1) > 0.0D0)
2391 bgStddev%data2d(:,2,bgStddev%nrep) = 1.0D0/stddevIn(:,1)
2392 elsewhere
2393 bgStddev%data2d(:,2,bgStddev%nrep) = 0.0D0
2394 end where
2395
2396 end subroutine bcsc_addBgStddev
2397
2398 !--------------------------------------------------------------------------
2399 ! bcsc_retrievebgStddev
2400 !--------------------------------------------------------------------------
2401 function bcsc_retrieveBgStddev(dim1,dim2,headerIndex) result(stddevOut)
2402 !
2403 !:Purpose: To retrieve previously saved background stddev profiles
2404 ! in bgStddev from the header index of the observation.
2405 !
2406 implicit none
2407
2408 ! Arguments:
2409 integer, intent(in) :: dim1 ! Dimensions of output array
2410 integer, intent(in) :: dim2 ! Dimensions of output array
2411 integer, intent(in) :: headerIndex ! Index of observation
2412 ! Result:
2413 real(8) :: stddevOut(dim1,dim2)
2414
2415 ! Locals:
2416 character(len=22) :: code
2417
2418 if (bgStddev%dim1 /= dim1 .or. bgStddev%dim2 /= dim2) then
2419 call utl_abort("bcsc_retrievebgStddev: Inconsitent dimensions")
2420 end if
2421
2422 write(code,'(I22)') headerIndex
2423
2424 stddevOut = oss_obsdata_get_array2d(bgStddev,code)
2425
2426 end function bcsc_retrieveBgStddev
2427
2428end module bCovarSetupChem_mod
