1module backgroundCheck_mod
2 ! MODULE backgroundCheck_mod (prefix='bgck' category='1. High-level functionality')
3 !
4 !:Purpose: Performs the background check on all conventional observations
5 !
6 use mathPhysConstants_mod
7 use bufr_mod
8 use obsSpaceData_mod
9 use gps_mod
10 use utilities_mod
11 use columnData_mod
12 use obsSpaceDiag_mod
13 use verticalCoord_mod
14 use horizontalCoord_mod
15 use obsSpaceErrorStdDev_mod
16 use obsFamilyList_mod
17 use codtyp_mod
18
19 implicit none
20 private
21
22 ! public procedures
23 public :: bgck_bgcheck_conv
24
25 integer, parameter :: numFamilyToProcess = 11
26 character(len=2), parameter :: familyListToProcess(numFamilyToProcess)= (/ &
27 'UA','AI','SF','SC','SW','PR','RO','GP','RA', &
28 'CH','AL' /)
29
30 contains
31
32 !--------------------------------------------------------------------------
33 ! bgck_bgcheck_conv
34 !--------------------------------------------------------------------------
35 subroutine bgck_bgcheck_conv( columnTrlOnAnlIncLev, columnTrlOnTrlLev, hco_anl, obsSpaceData )
36 !
37 !:Purpose: Do background check on all conventional observations
38 !
39 implicit none
40
41 ! Arguments:
42 type(struct_obs) , intent(inout) :: obsSpaceData ! Observation-related data
43 type(struct_columnData) , intent(inout) :: columnTrlOnAnlIncLev ! column data on analysis levels
44 type(struct_columnData) , intent(inout) :: columnTrlOnTrlLev ! column data on trial levels
45 type(struct_hco), pointer, intent(in) :: hco_anl ! horizontal grid structure
46
47 ! Locals:
48 integer :: familyIndex
49 integer :: nulNam, ier, fnom, fclos
50 logical :: noObsToProcess
51 character(len=*), parameter :: myName = 'bgck_bgcheck_conv'
52
53 ! namelist variables
54 logical :: new_bgck_sw ! choose to use the 'new' background check for SW obs
55 namelist /NAMBGCKCONV/ new_bgck_sw
56
57 write(*,*) myName//' begin conventional data background check'
58
59 ! Check if any observations are present for conventional background check
60 noObsToProcess = .true.
61 do familyIndex = 1, numFamilyToProcess
62 if (obs_famExist(obsSpaceData,familyListToProcess(familyIndex))) then
63 noObsToProcess = .false.
64 end if
65 end do
66 if (noObsToProcess) then
67 write(*,*) myName//': No observations to process'
68 return
69 end if
70
71 call utl_tmg_start(117,'--BgckConventional')
72
73 new_bgck_sw = .false.
74
75 nulNam=0
76 ier = fnom( nulNam, 'flnml', 'r/o', 0 )
77 read( nulNam, nml = NAMBGCKCONV, IOSTAT = ier )
78 if ( ier /= 0 ) write(*,*) myName//': no valid namelist NAMBGCKCONV found, default values will be taken:'
79 write(*,*) myName//': new_bgck_sw = ',new_bgck_sw
80 ier = fclos(nulNam)
81
82
83 ! Obtain or calc OmP-error std dev when requested and possible.
84 ! Otherwise calc HBHT contribution (sigma_B in observation space)
85 ! --------------------------------------------------------------------
86 call ose_computeStddev( columnTrlOnAnlIncLev, hco_anl, & ! IN
87 obsSpaceData ) ! INOUT
88
89 ! DO A BACKGROUND CHECK ON ALL THE OBSERVATIONS
90 ! ----------------------------------------------
91
92 do familyIndex = 1, ofl_numFamily
93 ! For SW only, old and new background check schemes controlled by "new_bgck_sw"
94 if ( obs_famExist( obsSpaceData, ofl_familyList( familyIndex )) ) then
95 call bgck_data( ofl_familyList(familyIndex), obsSpaceData, new_bgck_sw )
96 end if
97 end do
98
99 if (obs_famExist( obsSpaceData, 'RO' )) call bgck_gpsro( columnTrlOnTrlLev , obsSpaceData )
100
101 ! Conduct obs-space post-processing diagnostic tasks (some diagnostic
102 ! computations controlled by NAMOSD namelist in flnml)
103
104 call osd_ObsSpaceDiag( obsSpaceData, columnTrlOnAnlIncLev, hco_anl, analysisMode_opt = .false. )
105
106 call utl_tmg_stop(117)
107
108 end subroutine bgck_bgcheck_conv
109
110 !--------------------------------------------------------------------------
111 ! bgck_data
112 !--------------------------------------------------------------------------
113 subroutine bgck_data( obsFamily, obsData, new_bgck_sw )
114 !
115 !:Purpose: Calculate a background check for a data family and set the
116 ! appropriate quality-control flags in obsSpaceData
117 implicit none
118
119 ! NOTE 1: gps_gb_YSFERRWGT IN MODGPSZTD_MOD (FROM NML FILE) IS USED HERE FOR ERROR WEIGHTING
120 ! OF TIME SERIES (FGAT) GPS MET OBSERVATIONS PS, TS, DPDS. IT IS APPLIED
121 ! (ALONG WITH gps_gb_YZDERRWGT FOR ZTD) IN S/R SETERR AS A MULT. FACTOR TO ERRORS.
122 ! gps_gb_YZDERRWGT AND gps_gb_YSFERRWGT = 1 FOR NORMAL 3D-VAR (3D-THINNING).
123
124 ! Arguments:
125 character(len=2), intent(in) :: obsFamily ! current observation family
126 type(struct_obs), intent(inout) :: obsData ! obsSpaceData
127 logical , intent(in) :: new_bgck_sw
128
129 ! Locals:
130 real(8) :: averageFGE, averageOER
131 integer :: obsFlag, obsVarno, headerIndex, codeType, bodyIndex, obsCount
132 integer :: numberObs, numberObsRejected, INZOBS, INZREJ
133 integer :: INPOBS, INTOBS, INDOBS, INPREJ, INTREJ, INDREJ
134 real(8) :: ZOER,ZOMP,ZFGE,ZOMPER,ZBGCHK,ZVAR,ZLEV,ZLAT,ZLON,ZSOP
135 logical :: LLOK, LLZD
136 character(len=12) :: stnid
137 integer :: i_ass, i_vco, i_oth, bodyIndex_u, bodyIndex_v, bodyIndex_start
138 real(8) :: uu_d, uu_r, uu_f, vv_d, vv_r, vv_f, duv2, duv2_lim, zslev
139 logical :: found_u, found_v
140 character(len=*), parameter :: myName = 'bgck_data'
141
142 write(*,*)' '
143 write(*,*) ' ------------------------------'
144 write(*,*) myName//' background check for', obsFamily, ' data'
145 write(*,*) ' ------------------------------'
146 write(*,*) ' '
147 write(*,'(a55,a74)')' STNID LATITU LONGITU ID Elem Level ', &
148 ' Value Sigmao Sigmap O-P SigmaOP qcflag '
149 write(*,'(a55,a74)')' ----- ------ ------- -- ---- ----- ', &
150 ' ----- ------ ------ --- ------- ------ '
151
152 obsCount = 0
153 averageFGE = 0.d0
154 averageOER=0.d0
155
156 numberObs = 0
157 numberObsRejected = 0
158
159 ! Initialize counters for GP family observations
160
161 if ( obsFamily == 'GP' ) then
162 INZOBS=0
163 INPOBS=0
164 INTOBS=0
165 INDOBS=0
166 INZREJ=0
167 INPREJ=0
168 INTREJ=0
169 INDREJ=0
170 end if
171
172 if ( .not. new_bgck_sw .or. obsFamily /= 'SW' ) then
173
174 ! loop over all header indices of the current observation family
175 call obs_set_current_header_list( obsData, obsFamily )
176 HEADER: do
177 headerIndex = obs_getHeaderIndex( obsData )
178 if ( headerIndex < 0 ) exit HEADER
179
180 stnid = obs_elem_c( obsData, 'STID', headerIndex )
181
182 ! 1. Computation of (HX - Z)**2/(SIGMAo**2 +SIGMAp**2)
183 ! ----------------------------------------------------
184
185 ! loop over all body indices for the current headerIndex
186 call obs_set_current_body_list( obsData, headerIndex )
187 BODY: do
188 bodyIndex = obs_getBodyIndex( obsData )
189 if (bodyIndex < 0 ) exit BODY
190
191 obsVarno = obs_bodyElem_i( obsData, OBS_VNM, bodyIndex )
192 LLOK = ( obs_bodyElem_i( obsData, OBS_ASS, bodyIndex ) == obs_assimilated )
193 if ( LLOK ) then
194 numberObs = numberObs + 1
195 if ( obsFamily == 'GP' ) then
196 if ( obsVarno == BUFR_NEZD ) INZOBS = INZOBS+1
197 if ( obsVarno == BUFR_NEPS ) INPOBS = INPOBS+1
198 if ( obsVarno == BUFR_NETS ) INTOBS = INTOBS+1
199 if ( obsVarno == BUFR_NESS ) INDOBS = INDOBS+1
200 end if
201 ZVAR = obs_bodyElem_r( obsData, OBS_VAR, bodyIndex )
202 ZLEV = obs_bodyElem_r( obsData, OBS_PPP, bodyIndex )
203 ZLAT = obs_headElem_r( obsData, OBS_LAT, headerIndex ) * MPC_DEGREES_PER_RADIAN_R8
204 ZLON = obs_headElem_r( obsData, OBS_LON, headerIndex ) * MPC_DEGREES_PER_RADIAN_R8
205
206 ! BACKGROUND CHECK
207
208 ZOMP = obs_bodyElem_r( obsData, OBS_OMP, bodyIndex )
209
210 ! Get error std dev
211
212 ! std(O-P) (valid/available if > 0.0)
213 zomper = obs_bodyElem_r( obsData, OBS_OMPE, bodyIndex )
214
215 ! std(O)
216 zoer = obs_bodyElem_r( obsData, OBS_OER, bodyIndex )
217 ! std(P)
218 zfge = obs_bodyElem_r( obsData, OBS_HPHT, bodyIndex )
219 ! NOTE: For GB-GPS ZTD observations (GP family), ZFGE is not the FGE but
220 ! rather Std(O-P) set in s/r SETERRGPSGB.
221 if ( obsFamily == 'GP' .and. obsVarno == BUFR_NEZD ) then
222 if ( zomper <= 0.0d0 .and. zfge > 0.0d0 ) zomper = zfge
223 end if
224
225 ! Protect against error std dev values that are too small
226 if ( obsFamily == 'GP' ) then
227 if ( obsVarno == BUFR_NEZD ) then
228 ZOER = ZOER / gps_gb_YZDERRWGT
229 else
230 ZOER = ZOER / gps_gb_YSFERRWGT
231 end if
232 if ( ZOER < 1.0d-3 .and. obsVarno /= BUFR_NEZD ) then
233 write(*,*)' Problem for GP STNID ZOER= ' , stnid, ZOER
234 call utl_abort( myName//': PROBLEM WITH OER.')
235 end if
236 IF ( ZFGE < 1.0d-3 ) then
237 write(*,*)' Problem for GP STNID FGE= ', stnid, ZFGE
238 ZFGE=1.0d-3
239 end if
240 if ( zomper > 0.0d0 .and. zomper < 1.0d-3 ) zomper = 1.0d-3
241 else if ( obsFamily == 'CH' ) then
242 if ( ZFGE**2 + ZOER**2 < 1.0d-60) then
243 write(*,*) ' Problem for STNID FGE ZOER=', stnid, ZFGE, ZOER
244 ZFGE=1.0d30
245 ZOER=1.0d-30
246 end if
247 if ( zomper > 0.0 .and. zomper < 1.0d-30 ) zomper = 1.0d-30
248 else
249 if ( zfge < 0.0d0 ) zfge = 0.0d0
250 if ( ZFGE**2 + ZOER**2 < 1.0d-5)then
251 write(*,*) ' Problem for STNID FGE ZOER=', stnid, ZFGE, ZOER
252 ZFGE=1.0d-5
253 ZOER=1.0d-5
254 end if
255 if ( zomper > 0.0 .and. zomper < 1.0D-5 ) zomper = 1.0d-5
256 end if
257
258 ! Calculate zbgchk
259
260 if ( zomper > 0.0d0 ) then
261 zbgchk = (zomp**2)/(zomper**2)
262 zsop = zomper
263 else
264 ZBGCHK=(ZOMP)**2/(ZFGE**2 + ZOER**2)
265 ZSOP = SQRT(ZFGE**2 + ZOER**2)
266 end if
267
268 ! UPDATE QUALITY CONTROL FLAGS, based on zbgchk
269 ! ( ELEMENT FLAGS + GLOBAL HEADER FLAGS)
270 obsVarno = obs_bodyElem_i( obsData, OBS_VNM, bodyIndex )
271 if( obsVarno == bufr_nees .and. obs_bodyElem_i( obsData, OBS_XTR, bodyIndex ) == 0 ) then
272 averageFGE = averageFGE + zfge * zfge
273 averageOER = averageOER + zoer * zoer
274 obsCount = obsCount + 1
275 end if
276 codeType = obs_headElem_i( obsData, OBS_ITY, headerIndex )
277 obsFlag = ISETFLAG( obsFamily, codeType, obsVarno, ZBGCHK )
278
279 ! CONVERT ZTD VALUES FROM M TO MM FOR PRINTOUT
280
281 LLZD = .FALSE.
282 if ( obsFamily == 'GP' .and. obsVarno == BUFR_NEZD ) then
283 ZVAR = ZVAR * 1000.0D0
284 ZOER = ZOER * 1000.0D0
285 ZFGE = ZFGE * 1000.0D0
286 ZOMP = ZOMP * 1000.0D0
287 ZSOP = ZSOP * 1000.0D0
288 LLZD = .TRUE.
289 end if
290
291 stnid = obs_elem_c( obsData, 'STID', headerIndex )
292 if ( obsFlag >= 2 .or. ( LLZD .and. gps_gb_LTESTOP )) then
293 if ( obsFamily /= 'CH' ) then
294 write(*,122) &
295 stnid, zlat, zlon, codeType, obsVarno, ZLEV, ZVAR, ZOER, &
296 ZFGE, ZOMP, ZSOP, ZBGCHK, obsFlag
297 else
298 write(*,124) &
299 stnid, zlat, zlon, codeType, obsVarno, ZLEV, ZVAR, ZOER, &
300 ZFGE, ZOMP, ZSOP, ZBGCHK, obsFlag
301 end if
302 if ( obsFlag >= 2 ) numberObsRejected = numberObsRejected + 1
303 if ( obsFlag >= 2 .and. obsFamily == 'GP' ) then
304 if ( obsVarno == BUFR_NEZD ) INZREJ = INZREJ+1
305 if ( obsVarno == BUFR_NEPS ) INPREJ = INPREJ+1
306 if ( obsVarno == BUFR_NETS ) INTREJ = INTREJ+1
307 if ( obsVarno == BUFR_NESS ) INDREJ = INDREJ+1
308 end if
309 end if
310
311 if ( obsFlag == 1 ) then
312 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex ), 13 ))
313 else if ( obsFlag == 2 ) then
314 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex ), 14 ))
315 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex ), 16 ))
316 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex ), 09 ))
317 call obs_headSet_i( obsData, OBS_ST1, headerIndex, ibset( obs_headElem_i( obsData, OBS_ST1, headerIndex ), 06 ))
318
319 else if ( obsFlag == 3 ) then
320 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex), 15 ))
321 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex), 16 ))
322 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex), 09 ))
323 call obs_headSet_i( obsData, OBS_ST1, headerIndex, ibset( obs_headElem_i( obsData, OBS_ST1, headerIndex), 06 ))
324 end if
325 end if
326 end do BODY
327
328124 FORMAT(2x,a9,1x,f6.2,1x,f7.2,1x,I3,1x,I5,7(2x,G11.2),I3 )
329122 FORMAT(2x,a9,1x,f6.2,1x,f7.2,1x,I3,1x,I5,7(2x,F11.2),I3 )
330
331 end do HEADER
332
333 else !IF (.not. new_bgck_sw .or. obsFamily .ne. 'SW') then
334
335 call obs_set_current_body_list( obsData, 'SW' )
336 bodyuv: do
337 bodyIndex = obs_getBodyIndex(obsData)
338 if ( bodyIndex < 0 ) exit bodyuv
339
340 ! Only process pressure level observations flagged to be assimilated
341 i_ass = obs_bodyElem_i ( obsData, OBS_ASS, bodyIndex )
342 i_vco = obs_bodyElem_i ( obsData, OBS_VCO, bodyIndex )
343
344 if( i_ass /= obs_assimilated .or. i_vco /= 2 ) cycle bodyuv
345
346 headerIndex = obs_bodyElem_i( obsData, OBS_HIND, bodyIndex )
347 bodyIndex_start = obs_headElem_i( obsData, OBS_RLN, headerIndex )
348
349 obsVarno = obs_bodyElem_i( obsData, OBS_VNM, bodyIndex )
350 zlev = obs_bodyElem_r( obsData, OBS_PPP, bodyIndex )
351
352 found_u = .false.
353 found_v = .false.
354
355 if ( obsVarno == BUFR_NEUU ) then
356
357 bodyIndex_u = bodyIndex
358 found_u = .true.
359
360 do i_oth = bodyIndex_start, bodyIndex
361 obsVarno = obs_bodyElem_i( obsData, OBS_VNM, i_oth )
362 zslev = obs_bodyElem_r( obsData, OBS_PPP, i_oth )
363 if ( obsVarno == BUFR_NEVV .and. zslev == zlev ) then
364 bodyIndex_v = i_oth
365 found_v = .true.
366 end if
367 end do
368
369 else
370
371 bodyIndex_v = bodyIndex
372 found_v = .true.
373
374 do i_oth = bodyIndex_start, bodyIndex
375 obsVarno = obs_bodyElem_i( obsData, OBS_VNM, i_oth )
376 zslev = obs_bodyElem_r( obsData, OBS_PPP, i_oth )
377 if ( obsVarno == BUFR_NEUU .and. zslev == zlev ) then
378 bodyIndex_u = i_oth
379 found_u = .true.
380 end if
381 end do
382
383 end if
384
385 if ( found_u .and. found_v ) then
386
387 uu_d = obs_bodyElem_r( obsData, OBS_OMP, bodyIndex_u )
388 vv_d = obs_bodyElem_r( obsData, OBS_OMP, bodyIndex_v )
389 uu_r = obs_bodyElem_r( obsData, OBS_OER, bodyIndex_u )
390 vv_r = obs_bodyElem_r( obsData, OBS_OER, bodyIndex_v )
391 uu_f = obs_bodyElem_r( obsData, OBS_HPHT, bodyIndex_u )
392 vv_f = obs_bodyElem_r( obsData, OBS_HPHT, bodyIndex_v )
393
394 duv2 = uu_d**2 + vv_d**2
395
396 obsFlag = 0
397 duv2_lim = (uu_r**2 + uu_f**2 + vv_r**2 + vv_f**2)*1
398 if (duv2 > duv2_lim) then
399 obsFlag = 1
400 end if
401 duv2_lim = (uu_r**2 + uu_f**2 + vv_r**2 + vv_f**2)*2
402 if(duv2 > duv2_lim) then
403 obsFlag = 2
404 end if
405 duv2_lim = (uu_r**2 + uu_f**2 + vv_r**2 + vv_f**2)*4
406 if(duv2 > duv2_lim) then
407 obsFlag = 3
408 end if
409
410 numberObs = numberObs + 2
411 if ( obsFlag >= 2 ) numberObsRejected = numberObsRejected + 2
412
413 if ( obsFlag == 1 ) then
414 call obs_bodySet_i( obsData,OBS_FLG,bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 13 ))
415 call obs_bodySet_i( obsData,OBS_FLG,bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 13 ))
416
417 else if ( obsFlag == 2 ) then
418 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 14 ))
419 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 16 ))
420 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 09 ))
421 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 14 ))
422 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 16 ))
423 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 09 ))
424 call obs_headSet_i( obsData, OBS_ST1, headerIndex, ibset( obs_headElem_i( obsData, OBS_ST1, headerIndex ), 06 ))
425
426 else if ( obsFlag == 3 ) then
427 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 15 ))
428 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 16 ))
429 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_u, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_u ), 09 ))
430 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 15 ))
431 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 16 ))
432 call obs_bodySet_i( obsData, OBS_FLG, bodyIndex_v, ibset( obs_bodyElem_i( obsData, OBS_FLG, bodyIndex_v ), 09 ))
433 call obs_headSet_i( obsData, OBS_ST1, headerIndex, ibset( obs_headElem_i( obsData, OBS_ST1, headerIndex ), 06 ))
434 end if
435
436 end if !if(found_u .and. found_v)
437
438 end do bodyuv
439
440 end if !IF (.not. new_bgck_sw .or. obsFamily .ne. 'SW')
441
442 if ( numberObs > 0 ) then
443 write(*,*)' '
444 write(*,*) myName//': FINISHED BGCHECK OF ', obsFamily, ' DATA'
445 write(*,123) myName//': ', numberObsRejected, ' OBSERVATIONS REJECTED OUT OF ', numberObs
446 write(*,*)' '
447 end if
448
449 if ( numberObs > 0 .and. obsFamily == 'GP' ) then
450 write(*,*)' '
451 write(*,*) ' BGCDATA: REPORT FOR GP FAMILY OF OBSERVATIONS'
452 write(*,123) 'BGCDATA: ',INZREJ, ' ZTD OBSERVATIONS REJECTED OUT OF ', INZOBS
453 write(*,123) 'BGCDATA: ',INPREJ, ' PSFC OBSERVATIONS REJECTED OUT OF ', INPOBS
454 write(*,123) 'BGCDATA: ',INTREJ, ' TSFC OBSERVATIONS REJECTED OUT OF ', INTOBS
455 write(*,123) 'BGCDATA: ',INDREJ, ' DPDS OBSERVATIONS REJECTED OUT OF ', INDOBS
456 write(*,*)' '
457 end if
458
459123 FORMAT(2X,A,I0,A,I0)
460
461 write(*,*)' '
462 write(*,*)' ---------------------------'
463 write(*,*) myName//' DONE '
464 write(*,*)' ---------------------------'
465 write(*,*)' '
466
467 if ( obsCount > 0) then
468 write(*,*) myName//': obsCount: ', obsCount,'; mean(FGE): ', averageFGE / obsCount
469 write(*,*) myName//': obsCount: ', obsCount,'; mean(OER): ', averageOER / obsCount
470 end if
471
472 end subroutine bgck_data
473
474 !--------------------------------------------------------------------------
475 ! bgck_gpsro
476 !--------------------------------------------------------------------------
477 subroutine bgck_gpsro(columnTrlOnTrlLev,obsData)
478 !
479 !:Purpose: Set background-check flag on GPSRO data if ABS(O-P)/P is too
480 ! large
481 !
482 implicit none
483
484 ! Arguments:
485 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
486 type(struct_obs), intent(inout) :: obsData
487
488 ! Locals:
489 type(struct_vco), pointer :: vco_trl
490 real(8) :: HNH1, ZOBS, ZMHX, ZOMF, ZREF, ZOER, Rad
491 integer :: headerIndex
492 integer :: IDATYP, iProfile, varNum
493 integer :: IDATA , IDATEND, bodyIndex
494 integer :: NGPSLEV
495 integer :: iversion
496
497 vco_trl => col_getVco(columnTrlOnTrlLev)
498 iversion = vco_trl%vCode
499
500 write(*,*)'ENTER BGCSGPSRO'
501
502 ! 1. Initializations
503 ! ---------------
504
505 NGPSLEV=col_getNumLev(columnTrlOnTrlLev,'TH')
506
507 ! Loop over all files
508
509 ! loop over all header indices of the 'RO' family
510 call obs_set_current_header_list( obsData, 'RO' )
511 HEADER: do
512 headerIndex = obs_getHeaderIndex(obsData)
513 if (headerIndex < 0) exit HEADER
514
515 ! Process only refractivity data (codtyp 169)
516
517 IDATYP = obs_headElem_i(obsData,OBS_ITY,headerIndex)
518 IF ( IDATYP == 169 ) THEN
519 iProfile = gps_iprofile_from_index(headerIndex)
520 varNum = gps_vRO_IndexPrf(iProfile, 2)
521
522 ! Basic geometric variables of the profile
523
524 Rad = obs_headElem_r(obsData,OBS_TRAD,headerIndex)
525
526 ! Loops over data in the observation
527
528 IDATA = obs_headElem_i(obsData,OBS_RLN,headerIndex)
529 IDATEND = obs_headElem_i(obsData,OBS_NLV,headerIndex) + IDATA - 1
530
531 ! Scan for requested assimilations, and count them
532
533 do bodyIndex= IDATA, IDATEND
534 if ( obs_bodyElem_i( obsData, OBS_ASS, bodyIndex ) == obs_assimilated ) then
535 HNH1 = obs_bodyElem_r( obsData, OBS_PPP, bodyIndex )
536 if ( varNum == bufr_nebd ) HNH1 = HNH1-Rad
537
538 ! Increment OMF = Y - H(x)
539
540 ZOMF = obs_bodyElem_r( obsData, OBS_OMP, bodyIndex )
541
542 ! Observation value Y
543
544 ZOBS = obs_bodyElem_r( obsData, OBS_VAR, bodyIndex )
545 ZOER = obs_bodyElem_r( obsData, OBS_OER, bodyIndex )
546 ZMHX = ZOBS-ZOMF
547
548 ! Reference order of magnitude value:
549
550 if ( varNum == bufr_nebd ) then
551 ZREF = 0.025d0*exp(-HNH1/6500.d0)
552 else
553 if ( .not. gps_roBNorm ) then
554 ZREF = 300.d0*exp(-HNH1/6500.d0)
555 else
556 ZREF = ZMHX
557 end if
558 end if
559
560 ! OMF Tested criteria:
561
562 if ( .not. gps_roBNorm ) then
563 if (DABS(ZOMF)/ZREF > gps_BgckBand .or. DABS(ZOMF)/ZOER > 3.d0) then
564 call obs_bodySet_i(obsData,OBS_FLG,bodyIndex,ibset(obs_bodyElem_i(obsData,OBS_FLG,bodyIndex),16))
565 call obs_bodySet_i(obsData,OBS_FLG,bodyIndex,ibset(obs_bodyElem_i(obsData,OBS_FLG,bodyIndex),9))
566 end if
567 else
568 if ( DABS(ZOMF)/ZREF > gps_BgckBand .or. DABS(ZOMF)/ZOER > gps_roNsigma) then
569 call obs_bodySet_i(obsData,OBS_FLG,bodyIndex,ibset(obs_bodyElem_i(obsData,OBS_FLG,bodyIndex),16))
570 call obs_bodySet_i(obsData,OBS_FLG,bodyIndex,ibset(obs_bodyElem_i(obsData,OBS_FLG,bodyIndex),9))
571 write(*,'(A40,F10.0,3F12.4)') ' REJECT BGCSGPSRO H O P (O-P/ZREF) =',HNH1,ZOBS,ZMHX,(ZOMF)/ZREF
572 end if
573 end if
574
575 end if
576 end do
577 end if
578
579 end do HEADER
580
581 write(*,*)'EXIT BGCSGPSRO'
582 RETURN
583
584 end subroutine bgck_gpsro
585
586 !--------------------------------------------------------------------------
587 ! isetflag
588 !--------------------------------------------------------------------------
589 function isetflag( obsFamily, kodtyp, kvnam, zbgchk )
590 !
591 !:Purpose: Set BACKGROUND-CHECK FLAGS According to values set in a table.
592 ! Original values in table come from ecmwf.
593 !
594 implicit none
595
596 ! Arguments:
597 character(len=2), intent(in) :: obsFamily ! FAMILY NAME ( 'UA' , 'AI' ...etc.. )
598 integer, intent(in) :: kodtyp ! BURP CODE TYPE
599 integer, intent(in) :: kvnam ! VARIABLE NAME ( BURP )
600 real(8), intent(in) :: zbgchk ! NORMALIZED BACKGROUND DEPARTURE
601 ! Result:
602 integer :: isetflag
603
604 ! Locals:
605 real(8), parameter :: zsacrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
606 real(8), parameter :: zttcrit(3) = (/ 9.00D0, 16.00D0, 25.00D0 /)
607 real(8), parameter :: zalcrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
608 real(8), parameter :: zescrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
609 real(8), parameter :: zpscrit(3) = (/ 9.00D0, 16.00D0, 25.00D0 /)
610 real(8), parameter :: zpncrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
611 real(8), parameter :: zswcrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
612 real(8), parameter :: ztscrit(3) = (/ 5.00D0, 25.00D0, 30.00D0 /)
613 real(8), parameter :: zdzcrit(3) = (/ 2.25D0, 5.06D0, 7.56D0 /)
614 real(8), parameter :: zgzcrit(3) = (/ 12.25D0, 25.00D0, 36.00D0 /)
615 real(8), parameter :: zzdcrit(3) = (/ 9.00D0, 16.00D0, 25.00D0 /)
616 real(8), parameter :: zchcrit(3) = (/ 9.00D0, 16.00D0, 25.00D0 /)
617 real(8), parameter :: zLogViscrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
618 ! Temporary hardcoded values for radar Doppler velocity
619 real(8), parameter :: radvelcrit(3) = (/ 8.00D0, 20.00D0, 30.00D0 /)
620 !real(8) :: zuvcrit(3) = (/ 10.00D0, 20.00D0, 30.00D0 /)
621 real(8) :: zuvcrit(3)
622
623 zuvcrit(1) = 10.00D0
624 zuvcrit(3) = 30.00D0
625 if ( kodtyp == 37 ) then
626 zuvcrit(2)=25.00D0
627 else
628 zuvcrit(2)=20.00D0
629 end if
630
631 isetflag=0
632
633 if ( kvnam == BUFR_NEGZ ) then
634
635 ! SET FLAG FOR HEIGHTS
636
637 if ( zbgchk >= zgzcrit(1) .and. zbgchk < zgzcrit(2) ) then
638 isetflag=1
639 else if ( zbgchk >= zgzcrit(2) .and. zbgchk < zgzcrit(3) ) then
640 isetflag=2
641 else if ( zbgchk >= zgzcrit(3) )then
642 isetflag =3
643 end if
644
645 else if ( kvnam == BUFR_NETT ) then
646
647 ! SET FLAG FOR TEMPERATURE
648
649 if ( zbgchk >= zttcrit(1) .and. zbgchk < zttcrit(2) ) then
650 isetflag=1
651 else if ( zbgchk >= zttcrit(2) .and. zbgchk < zttcrit(3) ) then
652 isetflag=2
653 else if ( zbgchk >= zttcrit(3) )then
654 isetflag =3
655 end if
656
657 else if ( kvnam == BUFR_NEDZ ) then
658
659 ! SET FLAG FOR SATEMS
660
661 if ( zbgchk >= zdzcrit(1) .and. zbgchk < zdzcrit(2) ) then
662 isetflag=1
663 else if ( zbgchk >= zdzcrit(2) .and. zbgchk < zdzcrit(3) ) then
664 isetflag=2
665 else if ( zbgchk >= zdzcrit(3) )then
666 isetflag =3
667 end if
668
669 else if ( kvnam == BUFR_NEFS ) then
670
671 ! SET FLAG FOR WIND SPEED
672
673 if ( zbgchk >= zsacrit(1) .and. zbgchk < zsacrit(2) ) then
674 isetflag=1
675 else if ( zbgchk >= zsacrit(2) .and. zbgchk < zsacrit(3) ) then
676 isetflag=2
677 else if ( zbgchk >= zsacrit(3) )then
678 isetflag =3
679 end if
680
681 else if ( kvnam == BUFR_NEUU .or. kvnam == BUFR_NEVV ) then
682
683 ! SET FLAG FOR WIND COMPONENTS
684
685 if ( zbgchk >= zuvcrit(1) .and. zbgchk < zuvcrit(2) ) then
686 isetflag=1
687 else if ( zbgchk >= zuvcrit(2) .and. zbgchk < zuvcrit(3) ) then
688 isetflag=2
689 else if ( zbgchk >= zuvcrit(3) )then
690 isetflag =3
691 end if
692
693 else if ( kvnam == BUFR_NEAL ) then
694
695 ! SET FLAG FOR ALADIN HLOS WIND OBSERVATIONS
696
697 if ( zbgchk >= zalcrit(1) .and. zbgchk < zalcrit(2) ) then
698 isetflag=1
699 else if ( zbgchk >= zalcrit(2) .and. zbgchk < zalcrit(3) ) then
700 isetflag=2
701 else if ( zbgchk >= zalcrit(3) )then
702 isetflag=3
703 end if
704
705 else if ( kvnam == BUFR_NEUS .or. kvnam == BUFR_NEVS ) then
706
707 ! SET FLAG FOR SURFACE WIND COMPONENTS
708
709 if ( zbgchk >= zswcrit(1) .and. zbgchk < zswcrit(2) ) then
710 isetflag=1
711 else if ( zbgchk >= zswcrit(2) .and. zbgchk < zswcrit(3) ) then
712 isetflag=2
713 else if ( zbgchk >= zswcrit(3) )then
714 isetflag =3
715 end if
716
717 else if ( kvnam == bufr_gust ) then
718
719 ! SET FLAG FOR SURFACE WIND GUST
720
721 if ( zbgchk >= zswcrit(1) .and. zbgchk < zswcrit(2) ) then
722 isetflag=1
723 else if ( zbgchk >= zswcrit(2) .and. zbgchk < zswcrit(3) ) then
724 isetflag=2
725 else if ( zbgchk >= zswcrit(3) )then
726 isetflag =3
727 end if
728
729 else if ( kvnam == bufr_nees ) then
730
731 ! SET FLAG FOR DEW POINT DEPRESSION
732
733 if ( zbgchk >= zescrit(1) .and. zbgchk < zescrit(2) ) then
734 isetflag=1
735 else if ( zbgchk >= zescrit(2) .and. zbgchk < zescrit(3) ) then
736 isetflag=2
737 else if ( zbgchk >= zescrit(3) )then
738 isetflag =3
739 end if
740
741 else if ( kvnam == bufr_neps ) then
742
743 ! SET FLAG FOR SURFACE PRESSURE
744
745 if ( zbgchk >= zpscrit(1) .and. zbgchk < zpscrit(2) ) then
746 isetflag=1
747 else if ( zbgchk >= zpscrit(2) .and. zbgchk < zpscrit(3) ) then
748 isetflag=2
749 else if ( zbgchk >= zpscrit(3) )then
750 isetflag =3
751 end if
752
753 else if ( kvnam == bufr_nepn ) then
754
755 ! SET FLAG FOR MEAN SEA LEVEL PRESSURE
756
757 if ( zbgchk >= zpncrit(1) .and. zbgchk < zpncrit(2) ) then
758 isetflag=1
759 else if ( zbgchk >= zpncrit(2) .and. zbgchk < zpncrit(3) ) then
760 isetflag=2
761 else if ( zbgchk >= zpncrit(3) ) then
762 isetflag =3
763 end if
764
765 else if ( kvnam == bufr_nets ) then
766
767 ! SET FLAG FOR SURFACE TEMPERATURE
768
769 if ( zbgchk >= ztscrit(1) .and. zbgchk < ztscrit(2) ) then
770 isetflag=1
771 else if ( zbgchk >= ztscrit(2) .and. zbgchk < ztscrit(3) ) then
772 isetflag=2
773 else if ( zbgchk >= ztscrit(3) ) then
774 isetflag =3
775 end if
776
777 else if ( kvnam == BUFR_NESS ) then
778
779 ! SET FLAG FOR SURFACE DEW POINT DEPRESSION
780
781 if ( zbgchk >= zescrit(1) .and. zbgchk < zescrit(2) ) then
782 isetflag=1
783 else if ( zbgchk >= zescrit(2) .and. zbgchk < zescrit(3) ) then
784 isetflag=2
785 else if ( zbgchk >= zescrit(3) ) then
786 isetflag =3
787 end if
788
789 else if ( kvnam == bufr_vis ) then
790
791 ! SET FLAG FOR VISIBILITY
792
793 if ( zbgchk >= zLogViscrit(1) .and. zbgchk < zLogViscrit(2) ) then
794 isetflag=1
795 else if ( zbgchk >= zLogViscrit(2) .and. zbgchk < zLogViscrit(3) ) then
796 isetflag=2
797 else if ( zbgchk >= zLogViscrit(3) ) then
798 isetflag =3
799 end if
800
801 else if ( kvnam == bufr_nezd ) then
802
803 ! SET FLAG FOR GB-GPS ZENITH DELAY
804
805 if ( zbgchk >= zzdcrit(1) .and. zbgchk < zzdcrit(2) ) then
806 isetflag=1
807 else if ( zbgchk >= zzdcrit(2) .and. zbgchk < zzdcrit(3) ) then
808 isetflag=2
809 else if ( zbgchk >= zzdcrit(3) ) then
810 isetflag =3
811 end if
812
813 else if ( kvnam == bufr_radvel ) then
814
815 ! Set flag for Doppler Velocity (Radvel)
816 if ( zbgchk >= radvelcrit(1) .and. zbgchk < radvelcrit(2) ) then
817 isetflag = 1
818 else if ( zbgchk >= radvelcrit(2) .and. zbgchk < radvelcrit(3) ) then
819 isetflag = 2
820 else if ( zbgchk >= radvelcrit(3) ) then
821 isetflag = 3
822 end if
823
824 else if ( obsFamily == 'CH' ) then
825
826 ! SET FLAG FOR CHEMICAL CONSTITUENTS
827 if ( zbgchk >= zchcrit(1) .and. zbgchk < zchcrit(2) ) then
828 isetflag=1
829 else if ( zbgchk >= zchcrit(2) .and. zbgchk < zchcrit(3) ) then
830 isetflag=2
831 else if ( zbgchk >= zchcrit(3) ) then
832 isetflag =3
833 end if
834
835 end if
836
837 return
838
839 end function isetflag
840
841end module backgroundCheck_mod
