| blob | 32108a02c4f039e3c827a13a5043195b03e7328e |
1 subroutine da_qc_ahi (it, i, nchan, ob, iv)
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for ahi data.
5 ! Method: Assume cloud flag coming from GEOCAT processing
6 ! To be developed: built in cloud_detection method
7 ! HISTORY: 2020/03/01 - Add clear sky cloud detection procedures Dongmei Xu, NUIST, NCAR/MMM
8 !---------------------------------------------------------------------------
10 implicit none
12 integer, intent(in) :: it ! outer loop count
13 integer, intent(in) :: i ! sensor index.
14 integer, intent(in) :: nchan ! number of channel
15 type (y_type), intent(in) :: ob ! Observation structure.
16 type (iv_type), intent(inout) :: iv ! O-B structure.
17 ! local variables
18 logical :: lmix, cloud_detection
19 integer :: n,k,isflg,ios,fgat_rad_unit
20 integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
21 nrej_omb_std(nchan), &
22 nrej_clw(nchan),num_proc_domain, &
23 nrej_mixsurface,nrej_land
25 real :: inv_grosscheck, lowBTcheck
26 ! isflg: SEA(0),ICE(1),LAND(2),SNOW(3),MSEA(4),MICE(5),MLND(6),MSNO(7)
27 integer, parameter :: sea_flag = 0
28 integer, parameter :: ice_flag = 1
29 integer, parameter :: land_flag = 2
30 integer, parameter :: snow_flag = 3
31 integer, parameter :: msea_flag = 4
32 integer, parameter :: mice_flag = 5
33 integer, parameter :: mland_flag = 6
34 integer, parameter :: msnow_flag = 7
35 character(len=30) :: filename
36 logical :: print_cld_debug
38 !! Additional variables used by Zhuge and Zou (2017)
39 integer :: itest
40 logical :: reject_clddet
41 real :: crit_clddet
42 real :: rad_O14, rad_M14, rad_tropt
43 real :: rad_o_ch7, rad_b_ch7, rad_o_ch14, rad_b_ch14
44 real :: Relaz, Glintzen, tb_temp1
45 real :: wave_num(10)
46 real :: plbc1(10), plbc2(10)
47 real :: plfk1(10), plfk2(10)
48 integer, parameter :: num_clddet_tests = 10
49 integer, parameter :: num_clddet_cats = 4
50 real :: eps_clddet(num_clddet_tests+2,num_clddet_cats)
51 integer :: index_clddet(num_clddet_tests), offset_clddet
52 integer :: isflgs_clddet(num_clddet_cats)
53 logical :: qual_clddet(num_clddet_cats)
54 character(len=10) :: crit_names_clddet(num_clddet_tests)
55 integer :: nrej_clddet(nchan,num_clddet_tests+1)
56 integer*2 :: clddet_tests(iv%instid(i)%superob_width, &
57 iv%instid(i)%superob_width, &
58 num_clddet_tests)
59 integer :: superob_center
60 integer :: isuper, jsuper
61 real :: cm, co, ca ! variables for all-sky obs error
63 real, pointer :: tb_ob(:,:), tb_xb(:,:), tb_inv(:,:), tb_xb_clr(:,:), ca_mean(:,:)
64 integer :: tb_qc(nchan), tb_qc_clddet(nchan)
66 real :: big_num
67 real, parameter :: C1=1.19104276e-5 ! = 2 * h * c**2 mWm-2sr-1(cm-1)-4
68 real, parameter :: C2=1.43877516 ! = h * c / b = 1.43877 K(cm-1)-1
70 ! h = Planck's constant
71 ! b = Boltzmann constant
72 ! c = velocity of light
74 integer, parameter :: ch7 = 1
75 integer, parameter :: ch10 = 4
76 integer, parameter :: ch13 = 7
77 integer, parameter :: ch14 = 8
78 integer, parameter :: ch15 = 9
79 ! mmr or pf Cloud Detection Variables
80 integer :: kts_100hPa(1), kte_surf, ndim
81 integer :: numrad_local(nchan), numrad_global(nchan)
82 real :: tstore
83 real :: bias_local(nchan), bias_global(nchan)
84 integer :: kmin, kmax
85 integer, allocatable :: k_cloud_flag(:) ! cloud flags
86 if (trace_use) call da_trace_entry("da_qc_ahi")
87 ! These values can change as SRF (spectral response function) is updated
88 ! It is recommended to acquire these from L1B files, not copy them from GOES R PUG L1b Vol. 3
89 wave_num(1:10) = (/2570.373, 1620.528, 1443.554, 1363.228, 1184.220, &
90 1040.891, 968.001, 894.000, 815.294, 753.790/)
91 plbc1(1:10) = (/0.43361, 1.55228, 0.34427, 0.05651, 0.18733, &
92 0.09102, 0.07550, 0.22516, 0.21702, 0.06266/)
93 plbc2(1:10) = (/0.99939, 0.99667, 0.99918, 0.99986, 0.99948, &
94 0.99971, 0.99975, 0.99920, 0.99916, 0.99974/)
96 plfk1 = C1 * wave_num**3
97 plfk2 = C2 * wave_num
99 crit_names_clddet(1) = "rtct" !here
100 crit_names_clddet(2) = "etrop"
101 crit_names_clddet(3) = "pfmft"
102 crit_names_clddet(4) = "nfmft"
103 crit_names_clddet(5) = "rfmft" !here
104 crit_names_clddet(6) = "cirh2o" !here
105 crit_names_clddet(7) = "emiss4"
106 crit_names_clddet(8) = "ulst" !here
107 crit_names_clddet(9) = "notc"
108 crit_names_clddet(10) = "tempir" !here
110 big_num = huge(big_num)
111 !! Table 4 from Zhuge X. and Zou X. JAMC, 2016. [modified from ABI Cloud Mask Algorithm]
112 !ocean land snow ice (assume same as snow)
113 eps_clddet = transpose( reshape( (/ &
114 3.2, 4.1, big_num, big_num &
115 , 0.1, 0.3, 0.4, 0.4 &
116 , 0.8, 2.5, big_num, big_num &
117 , 1.0, 2.0, 5.0, 5.0 &
118 , 0.7, 1.0, big_num, big_num &
119 , 0.7, 0.7, 0.7, 0.7 &
120 , 0.1, 0.2, 0.3, 0.3 & ! Land values: 0.46 in ABI CM; 0.2 in ZZ16
121 , 2.86, big_num, big_num, big_num &
122 , 0.05, 0.1, 0.12, 0.12 &
123 , 15., 21., 10., 10. &
124 , 11., 15., 4.5, 4.5 &
125 , 2.0, 2.0, 2.0, 2.0 &
126 /), (/ size(eps_clddet, 2), size(eps_clddet, 1) /)) )
127 index_clddet = (/1, 2, 3, 4, 5, 6, 7, 9, 10, 12/)
128 isflgs_clddet = (/sea_flag, land_flag, snow_flag, ice_flag/)
129 print_cld_debug = .false.
130 ngood(:) = 0
131 nrej(:) = 0
132 nrej_omb_abs(:) = 0
133 nrej_omb_std(:) = 0
134 nrej_clw(:) = 0
135 nrej_mixsurface = 0
136 nrej_land = 0
137 nrej_clddet(:,:)= 0
138 num_proc_domain = 0
140 superob_center = ahi_superob_halfwidth + 1
142 tb_xb => iv%instid(i)%tb_xb
143 tb_inv => iv%instid(i)%tb_inv
144 AHIPixelQCLoop: do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
146 tb_ob => ob%instid(i)%tb
148 if (iv%instid(i)%info%proc_domain(1,n)) &
149 num_proc_domain = num_proc_domain + 1
151 ! 0.0 initialise QC by flags assuming good obs
152 !-----------------------------------------------------------------
153 tb_qc = qc_good
155 ! 1.0 reject all channels over mixture surface type
156 !------------------------------------------------------
157 isflg = iv%instid(i)%isflg(n)
158 lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
159 if (lmix) then
160 tb_qc = qc_bad
161 if (iv%instid(i)%info%proc_domain(1,n)) &
162 nrej_mixsurface = nrej_mixsurface + 1
163 end if
165 if ( isflg > 0 ) then
166 do k = 1, nchan
167 if ( k /= 2 .and. k /= 3 .and. k /= 4 ) then
168 if (only_sea_rad) then
169 tb_qc(k) = qc_bad
170 nrej_land = nrej_land + 1
171 end if
172 end if
173 end do
174 end if
175 ! 2.0 check iuse
176 !-----------------------------------------------------------------
177 do k = 1, nchan
178 if (satinfo(i)%iuse(k) .eq. -1) &
179 tb_qc(k) = qc_bad
180 end do
182 ! 3.0 check clw in fg
183 !-----------------------------------------------------------------
184 if (.not. crtm_cloud ) then
186 do k = 1, nchan
187 if (iv%instid(i)%clwp(n) >= 0.2) then
188 tb_qc = qc_bad
189 if (iv%instid(i)%info%proc_domain(1,n)) &
190 nrej_clw(k) = nrej_clw(k) + 1
191 end if
192 end do
193 end if
194 ! METHOD,Zhuge, X. and Zou, X., Test of a Modified Infrared-Only ABI Cloud Mask Algorithm for AHI Radiance Observations. J. Appl. Meteor. Climatol., 2016, 55: 2529–2546.
195 ahi_clddet: if ( use_clddet_zz &
196 .and. all(tb_inv( (/ch7,ch14,ch15/), n ) .gt. missing_r) &
197 .and. all(tb_ob( (/ch7,ch14,ch15/), n ) .gt. missing_r) &
198 .and. all(tb_xb( (/ch7,ch14,ch15/), n ) .gt. missing_r) &
199 ) then
200 !!===============================================================================
201 !!===============================================================================
202 !!
203 !! 4.0 AHI IR-only Cloud Mask Algorithm, combines:
204 !! (*) Heidinger A. and Straka W., ABI Cloud Mask, version 3.0, 11 JUN, 2013.
205 !! (*) Zhuge X. and Zou X. JAMC, 2016.
206 !!
207 !!===============================================================================
208 !!===============================================================================
210 !JJGDEBUG
211 if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG1: ', n, &
212 tb_inv(:,n)
213 if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG2: ', n, &
214 tb_xb(:,n)
215 if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG3: ', n, &
216 tb_ob(:,n)
217 if (crtm_cloud ) then
218 if (print_cld_debug) write(stdout,'(A,I8,*(2x,F10.4:))') 'PIXEL_DEBUG4: ', n, &
219 iv%instid(i)%tb_xb_clr(:,n)
220 end if
222 if (print_cld_debug) write(stdout,'(A,I8,8F12.4,2x,A)') 'PIXEL_DEBUG5: ', n, &
223 iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n), &
224 iv%instid(i)%satzen(n), iv%instid(i)%satazi(n), &
225 iv%instid(i)%solzen(n), iv%instid(i)%solazi(n), &
226 iv%instid(i)%tropt(n), iv%instid(i)%superob(superob_center,superob_center)%cld_qc(n)%terr_hgt, &
227 iv%instid(i)%info%date_char(n)
228 !JJGDEBUG
230 clddet_tests = 0
232 do jsuper = 1, iv%instid(i)%superob_width
233 do isuper = 1, iv%instid(i)%superob_width
235 tb_ob => iv%instid(i)%superob(isuper,jsuper)%tb_obs
237 if ( tb_xb(ch14,n) .gt. 0. .and. iv%instid(i)%tropt(n) .gt. 0. ) then
238 tb_temp1 = tb_ob(ch14,n)
239 rad_O14 = plfk1(ch14) / &
240 ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1 ) ) -1 )
241 tb_temp1 = tb_xb(ch14,n)
242 rad_M14 = plfk1(ch14) / &
243 ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1) ) -1 )
244 tb_temp1 = iv%instid(i)%tropt(n)
245 rad_tropt = plfk1(ch14) / &
246 ( exp( plfk2(ch14) / (plbc1(ch14) + plbc2(ch14) * tb_temp1) ) -1 )
247 else
248 rad_O14 = missing_r
249 rad_M14 = missing_r
250 rad_tropt = missing_r
251 end if
253 if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) then
254 tb_temp1 = tb_ob(ch7,n)
255 rad_o_ch7 = plfk1(ch7) / &
256 ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. )
257 tb_temp1 = tb_xb(ch7,n)
258 rad_b_ch7 = plfk1(ch7) / &
259 ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. )
260 tb_temp1 = tb_ob(ch14,n)
261 rad_o_ch14 = plfk1(ch7) / &
262 ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. )
263 tb_temp1 = tb_xb(ch14,n)
264 rad_b_ch14 = plfk1(ch7) / &
265 ( exp( plfk2(ch7) / ( plbc1(ch7) + plbc2(ch7) * tb_temp1 ) ) - 1. )
266 else
267 rad_o_ch7 = missing_r
268 rad_b_ch7 = missing_r
269 rad_o_ch14 = missing_r
270 rad_b_ch14 = missing_r
271 end if
273 tb_qc_clddet = tb_qc
275 AHICloudTestLoop: do itest = 1, num_clddet_tests
276 qual_clddet = .true.
277 offset_clddet = 0
278 crit_clddet = missing_r
279 select case (itest)
280 case (1)
281 !--------------------------------------------------------------------------
282 ! 4.1 Relative Thermal Contrast Test (RTCT)
283 !--------------------------------------------------------------------------
284 crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%RTCT
285 qual_clddet(3:4) = .false.
286 case (2)
287 !--------------------------------------------------------------------------
288 ! 4.2 Cloud check: step 1
289 ! Emissivity at Tropopause Test (ETROP)
290 !--------------------------------------------------------------------------
291 if ( tb_xb(ch14,n) .gt. 0. .and. iv%instid(i)%tropt(n) .gt. 0. ) &
292 crit_clddet = (rad_O14 - rad_M14) / (rad_tropt - rad_M14)
293 case (3)
294 !--------------------------------------------------------------------------
295 ! 4.3 Cloud check: step 2
296 ! Positive Fourteen Minus Fifteen Test (PFMFT)
297 !--------------------------------------------------------------------------
298 ! See ABI Cloud Mask Description for qual_clddet
299 qual_clddet = (tb_xb(ch14,n).ge.tb_xb(ch15,n))
300 if ( (tb_inv(ch14,n) + tb_xb(ch14,n)).le.310. .and. &
301 iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_14.ge.0.3 .and. &
302 tb_ob(ch14,n).gt.0. .and. tb_ob(ch15,n).gt.0. ) &
303 crit_clddet = ( tb_ob(ch14,n) - tb_ob(ch15,n) )
304 ! above using ob without VarBC
305 ! -------------------------------
306 ! crit_clddet = (tb_inv(ch14,n) + tb_xb(ch14,n) - &
307 ! (tb_inv(ch15,n) + tb_xb(ch15,n)) )-&
308 ! (tb_xb(ch14,n) - tb_xb(ch15,n)) * &
309 ! (tb_ob(ch14,n) - 260.) / (tb_xb(ch14,n) - 260.)
310 ! above using ob with VarBC
311 ! -------------------------------
312 if ( crit_clddet.gt.missing_r .and. &
313 (tb_inv(ch14,n) + tb_xb(ch14,n)).gt.270. .and. &
314 tb_xb(ch14,n).gt.270. ) &
315 crit_clddet = crit_clddet - &
316 (tb_xb(ch14,n) - tb_xb(ch15,n)) * &
317 (tb_ob(ch14,n) - 260.) / (tb_xb(ch14,n) - 260.)
318 ! above 1 line using ob without VarBC
319 ! (tb_inv(ch14,n) + tb_xb(ch14,n) - 260.)/ &
320 ! (tb_xb(ch14,n) - 260.)
321 ! above 2 lines using ob with VarBC
323 case (4)
324 !--------------------------------------------------------------------------
325 ! 4.4 Negative Fourteen Minus Fifteen Test (NFMFT)
326 !--------------------------------------------------------------------------
327 if (tb_ob(ch14,n) .gt. 0. .and. tb_ob(ch15,n) .gt. 0.) &
328 crit_clddet = tb_inv(ch15,n) - tb_inv(ch14,n)
329 case (5)
330 !--------------------------------------------------------------------------
331 ! 4.5 Relative Fourteen Minus Fifteen Test (RFMFT)
332 !--------------------------------------------------------------------------
333 ! See ABI Cloud Mask Description for qual_clddet
334 qual_clddet = ( tb_ob(ch14,n) - tb_ob(ch15,n) ) .lt. 1.0
335 qual_clddet(2) = qual_clddet(2) .and. tb_ob(ch14,n) .le. 300.
336 qual_clddet(3:4) = .false.
338 crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%RFMFT
339 case (6)
340 !--------------------------------------------------------------------------
341 ! 4.6 Cirrus Water Vapor Test (CIRH2O)
342 !--------------------------------------------------------------------------
343 ! See ABI Cloud Mask Description for qual_clddet
344 qual_clddet = &
345 iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%terr_hgt .le. 2000. &
346 .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%terr_hgt .ge. 0. &
347 .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_10 .gt. 0.5 &
348 .and. iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_14 .gt. 0.5
349 crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%CIRH2O
350 case (7)
351 !--------------------------------------------------------------------------
352 ! 4.7 Modified 4um Emissivity Test (M-4EMISS)
353 !--------------------------------------------------------------------------
354 ! Modify EMISS for sun glint area may be not work, because we are at north land
355 ! - compute relative azimuth
356 Relaz = RELATIVE_AZIMUTH(iv%instid(i)%solazi(n),iv%instid(i)%satazi(n))
358 ! - compute glint angle
359 Glintzen = GLINT_ANGLE(iv%instid(i)%solzen(n),iv%instid(i)%satzen(n),Relaz )
361 if ( Glintzen.lt.40.0 .and. isflg==sea_flag) then
362 crit_clddet = - tb_inv(ch7,n) ! (B_ch7 - O_ch7)
363 offset_clddet = 1
364 else
365 if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) &
366 crit_clddet = (rad_o_ch7/rad_o_ch14 - rad_b_ch7/rad_b_ch14)/ &
367 (rad_b_ch7 / rad_b_ch14)
368 end if
369 case (8)
370 !--------------------------------------------------------------------------
371 ! 4.8 Uniform low stratus Test (ULST)
372 !--------------------------------------------------------------------------
373 !JJG, AHI error: Changed this to solzen instead of solazi for night/day test
374 qual_clddet = iv%instid(i)%solzen(n) >= 85.0
375 if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch14,n) .gt. 0.) &
376 crit_clddet = rad_b_ch7/rad_b_ch14 - rad_o_ch7/rad_o_ch14
377 case (9)
378 !--------------------------------------------------------------------------
379 ! 4.9 New Optically Thin Cloud Test (N-OTC)
380 !--------------------------------------------------------------------------
381 !JJG, AHI error: Changed this to solzen instead of solazi for night/day test
382 if ( iv%instid(i)%solzen(n) .ge. 85.0 ) &
383 offset_clddet = 1 ! night time
385 if (tb_ob(ch7,n) .gt. 0. .and. tb_ob(ch15,n) .gt. 0.) &
386 ! using ob without VarBC
387 ! -------------------------------
388 crit_clddet = tb_ob(ch7,n) - tb_ob(ch15,n)
390 ! using ob with VarBC
391 ! -------------------------------
392 ! crit_clddet = tb_inv(ch7,n) + tb_xb(ch7,n) - &
393 ! (tb_inv(ch15,n) + tb_xb(ch15,n))
394 case (10)
395 !--------------------------------------------------------------------------
396 ! 4.10 Temporal Infrared Test (TEMPIR)
397 !--------------------------------------------------------------------------
398 crit_clddet = iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%TEMPIR
399 case default
400 cycle
401 end select
403 ! call evaluate_clddet_test ( &
404 ! isflg, isflgs_clddet, crit_clddet, eps_clddet(index_clddet(itest)+offset_clddet,:), qual_clddet, &
405 ! iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n), &
406 ! reject_clddet )
408 reject_clddet = crit_clddet.gt.missing_r .and. &
409 any( isflg.eq.isflgs_clddet .and. &
410 crit_clddet.gt.eps_clddet(index_clddet(itest)+offset_clddet,:) .and. &
411 qual_clddet )
413 if (reject_clddet) then
414 tb_qc_clddet = qc_bad ! CSS do we want to set it bad for a given pixel within a superob?
415 if (iv%instid(i)%info%proc_domain(1,n)) then
416 nrej_clddet(:,itest) = nrej_clddet(:,itest) + 1
417 ! write(stdout,"(A,F14.6,A,I4,2D12.4)") trim(crit_names_clddet(itest)), crit_clddet, " isflg", isflg, iv%instid(i)%info%lat(1,n), iv%instid(i)%info%lon(1,n)
418 end if
420 clddet_tests(isuper, jsuper, itest) = 1
421 end if
422 end do AHICloudTestLoop
423 end do ! isuper
424 end do ! jsuper
426 iv%instid(i)%cloud_frac(n) = &
427 real( count(sum(clddet_tests,3) > 0),8) / real( iv%instid(i)%superob_width**2,8)
429 if (.not. crtm_cloud ) tb_qc = tb_qc_clddet ! CSS logic here isn't quite right if superobbing
431 !JJGDEBUG
432 if (print_cld_debug) write(stdout,'(A,I8,*(2x,I1:))') 'PIXEL_DEBUG6: ', n, clddet_tests
433 !JJGDEBUG
434 else ! not clddet_zz
435 ! 4. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016)
436 !---------------------------------------------
437 if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
438 iv%instid(i)%cloud_flag(:,n) = qc_good
440 if (rtm_option == rtm_option_rttov) then
441 #ifdef RTTOV
442 kte_surf = iv%instid(i)%nlevels
443 kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), &
444 MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0)
445 do k=1,nchan
446 tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * &
447 (ob%instid(i)%tb(k,n) - bias_global(k))
448 iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / &
449 (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0)
450 end do
451 #endif
452 elseif (rtm_option == rtm_option_crtm) then
453 kte_surf = kte
454 kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), &
455 MASK = iv%instid(i)%pm(kts:kte,n) < 100.0)
457 do k = 1, nchan
458 CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), &
459 iv%instid(i)%rad_obs(k,n))
460 end do
462 end if
464 ndim = kte_surf - kts_100hPa(1) + 1
466 call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv)
467 end if
469 do k = 1, nchan
470 if (iv%instid(i)%cloud_flag(k,n) == qc_bad) tb_qc(k) = qc_bad
471 end do
472 end if ahi_clddet
474 tb_ob => ob%instid(i)%tb
476 ! ---------------------------calculate and save ca_mean for crtm_cloud and crtm_clr
477 ! 5.0 assigning obs errors
478 ! tb_xb_clr => iv%instid(i)%tb_xb_clr ! currently not used
479 if (.not. crtm_cloud ) then
480 do k = 1, nchan
481 if (use_error_factor_rad) then
482 iv%instid(i)%tb_error(k,n) = &
483 satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
484 else
485 iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
486 end if
487 end do ! nchan
489 else ! Added this else block...until now obs error = 500 if crtm_cloud = T...not good!
491 if ( ahi_use_symm_obs_err ) then
492 do k = 1, nchan
493 ! Okamato et al. (2014)
494 !cm = iv%instid(i)%tb_xb(k,n) - iv%instid(i)%tb_xb_clr(k,n)
495 !co = ob%instid(i)%tb(k,n) - iv%instid(i)%tb_xb_clr(k,n)
497 ! Harnisch et al. (2016)
498 cm = max(0.0, satinfo(i)%BTLim(k) - iv%instid(i)%tb_xb(k,n))
499 co = max(0.0, satinfo(i)%BTLim(k) - ob%instid(i)%tb(k,n))
501 ! Symmetric cloud amount
502 ca = 0.5*( abs(cm) + abs(co) )
504 ! Figure out observation error as a function of ca
505 if (ca.lt.satinfo(i)%ca1(k)) then
506 iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)
507 else if (ca.ge.satinfo(i)%ca1(k) .and. ca.lt.satinfo(i)%ca2(k)) then
508 iv%instid(i)%tb_error(k,n)= satinfo(i)%error_std(k)+ &
509 (ca-satinfo(i)%ca1(k))*(satinfo(i)%error_cld(k)-satinfo(i)%error_std(k))/(satinfo(i)%ca2(k)-satinfo(i)%ca1(k))
510 else
511 iv%instid(i)%tb_error(k,n)= satinfo(i)%error_cld(k)
512 end if
513 end do ! nchan
514 else
515 iv%instid(i)%tb_error(:,n)= 500.0 ! this is the default
516 endif
517 end if
519 ! 6.0 gross 8k check -clr,sdobs-clr
520 !-----------------------------------------------------------------
521 if (.not. crtm_cloud ) then
522 ! absolute departure check
523 do k = 1, nchan
524 inv_grosscheck = 8.0
525 if (use_satcv(2)) inv_grosscheck = 100.0
526 if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then
527 tb_qc(k) = qc_bad
528 if (iv%instid(i)%info%proc_domain(1,n)) &
529 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
530 end if
531 end do ! nchan
532 if (use_clddet_zz) then
533 ! SDob cloud inhomogeneous check
534 do isuper = 1, iv%instid(i)%superob_width
535 do jsuper = 1, iv%instid(i)%superob_width
536 if (iv%instid(i)%superob(isuper,jsuper)%cld_qc(n)%tb_stddev_13 >= 2) then ! only use abs clear pixel
537 tb_qc = qc_bad
538 if (iv%instid(i)%info%proc_domain(1,n)) &
539 nrej_clddet(:,11)= nrej_clddet(:,11)+1
540 end if
541 end do
542 end do
543 end if
544 end if
546 ! 7.0 3std check
547 !-----------------------------------------------------------------
548 do k = 1, nchan
549 ! relative departure check
550 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
551 tb_qc(k) = qc_bad
552 if (iv%instid(i)%info%proc_domain(1,n)) &
553 nrej_omb_std(k) = nrej_omb_std(k) + 1
555 end if
556 end do ! nchan
558 !final QC decsion
559 ! CSS comment out below. Dangerous, especially if satinfo(i)%iuse(k) = -1
560 ! CSS also okay to fail 3std check if using symmetric error model for obs errors
561 !if (crtm_cloud ) tb_qc(2:4) = qc_good ! no qc for crtm_cloud
563 iv%instid(i)%tb_qc(:,n) = tb_qc
564 do k = 1, nchan
565 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
566 iv%instid(i)%tb_error(k,n) = 500.0
567 if (iv%instid(i)%info%proc_domain(1,n)) &
568 nrej(k) = nrej(k) + 1
569 else
570 if (iv%instid(i)%info%proc_domain(1,n)) &
571 ngood(k) = ngood(k) + 1
572 end if
573 end do ! nchan
574 end do AHIPixelQCLoop
576 ! Do inter-processor communication to gather statistics.
577 call da_proc_sum_int (num_proc_domain)
578 call da_proc_sum_int (nrej_mixsurface)
579 call da_proc_sum_int (nrej_land)
580 call da_proc_sum_ints (nrej_omb_abs)
581 call da_proc_sum_ints (nrej_omb_std)
582 call da_proc_sum_ints (nrej_clw)
583 do itest = 1, num_clddet_tests+1
584 call da_proc_sum_ints (nrej_clddet(:,itest))
585 end do
586 call da_proc_sum_ints (nrej)
587 call da_proc_sum_ints (ngood)
589 if (rootproc) then
590 if (num_fgat_time > 1) then
591 write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
592 else
593 write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
594 end if
596 call da_get_unit(fgat_rad_unit)
597 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
598 if (ios /= 0) then
599 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
600 call da_error(__FILE__,__LINE__,message(1:1))
601 end if
603 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
604 if(num_proc_domain > 0) write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
605 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
606 write(fgat_rad_unit,'(a20,i7)') ' nrej_land = ', nrej_land
607 write(fgat_rad_unit,'(a20)') ' nrej_clw(:) = '
608 write(fgat_rad_unit,'(10i7)') nrej_clw(:)
609 do itest = 1, num_clddet_tests+1
610 write(fgat_rad_unit,'(a20,i2,a2)') ' nrej_clddet',itest,"="
611 write(fgat_rad_unit,'(10i7)') nrej_clddet(:,itest)
612 end do
613 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
614 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
615 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
616 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
617 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
618 write(fgat_rad_unit,'(10i7)') nrej(:)
619 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
620 write(fgat_rad_unit,'(10i7)') ngood(:)
622 close(fgat_rad_unit)
623 call da_free_unit(fgat_rad_unit)
624 end if
625 if (trace_use) call da_trace_exit("da_qc_ahi")
627 end subroutine da_qc_ahi
629 function relative_azimuth ( sol_az ,sen_az )
631 implicit none
633 real :: sol_az
634 real :: sen_az
635 real :: relative_azimuth
637 relative_azimuth = abs(sol_az - sen_az)
638 if (relative_azimuth > 180.0) then
639 relative_azimuth = 360.0 - relative_azimuth
640 endif
641 relative_azimuth = 180.0 - relative_azimuth
643 end function relative_azimuth
646 function glint_angle ( sol_zen , sat_zen , rel_az )
647 !------------------------------------------------------------------------------------
648 ! Glint angle (the angle difference between direct "specular" reflection off
649 ! the surface and actual reflection toward the satellite.)
650 !------------------------------------------------------------------------------------
652 implicit none
654 real :: sol_zen
655 real :: sat_zen
656 real :: rel_az
657 real :: glint_angle
659 glint_angle = cos(sol_zen * deg2rad) * cos(sat_zen * deg2rad) + &
660 sin(sol_zen * deg2rad) * sin(sat_zen * deg2rad) * cos(rel_az * deg2rad)
661 glint_angle = max(-1.0 , min( glint_angle ,1.0 ))
662 glint_angle = acos(glint_angle) / deg2rad
664 end function glint_angle