da_qc_amsua.inc
References to this file elsewhere.
1 subroutine da_qc_amsua (i, nchan, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for amsua data.
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 integer, intent(in) :: i ! sensor index.
10 integer, intent(in) :: nchan ! number of channel
11 type (y_type), intent(in) :: ob ! Observation structure.
12 type (ob_type), intent(inout) :: iv ! O-B structure.
13
14
15 ! local variables
16 integer :: n,scanpos,k,isflg,ios,fgat_rad_unit
17 logical :: lmix
18 real :: si
19 ! real :: satzen
20 integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
21 nrej_omb_std(nchan), &
22 nrej_mixsurface,nrej_windowchanl, nrej_si, &
23 nrej_clw,nrej_topo,num_rad, num_proc_domain, &
24 nrej_limb
25
26 character(len=30) :: filename
27
28 if (trace_use) call da_trace_entry("da_qc_amsua")
29
30 ngood(:) = 0
31 nrej(:) = 0
32 nrej_omb_abs(:) = 0
33 nrej_omb_std(:) = 0
34 nrej_mixsurface = 0
35 nrej_windowchanl= 0
36 nrej_si = 0
37 nrej_clw = 0
38 nrej_topo = 0
39 nrej_limb = 0
40 num_rad = iv%ob_numb(iv%current_ob_time)%radiance(i)- &
41 iv%ob_numb(iv%current_ob_time-1)%radiance(i)
42 ! num_rad = iv%instid(i)%num_rad
43 num_proc_domain = 0
44
45 if (num_rad > 0) then
46
47 ! do n= 1, iv%instid(i)%num_rad ! loop for pixel
48 do n= iv%ob_numb(iv%current_ob_time-1)%radiance(i)+1, iv%ob_numb(iv%current_ob_time)%radiance(i)
49
50 if (iv%instid(i)%proc_domain(n)) &
51 num_proc_domain = num_proc_domain + 1
52
53 ! 0.0 initialise QC by flags assuming good obs
54 !---------------------------------------------
55 iv%instid(i)%tb_qc(:,n) = qc_good
56
57 ! a. reject all channels over mixture surface type
58 !------------------------------------------------------
59 isflg = iv%instid(i)%isflg(n)
60 lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
61 if (lmix) then
62 iv%instid(i)%tb_qc(:,n) = qc_bad
63 if (iv%instid(i)%proc_domain(n)) &
64 nrej_mixsurface = nrej_mixsurface + 1
65 end if
66 ! b. reject channels 1~4 over land/sea-ice/snow
67 !------------------------------------------------------
68 if (isflg > 0) then
69 iv%instid(i)%tb_qc(1:4,n) = qc_bad
70 if (iv%instid(i)%proc_domain(n)) &
71 nrej_windowchanl = nrej_windowchanl + 1
72 ! reject whole pixel if not over sea for global case
73 if (global) iv%instid(i)%tb_qc(:,n) = qc_bad
74 if (only_sea_rad) iv%instid(i)%tb_qc(:,n) = qc_bad
75 end if
76
77 ! c. reject channels 13,14(above top model 10mb),15
78 !------------------------------------------------------
79 iv%instid(i)%tb_qc(13:15,n) = qc_bad
80
81 ! reject limb obs
82 !------------------------------------------------------
83 scanpos = iv%instid(i)%scanpos(n)
84 if (scanpos <= 3 .or. scanpos >= 28) then
85 iv%instid(i)%tb_qc(:,n) = qc_bad
86 if (iv%instid(i)%proc_domain(n)) &
87 nrej_limb = nrej_limb + 1
88 end if
89
90 ! satzen = rad%satzen
91 ! if (abs(satzen) > 45.) iv%instid(i)%tb_qc(:,n) = qc_bad
92
93 ! d. check precipitation
94 !-----------------------------------------------------------
95 if (ob%instid(i)%tb(1,n) > 0. .and. &
96 ob%instid(i)%tb(15,n) > 0.) then
97 si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(15,n)
98 if (si >= 3.) then
99 iv%instid(i)%tb_qc(:,n) = qc_bad
100 iv%instid(i)%cloud_flag(:,n) = qc_bad
101 if (iv%instid(i)%proc_domain(n)) &
102 nrej_si = nrej_si + 1
103 end if
104 end if
105
106 if (iv%instid(i)%clwp(n) >= 0.2) then
107 iv%instid(i)%tb_qc(:,n) = qc_bad
108 iv%instid(i)%cloud_flag(:,n) = qc_bad
109 if (iv%instid(i)%proc_domain(n)) &
110 nrej_clw = nrej_clw + 1
111 end if
112
113 ! 3.1 Estimate Cloud Liquid Water (CLW) in mm over sea
114 ! (Grody etal. 2001, JGR, Equation 5b,7c,7d,9)
115 !---------------------------------------------------------
116 ! if (isflg == 0) then
117 ! coszen = cos(iv%instid(i)%satzen(n))
118 ! d0 = 8.24-(2.622-1.846*coszen)*coszen
119 ! d1 = 0.754
120 ! d2 = -2.265
121 ! ts = iv%instid(i)%ts(n)
122 ! tb1 = ob%instid(i)%tb(1,n)
123 ! tb2 = ob%instid(i)%tb(2,n)
124 ! clw = coszen*(d0+d1*log(ts-tb1)+d2*log(ts-tb2))
125 ! clw = clw - 0.03
126 ! end if
127
128
129 ! e. check surface height/pressure
130 !-----------------------------------------------------------
131 ! sfchgt = ivrad%info(n)%elv
132 ! if (sfchgt >=) then
133 ! else
134 ! end if
135
136 if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.)) then
137 iv%instid(i)%tb_qc(5,n) = qc_bad
138 if (iv%instid(i)%proc_domain(n)) &
139 nrej_topo = nrej_topo + 1
140 end if
141
142 ! g. check iuse
143 !-----------------------------------------------------------
144 do k = 1, nchan
145 if (satinfo(i)%iuse(k) .eq. -1) &
146 iv%instid(i)%tb_qc(k,n) = qc_bad
147 end do
148
149 ! f. check innovation
150 !-----------------------------------------------------------
151 do k = 1, nchan
152
153 ! absolute departure check
154 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.) then
155 iv%instid(i)%tb_qc(k,n) = qc_bad
156 if (iv%instid(i)%proc_domain(n)) &
157 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
158 end if
159
160 ! relative departure check
161 if (use_error_factor_rad) then
162 iv%instid(i)%tb_error(k,n) = &
163 satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
164 else
165 iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
166 end if
167
168 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
169 iv%instid(i)%tb_qc(k,n) = qc_bad
170 if (iv%instid(i)%proc_domain(n)) &
171 nrej_omb_std(k) = nrej_omb_std(k) + 1
172 end if
173
174 ! final QC decsion
175 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
176 iv%instid(i)%tb_error(k,n) = 500.0
177 if (iv%instid(i)%proc_domain(n)) &
178 nrej(k) = nrej(k) + 1
179 else
180 if (iv%instid(i)%proc_domain(n)) &
181 ngood(k) = ngood(k) + 1
182 end if
183 end do ! chan
184 end do ! end loop pixel
185 end if
186
187 ! Do inter-processor communication to gather statistics.
188 call da_proc_sum_int (num_proc_domain)
189 call da_proc_sum_int (nrej_mixsurface)
190 call da_proc_sum_int (nrej_windowchanl)
191 call da_proc_sum_int (nrej_si )
192 call da_proc_sum_int (nrej_clw)
193 call da_proc_sum_int (nrej_topo)
194 call da_proc_sum_int (nrej_limb)
195 call da_proc_sum_ints (nrej_omb_abs(:))
196 call da_proc_sum_ints (nrej_omb_std(:))
197 call da_proc_sum_ints (nrej(:))
198 call da_proc_sum_ints (ngood(:))
199
200 if (rootproc) then
201 if (num_fgat_time > 1) then
202 write(filename,'(a,i2.2)') 'qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%current_ob_time
203 else
204 filename = 'qcstat_'//trim(iv%instid(i)%rttovid_string)
205 end if
206
207 call da_get_unit(fgat_rad_unit)
208 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
209 if (ios /= 0) Then
210 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
211 call da_error(__FILE__,__LINE__,message(1:1))
212 end if
213
214 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
215 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
216 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
217 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
218 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
219 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
220 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
221 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
222 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
223 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
224 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
225 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
226 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
227 write(fgat_rad_unit,'(10i7)') nrej(:)
228 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
229 write(fgat_rad_unit,'(10i7)') ngood(:)
230
231 close(fgat_rad_unit)
232 call da_free_unit(fgat_rad_unit)
233 end if
234
235 if (trace_use) call da_trace_exit("da_qc_amsua")
236
237 end subroutine da_qc_amsua
238
239