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