da_qc_airs.inc
References to this file elsewhere.
1 subroutine da_qc_airs (i, nchan, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for AQUA/EOS-2-AIRS 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 !!
21 !! TV change: added new diagnostic local varaibles
22 integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
23 nrej_omb_std(nchan),nrej_limb, &
24 nrej_landsurface,nrej_windowchshort,nrej_windowchlong, nrej_si, &
25 nrej_clw,nrej_sst,nrej_topo,num_rad, num_proc_domain
26 !
27 ! integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
28 ! nrej_omb_std(nchan), &
29 ! nrej_mixsurface,nrej_windowchanl, nrej_si, &
30 ! nrej_clw,nrej_topo,num_rad, num_proc_domain, &
31 ! nrej_limb
32 !
33 real :: SST_model, SST_airs, diffSST
34
35 !! TV end change
36
37
38 character(len=30) :: filename
39
40 if (trace_use) call da_trace_entry("da_qc_airs")
41
42 !! TV change: new diagnosic locals
43 ngood(:) = 0
44 nrej(:) = 0
45 nrej_omb_abs(:) = 0
46 nrej_omb_std(:) = 0
47 nrej_landsurface = 0
48 nrej_windowchshort= 0
49 nrej_windowchlong= 0
50 nrej_sst= 0
51 ! nrej_si = 0
52 nrej_clw = 0
53 nrej_topo = 0
54
55 ! nrej_mixsurface = 0
56 ! nrej_windowchanl= 0
57
58 !! TV end change
59
60 nrej_limb = 0
61 num_rad = iv%ob_numb(iv%current_ob_time)%radiance(i)- &
62 iv%ob_numb(iv%current_ob_time-1)%radiance(i)
63 ! num_rad = iv%instid(i)%num_rad
64 num_proc_domain = 0
65
66 if (num_rad > 0) then
67
68 ! do n= 1, iv%instid(i)%num_rad ! loop for pixel
69 do n= iv%ob_numb(iv%current_ob_time-1)%radiance(i)+1, iv%ob_numb(iv%current_ob_time)%radiance(i)
70
71 if (iv%instid(i)%proc_domain(n)) &
72 num_proc_domain = num_proc_domain + 1
73
74 ! 0.0 initialise QC by flags assuming good obs
75 !---------------------------------------------
76 iv%instid(i)%tb_qc(:,n) = qc_good
77
78 !! TV change: include sea locations only
79 ! a. reject all channels over land and mixture
80 !------------------------------------------------------
81 isflg = iv%instid(i)%isflg(n)
82 if (isflg > 0) then
83 iv%instid(i)%tb_qc(:,n) = -1
84 if (iv%instid(i)%proc_domain(n)) &
85 nrej_landsurface = nrej_landsurface + 1
86 end if
87
88 ! ! a. reject all channels over mixture surface type
89 ! !------------------------------------------------------
90 ! isflg = iv%instid(i)%isflg(n)
91 ! lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
92 ! if (lmix) then
93 ! iv%instid(i)%tb_qc(:,n) = qc_bad
94 ! if (iv%instid(i)%proc_domain(n)) &
95 ! nrej_mixsurface = nrej_mixsurface + 1
96 ! end if
97
98 ! ! b. reject all channels over land/sea-ice/snow
99 ! !------------------------------------------------------
100 ! if (isflg > 0) then
101 ! iv%instid(i)%tb_qc(:,n) = qc_bad
102 ! if (iv%instid(i)%proc_domain(n)) &
103 ! nrej_windowchanl = nrej_windowchanl + 1
104 ! end if
105
106 !! TV end change: include sea locations only
107
108
109 ! reject limb obs
110 !------------------------------------------------------
111 !scanpos = iv%instid(i)%scanpos(n)
112 !if (scanpos <= 3 .or. scanpos >= 88) then
113 ! iv%instid(i)%tb_qc(:,n) = qc_bad
114 ! if (iv%instid(i)%proc_domain(n)) &
115 ! nrej_limb = nrej_limb + 1
116 !end if
117
118
119 !! TV change: QC for clouds
120
121 ! c. Check for model clouds!!!!!!
122 !-----------------------------------------------------------
123
124 if (iv%instid(i)%clwp(n) > 0.05) then
125 iv%instid(i)%tb_qc(:,n) = qc_bad
126 iv%instid(i)%cloud_flag(:,n) = -1
127 if (iv%instid(i)%proc_domain(n)) &
128 nrej_clw = nrej_clw + 1
129 end if
130
131 ! d. Crude check for clouds in obs (assuming obs are used over ocean only)
132 ! Use long wave window channel #914 - 10.662 nm (965.43 cm^-1)
133 ! should be warmer than freezing temperature of the sea
134 !-----------------------------------------------------------
135 !
136 if(ob%instid(i)%tb(129,n) < 271.) then
137 iv%instid(i)%tb_qc(:,n) = qc_bad
138 if (iv%instid(i)%proc_domain(n)) &
139 nrej_windowchlong = nrej_windowchlong + 1
140 end if
141
142 ! e. Check for contaminated obs in warmest near-infrared: Sun contamination during day
143 !-----------------------------------------------------------
144 !
145 SST_airs=ob%instid(i)%tb(272,n) !! short wave window channel #2328 - 3.882 nm (2616.38 cm^-1)
146 if(SST_airs > 307.) then
147 iv%instid(i)%tb_qc(257:281,n) = qc_bad
148 if (iv%instid(i)%proc_domain(n)) &
149 nrej_windowchshort = nrej_windowchshort + 1
150 end if
151
152
153 ! f. Check for cloud free in obs (assuming obs are used over ocean only)
154 ! Criterion: model SST within 2 K of transparent (hottest) short wave window channel
155 ! includes check for contaminated near-infrared
156 !-----------------------------------------------------------
157 !
158 SST_model=iv%instid(i)%ts(n) !! SST in the model
159 diffSST=abs(SST_model-SST_airs)
160 if(SST_airs < 307 .and. diffSST > 2.) then
161 iv%instid(i)%tb_qc(:,n) = qc_bad
162 if (iv%instid(i)%proc_domain(n)) &
163 nrej_sst = nrej_sst + 1
164 end if
165
166 !! TV change
167
168 ! e. check surface height/pressure
169 !-----------------------------------------------------------
170 ! sfchgt = ivrad%info(n)%elv
171 ! if (sfchgt >=) then
172 ! else
173 ! end if
174
175 !if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.)) then
176 ! iv%instid(i)%tb_qc(5,n) = qc_bad
177 ! if (iv%instid(i)%proc_domain(n)) &
178 ! nrej_topo = nrej_topo + 1
179 !end if
180
181 ! g. check iuse from information file (channel selection)
182 !-----------------------------------------------------------
183 !do k = 1, nchan
184 ! if (satinfo(i)%iuse(k) .eq. -1) &
185 ! iv%instid(i)%tb_qc(k,n) = qc_bad
186 !end do
187
188 ! f. check innovation
189 !-----------------------------------------------------------
190 do k = 1, nchan
191
192 ! 1. check absolute value of innovation
193 !------------------------------------------------
194 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.) then
195 iv%instid(i)%tb_qc(k,n) = qc_bad
196 if (iv%instid(i)%proc_domain(n)) &
197 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
198 end if
199
200 ! 2. check relative value of innovation
201 ! and assign of the observation error (standard deviation)
202 !------------------------------------------------------------------------
203 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*satinfo(i)%error(k)) then
204 iv%instid(i)%tb_qc(k,n) = qc_bad
205 if (iv%instid(i)%proc_domain(n)) &
206 nrej_omb_std(k) = nrej_omb_std(k) + 1
207 end if
208
209 ! 3. Final QC decision
210 !---------------------------------------------
211 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then ! bad obs
212 iv%instid(i)%tb_error(k,n) = 500.0
213 if (iv%instid(i)%proc_domain(n)) &
214 nrej(k) = nrej(k) + 1
215 else ! good obs
216 if (iv%instid(i)%proc_domain(n)) &
217 ngood(k) = ngood(k) + 1
218 if (use_error_factor_rad) then ! if use error tuning factor
219 iv%instid(i)%tb_error(k,n) = &
220 satinfo(i)%error(k)*satinfo(i)%error_factor(k)
221 else
222 iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
223 end if
224 end if
225 end do ! chan
226 end do ! end loop pixel
227 end if
228
229 ! Do inter-processor communication to gather statistics.
230 call da_proc_sum_int (num_proc_domain)
231 !! TV change
232 call da_proc_sum_int ( nrej_landsurface )
233 call da_proc_sum_int ( nrej_windowchlong)
234 call da_proc_sum_int ( nrej_windowchshort)
235 call da_proc_sum_int ( nrej_sst)
236 ! call da_proc_sum_int ( nrej_si )
237 !! TV end change
238 call da_proc_sum_int ( nrej_clw )
239 call da_proc_sum_int ( nrej_topo )
240 call da_proc_sum_int (nrej_limb)
241 call da_proc_sum_ints (nrej_omb_abs(:))
242 call da_proc_sum_ints (nrej_omb_std(:))
243 call da_proc_sum_ints (nrej(:))
244 call da_proc_sum_ints (ngood(:))
245
246 if (rootproc) then
247 if (num_fgat_time > 1) then
248 write(filename,'(a,i2.2)') trim(iv%instid(i)%rttovid_string)//'.qcstat_',iv%current_ob_time
249 else
250 filename = trim(iv%instid(i)%rttovid_string)//'.qcstat'
251 end if
252
253 call da_get_unit(fgat_rad_unit)
254 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
255 if (ios /= 0) Then
256 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
257 call da_error(__FILE__,__LINE__,message(1:1))
258 end if
259
260 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
261 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
262 !!TV change
263 write(fgat_rad_unit,'(a20,i7)') ' nrej_landsurface = ', nrej_landsurface
264 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchlong = ', nrej_windowchlong
265 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchshort = ', nrej_windowchshort
266 write(fgat_rad_unit,'(a20,i7)') ' nrej_sst = ', nrej_sst
267 ! write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
268 ! write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
269 ! write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
270 !!TV end change
271 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
272 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
273 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
274 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
275 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
276 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
277 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
278 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
279 write(fgat_rad_unit,'(10i7)') nrej(:)
280 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
281 write(fgat_rad_unit,'(10i7)') ngood(:)
282
283 close(fgat_rad_unit)
284 call da_free_unit(fgat_rad_unit)
285 end if
286
287 if (trace_use) call da_trace_exit("da_qc_airs")
288
289 end subroutine da_qc_airs
290