da_qc_amsub.inc
References to this file elsewhere.
1 subroutine da_qc_amsub (i, nchan, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for amsub 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 ! local variables
15 integer :: n,scanpos,k,isflg,ios,fgat_rad_unit
16 logical :: lmix
17 real :: si
18 ! real :: satzen
19 integer :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
20 nrej_omb_std(nchan), &
21 nrej_mixsurface,nrej_windowchanl, nrej_si, &
22 nrej_clw,nrej_topo,num_rad, num_proc_domain, &
23 nrej_limb
24
25 character(len=30) :: filename
26
27 if (trace_use) call da_trace_entry("da_qc_amsub")
28
29 ngood(:) = 0
30 nrej(:) = 0
31 nrej_omb_abs(:) = 0
32 nrej_omb_std(:) = 0
33 nrej_mixsurface = 0
34 nrej_windowchanl= 0
35 nrej_si = 0
36 nrej_clw = 0
37 nrej_topo = 0
38 nrej_limb = 0
39 num_rad = iv%ob_numb(iv%current_ob_time)%radiance(i)- &
40 iv%ob_numb(iv%current_ob_time-1)%radiance(i)
41
42 num_proc_domain = 0
43
44 if (num_rad > 0) then
45 do n= iv%ob_numb(iv%current_ob_time-1)%radiance(i)+1, iv%ob_numb(iv%current_ob_time)%radiance(i)
46
47 if (iv%instid(i)%proc_domain(n)) &
48 num_proc_domain = num_proc_domain + 1
49
50 ! 0.0 initialise QC flags by assuming good obs
51 !---------------------------------------------
52 iv%instid(i)%tb_qc(:,n) = qc_good
53
54 ! a. reject all channels over mixture surface type
55 !------------------------------------------------------
56 isflg = iv%instid(i)%isflg(n)
57 lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
58 if (lmix) then
59 iv%instid(i)%tb_qc(:,n) = qc_bad
60 if (iv%instid(i)%proc_domain(n)) &
61 nrej_mixsurface = nrej_mixsurface + 1
62 end if
63
64 ! b. reject channels 1~2 over land/sea-ice/snow
65 !------------------------------------------------------
66 if (isflg > 0) then
67 iv%instid(i)%tb_qc(1:2,n) = qc_bad
68 if (iv%instid(i)%proc_domain(n)) &
69 nrej_windowchanl = nrej_windowchanl + 1
70 if (only_sea_rad) iv%instid(i)%tb_qc(:,n) = qc_bad
71 end if
72
73 ! reject limb obs
74 !------------------------------------------------------
75 scanpos = iv%instid(i)%scanpos(n)
76 if (scanpos <= 8 .or. scanpos >= 83) then
77 iv%instid(i)%tb_qc(:,n) = qc_bad
78 if (iv%instid(i)%proc_domain(n)) &
79 nrej_limb = nrej_limb + 1
80 end if
81
82 ! satzen = rad%satzen
83 ! if (abs(satzen) > 45.0) &
84 ! iv%instid(i)%tb_qc(:,n) = qc_bad
85
86 ! d. check cloud/precipitation
87 !-----------------------------------------------------------
88 if (ob%instid(i)%tb(1,n) > 0.0 .and. &
89 ob%instid(i)%tb(2,n) > 0.0) then
90 si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(2,n)
91 if (si >= 3.0) then
92 iv%instid(i)%tb_qc(:,n) = qc_bad
93 iv%instid(i)%cloud_flag(:,n) = qc_bad
94 if (iv%instid(i)%proc_domain(n)) &
95 nrej_si = nrej_si + 1
96 end if
97 end if
98
99 if (iv%instid(i)%clwp(n) >= 0.2) then
100 iv%instid(i)%tb_qc(:,n) = qc_bad
101 iv%instid(i)%cloud_flag(:,n) = qc_bad
102 if (iv%instid(i)%proc_domain(n)) &
103 nrej_clw = nrej_clw + 1
104 end if
105
106 ! e. check surface height/pressure
107 !------------------------------------------------------
108 ! sfchgt = ivrad%info(n)%elv
109 ! if (sfchgt >=) then
110 !
111 ! else
112 ! end if
113
114 if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 800.0)) then
115 iv%instid(i)%tb_qc(5,n) = qc_bad
116 if (iv%instid(i)%proc_domain(n)) &
117 nrej_topo = nrej_topo + 1
118 end if
119
120 ! g. check iuse (pre-rejected channels by .info files)
121 !------------------------------------------------------
122 do k = 1, nchan
123 if (satinfo(i)%iuse(k) .eq. -1) &
124 iv%instid(i)%tb_qc(k,n) = qc_bad
125 end do
126
127 ! f. check innovation
128 !-----------------------------------------------------------
129 do k = 1, nchan
130 ! absolute departure check
131 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
132 iv%instid(i)%tb_qc(k,n) = qc_bad
133 if (iv%instid(i)%proc_domain(n)) &
134 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
135 end if
136
137 ! relative departure check
138 if (use_error_factor_rad) then
139 iv%instid(i)%tb_error(k,n) = &
140 satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
141 else
142 iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
143 end if
144
145 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
146 iv%instid(i)%tb_qc(k,n) = qc_bad
147 if (iv%instid(i)%proc_domain(n)) &
148 nrej_omb_std(k) = nrej_omb_std(k) + 1
149 end if
150
151 ! final QC decision
152 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
153 iv%instid(i)%tb_error(k,n) = 500.0
154 if (iv%instid(i)%proc_domain(n)) &
155 nrej(k) = nrej(k) + 1
156 else
157 if (iv%instid(i)%proc_domain(n)) &
158 ngood(k) = ngood(k) + 1
159 end if
160 end do ! chan
161 end do ! end loop pixel
162 end if
163
164 ! Do inter-processor communication to gather statistics.
165
166 call da_proc_sum_int (num_proc_domain)
167 call da_proc_sum_int (nrej_mixsurface)
168 call da_proc_sum_int (nrej_windowchanl)
169 call da_proc_sum_int (nrej_si)
170 call da_proc_sum_int (nrej_clw)
171 call da_proc_sum_int (nrej_topo)
172 call da_proc_sum_int (nrej_limb)
173 call da_proc_sum_ints (nrej_omb_abs(:))
174 call da_proc_sum_ints (nrej_omb_std(:))
175 call da_proc_sum_ints (nrej(:))
176 call da_proc_sum_ints (ngood(:))
177
178 if (rootproc) then
179 if (num_fgat_time > 1) then
180 write(filename,'(a,i2.2)') 'qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%current_ob_time
181 else
182 filename = 'qcstat_'//trim(iv%instid(i)%rttovid_string)
183 end if
184 call da_get_unit(fgat_rad_unit)
185 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
186 if (ios /= 0) then
187 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
188 call da_error(__FILE__,__LINE__,message(1:1))
189 end if
190
191 write(fgat_rad_unit, fmt='(/a/)') &
192 'Quality Control Statistics for '//iv%instid(i)%rttovid_string
193 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
194 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
195 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
196 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
197 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
198 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
199 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
200 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
201 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
202 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
203 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
204 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
205 write(fgat_rad_unit,'(10i7)') nrej(:)
206 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
207 write(fgat_rad_unit,'(10i7)') ngood(:)
208
209 close(fgat_rad_unit)
210 call da_free_unit(fgat_rad_unit)
211 end if
212
213 if (trace_use) call da_trace_exit("da_qc_amsub")
214
215 end subroutine da_qc_amsub
216
217