1 subroutine da_qc_amsua (it, i, nchan, ob, iv)
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for amsua data.
5 !---------------------------------------------------------------------------
9 integer, intent(in) :: it ! outer loop count
10 integer, intent(in) :: i ! sensor index.
11 integer, intent(in) :: nchan ! number of channel
12 type (y_type), intent(in) :: ob ! Observation structure.
13 type (iv_type), intent(inout) :: iv ! O-B structure.
17 integer :: n,scanpos,k,isflg,ios,fgat_rad_unit
21 integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
22 nrej_omb_std(nchan), &
23 nrej_mixsurface,nrej_windowchanl, nrej_si, &
24 nrej_clw,nrej_topo, num_proc_domain, &
27 character(len=30) :: filename
29 if (trace_use) call da_trace_entry("da_qc_amsua")
43 do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
45 if (iv%instid(i)%info%proc_domain(1,n)) &
46 num_proc_domain = num_proc_domain + 1
48 ! 0.0 initialise QC by flags assuming good obs
49 !---------------------------------------------
50 iv%instid(i)%tb_qc(:,n) = qc_good
52 ! a. reject all channels over mixture surface type
53 !------------------------------------------------------
54 isflg = iv%instid(i)%isflg(n)
55 lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
57 iv%instid(i)%tb_qc(:,n) = qc_bad
58 if (iv%instid(i)%info%proc_domain(1,n)) &
59 nrej_mixsurface = nrej_mixsurface + 1
61 ! b. reject channels 1~4 over land/sea-ice/snow
62 !------------------------------------------------------
64 iv%instid(i)%tb_qc(1:4,n) = qc_bad
65 if (iv%instid(i)%info%proc_domain(1,n)) &
66 nrej_windowchanl = nrej_windowchanl + 1
67 ! reject whole pixel if not over sea for global case
68 if (global) iv%instid(i)%tb_qc(:,n) = qc_bad
69 if (only_sea_rad) iv%instid(i)%tb_qc(:,n) = qc_bad
72 ! c. reject channels 13,14(above top model 10mb),15
73 !------------------------------------------------------
74 iv%instid(i)%tb_qc(13:15,n) = qc_bad
77 !------------------------------------------------------
78 scanpos = iv%instid(i)%scanpos(n)
79 if (scanpos <= 3 .or. scanpos >= 28) then
80 iv%instid(i)%tb_qc(:,n) = qc_bad
81 if (iv%instid(i)%info%proc_domain(1,n)) &
82 nrej_limb = nrej_limb + 1
86 ! if (abs(satzen) > 45.0) iv%instid(i)%tb_qc(:,n) = qc_bad
88 ! d. check precipitation
89 !-----------------------------------------------------------
90 if (ob%instid(i)%tb(1,n) > 0.0 .and. &
91 ob%instid(i)%tb(15,n) > 0.0) then
92 si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(15,n)
94 iv%instid(i)%tb_qc(:,n) = qc_bad
95 iv%instid(i)%cloud_flag(:,n) = qc_bad
96 if (iv%instid(i)%info%proc_domain(1,n)) &
101 if (iv%instid(i)%clwp(n) >= 0.2) then
102 iv%instid(i)%tb_qc(:,n) = qc_bad
103 iv%instid(i)%cloud_flag(:,n) = qc_bad
104 if (iv%instid(i)%info%proc_domain(1,n)) &
105 nrej_clw = nrej_clw + 1
108 ! 3.1 Estimate Cloud Liquid Water (CLW) in mm over sea
109 ! (Grody etal. 2001, JGR, Equation 5b,7c,7d,9)
110 !---------------------------------------------------------
111 ! if (isflg == 0) then
112 ! coszen = cos(iv%instid(i)%satzen(n))
113 ! d0 = 8.24-(2.622-1.846*coszen)*coszen
116 ! ts = iv%instid(i)%ts(n)
117 ! tb1 = ob%instid(i)%tb(1,n)
118 ! tb2 = ob%instid(i)%tb(2,n)
119 ! clw = coszen*(d0+d1*log(ts-tb1)+d2*log(ts-tb2))
124 ! e. check surface height/pressure
125 !-----------------------------------------------------------
126 ! sfchgt = ivrad%info(n)%elv
127 ! if (sfchgt >=) then
131 if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.0)) then
132 iv%instid(i)%tb_qc(5,n) = qc_bad
133 if (iv%instid(i)%info%proc_domain(1,n)) &
134 nrej_topo = nrej_topo + 1
138 !-----------------------------------------------------------
140 if (satinfo(i)%iuse(k) .eq. -1) &
141 iv%instid(i)%tb_qc(k,n) = qc_bad
144 ! f. check innovation
145 !-----------------------------------------------------------
148 ! absolute departure check
149 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
150 iv%instid(i)%tb_qc(k,n) = qc_bad
151 if (iv%instid(i)%info%proc_domain(1,n)) &
152 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
155 ! relative departure check
156 if (use_error_factor_rad) then
157 iv%instid(i)%tb_error(k,n) = &
158 satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
160 iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
163 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
164 iv%instid(i)%tb_qc(k,n) = qc_bad
165 if (iv%instid(i)%info%proc_domain(1,n)) &
166 nrej_omb_std(k) = nrej_omb_std(k) + 1
170 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
171 iv%instid(i)%tb_error(k,n) = 500.0
172 if (iv%instid(i)%info%proc_domain(1,n)) &
173 nrej(k) = nrej(k) + 1
175 if (iv%instid(i)%info%proc_domain(1,n)) &
176 ngood(k) = ngood(k) + 1
180 end do ! end loop pixel
182 ! Do inter-processor communication to gather statistics.
183 call da_proc_sum_int (num_proc_domain)
184 call da_proc_sum_int (nrej_mixsurface)
185 call da_proc_sum_int (nrej_windowchanl)
186 call da_proc_sum_int (nrej_si )
187 call da_proc_sum_int (nrej_clw)
188 call da_proc_sum_int (nrej_topo)
189 call da_proc_sum_int (nrej_limb)
190 call da_proc_sum_ints (nrej_omb_abs(:))
191 call da_proc_sum_ints (nrej_omb_std(:))
192 call da_proc_sum_ints (nrej(:))
193 call da_proc_sum_ints (ngood(:))
196 if (num_fgat_time > 1) then
197 write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
199 write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
202 call da_get_unit(fgat_rad_unit)
203 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
205 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
206 call da_error(__FILE__,__LINE__,message(1:1))
209 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
210 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
211 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
212 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
213 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
214 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
215 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
216 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
217 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
218 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
219 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
220 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
221 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
222 write(fgat_rad_unit,'(10i7)') nrej(:)
223 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
224 write(fgat_rad_unit,'(10i7)') ngood(:)
227 call da_free_unit(fgat_rad_unit)
230 if (trace_use) call da_trace_exit("da_qc_amsua")
232 end subroutine da_qc_amsua