da_qc_hirs.inc
References to this file elsewhere.
1 subroutine da_qc_hirs (i, nchan, ob, iv)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: perform quality control for HIRS 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_hirs")
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
67 ! b. reject all channels over land/sea-ice/snow
68 !------------------------------------------------------
69 if (isflg > 0) then
70 iv%instid(i)%tb_qc(:,n) = qc_bad
71 if (iv%instid(i)%proc_domain(n)) &
72 nrej_windowchanl = nrej_windowchanl + 1
73 end if
74
75 ! c. reject channels 13,14(above top model 10mb),15
76 !------------------------------------------------------
77 !iv%instid(i)%tb_qc(13:15,n) = qc_bad
78
79 ! reject limb obs
80 !------------------------------------------------------
81 scanpos = iv%instid(i)%scanpos(n)
82 if (scanpos <= 3 .or. scanpos >= 54) then
83 iv%instid(i)%tb_qc(:,n) = qc_bad
84 if (iv%instid(i)%proc_domain(n)) &
85 nrej_limb = nrej_limb + 1
86 end if
87
88 ! d. cloud detection to be added
89 !-----------------------------------------------------------
90 if (iv%instid(i)%clwp(n) >= 0.2) 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_clw = nrej_clw + 1
95 end if
96
97 ! e. check surface height/pressure
98 !-----------------------------------------------------------
99 ! sfchgt = ivrad%info(n)%elv
100 ! if (sfchgt >=) then
101 ! else
102 ! end if
103
104 !if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.)) then
105 ! iv%instid(i)%tb_qc(5,n) = qc_bad
106 ! if (iv%instid(i)%proc_domain(n)) &
107 ! nrej_topo = nrej_topo + 1
108 !end if
109
110 ! g. check iuse from information file (channel selection)
111 !-----------------------------------------------------------
112 do k = 1, nchan
113 if (satinfo(i)%iuse(k) .eq. -1) &
114 iv%instid(i)%tb_qc(k,n) = qc_bad
115 end do
116
117 ! f. check innovation
118 !-----------------------------------------------------------
119 do k = 1, nchan
120
121 ! 1. check absolute value of innovation
122 !------------------------------------------------
123 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.) 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 ! 2. check relative value of innovation
130 ! and assign of the observation error (standard deviation)
131 !------------------------------------------------------------------------
132 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*satinfo(i)%error(k)) then
133 iv%instid(i)%tb_qc(k,n) = qc_bad
134 if (iv%instid(i)%proc_domain(n)) &
135 nrej_omb_std(k) = nrej_omb_std(k) + 1
136 end if
137
138 ! 3. Final QC decision
139 !---------------------------------------------
140 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then ! bad obs
141 iv%instid(i)%tb_error(k,n) = 500.0
142 if (iv%instid(i)%proc_domain(n)) &
143 nrej(k) = nrej(k) + 1
144 else ! good obs
145 if (iv%instid(i)%proc_domain(n)) &
146 ngood(k) = ngood(k) + 1
147 if (use_error_factor_rad) then ! if use error tuning factor
148 iv%instid(i)%tb_error(k,n) = &
149 satinfo(i)%error(k)*satinfo(i)%error_factor(k)
150 else
151 iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
152 end if
153 end if
154 end do ! chan
155 end do ! end loop pixel
156 end if
157
158 ! Do inter-processor communication to gather statistics.
159 call da_proc_sum_int (num_proc_domain)
160 call da_proc_sum_int (nrej_mixsurface)
161 call da_proc_sum_int (nrej_windowchanl)
162 call da_proc_sum_int (nrej_si )
163 call da_proc_sum_int (nrej_clw)
164 call da_proc_sum_int (nrej_topo)
165 call da_proc_sum_int (nrej_limb)
166 call da_proc_sum_ints (nrej_omb_abs(:))
167 call da_proc_sum_ints (nrej_omb_std(:))
168 call da_proc_sum_ints (nrej(:))
169 call da_proc_sum_ints (ngood(:))
170
171 if (rootproc) then
172 if (num_fgat_time > 1) then
173 write(filename,'(a,i2.2)') trim(iv%instid(i)%rttovid_string)//'.qcstat_',iv%current_ob_time
174 else
175 filename = trim(iv%instid(i)%rttovid_string)//'.qcstat'
176 end if
177
178 call da_get_unit(fgat_rad_unit)
179 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
180 if (ios /= 0) Then
181 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
182 call da_error(__FILE__,__LINE__,message(1:1))
183 end if
184
185 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
186 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
187 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
188 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
189 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
190 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
191 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
192 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
193 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
194 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
195 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
196 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
197 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
198 write(fgat_rad_unit,'(10i7)') nrej(:)
199 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
200 write(fgat_rad_unit,'(10i7)') ngood(:)
201
202 close(fgat_rad_unit)
203 call da_free_unit(fgat_rad_unit)
204 end if
205
206 if (trace_use) call da_trace_exit("da_qc_hirs")
207
208 end subroutine da_qc_hirs
209