1module obsSpaceDiag_mod
2 ! MODULE obsSpaceDiag_mod (prefix='osd' category='1. High-level functionality')
3 !
4 !:Purpose: Some experimental procedures for computing various diagnostics in
5 ! observation space.
6 !
7 use codePrecision_mod
8 use midasMpi_mod
9 use bufr_mod
10 use codtyp_mod
11 use MathPhysConstants_mod
12 use horizontalCoord_mod
13 use timeCoord_mod
14 use controlVector_mod
15 use obsSpaceData_mod
16 use columnData_mod
17 use verticalCoord_mod
18 use gridStateVector_mod
19 use bMatrix_mod
20 use bMatrixHi_mod
21 use bMatrixEnsemble_mod
22 use bCovarSetupChem_mod
23 use varNameList_mod
24 use stateToColumn_mod
25 use randomNumber_mod
26 use obsOperators_mod
27 use utilities_mod
28 use physicsfunctions_mod
29 use obsSubSpaceData_mod
30 use obsfiles_mod
31 use obsOperatorsChem_mod
32 use obsFamilyList_mod
33
34 implicit none
35 save
36 private
37
38 ! public procedures
39 !------------------
40 public :: osd_ObsSpaceDiag
41
42 type :: struct_osd_diagn
43
44 ! Structure for storing observation-space diagnostic arrays
45 !
46 ! Real arrays indexed by (lat,lon,lev,stat), where stat=1 for RMS and stat=2 for mean for all except Jo_stats.
47 ! For Jo_stats, stat=1 is Jo for x=x_analysis and stat=2 is Jo for x=x_background.
48 ! The array counts is indexed by (lat,lon,lev).
49 ! The array status_count is indexed (lat,lon,lev,status), where the status number ranges from 0 to 2.
50 !
51 ! Variable Description
52 ! -------- -----------
53 ! OmP_stats obs - background statistics
54 ! OmA_stats obs - analysis statistics
55 ! obs_stats observation statistics
56 ! Jo_stats cost function statistics
57 ! - First four elements (of last dimension) apply only prescribed obs error variances
58 ! as normalizing-scaling denomicators.
59 ! - Fifth element (of last dimension) applies the sum of the prescribed obs error variances
60 ! and the diag(HPHT) as normalizing-scaling denominator. This does not include consideration
61 ! of spatial correlations of (O-P) between obs points
62 ! associated to HPHT in the normalizing denominator.
63 ! Jpa_stats statistics (P-A) in observation space.
64 ! - Not exactly equivalent to Jb of the cost function
65 ! - Applies the prescribed diag(HPHT) as normalizing-scaling denomiator.
66 ! - Does not include consideration of spatial correlations of (P-A) between obs points
67 ! associated to HPHT in normalizing denominator
68 ! - Vert. coordinate binning included but not currently output.
69 ! diagR_stats Elements for the calc of mean{[(O-P)-mean(O-P)][(O-A)-mean(O-A)]}
70 ! (with each O-P and O-A difference divided by sigma_obs)
71 ! for scaling factor adjustments of obs std. dev via the Desroziers approach.
72 ! diagHPHT_stats Elements for the calc of mean{[(O-P)-mean(O-P)][[(O-P)-mean(O-P)]-[(O-A)-mean(O-A)]]}
73 ! = mean{[(O-P)-mean(O-P)][(A-P)-mean(A-P)]}(with each O-P and O-A difference divided by sqrtHPHT)
74 ! for scaling factor adjustments of background std. dev. in obs space via the Desroziers approach.
75 ! counts number of observations in a (lat,lon,lev) bin
76 ! nlat number of latitude levels
77 ! nlon number of longitude levels
78 ! nlev number of vertical levels
79 ! nbin total number of (lat,lon,lev) bins
80 ! nstat number of different statistics types
81 ! nstatus indicates the number of observations with a certain status,
82 ! with status values denoting:
83 ! 0 - observation has been rejected and not included for diagnostics
84 ! 1 - observation has been assimilated
85 ! 2 - observation has been used for diagnostics only (not assimilated)
86 ! deltaLat latitude bin size (deg)
87 ! deltaLon longitude bin size (deg)
88 ! deltaLogPressure vertical bin size in log(pressure)
89 ! allow_print_summary indicates fs printing of summary diagnostics is allowed
90 ! assim_mode indicates if assimilation was performed for this dataset
91
92 real(8), allocatable :: OmP_stats(:,:,:,:),OmA_stats(:,:,:,:),obs_stats(:,:,:,:),Jo_stats(:,:,:,:)
93 real(8), allocatable :: diagR_stats(:,:,:,:),diagHPHT_stats(:,:,:,:), Jpa_stats(:,:,:)
94 integer, allocatable :: counts(:,:,:),nstatus(:,:,:,:)
95
96 integer :: nlev,nlat,nlon,nbin,nstat
97 real(8) :: deltaLat
98 real(8) :: deltaLon
99 real(8) :: deltaLogPressure
100 logical :: allow_print_summary,assim_mode
101
102 end type struct_osd_diagn
103
104 ! Parameters identifying obs sets and related actions for diagnostic calcs of each family
105 !
106 ! diagn_num Prescribed (starting) number of (stnid, bufr, nlev) for the diagnostics calc
107 !
108 ! diagn_stnid Prescribed (starting) list of stnid (with *s as needed) for the diagnostics calc
109 ! with '*' denoting wild cards
110 !
111 ! diagn_varno Prescribed (starting) list of bufr elements for the diagnostics calc
112 !
113 ! diagn_unilev Prescribed (starting) list of logicals indicating uni-level obs for the diagnostics calc
114 !
115 ! diagn_pressmin Bottom of top layer for diagnostics (in Pa).
116 !
117 ! diagn_save Logical indicating gridded diagnostics are to be saved
118 ! in an ascii file in addition to overall diagnostics.
119 !
120 ! diagn_nset Integer indicating grouping of diagnostics with
121 ! 1: group by stnid
122 ! 2: group by (stnid,bufr)
123 ! 3: group by (stnid,bufr,nlev)
124 !
125 ! diagn_all Logical indicating if all combinations specified by diagn_nset are to be
126 ! used in diagnostics or only those specified by the diagn_* arrays
127 !
128 ! obsspace_diagn_filename
129 ! File name for file containing obs space diagnostics related to chemical constituents.
130
131 character(len=100) :: obsspace_diagn_filename(ofl_numFamily)
132 integer, parameter :: max_cfg_size=100
133
134 ! Namelist variables
135 real(8) :: deltaLat ! Size of latitude bins for diagnostics (degrees)
136 real(8) :: deltaLon ! Size of longitude bins for diagnostics (degrees)
137 real(8) :: deltaPressure ! Size of vertical bins for diagnostics (Pascal)
138 real(8) :: deltaHeight ! Size of vertical bins for diagnostics (meters)
139 integer :: numFamily ! MUST NOT BE INCLUDED IN NAMELIST!
140 character(len=2) :: familyList(ofl_numFamily) ! List of family names to consider
141 integer :: numElement ! MUST NOT BE INCLUDED IN NAMELIST!
142 integer :: elementList(ofl_numFamily) ! List of bufr element ids to consider
143 integer :: nrandseed ! Random seed value
144 logical :: lrandom ! Perform diagnostics from random perturbations
145 integer :: diagn_num(ofl_numFamily) ! Prescribed (starting) number of (stnid, bufr, nlev) for diag calc
146 integer :: diagn_nset(ofl_numFamily) ! Choose to group by 1: stnid, 2: stnid,bufr, 3: stnid,bufr,nlev
147 integer :: diagn_varno(ofl_numFamily,max_cfg_size) ! List of bufr element ids for diag calc
148 character(len=9) :: diagn_stnid(ofl_numFamily,max_cfg_size) ! List of stnid for diag calc
149 logical :: diagn_save(ofl_numFamily) ! Choose to save gridded info in ascii file
150 logical :: diagn_all(ofl_numFamily) ! Choose to use all combinations specified by diagn_nset
151 logical :: diagn_unilev(ofl_numFamily,max_cfg_size) ! List indicating if uni-lev obs for diag calc
152 real(8) :: diagn_pressmin(ofl_numFamily) ! Bottom of top layer for diagnostics (in Pa)
153
154contains
155
156 !--------------------------------------------------------------------------
157 ! osd_ObsSpaceDiag
158 !--------------------------------------------------------------------------
159 subroutine osd_ObsSpaceDiag( obsSpaceData, columnTrlOnAnlIncLev, hco_anl, analysisMode_opt )
160 !
161 !:Purpose: Calls routines to perform observation-space diagnostic tasks
162 !
163 !:Arguments:
164 ! :obsSpaceData: Obs space data structure
165 ! :columnTrlOnAnlIncLev: Structure of vertical columns at obs locations.
166 ! Expected to be for analysis vertical levels if to be used.
167 ! :analysisMode: logical indicating if following analysis mode or not (optional)
168 ! Assumed .true. if not present.
169 !
170 implicit none
171
172 ! Arguments:
173 type(struct_obs) , intent(inout) :: obsSpaceData
174 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
175 type(struct_hco), pointer, intent(in) :: hco_anl
176 logical, optional, intent(in) :: analysisMode_opt
177
178 ! Locals:
179 logical :: nmlExists,anlm_mod
180 integer :: ierr
181 integer :: dateprnt,timeprnt,newdate
182
183 if (present(analysisMode_opt)) then
184 anlm_mod = analysisMode_opt
185 else
186 anlm_mod = .true.
187 end if
188
189 call osd_setup(nmlExists)
190 ierr = newdate(tim_getDatestamp(),dateprnt,timeprnt,-3)
191 dateprnt=dateprnt*100+timeprnt/1000000
192
193 ! Perform diagnostics based on OmP (and OmA if available)
194
195 call osd_obsPostProc(obsSpaceData,columnTrlOnAnlIncLev,deltaLat,deltaLon,deltaPressure,anlm_mod)
196
197 if ((.not. anlm_mod) .or. (.not.lrandom) .or. (.not.nmlExists)) return
198
199 ! Perform diagnostics from random perturbations
200
201 call osd_calcInflation(obsSpaceData,columnTrlOnAnlIncLev,hco_anl,dateprnt)
202
203 end subroutine osd_ObsSpaceDiag
204
205 !--------------------------------------------------------------------------
206 ! osd_calcInflation
207 !--------------------------------------------------------------------------
208 subroutine osd_calcInflation( obsSpaceData, columnTrlOnAnlIncLev, hco_anl, dateprnt )
209 !
210 !:Purpose: Calculates observation-space diagnostics from random perturbations
211 !
212
213 implicit none
214
215 ! Arguments:
216 type(struct_obs) , intent(inout) :: obsSpaceData
217 type(struct_columnData) , intent(inout) :: columnTrlOnAnlIncLev
218 type(struct_hco), pointer, intent(in) :: hco_anl
219 integer , intent(in) :: dateprnt
220
221 ! Locals:
222 type(struct_gsv) :: statevector
223 type(struct_columnData) :: column
224 type(struct_vco), pointer :: vco_anl
225 real(8), allocatable,target :: controlVector(:)
226 integer :: familyIndex,elementIndex,bodyIndex,headerIndex,latIndex,lonIndex,verticalIndex
227 integer :: maxLat,maxLon,maxVertical
228 real(8), allocatable :: innovStd(:,:,:),bmatHiStd(:,:,:),bmatEnStd(:,:,:)
229 integer, allocatable :: counts(:,:,:)
230 real(8), allocatable :: my_innovStd(:,:,:),my_bmatHiStd(:,:,:),my_bmatEnStd(:,:,:)
231 integer, allocatable :: my_counts(:,:,:)
232 integer :: ierr,nulinnov,nulBmatHi,nulBmatEn,nulcount,fnom,fclos,ivco,ivco_recv,iseed,jj,jlev,jvar
233 integer :: ivar_count,nlev_max
234 logical :: lpert_static, lpert_ens
235 real(8), pointer :: cvBhi(:), cvBen(:), cvBchm(:)
236 real(pre_incrReal), pointer :: field(:,:,:,:)
237 real(8), allocatable :: HxBhi(:), HxBen(:)
238 real(8), allocatable :: scaleFactor(:),scaleFactorChm(:,:)
239 character(len=128) :: innovFileName,bmatHiFileName,bmatEnFileName,countFileName
240 character(len=6) :: elementStr
241 character(len=10) :: dateStr
242
243 write(*,*) 'osd_calcInflation: Starting'
244
245 if( nrandseed == 999 ) nrandseed=dateprnt ! if seed not set by namelist, use valid date/time
246 write(*,*) 'osd_calcInflation: random seed set to ',nrandseed
247
248 maxLat = nint(180.0d0/deltaLat)
249 maxLon = nint(360.0d0/deltaLon)
250 maxVertical = max(1+nint(110000.0d0/deltaPressure),1+nint(80000.0d0/deltaHeight),200)
251
252 write(*,*) 'osd_calcInflation: Compute random realization of background error'
253
254 ! allocate vectors to store Hx for the static and ensemble-based covariance matrices
255 allocate(HxBhi(obs_numbody(obsSpaceData)))
256 allocate(HxBen(obs_numbody(obsSpaceData)))
257
258 ! initialize columnData object for increment
259 call col_setVco(column,col_getVco(columnTrlOnAnlIncLev))
260 call col_allocate(column,col_getNumCol(columnTrlOnAnlIncLev),mpiLocal_opt=.true.)
261
262 ! initialize gridStateVector object for increment
263 vco_anl => col_getVco(columnTrlOnAnlIncLev)
264 call gsv_allocate(statevector, tim_nstepobsinc, hco_anl, vco_anl, &
265 dataKind_opt=pre_incrReal, mpi_local_opt=.true.)
266
267 nlev_max=max(col_getNumLev(columnTrlOnAnlIncLev,'MM'),col_getNumLev(columnTrlOnAnlIncLev,'TH'))
268
269 allocate(controlVector(cvm_nvadim))
270 allocate(scaleFactor(nlev_max))
271 allocate(scaleFactorChm(100,nlev_max))
272
273 ! COMPUTE BMATRIX PERTURBATION FOR THE STATIC COVARIANCES CASE; from Bhi and or BChm
274
275 if (all(familyList(1:numFamily) == 'CH').or.(.not.cvm_subVectorExists('B_HI'))) then
276 cvBhi => null()
277 else
278 cvBhi => cvm_getSubVector(controlVector,'B_HI')
279 end if
280 if (any(familyList(1:numFamily) == 'CH').and.obs_famExist(obsSpaceData,'CH').and.cvm_subVectorExists('B_CHM')) then
281 cvBChm => cvm_getSubVector(controlVector,'B_CHM')
282 else
283 cvBChm => null()
284 end if
285
286 HxBhi(:) = 0.0d0
287
288 iseed = abs(nrandseed)
289 call rng_setup(iseed)
290
291 if (cvm_subVectorExists('B_HI').or.cvm_subVectorExists('B_CHM')) then
292
293 ! compute random control vector
294
295 controlVector(:) = 0.0d0
296
297 if (cvm_subVectorExists('B_HI')) then
298 do jj = 1,size(cvBhi)
299 cvBhi(jj)=rng_gaussian()
300 enddo
301 ! initialize vector of scaleFactors
302 call bhi_getScaleFactor(scaleFactor)
303 else
304 scaleFactor(:)=1.0
305 end if
306
307 if (cvm_subVectorExists('B_CHM')) then
308 do jj = 1,size(cvBChm)
309 cvBChm(jj)=rng_gaussian()
310 enddo
311 ! initialize vector of scaleFactors
312 call bcsc_getScaleFactor(scaleFactorChm)
313 else
314 scaleFactorChm(:,:)=1.0
315 end if
316
317 ! multiply vector by B^1/2
318 call bmat_sqrtB(controlVector,cvm_nvadim,statevector)
319
320 ! undo the scaleFactor (THIS IS NOT CORRECT FOR 2D VARIABLES!!! (P0 and TG) )
321
322 ivar_count=0
323 do jvar=1,vnl_numvarmax
324 if(gsv_varExist(statevector,vnl_varNameList(jvar))) then
325
326 call gsv_getField(statevector,field,vnl_varNameList(jvar))
327
328 if (vnl_varKindFromVarname(vnl_varNameList(jvar)) == 'MT') then
329 do jlev = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
330 if(scaleFactor(jlev) > 0.0d0) field(:,:,jlev,:)=field(:,:,jlev,:)/scaleFactor(jlev)
331 enddo
332 else if (vnl_varKindFromVarname(vnl_varNameList(jvar)) == 'CH') then
333 ivar_count=ivar_count+1
334 do jlev = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
335 if(scaleFactorChm(ivar_count,jlev) > 0.0d0) field(:,:,jlev,:)=field(:,:,jlev,:) &
336 /scaleFactorChm(ivar_count,jlev)
337 end do
338 end if
339 endif
340 enddo
341
342 ! multiply by H
343 call s2c_tl(statevector,column,columnTrlOnAnlIncLev,obsSpaceData) ! put in column H_horiz dx
344 call oop_Htl(column,columnTrlOnAnlIncLev,obsSpaceData,1) ! Save as OBS_WORK: H_vert H_horiz dx = Hdx
345 do bodyIndex=1,obs_numBody(obsSpaceData)
346 HxBhi(bodyIndex) = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
347 enddo
348
349 end if
350
351 ! COMPUTE BMATRIX PERTURBATION FOR THE ENSEMBLE COVARIANCES CASE; from Ben
352
353 if (cvm_subVectorExists('B_ENS')) then
354
355 cvBen => cvm_getSubVector(controlVector,'B_ENS')
356 HxBen(:) = 0.0d0
357
358 ! compute random control vector
359
360 controlVector(:) = 0.0d0
361
362 do jj = 1,size(cvBen)
363 cvBen(jj)=rng_gaussian()
364 enddo
365
366 ! initialize vector of scaleFactors
367 call ben_getScaleFactor(scaleFactor)
368
369 scaleFactorChm(:,:)=1.0
370
371 ! multiply vector by B^1/2
372 call bmat_sqrtB(controlVector,cvm_nvadim,statevector)
373
374 ! undo the scaleFactor
375
376 ivar_count=0
377 do jvar=1,vnl_numvarmax
378 if(gsv_varExist(statevector,vnl_varNameList(jvar))) then
379
380 call gsv_getField(statevector,field,vnl_varNameList(jvar))
381
382 if (vnl_varKindFromVarname(vnl_varNameList(jvar)) == 'MT') then
383 do jlev = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
384 if(scaleFactor(jlev) > 0.0d0) field(:,:,jlev,:)=field(:,:,jlev,:)/scaleFactor(jlev)
385 enddo
386 else if (vnl_varKindFromVarname(vnl_varNameList(jvar)) == 'CH') then
387 ivar_count=ivar_count+1
388 do jlev = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname(vnl_varNameList(jvar)))
389 if(scaleFactorChm(ivar_count,jlev) > 0.0d0) field(:,:,jlev,:)=field(:,:,jlev,:) &
390 /scaleFactorChm(ivar_count,jlev)
391 end do
392 end if
393 endif
394 enddo
395
396 ! multiply vector by H
397 call s2c_tl(statevector,column,columnTrlOnAnlIncLev,obsSpaceData) ! put in column H_horiz dx
398 call oop_Htl(column,columnTrlOnAnlIncLev,obsSpaceData,1) ! Save as OBS_WORK: H_vert H_horiz dx = Hdx
399 do bodyIndex=1,obs_numBody(obsSpaceData)
400 HxBen(bodyIndex) = obs_bodyElem_r(obsSpaceData,OBS_WORK,bodyIndex)
401 enddo
402 else
403 cvBen => null()
404 HxBen(:) = 0.0d0
405 end if
406
407 deallocate(scaleFactor)
408 deallocate(scaleFactorChm)
409 call col_deallocate(column)
410 call gsv_deallocate(statevector)
411
412 allocate(my_innovStd(maxLat,maxLon,maxVertical))
413 allocate(my_bmatHiStd(maxLat,maxLon,maxVertical))
414 allocate(my_bmatEnStd(maxLat,maxLon,maxVertical))
415 allocate(my_counts(maxLat,maxLon,maxVertical))
416
417 allocate(innovStd(maxLat,maxLon,maxVertical))
418 allocate(bmatHiStd(maxLat,maxLon,maxVertical))
419 allocate(bmatEnStd(maxLat,maxLon,maxVertical))
420 allocate(counts(maxLat,maxLon,maxVertical))
421
422 FAMILY: do familyIndex = 1, numFamily
423
424 ELEMENT: do elementIndex = 1, numElement
425
426 ! Initialize logicals for calc of perturbation diagnostics.
427
428 lpert_static=.false.
429 lpert_ens=.false.
430
431 if (familyList(familyIndex) /= 'CH') then
432 if (cvm_subVectorExists('B_HI')) lpert_static=.true.
433 else
434 if (cvm_subVectorExists('B_CHM')) lpert_static=.true.
435 end if
436 if (cvm_subVectorExists('B_ENS')) lpert_ens=.true.
437
438 ivco = -999
439 my_innovStd(:,:,:) = 0.0d0
440 my_bmatHiStd(:,:,:) = 0.0d0
441 my_bmatEnStd(:,:,:) = 0.0d0
442 my_counts(:,:,:) = 0
443
444 call obs_set_current_body_list(obsSpaceData,familyList(familyIndex))
445 BODY: do
446 bodyIndex = obs_getBodyIndex(obsSpaceData)
447 if (bodyIndex < 0) exit BODY
448
449 if(obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == elementList(elementIndex) .and. &
450 obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated) then
451
452 call osd_getIndices(obsSpaceData,bodyIndex,latIndex,lonIndex,verticalIndex)
453
454 if(verticalIndex == -1) then
455 ! skip this obs for whatever reason
456 cycle BODY
457 else if(latIndex > maxLat .or. lonIndex > maxLon .or. verticalIndex > maxVertical) then
458 write(*,*) 'osd_calcInflation: index too big: lat,lon,vertical=',latIndex,lonIndex,verticalIndex, &
459 ' lat_max,lon_max,vertical_max=',maxlat,maxlon,maxvertical
460 call utl_abort('osd_calcInflation')
461 endif
462
463 ivco = obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex)
464 my_counts(latIndex,lonIndex,verticalIndex) = my_counts(latIndex,lonIndex,verticalIndex) + 1
465 my_innovStd(latIndex,lonIndex,verticalIndex) = my_innovStd(latIndex,lonIndex,verticalIndex) + &
466 obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)* &
467 obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)
468 if (lpert_static) my_bmatHiStd(latIndex,lonIndex,verticalIndex) = my_bmatHiStd(latIndex,lonIndex,verticalIndex) + &
469 HxBhi(bodyIndex)*HxBhi(bodyIndex)
470 if (lpert_ens) my_bmatEnStd(latIndex,lonIndex,verticalIndex) = my_bmatEnStd(latIndex,lonIndex,verticalIndex) + &
471 HxBen(bodyIndex)*HxBen(bodyIndex)
472
473 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
474
475 endif
476 enddo BODY
477
478 call rpn_comm_allreduce(ivco,ivco_recv,1,"MPI_INTEGER","MPI_MAX","GRID",ierr)
479 ivco = ivco_recv
480
481 call rpn_comm_allreduce(my_counts,counts,maxLat*maxLon*maxVertical,"MPI_INTEGER","MPI_SUM","GRID",ierr)
482 call rpn_comm_allreduce(my_innovStd,innovStd,maxLat*maxLon*maxVertical,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
483 if (lpert_static) call rpn_comm_allreduce(my_bmatHiStd,bmatHiStd,maxLat*maxLon*maxVertical,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
484 if (lpert_ens) call rpn_comm_allreduce(my_bmatEnStd,bmatEnStd,maxLat*maxLon*maxVertical,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
485
486 where (counts > 0) innovStd = sqrt(innovStd/counts)
487 if (lpert_static) then
488 where (counts > 0) bmatHiStd = sqrt(bmatHiStd/counts)
489 end if
490 if (lpert_ens) then
491 where (counts > 0) bmatEnStd = sqrt(bmatEnStd/counts)
492 end if
493
494 if(mmpi_myid == 0 .and. sum(counts(:,:,:)) > 0) then
495
496 ! determine file names
497 write(dateStr,'(i10.10)') dateprnt
498 write(elementStr,'(i6.6)') elementList(elementIndex)
499 innovFileName = 'innov' // dateStr // '_' // trim(familyList(familyIndex)) // '_' // trim(elementStr) // '.dat'
500 if (lpert_static) bmatHiFileName = 'bmathi' // dateStr // '_' // trim(familyList(familyIndex)) // '_' // trim(elementStr) // '.dat'
501 if (lpert_ens) bmatEnFileName = 'bmaten' // dateStr // '_' // trim(familyList(familyIndex)) // '_' // trim(elementStr) // '.dat'
502 countFileName = 'count' // dateStr // '_' // trim(familyList(familyIndex)) // '_' // trim(elementStr) // '.dat'
503
504 ! open files
505 nulinnov=0
506 nulBmatHi =0
507 nulBmatEn =0
508 nulcount=0
509 ierr = fnom(nulinnov,innovFileName,'FMT+R/W',0)
510 if (lpert_static) ierr = fnom(nulBmatHi ,bmatHiFileName ,'FMT+R/W',0)
511 if (lpert_ens) ierr = fnom(nulBmatEn ,bmatEnFileName ,'FMT+R/W',0)
512 ierr = fnom(nulcount,countFileName,'FMT+R/W',0)
513
514 ! write data for this family/element
515 write(nulinnov,*) '***maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight='
516 write(nulinnov,*) maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight
517
518 if (lpert_static) then
519 write(nulBmatHi,*) '***maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight='
520 write(nulBmatHi,*) maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight
521 end if
522 if (lpert_ens) then
523 write(nulBmatEn,*) '***maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight='
524 write(nulBmatEn,*) maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight
525 end if
526
527 write(nulcount,*) '***maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight='
528 write(nulcount,*) maxLon,maxLat,deltaLon,deltaLat,deltaPressure,deltaHeight
529 do verticalIndex = 1,maxVertical
530 if(sum(counts(:,:,verticalIndex)) > 0) then
531 write(nulinnov,*) '***verticalIndex,vco='
532 write(nulinnov,*) verticalIndex,ivco
533
534 if (lpert_static) then
535 write(nulBmatHi,*) '***verticalIndex,vco='
536 write(nulBmatHi,*) verticalIndex,ivco
537 end if
538 if (lpert_ens) then
539 write(nulBmatEn,*) '***verticalIndex,vco='
540 write(nulBmatEn,*) verticalIndex,ivco
541 end if
542
543 write(nulcount,*) '***verticalIndex,vco='
544 write(nulcount,*) verticalIndex,ivco
545 do latIndex = 1,maxLat
546 write(nulinnov,*) innovStd(latIndex,:,verticalIndex)
547 write(nulcount,*) counts(latIndex,:,verticalIndex)
548 enddo
549
550 if (lpert_static) then
551 do latIndex = 1,maxLat
552 write(nulBmatHi ,*) bmatHiStd(latIndex,:,verticalIndex)
553 enddo
554 end if
555 if (lpert_ens) then
556 do latIndex = 1,maxLat
557 write(nulBmatEn ,*) bmatEnStd(latIndex,:,verticalIndex)
558 enddo
559 end if
560
561 endif
562 enddo
563
564 ! close files
565 ierr = fclos(nulinnov)
566 if (lpert_static) ierr = fclos(nulBmatHi)
567 if (lpert_ens) ierr = fclos(nulBmatEn)
568 ierr = fclos(nulcount)
569 endif
570 enddo ELEMENT
571
572 enddo FAMILY
573
574 deallocate(my_counts)
575 deallocate(my_innovStd)
576 deallocate(my_bmatHiStd)
577 deallocate(my_bmatEnStd)
578
579 deallocate(innovStd)
580 deallocate(bmatHiStd)
581 deallocate(bmatEnStd)
582 deallocate(counts)
583
584 deallocate(HxBhi)
585 deallocate(HxBen)
586
587 write(*,*) 'osd_calcInflation: Finished'
588
589 end subroutine osd_calcInflation
590
591 !--------------------------------------------------------------------------
592 ! osd_getIndices
593 !--------------------------------------------------------------------------
594 subroutine osd_getIndices( obsSpaceData, bodyIndex, latIndex, lonIndex, verticalIndex )
595 !
596 implicit none
597
598 ! Arguments:
599 type(struct_obs), intent(in) :: obsSpaceData
600 integer , intent(in) :: bodyIndex
601 integer , intent(out) :: latIndex
602 integer , intent(out) :: lonIndex
603 integer , intent(out) :: verticalIndex
604
605 ! Locals:
606 real(8), parameter :: epsilon=0.001
607 integer :: headerIndex, bodyElem_i
608
609 ! codtypes for tovs: 164(AMSUA) 168 180 181 182 183 185 186 192 193
610
611 ! epsilon is added below to handle case where lon/dlon~1 or lat/dlat~1
612 headerIndex = obs_bodyElem_i(obsSpaceData,OBS_HIND,bodyIndex)
613 latIndex = 1 + floor((90.0d0 + obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)*MPC_DEGREES_PER_RADIAN_R8)/deltaLat - epsilon)
614 if (latIndex == 0) latIndex=1
615 lonIndex = 1 + floor(obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)*MPC_DEGREES_PER_RADIAN_R8/deltaLon - epsilon)
616 if (lonIndex == 0) lonIndex=1
617
618 select case(obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex))
619 case(1)
620 ! height coordinate
621 verticalIndex = 1 + nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)/deltaHeight)
622 case(2)
623 ! pressure coordinate
624 verticalIndex = 1 + nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)/deltaPressure)
625 case(3)
626 ! channel number
627 verticalIndex = nint(obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex))
628 if(obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex) == codtyp_get_codtyp("AMSUA")) then
629 ! amsu-a
630 verticalIndex = verticalIndex - 27
631 else
632 ! ignore other types of TOVS for now
633 verticalIndex = -1
634 endif
635 case(4,5)
636 ! Integrated column or surface value - assign to first level
637 verticalIndex = 1
638 case default
639 ! unknown vertical coordinate
640 bodyElem_i = obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex)
641 write(*,*) 'osd_getIndices: Unknown VCO! ', bodyElem_i
642 verticalIndex = -1
643 end select
644
645 end subroutine osd_getIndices
646
647 !--------------------------------------------------------------------------
648 ! osd_setup
649 !--------------------------------------------------------------------------
650 subroutine osd_setup(nmlExists)
651 !
652 implicit none
653
654 ! Arguments:
655 logical, intent(out) :: nmlExists
656
657 ! Locals:
658 integer :: nulnam,ierr,fnom,fclos
659 integer :: familyIndex, elementIndex
660
661 namelist /namosd/nrandseed,deltaLat,deltaLon,deltaPressure,deltaHeight, &
662 numFamily,familyList,numElement,elementList,lrandom, &
663 diagn_save,diagn_all,diagn_num,diagn_stnid,diagn_varno,diagn_unilev, &
664 diagn_nset,diagn_pressmin
665
666 ! set default values for namelist variables
667 nrandseed = 999
668 lrandom=.true.
669 deltaLat = 10.0d0
670 deltaLon = 10.0d0
671 deltaPressure = 10000.0d0
672 deltaHeight = 5000.0d0
673
674 !numFamily = ofl_numFamily
675 !familyList(:) = ofl_familyList(:)
676 numFamily = MPC_missingValue_INT
677 familyList(:) = ' '
678
679 numElement = MPC_missingValue_INT
680 elementList(:) = -1
681
682 diagn_save(:)=.false.
683 diagn_all(:)=.true.
684 diagn_pressmin(:)=10. ! 0.1 hPa
685 diagn_nset(:)=2
686 diagn_num(:)=0
687 diagn_stnid(:,:)='*********'
688 diagn_varno(:,:)=0
689 diagn_unilev(:,:)=.false.
690
691 nulnam = 0
692 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
693 read(nulnam,nml=namosd,iostat=ierr)
694 if(ierr /= 0) then
695 write(*,*) 'osd_setup: No valid namelist NAMOSD found, skipping some diagnostics'
696 nmlExists = .false.
697 ierr = fclos(nulnam)
698 return
699 else
700 nmlExists = .true.
701 endif
702 if(mmpi_myid == 0) write(*,nml=namosd)
703 ierr = fclos(nulnam)
704 if (numFamily /= MPC_missingValue_INT) then
705 call utl_abort('osd_setup: check NAMOSD namelist section: numFamily should be removed')
706 end if
707 numFamily = 0
708 do familyIndex = 1, ofl_numFamily
709 if (familyList(familyIndex) == ' ') exit
710 numFamily = numFamily + 1
711 end do
712 if (numElement /= MPC_missingValue_INT) then
713 call utl_abort('osd_setup: check NAMOSD namelist section: numElement should be removed')
714 end if
715 numElement = 0
716 do elementIndex = 1, ofl_numFamily
717 if (elementList(elementIndex) == -1) exit
718 numElement = numElement + 1
719 end do
720 do familyIndex = 1, numFamily
721 if ( .not. ofl_isFamilyTypeInList(familyList(familyIndex)) ) &
722 call utl_abort('osd_setup: Specified family '//familyList(familyIndex)//' was not recognized')
723 obsspace_diagn_filename(familyIndex) ='obsspace_diag_'//familyList(familyIndex)//'_'
724 if (diagn_num(familyIndex) > max_cfg_size) call utl_abort('osd_setup: Number exceeds allowed size of max_cfg_size')
725 end do
726
727 end subroutine osd_setup
728
729 !--------------------------------------------------------------------------
730 ! osd_update_obsfile
731 !--------------------------------------------------------------------------
732 subroutine osd_update_obsfile( obsSpaceData )
733 !
734 !:Purpose: Update of obs file(s) for content other
735 ! than OBS,OMA,OMP,OER,FGE,MRK in obsSpaceData
736 ! Content can be augmented as needed.
737 !
738 implicit none
739
740 ! Arguments:
741 type (struct_obs), intent(inout) :: obsSpaceData
742
743 ! If needed, add effective temperature values in CH family obs file
744 ! for total column measurements
745
746 if (obs_famExist(obsSpaceData,'CH')) call oopc_addEfftempObsfile()
747
748 end subroutine osd_update_obsfile
749
750 !-----------------------------------------------------------------------------------------
751 !------------------- Observation-space diagnostic functions and routines -----------------
752
753 !--------------------------------------------------------------------------
754 ! osd_obsPostProc
755 !--------------------------------------------------------------------------
756 subroutine osd_obsPostProc( obsSpaceData, columnTrlOnAnlIncLev, deltaLat, deltaLon, deltaPressure, anlm_mode )
757 !
758 !:Purpose: Interface for observation-space post-processing procedures.
759 !
760 !:Arguments:
761 ! :obsSpaceData: Obs space data structure
762 ! :obsfam: Target obs family (e.g. CH)
763 ! :codtyplist: Code type list asscoiated to obsfam.
764 ! :columnTrlOnAnlIncLev: Columns from analysis vertical coordinate in obs space (at obs location)
765 ! :date: YYYYMMDDHH
766 ! :deltaLat: Size of latitude bins for diagnostics (degrees)
767 ! :deltaLon: Size of longitude bins for diagnostics (degrees)
768 ! :deltaPressure: Size of vertical bins for diagnostics (Pascal)
769 ! :anlm_mode: Logical indicating if OmA (and Jo) diagnostics to be generated.
770 !
771 implicit none
772
773 ! Arguments:
774 type(struct_obs), intent(inout) :: obsSpaceData
775 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
776 real(8), intent(in) :: deltaLat
777 real(8), intent(in) :: deltaLon
778 real(8), intent(in) :: deltaPressure
779 logical, intent(in) :: anlm_mode
780
781 ! Locals:
782 integer, allocatable :: codtyplist(:)
783 integer :: jelm,ifam
784
785 if (mmpi_myid == 0) then
786 write(*,*)
787 write(*,*) "Enter osd_obsPostProc: Observation-space post-processing tasks for chemical constituents"
788 write(*,*)
789 end if
790
791 ! Generate and output cost function, OmP, and OmA diagnostics
792
793 if (obs_famExist(obsSpaceData,'CH')) then
794
795 ifam=0
796 do jelm=1,numFamily
797 if (familyList(jelm) == 'CH') then
798 ifam=jelm
799 exit
800 end if
801 end do
802 if (ifam == 0) then
803
804 write(*,*) 'osd_obsPostProc: Warning - No post-processing requested for CH family.'
805
806 else
807
808 ! Initialize oss_comboIDlist and add (stnid,varno) pairs from the namelist
809 ! Sets list of identifiers for observations to be processed in osd_obsDiagnostics within the CH family
810
811 call oss_comboIdlist(initialize_opt=.true., nset_opt=diagn_nset(ifam), all_combos_opt=diagn_all(ifam))
812 do jelm=1,diagn_num(ifam)
813 call oss_comboIdlist(stnid_add_opt=diagn_stnid(ifam,jelm), varno_add_opt=diagn_varno(ifam,jelm), unilev_add_opt=diagn_unilev(ifam,jelm))
814 end do
815
816 ! Diagnostics for retrievd chemical constituents (CH family)
817
818 allocate(codtyplist(2))
819
820 codtyplist(1)=codtyp_get_codtyp('CHEMREMOTE')
821 codtyplist(2)=codtyp_get_codtyp('CHEMINSITU')
822
823 call osd_obsDiagnostics(obsSpaceData,columnTrlOnAnlIncLev,'CH',codtyplist,trim(obsspace_diagn_filename(ifam)), &
824 diagn_save(ifam),deltaLat,deltaLon,deltaPressure,diagn_pressmin(ifam),anlm_mode)
825
826 deallocate(codtyplist)
827
828 end if
829
830 end if
831
832 ! Generate and output cost function, OmP, and OmA diagnostics for
833 ! channels/instruments of the TO family (when processed with accompanying
834 ! CH obs).
835
836 ! call osd_TO_obsDiagnostics(obsSpaceData,columnTrlOnAnlIncLev,date,deltaLat,deltaLon,deltaPressure,anlm_mode)
837
838 ! Apply any required obs file update
839
840 call osd_update_obsfile(obsSpaceData)
841
842 if (mmpi_myid == 0) then
843 write(*,*)
844 write(*,*) "Exit osd_obsPostProc"
845 write(*,*)
846 end if
847
848 end subroutine osd_obsPostProc
849
850 !--------------------------------------------------------------------------
851 ! osd_obsDiagnostics
852 !--------------------------------------------------------------------------
853 subroutine osd_obsDiagnostics( obsSpaceData, columnTrlOnAnlIncLev, obsfam, codtyplist, filename, save_diagn, &
854 deltaLat, deltaLon, deltaPressure, pressmin, anlm_mode )
855 !
856 !:Purpose: Calculates and prints observation-space diagnostics for chemical constituents
857 !
858 !:Arguments:
859 ! :obsSpaceData: Obs space data structure
860 ! :columnTrlOnAnlIncLev: Columns from analysis vertical coordinate in obs space (at obs location)
861 ! :obsfam: Obs family (e.g. 'CH'
862 ! :codtypelist: Code type list
863 ! :filename: Output file name
864 ! :save_diagn: Logical indicating gridded diagnostics are to be save
865 ! :date: YYYYMMDDHH
866 ! :deltaLat: Size of latitude bins for diagnostics (degrees)
867 ! :deltaLon: Size of longitude bins for diagnostics (degrees)
868 ! :deltaPressure: Size of vertical bins for diagnostics (Pascal)
869 ! :pressmin: bottom of top layer for diagnostics (in Pa).
870 ! :anlm_mode: Logical indicating if OmA diagnostics are to be generated.
871 !
872 !:Output: Content of ascii file with obs space diagnostics
873 !
874 !:Comments:
875 !
876 ! - Although Jo_analysis is already calculated in OBS_JOBS obsSpaceData and can be passed
877 ! to osd_obsspace_diagn_add, it is recalculated in osd_obsspace_diagn_add since OBS_JOBS
878 ! will be set to zero for diagnostic-only observations.
879 !
880 implicit none
881
882 ! Arguments:
883 type(struct_obs) , intent(inout) :: obsSpaceData
884 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
885 character(len=*) , intent(in) :: obsfam
886 character(len=*) , intent(in) :: filename
887 integer, intent(in) :: codtyplist(:)
888 real(8), intent(in) :: deltaLat
889 real(8), intent(in) :: deltaLon
890 real(8), intent(in) :: deltaPressure
891 real(8), intent(in) :: pressmin
892 logical, intent(in) :: anlm_mode
893 logical, intent(in) :: save_diagn
894
895 ! Locals:
896 type(struct_osd_diagn) :: obs_diagn
897 integer :: headerIndex,bodyIndex,vco,nlev_obs,ilev_obs,nlev_mod,ilev_mod
898 integer, parameter :: nmax=100
899 integer :: varno,varno_elemID(nmax)
900 integer :: elemID,num_elemID,nset,iass
901 character(len=9) :: stnid_elemID(nmax)
902 logical :: unilev_elemID(nmax),unilevel,diagn_only,assim_obs,status_hpht
903 character(len=256) :: label
904 real(8) :: lat,lon
905 character(len=12) :: stnid
906 real(8), allocatable :: lev(:), omp(:), oma(:), obs(:)
907 real(8), allocatable :: pres_mod(:), sigma_obs(:), sqrtHPHT(:)
908 logical, allocatable :: success(:)
909 integer, allocatable :: status(:)
910 real(8), pointer :: height_mod(:)
911
912 ! Get combination lists to group diagnostics by
913 call oss_get_comboIdlist(obsSpaceData,stnid_elemID,varno_elemID,unilev_elemID,num_elemID,nset)
914
915 if (num_elemID == 0) return
916
917 if (mmpi_myid == 0) then
918 write(*,*)
919 write(*,*) "osd_obsDiagnostics: Observation-space diagnostics for chemical constituents"
920 write(*,*)
921 end if
922
923 ! Read forecast error std. dev. at obs locations from obs files (saved in OBS_HPHT in obsSpaceData)
924 !status_hpht=.false.
925 if (anlm_mode) call osd_ReadSqrtHPHT(obsSpaceData,obsfam,codtyplist,status_hpht)
926
927 ! Allocate memory for diagnostic arrays in obs_diagn
928 call osd_obsspace_diagn_alloc(obs_diagn,deltaLat,deltaLon,deltaPressure,pressmin)
929
930 ! Loop over all pairs in *_elemID lists
931
932 do elemID=1,num_elemID
933
934 ! Initialize the diagnostic arrays
935 call osd_obsspace_diagn_init(obs_diagn)
936
937 call obs_set_current_header_list(obsSpaceData,obsfam)
938 HEADER: do
939
940 headerIndex = obs_getHeaderIndex(obsSpaceData)
941 if (headerIndex < 0) exit HEADER
942
943 ! Body info that we only need for first point in the profile
944 bodyIndex = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
945 vco = obs_bodyElem_i(obsSpaceData,OBS_VCO,bodyIndex)
946
947 if (vco /= 1 .and. vco /= 2 .and. vco /= 4 .and. vco /= 5) then
948 ! Vertical coordinate not handled
949 write(*,*) 'osd_obsDiagnostics: Currently unaccounted VCO = ',vco
950 cycle HEADER
951 end if
952
953 varno = obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
954 nlev_obs = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)
955 stnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
956
957 ! Identify max number of profile points in the profile (exclude BUFR_SCALE_EXPONENT elements)
958 call obs_set_current_body_list(obsSpaceData,headerIndex)
959 do
960 bodyIndex = obs_getBodyIndex(obsSpaceData)
961 if (bodyIndex < 0) exit
962 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == BUFR_SCALE_EXPONENT) then
963 nlev_obs = nlev_obs-1
964 else
965 varno=obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
966 end if
967 end do
968
969 ! Determine if this observation should be added to this group (as specified by nset)
970 if (.not.utl_stnid_equal(stnid_elemID(elemID),stnid)) cycle HEADER
971 if (nset >= 2.and.varno /= varno_elemID(elemID)) cycle HEADER
972 if (nset >= 3.and..not.(( nlev_obs == 1 .and. vco == 4 ).eqv.unilev_elemID(elemID))) cycle HEADER
973
974 ! Accumulate for this combo
975
976 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)*MPC_DEGREES_PER_RADIAN_R8
977 lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)*MPC_DEGREES_PER_RADIAN_R8
978
979 allocate(lev(nlev_obs), omp(nlev_obs), oma(nlev_obs), obs(nlev_obs))
980 allocate(sigma_obs(nlev_obs), success(nlev_obs), status(nlev_obs))
981 if (anlm_mode) allocate(sqrtHPHT(nlev_obs))
982
983 lev(:) = 0.0d0
984 omp(:) = 0.0d0
985 oma(:) = 0.0d0
986 obs(:) = 0.0d0
987 sigma_obs(:) = 0.0d0
988 if (anlm_mode) sqrtHPHT(:) = 0.0d0
989 status(:) = 0
990 assim_obs = .false.
991 ilev_obs = 0
992
993 call obs_set_current_body_list(obsSpaceData,headerIndex)
994 BODY: do
995
996 bodyIndex = obs_getBodyIndex(obsSpaceData)
997 if (bodyIndex < 0) exit BODY
998
999 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= varno) cycle BODY
1000
1001 ilev_obs = ilev_obs+1
1002
1003 iass = obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex)
1004
1005 ! Indicates if diagnostics are to be calculated but observation not assimilated
1006 diagn_only = oopc_diagnOnly(obsfam,stnid,varno,nlev_obs,obs_bodyElem_i(obsSpaceData,OBS_FLG,bodyIndex))
1007
1008 assim_obs = ((.not.diagn_only).and.anlm_mode) .or. assim_obs
1009
1010 if (iass == obs_assimilated) then
1011 status(ilev_obs) = 1
1012 else if (diagn_only) then
1013 status(ilev_obs) = 2
1014 end if
1015
1016 lev(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_PPP,bodyIndex)
1017
1018 ! Include in the sum assimilated data and diagnostic only data
1019 if (status(ilev_obs) > 0) then
1020
1021 omp(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_OMP,bodyIndex)
1022 obs(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_VAR,bodyIndex)
1023 sigma_obs(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_OER,bodyIndex)
1024 if (anlm_mode) then
1025 oma(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_OMA,bodyIndex)
1026 if (status_hpht) then
1027 sqrtHPHT(ilev_obs) = obs_bodyElem_r(obsSpaceData,OBS_HPHT,bodyIndex)
1028 if (sqrtHPHT(ilev_obs) < 1.D-30) then
1029 write(*,*) 'osd_obsDiagnostics: WARNING. sqrtHPHT not found for all obs'
1030 write(*,*) 'Will not be used in Desroziers-based diagnostics.'
1031 status_hpht=.false.
1032 end if
1033 end if
1034 end if
1035
1036 end if
1037
1038 end do BODY
1039
1040 ! Convert to pressure if needed and identify unilevel observations
1041 unilevel = .false.
1042 select case(vco)
1043 case(1)
1044 ! Height coordinate
1045
1046 nlev_mod = col_getNumLev(columnTrlOnAnlIncLev,'TH') ! number of model levels
1047 height_mod => col_getColumn(columnTrlOnAnlIncLev,headerIndex,'Z_T') ! geopotential
1048
1049 allocate(pres_mod(nlev_mod))
1050
1051 do ilev_mod=1,nlev_mod
1052 pres_mod(ilev_mod) = col_getPressure(columnTrlOnAnlIncLev,ilev_mod,headerIndex,'TH') ! model pressure
1053 end do
1054
1055 ! Convert altidudes to pressure
1056 success = status > 0
1057 lev = phf_convert_z_to_pressure(lev,height_mod,pres_mod,nlev_obs,nlev_mod,lat/MPC_DEGREES_PER_RADIAN_R8,success)
1058
1059 deallocate(pres_mod)
1060
1061 case(4,5)
1062 ! Uni-level observations
1063 unilevel = .true.
1064 end select
1065
1066 ! Add observation to diagnostic arrays
1067 if (anlm_mode) then
1068 if (status_hpht) then
1069 call osd_obsspace_diagn_add(obs_diagn, lat, lon, lev, &
1070 pressmin, omp, obs, sigma_obs, &
1071 nlev_obs, unilevel, assim_obs, status, &
1072 oma_opt=oma, sqrtHPHT_opt=sqrtHPHT)
1073 else
1074 call osd_obsspace_diagn_add(obs_diagn, lat, lon, lev, &
1075 pressmin, omp, obs, sigma_obs, &
1076 nlev_obs, unilevel, assim_obs, status,oma_opt=oma)
1077 end if
1078 else
1079 call osd_obsspace_diagn_add(obs_diagn, lat, lon, lev, &
1080 pressmin, omp, obs, sigma_obs, &
1081 nlev_obs, unilevel, assim_obs, status)
1082 end if
1083
1084 deallocate(lev,omp,oma,obs,sigma_obs,success,status)
1085 if (anlm_mode) deallocate(sqrtHPHT)
1086
1087 end do HEADER
1088
1089 ! Prepare output identification
1090 select case(nset)
1091 case(1)
1092 write(label,*) ">>> Statistics for STNID ",stnid_elemID(elemID)
1093 case(2)
1094 write(label,*) ">>> Statistics for STNID ",stnid_elemID(elemID)," and BUFR # ",varno_elemID(elemID)
1095 case(3)
1096 write(label,*) ">>> Statistics for STNID ",stnid_elemID(elemID)," and BUFR # ",varno_elemID(elemID)," and UNILEV = ",unilev_elemID(elemID)
1097 end select
1098
1099 ! Sum over different processors
1100 call osd_obsspace_diagn_MPIreduce(obs_diagn)
1101
1102 ! Output, and deallocate diagnostic arrays
1103 if (mmpi_myid == 0) call osd_obsspace_diagn_print(obs_diagn,filename, save_diagn, &
1104 'stats', pressmin, status_hpht, label_opt=label, openfile_opt=.true.)
1105
1106 end do
1107
1108 ! Output diagnostics summary (over all CH observations)
1109 if (mmpi_myid == 0) call osd_obsspace_diagn_print(obs_diagn,filename, save_diagn, &
1110 'summary', pressmin, status_hpht, openfile_opt=.true.)
1111
1112 ! Deallocate arrays in obs_diagn
1113 call osd_obsspace_diagn_dealloc(obs_diagn)
1114
1115 if (mmpi_myid == 0) then
1116 write(*,*)
1117 write(*,*) "End osd_obsDiagnostics"
1118 write(*,*)
1119 end if
1120
1121 end subroutine osd_obsDiagnostics
1122
1123 !--------------------------------------------------------------------------
1124 ! osd_ReadSqrtHPHT
1125 !--------------------------------------------------------------------------
1126 subroutine osd_ReadSqrtHPHT( obsSpaceData, obsfam, codtyplist, status_hpht )
1127 !
1128 !:Purpose: Read background error std. dev. at obs locations from the obs files and store
1129 ! under OBS_HPHT in obsSpaceData
1130 !
1131 !:Arguments:
1132 ! :obsSpaceData: Observation space data
1133 ! :obsfam: Obs family. e.g. 'CH'
1134 ! :codtyplist: Code type list associated to obsfam
1135 ! :status_hpht: logical indicating if successfully retrieved sqrtHPHT from obs file
1136 !
1137 implicit none
1138
1139 ! Arguments:
1140 logical, intent(out) :: status_hpht
1141 type(struct_obs), intent(inout) :: obsSpaceData
1142 character(len=*), intent(in) :: obsfam
1143 integer, intent(in) :: codtyplist(:)
1144
1145 ! Locals:
1146 integer :: bodyIndex,headerIndex,rln,nlv,kk
1147 integer :: stat,varno,icodtyp
1148 integer, parameter :: max_nlev=500
1149 integer, parameter :: ndim=1
1150 real(8) :: array(max_nlev)
1151 character(len=12), parameter :: stnid='************'
1152 type(struct_oss_obsdata) :: SqrtHPHT_struct
1153
1154 write(*,*) 'osd_ReadSqrtHPHT: STARTING'
1155
1156 ! Retrieve data from FGE blocks, i.e. sqrt(HPHT)
1157 ! (with ndim=1, bkstp=15 and block_type='DATA')
1158
1159 status_hpht=.true.
1160 SqrtHPHT_struct = obsf_obsSub_read(obsfam,stnid,-1,max_nlev,ndim,bkstp_opt=15,block_opt='DATA', &
1161 match_nlev_opt=.false.,codtyp_opt=codtyplist)
1162
1163 if (SqrtHPHT_struct%nrep == 0) then
1164 write(*,*) 'osd_ReadSqrtHPHT: WARNING. sqrtHPHT not found in obs file(s).'
1165 write(*,*) 'Will not be used in Desroziers-based diagnostics.'
1166 status_hpht=.false.
1167 call oss_obsdata_dealloc(SqrtHPHT_struct)
1168 return
1169 end if
1170
1171 ! Save in OBS_HPHT of obsSpaceData
1172
1173 call obs_set_current_header_list(obsSpaceData,obsfam)
1174 HEADER: do
1175
1176 headerIndex = obs_getHeaderIndex(obsSpaceData)
1177 if (headerIndex < 0) exit HEADER
1178
1179 icodtyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
1180 if (all(icodtyp /= codtyplist)) cycle HEADER
1181
1182 ! Search for corresponding HPHT profile/element
1183
1184 array(:)=0.0D0
1185 array=oss_obsdata_get_data1d(SqrtHPHT_struct, &
1186 obs_headElem_r(obsSpaceData,OBS_LON,headerIndex),obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex), &
1187 obs_headElem_i(obsSpaceData,OBS_DAT,headerIndex),obs_headElem_i(obsSpaceData,OBS_ETM,headerIndex), &
1188 obs_elem_c(obsSpaceData,'STID',headerIndex),stat_opt=stat)
1189
1190 if (stat == 0) then
1191
1192 ! Store OBS_HPHT profile
1193
1194 rln=obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
1195 nlv=obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)
1196 kk=0
1197 do bodyIndex = rln, nlv + rln -1
1198 varno=obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex)
1199 if ( varno > 10000 .and. varno < 16000 ) then
1200 kk=kk+1
1201 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,array(kk))
1202 else
1203 call obs_bodySet_r(obsSpaceData,OBS_HPHT,bodyIndex,0.0D0)
1204 end if
1205 end do
1206 end if
1207
1208 enddo HEADER
1209
1210 call oss_obsdata_dealloc(SqrtHPHT_struct)
1211
1212 write(*,*) 'osd_ReadSqrtHPHT: DONE'
1213
1214 end subroutine osd_ReadSqrtHPHT
1215
1216 !--------------------------------------------------------------------------
1217 ! osd_obsspace_diagn_alloc
1218 !--------------------------------------------------------------------------
1219 subroutine osd_obsspace_diagn_alloc( obs_diagn, deltaLat, deltaLon, deltaPressure, pressmin )
1220 !
1221 !:Purpose: Allocates diagnostic arrays in obs_diagn.
1222 !
1223 !:Arguments:
1224 ! :deltaLat: latitude bin size in degrees
1225 ! :deltaLon: longitutde in degrees
1226 ! :deltaPressure: pressures bin size in Pa (approximate)
1227 ! :pressmin: bottom of top layer for diagnostics (in Pa).
1228 !
1229 implicit none
1230
1231 ! Arguments:
1232 type(struct_osd_diagn), intent(inout) :: obs_diagn
1233 real(8) , intent(in) :: deltaLat
1234 real(8) , intent(in) :: deltaLon
1235 real(8) , intent(in) :: deltaPressure
1236 real(8) , intent(in) :: pressmin
1237
1238 ! Locals:
1239 integer :: nlev,nlat,nlon,nbin,nstat
1240
1241 obs_diagn%deltaLat = deltaLat
1242 obs_diagn%deltaLon = deltaLon
1243 obs_diagn%deltaLogPressure = deltaPressure/1.0d5 ! set constant delta ln(P) bin
1244
1245 nlat = floor(180.0d0/deltaLat)
1246 nlon = floor(360.0d0/deltaLon)
1247
1248 ! Add a last unequal size bin if remainder is larger than one degree
1249 if (180.0d0-nlat*deltaLat > 1.) nlat = nlat+1
1250 if (360.0d0-nlon*deltaLon > 1.) nlon = nlon+1
1251
1252 ! Set number of levels for a pressure coordinate in hPa to cover the range
1253 ! of 0.01*pressmin (in hPa) to 1000 hPa for layers 2 to nlev-1.
1254 ! Layer 1 covers all pressure levels <= 0.01*pressmin hPa = top of layer 2,
1255 ! The bottom of layer nlev-1 is provided by press_bins(nlev).
1256 nlev = 2 + nint(log(1.0d5/pressmin)/obs_diagn%deltaLogPressure) ! last index is for unilevel observations
1257 nbin = nlev*nlat*nlon
1258
1259 ! Two different statistics held, stat=1 for RMS and stat=2 for mean
1260 nstat = 2
1261
1262 obs_diagn%nlat = nlat
1263 obs_diagn%nlon = nlon
1264 obs_diagn%nlev = nlev
1265 obs_diagn%nbin = nbin
1266 obs_diagn%nstat = nstat
1267
1268 if (.not. allocated(obs_diagn%OmP_stats)) allocate(obs_diagn%OmP_stats(nlat,nlon,nlev,nstat))
1269 if (.not. allocated(obs_diagn%OmA_stats)) allocate(obs_diagn%OmA_stats(nlat,nlon,nlev,nstat))
1270 if (.not. allocated(obs_diagn%obs_stats)) allocate(obs_diagn%obs_stats(nlat,nlon,nlev,nstat))
1271 if (.not. allocated(obs_diagn%Jo_stats)) allocate(obs_diagn%Jo_stats(nlat,nlon,nlev,nstat*2+1))
1272 if (.not. allocated(obs_diagn%Jpa_stats)) allocate(obs_diagn%Jpa_stats(nlat,nlon,nlev))
1273 if (.not. allocated(obs_diagn%counts)) allocate(obs_diagn%counts(nlat,nlon,nlev))
1274 if (.not. allocated(obs_diagn%nstatus)) allocate(obs_diagn%nstatus(nlat,nlon,nlev,0:2))
1275 if (.not. allocated(obs_diagn%diagR_stats)) allocate(obs_diagn%diagR_stats(nlat,nlon,nlev,3))
1276 if (.not. allocated(obs_diagn%diagHPHT_stats)) allocate(obs_diagn%diagHPHT_stats(nlat,nlon,nlev,3))
1277
1278 end subroutine osd_obsspace_diagn_alloc
1279
1280 !--------------------------------------------------------------------------
1281 ! osd_obsspace_diagn_init
1282 !--------------------------------------------------------------------------
1283 subroutine osd_obsspace_diagn_init(obs_diagn)
1284 !
1285 !:Purpose: Initializes diagnostic arrays in obs_diagn.
1286 !
1287 implicit none
1288
1289 ! Arguments:
1290 type(struct_osd_diagn), intent(inout) :: obs_diagn
1291
1292 obs_diagn%OmP_stats(:,:,:,:) = 0.0d0
1293 obs_diagn%OmA_stats(:,:,:,:) = 0.0d0
1294 obs_diagn%obs_stats(:,:,:,:) = 0.0d0
1295 obs_diagn%Jo_stats(:,:,:,:) = 0.0d0
1296 obs_diagn%Jpa_stats(:,:,:) = 0.0d0
1297 obs_diagn%counts(:,:,:) = 0
1298 obs_diagn%nstatus(:,:,:,:) = 0
1299 obs_diagn%diagR_stats(:,:,:,:) = 0.0d0
1300 obs_diagn%diagHPHT_stats(:,:,:,:) = 0.0d0
1301
1302 obs_diagn%allow_print_summary = .false.
1303 obs_diagn%assim_mode = .false.
1304
1305 end subroutine osd_obsspace_diagn_init
1306
1307 !--------------------------------------------------------------------------
1308 ! osd_obsspace_diagn_dealloc
1309 !--------------------------------------------------------------------------
1310 subroutine osd_obsspace_diagn_dealloc(obs_diagn)
1311 !
1312 !:Purpose: Deallocates diagnostic arrays in obs_diagn.
1313 !
1314 implicit none
1315
1316 ! Arguments:
1317 type(struct_osd_diagn), intent(inout) :: obs_diagn
1318
1319 deallocate(obs_diagn%OmP_stats,obs_diagn%OmA_stats,obs_diagn%obs_stats)
1320 deallocate(obs_diagn%Jo_stats,obs_diagn%Jpa_stats,obs_diagn%counts,obs_diagn%nstatus)
1321 deallocate(obs_diagn%diagR_stats,obs_diagn%diagHPHT_stats)
1322
1323 end subroutine osd_obsspace_diagn_dealloc
1324
1325 !--------------------------------------------------------------------------
1326 ! osd_obsspace_diagn_add
1327 !--------------------------------------------------------------------------
1328 subroutine osd_obsspace_diagn_add( obs_diagn, lat, lon, pressure, pressmin, OmP, obs, sigma_obs, nlev_obs, &
1329 unilevel, assim_obs, status, OmA_opt, sqrtHPHT_opt)
1330 !
1331 !:Purpose: Adds an observation to the diagnostic arrays in obs_diagn.
1332 !
1333 !:Arguments:
1334 ! :lat: latitude in degrees
1335 ! :lon: longitutde in degrees
1336 ! :pressure: pressures of the profile (Pa)
1337 ! :pressmin: bottom of top layer for diagnostics (in Pa).
1338 ! :OmP: obs - background
1339 ! :OmA_opt: obs - analysis
1340 ! :obs: observations
1341 ! :Jo: cost function
1342 ! :sigma_obs: observation error standard deviation
1343 ! :sqrtHPHT_opt: forecast error standard deviation in obs space
1344 ! :assim_obs: indicates if the profile belongs to an assimilated data set
1345 ! :status: indicates status of the observations, with values denoting:
1346 !
1347 ! - 0 - observation has been rejected and not included in diagnostics
1348 ! - 1 - observation has been assimilated
1349 ! - 2 - observation has been used for diagnostics only (not assimilated)
1350 ! only observations with status=1,2 will be added to the statistic arrays
1351 ! :nlev_obs: number of observations in the profile
1352 ! :unilevel: if the observation does not have a defined height coordinate
1353 !
1354 implicit none
1355
1356 ! Arguments:
1357 type(struct_osd_diagn), intent(inout) :: obs_diagn
1358 real(8) , intent(in) :: lat
1359 real(8) , intent(in) :: lon
1360 integer , intent(in) :: nlev_obs
1361 integer , intent(in) :: status(nlev_obs)
1362 real(8) , intent(in) :: pressure(nlev_obs)
1363 real(8) , intent(in) :: OmP(nlev_obs)
1364 real(8) , intent(in) :: obs(nlev_obs)
1365 real(8) , intent(in) :: sigma_obs(nlev_obs)
1366 real(8) , intent(in) :: pressmin
1367 logical , intent(in) :: unilevel
1368 logical , intent(in) :: assim_obs
1369 real(8) , optional, intent(in) :: OmA_opt(nlev_obs)
1370 real(8) , optional, intent(in) :: sqrtHPHT_opt(nlev_obs)
1371
1372 ! Locals:
1373 integer :: ilat,ilon,ilev,ilev_obs
1374
1375 if (assim_obs) obs_diagn%assim_mode = .true.
1376
1377 ! Put in first/list bin if lat,lon,level lower/higher than diagnostic range
1378
1379 ilat = max(min(1 + floor((90.0d0+lat)/obs_diagn%deltaLat), obs_diagn%nlat), 1)
1380 ilon = max(min(1 + floor(lon/obs_diagn%deltaLon), obs_diagn%nlon), 1)
1381
1382 LEVELS: do ilev_obs=1,nlev_obs
1383
1384 if (unilevel) then
1385 ilev = obs_diagn%nlev ! put unilevel data in last level index
1386 else
1387 ilev = max(min(2 + floor(log(pressure(ilev_obs)/pressmin)/obs_diagn%deltaLogPressure), obs_diagn%nlev-1), 1)
1388 end if
1389
1390 obs_diagn%nstatus(ilat,ilon,ilev,status(ilev_obs)) = obs_diagn%nstatus(ilat,ilon,ilev,status(ilev_obs)) + 1
1391
1392 if (status(ilev_obs) == 0) cycle LEVELS ! skip adding of stats if the observation was rejected
1393
1394 obs_diagn%counts(ilat,ilon,ilev) = obs_diagn%counts(ilat,ilon,ilev) + 1
1395
1396 obs_diagn%OmP_stats(ilat,ilon,ilev,1) = obs_diagn%OmP_stats(ilat,ilon,ilev,1) + OmP(ilev_obs)**2
1397 obs_diagn%OmP_stats(ilat,ilon,ilev,2) = obs_diagn%OmP_stats(ilat,ilon,ilev,2) + OmP(ilev_obs)
1398
1399 obs_diagn%obs_stats(ilat,ilon,ilev,1) = obs_diagn%obs_stats(ilat,ilon,ilev,1) + obs(ilev_obs)**2
1400 obs_diagn%obs_stats(ilat,ilon,ilev,2) = obs_diagn%obs_stats(ilat,ilon,ilev,2) + obs(ilev_obs)
1401
1402 obs_diagn%Jo_stats(ilat,ilon,ilev,2) = obs_diagn%Jo_stats(ilat,ilon,ilev,2) + 0.5 * OmP(ilev_obs)**2 / sigma_obs(ilev_obs)**2
1403 if (status(ilev_obs) == 1) obs_diagn%Jo_stats(ilat,ilon,ilev,4) = obs_diagn%Jo_stats(ilat,ilon,ilev,4) + 0.5 * OmP(ilev_obs)**2 / sigma_obs(ilev_obs)**2
1404
1405 if (present(OmA_opt)) then
1406
1407 obs_diagn%OmA_stats(ilat,ilon,ilev,1) = obs_diagn%OmA_stats(ilat,ilon,ilev,1) + OmA_opt(ilev_obs)**2
1408 obs_diagn%OmA_stats(ilat,ilon,ilev,2) = obs_diagn%OmA_stats(ilat,ilon,ilev,2) + OmA_opt(ilev_obs)
1409 obs_diagn%Jo_stats(ilat,ilon,ilev,1) = obs_diagn%Jo_stats(ilat,ilon,ilev,1) + 0.5 * OmA_opt(ilev_obs)**2 / sigma_obs(ilev_obs)**2
1410
1411 if (status(ilev_obs) == 1) then
1412
1413 obs_diagn%diagR_stats(ilat,ilon,ilev,1) = obs_diagn%diagR_stats(ilat,ilon,ilev,1) &
1414 + OmP(ilev_obs)*OmA_opt(ilev_obs)/sigma_obs(ilev_obs)**2
1415 obs_diagn%diagR_stats(ilat,ilon,ilev,2) = obs_diagn%diagR_stats(ilat,ilon,ilev,2) &
1416 + OmP(ilev_obs)/sigma_obs(ilev_obs)
1417 obs_diagn%diagR_stats(ilat,ilon,ilev,3) = obs_diagn%diagR_stats(ilat,ilon,ilev,3) &
1418 + OmA_opt(ilev_obs)/sigma_obs(ilev_obs)
1419
1420 if (present(sqrtHPHT_opt)) then
1421 if (sqrtHPHT_opt(ilev_obs) > 0.0) then
1422 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,1) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,1) &
1423 + OmP(ilev_obs)*(OmP(ilev_obs)-OmA_opt(ilev_obs))/sqrtHPHT_opt(ilev_obs)**2
1424 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,2) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,2) &
1425 + OmP(ilev_obs)/sqrtHPHT_opt(ilev_obs)
1426 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,3) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,3) &
1427 + (OmP(ilev_obs)-OmA_opt(ilev_obs))/sqrtHPHT_opt(ilev_obs)
1428 ! Following two diagnostics do not account for spatial correlations (spatial correlations of HPHT)!
1429 ! As a consequence, Jb likely overestimated.
1430 obs_diagn%Jo_stats(ilat,ilon,ilev,5) = obs_diagn%Jo_stats(ilat,ilon,ilev,5) + 0.5 * (OmP(ilev_obs))**2 &
1431 / (sigma_obs(ilev_obs)**2 + sqrtHPHT_opt(ilev_obs)**2)
1432 obs_diagn%Jpa_stats(ilat,ilon,ilev) = obs_diagn%Jpa_stats(ilat,ilon,ilev) + 0.5 * (OmA_opt(ilev_obs)-OmP(ilev_obs))**2 &
1433 / sqrtHPHT_opt(ilev_obs)**2
1434 end if
1435 else
1436
1437 ! No division by sqrtHPHT
1438
1439 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,1) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,1) + OmP(ilev_obs)*(OmP(ilev_obs)-OmA_opt(ilev_obs))
1440 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,2) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,2) + OmP(ilev_obs)
1441 obs_diagn%diagHPHT_stats(ilat,ilon,ilev,3) = obs_diagn%diagHPHT_stats(ilat,ilon,ilev,3) + OmP(ilev_obs)-OmA_opt(ilev_obs)
1442 end if
1443
1444 end if
1445
1446 end if
1447
1448 end do LEVELS
1449
1450 end subroutine osd_obsspace_diagn_add
1451
1452 !--------------------------------------------------------------------------
1453 ! osd_obsspace_diagn_MPIreduce
1454 !--------------------------------------------------------------------------
1455 subroutine osd_obsspace_diagn_MPIreduce(obs_diagn)
1456 !
1457 !:Purpose: Performs a MPI allreduce on diagnostic arrays in obs_diagn.
1458 !
1459 implicit none
1460
1461 ! Arguments:
1462 type(struct_osd_diagn), intent(inout) :: obs_diagn
1463
1464 ! Locals:
1465 real(8), allocatable :: OmP_global(:,:,:,:), OmA_global(:,:,:,:), obs_global(:,:,:,:), Jo_global(:,:,:,:)
1466 real(8), allocatable :: Jpa_global(:,:,:), diagR_global(:,:,:,:), diagHPHT_global(:,:,:,:)
1467 integer, allocatable :: counts_global(:,:,:),nstatus_global(:,:,:,:)
1468 logical :: assim_global
1469 integer :: nlat,nlon,nlev,nstat,nbin,ierr
1470
1471 nlat = obs_diagn%nlat
1472 nlon = obs_diagn%nlon
1473 nlev = obs_diagn%nlev
1474 nstat = obs_diagn%nstat
1475 nbin = obs_diagn%nbin
1476
1477 ! Allocate memory for mpi global arrays
1478 allocate(OmP_global(nlat,nlon,nlev,nstat))
1479 allocate(OmA_global(nlat,nlon,nlev,nstat))
1480 allocate(obs_global(nlat,nlon,nlev,nstat))
1481 allocate(Jo_global(nlat,nlon,nlev,nstat*2+1))
1482 allocate(Jpa_global(nlat,nlon,nlev))
1483 allocate(diagR_global(nlat,nlon,nlev,3))
1484 allocate(diagHPHT_global(nlat,nlon,nlev,3))
1485 allocate(counts_global(nlat,nlon,nlev))
1486 allocate(nstatus_global(nlat,nlon,nlev,0:2))
1487
1488 ! Reduce from all mpi processes
1489 call rpn_comm_allreduce(obs_diagn%OmP_stats,OmP_global,nbin*nstat,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1490 call rpn_comm_allreduce(obs_diagn%OmA_stats,OmA_global,nbin*nstat,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1491 call rpn_comm_allreduce(obs_diagn%obs_stats,obs_global,nbin*nstat,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1492 call rpn_comm_allreduce(obs_diagn%Jo_stats,Jo_global,nbin*(nstat*2+1),"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1493 call rpn_comm_allreduce(obs_diagn%Jpa_stats,Jpa_global,nbin,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1494 call rpn_comm_allreduce(obs_diagn%diagR_stats,diagR_global,nbin*3,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1495 call rpn_comm_allreduce(obs_diagn%diagHPHT_stats,diagHPHT_global,nbin*3,"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
1496 call rpn_comm_allreduce(obs_diagn%counts,counts_global,nbin,"MPI_INTEGER","MPI_SUM","GRID",ierr)
1497 call rpn_comm_allreduce(obs_diagn%nstatus,nstatus_global,nbin*3,"MPI_INTEGER","MPI_SUM","GRID",ierr)
1498 call rpn_comm_allreduce(obs_diagn%assim_mode,assim_global,1,"MPI_LOGICAL","MPI_LOR","GRID",ierr)
1499
1500 ! save in struct_osd_diagn
1501 obs_diagn%OmP_stats = OmP_global
1502 obs_diagn%OmA_stats = OmA_global
1503 obs_diagn%obs_stats = obs_global
1504 obs_diagn%Jo_stats = Jo_global
1505 obs_diagn%Jpa_stats = Jpa_global
1506 obs_diagn%diagR_stats = diagR_global
1507 obs_diagn%diagHPHT_stats = diagHPHT_global
1508 obs_diagn%counts = counts_global
1509 obs_diagn%nstatus = nstatus_global
1510 obs_diagn%assim_mode = assim_global
1511
1512 deallocate(OmP_global,OmA_global,obs_global,Jo_global,Jpa_global,diagR_global,diagHPHT_global, &
1513 counts_global,nstatus_global)
1514
1515 end subroutine osd_obsspace_diagn_MPIreduce
1516
1517 !--------------------------------------------------------------------------
1518 ! osd_obsspace_diagn_print
1519 !--------------------------------------------------------------------------
1520 subroutine osd_obsspace_diagn_print(obs_diagn, filename, save_diagn, print_type, pressmin, status_hpht, label_opt, openfile_opt )
1521 !
1522 !:Purpose: Prints observation space diagnostics. If called with print_type = 'stats', the
1523 ! printed statistics will be added to the total diagnostic arrays.
1524 !
1525 !:Arguments:
1526 ! :filename: output file name
1527 ! :save_diagn: Logical indicating gridded diagnostics are to be save
1528 ! :print_type: Specifies which statistics to print, with possible values:
1529 !
1530 ! - 'stats' - prints statistics for the the arrays within obs_diagn
1531 ! - 'summary'- prints total statistics held in the saved variables
1532 ! within this subrouine
1533 ! :pressmin: min pressure level for output
1534 ! :label_opt: label to print (only relevant if print_type = 'stats')
1535 ! :openfile_opt: logical indicating if file filename is to be opened.
1536 ! :status_hpht: logical indicating if sqrtHPHT were available.
1537 !
1538 implicit none
1539
1540 ! Arguments:
1541 type(struct_osd_diagn) , intent(inout) :: obs_diagn
1542 character(len=*) , intent(in) :: print_type
1543 character(len=*) , intent(in) :: filename
1544 real(8) , intent(in) :: pressmin
1545 logical , intent(in) :: status_hpht
1546 logical , intent(in) :: save_diagn
1547 logical , optional, intent(in) :: openfile_opt
1548 character(len=256), optional, intent(in) :: label_opt
1549
1550 ! Locals:
1551 integer, external :: fnom, fclos
1552 real(8) :: Jo_a,Jo_b,Jpa_assim, Jo_a_assim,Jo_b_assim, Jo_p_assim
1553 real(8), save :: Jo_a_total=0.0d0, Jo_b_total=0.0d0, Jpa_total_assim=0.0d0
1554 real(8), save :: Jo_a_total_assim=0.0d0, Jo_b_total_assim=0.0d0, Jo_p_total_assim=0.0d0
1555 integer, save :: counts_total=0, counts_total_assim=0
1556 integer :: ierr,unit,icount,nlat,nlon,nlev,ilev,ilat,ilon,icount_assim
1557 integer, allocatable :: ncounts(:), ncounts_assim(:)
1558 real(8), allocatable :: press_bins(:)
1559 logical :: fileout_exist,multilevel,unilevel
1560
1561 select case(trim(print_type))
1562 case('stats','STATS')
1563
1564 ! Print observation-space statistics to listing file or to output file obspace_diag_filename.
1565
1566 ! Open and append to output file if requested
1567
1568 if (present(openfile_opt)) then
1569 if (openfile_opt) then
1570 inquire(file=filename, exist=fileout_exist)
1571 call utl_open_asciifile(filename,unit)
1572 else
1573 unit=6
1574 end if
1575 else
1576 unit=6
1577 end if
1578
1579 if (present(label_opt)) then
1580 write(unit,*)
1581 write(unit,*) trim(label_opt)
1582 end if
1583
1584 if (any(obs_diagn%counts > 0)) then
1585
1586 nlat = obs_diagn%nlat
1587 nlon = obs_diagn%nlon
1588 nlev = obs_diagn%nlev
1589
1590 allocate(ncounts(nlev), ncounts_assim(nlev), press_bins(nlev+1))
1591
1592 ! Pressure boundaries for each bin starting from a top layer with lower boundary of 0.01*pressmin (i=2) in hPa
1593 ! and extending to the surface.
1594
1595 press_bins(2:nlev-1) = (/ (0.01*pressmin*exp((ilev-2)*obs_diagn%deltaLogPressure), ilev=2,nlev-1) /)
1596 press_bins(1) = 0.0d0
1597 press_bins(nlev) = 1200.0d0 ! set to pressure larger than the largest expected surface pressure.
1598 press_bins(nlev+1) = 0.0d0 ! set to zero for unilevel
1599
1600 ! Total counts for each level
1601 ncounts = sum(sum(obs_diagn%counts,dim=1),dim=1)
1602 counts_total = counts_total + sum(ncounts)
1603 ncounts_assim = sum(sum(obs_diagn%nstatus(:,:,:,1),dim=1),dim=1)
1604 counts_total_assim = counts_total_assim + sum(ncounts_assim)
1605
1606 ! Indicates if any multilevel or unilevel observations exist
1607 multilevel = any(obs_diagn%nstatus(:,:,1:nlev-1,:) > 0)
1608 unilevel = any(obs_diagn%nstatus(:,:,nlev,:) > 0)
1609
1610 if (obs_diagn%assim_mode) then
1611 write(unit,*)
1612 write(unit,'(A)') " Elements for calc of obs and background error standard deviation scaling factors via"
1613 write(unit,'(A)') " the Desroziers approach can be found in the third and fourth blocks of statistics."
1614 write(unit,*)
1615 end if
1616
1617 write(unit,*)
1618 write(unit,'(A)') " Global statistics"
1619 write(unit,'(A)') " -----------------"
1620
1621 ! Multi-level data
1622
1623 if (multilevel) then
1624
1625 icount = sum(ncounts(1:nlev-1))
1626 icount_assim = sum(ncounts_assim(1:nlev-1))
1627
1628 if (icount > 0) then
1629 Jo_a = sum(obs_diagn%Jo_stats(:,:,1:nlev-1,1))
1630 Jo_b = sum(obs_diagn%Jo_stats(:,:,1:nlev-1,2))
1631 Jo_a_assim = sum(obs_diagn%Jo_stats(:,:,1:nlev-1,3))
1632 Jo_b_assim = sum(obs_diagn%Jo_stats(:,:,1:nlev-1,4))
1633 Jo_p_assim = sum(obs_diagn%Jo_stats(:,:,1:nlev-1,5))
1634 Jpa_assim = sum(obs_diagn%Jpa_stats(:,:,1:nlev-1))
1635 Jo_a_total = Jo_a_total + Jo_a
1636 Jo_b_total = Jo_b_total + Jo_b
1637 Jo_a_total_assim = Jo_a_total_assim + Jo_a_assim
1638 Jo_b_total_assim = Jo_b_total_assim + Jo_b_assim
1639 Jo_p_total_assim = Jo_p_total_assim + Jo_p_assim
1640 Jpa_total_assim = Jpa_total_assim + Jpa_assim
1641 else
1642 Jo_a = 0.0d0
1643 Jo_a_assim = 0.0d0
1644 Jo_b_assim = 0.0d0
1645 Jo_p_assim = 0.0d0
1646 Jpa_assim = 0.0d0
1647 end if
1648
1649 obs_diagn%allow_print_summary = .true.
1650
1651 call print_J(unit,"Multi-level data:",Jo_a,Jo_b,icount,Jo_a_assim,Jo_b_assim,Jpa_assim,Jo_p_assim,icount_assim)
1652
1653 call print_stats(unit,obs_diagn,press_bins,1,nlat,1,nlon,1,nlev-1)
1654 if (obs_diagn%assim_mode) call print_Desroziers(unit,obs_diagn,press_bins,1,nlat,1,nlon,1,nlev-1,status_hpht)
1655
1656 end if
1657
1658 ! Uni-level data
1659
1660 if (unilevel) then
1661
1662 if (ncounts(nlev) > 0) then
1663 Jo_a = sum(obs_diagn%Jo_stats(:,:,nlev,1))
1664 Jo_b = sum(obs_diagn%Jo_stats(:,:,nlev,2))
1665 Jo_a_assim = sum(obs_diagn%Jo_stats(:,:,nlev,3))
1666 Jo_b_assim = sum(obs_diagn%Jo_stats(:,:,nlev,4))
1667 Jo_p_assim = sum(obs_diagn%Jo_stats(:,:,nlev,5))
1668 Jpa_assim = sum(obs_diagn%Jpa_stats(:,:,nlev))
1669 Jo_a_total = Jo_a_total + Jo_a
1670 Jo_b_total = Jo_b_total + Jo_b
1671 Jo_a_total_assim = Jo_a_total_assim + Jo_a_assim
1672 Jo_b_total_assim = Jo_b_total_assim + Jo_b_assim
1673 Jo_p_total_assim = Jo_p_total_assim + Jo_p_assim
1674 Jpa_total_assim = Jpa_total_assim + Jpa_assim
1675 else
1676 Jo_a = 0.0d0
1677 Jo_b = 0.0d0
1678 Jo_a_assim = 0.0d0
1679 Jo_b_assim = 0.0d0
1680 Jo_p_assim = 0.0d0
1681 Jpa_assim = 0.0d0
1682 end if
1683
1684 obs_diagn%allow_print_summary = .true.
1685
1686 call print_J(unit,"Uni-level data:",Jo_a,Jo_b,ncounts(nlev),Jo_a_assim,Jo_b_assim,Jpa_assim,Jo_p_assim,ncounts_assim(nlev))
1687
1688 call print_stats(unit,obs_diagn,press_bins,1,nlat,1,nlon,nlev,nlev)
1689 if (obs_diagn%assim_mode) call print_Desroziers(unit,obs_diagn,press_bins,1,nlat,1,nlon,nlev,nlev,status_hpht)
1690
1691 end if
1692
1693 deallocate(ncounts,ncounts_assim)
1694
1695 ! Output lat,lon dependent averages to file if obsspace_diagn_filename is provided
1696 if (present(openfile_opt)) then
1697 if (openfile_opt .and. save_diagn .and. (nlat > 1 .or. nlon > 1)) then
1698
1699 write(unit,*)
1700 write(unit,'(A)') " Lat-lon gridded statistics"
1701 write(unit,'(A)') " --------------------------"
1702 write(unit,*)
1703 write(unit,'(2X,3(A,I4))') "nlat = ",nlat," , nlon = ",nlon," , nlev = ",nlev
1704 write(unit,*)
1705
1706 do ilat=1,nlat
1707 do ilon=1,nlon
1708 write(unit,'(2X,2(A,I6),3X,2(F8.1,A,F8.1))') "ilat = ",ilat," , ilon = ",ilon, &
1709 (ilat-1.)*obs_diagn%deltaLat-90.," < lat < ",ilat*obs_diagn%deltaLat-90., &
1710 (ilon-1.)*obs_diagn%deltaLon," < lon < ",ilon*obs_diagn%deltaLon
1711 if (any(obs_diagn%nstatus(ilat,ilon,1:nlev-1,:) > 0)) then
1712 write(unit,*)
1713 write(unit,'(A)') " Multi-level data:"
1714 write(unit,*)
1715 call print_stats(unit,obs_diagn,press_bins,ilat,ilat,ilon,ilon,1,nlev-1)
1716 if (obs_diagn%assim_mode) call print_Desroziers(unit,obs_diagn,press_bins,ilat,ilat,ilon,ilon,1,nlev-1,status_hpht)
1717 else if (multilevel) then
1718 write(unit,*)
1719 write(unit,'(A)') " No multi-level data."
1720 write(unit,*)
1721 end if
1722 if (any(obs_diagn%nstatus(ilat,ilon,nlev,:) > 0)) then
1723 write(unit,*)
1724 write(unit,'(A)') " Uni-level data:"
1725 write(unit,*)
1726 call print_stats(unit,obs_diagn,press_bins,ilat,ilat,ilon,ilon,nlev,nlev)
1727 if (obs_diagn%assim_mode) call print_Desroziers(unit,obs_diagn,press_bins,ilat,ilat,ilon,ilon,nlev,nlev,status_hpht)
1728 else if (unilevel) then
1729 write(unit,*)
1730 write(unit,'(A)') " No uni-level data."
1731 write(unit,*)
1732 end if
1733 end do
1734 end do
1735
1736 end if
1737 end if
1738
1739 deallocate(press_bins)
1740
1741 else
1742 write(unit,*)
1743 write(unit,*) "No data found for this combination."
1744 write(unit,*)
1745 end if
1746
1747 if (present(openfile_opt)) then
1748 if (openfile_opt) ierr=fclos(unit)
1749 end if
1750
1751 case('summary','SUMMARY')
1752
1753 if (.not.obs_diagn%allow_print_summary) then
1754 write(*,*) "osd_obsspace_diagn_print: allow_print_summary is set to false, no summary will be printed."
1755 return
1756 end if
1757
1758 if (present(openfile_opt)) then
1759 if (openfile_opt) then
1760 inquire(file=filename, exist=fileout_exist)
1761 call utl_open_asciifile(filename,unit)
1762 else
1763 unit=6
1764 end if
1765 else
1766 unit=6
1767 end if
1768
1769 call print_J(unit,"Total cost function for CH observations:",Jo_a_total,Jo_b_total,counts_total, &
1770 Jo_a_total_assim,Jo_b_total_assim,Jpa_total_assim,Jo_a_total_assim,counts_total_assim)
1771
1772 if (present(openfile_opt)) then
1773 if (openfile_opt) ierr=fclos(unit)
1774 end if
1775
1776 case default
1777 call utl_abort("osd_obsspace_diagn_print: Invalid print_type select of " // trim(print_type) )
1778 end select
1779
1780
1781 contains
1782
1783 !--------------------------------------------------------------------------
1784 ! print_J
1785 !--------------------------------------------------------------------------
1786 subroutine print_J( unit, title, Jo_analysis, Jo_backgrnd, nobs, Jo_anal_assim, Jo_bgck_assim, Jpa_assim, Jo_p_assim, nobs_assim )
1787 !
1788 implicit none
1789
1790 ! Arguments:
1791 integer , intent(in) :: unit
1792 integer , intent(in) :: nobs
1793 integer , intent(in) :: nobs_assim
1794 character(len=*), intent(in) :: title
1795 real(8) , intent(in) :: Jo_analysis
1796 real(8) , intent(in) :: Jo_backgrnd
1797 real(8) , intent(in) :: Jo_anal_assim
1798 real(8) , intent(in) :: Jo_bgck_assim
1799 real(8) , intent(in) :: Jpa_assim
1800 real(8) , intent(in) :: Jo_p_assim
1801
1802 ! Locals:
1803 real(8) :: Jo_analysis_norm,Jo_backgrnd_norm,Jt_assim
1804 real(8) :: Jpa_norm_assim,Jt_norm_assim,Jo_p_norm_assim
1805 character(len=100) :: fmt
1806
1807 if (Jo_analysis < 1.0d6 .and. Jo_backgrnd < 1.0d6) then
1808 fmt = '(A,F24.8,A,F24.8)'
1809 else
1810 fmt = '(A,ES24.8,A,ES24.8)'
1811 end if
1812
1813 if (nobs > 0) then
1814 Jo_analysis_norm = 2.*Jo_analysis/nobs
1815 Jo_backgrnd_norm = 2.*Jo_backgrnd/nobs
1816 else
1817 Jo_analysis_norm = 0.0d0
1818 Jo_backgrnd_norm = 0.0d0
1819 end if
1820
1821 if ((nobs > 0 .and. nobs_assim > 0) .or. nobs_assim == 0) then
1822 write(unit,*)
1823 write(unit,'(A)') " " // title // " totals "
1824 write(unit,*)
1825 write(unit,fmt) " Jo(x_analysis) = ",Jo_analysis," , 2*Jo(x_analysis)/N = ",Jo_analysis_norm
1826 write(unit,fmt) " Jo(x_backgrnd) = ",Jo_backgrnd," , 2*Jo(x_backgrnd)/N = ",Jo_backgrnd_norm
1827 write(unit,'(A,I24)') " N = ",nobs
1828 write(unit,*)
1829 end if
1830
1831 if (nobs_assim > 0) then
1832 Jo_analysis_norm = 2.*Jo_anal_assim/nobs_assim
1833 Jo_backgrnd_norm = 2.*Jo_bgck_assim/nobs_assim
1834 Jpa_norm_assim = 2.*Jpa_assim/nobs_assim
1835 Jt_norm_assim = Jpa_norm_assim+Jo_analysis_norm
1836 Jo_p_norm_assim = 2.*Jo_p_assim/nobs_assim
1837
1838 write(unit,*)
1839 write(unit,'(A)') " " // title // " assimilated data "
1840 write(unit,*)
1841 write(unit,fmt) " J(A-P)+Jo(x_a) = ",Jt_assim, " , 2*[J(A-P)+Jo(x_a)]/N = ",Jt_norm_assim
1842 write(unit,fmt) " J(A-P) = ",Jpa_assim, " , 2*J(A-P)/N = ",Jpa_norm_assim
1843 write(unit,fmt) " Jo(x_analysis) = ",Jo_anal_assim," , 2*Jo(x_analysis)/N = ",Jo_analysis_norm
1844 write(unit,fmt) " Jo(x_backgrnd) = ",Jo_bgck_assim," , 2*Jo(x_backgrnd)/N = ",Jo_backgrnd_norm
1845 write(unit,fmt) " J(O-P) = ",Jo_p_assim, " , 2*J(O-P)/N = ",Jo_p_norm_assim
1846 write(unit,'(A,I24)') " N = ",nobs_assim
1847 write(unit,*)
1848 end if
1849
1850 end subroutine print_J
1851
1852 !--------------------------------------------------------------------------
1853 ! print_stats
1854 !--------------------------------------------------------------------------
1855 subroutine print_stats( unit, obs_diagn, pressure, ilat_start, ilat_end, ilon_start, ilon_end, ilev_start, ilev_end )
1856 implicit none
1857
1858 ! Arguments:
1859 type(struct_osd_diagn), intent(in) :: obs_diagn
1860 real(8) , intent(in) :: pressure(obs_diagn%nlev+1)
1861 integer , intent(in) :: unit
1862 integer , intent(in) :: ilat_start
1863 integer , intent(in) :: ilat_end
1864 integer , intent(in) :: ilon_start
1865 integer , intent(in) :: ilon_end
1866 integer , intent(in) :: ilev_start
1867 integer , intent(in) :: ilev_end
1868
1869 ! Locals:
1870 integer :: ilev,level,counts(obs_diagn%nlev),N_assim,N_diagn,N_rej
1871 real(8) :: pres1,pres2,jo_a,jo_b,jo_a_norm,jo_b_norm,obs_sum,obs_mean,obs_std,OmP_mean,OmP_rms,OmA_mean,OmA_rms
1872 logical :: skip(obs_diagn%nlev)
1873
1874 skip(:) = .false.
1875
1876 write(unit,'(A)') " Layer Pressure (hPa) Counts (N) N_assim N_diagn N_rej Jo(O-A) 2*Jo(O-A)/N Jo(O-P) 2*Jo(O-P)/N"
1877 write(unit,'(A)') " ----- -------------- ---------- ------- ------- ----- ------- ----------- ------- -----------"
1878
1879 do ilev = ilev_start, ilev_end
1880
1881 counts(ilev) = sum(obs_diagn%counts(ilat_start:ilat_end,ilon_start:ilon_end,ilev))
1882
1883 N_rej = sum(obs_diagn%nstatus(ilat_start:ilat_end,ilon_start:ilon_end,ilev,0))
1884 N_assim = sum(obs_diagn%nstatus(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))
1885 N_diagn = sum(obs_diagn%nstatus(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2))
1886
1887 skip(ilev) = counts(ilev) == 0 .and. N_rej == 0
1888
1889 if (skip(ilev)) cycle
1890
1891 if (ilev < nlev) then
1892 pres1 = pressure(ilev)
1893 pres2 = pressure(ilev+1)
1894 level = ilev
1895 else
1896 pres1 = 0.0d0
1897 pres2 = 0.0d0
1898 level = 0
1899 end if
1900
1901 jo_a = sum(obs_diagn%Jo_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))
1902 jo_b = sum(obs_diagn%Jo_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2))
1903
1904 if ( counts(ilev) == 0 ) then
1905 jo_a_norm = 0.0d0
1906 jo_b_norm = 0.0d0
1907 else
1908 jo_a_norm = 2.*jo_a/counts(ilev)
1909 jo_b_norm = 2.*jo_b/counts(ilev)
1910 end if
1911
1912 write(unit,'(2X, I3, 2(2X, F11.4), 4(2X,I9), 4(2X, ES11.4))') &
1913 level,pres1,pres2,counts(ilev),N_assim,N_diagn,N_rej,jo_a,jo_a_norm,jo_b,jo_b_norm
1914
1915 end do
1916
1917 write(unit,*)
1918 ! Division by <O> removed to prevent unlikely division by zero (or <= zero)
1919 ! write(unit,'(A)') " Layer Pressure (hPa) Counts obs (mean) obs (std) rms(O-P)/<O> <O-P>/<O> rms(O-A)/<O> <O-A>/<O>"
1920 ! write(unit,'(A)') " ----- -------------- ------ ---------- --------- ------------ --------- ------------ ---------"
1921 write(unit,'(A)') " Layer Pressure (hPa) Counts obs (mean) obs (std) rms(O-P) <O-P> rms(O-A) <O-A>"
1922 write(unit,'(A)') " ----- -------------- ------ ---------- --------- -------- ----- -------- -----"
1923
1924
1925 do ilev=ilev_start,ilev_end
1926
1927 if (skip(ilev)) cycle
1928
1929 if (ilev < nlev) then
1930 pres1 = pressure(ilev)
1931 pres2 = pressure(ilev+1)
1932 level = ilev
1933 else
1934 pres1 = 0.0d0
1935 pres2 = 0.0d0
1936 level = 0
1937 end if
1938
1939 if (counts(ilev) == 0) then
1940 obs_mean = 0.0d0
1941 obs_std = 0.0d0
1942 OmP_rms = 0.0d0
1943 OmP_mean = 0.0d0
1944 OmA_rms = 0.0d0
1945 OmA_mean = 0.0d0
1946
1947 write(unit,'(2X, I3, 2(2X, F11.4), 2X, I6, 6(2X, ES11.4))') &
1948 level,pres1,pres2,counts(ilev),obs_mean,obs_std,OmP_rms,OmP_mean,OmA_rms,OmA_mean
1949 else
1950 obs_sum = sum(obs_diagn%obs_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2))
1951 obs_mean = obs_sum / counts(ilev)
1952 obs_std = sqrt(max(0.0D0, sum(obs_diagn%obs_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))/counts(ilev) - (obs_sum/counts(ilev))**2 ))
1953
1954 OmP_rms = sqrt(max(0.0D0, sum(obs_diagn%OmP_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1)) / counts(ilev) ))
1955 OmP_mean = sum(obs_diagn%OmP_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2)) / obs_sum
1956
1957 OmA_rms = sqrt(max(0.0D0, sum(obs_diagn%OmA_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1)) / counts(ilev) ))
1958 OmA_mean = sum(obs_diagn%OmA_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2)) / counts(ilev)
1959
1960 ! write(unit,'(2X, I3, 2(2X, F11.4), 6(2X, ES11.4))') &
1961 ! level,pres1,pres2,obs_mean,obs_std,OmP_rms/obs_mean,OmP_mean/obs_sum,OmA_rms/obs_mean,OmA_mean/obs_sum
1962
1963 write(unit,'(2X, I3, 2(2X, F11.4), 2X, I6, 6(2X, ES11.4))') &
1964 level,pres1,pres2,counts(ilev),obs_mean,obs_std,OmP_rms,OmP_mean,OmA_rms,OmA_mean
1965
1966 end if
1967
1968 end do
1969
1970 write(unit,*)
1971
1972 end subroutine print_stats
1973
1974 !--------------------------------------------------------------------------
1975 ! print_Desroziers
1976 !--------------------------------------------------------------------------
1977 subroutine print_Desroziers( unit, obs_diagn, pressure, ilat_start, ilat_end, ilon_start, ilon_end, ilev_start, ilev_end, status_hpht )
1978 !
1979 !:Purpose: Prints elements contributing to the calc of scaling factors for observation and background
1980 ! error std. dev. based on the Desroziers approach.
1981 !
1982 implicit none
1983
1984 ! Arguments:
1985 type(struct_osd_diagn), intent(in) :: obs_diagn
1986 real(8) , intent(in) :: pressure(obs_diagn%nlev+1)
1987 integer , intent(in) :: unit
1988 integer , intent(in) :: ilat_start
1989 integer , intent(in) :: ilat_end
1990 integer , intent(in) :: ilon_start
1991 integer , intent(in) :: ilon_end
1992 integer , intent(in) :: ilev_start
1993 integer , intent(in) :: ilev_end
1994 logical , intent(in) :: status_hpht
1995
1996 ! Locals:
1997 integer :: ilev,level,N_assim(obs_diagn%nlev)
1998 real(8) :: pres1,pres2,sum_prod,sum_OmP,sum_OmA,sum_AmP,scaling
1999
2000 N_assim(:) = 0
2001
2002 write(unit,'(A)') " Layer Pressure (hPa) N_assim sum[(O-P)(O-A)/var(O)] sum[(O-P)/sig(O)] sum[(O-A)/sig(O)] mean[(O-P)(O-A)/var(O)]-mean[(O-P)/sig(O)]*mean[(O-A)/sig(O)]"
2003 write(unit,'(A)') " ----- -------------- ------- ---------------------- ----------------- ----------------- -------------------------------------------------------------"
2004
2005 do ilev=ilev_start,ilev_end
2006
2007 N_assim(ilev) = sum(obs_diagn%nstatus(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))
2008
2009 if (N_assim(ilev) == 0) cycle
2010
2011 if (ilev < nlev) then
2012 pres1 = pressure(ilev)
2013 pres2 = pressure(ilev+1)
2014 level = ilev
2015 else
2016 pres1 = 0.0d0
2017 pres2 = 0.0d0
2018 level = 0
2019 end if
2020
2021 sum_prod = sum(obs_diagn%diagR_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))
2022 sum_OmP = sum(obs_diagn%diagR_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2))
2023 sum_OmA = sum(obs_diagn%diagR_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,3))
2024
2025 scaling = sum_prod/N_assim(ilev) - sum_OmP*sum_OmA/N_assim(ilev)**2
2026
2027 write(unit,'(2X, I3, 2(2X, F11.4), 2X, I9, 11X, ES11.4, 2(8X,ES11.4), 10X, ES11.4)') level,pres1,pres2,N_assim(ilev),sum_prod,sum_OmP,sum_OmA,scaling
2028
2029 end do
2030
2031 write(unit,*)
2032 if (status_hpht) then
2033 write(unit,'(A)') " Layer Pressure (hPa) N_assim sum[(O-P)(A-P)/var(P)] sum[(O-P)/sig(P)] sum[(A-P)/sig(P)] mean[(O-P)(A-P)/var(P)]-mean[(O-P)/sig(P)]*mean[(A-P)/sig(P)]"
2034 else
2035 write(unit,'(A)') " Layer Pressure (hPa) N_assim sum[(O-P)(A-P)] sum[O-P] sum[A-P] mean[(O-P)(A-P)]-mean[O-P]*mean[A-P]"
2036 end if
2037 write(unit,'(A)') " ----- -------------- ------- ---------------------- ----------------- ----------------- -------------------------------------------------------------"
2038
2039 do ilev=ilev_start,ilev_end
2040
2041 if (N_assim(ilev) == 0) cycle
2042
2043 if (ilev < nlev) then
2044 pres1 = pressure(ilev)
2045 pres2 = pressure(ilev+1)
2046 level = ilev
2047 else
2048 pres1 = 0.0d0
2049 pres2 = 0.0d0
2050 level = 0
2051 end if
2052
2053 sum_prod = sum(obs_diagn%diagHPHT_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,1))
2054 sum_OmP = sum(obs_diagn%diagHPHT_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,2))
2055 sum_AmP = sum(obs_diagn%diagHPHT_stats(ilat_start:ilat_end,ilon_start:ilon_end,ilev,3))
2056
2057 scaling = sum_prod/N_assim(ilev) - sum_OmP*sum_AmP/N_assim(ilev)**2
2058
2059 write(unit,'(2X, I3, 2(2X, F11.4), 2X, I9, 11X, ES11.4, 2(8X,ES11.4), 10X, ES11.4)') level,pres1,pres2,N_assim(ilev),sum_prod,sum_OmP,sum_AmP,scaling
2060
2061 end do
2062
2063 write(unit,*)
2064
2065 end subroutine print_Desroziers
2066
2067 end subroutine osd_obsspace_diagn_print
2068
2069end module obsSpaceDiag_mod
