1module obsSpaceErrorStdDev_mod
2 ! MODULE obsSpaceErrorStdDev_mod (prefix='ose' category='1. High-level functionality')
3 !
4 !:Purpose: Contains subroutines for computing background-error and OmP-error
5 ! standard deviations in observation space
6 !
7 use midasMpi_mod
8 use obsSpaceData_mod
9 use columnData_mod
10 use bufr_mod
11 use utilities_mod
12 use earthConstants_mod
13 use MathPhysConstants_mod
14 use stateToColumn_mod
15 use gridStateVector_mod
16 use verticalCoord_mod
17 use horizontalCoord_mod
18 use bmatrixhi_mod
19 use obsOperators_mod
20 use gps_mod
21 use codtyp_mod
22 use bCovarSetupChem_mod
23 use timeCoord_mod
24 use obsTimeInterp_mod
25 use bMatrixEnsemble_mod
26 use varNameList_mod
27 use obsOperatorsChem_mod
28 use obsFamilyList_mod
29 use calcHeightAndPressure_mod
30
31 implicit none
32 private
33
34 ! public procedures
35 public :: ose_computeStddev
36
37 ! module structures
38 ! -----------------
39
40 type :: struct_OmPStdDevCH
41 !
42 ! Structure containing information retrieved from auxiliary file for holding
43 ! OmP std dev information
44 !
45 ! Variable Description
46 ! -------- -----------
47 ! n_stnid Number of sub-families (identified via STNIDs)
48 ! stnids Sub-families (STNIDs; * are wild cards)
49 ! element BUFR element in data block
50 ! source 0: Set entirely from the auxiliary file being read
51 ! as a function latitude, vertical level, and month.
52 ! 1: Values to be calculated from OmP differences
53 ! (requires a large data set size for each assim)
54 ! using read in latitudes and vertical levels
55 ! 2: Values to be calculated from OBS_OER and estimated
56 ! OBS_HPHT (default in absence of stnid and element)
57 ! std_type Index of std dev being read for source=0 or calculated for source=1
58 ! 0: absolute value
59 ! 1: fraction of measurement value (% / 100)
60 ! ibegin Position index of start of data for given
61 ! sub-family in the arrays OmPstd,levels,lat,month
62 ! n_lvl Number of vertical levels (max number when source=2)
63 ! levels Vertical levels (in coordinate of sub-family data)
64 ! n_lat Number of latitudes
65 ! lat Latitudes (degrees; ordered in increasing size)
66 ! n_month Number of months (requires n_lat>1 to be >1)
67 ! month Months numbered between 1 and 12.
68 ! std OmP error std dev (or fraction)
69
70 integer :: n_stnid
71 character(len=12), allocatable :: stnids(:)
72 integer, allocatable :: element(:),n_lat(:),n_month(:),std_type(:)
73 integer, allocatable :: source(:),ibegin(:),n_lvl(:)
74 real(8), allocatable :: std(:)
75 real(8), allocatable :: levels(:),lat(:),month(:)
76
77 end type struct_OmPStdDevCH
78
79 type(struct_OmPStdDevCH) :: OmPstdCH
80 type(struct_hco), pointer :: hco_anl => null()
81 real*8 , allocatable :: ose_vRO_Jacobian(:,:,:)
82 logical, allocatable :: ose_vRO_lJac(:)
83
84 contains
85
86 !--------------------------------------------------------------------------
87 ! ose_computeStddev
88 !--------------------------------------------------------------------------
89 subroutine ose_computeStddev(columnTrlOnAnlIncLev,hco_anl_in,obsSpaceData)
90 !
91 !:Purpose: To set OmP-error std dev when possible. Otherwise
92 ! compute background-error stddev in observation space to
93 ! estimate OmP-error std dev.
94 !
95 implicit none
96
97 ! Arguments:
98 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev ! Background columns on anl levels, obs horiz locations
99 type(struct_hco), pointer, intent(in) :: hco_anl_in
100 type(struct_obs), intent(inout) :: obsSpaceData ! Observation-related data
101
102 ! Locals:
103 real(8) :: HBHT_static, HBHT_ensemble, HBHT_hybrid
104 logical :: staticHBHT = .false.
105 logical :: staticOMPE = .false.
106 logical :: staticHBHT_ch = .false.
107 logical :: staticOMPE_ch = .false.
108 logical :: ensemble = .false.
109 integer :: index_body
110 integer :: fnom, fclos, ierr, nulnam
111
112 ! Namelist variables:
113 character(len=12) :: hybrid_mode ! can be 'WEIGHTED_SUM' or 'MAX_VALUE'
114 NAMELIST /NAMHBHT/hybrid_mode
115
116 ! set the module variable hco_anl
117 hco_anl => hco_anl_in
118
119 ! return if there is no observation.
120 if ( obs_numheader(obsSpaceData) == 0 ) return
121
122 !- 1. Compute HBHT (sigma_b in observation space) or OmP error std dev
123
124 !- 1.1 Compute error std dev from static statistics
125 ! OmP error std dev (OBS_OMPE) when possible and required.
126 ! Otherwise, compute HBHT
127 ! obsSpaceData - INOUT ( OmP error std dev outputted in OBS_OMPE or/and
128 ! sqrt(diag(H*B*H^T)) with B_static_chm outputted in OBS_HPHT )
129
130 call ose_setStaticErrorStddev( columnTrlOnAnlIncLev, obsSpaceData, staticHBHT, staticHBHT_ch, staticOMPE_ch )
131
132 !- 1.2 HBHT from the Bens
133 ! obsSpaceData - INOUT (HBensHT std. dev. outputted in OBS_WORK)
134
135 call ose_compute_hbht_ensemble( columnTrlOnAnlIncLev, &
136 obsSpaceData, &
137 ensemble )
138
139 if ( .not. staticHBHT .and. .not. staticOMPE .and. .not. ensemble .and. &
140 .not. staticHBHT_ch .and. .not. staticOMPE_ch ) &
141 call utl_abort('ose_computeStddev: no OmP std dev or sqrt(HBHT) was initialized')
142
143 !- 2. Select/Blend HBHT.
144
145 if ( staticHBHT .and. .not. ensemble ) then
146 ! Bnmc only
147 write(*,*)
148 write(*,*) 'ose_computeStddev: Using B_static ONLY'
149 ! HBnmcHT std. dev. already in OBS_HPHT
150 else if ( .not. staticHBHT .and. ensemble ) then
151 write(*,*)
152 write(*,*) 'ose_computeStddev: Using B_ensemble ONLY'
153 ! Transfer HBensHT std. dev. values in OBS_WORK to OBS_HPHT
154 do index_body = 1, obs_numBody(obsSpaceData)
155 HBHT_ensemble = obs_bodyElem_r(obsSpaceData,OBS_WORK,index_body)
156 call obs_bodySet_r(obsSpaceData,OBS_HPHT,index_body, HBHT_ensemble)
157 end do
158 else if ( staticHBHT .and. ensemble ) then
159 ! Read Namelist first
160 hybrid_mode = 'WEIGHTED_SUM' ! default value
161 nulnam = 0
162 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
163 read(nulnam,nml=namhbht,iostat=ierr)
164 if ( ierr /= 0) call utl_abort('ose_computeStddev: Error reading namelist')
165 if ( mmpi_myid == 0 ) write(*,nml=namhbht)
166 ierr = fclos(nulnam)
167
168 write(*,*)
169 write(*,*) 'ose_computeStddev: Using hybrid approach (blend of B_static and B_ensemble) in mode = ', trim(hybrid_mode)
170
171 do index_body = 1, obs_numBody(obsSpaceData)
172 HBHT_static = obs_bodyElem_r(obsSpaceData,OBS_HPHT,index_body)
173 HBHT_ensemble = obs_bodyElem_r(obsSpaceData,OBS_WORK,index_body)
174 if ( HBHT_static <= 0.0d0 ) then
175 call obs_bodySet_r(obsSpaceData,OBS_HPHT,index_body, HBHT_ensemble)
176 else
177 select case ( trim(hybrid_mode) )
178 case ('WEIGHTED_SUM')
179 HBHT_hybrid = sqrt(HBHT_static**2 + HBHT_ensemble**2)
180 case ('MAX_VALUE')
181 HBHT_hybrid = max(HBHT_static,HBHT_ensemble)
182 case default
183 write(*,*)
184 write(*,*) 'ose_computeStddev: Unknown hybrid_mode ', trim(hybrid_mode)
185 call utl_abort('ose_compute_HBHT')
186 end select
187 call obs_bodySet_r(obsSpaceData,OBS_HPHT,index_body, HBHT_hybrid)
188 end if
189 end do
190
191 end if
192
193 end subroutine ose_computeStddev
194
195 !--------------------------------------------------------------------------
196 ! ose_setStaticErrorStddev
197 !--------------------------------------------------------------------------
198 subroutine ose_setStaticErrorStddev( columnTrlOnAnlIncLev, obsSpaceData, statusHBHT, statusHBHT_ch, statusOMPE_ch )
199 !
200 !:Purpose: To assign or compute the OmP error standard deviations in
201 ! observation space where requested. If not possible or available,
202 ! compute background-error standard deviation.
203 !
204 !:Note: OmP error std dev assigment currently only available for the CH obs family.
205 !
206 !-------------------------------------------------------------------------
207 implicit none
208
209 ! Arguments:
210 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
211 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data, output saved in OBS_HPHT column
212 logical, intent(inout) :: statusHBHT
213 logical, intent(inout) :: statusHBHT_ch
214 logical, intent(inout) :: statusOMPE_ch
215
216 ! Locals:
217 integer :: famIndex
218 character(len=4), allocatable :: availableOMPE(:)
219
220 allocate(availableOMPE(ofl_numFamily))
221 availableOMPE(:) = ' '
222
223 ! Assignment of OBS_OMPE in obsSpaceData according to obs family where possible and required
224
225 do famIndex=1,ofl_numFamily
226 if ( .not. obs_famExist(obsSpaceData,ofl_familyList(famIndex),localMPI_opt=.true.) ) cycle
227
228 if ( ofl_familyList(famIndex) == 'CH' ) then
229 write(*,*)
230 write(*,*) 'ose_setStaticErrorStddev: Setting of OmP error std dev begins for family ',ofl_familyList(famIndex)
231 availableOMPE(famIndex) = ose_setOmPstddevCH(obsSpaceData)
232 else
233 ! Not available
234 availableOMPE(famIndex) = 'None'
235 end if
236 if ( trim(availableOMPE(famIndex)) == 'Some' .or. &
237 trim(availableOMPE(famIndex)) == 'All' ) statusOMPE_ch = .true.
238
239 if ( trim(availableOMPE(famIndex)) == '' ) then
240 write(*,*) 'ose_setStaticErrorStddev: No ',ofl_familyList(famIndex),' obs. Setting of error std dev not required.'
241 else if ( trim(availableOMPE(famIndex)) == 'None' ) then
242 if ( ofl_familyList(famIndex) == 'CH' ) then
243 write(*,*) 'ose_setStaticErrorStddev: Setting of ',ofl_familyList(famIndex),' OmP error std dev to be estimated via HBHT calc for all obs'
244 end if
245 else if ( trim(availableOMPE(famIndex)) == 'Some' ) then
246 write(*,*) 'ose_setStaticErrorStddev: Setting of ',ofl_familyList(famIndex),' OmP error std dev partially completed (some HBHT calc needed to complete)'
247 else
248 write(*,*) 'ose_setStaticErrorStddev: Setting of ',ofl_familyList(famIndex),' OmP error std dev completed (no HBHT calc needed)'
249 end if
250
251 end do
252
253 ! Computation of background-error standard deviation for obs families or groups of families when required
254
255 ! HBHT from the Bnmc matrix for weather
256 ! obsSpaceData - INOUT (HBnmcHT std. dev. output in OBS_HPHT)
257
258 if ( any(ofl_familyList /= 'CH' .and. ofl_familyList /= 'TO' .and. &
259 ( availableOMPE == 'Some' .or. availableOMPE == 'None' ) ) ) &
260 call ose_compute_hbht_static( columnTrlOnAnlIncLev, obsSpaceData, statusHBHT )
261
262 ! HBHT from the B matrix for constituents
263
264 if ( any(ofl_familyList == 'CH' .and. ( availableOMPE == 'Some' .or. availableOMPE == 'None' ) ) ) &
265 call ose_compute_hbht_static_chem( columnTrlOnAnlIncLev, obsSpaceData, statusHBHT_ch )
266
267 deallocate(availableOMPE)
268
269 end subroutine ose_setStaticErrorStddev
270
271 !--------------------------------------------------------------------------
272 ! ose_compute_hbht_static
273 !--------------------------------------------------------------------------
274 subroutine ose_compute_hbht_static(columnTrlOnAnlIncLev,lobsSpaceData,active)
275 !
276 !:Purpose: To compute background-error stddev in observation space using
277 ! fixed statistics specific in stats file.
278 !
279 implicit none
280
281 ! Arguments:
282 type(struct_obs), intent(inout) :: lobsSpaceData
283 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
284 logical, intent(out) :: active
285
286 ! Locals:
287 type(struct_vco), pointer :: vco_anl
288 type(struct_columnData) :: column
289 type(struct_gsv) :: statevector
290 INTEGER JLAT, JLON, JLEV, JOBS
291 CHARACTER*12 CLETIKET
292 CHARACTER*2 CLTYPVAR
293 CHARACTER*1 CLGRTYP
294 CHARACTER*4 CLNOMVAR
295 INTEGER IULSSF,IDATEO
296 INTEGER FSTPRM,FNOM,FSTOUV
297 INTEGER IKEY,IERR,IDATE
298 REAL*8, allocatable :: ZBUFFER(:,:)
299 real*8, pointer :: height_column(:), tt_column(:), field_ptr(:,:,:)
300 INTEGER INI,INJ,INK, INPAS, INBITS, IDATYP, IDEET
301 INTEGER IP1,IP2,IP3,IG1,IG2,IG3,IG4,ISWA,ILENGTH,IDLTF
302 INTEGER IUBC,IEXTR1,IEXTR2,IEXTR3
303 integer :: nLev_M,nLev_T,shift_level,Vcode_anl
304 integer :: cvdim
305 real(8), allocatable :: scaleFactor(:)
306
307 !- Get the appropriate Vertical Coordinate
308 vco_anl => col_getVco(columnTrlOnAnlIncLev)
309
310 ! Note : Here we can use the global B_hi even if we are in LAM mode since,
311 ! in BackgroundCheck mode, the only purpose of bhi_setup is to read
312 ! the scaleFactor and to check the consistency between the vco from
313 ! analysisgrid and cov file
314 call bhi_setup( hco_anl,vco_anl, & ! IN
315 cvdim, & ! OUT
316 mode_opt='BackgroundCheck' ) ! IN
317
318 if ( cvdim > 0 ) then
319 write(*,*)
320 write(*,*) 'Computing HBHT from Bnmc - Start'
321 active = .true.
322 else
323 if ( mmpi_myid == 0 ) write(*,*) 'ose_compute_hbht_static: option NOT ACTIVATED'
324 active = .false.
325 return
326 end if
327
328 nlev_T = col_getNumLev(columnTrlOnAnlIncLev,'TH')
329 nlev_M = col_getNumLev(columnTrlOnAnlIncLev,'MM')
330
331 allocate(scaleFactor(max(nLev_M,nLev_T)))
332 call bhi_getScaleFactor(scaleFactor)
333
334 Vcode_anl = vco_anl%Vcode
335 if ( Vcode_anl == 5001 ) then
336 shift_level = 0
337 else if ( Vcode_anl == 5002) then
338 shift_level = 1
339 else if ( Vcode_anl == 5005 ) then
340 shift_level = 0
341 else
342 write(*,*) 'Vcode_anl = ',Vcode_anl
343 call utl_abort('ose_compute_hbht_static: unknown vertical coordinate type!')
344 end if
345
346 allocate(ZBUFFER(HCO_ANL%NI,HCO_ANL%NJ))
347
348 call gsv_allocate(statevector, 1, hco_anl, vco_anl, &
349 allocHeight_opt=.false., allocPressure_opt=.false.)
350 call gsv_zero(statevector)
351
352 call col_setVco(column,col_getVco(columnTrlOnAnlIncLev))
353 call col_allocate(column,col_getNumCol(columnTrlOnAnlIncLev))
354
355 ! Set the value of OBS_LYR required by setfge routines
356
357 call oop_vobslyrs(columnTrlOnAnlIncLev,lobsSpaceData,beSilent=.false.)
358
359 ! 1. Opening the statistics file
360
361 IULSSF=0
362 IERR=FNOM(iulssf,'./bgcov','RND+OLD+R/O',0)
363 IF ( IERR .EQ. 0 ) THEN
364 write(*,*) 'IBGST - File : ./bgcov'
365 write(*,*) ' opened as unit file ',iulssf
366 ierr = fstouv(iulssf,'RND+OLD')
367 ELSE
368 CALL utl_abort('HBHT_static:NO BACKGROUND STAT FILE!!')
369 ENDIF
370
371 ! 2.1 Background error standard deviations
372
373 CLETIKET = 'BGCK_STDDEV'
374 write(*,*) 'HBHT_static: CLETIKET = ',CLETIKET
375 IDATE = -1
376 IP2 = -1
377 IP3 = -1
378 CLTYPVAR =' '
379
380 ! READ IN STANDARD DEVIATION FOR EACH OBSERVATION TYPE
381
382 clnomvar = 'UU'
383 write(*,*) clnomvar
384 call gsv_getField(statevector,field_ptr,'UU')
385 do jlev = 1, nlev_M
386 ip1 = vco_anl%ip1_M(jlev)
387 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
388 ierr = fstprm(ikey,idateo,ideet,inpas &
389 ,ini,inj,ink, inbits, idatyp &
390 ,ip1,ip2,ip3,cltypvar,clnomvar,cletiket,clgrtyp &
391 ,ig1,ig2,ig3,ig4,iswa,ilength,idltf &
392 ,iubc,iextr1,iextr2,iextr3)
393 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
394 write(*,*)
395 write(*,*) 'HBHT_static: Invalid dimensions for...'
396 write(*,*) 'nomvar =', trim(clnomvar)
397 write(*,*) 'etiket =', trim(cletiket)
398 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
399 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_M
400 call utl_abort('ose_compute_hbht_static')
401 end if
402 do jlat = 1, hco_anl%nj
403 do jlon = 1, hco_anl%ni
404 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev+shift_level)*zbuffer(jlon,jlat)*MPC_M_PER_S_PER_KNOT_R8
405 end do
406 end do
407 end do
408
409 clnomvar = 'VV'
410 write(*,*) clnomvar
411 call gsv_getField(statevector,field_ptr,'VV')
412 do jlev = 1, nlev_M
413 ip1 = vco_anl%ip1_M(jlev)
414 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
415 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
416 write(*,*)
417 write(*,*) 'HBHT_static: Invalid dimensions for...'
418 write(*,*) 'nomvar =', trim(clnomvar)
419 write(*,*) 'etiket =', trim(cletiket)
420 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
421 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_M
422 call utl_abort('ose_compute_hbht_static')
423 end if
424
425 do jlat = 1, hco_anl%nj
426 do jlon = 1, hco_anl%ni
427 field_ptr(jlon,jlat,jlev) = scalefactor(jlev+shift_level)*zbuffer(jlon,jlat)*MPC_M_PER_S_PER_KNOT_R8
428 end do
429 end do
430 end do
431
432 clnomvar = 'ES'
433 write(*,*) clnomvar
434 call gsv_getField(statevector,field_ptr,'HU')
435 do jlev = 1, nlev_T
436 ip1 = vco_anl%ip1_T(jlev)
437 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
438 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
439 write(*,*)
440 write(*,*) 'HBHT_static: Invalid dimensions for...'
441 write(*,*) 'nomvar =', trim(clnomvar)
442 write(*,*) 'etiket =', trim(cletiket)
443 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
444 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
445 call utl_abort('ose_compute_hbht_static')
446 end if
447 do jlat = 1, hco_anl%nj
448 do jlon = 1, hco_anl%ni
449 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
450 end do
451 end do
452 end do
453
454 ! height is put into TT slot in gridStateVector
455 clnomvar = 'GZ'
456 write(*,*) clnomvar
457 call gsv_getField(statevector,field_ptr,'TT')
458 do jlev = 1, nlev_T
459 ip1 = vco_anl%ip1_T(jlev)
460 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
461 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
462 write(*,*)
463 write(*,*) 'HBHT_static: Invalid dimensions for...'
464 write(*,*) 'nomvar =', trim(clnomvar)
465 write(*,*) 'etiket =', trim(cletiket)
466 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
467 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
468 call utl_abort('ose_compute_hbht_static')
469 end if
470 do jlat = 1, hco_anl%nj
471 do jlon = 1, hco_anl%ni
472 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)*ec_rg*10.d0
473 end do
474 end do
475 end do
476
477 call s2c_bgcheck_bilin(column,statevector,lobsSpaceData)
478
479 ! copy height data from TT to height slot in columnData
480 do jobs= 1, col_getNumCol(column)
481 height_column => col_getColumn(column,jobs,'Z_T')
482 tt_column => col_getColumn(column,jobs,'TT')
483 do jlev = 1,col_getNumLev(column,'TH')
484 height_column(jlev)=tt_column(jlev)
485 enddo
486 enddo
487
488 ! SET THE FIRST-GUESS ERRORS FOR CONVENTIONAL DATA ON PRESSURE LEVELS
489 ! --------------------------------------------------------------------
490
491 call setfgefam('AI',column,columnTrlOnAnlIncLev,lobsSpaceData)
492 call setfgefam('SW',column,columnTrlOnAnlIncLev,lobsSpaceData)
493 call setfgefam('UA',column,columnTrlOnAnlIncLev,lobsSpaceData)
494 call setfgefam('SF',column,columnTrlOnAnlIncLev,lobsSpaceData)
495 call setfgefam('HU',column,columnTrlOnAnlIncLev,lobsSpaceData)
496 call setfgefamz('PR',column,columnTrlOnAnlIncLev,lobsSpaceData)
497 call setfgefamz('AL',column,columnTrlOnAnlIncLev,lobsSpaceData)
498 call setfgefamz('RA',column,columnTrlOnAnlIncLev,lobsSpaceData)
499
500 ! SET THE FIRST-GUESS ERRORS FOR RADIO OCCULTATION DATA
501 ! -----------------------------------------------------
502
503 call setfgedif('RO',columnTrlOnAnlIncLev,lobsSpaceData)
504
505 ! DO TEMPERATURE FIRST-GUESS ERROR
506 ! ---------------------------------
507
508 clnomvar = 'TT'
509 write(*,*) clnomvar
510 call gsv_getField(statevector,field_ptr,'TT')
511 do jlev = 1, nlev_T
512 ip1 = vco_anl%ip1_T(jlev)
513 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
514 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
515 write(*,*)
516 write(*,*) 'HBHT_static: Invalid dimensions for...'
517 write(*,*) 'nomvar =', trim(clnomvar)
518 write(*,*) 'etiket =', trim(cletiket)
519 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
520 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
521 call utl_abort('ose_compute_hbht_static')
522 end if
523 do jlat = 1, hco_anl%nj
524 do jlon = 1, hco_anl%ni
525 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
526 end do
527 end do
528 end do
529
530 call s2c_bgcheck_bilin(column,statevector,lobsSpaceData)
531 call setfgett(column,columnTrlOnAnlIncLev,lobsSpaceData)
532
533 ! RELOAD DATA TO DO SURFACE FIRST-GUESS ERRORS
534 ! ---------------------------------------------
535
536 clnomvar = 'P0'
537 write(*,*) clnomvar
538 ip1 = -1
539 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
540 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
541 write(*,*)
542 write(*,*) 'HBHT_static: Invalid dimensions for...'
543 write(*,*) 'nomvar =', trim(CLNOMVAR)
544 write(*,*) 'etiket =', trim(CLETIKET)
545 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
546 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, 1
547 call utl_abort('ose_compute_hbht_static')
548 end if
549 call gsv_getField(statevector,field_ptr,'P0')
550 do jlev = 1, ink
551 do jlat = 1, hco_anl%nj
552 do jlon = 1, hco_anl%ni
553 field_ptr(jlon,jlat,jlev) = scaleFactor(max(nLev_M,nLev_T))*zbuffer(jlon,jlat)*MPC_PA_PER_MBAR_R8
554 end do
555 end do
556 end do
557
558 clnomvar = 'UU'
559 write(*,*) clnomvar
560 call gsv_getField(statevector,field_ptr,'UU')
561 do jlev = 1, nlev_M
562 ip1 = vco_anl%ip1_M(jlev)
563 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
564 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
565 write(*,*)
566 write(*,*) 'HBHT_static: Invalid dimensions for...'
567 write(*,*) 'nomvar =', trim(clnomvar)
568 write(*,*) 'etiket =', trim(cletiket)
569 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
570 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_M
571 call utl_abort('ose_compute_hbht_static')
572 end if
573 do jlat = 1, hco_anl%nj
574 do jlon = 1, hco_anl%ni
575 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev+shift_level)*zbuffer(jlon,jlat)*MPC_M_PER_S_PER_KNOT_R8
576 end do
577 end do
578 end do
579
580 clnomvar = 'VV'
581 write(*,*) clnomvar
582 call gsv_getField(statevector,field_ptr,'VV')
583 do jlev = 1, nlev_M
584 ip1 = vco_anl%ip1_M(jlev)
585 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
586 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
587 write(*,*)
588 write(*,*) 'HBHT_static: Invalid dimensions for...'
589 write(*,*) 'nomvar =', trim(clnomvar)
590 write(*,*) 'etiket =', trim(cletiket)
591 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
592 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_M
593 call utl_abort('ose_compute_hbht_static')
594 end if
595 do jlat = 1, hco_anl%nj
596 do jlon = 1, hco_anl%ni
597 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev+shift_level)*zbuffer(jlon,jlat)*MPC_M_PER_S_PER_KNOT_R8
598 end do
599 end do
600 end do
601
602 clnomvar = 'TT'
603 write(*,*) clnomvar
604 call gsv_getField(statevector,field_ptr,'TT')
605 do jlev = 1, nlev_T
606 ip1 = vco_anl%ip1_T(jlev)
607 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
608 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
609 write(*,*)
610 write(*,*) 'HBHT_static: Invalid dimensions for...'
611 write(*,*) 'nomvar =', trim(clnomvar)
612 write(*,*) 'etiket =', trim(cletiket)
613 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
614 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
615 call utl_abort('ose_compute_hbht_static')
616 end if
617 do jlat = 1, hco_anl%nj
618 do jlon = 1, hco_anl%ni
619 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
620 end do
621 end do
622 end do
623
624 clnomvar = 'ES'
625 write(*,*) clnomvar
626 call gsv_getField(statevector,field_ptr,'HU')
627 do jlev = 1, nlev_T
628 ip1 = vco_anl%ip1_T(jlev)
629 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
630 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
631 write(*,*)
632 write(*,*) 'HBHT_static: Invalid dimensions for...'
633 write(*,*) 'nomvar =', trim(clnomvar)
634 write(*,*) 'etiket =', trim(cletiket)
635 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
636 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
637 call utl_abort('ose_compute_hbht_static')
638 end if
639 do jlat = 1, hco_anl%nj
640 do jlon = 1, hco_anl%ni
641 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
642 end do
643 end do
644 end do
645
646 clnomvar = 'LVIS'
647 if (gsv_varExist(statevector,clnomvar)) then
648 write(*,*) clnomvar
649 call gsv_getField(statevector,field_ptr,clnomvar)
650 do jlev = 1, nlev_T
651 ip1 = vco_anl%ip1_T(jlev)
652 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
653 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
654 write(*,*)
655 write(*,*) 'HBHT_static: Invalid dimensions for...'
656 write(*,*) 'nomvar =', trim(clnomvar)
657 write(*,*) 'etiket =', trim(cletiket)
658 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
659 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
660 call utl_abort('ose_compute_hbht_static')
661 end if
662 do jlat = 1, hco_anl%nj
663 do jlon = 1, hco_anl%ni
664 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
665 end do
666 end do
667 end do
668 end if
669
670 clnomvar = 'WGE'
671 if (gsv_varExist(statevector,clnomvar)) then
672 write(*,*) clnomvar
673 ip1 = -1
674 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
675 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
676 write(*,*)
677 write(*,*) 'HBHT_static: Invalid dimensions for...'
678 write(*,*) 'nomvar =', trim(CLNOMVAR)
679 write(*,*) 'etiket =', trim(CLETIKET)
680 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
681 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, 1
682 call utl_abort('ose_compute_hbht_static')
683 end if
684 call gsv_getField(statevector,field_ptr,clnomvar)
685 do jlev = 1, ink
686 do jlat = 1, hco_anl%nj
687 do jlon = 1, hco_anl%ni
688 field_ptr(jlon,jlat,jlev) = scaleFactor(max(nLev_M,nLev_T))*zbuffer(jlon,jlat)
689 end do
690 end do
691 end do
692 end if
693
694 call s2c_bgcheck_bilin(column,statevector,lobsSpaceData)
695
696 ! SET THE FIRST-GUESS ERRORS FOR THE SURFACE DATA
697 ! ------------------------------------------------
698
699 call setfgesurf(column,columnTrlOnAnlIncLev,lobsSpaceData)
700
701 ! READ IN LN Q FIRST-GUESS ERRORS FOR SETFGEGPS
702 ! ---------------------------------------------
703
704 clnomvar = 'LQ'
705 write(*,*) clnomvar
706 call gsv_getField(statevector,field_ptr,'HU')
707 do jlev = 1, nlev_T
708 ip1 = vco_anl%ip1_T(jlev)
709 ikey = utl_fstlir(zbuffer,iulssf,ini,inj,ink,idate,cletiket,ip1,ip2,ip3,cltypvar,clnomvar)
710 if (ini /= hco_anl%ni .or. inj /= hco_anl%nj .or. ink /= 1) then
711 write(*,*)
712 write(*,*) 'HBHT_static: Invalid dimensions for...'
713 write(*,*) 'nomvar =', trim(clnomvar)
714 write(*,*) 'etiket =', trim(cletiket)
715 write(*,*) 'Found ni,nj,nk =', ini, inj, ink
716 write(*,*) 'Should be =', hco_anl%ni, hco_anl%nj, nlev_T
717 call utl_abort('ose_compute_hbht_static')
718 end if
719 do jlat = 1, hco_anl%nj
720 do jlon = 1, hco_anl%ni
721 field_ptr(jlon,jlat,jlev) = scaleFactor(jlev)*zbuffer(jlon,jlat)
722 end do
723 end do
724 end do
725
726 call s2c_bgcheck_bilin(column,statevector,lobsSpaceData)
727
728 ! OPTIONAL TEST OF THE GB-GPS ZTD OPERATOR JACOBIAN
729 ! -------------------------------------------------
730
731 if (gps_gb_ltestop) call setfgegps(column,columnTrlOnAnlIncLev,lobsSpaceData)
732
733 deallocate(scaleFactor)
734 call col_deallocate(column)
735 deallocate(zbuffer)
736
737 write(*,*)
738 write(*,*) 'Computing HBHT from Bnmc - END'
739
740 end subroutine ose_compute_hbht_static
741
742 !--------------------------------------------------------------------------
743 ! ose_compute_hbht_static_chem
744 !--------------------------------------------------------------------------
745 subroutine ose_compute_hbht_static_chem(columnTrlOnAnlIncLev,obsSpaceData,active)
746 !
747 !:Purpose: To compute the background error standard deviations in
748 ! observation space, sqrt(diag(H*B_static*H^T)).
749 !
750 implicit none
751
752 ! Arguments:
753 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev ! column at observation location
754 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data, output saved in OBS_HPHT column
755 logical, intent(out) :: active ! flag to indicate if chemical constituents are to be used
756
757 ! Locals:
758 type(struct_vco), pointer :: vco_anl
759
760 !- Get the appropriate Vertical Coordinate
761 vco_anl => col_getVco(columnTrlOnAnlIncLev)
762
763 call bcsc_setupCH( hco_anl,vco_anl,active,'BackgroundCheck' )
764
765 if (active) then
766 write(*,*)
767 write(*,*) 'Computing H*B*H^T using B_static_chm - Start'
768 else
769 if ( mmpi_myid == 0 ) write(*,*) 'ose_compute_HBHT_static_chem: option NOT ACTIVATED'
770 return
771 end if
772
773 call oopc_CHobsoperators(columnTrlOnAnlIncLev,obsSpaceData,'HBHT')
774
775 write(*,*)
776 write(*,*) 'Computing H*B*H^T using B_static_chm - End'
777
778 RETURN
779 end subroutine ose_compute_hbht_static_chem
780
781 !--------------------------------------------------------------------------
782 ! ose_compute_hbht_ensemble
783 !--------------------------------------------------------------------------
784 subroutine ose_compute_hbht_ensemble(columnTrlOnAnlIncLev,obsSpaceData,active)
785 !
786 !:Purpose: To compute background-error stddev in observation space using
787 ! ensemble-based statistics.
788 !
789 implicit none
790
791 ! Arguments:
792 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev ! Background columns interpolated to anl levels, obs locations
793 type(struct_obs), intent(inout) :: obsSpaceData ! Observation-related data
794 logical, intent(out) :: active
795
796 ! Locals:
797 type(struct_columnData) :: column
798 type(struct_gsv) :: statevector
799 type(struct_vco), pointer :: vco_anl
800 real(8), allocatable :: HBHT_ens(:)
801 integer :: memberIndex, index_body
802 integer, allocatable :: cvdim(:)
803
804 !
805 !- 1. Initialization
806 !
807
808 !- 1.1 Get vertical analysis grid attributes
809 vco_anl => col_getVco(columnTrlOnAnlIncLev)
810
811 !- 1.2 Initialize/Read the flow-dependent ensemble perturbations
812 call ben_Setup( hco_anl, hco_anl, & ! IN
813 vco_anl, & ! IN
814 cvdim, & ! OUT
815 'BackgroundCheck' ) ! IN
816
817 if ( cvdim(1) > 0 ) then
818 write(*,*)
819 write(*,*) 'Computing HBHT from ensemble perturbations - START'
820 active = .true.
821 else
822 if ( mmpi_myid == 0 ) write(*,*) 'ose_compute_hbht_ensemble: option NOT ACTIVATED'
823 active = .false.
824 return
825 end if
826
827 !- 1.3 Create a gridstatevector to store the ensemble perturbations
828 call gsv_allocate(statevector, tim_nstepobsinc, hco_anl, vco_anl, &
829 mpi_local_opt=.true., allocHeight_opt=.false., allocPressure_opt=.false.)
830
831 !- 1.4 Create column vectors to store the ens perturbation interpolated to obs horizontal locations
832 call col_setVco(column,vco_anl)
833 call col_allocate(column,col_getNumCol(columnTrlOnAnlIncLev),mpiLocal_opt=.true.)
834
835 !- 1.5 Create a working a array to sum H ensPert HT
836 allocate(HBHT_ens(obs_numBody(obsSpaceData)))
837 HBHT_ens(:) = 0.d0
838
839 !- 1.6
840 call oti_timeBinning(obsSpaceData,tim_nstepobsinc)
841
842 !
843 !- 2. Compute HBHT from the ensemble perturbations
844 !
845 do memberIndex = 1, ben_getnEns()
846
847 !- 2.1 Extract perturbations from the current memberIndex
848 write(*,*)
849 write(*,*) 'Reading ensemble perturbation from member = ', memberIndex
850 call ben_getPerturbation( statevector, &
851 memberIndex, &
852 'ConstantValue' )
853
854 !- 2.2 Interpolation to the observation horizontal locations
855 ! column - OUT (H_horiz EnsPert)
856 call s2c_tl( statevector, &
857 column, &
858 columnTrlOnAnlIncLev, obsSpaceData )
859
860 !- 2.3 Interpolation to observation space
861 ! obsSpaceData - OUT (Save as OBS_WORK: H_vert H_horiz EnsPert = H EnsPert)
862 call oop_Htl( column, columnTrlOnAnlIncLev, &
863 obsSpaceData, &
864 1 )
865
866 !- 2.4 alpha * HBH^T = sum(OBS_WORK^2)
867 do index_body = 1, obs_numBody(obsSpaceData)
868 HBHT_ens(index_body) = HBHT_ens(index_body) + &
869 (obs_bodyElem_r(obsSpaceData,OBS_WORK,index_body))**2
870 end do
871
872 end do
873
874 !- 2.5 Insert the standard deviations in OBS_WORK
875 do index_body = 1, obs_numBody(obsSpaceData)
876 call obs_bodySet_r(obsSpaceData,OBS_WORK,index_body,sqrt(HBHT_ens(index_body)))
877 end do
878
879 !
880 !- 3. Ending/Deallocation
881 !
882 deallocate(HBHT_ens)
883 call col_deallocate(column)
884 call gsv_deallocate(statevector)
885
886 write(*,*)
887 write(*,*) 'Computing HBHT from ensemble perturbations - END'
888
889 end subroutine ose_compute_hbht_ensemble
890
891 !-------------------- Weather obs FGE std dev routines --------------------
892
893 !--------------------------------------------------------------------------
894 ! setfgefam
895 !--------------------------------------------------------------------------
896 subroutine setfgefam(cdfam,column,columnTrlOnAnlIncLev,lobsSpaceData)
897 !
898 !:Purpose: To interpolate vertically the contents of "column" to
899 ! the pressure levels of the observations. Then to compute
900 ! THE FIRST GUESS ERROR VARIANCES. A linear interpolation in ln(p)
901 ! is performed.
902 !
903 implicit none
904
905 ! Arguments:
906 character(len=2), intent(in) :: cdfam
907 type(struct_columnData), intent(in) :: column
908 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
909 type(struct_obs), intent(inout) :: lobsSpaceData
910
911 ! Locals:
912 INTEGER IPB,IPT
913 INTEGER INDEX_HEADER,ITYP,IK
914 INTEGER INDEX_BODY
915 REAL*8 ZWB,ZWT, columnElem
916 REAL*8 ZLEV,ZPB,ZPT
917 character(len=4) :: varLevel
918
919 ! loop over all header indices of the CDFAM family
920 call obs_set_current_header_list(lobsSpaceData,CDFAM)
921 HEADER: do
922 index_header = obs_getHeaderIndex(lobsSpaceData)
923 if (index_header < 0) exit HEADER
924
925 ! loop over all body indices for this index_header
926 call obs_set_current_body_list(lobsSpaceData, index_header)
927 BODY: do
928 index_body = obs_getBodyIndex(lobsSpaceData)
929 if (index_body < 0) exit BODY
930 !
931 !* 1. Computation of sigmap
932 ! . -----------------------------
933 !
934 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,index_body) == obs_assimilated .AND. &
935 obs_bodyElem_i(lobsSpaceData,OBS_VCO,index_body) .EQ. 2 ) then
936 IF (obs_bodyElem_i(lobsSpaceData,OBS_XTR,index_body) .NE. 0) THEN
937 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,index_body)
938 varLevel = vnl_varLevelFromVarnum(ityp)
939 IK = col_getNumLev(columnTrlOnAnlIncLev,varLevel)
940 IPB = IK + col_getOffsetFromVarno(columnTrlOnAnlIncLev,ityp)
941 if(ITYP .ne. BUFR_NEGZ) then
942 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body,col_getElem(column,IPB,INDEX_HEADER))
943 else
944 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body,col_getHeight(column,IK,INDEX_HEADER,'TH'))
945 endif
946 ELSE
947 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,index_body)
948 varLevel = vnl_varLevelFromVarnum(ityp)
949 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,index_body)
950 IK = obs_bodyElem_i(lobsSpaceData,OBS_LYR,index_body)
951 IPT = IK + col_getOffsetFromVarno(columnTrlOnAnlIncLev,ityp)
952 IPB = IPT+1
953 ZPT = col_getPressure(columnTrlOnAnlIncLev,IK,INDEX_HEADER,varLevel)
954 ZPB = col_getPressure(columnTrlOnAnlIncLev,IK+1,INDEX_HEADER,varLevel)
955 ZWB = LOG(ZLEV/ZPT)/LOG(ZPB/ZPT)
956 ZWT = 1.0D0 - ZWB
957
958 ! FIRST GUESS ERROR VARIANCE
959
960 if(ITYP .ne. BUFR_NEGZ) then
961 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body, &
962 (ZWB*col_getElem(column,IPB,INDEX_HEADER) + ZWT*col_getElem(column,IPT,INDEX_HEADER)))
963 else
964 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body, &
965 (ZWB*col_getHeight(column,IK+1,INDEX_HEADER,'TH') + ZWT*col_getHeight(column,IK,INDEX_HEADER,'TH')))
966 endif
967 if(obs_bodyElem_r(lobsSpaceData,OBS_HPHT,index_body).le.0.d0) then
968 write(*,*) 'SETFGEFAM: CDFAM = ',CDFAM
969 write(*,*) 'SETFGEFAM: IPB,IPT,ZWB,ZWT,ITYP,ZLEV=',IPB,IPT,ZWB,ZWT,ITYP,ZLEV
970 columnElem = col_getElem(column,IPB,INDEX_HEADER)
971 write(*,*) 'SETFGEFAM: column_all(IPB,INDEX_HEADER)=',columnElem
972 columnElem = col_getElem(column,IPT,INDEX_HEADER)
973 write(*,*) 'SETFGEFAM: column_all(IPT,INDEX_HEADER)=',columnElem
974 columnElem = col_getHeight(column,IK+1,INDEX_HEADER,'TH')
975 write(*,*) 'SETFGEFAM: get_height(IK+1,INDEX_HEADER)=',columnElem
976 columnElem = col_getHeight(column,IK ,INDEX_HEADER,'TH')
977 write(*,*) 'SETFGEFAM: get_height(IK ,INDEX_HEADER)=',columnElem
978 CALL utl_abort('SETFGEFAM: First-guess stdev bad value')
979 endif
980 ENDIF
981 ENDIF
982
983 END DO BODY
984
985 END DO HEADER
986
987 end subroutine setfgefam
988
989 !--------------------------------------------------------------------------
990 ! setfgefamz
991 !--------------------------------------------------------------------------
992 subroutine setfgefamz(cdfam, column, columnTrlOnAnlIncLev, obsSpaceData)
993 !
994 !:Purpose: To interpolate vertically the contents of "column" to the levels
995 ! of the observations (in meters). Then to compute THE FIRST GUESS
996 ! ERROR VARIANCES. A linear interpolation in z is performed.
997 !
998 implicit none
999
1000 ! Arguments:
1001 character(len=2), intent(in) :: cdfam
1002 type(struct_columnData), intent(in) :: column
1003 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1004 type(struct_obs), intent(inout) :: obsSpaceData
1005
1006 ! Locals:
1007 integer :: IPB, IPT
1008 integer :: headerIndex, ITYP, IK
1009 integer :: bodyIndex
1010 integer :: bodyIndexStart, bodyIndexEnd, bodyIndex2
1011 real*8 :: ZWB, ZWT, sigmap_uuT, sigmap_vvT , sigmap_uuB, sigmap_vvB
1012 real*8 :: ZLEV, ZPB, ZPT
1013 real(8) :: azimuth, fge_uu, fge_vv, fge_fam
1014 character(len=4) :: varLevel
1015
1016 ! loop over all header indices of the CDFAM family
1017 call obs_set_current_header_list(obsSpaceData, CDFAM)
1018 HEADER: do
1019 headerIndex = obs_getHeaderIndex(obsSpaceData)
1020 if ( headerIndex < 0 ) exit HEADER
1021
1022 if ( cdfam == 'RA' ) then
1023 ! Azimuth of the radar beam
1024 azimuth = obs_headElem_r(obsSpaceData, OBS_RZAM, headerIndex )
1025 end if
1026
1027 ! loop over all body indices for this headerIndex
1028 call obs_set_current_body_list(obsSpaceData, headerIndex)
1029 BODY: do
1030 bodyIndex = obs_getBodyIndex(obsSpaceData)
1031 if ( bodyIndex < 0 ) exit BODY
1032
1033 !* 1. Computation of sigmap
1034 ! . -----------------------------
1035 if ( obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated .and. &
1036 obs_bodyElem_i(obsSpaceData, OBS_VCO, bodyIndex) == 1 )then
1037 ITYP = obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex)
1038 if ( ITYP == bufr_radarPrecip ) cycle BODY
1039 varLevel = vnl_varLevelFromVarnum(ityp)
1040
1041 ! Interpolate the background-covariance statistics
1042 if ( obs_bodyElem_i(obsSpaceData, OBS_XTR, bodyIndex) /= 0 ) then
1043 IK=col_getNumLev(columnTrlOnAnlIncLev, varLevel)-1
1044 IPT = IK + col_getOffsetFromVarno(columnTrlOnAnlIncLev, ityp)
1045 IPB = IPT +1
1046 fge_uu = col_getElem(column, IPB, headerIndex, 'UU')
1047 fge_vv = col_getElem(column, IPB, headerIndex, 'VV')
1048
1049 else
1050 ZLEV = obs_bodyElem_r(obsSpaceData, OBS_PPP, bodyIndex)
1051 IK = obs_bodyElem_i(obsSpaceData, OBS_LYR, bodyIndex)
1052 IPT = IK + col_getOffsetFromVarno(columnTrlOnAnlIncLev, ityp)
1053 IPB = IPT+1
1054 ZPT = col_getHeight(columnTrlOnAnlIncLev, IK, headerIndex, varLevel)
1055 ZPB = col_getHeight(columnTrlOnAnlIncLev, IK+1, headerIndex, varLevel)
1056 ZWB = (ZPT-ZLEV)/(ZPT-ZPB)
1057 ZWT = 1.d0 - ZWB
1058 fge_uu = ZWB*col_getElem(column, IPB, headerIndex, 'UU') &
1059 + ZWT*col_getElem(column, IPT, headerIndex, 'UU')
1060 fge_vv = ZWB*col_getElem(column, IPB, headerIndex, 'VV') &
1061 + ZWT*col_getElem(column, IPT, headerIndex, 'VV')
1062 end if
1063
1064 ! First-Guess Error Variance
1065 if ( cdfam == 'AL' )then
1066 ! Scan body indices for the azimuth
1067 bodyIndexStart = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
1068 bodyIndexEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) &
1069 + bodyIndexStart - 1
1070 BODY_SUPP: do bodyIndex2 = bodyIndexStart, bodyIndexEnd
1071 if(obs_bodyElem_i(obsSpaceData, OBS_VNM, bodyIndex2) == 5021)then
1072 azimuth=obs_bodyElem_r(obsSpaceData, OBS_VAR, bodyIndex2) * MPC_RADIANS_PER_DEGREE_R8
1073 exit BODY_SUPP
1074 end if
1075 end do BODY_SUPP
1076
1077 fge_fam = sqrt((fge_vv*cos(azimuth))**2 + (fge_uu*sin(azimuth))**2)
1078
1079 else if( cdfam == 'PR' )then
1080 fge_fam = ZWB*col_getElem(column, IPB, headerIndex) &
1081 + ZWT*col_getElem(column, IPT, headerIndex)
1082
1083 else if( cdfam == 'RA' .and. ITYP == bufr_radvel ) then
1084
1085 ! Calculation of sigmap^2 = diag(H*P*H^t) to save sigmap in OBS_HPHT
1086 !
1087 ! H includes vertical interpolation
1088 ! and projection of U and V wind components along the direction of the beam
1089 !
1090 ! P forecast error covariance
1091 !
1092 ! H = [ ZWT*sin(az) ZWT*cos(az) ZWB*sin(az) ZWB*cos(az) ]
1093 ! diag(P) = [ sigmap_uuT^2 sigmap_vvT^2 sigmap_uuB^2 sigmap_vvB^2 ]
1094 sigmap_uuT = col_getElem(column, IPT, headerIndex, 'UU')
1095 sigmap_uuB = col_getElem(column, IPB, headerIndex, 'UU')
1096 sigmap_vvT = col_getElem(column, IPT, headerIndex, 'VV')
1097 sigmap_vvB = col_getElem(column, IPB, headerIndex, 'VV')
1098
1099 fge_fam = sqrt((sigmap_uuT*ZWT*sin(azimuth))**2 + (sigmap_vvT*ZWT*cos(azimuth))**2 + &
1100 (sigmap_uuB*ZWB*sin(azimuth))**2 + (sigmap_vvB*ZWB*cos(azimuth))**2)
1101
1102 else
1103
1104 write(*,*)"ERROR: The family", cdfam, " is not supported by setfgefamz"
1105 call utl_abort('setfgefamz')
1106
1107 end if
1108
1109 ! Store fge_fam in OBS_HPHT
1110 call obs_bodySet_r(obsSpaceData, OBS_HPHT, bodyIndex, fge_fam)
1111
1112 end if
1113
1114 end do BODY
1115
1116 end do HEADER
1117
1118 end subroutine setfgefamz
1119
1120 !--------------------------------------------------------------------------
1121 ! setfgett
1122 !--------------------------------------------------------------------------
1123 subroutine setfgett(column,columnTrlOnAnlIncLev,lobsSpaceData)
1124 !
1125 !:Purpose: To interpolate vertically the contents of "column" to the
1126 ! pressure levels of the observations. Then to compute THE FIRST
1127 ! GUESS ERROR VARIANCES. A linear interpolation in ln(p) is
1128 ! performed.
1129 !
1130 implicit none
1131
1132 ! Arguments:
1133 type(struct_columnData), intent(in) :: column
1134 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1135 type(struct_obs), intent(inout) :: lobsSpaceData
1136
1137 ! Locals:
1138 INTEGER IPB,IPT
1139 INTEGER INDEX_HEADER,ITYP,IK
1140 INTEGER INDEX_BODY
1141 REAL*8 ZWB,ZWT
1142 REAL*8 ZLEV,ZPB,ZPT
1143 character(len=4) :: varLevel
1144
1145 ! loop over all body rows
1146 BODY: do index_body=1,obs_numbody(lobsSpaceData)
1147
1148 ! 1. Computation of sigmap
1149 ! ---------------------
1150
1151 IF ( (obs_bodyElem_i(lobsSpaceData,OBS_ASS,index_body) == obs_assimilated) .and. &
1152 (obs_bodyElem_i(lobsSpaceData,OBS_VNM,index_body).EQ. BUFR_NETT) ) THEN
1153
1154 IF ( (obs_bodyElem_i(lobsSpaceData,OBS_XTR,index_body) .NE. 0) .and. &
1155 (obs_bodyElem_i(lobsSpaceData,OBS_VCO,index_body) .EQ. 2) ) THEN
1156 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,index_body)
1157 varLevel = vnl_varLevelFromVarnum(ityp)
1158 IK=col_getNumLev(columnTrlOnAnlIncLev,varLevel)-1
1159 INDEX_HEADER = obs_bodyElem_i(lobsSpaceData,OBS_HIND,index_body)
1160 IPT = IK + col_getOffsetFromVarno(columnTrlOnAnlIncLev,ityp)
1161 IPB = IPT +1
1162 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body,col_getElem(column,IPB,INDEX_HEADER))
1163 ELSE
1164 INDEX_HEADER = obs_bodyElem_i(lobsSpaceData,OBS_HIND,index_body)
1165 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,index_body)
1166 IK = obs_bodyElem_i(lobsSpaceData,OBS_LYR,index_body)
1167 IPT = IK
1168 IPB = IPT+1
1169 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,index_body)
1170 varLevel = vnl_varLevelFromVarnum(ityp)
1171 ZPT = col_getPressure(columnTrlOnAnlIncLev,IK,INDEX_HEADER,varLevel)
1172 ZPB = col_getPressure(columnTrlOnAnlIncLev,IK+1,INDEX_HEADER,varLevel)
1173 ZWB = LOG(ZLEV/ZPT)/LOG(ZPB/ZPT)
1174 ZWT = 1.0D0 - ZWB
1175
1176 ! FIRST GUESS ERROR VARIANCE
1177
1178 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body, &
1179 (ZWB*col_getElem(column,IPB,INDEX_HEADER,'TT') + ZWT*col_getElem(column,IPT,INDEX_HEADER,'TT')))
1180 ENDIF
1181
1182 ENDIF
1183
1184 END DO BODY
1185
1186 end subroutine setfgett
1187
1188 !--------------------------------------------------------------------------
1189 ! setfgeSurf
1190 !--------------------------------------------------------------------------
1191 subroutine setfgeSurf(column, columnTrlOnAnlIncLev, lobsSpaceData)
1192 !
1193 !:Purpose: To interpolate vertically the contents of "column" to the
1194 ! pressure levels of the observations. A linear interpolation in
1195 ! ln(p) is performed.
1196 !
1197 implicit none
1198
1199 ! Arguments:
1200 type(struct_columnData), intent(in) :: column
1201 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1202 type(struct_obs) , intent(inout) :: lobsSpaceData
1203
1204 ! Locals:
1205 integer :: ipb, ipt, idim, headerIndex, ik, bodyIndex, ityp, bodyElem_i
1206 real(8) :: zwb, zwt, zlev, zpt, zpb, zhhh, bodyElem_r, colElem1, colElem2
1207 character(len=2) :: cfam
1208 character(len=4) :: varLevel
1209 character(len=12) :: stnid
1210 logical :: ok
1211
1212 ! loop over all body rows
1213 BODY: do bodyIndex = 1, obs_numbody( lobsSpaceData )
1214
1215 cfam = obs_getFamily( lobsSpaceData, bodyIndex_opt = bodyIndex )
1216 if( cfam == 'SF'.or. cfam == 'TM' .or. cfam == 'UA' .or. cfam == 'SC' .or. cfam == 'GP' .or. cfam == 'GL' ) then
1217
1218 ! Process all data within the domain of the model (excluding GB-GPS ZTD data)
1219 ok = .false.
1220
1221 if ( obs_bodyElem_i( lobsSpaceData, OBS_VCO, bodyIndex ) == 1 ) then
1222
1223 ityp = obs_bodyElem_i( lobsSpaceData, OBS_VNM, bodyIndex )
1224 if ( ityp == BUFR_NETS .or. ityp == BUFR_NEPS .or. ityp == BUFR_NEPN .or. ityp == BUFR_NESS .or. &
1225 ityp == BUFR_NEUS .or. ityp == BUFR_NEVS .or. ityp == BUFR_NEFS .or. ityp == BUFR_NEDS .or. &
1226 ityp == bufr_sst .or. ityp == BUFR_ICEC .or. ityp == bufr_logVis .or. ityp == bufr_gust .or. &
1227 ityp == bufr_riverFlow) then
1228
1229 ok = ( obs_bodyElem_i( lobsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated )
1230
1231 else if ( ityp == BUFR_NEZD ) then
1232
1233 ! make sure total zenith delay (from ground-based GPS) not treated
1234 ok=.false.
1235
1236 else
1237
1238 ok=(obs_bodyElem_i( lobsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated .and. &
1239 obs_bodyElem_i( lobsSpaceData, OBS_XTR, bodyIndex ) >= 0)
1240 if ( ok ) write(*,*) 'setfgesurf: WARNING!!! unknown obs seen'
1241 if ( ok ) write(*,*) 'setfgesurf: ityp=',ityp,', cfam=',cfam
1242
1243 end if
1244
1245 if ( ok ) then
1246
1247 headerIndex = obs_bodyElem_i( lobsSpaceData, OBS_HIND, bodyIndex )
1248 ityp = obs_bodyElem_i( lobsSpaceData, OBS_VNM , bodyIndex )
1249 varLevel = vnl_varLevelFromVarnum( ityp )
1250 idim = 1
1251 if ( varLevel == 'SF') idim = 0
1252 ik = obs_bodyElem_i( lobsSpaceData, OBS_LYR, bodyIndex )
1253 zlev = obs_bodyElem_r( lobsSpaceData, OBS_PPP, bodyIndex )
1254 zhhh = zlev
1255
1256 if ( ityp == BUFR_NETS .or. ityp == BUFR_NEPS .or. ityp == BUFR_NEPN .or. &
1257 ityp == BUFR_NESS .or. ityp == BUFR_NEUS .or. ityp == BUFR_NEVS .or. &
1258 ityp == bufr_logVis .or. ityp == bufr_gust) then
1259
1260 ipt = ik + col_getOffsetFromVarno(columnTrlOnAnlIncLev,ityp)
1261 ipb = ipt+1
1262 call obs_bodySet_r( lobsSpaceData, OBS_HPHT, bodyIndex, col_getElem( column, ipb, headerIndex ) )
1263
1264 else
1265
1266 ipt = ik + col_getOffsetFromVarno(columnTrlOnAnlIncLev,ityp)
1267 ipb = ipt+1
1268 zpt = col_getHeight(columnTrlOnAnlIncLev,ik,headerIndex,varLevel)
1269 zpb = col_getHeight(columnTrlOnAnlIncLev,ik+1,headerIndex,varLevel)
1270 zwb = idim*(zpt-zhhh)/(zpt-zpb)
1271 zwt = 1.d0 - zwb
1272
1273 if ( obs_bodyElem_i( lobsSpaceData, OBS_XTR, bodyIndex ) == 0 ) then
1274
1275 call obs_bodySet_r( lobsSpaceData, OBS_HPHT, bodyIndex, &
1276 zwb * col_getElem( column, ipb, headerIndex ) + zwt * col_getElem( column, ipt, headerIndex ))
1277
1278 else
1279
1280 call obs_bodySet_r( lobsSpaceData, OBS_HPHT, bodyIndex, &
1281 col_getElem( column, ik + col_getOffsetFromVarno( columnTrlOnAnlIncLev, ityp ), headerIndex ))
1282
1283 end if
1284
1285 stnid = obs_elem_c( lobsSpaceData, 'STID', headerIndex )
1286 if(stnid == '99999999' ) then
1287
1288 bodyElem_i = obs_bodyElem_i( lobsSpaceData, OBS_XTR, bodyIndex )
1289 write(*,*) 'setfgesurf: stn, ityp, xtr, ipt, ipb, zwt, zwb', &
1290 stnid, ityp, &
1291 bodyElem_i, ipt, ipb, zwt, zwb
1292 bodyElem_r = obs_bodyElem_i( lobsSpaceData, OBS_HPHT, bodyIndex )
1293 colElem1 = col_getElem( column, ipb, headerIndex )
1294 colElem2 = col_getElem( column, ipt, headerIndex )
1295 write(*,*) 'setfgesurf: gobs(ipb), gobs(ipt), fge', &
1296 colElem1, colElem2, bodyElem_r
1297
1298 endif
1299
1300 end if
1301 end if
1302 end if
1303
1304 end if
1305
1306 end do BODY
1307
1308 end subroutine setfgeSurf
1309
1310 !--------------------------------------------------------------------------
1311 ! setfgedif
1312 !--------------------------------------------------------------------------
1313 subroutine setfgedif(cdfam,columnTrlOnAnlIncLev,lobsSpaceData)
1314 !
1315 !:Purpose: To construct the FIRST GUESS ERROR VARIANCES from the
1316 ! diff-calculated dependencies and the primary errors.
1317 !
1318 implicit none
1319
1320 ! Arguments:
1321 character(len=2), intent(in) :: cdfam
1322 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1323 type(struct_obs), intent(inout) :: lobsSpaceData
1324
1325 ! Locals:
1326 INTEGER INDEX_HEADER, IDATYP, INDEX_BODY, iProfile, varNum
1327 REAL*8 zLat, Lat, sLat
1328 REAL*8 zLon, Lon
1329 REAL*8 zAzm !, Azm
1330 INTEGER ISAT
1331 REAL*8 Rad, Geo
1332 REAL*8, allocatable :: zPP(:)
1333 REAL*8, allocatable :: zDP(:)
1334 REAL*8, allocatable :: zTT(:)
1335 REAL*8, allocatable :: zHU(:)
1336 REAL*8, allocatable :: zHT(:)
1337 REAL*8, allocatable :: zUU(:)
1338 REAL*8, allocatable :: zVV(:)
1339 INTEGER JL
1340 REAL*8 ZP0, ZMT
1341 REAL*8 ZFGE, ZERR
1342 INTEGER JV, NGPSLEV, NWNDLEV
1343 LOGICAL ASSIM, LFIRST, FIRSTHEADER
1344 INTEGER NH, NH1
1345 REAL*8 DV (gps_ncvmx)
1346 TYPE(GPS_PROFILE) :: PRF
1347 REAL*8 , allocatable :: H (:),AZMV(:)
1348 TYPE(GPS_DIFF), allocatable :: RSTV(:)
1349 type(struct_vco), pointer :: vco_anl
1350 real*8, dimension(:), pointer :: dPdPs
1351
1352 WRITE(*,*)'ENTER SETFGEDIFF'
1353
1354 ! Initializations
1355
1356 nullify(dPdPs)
1357
1358 NGPSLEV=col_getNumLev(columnTrlOnAnlIncLev,'TH')
1359 NWNDLEV=col_getNumLev(columnTrlOnAnlIncLev,'MM')
1360 LFIRST=.FALSE.
1361 if ( .NOT.allocated(ose_vRO_Jacobian) ) then
1362 LFIRST = .TRUE.
1363 allocate(zPP (NGPSLEV))
1364 allocate(zDP (NGPSLEV))
1365 allocate(zTT (NGPSLEV))
1366 allocate(zHU (NGPSLEV))
1367 allocate(zHT (NGPSLEV))
1368 allocate(zUU (NGPSLEV))
1369 allocate(zVV (NGPSLEV))
1370
1371 allocate(ose_vRO_Jacobian(gps_numROProfiles,gps_RO_MAXPRFSIZE,2*NGPSLEV+1))
1372 allocate(ose_vRO_lJac (gps_numROProfiles))
1373 ose_vRO_Jacobian = 0.d0
1374 ose_vRO_lJac = .False.
1375
1376 allocate( H (gps_RO_MAXPRFSIZE) )
1377 allocate( AZMV (gps_RO_MAXPRFSIZE) )
1378 allocate( RSTV (gps_RO_MAXPRFSIZE) )
1379 endif
1380
1381 vco_anl => col_getVco(columnTrlOnAnlIncLev)
1382
1383 ! Loop over all header indices of the 'RO' family:
1384
1385 call obs_set_current_header_list(lobsSpaceData,CDFAM)
1386 FIRSTHEADER=.TRUE.
1387 HEADER: do
1388 index_header = obs_getHeaderIndex(lobsSpaceData)
1389 if (index_header < 0) exit HEADER
1390
1391 ! Process only refractivity data (codtyp 169)
1392
1393 IDATYP = obs_headElem_i(lobsSpaceData,OBS_ITY,INDEX_HEADER)
1394 IF ( IDATYP .EQ. 169 ) THEN
1395 iProfile = gps_iprofile_from_index(index_header)
1396 varNum = gps_vRO_IndexPrf(iProfile, 2)
1397
1398 ! Scan for requested data values of the profile, and count them
1399
1400 ASSIM = .FALSE.
1401 NH = 0
1402 call obs_set_current_body_list(lobsSpaceData, INDEX_HEADER)
1403 BODY: do
1404 INDEX_BODY = obs_getBodyIndex(lobsSpaceData)
1405 if (INDEX_BODY < 0) exit BODY
1406 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated ) THEN
1407 ASSIM = .TRUE.
1408 NH = NH + 1
1409 ENDIF
1410 ENDDO BODY
1411
1412 ! If assimilations are requested, prepare and apply the observation operator
1413
1414 IF (ASSIM) THEN
1415 iProfile=gps_iprofile_from_index(INDEX_HEADER)
1416
1417 ! Profile at the observation location:
1418
1419 if (.not.ose_vRO_lJac(iProfile)) then
1420
1421 ! Basic geometric variables of the profile:
1422
1423 zLat = obs_headElem_r(lobsSpaceData,OBS_LAT,INDEX_HEADER)
1424 zLon = obs_headElem_r(lobsSpaceData,OBS_LON,INDEX_HEADER)
1425 ISAT = obs_headElem_i(lobsSpaceData,OBS_SAT,INDEX_HEADER)
1426 Rad = obs_headElem_r(lobsSpaceData,OBS_TRAD,INDEX_HEADER)
1427 Geo = obs_headElem_r(lobsSpaceData,OBS_GEOI,INDEX_HEADER)
1428 zAzm = obs_headElem_r(lobsSpaceData,OBS_AZA,INDEX_HEADER) / MPC_DEGREES_PER_RADIAN_R8
1429 zMT = col_getHeight(columnTrlOnAnlIncLev,NGPSLEV,INDEX_HEADER,'TH')
1430 Lat = zLat * MPC_DEGREES_PER_RADIAN_R8
1431 Lon = zLon * MPC_DEGREES_PER_RADIAN_R8
1432 !Azm = zAzm * MPC_DEGREES_PER_RADIAN_R8
1433 sLat = sin(zLat)
1434 zMT = zMT * ec_rg / gps_gravitysrf(sLat)
1435 zP0 = col_getElem(columnTrlOnAnlIncLev,1,INDEX_HEADER,'P0')
1436
1437 ! approximation for dPdPs
1438 if (associated(dPdPs)) then
1439 deallocate(dPdPs)
1440 end if
1441 call czp_fetch1DdPdPs(vco_anl, zP0, profT_opt=dPdPs)
1442 zDP(1:NGPSLEV) = dPdPs(1:NGPSLEV)
1443
1444 DO JL = 1, NGPSLEV
1445
1446 ! Profile x
1447
1448 zPP(JL) = col_getPressure(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1449 zTT(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TT') - MPC_K_C_DEGREE_OFFSET_R8
1450 zHU(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'HU')
1451 zHT(JL) = col_getHeight(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1452 zUU(JL) = 0.d0
1453 zVV(JL) = 0.d0
1454 ENDDO
1455 DO JL = 1, NWNDLEV
1456 zUU(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'UU')
1457 zVV(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'VV')
1458 ENDDO
1459 zUU(NGPSLEV) = zUU(NWNDLEV)
1460 zVV(NGPSLEV) = zUU(NWNDLEV)
1461
1462 ! GPS profile structure:
1463
1464 call gps_struct1sw(ngpslev,zLat,zLon,zAzm,zMT,Rad,geo,zP0,zPP,zDP,zTT,zHU,zUU,zVV,prf)
1465 !call gps_struct1sw_v2(ngpslev,zLat,zLon,zAzm,zMT,Rad,geo,zP0,zPP,zDP,zTT,zHU,zUU,zVV,zHT,prf)
1466
1467 ! Prepare the vector of all the observations:
1468
1469 NH1 = 0
1470 call obs_set_current_body_list(lobsSpaceData, INDEX_HEADER)
1471 BODY_2: do
1472 INDEX_BODY = obs_getBodyIndex(lobsSpaceData)
1473 if (INDEX_BODY < 0) exit BODY_2
1474 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated ) THEN
1475 NH1 = NH1 + 1
1476 H(NH1) = obs_bodyElem_r(lobsSpaceData,OBS_PPP,INDEX_BODY)
1477 AZMV(NH1)= zAzm
1478 ENDIF
1479 ENDDO BODY_2
1480
1481 ! Apply the observation operator:
1482
1483 IF (varNum == bufr_nebd) THEN
1484 CALL GPS_BNDOPV1(H, AZMV, NH, PRF, RSTV)
1485 ELSE
1486 CALL GPS_REFOPV (H, NH, PRF, RSTV)
1487 ENDIF
1488 DO NH1=1,NH
1489 ose_vRO_Jacobian(iProfile,NH1,:)= RSTV(NH1)%DVAR(1:2*NGPSLEV+1)
1490 ENDDO
1491 ose_vRO_lJac(iProfile) = .True.
1492 endif
1493
1494 ! Local error
1495
1496 DO JL = 1, NGPSLEV
1497 DV ( JL) = 1.d0
1498 DV (NGPSLEV+JL) = 1.d0
1499 ENDDO
1500 DV (2*NGPSLEV+1) = 2.d0
1501
1502 ! Perform the H(xb)DV operation:
1503
1504 NH1 = 0
1505 call obs_set_current_body_list(lobsSpaceData, INDEX_HEADER)
1506 BODY_3: do
1507 INDEX_BODY = obs_getBodyIndex(lobsSpaceData)
1508 if (INDEX_BODY < 0) exit BODY_3
1509 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated ) THEN
1510 NH1 = NH1 + 1
1511
1512 ! Observation jacobian
1513
1514 ! JAC = RSTV(NH1)%DVAR
1515
1516 ! Evaluate sqrt( H(xb)DV **2 )
1517
1518 ZFGE = 0.d0
1519 DO JV = 1, 2*PRF%NGPSLEV+1
1520 ZFGE = ZFGE + (ose_vRO_Jacobian(iProfile,NH1,JV) * DV(JV))**2
1521 ENDDO
1522 ZFGE = SQRT(ZFGE)
1523 ZERR = obs_bodyElem_r(lobsSpaceData,OBS_OER,INDEX_BODY)
1524
1525 ! FIRST GUESS ERROR VARIANCE
1526
1527 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,INDEX_BODY,ZFGE)
1528 IF (FIRSTHEADER) THEN
152911 FORMAT(A12,2I5,F12.2,3F16.8)
1530 WRITE(*,11)'SETFGEDIFFGE',NH1,NH,H(NH1),RSTV(NH1)%VAR,ZFGE,ZERR
1531 ENDIF
1532 ENDIF
1533 ENDDO BODY_3
1534 ENDIF
1535 ENDIF
1536 FIRSTHEADER = .FALSE.
1537 ENDDO HEADER
1538
1539 IF (LFIRST) THEN
1540 deallocate( RSTV )
1541 deallocate( AZMV )
1542 deallocate( H )
1543
1544 deallocate(zVV)
1545 deallocate(zUU)
1546 deallocate(zHT)
1547 deallocate(zHU)
1548 deallocate(zTT)
1549 deallocate(zDP)
1550 deallocate(zPP)
1551 ENDIF
1552
1553 WRITE(*,*)'EXIT SETFGEDIFF'
1554 end subroutine setfgedif
1555
1556 !--------------------------------------------------------------------------
1557 ! setfgegps
1558 !--------------------------------------------------------------------------
1559 subroutine setfgegps(column,columnTrlOnAnlIncLev,lobsSpaceData)
1560 !
1561 !:Purpose: To set FGE for all GPS ZTD observations using Jacobians from ZTD
1562 ! observation operator
1563 !
1564 !:Option: Test ZTD operators (compares H(x+dx)-H(x) with (dH/dx)*dx
1565 ! when gps_gb_LTESTOP = .true.)
1566 !
1567 !:Note:
1568 ! _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
1569 ! 9 October 2015
1570 !
1571 ! :NOTE: Effective Rev644M, this routine is no longer used!
1572 ! FGE for ZTD is no longer needed for background check.
1573 ! Routine is called only when gps_gb_LTESTOP=.true., in which
1574 ! case the operator test only is done.
1575 !
1576 ! _/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/
1577 !
1578 implicit none
1579
1580 ! Arguments:
1581 type(struct_columnData), intent(in) :: column
1582 type(struct_columnData), intent(in) :: columnTrlOnAnlIncLev
1583 type(struct_obs), intent(inout) :: lobsSpaceData
1584
1585 ! Locals:
1586 ! column contains background errors for control variables on model levels
1587 ! columnTrlOnAnlIncLev contains lo-res first guess profiles at obs locations
1588 type(struct_vco), pointer :: vco_anl
1589 REAL*8 ZLAT, Lat
1590 REAL*8 ZLON, Lon
1591 REAL*8, allocatable :: ZPP(:)
1592 REAL*8, allocatable :: ZDP(:)
1593 REAL*8, allocatable :: ZTT(:)
1594 REAL*8, allocatable :: ZHU(:)
1595 REAL*8, allocatable :: zHeight(:)
1596 REAL*8, allocatable :: zHeight2(:)
1597 REAL*8, allocatable :: ZTTB(:)
1598 REAL*8, allocatable :: ZHUB(:)
1599 REAL*8, allocatable :: ZQQB(:)
1600 REAL*8, allocatable :: ZQQ(:)
1601 REAL*8, allocatable :: ZTTB_P(:)
1602 REAL*8, allocatable :: ZQQB_P(:)
1603 REAL*8, allocatable :: zHeight_P(:)
1604 REAL*8, allocatable :: ZPP_P(:)
1605 REAL*8 ZP0
1606 REAL*8 ZP0B, ZP0B_P
1607 REAL*8 ZMT, ZTOP, ZBOT
1608 REAL*8 JAC(gps_ncvmx)
1609 REAL*8 DX (gps_ncvmx)
1610 REAL*8 ZLEV, ZTDOBS, ZPSMOD
1611 REAL*8 ZLSUM
1612 REAL*8 DELTAH_NL, DELTAH_TL
1613 REAL*8 PERTFAC, ZTDM
1614 REAL*8 ZDZMIN, ZSUMTEST
1615 INTEGER INDEX_HEADER, FIRST_HEADER
1616 INTEGER IDATYP, ITYP
1617 INTEGER IDATA, IDATEND, INDEX_BODY
1618 INTEGER JL, NFLEV_T, ILYR, IOBS
1619 INTEGER INOBS_OPT, INOBS_JAC, icount, iversion
1620 LOGICAL ASSIM, OK, LSTAG
1621 CHARACTER*9 STN_JAC
1622 character(len=12) :: stnid
1623 CHARACTER(len=4) :: varLevel
1624 TYPE(GPS_PROFILEZD) :: PRF, PRFP
1625 TYPE(GPS_DIFF) :: ZTDopv, ZTDopvP
1626 real*8, dimension(:), pointer :: dPdPs
1627
1628 IF (gps_gb_numZTD .EQ. 0) RETURN
1629
1630 ! Initializations
1631
1632 nullify(dPdPs)
1633
1634 NFLEV_T = col_getNumLev(columnTrlOnAnlIncLev,'TH')
1635 allocate(ZPP(NFLEV_T))
1636 allocate(ZDP(NFLEV_T))
1637 allocate(ZTT(NFLEV_T))
1638 allocate(ZHU(NFLEV_T))
1639 allocate(zHeight(NFLEV_T))
1640 allocate(ZTTB(NFLEV_T))
1641 allocate(ZHUB(NFLEV_T))
1642 allocate(ZQQB(NFLEV_T))
1643 allocate(ZQQ(NFLEV_T))
1644
1645 ! Number of locations/sites for observation operator test
1646 INOBS_OPT = 50
1647 ! Number of locations/sites for Jacobian printout
1648 INOBS_JAC = 5
1649 ! Factor to multiply background errors for perturbation vector
1650 PERTFAC = 0.75d0
1651
1652 STN_JAC = 'FSL_BRFT '
1653
1654 ZDZMIN = gps_gb_DZMIN
1655
1656 vco_anl => col_getVco(columnTrlOnAnlIncLev)
1657 iversion = vco_anl%Vcode
1658 if (iversion .eq. 5002) then
1659 LSTAG = .TRUE.
1660 WRITE(*,*)'VERTICAL COORD OF ANALYSIS FIELDS IS STAGGERED'
1661 WRITE(*,*)'VCODE= ',iversion,' LSTAG= ',LSTAG
1662 else
1663 LSTAG = .FALSE.
1664 WRITE(*,*)'VERTICAL COORD OF ANALYSIS FIELDS IS NOT STAGGERED'
1665 WRITE(*,*)'VCODE= ',iversion,' LSTAG= ',LSTAG
1666 endif
1667
1668 IF ( .NOT.gps_gb_LTESTOP ) THEN
1669
1670 first_header=-1
1671 icount = 0
1672
1673 ! loop over all header indices of the 'GP' family
1674 call obs_set_current_header_list(lobsSpaceData,'GP')
1675 HEADER: do
1676 index_header = obs_getHeaderIndex(lobsSpaceData)
1677 if (index_header < 0) exit HEADER
1678 if (first_header .eq. -1) first_header = index_header
1679
1680 ! Process only zenith delay data (codtyp 189 and BUFR_NEZD)
1681
1682 IDATYP = obs_headElem_i(lobsSpaceData,OBS_ITY,INDEX_HEADER)
1683 IF ( IDATYP .EQ. 189 ) THEN
1684
1685 ! Loop over data in the observations
1686
1687 IDATA = obs_headElem_i(lobsSpaceData,OBS_RLN,INDEX_HEADER)
1688 IDATEND = obs_headElem_i(lobsSpaceData,OBS_NLV,INDEX_HEADER) + IDATA - 1
1689 ASSIM = .FALSE.
1690
1691 ! Scan for requested assimilations, and count them.
1692
1693 DO INDEX_BODY= IDATA, IDATEND
1694 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,INDEX_BODY)
1695 OK = ( (ITYP .EQ. BUFR_NEZD) .AND. (obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated) )
1696 IF ( OK ) THEN
1697 ASSIM = .TRUE.
1698 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,INDEX_BODY)
1699 icount = icount + 1
1700 ENDIF
1701 ENDDO
1702
1703 ! If assimilations are requested, apply the AD observation operator
1704
1705 IF (ASSIM) THEN
1706
1707 ! LR background profile and background errors at the observation location x :
1708
1709 Lat = obs_headElem_r(lobsSpaceData,OBS_LAT,INDEX_HEADER)
1710 Lon = obs_headElem_r(lobsSpaceData,OBS_LON,INDEX_HEADER)
1711 ZLAT = Lat * MPC_DEGREES_PER_RADIAN_R8
1712 ZLON = Lon * MPC_DEGREES_PER_RADIAN_R8
1713 ZP0B = col_getElem(columnTrlOnAnlIncLev,1,INDEX_HEADER,'P0')
1714 DO JL = 1, NFLEV_T
1715 ZPP(JL) = col_getPressure(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1716 ZTTB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TT')- 273.15d0
1717 ZTT(JL) = col_getElem(column,JL,INDEX_HEADER,'TT')
1718 DX(JL) = ZTT(JL)
1719 ZHUB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'HU')
1720 ZQQB(JL) = ZHUB(JL)
1721 ZHU(JL) = col_getElem(column,JL,INDEX_HEADER,'HU')
1722 DX(NFLEV_T+JL) = ZHU(JL)
1723 zHeight(JL) = col_getHeight(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1724 DX(2*NFLEV_T+JL) = col_getHeight(column,JL,INDEX_HEADER,'TH')
1725 ENDDO
1726 ZP0 = col_getElem(column,1,INDEX_HEADER,'P0')
1727 DX(3*NFLEV_T+1) = ZP0
1728 ZMT = zHeight(NFLEV_T)
1729 CALL gps_structztd_v2(NFLEV_T,Lat,Lon,ZMT,ZP0B,ZPP,ZTTB,ZHUB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF)
1730 CALL gps_ztdopv(ZLEV,PRF,gps_gb_LBEVIS,ZDZMIN,ZTDopv,ZPSMOD,gps_gb_IZTDOP)
1731 JAC = ZTDopv%DVar
1732
1733 DO INDEX_BODY= IDATA, IDATEND
1734 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,INDEX_BODY)
1735 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated .AND. ITYP.EQ.BUFR_NEZD ) THEN
1736
1737 ! Observation error SDERR
1738 ! ZOER = obs_bodyElem_r(lobsSpaceData,OBS_PPP,INDEX_BODY)
1739
1740 ! Observation height (m)
1741 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,INDEX_BODY)
1742
1743 ZLSUM = 0.0d0
1744
1745 DO JL = 1, 3*NFLEV_T+1
1746 ZLSUM = ZLSUM + (JAC(JL)*DX(JL))**2
1747 ENDDO
1748 call obs_bodySet_r(lobsSpaceData,OBS_HPHT,index_body,SQRT(ZLSUM))
1749
1750 IF (icount .LE. INOBS_JAC) THEN
1751 ! IF ( obs_elem_c(lobsSpaceData,'STID',INDEX_HEADER) .EQ. STN_JAC ) THEN
1752 stnid = obs_elem_c(lobsSpaceData,'STID',INDEX_HEADER)
1753 WRITE(*,'(A11,A9)') 'SETFGEGPS: ', stnid
1754 WRITE(*,*) ' ZTD, ZTD FGE = ', ZTDopv%Var, SQRT(ZLSUM)
1755 WRITE(*,'(A11,A9,3(1x,f7.2))') &
1756 'SETFGEGPS: ',stnid,ZLAT,ZLON,ZLEV
1757 WRITE(*,*) 'JL JACT JACQ FGE_T FGE_LQ QQ'
1758 DO JL = 1, NFLEV_T
1759 WRITE(*,'(1X,I2,5(1x,E13.6))') JL,JAC(JL),JAC(JL+NFLEV_T)/ZQQB(JL),ZTT(JL),ZHU(JL),ZQQB(JL)
1760 ENDDO
1761 WRITE(*,*) 'JACPS FGE_PS'
1762 WRITE(*,'(2(1x,E13.6))') JAC(3*NFLEV_T+1), ZP0
1763 ENDIF
1764
1765 ENDIF
1766 ENDDO
1767 ENDIF
1768 ENDIF
1769
1770 ENDDO HEADER
1771
1772 ENDIF
1773
1774 !-------------------------------------------------------------------------
1775
1776 IF ( gps_gb_LTESTOP ) THEN
1777
1778 allocate(ZTTB_P(NFLEV_T))
1779 allocate(ZQQB_P(NFLEV_T))
1780 allocate(zHeight2(NFLEV_T))
1781 allocate(zHeight_P(NFLEV_T))
1782 allocate(ZPP_P(NFLEV_T))
1783
1784 icount = 0
1785 ZSUMTEST = 0
1786
1787 ! loop over all header indices of the 'GP' family
1788 call obs_set_current_header_list(lobsSpaceData,'GP')
1789 HEADER2: DO
1790 index_header = obs_getHeaderIndex(lobsSpaceData)
1791 if (index_header < 0) exit HEADER2
1792 if (icount > INOBS_OPT ) exit HEADER2
1793
1794 ! Loop over data in the observations
1795
1796 IDATA = obs_headElem_i(lobsSpaceData,OBS_RLN,INDEX_HEADER)
1797 IDATEND = obs_headElem_i(lobsSpaceData,OBS_NLV,INDEX_HEADER) + IDATA - 1
1798
1799 ! LR background profile and background errors at the observation location x :
1800
1801 Lat = obs_headElem_r(lobsSpaceData,OBS_LAT,INDEX_HEADER)
1802 Lon = obs_headElem_r(lobsSpaceData,OBS_LON,INDEX_HEADER)
1803 ZLAT = Lat * MPC_DEGREES_PER_RADIAN_R8
1804 ZLON = Lon * MPC_DEGREES_PER_RADIAN_R8
1805 ZP0B = col_getElem(columnTrlOnAnlIncLev,1,INDEX_HEADER,'P0')
1806
1807 ! approximation for dPdPs
1808 if (associated(dPdPs)) then
1809 deallocate(dPdPs)
1810 end if
1811 call czp_fetch1DdPdPs(vco_anl, ZP0B, profT_opt=dPdPs)
1812 zDP(1:NFLEV_T) = dPdPs(1:NFLEV_T)
1813
1814 DO JL = 1, NFLEV_T
1815 ZPP(JL) = col_getPressure(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1816 ZTTB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TT')- 273.15d0
1817 ZTT(JL) = col_getElem(column,JL,INDEX_HEADER,'TT') * PERTFAC
1818 ZQQB(JL) = col_getElem(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'HU')
1819 ZQQ(JL) = col_getElem(column,JL,INDEX_HEADER,'HU') * PERTFAC
1820 zHeight(JL) = col_getHeight(columnTrlOnAnlIncLev,JL,INDEX_HEADER,'TH')
1821 zHeight2(JL) = col_getHeight(column,JL,INDEX_HEADER,'TH') * PERTFAC
1822 ENDDO
1823 ZP0 = col_getElem(column,1,INDEX_HEADER,'P0') * PERTFAC
1824 ZMT = zHeight(NFLEV_T)
1825
1826 DO JL = 1, NFLEV_T
1827 DX ( JL) = ZTT(JL)
1828 DX (NFLEV_T+JL) = ZQQ(JL)
1829 DX (2*NFLEV_T+JL) = zHeight2(JL)
1830 ENDDO
1831 DX (3*NFLEV_T+1) = ZP0
1832
1833 ZTDOBS = -1.0d0
1834 DO INDEX_BODY = IDATA, IDATEND
1835 ITYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,INDEX_BODY)
1836 IOBS = obs_bodyElem_i(lobsSpaceData,OBS_HIND,INDEX_BODY)
1837 IF ( obs_bodyElem_i(lobsSpaceData,OBS_ASS,INDEX_BODY) == obs_assimilated .AND. ITYP .EQ. BUFR_NEZD ) THEN
1838 varLevel = vnl_varLevelFromVarnum(ITYP)
1839 ZTDOBS = obs_bodyElem_r(lobsSpaceData,OBS_VAR,INDEX_BODY)
1840 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,INDEX_BODY)
1841 ILYR = obs_bodyElem_i(lobsSpaceData,OBS_LYR,INDEX_BODY)
1842 ZTOP = col_getHeight(columnTrlOnAnlIncLev,ILYR,IOBS,varLevel)
1843 if ( ILYR .LT. NFLEV_T ) then
1844 ZBOT = col_getHeight(columnTrlOnAnlIncLev,ILYR+1,IOBS,varLevel)
1845 else
1846 ZBOT = ZTOP
1847 endif
1848 icount = icount + 1
1849 ENDIF
1850 ENDDO
1851
1852 IF ( ZTDOBS .GT. 0.d0 ) THEN
1853 ! Create the pertubation control vector
1854 DO JL = 1, NFLEV_T
1855 ZPP_P(JL) = ZPP(JL) + ZDP(JL)*ZP0
1856 ZTTB_P(JL) = ZTTB(JL) + ZTT(JL)
1857 ZQQB_P(JL) = ZQQB(JL) + ZQQ(JL)
1858 zHeight_P(JL) = zHeight(JL) + zHeight2(JL)
1859 ENDDO
1860 ZP0B_P = ZP0B + ZP0
1861
1862 ! Non-linear observation operator --> delta_H = H(x+delta_x) - H(x)
1863
1864 CALL gps_structztd_v2(NFLEV_T,Lat,Lon,ZMT,ZP0B,ZPP,ZTTB,ZQQB,zHeight,gps_gb_LBEVIS,gps_gb_IREFOPT,PRF)
1865 CALL gps_structztd_v2(NFLEV_T,Lat,Lon,ZMT,ZP0B_P,ZPP_P,ZTTB_P,ZQQB_P,zHeight_P,gps_gb_LBEVIS,gps_gb_IREFOPT,PRFP)
1866 CALL gps_ztdopv(ZLEV,PRF,gps_gb_LBEVIS,ZDZMIN,ZTDopv,ZPSMOD,gps_gb_IZTDOP)
1867 JAC = ZTDopv%DVar
1868 ZTDM = ZTDopv%Var
1869 CALL gps_ztdopv(ZLEV,PRFP,gps_gb_LBEVIS,ZDZMIN,ZTDopvP,ZPSMOD,gps_gb_IZTDOP)
1870 DELTAH_NL = ZTDopvP%Var - ZTDopv%Var
1871
1872 ! Linear --> delta_H = dH/dx * delta_x
1873
1874 DELTAH_TL = 0.0d0
1875 DO JL = 1, 3*NFLEV_T+1
1876 DELTAH_TL = DELTAH_TL + JAC(JL)*DX(JL)
1877 ENDDO
1878
1879 stnid = obs_elem_c(lobsSpaceData,'STID',INDEX_HEADER)
1880 WRITE(*,*) 'SETFGEGPS: GPS ZTD OBSOP TEST FOR SITE ', stnid
1881 WRITE(*,*) ' '
1882 WRITE(*,*) ' DZ (M), MODEL LEVEL ABOVE = ', ZLEV-ZMT, ILYR
1883 WRITE(*,*) ' ZLEV (M), ZTOP (M), ZBOT (M) = ', ZLEV, ZTOP, ZBOT
1884 WRITE(*,*) ' ZTD OBS (MM) = ', ZTDOBS*1000.d0
1885 WRITE(*,*) ' ZTD_MOD = ', ZTDM*1000.d0
1886 WRITE(*,*) ' DELTAH_NL, DELTAH_TL = ', DELTAH_NL*1000.d0, DELTAH_TL*1000.d0
1887 WRITE(*,*) ' '
1888 WRITE(*,*) ' DELTAH_TL/DELTAH_NL = ', DELTAH_TL/DELTAH_NL
1889 WRITE(*,*) ' '
1890
1891 ZSUMTEST = ZSUMTEST + (DELTAH_TL/DELTAH_NL)
1892
1893 ENDIF
1894
1895 ENDDO HEADER2
1896
1897 WRITE(*,*) ' '
1898 WRITE(*,*) 'SETFGEGPS: ----- GPS ZTD OBSOP TEST SUMMARY -----'
1899 WRITE(*,*) ' NUMBER OF TESTS (sites) = ', icount
1900 WRITE(*,*) ' AVG DELTAH_TL/DELTAH_NL = ', ZSUMTEST/FLOAT(icount)
1901 WRITE(*,*) ' '
1902
1903 deallocate(ZTTB_P)
1904 deallocate(ZQQB_P)
1905 deallocate(zHeight2)
1906 deallocate(zHeight_P)
1907 deallocate(ZPP_P)
1908
1909 ENDIF
1910 !-------------------------------------------------------------------------
1911
1912 deallocate(ZPP)
1913 deallocate(ZDP)
1914 deallocate(ZTT)
1915 deallocate(ZHU)
1916 deallocate(zHeight)
1917 deallocate(ZTTB)
1918 deallocate(ZHUB)
1919 deallocate(ZQQB)
1920 deallocate(ZQQ)
1921
1922 end subroutine setfgegps
1923
1924 !------------------ CH obs family OmP error std dev routines --------------
1925
1926 !--------------------------------------------------------------------------
1927 ! ose_setOmPstddevCH
1928 !--------------------------------------------------------------------------
1929 character(len=4) function ose_setOmPstddevCH(obsSpaceData) result(availableOMPE)
1930 !
1931 !:Purpose: To read OmP error std dev from auxiliary file or calculate from OmP.
1932 !
1933 implicit none
1934
1935 ! Arguments:
1936 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data; output saved in OBS_OMPE column
1937
1938 ! Locals:
1939 integer, parameter :: ndim=1
1940
1941 ! Check for the presence of CH observations
1942 availableOMPE = ' '
1943 if ( .not.obs_famExist(obsSpaceData,'CH',localMPI_opt=.true.) ) return
1944
1945 ! read the OmP error std. dev. information from the auxiliary file
1946 call ose_readOmPstddev_auxfileCH
1947
1948 availableOMPE='None'
1949 if ( OmPstdCH%n_stnid == 0 ) return ! All CH family OBS_OMPE to be estimated elsewhere via use of OBS_HPHT and OBS_OER
1950
1951 ! Calc from the OmP dataset if requested (and there is enough data)
1952 if ( any(OmPstdCH%source(1:OmPstdCH%n_stnid) == 1) ) call ose_calcOmPstddevCH(obsSpaceData)
1953
1954 ! Assign OmP error std dev values to OBS_OMPE for CH family
1955 if ( any(OmPstdCH%source(1:OmPstdCH%n_stnid) <= 1) ) call ose_fillOmPstddevCH(obsSpaceData)
1956
1957 ! Deallocate OmPstdCH space
1958 call ose_deallocOmPstddevCH
1959
1960 ! Check if ALL CH family obs have been assigned usable values in OBS_OMPE
1961 if ( .not. ose_OmPstddevExistsForAllCH(ObsSpaceData) ) then
1962 ! For some obs, need to calc OBS_HPHT to get OBS_OMPE
1963 availableOMPE='Some'
1964 else
1965 ! All OBS_OMPE are available
1966 availableOMPE='All '
1967 end if
1968
1969 end function ose_setOmPstddevCH
1970
1971 !--------------------------------------------------------------------------
1972 ! ose_readOmPstddev_auxfileCH
1973 !--------------------------------------------------------------------------
1974 subroutine ose_readOmPstddev_auxfileCH
1975 !
1976 !:Purpose: To read and store OmP error std. dev. as needed for CH
1977 ! family obs - if/when available.
1978 !
1979 implicit none
1980
1981 ! Locals:
1982 integer, external :: fnom, fclos
1983 integer :: ierr, nulstat
1984 logical :: LnewExists
1985 character (len=128) :: ligne
1986 character(len=11) :: AuxObsDataFileCH = 'obsinfo_chm'
1987 integer :: ipos, stnidIndex, monthIndex, levIndex, ios, isize, icount
1988 character(len=20) :: abortText
1989
1990 ! Initialization
1991
1992 OmPstdCH%n_stnid=0
1993
1994 ! Check the existence of the text file with statistics
1995
1996 INQUIRE(FILE=trim(AuxObsDataFileCH),EXIST=LnewExists)
1997 if (.not.LnewExists) then
1998 WRITE(*,*) '-----------------------------------------------------------------------------'
1999 WRITE(*,*) 'WARNING! ose_readOmPstddev_auxfileCH: auxiliary file ' // trim(AuxObsDataFileCH)
2000 WRITE(*,*) 'WARNING! not available. Default CH family OmP stddev to be applied if needed.'
2001 WRITE(*,*) '-----------------------------------------------------------------------------'
2002 return
2003 end if
2004
2005 ! Read error std dev. from file auxiliary file for constituent data
2006
2007 nulstat=0
2008 ierr=fnom(nulstat,trim(AuxObsDataFileCH),'SEQ',0)
2009 if ( ierr == 0 ) then
2010 open(unit=nulstat, file=trim(AuxObsDataFileCH), status='OLD')
2011 else
2012 call utl_abort('ose_readOmPstddev_auxfileCH: COULD NOT OPEN AUXILIARY FILE ' // trim(AuxObsDataFileCH) )
2013 end if
2014
2015 ! Read OmP error standard deviations for constituents or related directives if available.
2016
2017 ios=0
2018 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
2019 do while (trim(adjustl(ligne(1:12))) /= 'SECTION V:')
2020 read(nulstat,'(A)',iostat=ios,err=10,end=15) ligne
2021 end do
2022
2023 ! Read number of observation set sub-families (STNIDs and ...) and allocate space
2024
2025 read(nulstat,*,iostat=ios,err=10,end=10) OmPstdCH%n_stnid
2026 read(nulstat,*,iostat=ios,err=10,end=10) isize
2027
2028 allocate(OmPstdCH%stnids(OmPstdCH%n_stnid),OmPstdCH%std_type(OmPstdCH%n_stnid))
2029 allocate(OmPstdCH%n_month(OmPstdCH%n_stnid),OmPstdCH%n_lat(OmPstdCH%n_stnid))
2030 allocate(OmPstdCH%source(OmPstdCH%n_stnid),OmPstdCH%ibegin(OmPstdCH%n_stnid))
2031 allocate(OmPstdCH%element(OmPstdCH%n_stnid),OmPstdCH%n_lvl(OmPstdCH%n_stnid))
2032 allocate(OmPstdCH%std(isize))
2033 allocate(OmPstdCH%levels(isize),OmPstdCH%lat(isize),OmPstdCH%month(isize))
2034
2035 OmPstdCH%element(:)=0
2036 OmPstdCH%source(:)=0
2037 OmPstdCH%std_type(:)=0
2038 OmPstdCH%n_lvl(:)=1
2039 OmPstdCH%n_lat(:)=1
2040 OmPstdCH%n_month(:)=1
2041
2042 ! Begin reading for each sub-family
2043 ! Important: Combination of STNID, BUFR element and number of vertical levels
2044 ! to determine association to the observations.
2045
2046 icount=0
2047 STNIDLOOP: do stnidIndex=1,OmPstdCH%n_stnid
2048
2049 ! disregard line of dashes
2050 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
2051
2052 ! Read STNID (* as wildcard)
2053 read(nulstat,'(2X,A9)',iostat=ios,err=10,end=10) OmPstdCH%stnids(stnidIndex)
2054
2055 ! Read (1) BUFR element
2056 ! (2) Flag indication if OmP error std dev provided from this auxiliary file or
2057 ! to be calculated via OmP differences or HPHT
2058 ! (3) Type of input/set std dev
2059 ! (4) Number of vertical levels
2060 ! (5) Number of latitudes
2061 ! (6) Number of months
2062 !
2063 ! Important: Combination of STNID and BUFR element
2064 ! to determine association to the observations.
2065
2066 read(nulstat,*,iostat=ios,err=10,end=10) &
2067 OmPstdCH%element(stnidIndex),OmPstdCH%source(stnidIndex), &
2068 OmPstdCH%std_type(stnidIndex),OmPstdCH%n_lvl(stnidIndex), &
2069 OmPstdCH%n_lat(stnidIndex),OmPstdCH%n_month(stnidIndex)
2070
2071 if ( OmPstdCH%n_lvl(stnidIndex) < 1 ) OmPstdCH%n_lvl(stnidIndex)=1
2072 if ( OmPstdCH%n_lat(stnidIndex) < 1 ) OmPstdCH%n_lat(stnidIndex)=1
2073 if ( OmPstdCH%n_month(stnidIndex) < 1 ) OmPstdCH%n_month(stnidIndex)=1
2074
2075 if ( icount+OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex)*OmPstdCH%n_month(stnidIndex) > isize ) then
2076 write(*,'(10X,"Max array size exceeded: ",I6)') isize
2077 CALL utl_abort('ose_readOmPstddev_auxfileHCHin: PROBLEM READING OBSERR STD DEV.')
2078 else if ( OmPstdCH%n_lat(stnidIndex) == 1 .and. OmPstdCH%n_month(stnidIndex) > 1 ) then
2079 write(*,'(10X,"Fails for stnid number: ",I6)') stnidIndex
2080 CALL utl_abort('ose_readOmPstddev_auxfileHCHin: Cannot depend on month if not dependent on latitude')
2081 else if ( OmPstdCH%n_month(stnidIndex) /= 1 .and. OmPstdCH%n_month(stnidIndex) /= 12 ) then
2082 write(*,'(10X,"Fails for stnid number: ",I6)') stnidIndex
2083 CALL utl_abort('ose_readOmPstddev_auxfileHCHin: Number of months must be 1 or 12')
2084 end if
2085
2086 ! disregard line of dashes
2087 read(nulstat,'(A)',iostat=ios,err=10,end=10) ligne
2088
2089 if ( OmPstdCH%source(stnidIndex) > 1 ) then
2090 ! Disregard data section
2091 OmPstdCH%ibegin(stnidIndex)=icount
2092 cycle STNIDLOOP
2093 else if ( OmPstdCH%source(stnidIndex) == 1 ) then
2094
2095 OmPstdCH%ibegin(stnidIndex)=icount+1
2096
2097 ! Only reference vertical levels and latitudes to be determined
2098 ! For use in determining error std dev from OmPs of current processing window
2099
2100 OmPstdCH%n_month(stnidIndex) = 1
2101 if ( OmPstdCH%n_lvl(stnidIndex) == 1 .and. OmPstdCH%n_lat(stnidIndex) == 1 ) then
2102
2103 icount=icount+1
2104
2105 else if ( OmPstdCH%n_lvl(stnidIndex) == 1 .and. OmPstdCH%n_lat(stnidIndex) > 1 ) then
2106
2107 ! Read reference latitudes (must be in order of increasing size)
2108
2109 read(nulstat,*,iostat=ios,err=10,end=10) &
2110 OmPstdCH%lat(icount+1:icount+OmPstdCH%n_lat(stnidIndex))
2111 icount=icount+OmPstdCH%n_lat(stnidIndex)
2112
2113 else if ( OmPstdCH%n_lvl(stnidIndex) > 1 .and. OmPstdCH%n_lat(stnidIndex) == 1 ) then
2114
2115 ! Read reference vertical levels (must be in order of increasing size)
2116
2117 read(nulstat,*,iostat=ios,err=10,end=10) &
2118 OmPstdCH%levels(icount+1:icount+OmPstdCH%n_lvl(stnidIndex))
2119 icount=icount+OmPstdCH%n_lvl(stnidIndex)
2120
2121 else
2122
2123 ! Read reference latitudes (must be in order of increasing size)
2124 read(nulstat,*,iostat=ios,err=10,end=10) &
2125 OmPstdCH%lat(icount+1:icount+OmPstdCH%n_lat(stnidIndex))
2126
2127 ! Read reference vertical levels (must be in order of increasing size)
2128 read(nulstat,*,iostat=ios,err=10,end=10) &
2129 OmPstdCH%levels(icount+1:icount+OmPstdCH%n_lvl(stnidIndex))
2130
2131 icount=icount+OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex)
2132
2133 end if
2134
2135 cycle STNIDLOOP
2136
2137 end if
2138
2139 ! For OmPstdCH%source(stnidIndex) == 0
2140
2141 OmPstdCH%ibegin(stnidIndex)=icount+1
2142 if ( OmPstdCH%n_lvl(stnidIndex) == 1 .and. OmPstdCH%n_lat(stnidIndex) == 1 ) then
2143
2144 ! Read one value only (independent of level and latitude)
2145 ! Assumes one month only.
2146
2147 read(nulstat,*,iostat=ios,err=10,end=10) OmPstdCH%std(icount+1)
2148 icount=icount+1
2149
2150 else if ( OmPstdCH%n_lvl(stnidIndex) == 1 .and. OmPstdCH%n_lat(stnidIndex) > 1 ) then
2151
2152 ! Value dependent on latitude
2153
2154 ! Read reference latitudes (must be in order of increasing size)
2155
2156 read(nulstat,*,iostat=ios,err=10,end=10) &
2157 OmPstdCH%lat(icount+1:icount+OmPstdCH%n_lat(stnidIndex))
2158
2159 ! Read OMP error std dev related values
2160
2161 if ( OmPstdCH%n_month(stnidIndex) == 1 ) then
2162
2163 read(nulstat,*,iostat=ios,err=10,end=10) &
2164 OmPstdCH%std(icount+1:icount+OmPstdCH%n_lat(stnidIndex))
2165
2166 else
2167
2168 do monthIndex=1,OmPstdCH%n_month(stnidIndex)
2169
2170 ! Read reference month (must be in order of increasing size)
2171 read(nulstat,*,iostat=ios,err=10,end=10) OmPstdCH%month(icount+monthIndex)
2172
2173 ipos=icount+(monthIndex-1)*OmPstdCH%n_lat(stnidIndex)
2174
2175 read(nulstat,*,iostat=ios,err=10,end=10) &
2176 OmPstdCH%std(ipos+1:ipos+OmPstdCH%n_lat(stnidIndex))
2177 end do
2178
2179 end if
2180
2181 icount = icount + OmPstdCH%n_lat(stnidIndex)*OmPstdCH%n_month(stnidIndex)
2182
2183 else if ( OmPstdCH%n_lvl(stnidIndex) > 1 .and. OmPstdCH%n_lat(stnidIndex) == 1 ) then
2184
2185 ! Value dependent on vertical level and not latitude
2186
2187 if ( OmPstdCH%n_month(stnidIndex) == 1 ) then
2188
2189 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)
2190 icount=icount+1
2191
2192 ! Read vertical level and OMP error std dev related value.
2193
2194 read(nulstat,*,iostat=ios,err=10,end=10) &
2195 OmPstdCH%levels(icount),OmPstdCH%std(icount)
2196
2197 end do
2198
2199 else
2200
2201 do monthIndex=1,OmPstdCH%n_month(stnidIndex)
2202
2203 ! Read reference month (must be in order of increasing size)
2204 read(nulstat,*,iostat=ios,err=10,end=10) OmPstdCH%month(icount+monthIndex)
2205
2206 ipos=icount+(monthIndex-1)*OmPstdCH%n_lvl(stnidIndex)
2207
2208 read(nulstat,*,iostat=ios,err=10,end=10) &
2209 OmPstdCH%std(ipos+1:ipos+OmPstdCH%n_lvl(stnidIndex))
2210 end do
2211
2212 end if
2213
2214 icount = icount + OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_month(stnidIndex)
2215
2216 else if ( OmPstdCH%n_lvl(stnidIndex) > 1 .and. OmPstdCH%n_lat(stnidIndex) > 1 ) then
2217
2218 ! Value dependent on vertical level and latitude
2219
2220 ! Read reference latitudes (must be in order of increasing size)
2221 read(nulstat,*,iostat=ios,err=10,end=10) &
2222 OmPstdCH%lat(icount+1:icount+OmPstdCH%n_lat(stnidIndex))
2223
2224 if ( OmPstdCH%n_month(stnidIndex) == 1 ) then
2225
2226 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)
2227
2228 ! Read vertical level and OMP error std dev related lat-dependent values.
2229
2230 read(nulstat,*,iostat=ios,err=10,end=10) &
2231 OmPstdCH%levels(icount+levIndex), &
2232 OmPstdCH%std(icount+(levIndex-1)* &
2233 OmPstdCH%n_lat(stnidIndex)+1:icount+levIndex*OmPstdCH%n_lat(stnidIndex))
2234
2235 end do
2236
2237 else
2238
2239 do monthIndex=1,OmPstdCH%n_month(stnidIndex)
2240
2241 ! Read reference month (must be in order of increasing size)
2242 read(nulstat,*,iostat=ios,err=10,end=10) OmPstdCH%month(icount+monthIndex)
2243
2244 ipos=icount+(monthIndex-1)*OmPstdCH%n_lat(stnidIndex)*OmPstdCH%n_lvl(stnidIndex)
2245
2246 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)
2247
2248 ! Read vertical level and OMP error std dev related lat-dependent values.
2249
2250 read(nulstat,*,iostat=ios,err=10,end=10) &
2251 OmPstdCH%levels(icount+levIndex), &
2252 OmPstdCH%std(ipos+(levIndex-1)* &
2253 OmPstdCH%n_lat(stnidIndex)+1:ipos+levIndex*OmPstdCH%n_lat(stnidIndex))
2254
2255 end do
2256
2257 end do
2258
2259 end if
2260
2261 icount=icount+OmPstdCH%n_lat(stnidIndex)*OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_month(stnidIndex)
2262 end if
2263 end do STNIDLOOP
2264
2265 10 if (ios.gt.0) then
2266 write(abortText,*) ios
2267 call utl_abort('ose_readOmPstddev_auxfileCHin: PROBLEM READING OMP ERR STD DEV. - ' // &
2268 'File read error message number: ' // trim(abortText) )
2269 end if
2270
2271 return
2272
2273 ! Reached end of file and no related section found.
2274
2275 15 OmPstdCH%n_stnid = 0
2276
2277 close(unit=nulstat)
2278 ierr=fclos(nulstat)
2279
2280 end subroutine ose_readOmPstddev_auxfileCH
2281
2282 !--------------------------------------------------------------------------
2283 ! ose_calcOmPstddevCH
2284 !--------------------------------------------------------------------------
2285 subroutine ose_calcOmPstddevCH(obsSpaceData)
2286 !
2287 !:Purpose: To calc OmP error std dev for some obs sets of the CH family
2288 !
2289 implicit none
2290
2291 ! Arguments:
2292 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data; output saved in OBS_OMPE column
2293
2294 ! Locals:
2295 logical :: availableOmP
2296 integer :: stnidIndex, headerIndex, bodyIndex, bodyIndex_start, bodyIndex_end, icodtyp, ierr
2297 integer :: idate, itime, iass, latIndex, levIndex, monthIndex, ibegin, loopIndex, posIndex
2298 real(8) :: zlat, zval, zlev, lat, sumOmP, sumSqrOmP, varOmP, maxOmP,meanOmP,medianOmP
2299 character(len=12) :: stnid
2300 real(8), allocatable :: series(:,:,:),sumOmP2d(:,:),sumSqrOmP2d(:,:)
2301 real(8), allocatable :: sumOmP2dt(:,:),sumSqrOmP2dt(:,:)
2302 integer, allocatable :: nSeries(:,:),nSeriest(:,:)
2303 integer, parameter :: maxCount = 10000
2304 integer, parameter :: minCount = 5
2305 real(8), parameter :: stdScale = 2.0
2306 real(8), parameter :: minVal = 1.0d-20
2307 real(4) :: rseries(maxCount)
2308 integer :: ip(maxCount)
2309
2310 ! Loop over all obs types
2311
2312 do stnidIndex=1,OmPstdCH%n_stnid
2313 ! Check if OBS_OMPE is to be estimated elsewhere via use of OBS_HPHT and OBS_OER
2314 ! instead of via ose_fillOmPstddevCH
2315 if ( OmPstdCH%source(stnidIndex) /= 1 ) cycle
2316
2317 write(*,*)
2318 write(*,*) 'ose_calcOmPstddevCH: Online calculation of OmP error std dev of CH family stnid ',OmPstdCH%stnids(stnidIndex)
2319 write(*,*)
2320
2321 ! Identify start position index (-1) in OmPstdCH
2322 ibegin=OmPstdCH%ibegin(stnidIndex)-1
2323
2324 allocate(series(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex),maxCount))
2325 allocate(nSeries(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2326 allocate(sumOmP2d(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2327 allocate(sumSqrOmP2d(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2328 nSeries(:,:)=0
2329 sumOmP2d(:,:)=0.0d0
2330 sumSqrOmP2d(:,:)=0.0d0
2331
2332 allocate(nSeriest(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2333 allocate(sumOmP2dt(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2334 allocate(sumSqrOmP2dt(OmPstdCH%n_lvl(stnidIndex),OmPstdCH%n_lat(stnidIndex)))
2335 nSeriest(:,:)=0
2336
2337 ! Loop over all header indices of the 'CH' family:
2338
2339 call obs_set_current_header_list(obsSpaceData,'CH')
2340 HEADER: do
2341
2342 headerIndex = obs_getHeaderIndex(obsSpaceData)
2343 if (headerIndex < 0) exit HEADER
2344
2345 icodtyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
2346 if (icodtyp.ne.codtyp_get_codtyp('CHEMREMOTE').and.icodtyp.ne.codtyp_get_codtyp('CHEMINSITU')) cycle HEADER
2347
2348 stnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
2349 if ( .not. utl_stnid_equal(OmPstdCH%stnids(stnidIndex),stnid) ) cycle HEADER
2350
2351 bodyIndex_start = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2352 bodyIndex_end = bodyIndex_start+obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
2353
2354 ! Check OBS_VNM value.
2355 do bodyIndex=bodyIndex_start,bodyIndex_end
2356 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= BUFR_SCALE_EXPONENT) then
2357 if ( obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= OmPstdCH%element(stnidIndex) ) cycle HEADER
2358 end if
2359 end do
2360
2361 zlat = obs_headElem_r( obsSpaceData, OBS_LAT, headerIndex )
2362 idate = obs_headElem_i( obsSpaceData, OBS_DAT, headerIndex )
2363 itime = obs_headElem_i( obsSpaceData, OBS_ETM, headerIndex )
2364
2365 ! Identify lat and month index
2366 if (OmPstdCH%n_lat(stnidIndex) == 1) then
2367 latIndex = 1
2368 else
2369
2370 ! Find latitude index for specifying the latitude bin
2371 ! Assuming increasing latitudes in OmPstdCH%lat
2372
2373 lat = zlat / MPC_RADIANS_PER_DEGREE_R8 ! radians to degrees
2374
2375 if (lat >= 0.5*(OmPstdCH%lat(ibegin+OmPstdCH%n_lat(stnidIndex)) &
2376 +OmPstdCH%lat(ibegin+OmPstdCH%n_lat(stnidIndex)-1)) ) then
2377 latIndex = OmPstdCH%n_lat(stnidIndex)
2378 else
2379 do latIndex=1,OmPstdCH%n_lat(stnidIndex)-1
2380 if (lat <= 0.5*(OmPstdCH%lat(ibegin+latIndex)+OmPstdCH%lat(ibegin+latIndex+1)) ) exit
2381 end do
2382 end if
2383 end if
2384
2385 if (OmPstdCH%n_month(stnidIndex) == 1) then
2386 monthIndex = 1
2387 else
2388 ! Find month index
2389 monthIndex=(idate-(idate/10000)*10000)/100
2390 end if
2391
2392 do bodyIndex=bodyIndex_start,bodyIndex_end
2393 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == BUFR_SCALE_EXPONENT) cycle
2394
2395 iass = obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex )
2396 if ( iass /= obs_assimilated ) cycle
2397
2398 ! Identify vertical level
2399 zlev = obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex )
2400 if ( OmPstdCH%n_lvl(stnidIndex) == 1 ) then
2401 levIndex = 1
2402 else if ( OmPstdCH%levels(ibegin+1) < OmPstdCH%levels(ibegin+2) ) then
2403 if (zlev <= 0.5*(OmPstdCH%levels(ibegin+1)+OmPstdCH%levels(ibegin+2)) ) then
2404 levIndex = 1
2405 else
2406 do levIndex=OmPstdCH%n_lvl(stnidIndex),2,-1
2407 if ( zlev >= 0.5*(OmPstdCH%levels(ibegin+levIndex)+OmPstdCH%levels(ibegin+levIndex-1)) ) exit
2408 end do
2409 end if
2410 else
2411 if (zlev <= 0.5*(OmPstdCH%levels(ibegin+OmPstdCH%n_lvl(stnidIndex)) &
2412 +OmPstdCH%levels(ibegin+OmPstdCH%n_lvl(stnidIndex)-1)) ) then
2413 levIndex = OmPstdCH%n_lvl(stnidIndex)
2414 else
2415 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)-1
2416 if ( zlev >= 0.5*(OmPstdCH%levels(ibegin+levIndex)+OmPstdCH%levels(ibegin+levIndex+1)) ) exit
2417 end do
2418 end if
2419 end if
2420
2421 ! Store OmP in work array
2422
2423 if ( OmPstdCH%std_type(stnidIndex) == 1 ) then
2424 zval = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
2425 if ( zval > minVal ) then
2426 nSeries(levIndex,latIndex) = nSeries(levIndex,latIndex) + 1
2427 series(levIndex,latIndex,nSeries(levIndex,latIndex)) = obs_bodyElem_r( obsSpaceData, OBS_OMP, bodyIndex )/zval
2428 end if
2429 else
2430 nSeries(levIndex,latIndex) = nSeries(levIndex,latIndex) + 1
2431 series(levIndex,latIndex,nSeries(levIndex,latIndex)) = obs_bodyElem_r( obsSpaceData, OBS_OMP, bodyIndex )
2432 end if
2433 ! Check if enough data
2434 if ( nSeries(levIndex,latIndex) == maxCount ) exit
2435
2436 end do
2437
2438 end do HEADER
2439
2440 if ( any(nSeries > 0) ) then
2441
2442 ! Calc OmP error std dev
2443
2444 !write(*,*) 'latbin levbin Npts AvgOmP OmPStddev'
2445 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)
2446 do latIndex=1,OmPstdCH%n_lat(stnidIndex)
2447 availableOmP=.true.
2448 if (nSeries(levIndex,latIndex) > minCount ) then
2449 rseries(1:nSeries(levIndex,latIndex))=series(levIndex,latIndex,1:nSeries(levIndex,latIndex))
2450 call ipsort(ip,rseries(1:nSeries(levIndex,latIndex)),nSeries(levIndex,latIndex))
2451 medianOmP=series(levIndex,latIndex,ip(nSeries(levIndex,latIndex)/2))
2452
2453 sumOmP = sum(series(levIndex,latIndex,1:nSeries(levIndex,latIndex)))
2454 meanOmP=sumOmP/nSeries(levIndex,latIndex)
2455 sumSqrOmP = sum(series(levIndex,latIndex,1:nSeries(levIndex,latIndex))*series(levIndex,latIndex,1:nSeries(levIndex,latIndex)))
2456 varOmP = sumSqrOmP/nSeries(levIndex,latIndex) - (sumOmP/nSeries(levIndex,latIndex))**2
2457 if (varOmP > minVal*minVal ) then
2458 do loopIndex=1,nSeries(levIndex,latIndex)
2459 maxOmP = stdScale*sqrt(varOmP)
2460 ! if ( abs(series(levIndex,latIndex,loopIndex)-meanOmP) > maxOmP ) then
2461 if ( abs(series(levIndex,latIndex,loopIndex)-medianOmP) > maxOmP ) then
2462 sumOmP = sumOmP - series(levIndex,latIndex,loopIndex)
2463 sumSqrOmP = sumSqrOmP - series(levIndex,latIndex,loopIndex)*series(levIndex,latIndex,loopIndex)
2464 nSeries(levIndex,latIndex) = nSeries(levIndex,latIndex) - 1
2465 end if
2466 end do
2467 varOmP = sumSqrOmP/nSeries(levIndex,latIndex) - (sumOmP/nSeries(levIndex,latIndex))**2
2468 !write(*,110) latIndex,levIndex,nSeries(levIndex,latIndex),sumOmP/nSeries(levIndex,latIndex),sqrt(varOmP)
2469 sumOmP2d(levIndex,latIndex)=sumOmP
2470 sumSqrOmP2d(levIndex,latIndex)=sumSqrOmP
2471 if (varOmP < minVal*minVal ) availableOmP = .false.
2472 else
2473 availableOmP = .false.
2474 end if
2475 else
2476 availableOmP = .false.
2477 end if
2478 if (.not.availableOmP) nSeries(levIndex,latIndex) = 0
2479 end do
2480 end do
2481 end if
2482
2483 ! Combine from all processors
2484
2485 call rpn_comm_allreduce(nSeries,nSeriest,OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex),"MPI_INTEGER","MPI_SUM","GRID",ierr)
2486 call rpn_comm_allreduce(sumOmP2d,sumOmP2dt,OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex),"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
2487 call rpn_comm_allreduce(sumSqrOmP2d,sumSqrOmP2dt,OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex),"MPI_DOUBLE_PRECISION","MPI_SUM","GRID",ierr)
2488
2489 if (any(nSeriest > 5*minCount)) then
2490
2491 where ( nSeriest > 5*minCount ) sumSqrOmP2dt = sumSqrOmP2dt/nSeriest - (sumOmP2dt/nSeriest)**2
2492
2493 if (mmpi_myid == 0) then
2494 write(*,*) 'Calculated OmP error std dev'
2495 write(*,*) 'latbin levbin Npts AvgOmP OmPStddev'
2496 end if
2497 do levIndex=1,OmPstdCH%n_lvl(stnidIndex)
2498 do latIndex=1,OmPstdCH%n_lat(stnidIndex)
2499 availableOmP=.true.
2500 if (nSeriest(levIndex,latIndex) > 5*minCount ) then
2501 varOmP = sumSqrOmP2dt(levIndex,latIndex)
2502 if (varOmP < minVal*minVal ) then
2503 availableOmP = .false.
2504 else
2505 if (mmpi_myid == 0 ) write(*,110) latIndex,levIndex,nSeriest(levIndex,latIndex),sumOmP2dt(levIndex,latIndex)/nSeriest(levIndex,latIndex),sqrt(varOmP)
2506 end if
2507 else
2508 availableOmP = .false.
2509 end if
2510
2511 ! Store in structure
2512
2513 ! Indentify position in structure
2514 if (OmPstdCH%n_lvl(stnidIndex) > 1) then
2515 posIndex=ibegin+(levIndex-1)*OmPstdCH%n_lat(stnidIndex)
2516 else
2517 posIndex=ibegin
2518 end if
2519
2520 if ( OmPstdCH%n_month(stnidIndex) == 1 ) then
2521 posIndex = posIndex + latIndex
2522 else
2523 posIndex = posIndex + (monthIndex-1)*OmPstdCH%n_lvl(stnidIndex)*OmPstdCH%n_lat(stnidIndex) + latIndex
2524 end if
2525
2526 if (availableOmP) then
2527 OmPstdCH%std(posIndex) = sqrt(varOmP)
2528 else
2529 OmPstdCH%std(posIndex) = MPC_missingValue_R8
2530 end if
2531
2532 end do
2533 end do
2534 end if
2535110 format(I5,I7,I7,3x,G11.3,G11.3)
2536
2537 deallocate(series,nSeries,sumOmP2d,sumSqrOmP2d)
2538 deallocate(nSeriest,sumOmP2dt,sumSqrOmP2dt)
2539 end do
2540
2541 end subroutine ose_calcOmPstddevCH
2542
2543 !--------------------------------------------------------------------------
2544 ! ose_fillOmPstddevCH
2545 !--------------------------------------------------------------------------
2546 subroutine ose_fillOmPstddevCH(obsSpaceData)
2547 !
2548 !:Purpose: To assign the Omp error std dev where possible for the obs
2549 ! of the CH obs family.
2550 !
2551 implicit none
2552
2553 ! Arguments:
2554 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data; output saved in OBS_OMPE column
2555
2556 ! Locals:
2557 integer :: stnidIndex, headerIndex, bodyIndex, bodyIndex_start, bodyIndex_end, icodtyp
2558 integer :: idate, itime, iass, latIndex, monthIndex, ibegin
2559 real(8) :: zlat, zlev, OmP_err_stddev, lat
2560 character(len=12) :: stnid
2561
2562 ! Loop over all obs types
2563
2564 do stnidIndex=1,OmPstdCH%n_stnid
2565 ! Check if OBS_OMPE is to be estimated elsewhere via use of OBS_HPHT and OBS_OER
2566 ! instead of via ose_fillOmPstddevCH
2567 if ( OmPstdCH%source(stnidIndex) > 1 ) cycle
2568
2569 ! Loop over all header indices of the 'CH' family:
2570
2571 call obs_set_current_header_list(obsSpaceData,'CH')
2572 HEADER: do
2573
2574 headerIndex = obs_getHeaderIndex(obsSpaceData)
2575 if (headerIndex < 0) exit HEADER
2576
2577 icodtyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
2578 if (icodtyp.ne.codtyp_get_codtyp('CHEMREMOTE').and.icodtyp.ne.codtyp_get_codtyp('CHEMINSITU')) cycle HEADER
2579
2580 stnid = obs_elem_c(obsSpaceData,'STID',headerIndex)
2581 if ( .not. utl_stnid_equal(OmPstdCH%stnids(stnidIndex),stnid) ) cycle HEADER
2582
2583 bodyIndex_start = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2584 bodyIndex_end = bodyIndex_start+obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
2585
2586 ! Check OBS_VNM value.
2587 do bodyIndex=bodyIndex_start,bodyIndex_end
2588 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= BUFR_SCALE_EXPONENT) then
2589 if ( obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) /= OmPstdCH%element(stnidIndex) ) cycle HEADER
2590 end if
2591 end do
2592
2593 zlat = obs_headElem_r( obsSpaceData, OBS_LAT, headerIndex )
2594 idate = obs_headElem_i( obsSpaceData, OBS_DAT, headerIndex )
2595 itime = obs_headElem_i( obsSpaceData, OBS_ETM, headerIndex )
2596
2597 ! Identify lat and month index
2598 if (OmPstdCH%n_lat(stnidIndex) == 1) then
2599 latIndex = 0 ! no interpolation
2600 else
2601
2602 ! Find latitude index for interpolation.
2603 ! Assuming increasing latitudes in OmPstdCH%lat
2604 ! Interpolation to between latIndex and latIndex-1
2605
2606 lat = zlat / MPC_RADIANS_PER_DEGREE_R8 ! radians to degrees
2607
2608 ibegin=OmPstdCH%ibegin(stnidIndex)-1
2609 if ( lat <= OmPstdCH%lat(ibegin+1) ) then
2610 latIndex=2
2611 else if ( lat >= OmPstdCH%lat(ibegin+OmPstdCH%n_lat(stnidIndex)) ) then
2612 latIndex=OmPstdCH%n_lat(stnidIndex)
2613 else
2614 do latIndex=2,OmPstdCH%n_lat(stnidIndex)
2615 if (lat >= OmPstdCH%lat(ibegin+latIndex-1) .and. &
2616 lat <= OmPstdCH%lat(ibegin+latIndex)) exit
2617 end do
2618 end if
2619 end if
2620
2621 if (OmPstdCH%n_month(stnidIndex) == 1) then
2622 monthIndex = 1
2623 else
2624 ! Find month index
2625 monthIndex=(idate-(idate/10000)*10000)/100
2626 end if
2627
2628 do bodyIndex=bodyIndex_start,bodyIndex_end
2629 if (obs_bodyElem_i(obsSpaceData,OBS_VNM,bodyIndex) == BUFR_SCALE_EXPONENT) cycle
2630
2631 zlev = obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex )
2632
2633 ! Get OmP error std dev
2634 OmP_err_stddev = ose_getOmPstddevCH( lat, zlev, stnidIndex, latIndex, monthIndex )
2635 if ( OmPstdCH%std_type(stnidIndex) == 1 ) then
2636 iass = obs_bodyElem_i( obsSpaceData, OBS_ASS, bodyIndex )
2637 if ( iass == obs_assimilated ) then
2638 OmP_err_stddev = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )*OmP_err_stddev
2639 call obs_bodySet_r( obsSpaceData, OBS_OMPE, bodyIndex, OmP_err_stddev )
2640 end if
2641 else
2642 call obs_bodySet_r( obsSpaceData, OBS_OMPE, bodyIndex, OmP_err_stddev )
2643 end if
2644
2645 end do
2646 end do HEADER
2647 end do
2648
2649 end subroutine ose_fillOmPstddevCH
2650
2651 !--------------------------------------------------------------------------
2652 ! ose_getOmPstddevCH
2653 !--------------------------------------------------------------------------
2654 real(8) function ose_getOmPstddevCH(zlat,zlev,stnidIndex,latIndex,monthIndex) result(OmP_err_stddev)
2655 !
2656 !:Purpose: To return the OmP error std dev for a CH family measurement
2657 !
2658 implicit none
2659
2660 ! Arguments:
2661 real(8), intent(in) :: zlat ! latitude (radians)
2662 real(8), intent(in) :: zlev ! vertical coordinate value
2663 integer, intent(in) :: stnidIndex ! station and obs type index
2664 integer, intent(in) :: latIndex ! reference lat for interpolation
2665 integer, intent(in) :: monthIndex ! month index
2666
2667 ! Locals:
2668 real(8) :: levDiff
2669 integer :: ibegin,levIndex,loopIndex,posIndex
2670 real(8), parameter :: minDiff = 0.01*abs(MPC_missingValue_R8)
2671
2672 OmP_err_stddev = MPC_missingValue_R8
2673
2674 ! Get OmP error std. dev.
2675
2676 ibegin=OmPstdCH%ibegin(stnidIndex)-1
2677
2678 if (OmPstdCH%n_lvl(stnidIndex) > 1) then
2679
2680 ! Find nearest vertical level (no vertical interpolation in this version)
2681
2682 levDiff=1.0d10
2683 do loopIndex=1,OmPstdCH%n_lvl(stnidIndex)
2684 if ( levDiff > abs(zlev-OmPstdCH%levels(ibegin+loopIndex)) ) THEN
2685 levIndex=loopIndex
2686 levDiff=abs(zlev-OmPstdCH%levels(ibegin+loopIndex))
2687 END IF
2688 END DO
2689 posIndex=ibegin+(levIndex-1)*OmPstdCH%n_lat(stnidIndex)
2690 else
2691 posIndex=ibegin
2692 end if
2693
2694 if ( OmPstdCH%n_month(stnidIndex) == 1 ) then
2695 posIndex = posIndex + latIndex
2696 else
2697 posIndex = posIndex + (monthIndex-1)*OmPstdCH%n_lat(stnidIndex)*OmPstdCH%n_lvl(stnidIndex) + latIndex
2698 end if
2699
2700 if ( OmPstdCH%n_lat(stnidIndex) > 1 ) then
2701 ! Apply interpolation in latitude
2702
2703 if ( latIndex == 1 .or. latIndex > OmPstdCH%n_lat(stnidIndex)) then
2704 if ( abs(OmPstdCH%std(posIndex) - MPC_missingValue_R8) > minDiff ) then
2705 OmP_err_stddev = MPC_missingValue_R8
2706 else
2707 OmP_err_stddev = OmPstdCH%std(posIndex)
2708 end if
2709 else
2710 if ( abs(OmPstdCH%std(posIndex-1) - MPC_missingValue_R8) < minDiff .and. &
2711 abs(OmPstdCH%std(posIndex) - MPC_missingValue_R8) < minDiff ) then
2712 OmP_err_stddev = MPC_missingValue_R8
2713 else if ( abs(OmPstdCH%std(posIndex-1) - MPC_missingValue_R8) < minDiff ) then
2714 OmP_err_stddev = OmPstdCH%std(posIndex)
2715 else if ( abs(OmPstdCH%std(posIndex) - MPC_missingValue_R8) < minDiff ) then
2716 OmP_err_stddev = OmPstdCH%std(posIndex-1)
2717 else
2718 OmP_err_stddev = (OmPstdCH%std(posIndex-1)*(OmPstdCH%lat(ibegin+latIndex)-zlat)+ &
2719 OmPstdCH%std(posIndex)*(zlat-OmPstdCH%lat(ibegin+latIndex-1)))/ &
2720 (OmPstdCH%lat(ibegin+latIndex)-OmPstdCH%lat(ibegin+latIndex-1))
2721 end if
2722 end if
2723 else
2724 if ( abs(OmPstdCH%std(posIndex) - MPC_missingValue_R8) > minDiff ) then
2725 OmP_err_stddev = MPC_missingValue_R8
2726 else
2727 OmP_err_stddev = OmPstdCH%std(posIndex)
2728 end if
2729 end if
2730
2731 end function ose_getOmPstddevCH
2732
2733 !--------------------------------------------------------------------------
2734 ! ose_OmPstddevExistsForAllCH
2735 !--------------------------------------------------------------------------
2736 logical function ose_OmPstddevExistsForAllCH(obsSpaceData) result(allOMPE)
2737 !
2738 !:Purpose: To determine if all obs to be processed have usable OBS_OMPE values
2739 ! for the CH obs family.
2740 !
2741 implicit none
2742
2743 ! Arguments:
2744 type(struct_obs), intent(inout) :: obsSpaceData ! observation-space data; output saved in OBS_OMPE column
2745
2746 ! Local:
2747 integer :: headerIndex, bodyIndex, bodyIndex_start, bodyIndex_end, icodtyp
2748
2749 ! Loop over all header indices of the 'CH' family:
2750
2751 allOMPE = .true.
2752
2753 call obs_set_current_header_list(obsSpaceData,'CH')
2754 HEADER: do
2755
2756 headerIndex = obs_getHeaderIndex(obsSpaceData)
2757 if (headerIndex < 0) exit HEADER
2758
2759 icodtyp = obs_headElem_i(obsSpaceData,OBS_ITY,headerIndex)
2760 if (icodtyp.ne.codtyp_get_codtyp('CHEMREMOTE').and.icodtyp.ne.codtyp_get_codtyp('CHEMINSITU')) cycle HEADER
2761
2762 bodyIndex_start = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
2763 bodyIndex_end = bodyIndex_start+obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex)-1
2764
2765 ! Check for cases were OmP error std dev is not available
2766
2767 do bodyIndex=bodyIndex_start,bodyIndex_end
2768 if ( obs_bodyElem_r(obsSpaceData,OBS_OMPE,bodyIndex) <= 0.0d0 .and. &
2769 obs_bodyElem_i(obsSpaceData,OBS_ASS,bodyIndex) == obs_assimilated ) then
2770 allOMPE = .false.
2771 write(*,*)
2772 write(*,*) 'ose_OmPstddevExistsForAllCH: Not all CH obs to be processed have available OmP error std dev.'
2773 write(*,*) 'HPHT will be calculated for the obs that do not have OmP error std dev available.'
2774 write(*,*)
2775 return
2776 end if
2777 end do
2778 end do HEADER
2779
2780 end function ose_OmPstddevExistsForAllCH
2781
2782 !--------------------------------------------------------------------------
2783 ! ose_deallocOmPstddevCH
2784 !--------------------------------------------------------------------------
2785 subroutine ose_deallocOmPstddevCH
2786 !
2787 !:Purpose: To deallocate temporary storage space used for OmP error std dev
2788 ! for the CH family.
2789 !
2790 implicit none
2791
2792 if (OmPstdCH%n_stnid.eq.0) return
2793
2794 if (allocated(OmPstdCH%stnids)) deallocate(OmPstdCH%stnids)
2795 if (allocated(OmPstdCH%n_lvl)) deallocate(OmPstdCH%n_lvl)
2796 if (allocated(OmPstdCH%ibegin)) deallocate(OmPstdCH%ibegin)
2797 if (allocated(OmPstdCH%element)) deallocate(OmPstdCH%element)
2798 if (allocated(OmPstdCH%source)) deallocate(OmPstdCH%source)
2799 if (allocated(OmPstdCH%std_type)) deallocate(OmPstdCH%std_type)
2800 if (allocated(OmPstdCH%n_lat)) deallocate(OmPstdCH%n_lat)
2801 if (allocated(OmPstdCH%std)) deallocate(OmPstdCH%std)
2802 if (allocated(OmPstdCH%levels)) deallocate(OmPstdCH%levels)
2803 if (allocated(OmPstdCH%lat)) deallocate(OmPstdCH%lat)
2804 if (allocated(OmPstdCH%n_month)) deallocate(OmPstdCH%n_month)
2805 if (allocated(OmPstdCH%month)) deallocate(OmPstdCH%month)
2806
2807 end subroutine ose_deallocOmPstddevCH
2808
2809end module obsSpaceErrorStdDev_mod
