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 (iv_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 :: 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_proc_domain, &
23 nrej_limb
24
25 character(len=30) :: filename
26
27 if (trace_use) call da_trace_entry("da_qc_hirs")
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_proc_domain = 0
40
41 do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
42
43 if (iv%instid(i)%info%proc_domain(1,n)) num_proc_domain = num_proc_domain + 1
44
45 ! 0.0 initialise QC by flags assuming good obs
46 !---------------------------------------------
47 iv%instid(i)%tb_qc(:,n) = qc_good
48
49 ! a. reject all channels over mixture surface type
50 !------------------------------------------------------
51 isflg = iv%instid(i)%isflg(n)
52 lmix = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
53 if (lmix) then
54 iv%instid(i)%tb_qc(:,n) = qc_bad
55 if (iv%instid(i)%info%proc_domain(1,n)) &
56 nrej_mixsurface = nrej_mixsurface + 1
57 end if
58
59 ! b. reject all channels over land/sea-ice/snow
60 !------------------------------------------------------
61 if (isflg > 0) then
62 iv%instid(i)%tb_qc(:,n) = qc_bad
63 if (iv%instid(i)%info%proc_domain(1,n)) &
64 nrej_windowchanl = nrej_windowchanl + 1
65 end if
66
67 ! c. reject channels 13,14(above top model 10mb),15
68 !------------------------------------------------------
69 !iv%instid(i)%tb_qc(13:15,n) = qc_bad
70
71 ! reject limb obs
72 !------------------------------------------------------
73 scanpos = iv%instid(i)%scanpos(n)
74 if (scanpos <= 3 .or. scanpos >= 54) then
75 iv%instid(i)%tb_qc(:,n) = qc_bad
76 if (iv%instid(i)%info%proc_domain(1,n)) &
77 nrej_limb = nrej_limb + 1
78 end if
79
80 ! d. cloud detection to be added
81 !-----------------------------------------------------------
82 if (iv%instid(i)%clwp(n) >= 0.2) then
83 iv%instid(i)%tb_qc(:,n) = qc_bad
84 iv%instid(i)%cloud_flag(:,n) = qc_bad
85 if (iv%instid(i)%info%proc_domain(1,n)) &
86 nrej_clw = nrej_clw + 1
87 end if
88
89 ! e. check surface height/pressure
90 !-----------------------------------------------------------
91 ! sfchgt = ivrad%info(n)%elv
92 ! if (sfchgt >=) then
93 ! else
94 ! end if
95
96 !if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.0)) then
97 ! iv%instid(i)%tb_qc(5,n) = qc_bad
98 ! if (iv%instid(i)%info%proc_domain(1,n)) &
99 ! nrej_topo = nrej_topo + 1
100 !end if
101
102 ! g. check iuse from information file (channel selection)
103 !-----------------------------------------------------------
104 do k = 1, nchan
105 if (satinfo(i)%iuse(k) .eq. -1) &
106 iv%instid(i)%tb_qc(k,n) = qc_bad
107 end do
108
109 ! f. check innovation
110 !-----------------------------------------------------------
111 do k = 1, nchan
112
113 ! 1. check absolute value of innovation
114 !------------------------------------------------
115 if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
116 iv%instid(i)%tb_qc(k,n) = qc_bad
117 if (iv%instid(i)%info%proc_domain(1,n)) &
118 nrej_omb_abs(k) = nrej_omb_abs(k) + 1
119 end if
120
121 ! 2. check relative value of innovation
122 ! and assign of the observation error (standard deviation)
123 !------------------------------------------------------------------------
124 if (use_error_factor_rad) then ! if use error tuning factor
125 iv%instid(i)%tb_error(k,n) = &
126 satinfo(i)%error(k)*satinfo(i)%error_factor(k)
127 else
128 iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
129 end if
130
131 if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
132 iv%instid(i)%tb_qc(k,n) = qc_bad
133 if (iv%instid(i)%info%proc_domain(1,n)) &
134 nrej_omb_std(k) = nrej_omb_std(k) + 1
135 end if
136
137 ! 3. Final QC decision
138 !---------------------------------------------
139 if (iv%instid(i)%tb_qc(k,n) == qc_bad) then ! bad obs
140 iv%instid(i)%tb_error(k,n) = 500.0
141 if (iv%instid(i)%info%proc_domain(1,n)) &
142 nrej(k) = nrej(k) + 1
143 else ! good obs
144 if (iv%instid(i)%info%proc_domain(1,n)) &
145 ngood(k) = ngood(k) + 1
146 end if
147 end do ! chan
148 end do ! end loop pixel
149
150 ! Do inter-processor communication to gather statistics.
151 call da_proc_sum_int (num_proc_domain)
152 call da_proc_sum_int (nrej_mixsurface)
153 call da_proc_sum_int (nrej_windowchanl)
154 call da_proc_sum_int (nrej_si )
155 call da_proc_sum_int (nrej_clw)
156 call da_proc_sum_int (nrej_topo)
157 call da_proc_sum_int (nrej_limb)
158 call da_proc_sum_ints (nrej_omb_abs(:))
159 call da_proc_sum_ints (nrej_omb_std(:))
160 call da_proc_sum_ints (nrej(:))
161 call da_proc_sum_ints (ngood(:))
162
163 if (rootproc) then
164 if (num_fgat_time > 1) then
165 write(filename,'(a,i2.2)') 'qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
166 else
167 filename = 'qcstat_'//trim(iv%instid(i)%rttovid_string)
168 end if
169
170 call da_get_unit(fgat_rad_unit)
171 open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
172 if (ios /= 0) then
173 write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
174 call da_error(__FILE__,__LINE__,message(1:1))
175 end if
176
177 write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
178 write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain = ', num_proc_domain
179 write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface = ', nrej_mixsurface
180 write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
181 write(fgat_rad_unit,'(a20,i7)') ' nrej_si = ', nrej_si
182 write(fgat_rad_unit,'(a20,i7)') ' nrej_clw = ', nrej_clw
183 write(fgat_rad_unit,'(a20,i7)') ' nrej_topo = ', nrej_topo
184 write(fgat_rad_unit,'(a20,i7)') ' nrej_limb = ', nrej_limb
185 write(fgat_rad_unit,'(a20)') ' nrej_omb_abs(:) = '
186 write(fgat_rad_unit,'(10i7)') nrej_omb_abs(:)
187 write(fgat_rad_unit,'(a20)') ' nrej_omb_std(:) = '
188 write(fgat_rad_unit,'(10i7)') nrej_omb_std(:)
189 write(fgat_rad_unit,'(a20)') ' nrej(:) = '
190 write(fgat_rad_unit,'(10i7)') nrej(:)
191 write(fgat_rad_unit,'(a20)') ' ngood(:) = '
192 write(fgat_rad_unit,'(10i7)') ngood(:)
193
194 close(fgat_rad_unit)
195 call da_free_unit(fgat_rad_unit)
196 end if
197
198 if (trace_use) call da_trace_exit("da_qc_hirs")
199
200 end subroutine da_qc_hirs
201