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 end if
71
72 ! reject limb obs
73 !------------------------------------------------------
74 scanpos = iv%instid(i)%scanpos(n)
75 if (scanpos <= 8 .or. scanpos >= 83) then
76 iv%instid(i)%tb_qc(:,n) = qc_bad
77 if (iv%instid(i)%proc_domain(n)) &
78 nrej_limb = nrej_limb + 1
79 end if
80
81 ! satzen = rad%satzen
82 ! if (abs(satzen) > 45.0) &
83 ! iv%instid(i)%tb_qc(:,n) = qc_bad
84
85 ! d. check cloud/precipitation
86 !-----------------------------------------------------------
87 if (ob%instid(i)%tb(1,n) > 0.0 .and. &
88 ob%instid(i)%tb(2,n) > 0.0) then
89 si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(2,n)
90 if (si >= 3.0) then
91 iv%instid(i)%tb_qc(:,n) = qc_bad
92 iv%instid(i)%cloud_flag(:,n) = qc_bad
93 if (iv%instid(i)%proc_domain(n)) &
94 nrej_si = nrej_si + 1
95 end if
96 end if
97
98 if (iv%instid(i)%clwp(n) >= 0.2) 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_clw = nrej_clw + 1
103 end if
104
105 ! e. check surface height/pressure
106 !------------------------------------------------------
107 ! sfchgt = ivrad%info(n)%elv
108 ! if (sfchgt >=) then
109 !
110 ! else
111 ! end if
112
113 if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 800.0)) then
114 iv%instid(i)%tb_qc(5,n) = qc_bad
115 if (iv%instid(i)%proc_domain(n)) &
116 nrej_topo = nrej_topo + 1
117 end if
118
119 ! f. check innovation
120 !-----------------------------------------------------------
121 do k = 1, nchan
122 ! absolute departure check
123 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
124 iv%instid(i)%tb_qc(k,n) = qc_bad
125 if (iv%instid(i)%proc_domain(n)) &
126 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
127 end if
128
129 ! relative departure check
130 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*satinfo(i)%error_std(k)) then
131 iv%instid(i)%tb_qc(k,n) = qc_bad
132 if (iv%instid(i)%proc_domain(n)) &
133 nrej_omb_std(k) = nrej_omb_std(k) + 1
134 end if
135
136 ! final QC decision
137 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
138 iv%instid(i)%tb_error(k,n) = 500.0
139 if (iv%instid(i)%proc_domain(n)) &
140 nrej(k) = nrej(k) + 1
141 else
142 if (iv%instid(i)%proc_domain(n)) &
143 ngood(k) = ngood(k) + 1
144 if (use_error_factor_rad) then
145 iv%instid(i)%tb_error(k,n) = &
146 satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
147 else
148 iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
149 end if
150 end if
151 end do ! chan
152 end do ! end loop pixel
153 end if
154
155 ! Do inter-processor communication to gather statistics.
156
157 call da_proc_sum_int (num_proc_domain)
158 call da_proc_sum_int (nrej_mixsurface)
159 call da_proc_sum_int (nrej_windowchanl)
160 call da_proc_sum_int (nrej_si)
161 call da_proc_sum_int (nrej_clw)
162 call da_proc_sum_int (nrej_topo)
163 call da_proc_sum_int (nrej_limb)
164 call da_proc_sum_ints (nrej_omb_abs(:))
165 call da_proc_sum_ints (nrej_omb_std(:))
166 call da_proc_sum_ints (nrej(:))
167 call da_proc_sum_ints (ngood(:))
168
169 if (rootproc) then
170 if (num_fgat_time > 1) then
171 write(filename,'(a,i2.2)') trim(iv%instid(i)%rttovid_string)//'.qcstat',iv%current_ob_time
172 else
173 filename = trim(iv%instid(i)%rttovid_string)//'.qcstat'
174 end if
175 call da_get_unit(fgat_rad_unit)
176 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
177 if (ios /= 0) then
178 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
179 call da_error(__FILE__,__LINE__,message(1:1))
180 end if
181
182 write(fgat_rad_unit, fmt='(/a/)') &
183 'Quality Control Statistics for '//iv%instid(i)%rttovid_string
184 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
185 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
186 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
187 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
188 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
189 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
190 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
191 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
192 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
193 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
194 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
195 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
196 write(fgat_rad_unit,'(10i7)') nrej(:)
197 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
198 write(fgat_rad_unit,'(10i7)') ngood(:)
199
200 close(fgat_rad_unit)
201 call da_free_unit(fgat_rad_unit)
202 end if
203
204 if (trace_use) call da_trace_exit("da_qc_amsub")
205
206 end subroutine da_qc_amsub
207
208