1module varQC_mod
2 ! MODULE varQC_mod (prefix='vqc' category='1. High-level functionality')
3 !
4 !:Purpose: Procedures related to variational quality control including
5 ! hard-coded values that determine how quickly the observation
6 ! weight begins to be reduced.
7 !
8 use MathPhysConstants_mod
9 use earthConstants_mod
10 use codtyp_mod
11 use bufr_mod
12 use obsSpaceData_mod
13 use columnData_mod
14 use rmatrix_mod ,only : rmat_lnondiagr
15 use varNameList_mod
16 use obsFamilyList_mod
17
18 implicit none
19 save
20 private
21
22 ! public procedures
23 public :: vqc_setup, vqc_NlTl, vqc_ad, vqc_listrej
24
25
26 contains
27
28 subroutine vqc_setup(obsSpaceData)
29 !
30 !:Purpose: To set certain parameters for the asymmetric check
31 ! and for variational quality control
32 !
33 implicit none
34
35 ! Arguments:
36 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
37
38 ! Locals:
39 integer jdata, jjo, idata, idatend, idburp
40 integer ityp, iass, iother, jj, istyp, ilev
41 real(8) zagz, zahu, zatt, zauv, zabt, zabtb, zapn, zaps, zazd, zach, zatm, zaice, zapr
42 real(8) zdgz, zdhu, zdtt, zduv, zdbt, zdbtb, zdpn, zdps, zdzd, zdch, zdtm, zdice, zdpr
43 real(8) zattra, zauvra, zattsym, zdvis, zavis
44 real(8) zslev, zlev, zval, zspdo, zspdf, zofcst, zoval, zdiff, zaasym, zoer
45 real(8) zfcst, zlat, zlon, zprior
46 logical llok
47
48 !
49 !_____prior probabilities for winds:ZAUV
50 ! prior probabilities for scalar variables: zagz and zahu
51 ! standard deviation multiple for background check for winds: zduv
52 ! standard deviation multiple for background check for heights: zdgz
53 ! standard deviation multiple for background check for humidity: zdhu
54 !
55 zagz = 1.0d-12
56! zahu = 1.0d-2
57 zahu = 5.0d-2
58 zatt = 5.0d-2
59! zauv = 1.0d-2
60 zauv = 2.0d-2
61 zabt = 1.0d-1
62 zabtb = 1.0d-1
63 zapn = 1.0d-4
64 zaps = 1.0d-4
65 zazd = 2.0d-2
66 zach = 1.0d-3
67 zatm = 1.0d-12
68 zaice = 0.0d0
69 zapr = 0.0d0
70
71 zattra = 0.005d0
72 zattsym = 1.d-1
73 zauvra = 1.d-5
74 zavis= 1.0d-3
75
76 zdgz = 5.d0
77 zdhu = 5.d0
78 zdtt = 5.d0
79 zduv = 5.d0
80 zdbt = 3.d0
81 zdbtb = 3.d0
82 zdpn = 5.d0
83 zdps = 5.d0
84 zdzd = 3.d0
85 zdch = 1.d1
86 zdtm = 5.d0
87 zdice = 5.d2
88 zdpr = 5.d2
89 zdvis= 5.d0
90
91 do jjo = 1, obs_numheader(obsSpaceData)
92
93 IDATA = obs_headElem_i( obsSpaceData, OBS_RLN, jjo )
94 IDATEND = obs_headElem_i( obsSpaceData, OBS_NLV, jjo ) + IDATA - 1
95 IDBURP = obs_headElem_i( obsSpaceData, OBS_ITY, jjo )
96 ZLAT = obs_headElem_r( obsSpaceData, OBS_LAT, jjo ) * MPC_DEGREES_PER_RADIAN_R8
97 ZLON = obs_headElem_r( obsSpaceData, OBS_LON, jjo ) * MPC_DEGREES_PER_RADIAN_R8
98
99 do JDATA = IDATA, IDATEND
100
101 ityp = obs_bodyElem_i( obsSpaceData, OBS_VNM, JDATA )
102 IASS = obs_bodyElem_i( obsSpaceData, OBS_ASS, JDATA )
103 ZLEV = obs_bodyElem_r( obsSpaceData, OBS_PPP, JDATA ) * MPC_MBAR_PER_PA_R8
104 ZOER = obs_bodyElem_r( obsSpaceData, OBS_OER, JDATA )
105 ZVAL = obs_bodyElem_r( obsSpaceData, OBS_VAR, JDATA )
106
107 ZFCST = ZVAL - obs_bodyElem_r( obsSpaceData, OBS_OMP,JDATA)
108
109 if (ityp == BUFR_NETS .or. ityp == BUFR_NEPS .or. &
110 ityp == BUFR_NEPN .or. ityp == BUFR_NESS .or. &
111 ityp == BUFR_NEUS .or. ityp == BUFR_NEVS .or. &
112 ityp == BUFR_NEZD .or. ityp == bufr_vis .or. &
113 ityp == bufr_logVis .or. ityp == bufr_gust .or. &
114 ityp == bufr_radarPrecip .or. ityp == bufr_logRadarPrecip ) then
115 LLOK = (obs_bodyElem_i(obsSpaceData,OBS_ASS,JDATA) == obs_assimilated)
116 else
117 LLOK = (IASS == obs_assimilated) .and. ((obs_bodyElem_i(obsSpaceData,OBS_XTR,JDATA) ==0) &
118 .or. ((obs_bodyElem_i(obsSpaceData,OBS_XTR,JDATA) == 2) .and. &
119 (obs_bodyElem_i(obsSpaceData,OBS_VNM,JDATA) == BUFR_NEGZ)))
120 end if
121
122 if (LLOK) then
123 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEVV .or. &
124 ityp == BUFR_NEUS .or. ityp == BUFR_NEVS) then
125 ZAASYM = 1.0d0
126 IOTHER = -1
127 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEUS) then
128 do JJ = IDATA, JDATA
129 ISTYP = obs_bodyElem_i( obsSpaceData, OBS_VNM, JJ )
130 ZSLEV = obs_bodyElem_r( obsSpaceData, OBS_PPP, JJ ) * MPC_MBAR_PER_PA_R8
131 if ((ISTYP == BUFR_NEVV .or. ISTYP == BUFR_NEVS) &
132 .and. ZLEV == ZSLEV) then
133 IOTHER = JJ
134 end if
135 end do
136 else
137 do JJ=IDATA,JDATA
138 ISTYP = obs_bodyElem_i(obsSpaceData,OBS_VNM,JJ)
139 ZSLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JJ)*MPC_MBAR_PER_PA_R8
140 if ((ISTYP == BUFR_NEUU .or. ISTYP == BUFR_NEUS) &
141 .and. ZLEV == ZSLEV) then
142 IOTHER = JJ
143 end if
144 end do
145 end if
146 if (IOTHER /= -1) then
147 ZOER = obs_bodyElem_r(obsSpaceData,OBS_OER,JDATA)
148 ZOVAL = obs_bodyElem_r(obsSpaceData,OBS_VAR,IOTHER)
149 ZOFCST = ZOVAL-obs_bodyElem_r(obsSpaceData,OBS_OMP,IOTHER)
150 ZSPDO = SQRT(ZOVAL*ZOVAL + ZVAL*ZVAL)
151 ZSPDF = SQRT(ZOFCST*ZOFCST + ZFCST*ZFCST)
152 ZDIFF = ZSPDO - ZSPDF
153 ILEV = NINT(ZLEV)
154 !
155 !___tighten rejection criterion for satob winds
156 !
157 if (IDBURP == 88 .or. IDBURP == 188) then
158 if (ZDIFF < 0.0d0 .and. ABS(ZLAT) > 20.d0 .and. &
159 ILEV < 550) then
160 ZAASYM = 0.7d0*MAX(ZSPDF,1.0d0)
161 ZPRIOR = ZAASYM*ZAUV
162 if (ZPRIOR > 0.99d0) ZAASYM = 0.99d0/ZAUV
163 else
164 !
165 !___zaasym used to specify new criterion for satwinds
166 ! than are not included in the asymmetric test
167 !
168 ZAASYM = 10.d0
169 end if
170 end if
171
172 end if
173 !
174 !__INITIAL VALUE OF GAMMA FOR QCVAR (WINDS)
175 !
176 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(1.d0 - &
177 (1.d0-(ZAUV*ZAASYM))*(1.d0-(ZAUV*ZAASYM)))*(2.d0*MPC_PI_r8)/ &
178 ((1.d0-(ZAUV*ZAASYM))*(1.d0-(ZAUV*ZAASYM))* &
179 (2.d0*ZDUV)*(2.d0*ZDUV)))
180 if ((IDBURP >= 32 .and. IDBURP <= 38) .or. &
181 (IDBURP >= 135 .and. IDBURP <= 142) .or. &
182 (IDBURP == 132) ) &
183 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA, &
184 (1.d0 - (1.d0-ZAUVRA)*(1.d0-ZAUVRA)) &
185 *(2.d0*MPC_PI_R8)/((1.d0-ZAUVRA)*(1.d0-ZAUVRA)* &
186 (2.d0*ZDUV)*(2.d0*ZDUV)))
187 else if (ityp == BUFR_NEGZ) then
188 !
189 ! INITIAL VALUE OF GAMMA FOR QCVAR (HEIGHTS)
190 !
191 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZAGZ* &
192 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZAGZ)*(2.d0*ZDGZ)))
193 else if (ityp == BUFR_NETT) then
194 !
195 ! INITIAL VALUE OF GAMMA FOR QCVAR (TEMPERATURES)
196 !
197 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZATT* &
198 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZATT)*(2.d0*ZDTT)))
199 if ((IDBURP >= 32 .and. IDBURP <= 38) .or. &
200 (IDBURP >= 135 .and. IDBURP <= 142) .or. &
201 (IDBURP == 132) ) &
202 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZATTRA* &
203 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZATTRA)*(2.d0*ZDTT)))
204 else if (ityp == BUFR_NETS) then
205 !
206 ! INITIAL VALUE OF GAMMA FOR QCVAR (SCREEN-LEVEL TEMPERATURES)
207 !
208 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZATT* &
209 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZATT)*(2.d0*ZDTT)))
210 !
211 ! ASYMMETRIC TEST FOR SHIP TEMPERATURES
212 !
213 if ((IDBURP == 13) .and. (ZVAL > ZFCST)) then
214 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZATTSYM* &
215 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZATTSYM)*(2.d0*ZDTT)))
216 end if
217 else if (ityp == BUFR_NEPN) then
218 !
219 ! INITIAL VALUE OF GAMMA FOR QCVAR (MSL PRESSURE)
220 !
221 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZAPN* &
222 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZAPN)*(2.d0*ZDPN)))
223 else if (ityp == BUFR_NEPS) then
224 !
225 ! INITIAL VALUE OF GAMMA FOR QCVAR (STATION PRESSURE)
226 !
227 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZAPS* &
228 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZAPS)*(2.d0*ZDPS)))
229 else if ( ityp == 12062 .or. ityp == 12063 .or. &
230 ityp == 12163) then
231 !
232 ! INITIAL VALUE OF GAMMA FOR BRIGHTNESS TEMPERATURES
233 ! TOVS AMSU-A + AIRS + IASI + GEORAD !!!!!
234 !
235 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZABT* &
236 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZABT)*(2.d0*ZDBT)))
237 if (IDBURP == 181 .or. IDBURP == 182 .or. &
238 IDBURP == 200) then
239 !
240 ! INITIAL VALUE OF GAMMA FOR TOVS AMSU-B (181) AND MHS (182)
241 ! AND MWHS2(200)
242 !
243 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZABTB* &
244 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZABTB)*(2.d0*ZDBTB)))
245 end if
246 else if (ityp == BUFR_NEZD) then
247 !
248 ! INITIAL VALUE OF GAMMA FOR QCVAR (GPS ZENITH DELAY)
249 !
250 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZAZD* &
251 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZAZD)*(2.d0*ZDZD)))
252 else if (bufr_IsAtmosConstituent(ityp)) then
253 !
254 ! INITIAL VALUE OF GAMMA FOR THE CH (constituents) family
255 !
256 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZACH* &
257 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZACH)*(2.d0*ZDCH)))
258 else if (ityp == bufr_sst ) then
259 !
260 ! INITIAL VALUE OF GAMMA FOR SST
261 !
262 call obs_bodySet_r( obsSpaceData, OBS_POB, JDATA, ( zatm * &
263 sqrt( 2.d0 * MPC_PI_R8 )) / (( 1.d0 - zatm ) * ( 2.d0 * zdtm )))
264 else if (ityp == BUFR_ICEC ) then
265 !
266 ! INITIAL VALUE OF GAMMA FOR SEA ICE
267 !
268 call obs_bodySet_r( obsSpaceData, OBS_POB, JDATA, ( zaice * &
269 sqrt( 2.d0 * MPC_PI_R8 )) / (( 1.d0 - zaice ) * ( 2.d0 * zdice )))
270 else if (ityp == BUFR_radarPrecip .or. ityp == BUFR_logRadarPrecip ) then
271 !
272 ! INITIAL VALUE OF GAMMA FOR RADAR PRECIPITATION
273 !
274 call obs_bodySet_r( obsSpaceData, OBS_POB, JDATA, ( zapr * &
275 sqrt( 2.d0 * MPC_PI_R8 )) / (( 1.d0 - zapr ) * ( 2.d0 * zdpr )))
276 else if (ityp == bufr_gust ) then
277 !
278 ! INITIAL VALUE OF GAMMA FOR WIND GUST
279 !
280 call obs_bodySet_r( obsSpaceData, OBS_POB, JDATA, ( zauv * &
281 sqrt( 2.d0 * MPC_PI_R8 )) / (( 1.d0 - zauv ) * ( 2.d0 * zduv )))
282 else if (ityp == bufr_vis .or. ityp == bufr_logVis) then
283 !
284 ! INITIAL VALUE OF GAMMA FOR VISIBILITY
285 !
286 call obs_bodySet_r( obsSpaceData, OBS_POB, JDATA, ( zavis * &
287 sqrt( 2.d0 * MPC_PI_R8 )) / (( 1.d0 - zavis ) * ( 2.d0 * zdvis )))
288 else
289 !
290 ! INITIAL VALUE OF GAMMA FOR QCVAR (OTHERS)
291 !
292 call obs_bodySet_r(obsSpaceData,OBS_POB,JDATA,(ZAHU* &
293 SQRT(2.d0*MPC_PI_R8))/((1.d0-ZAHU)*(2.d0*ZDHU)))
294 end if
295 end if
296 end do
297
298 end do
299
300 end subroutine vqc_setup
301
302
303 subroutine vqc_NlTl(obsSpaceData)
304 !
305 !:Purpose: 1) Modify Jo [OBS_JOBS] according to
306 ! Andersson and Jarvinen 1999, Variational quality control,
307 ! Q.J.R., 125, pp. 697-722.
308 ! 2) Save the values of (1-Wqc) in OBS_QCV
309 ! for gradient factorization and postalt flag criterion.
310 !
311 implicit none
312
313 ! Arguments:
314 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
315
316 ! Locals:
317 integer :: index_body,istyp,jj,index_header,ityp,index_body_start
318 real*8 :: zgami,zjon,zqcarg,zppost,zlev,zslev
319 logical :: lluv
320 logical :: includeFlag
321
322 BODY: do index_body = 1, obs_numbody(obsSpaceData)
323 includeFlag = (obs_bodyElem_i(obsSpaceData,OBS_ASS,index_body) == obs_assimilated).and. &
324 (obs_getFamily(obsSpaceData,bodyIndex_opt=index_body).ne.'RO')
325 ! pas de qcvar pour les radiances en mode matrice R non diagonale
326 if (rmat_lnondiagr) includeFlag = includeFlag .and. &
327 (obs_getFamily(obsSpaceData,bodyIndex_opt=index_body) /= 'TO')
328
329 if (includeFlag) then
330 index_header = obs_bodyElem_i(obsSpaceData,OBS_HIND,index_body)
331 index_body_start = obs_headElem_i(obsSpaceData,OBS_RLN,INDEX_HEADER)
332 ZLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,INDEX_BODY)
333 zgami = obs_bodyElem_r(obsSpaceData,OBS_POB,index_body)
334 ityp = obs_bodyElem_i(obsSpaceData,OBS_VNM,INDEX_BODY)
335 LLUV = ((ityp == BUFR_NEUU .or. ityp == BUFR_NEUS) .and. &
336 col_varExist(varName='UU')) .or. ((ityp == BUFR_NEVV .or. &
337 ityp == BUFR_NEVS).and.col_varExist(varName='VV'))
338 if (LLUV) then
339 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEUS)then
340 !
341 ! In order to calculate the contribution to Jo from a wind, the o-a
342 ! must be available for both u and v components. Hence, loop over only
343 ! data for which o-a has already been calculated
344 !
345 do JJ=INDEX_BODY_START, INDEX_BODY
346 ISTYP = obs_bodyElem_i(obsSpaceData,OBS_VNM,JJ)
347 ZSLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JJ)
348 if ((ISTYP == BUFR_NEVV .or. &
349 ISTYP == BUFR_NEVS) .and. &
350 ZSLEV == ZLEV) then
351 ZJON=obs_bodyElem_r(obsSpaceData,OBS_JOBS,INDEX_BODY)+ &
352 obs_bodyElem_r(obsSpaceData,OBS_JOBS,JJ)
353 ZQCARG = ZGAMI + EXP(-1.0D0*ZJON)
354 ZPPOST = ZGAMI/ZQCARG
355 !
356 ! Store the value of o-a multiplied by one minus the posterior
357 ! probability of gross error (needed for the adjoint calculations)
358 !
359 call obs_bodySet_r(obsSpaceData,OBS_QCV,INDEX_BODY, ZPPOST)
360 call obs_bodySet_r(obsSpaceData,OBS_QCV,JJ, ZPPOST)
361
362 call obs_bodySet_r(obsSpaceData,OBS_JOBS,INDEX_BODY,-LOG(ZQCARG/(ZGAMI+1.D0))/2.D0)
363 call obs_bodySet_r(obsSpaceData,OBS_JOBS,JJ, -LOG(ZQCARG/(ZGAMI+1.D0))/2.D0)
364 end if
365 end do
366 else ! ityp
367 do JJ=INDEX_BODY_START, INDEX_BODY
368 ISTYP = obs_bodyElem_i(obsSpaceData,OBS_VNM,JJ)
369 ZSLEV = obs_bodyElem_r(obsSpaceData,OBS_PPP,JJ)
370 if ((ISTYP == BUFR_NEUU .or. &
371 ISTYP == BUFR_NEUS) .and. &
372 ZSLEV == ZLEV) then
373 ZJON=obs_bodyElem_r(obsSpaceData,OBS_JOBS,INDEX_BODY)+ &
374 obs_bodyElem_r(obsSpaceData,OBS_JOBS,JJ)
375 ZQCARG = ZGAMI + EXP(-1.0D0*ZJON)
376 ZPPOST = ZGAMI/ZQCARG
377 call obs_bodySet_r(obsSpaceData,OBS_QCV,INDEX_BODY, ZPPOST)
378 call obs_bodySet_r(obsSpaceData,OBS_QCV,JJ, ZPPOST)
379 call obs_bodySet_r(obsSpaceData,OBS_JOBS,INDEX_BODY,-LOG(ZQCARG/(ZGAMI+1.D0))/2.D0)
380 call obs_bodySet_r(obsSpaceData,OBS_JOBS,JJ, -LOG(ZQCARG/(ZGAMI+1.D0))/2.D0)
381 end if
382 enddo
383 endif !ityp
384 else ! LLUV
385 zjon = obs_bodyElem_r(obsSpaceData,OBS_JOBS,index_body)
386 zqcarg = zgami + exp(-1.0D0*zjon)
387 zppost = zgami/zqcarg
388 call obs_bodySet_r(obsSpaceData,OBS_QCV,index_body, zppost)
389 call obs_bodySet_r(obsSpaceData,OBS_JOBS,INDEX_BODY, - log(zqcarg/(zgami+1.D0)))
390 endif ! LLUV
391
392 endif ! includeFlag
393
394 enddo BODY
395
396 end subroutine vqc_NlTl
397
398
399 subroutine vqc_ad(obsSpaceData)
400 !
401 !:Purpose: Factorizes Grad(Jo) according to Andersson and Jarvinen
402 ! 1999, Variational quality control, Q.J.R., 125, pp. 697-722.
403 ! It uses the value of (1-Wqc) saved in OBS_QCV in vqc_NlTl
404 !
405 implicit none
406
407 ! Arguments:
408 type(struct_obs), intent(inout) :: obsSpaceData ! obsSpaceData object
409
410 ! Locals:
411 integer :: index_body
412 logical :: includeFlag
413
414 do index_body=1,obs_numbody(obsSpaceData)
415 includeFlag = (obs_bodyElem_i(obsSpaceData,OBS_ASS,index_body) == obs_assimilated) .and. &
416 (obs_getFamily(obsSpaceData,bodyIndex_opt=index_body) /= 'RO')
417 ! pas de qcvar pour les radiances en mode matrice R non diagonale
418 if (rmat_lnondiagr) includeFlag = includeFlag .and. &
419 (obs_getFamily(obsSpaceData,bodyIndex_opt=index_body) /= 'TO')
420
421 if (includeFlag) then
422 call obs_bodySet_r(obsSpaceData,OBS_WORK,index_body, &
423 obs_bodyElem_r(obsSpaceData,OBS_WORK,index_body) &
424 *(1.d0 - obs_bodyElem_r(obsSpaceData,OBS_QCV,index_body)))
425 endif
426 enddo
427
428 end subroutine vqc_ad
429
430
431 subroutine vqc_listrej(lobsSpaceData)
432 !
433 !:Purpose: List all observations rejected by the variational QC
434 ! Set QC flags consistent with VARQC decisions
435 ! Set global flag indicating report contains rejected observations
436 ! as required
437 !
438 implicit none
439
440 ! Arguments:
441 type(struct_obs), intent(inout) :: lobsSpaceData
442
443 ! Locals:
444 integer, parameter :: numitem = 16
445 integer :: icount( numitem, ofl_numFamily )
446 integer :: jfam, jitem, headerIndex
447 integer :: bodyIndex, bodyIndex2, ISTART,ityp,ISTYP
448 integer :: ISPDO,IDIRO,ISPDF,IDIRF,ISPDA,IDIRA,ILEV
449 integer :: IDBURP,IOTHER,obsONM
450 real(8) :: ZVAR,ZFCST,ZANA,ZPOST,ZLAT,ZLON,ZUU,ZVV
451 real(8) :: SPD,DEG,ZLEV,ZSLEV,ZCUT
452 character(len=4) :: CLITM(NUMITEM),CLDESC
453 character(len=2) :: CLUNITS
454 character(len=21) :: CODTYPNAME
455 character(len=12) :: stnId
456 logical :: LLOK,LLELREJ
457
458 ! ------NOTE----------
459 ! CURRENTLY SUPPORTED FAMILIES OF DATA 'UA' 'AI' 'SF' 'HU' 'TO' 'GO' 'GP'
460
461 DATA ZCUT /0.75D0/
462 DATA CLITM(1), CLITM(2), CLITM(3), CLITM(4), CLITM(5), CLITM(6) &
463 /'WND', 'WND', 'HGT', 'TMP', 'DPD', 'STNP'/
464 DATA CLITM(7), CLITM(8), CLITM(9), CLITM(10), CLITM(11), CLITM(12) &
465 /'MSLP', 'TSFC', 'SDPD', 'SWND', 'SWND', 'BTMP'/
466 DATA CLITM(13), CLITM(14), CLITM(15) , CLITM(16) &
467 /'ZTD', 'CHM', 'SST', 'ICEC'/
468
469
470 do jfam = 1, ofl_numFamily
471 do jitem = 1, numitem
472 icount( jitem, jfam ) = 0
473 end do
474 end do
475 write(*,*) 'LIST OF DATA REJECTED BY VARIATIONAL QUALITY CONTROL'
476 !
477 !_____LOOP OVER ALL REPORTS, ONE FAMILY AT A TIME
478 !
479 FAMILY: do jfam = 1, ofl_numFamily
480
481 write(*,*) ' '
482 if (ofl_familyList(jfam) /= 'TO' ) then
483 write(*,*) 'IDENT TYPE DESCRIPTION LAT LONG LEVEL ITEM OBSVD BKGND ANAL POST PROB REPORT'
484 else
485 write(*,*) 'IDENT TYPE DESCRIPTION LAT LONG CHNL ITEM OBSVD BKGND ANAL POST PROB REPORT'
486 end if
487
488 ! loop over all header indices of the family
489 call obs_set_current_header_list( lobsSpaceData, ofl_familyList( jfam ) )
490
491 HEADER: do
492 headerIndex = obs_getHeaderIndex( lobsSpaceData )
493 if (headerIndex < 0) exit HEADER
494
495 llelrej = .false.
496 idburp = obs_headElem_i( lobsSpaceData, OBS_ITY, headerIndex )
497 zlat = obs_headElem_r( lobsSpaceData, OBS_LAT, headerIndex ) * MPC_DEGREES_PER_RADIAN_R8
498 zlon = obs_headElem_r( lobsSpaceData, OBS_LON, headerIndex ) * MPC_DEGREES_PER_RADIAN_R8
499
500 ! loop over all body indices for this headerIndex
501 call obs_set_current_body_list( lobsSpaceData, headerIndex )
502
503 BODY: do
504 bodyIndex = obs_getBodyIndex( lobsSpaceData )
505 if ( bodyIndex < 0) exit BODY
506 ityp = obs_bodyElem_i( lobsSpaceData, OBS_VNM, bodyIndex )
507 if (ityp == BUFR_NETS .or. ityp == BUFR_NEPS .or. &
508 ityp == BUFR_NEPN .or. ityp == BUFR_NESS .or. &
509 ityp == BUFR_NEUS .or. ityp == BUFR_NEVS .or. &
510 ityp == BUFR_NEZD .or. ityp == bufr_sst .or. ityp == BUFR_ICEC .or. &
511 ityp == bufr_vis .or. ityp == bufr_logVis .or. ityp == bufr_gust .or. &
512 ityp == bufr_radarPrecip .or. ityp == bufr_logRadarPrecip ) then
513 llok = (obs_bodyElem_i( lobsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated )
514 else
515 llok = (obs_bodyElem_i( lobsSpaceData, OBS_ASS, bodyIndex ) == obs_assimilated .and. &
516 obs_bodyElem_i( lobsSpaceData, OBS_XTR, bodyIndex ) == 0) .or. &
517 (obs_bodyElem_i( lobsSpaceData, OBS_XTR, bodyIndex ) == 2 .and. &
518 obs_bodyElem_i( lobsSpaceData, OBS_VNM, bodyIndex ) == BUFR_NEGZ)
519 end if
520 if ( llok ) then
521 zpost = obs_bodyElem_r( lobsSpaceData, OBS_QCV, bodyIndex )
522 zlev = obs_bodyElem_r( lobsSpaceData, OBS_PPP, bodyIndex ) * MPC_MBAR_PER_PA_R8
523 if ( obs_bodyElem_i( lobsSpaceData, OBS_VCO, bodyIndex ) == 2) then
524 zlev = obs_bodyElem_r( lobsSpaceData, OBS_PPP, bodyIndex ) * MPC_MBAR_PER_PA_R8
525 CLUNITS = 'MB'
526 else if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex) == 1) then
527 !
528 ! VERTICAL COORDINATE IS HEIGHT
529 !
530 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex)
531 CLUNITS = ' M'
532 if (ofl_familyList(JFAM)=='CH') then
533 ZLEV = ZLEV*0.01
534 CLUNITS = 'HM'
535 end if
536 else if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex) == -1) then
537 !
538 ! TOVS CHANNEL NUMBER
539 !
540 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex)
541 CLUNITS = ' '
542 else
543 ZLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex)
544 CLUNITS = ' '
545 end if
546
547 ZVAR = obs_bodyElem_r(lobsSpaceData,OBS_VAR,bodyIndex)
548
549 ZFCST= ZVAR - obs_bodyElem_r(lobsSpaceData,OBS_OMP,bodyIndex)
550 ZANA = ZVAR - obs_bodyElem_r(lobsSpaceData,OBS_OMA,bodyIndex)
551 !
552 !_____________TREAT WINDS AS SPECIAL CASE
553 ! BUFR_NEUU = 011003 (U COMPONENT) (m/s)
554 ! BUFR_NEVV = 011004 (V COMPONENT) (m/s)
555 ! BUFR_NEUS = 011215 (U COMPONENT AT 10 M) (m/s)
556 ! BUFR_NEVS = 011216 (V COMPONENT AT 10 M) (m/s)
557 !
558 if (((ityp==BUFR_NEUU .or. ityp == BUFR_NEUS) .and. &
559 col_varExist(varName='UU')).or. &
560 ((ityp==BUFR_NEVV .or. ityp == BUFR_NEVS) .and. &
561 col_varExist(varName='VV'))) then
562 IOTHER = -1
563 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEUS) then
564 ISTART=obs_headElem_i(lobsSpaceData,OBS_RLN,headerIndex)
565 do bodyIndex2=ISTART,bodyIndex
566 ISTYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,bodyIndex2)
567 if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex2) == 2) then
568 ZSLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex2)*MPC_MBAR_PER_PA_R8
569 end if
570 if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex2) == 1) then
571 ZSLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex2)
572 end if
573 if ((ISTYP == BUFR_NEVV .or. ISTYP == BUFR_NEVS) .and. &
574 ZLEV == ZSLEV) then
575 IOTHER = bodyIndex2
576 end if
577 end do
578 else
579 ISTART=obs_headElem_i(lobsSpaceData,OBS_RLN,headerIndex)
580 do bodyIndex2=ISTART,bodyIndex
581 ISTYP = obs_bodyElem_i(lobsSpaceData,OBS_VNM,bodyIndex2)
582 if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex2) == 2) then
583 ZSLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex2)*MPC_MBAR_PER_PA_R8
584 end if
585 if (obs_bodyElem_i(lobsSpaceData,OBS_VCO,bodyIndex2) == 1) then
586 ZSLEV = obs_bodyElem_r(lobsSpaceData,OBS_PPP,bodyIndex2)
587 end if
588 if ((ISTYP == BUFR_NEUU .or. ISTYP == BUFR_NEUS) .and. &
589 ZLEV == ZSLEV) then
590 IOTHER = bodyIndex2
591 end if
592 end do
593 end if
594 if (IOTHER /= -1) then
595 if (ityp == BUFR_NEVV .or. ityp == BUFR_NEUS) then
596 ZUU = ZVAR
597 ZVV = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER)
598 else
599 ZVV = ZVAR
600 ZUU = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER)
601 end if
602 SPD = SQRT(ZUU*ZUU + ZVV*ZVV)
603 ISPDO = NINT(SPD)
604 if (ZUU==0.D0 .and. ZVV==0.D0)then
605 IDIRO = 999
606 else
607 DEG = 270. - ATAN2(ZVV,ZUU)*MPC_DEGREES_PER_RADIAN_R8
608 if (DEG > 360.D0)DEG = DEG - 360.D0
609 if (DEG < 0.D0) DEG = DEG + 360.D0
610 IDIRO = NINT(DEG)
611 end if
612 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEUS) then
613 ZUU = ZFCST
614 ZVV = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER) - &
615 obs_bodyElem_r(lobsSpaceData,OBS_OMP,IOTHER)
616 else
617 ZVV = ZFCST
618 ZUU = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER) - &
619 obs_bodyElem_r(lobsSpaceData,OBS_OMP,IOTHER)
620 end if
621 SPD=SQRT(ZUU*ZUU + ZVV*ZVV)
622 ISPDF = NINT(SPD)
623 if (ZUU == 0.D0 .and. ZVV == 0.D0) then
624 IDIRF = 999
625 else
626 DEG = 270.D0 - ATAN2(ZVV,ZUU)*MPC_DEGREES_PER_RADIAN_R8
627 if (DEG > 360.D0)DEG = DEG - 360.D0
628 if (DEG < 0.D0) DEG = DEG + 360.D0
629 IDIRF = NINT(DEG)
630 end if
631 if (ityp == BUFR_NEUU .or. ityp == BUFR_NEUS) then
632 ZUU = ZANA
633 ZVV = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER) - &
634 obs_bodyElem_r(lobsSpaceData,OBS_OMA,IOTHER)
635 else
636 ZVV = ZANA
637 ZUU = obs_bodyElem_r(lobsSpaceData,OBS_VAR,IOTHER) - &
638 obs_bodyElem_r(lobsSpaceData,OBS_OMA,IOTHER)
639 end if
640 SPD=SQRT(ZUU*ZUU + ZVV*ZVV)
641 ISPDA = NINT(SPD)
642 if (ZUU==0.D0 .and. ZVV==0.D0) then
643 IDIRA = 999
644 else
645 DEG = 270.D0 - ATAN2(ZVV,ZUU)*MPC_DEGREES_PER_RADIAN_R8
646 if (DEG > 360.D0)DEG = DEG - 360.D0
647 if (DEG < 0.D0) DEG = DEG + 360.D0
648 IDIRA = NINT(DEG)
649 end if
650 ILEV = NINT(ZLEV)
651 if (ZPOST > ZCUT) then
652 LLELREJ = .TRUE.
653 call obs_bodySet_i(lobsSpaceData,OBS_FLG,bodyIndex, &
654 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,bodyIndex),9))
655 call obs_bodySet_i(lobsSpaceData,OBS_FLG,IOTHER, &
656 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,IOTHER),9))
657 call obs_bodySet_i(lobsSpaceData,OBS_FLG,bodyIndex, &
658 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,bodyIndex),17))
659 call obs_bodySet_i(lobsSpaceData,OBS_FLG,IOTHER, &
660 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,IOTHER),17))
661 if (ityp == BUFR_NEUU .or. &
662 ityp == BUFR_NEVV) then
663 ICOUNT(1,JFAM) = ICOUNT(1,JFAM) + 1
664 end if
665 if (ityp == BUFR_NEUS .or. ityp == BUFR_NEVS) then
666 ICOUNT(10,JFAM) = ICOUNT(10,JFAM) + 1
667 end if
668 codtypname=codtyp_get_name(IDBURP)
669 stnId = obs_elem_c(lobsSpaceData,'STID',headerIndex)
670 obsONM = obs_headElem_i(lobsSpaceData,OBS_ONM,headerIndex)
671 write(*,620) stnId,IDBURP,codtypname,ZLAT,ZLON, &
672 ILEV,CLUNITS,IDIRO,ISPDO,IDIRF,ISPDF,IDIRA, &
673 ISPDA,ZPOST,obsONM
674 620 FORMAT(A9,1X,i3,1x,A21,1X,F5.1,2X,F5.1,2X,I4,A2,2X, &
675 'WND',3X,I3,'/',I3,3X,I3,'/',I3,3X,I3, &
676 '/',I3,2X,F7.4,1X,I8)
677
678 end if
679 end if
680 else
681 ILEV = NINT(ZLEV)
682 if (ZPOST > ZCUT) then
683 LLELREJ = .TRUE.
684 if (ofl_familyList(JFAM)=='CH') then
685 ICOUNT(14,JFAM)=ICOUNT(14,JFAM)+1
686 CLDESC=CLITM(14)
687 if (obs_headElem_i(lobsSpaceData,OBS_CHM,headerIndex).ge.0) then
688
689 ! CONVERT TO PERCENTAGE DIFFERENCE FROM THE FORECAST
690
691 if (abs(ZFCST).lt.1.D-10*max(1.D-30,abs(ZVAR-ZFCST))) then
692 write(*,*) 'vqc_listrej: WARNING for CH obs. BKGND and or OBSVD value out of bounds for output.'
693 write(*,*) 'OBSVD and ANAL set to 0.0 below in addition to BKGND.'
694 ZVAR=0.0D0
695 ZANA=0.0D0
696 else
697 ZVAR=(ZVAR-ZFCST)/ZFCST*100.
698 ZANA=(ZANA-ZFCST)/ZFCST*100.
699 end if
700 ZFCST=0.0
701 CLDESC=vnl_varnameFromVarnum(ityp,obs_headElem_i(lobsSpaceData,OBS_CHM,headerIndex))
702
703 else if (ityp == BUFR_NETT) then
704
705 ! CONVERT FROM KELVIN TO CELCIUS
706
707 CLDESC = CLITM(4)
708 ZVAR = ZVAR - MPC_K_C_DEGREE_OFFSET_R8
709 ZFCST = ZFCST - MPC_K_C_DEGREE_OFFSET_R8
710 ZANA = ZANA - MPC_K_C_DEGREE_OFFSET_R8
711 end if
712 else if (ityp == BUFR_NEGZ) then
713 CLDESC = CLITM(3)
714 ICOUNT(3,JFAM) = ICOUNT(3,JFAM) + 1
715 !
716 ! CONVERT FROM GEOPOTENTIAL TO GEOPOTENTIAL HEIGHT
717 !
718 ZVAR = ZVAR/ec_rg
719 ZFCST = ZFCST/ec_rg
720 ZANA = ZANA/ec_rg
721 else if (ityp == BUFR_NETT) then
722 CLDESC = CLITM(4)
723 ICOUNT(4,JFAM) = ICOUNT(4,JFAM) + 1
724 !
725 ! CONVERT FROM KELVIN TO CELCIUS
726 !
727 ZVAR = ZVAR - MPC_K_C_DEGREE_OFFSET_R8
728 ZFCST = ZFCST - MPC_K_C_DEGREE_OFFSET_R8
729 ZANA = ZANA - MPC_K_C_DEGREE_OFFSET_R8
730 else if (ityp == BUFR_NEES) then
731 CLDESC = CLITM(5)
732 ICOUNT(5,JFAM) = ICOUNT(5,JFAM) + 1
733 else if (ityp == BUFR_NEPS) then
734 CLDESC = CLITM(6)
735 ICOUNT(6,JFAM) = ICOUNT(6,JFAM) + 1
736 !
737 ! CONVERT FROM PASCALS TO MILLIBARS
738 !
739 ZVAR = ZVAR*MPC_MBAR_PER_PA_R8
740 ZANA = ZANA*MPC_MBAR_PER_PA_R8
741 ZFCST = ZFCST*MPC_MBAR_PER_PA_R8
742 else if (ityp == BUFR_NEPN) then
743 CLDESC = CLITM(7)
744 ICOUNT(7,JFAM) = ICOUNT(7,JFAM) + 1
745 !
746 ! CONVERT FROM PASCALS TO MILLIBARS
747 !
748 ZVAR = ZVAR*MPC_MBAR_PER_PA_R8
749 ZANA = ZANA*MPC_MBAR_PER_PA_R8
750 ZFCST = ZFCST*MPC_MBAR_PER_PA_R8
751 else if (ityp == BUFR_NETS) then
752 CLDESC = CLITM(8)
753 ICOUNT(8,JFAM) = ICOUNT(8,JFAM) + 1
754 !
755 ! CONVERT FROM KELVIN TO CELCIUS
756 !
757 ZVAR = ZVAR - MPC_K_C_DEGREE_OFFSET_R8
758 ZFCST = ZFCST - MPC_K_C_DEGREE_OFFSET_R8
759 ZANA = ZANA - MPC_K_C_DEGREE_OFFSET_R8
760 else if (ityp == BUFR_NESS) then
761 CLDESC = CLITM(9)
762 ICOUNT(9,JFAM) = ICOUNT(9,JFAM) + 1
763 else if ( ityp==BUFR_NBT1 .or. ityp==BUFR_NBT2 .or. &
764 ityp==BUFR_NBT3 ) then
765 CLDESC = CLITM(12)
766 ICOUNT(12,JFAM) = ICOUNT(12,JFAM) + 1
767 else if (ityp == BUFR_NEZD) then
768 CLDESC = CLITM(13)
769 ICOUNT(13,JFAM) = ICOUNT(13,JFAM) + 1
770 !
771 ! CONVERT ZTD FROM METRES TO MILLIMETRES
772 !
773 ZVAR = ZVAR * 1000.D0
774 ZFCST = ZFCST * 1000.D0
775 ZANA = ZANA * 1000.D0
776 else if (ityp == bufr_sst) then
777 CLDESC = CLITM(15)
778 ICOUNT(15,JFAM) = ICOUNT(15,JFAM) + 1
779 else if (ityp == BUFR_ICEC) then
780 CLDESC = CLITM(16)
781 ICOUNT(16,JFAM) = ICOUNT(16,JFAM) + 1
782 end if
783
784 call obs_bodySet_i(lobsSpaceData,OBS_FLG,bodyIndex, &
785 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,bodyIndex),9))
786 call obs_bodySet_i(lobsSpaceData,OBS_FLG,bodyIndex, &
787 IBSET(obs_bodyElem_i(lobsSpaceData,OBS_FLG,bodyIndex),17))
788 codtypname=codtyp_get_name(IDBURP)
789 stnId = obs_elem_c(lobsSpaceData,'STID',headerIndex)
790 obsONM = obs_headElem_i(lobsSpaceData,OBS_ONM,headerIndex)
791 write(*,630) stnId,IDBURP, &
792 codtypname,ZLAT,ZLON,ILEV,CLUNITS,CLDESC, &
793 ZVAR,ZFCST,ZANA,ZPOST,obsONM
794 630 FORMAT(A9,1X,I3,1X,A21,1X,F5.1,2X,F5.1,2X,I4,A2,2X, &
795 A4,2X,F7.1,3X,F7.1,3X,F7.1,2X,F7.4,1X,I8,1X)
796 end if
797 end if
798 end if
799 end do BODY
800 !
801 ! NOW SET THE GLOBAL FLAGS IN THE BURP RECORD HEADER
802 !
803 if (LLELREJ) then
804 call obs_headSet_i(lobsSpaceData,OBS_ST1,headerIndex, IBSET( obs_headElem_i(lobsSpaceData,OBS_ST1,headerIndex), 06))
805 end if
806 end do HEADER
807 end do FAMILY
808
809
810 write(*,640)
811 640 FORMAT(//)
812 write(*,*) ' REJECTED DATA ACCORDING TO FAMILY OF REPORT.'
813 write(*,670)(ofl_familyList(JFAM),JFAM=1,ofl_numFamily)
814 670 FORMAT(5X,11(7X,A2))
815 do JITEM=1,NUMITEM
816 if ( .NOT. (JITEM == 2 .or. JITEM == 11)) then
817 write(*,680)CLITM(JITEM),(ICOUNT(JITEM,JFAM),JFAM=1,ofl_numFamily)
818 680 FORMAT(1X,A4,11(4X,I5))
819 end if
820 end do
821 write(*,640)
822
823 end subroutine vqc_listrej
824
825
826end module varQC_mod
