wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / var / da / da_radiance / da_qc_amsua.inc
blobbdc321e88ca3fc321ee45a64f08333e4186e68c9
1 subroutine da_qc_amsua (it, i, nchan, ob, iv)
3    !---------------------------------------------------------------------------
4    ! Purpose: perform quality control for amsua data.
5    !---------------------------------------------------------------------------
7    implicit none
9    integer, intent(in)             :: it         ! outer loop count
10    integer, intent(in)             :: i          ! sensor index.
11    integer, intent(in)             :: nchan      ! number of channel
12    type (y_type),  intent(in)      :: ob         ! Observation structure.
13    type (iv_type), intent(inout)   :: iv         ! O-B structure.
16    ! local variables
17    integer   :: n,scanpos,k,isflg,ios,fgat_rad_unit
18    logical   :: lmix
19    real      :: si
20    ! real    :: satzen
21    integer   :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
22                 nrej_omb_std(nchan),      &
23                 nrej_mixsurface,nrej_windowchanl, nrej_si,    &
24                 nrej_clw,nrej_topo, num_proc_domain,  &
25                 nrej_limb
27    character(len=30)  :: filename
29    if (trace_use) call da_trace_entry("da_qc_amsua")
31    ngood(:)        = 0
32    nrej(:)         = 0
33    nrej_omb_abs(:) = 0
34    nrej_omb_std(:) = 0
35    nrej_mixsurface = 0
36    nrej_windowchanl= 0
37    nrej_si         = 0
38    nrej_clw        = 0
39    nrej_topo       = 0
40    nrej_limb       = 0
41    num_proc_domain = 0
43       do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
45          if (iv%instid(i)%info%proc_domain(1,n)) &
46                num_proc_domain = num_proc_domain + 1
48          !  0.0  initialise QC by flags assuming good obs
49          !---------------------------------------------
50          iv%instid(i)%tb_qc(:,n) = qc_good
52          !  a.  reject all channels over mixture surface type
53          !------------------------------------------------------
54          isflg = iv%instid(i)%isflg(n)
55          lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
56          if (lmix) then
57             iv%instid(i)%tb_qc(:,n)  =  qc_bad
58             if (iv%instid(i)%info%proc_domain(1,n)) &
59                nrej_mixsurface = nrej_mixsurface + 1
60          end if
61          !  b.  reject channels 1~4 over land/sea-ice/snow
62          !------------------------------------------------------
63          if (isflg > 0) then 
64             iv%instid(i)%tb_qc(1:4,n)  = qc_bad
65             if (iv%instid(i)%info%proc_domain(1,n)) &
66                nrej_windowchanl = nrej_windowchanl + 1
67            ! reject whole pixel if not over sea for global case
68             if (global) iv%instid(i)%tb_qc(:,n)  = qc_bad
69             if (only_sea_rad) iv%instid(i)%tb_qc(:,n)  = qc_bad
70          end if
72          !  c.  reject channels 13,14(above top model 10mb),15 
73          !------------------------------------------------------
74          iv%instid(i)%tb_qc(13:15,n)  = qc_bad
76          !    reject limb obs 
77          !------------------------------------------------------
78          scanpos = iv%instid(i)%scanpos(n)
79          if (scanpos <= 3 .or. scanpos >= 28) then
80             iv%instid(i)%tb_qc(:,n)  =  qc_bad
81             if (iv%instid(i)%info%proc_domain(1,n)) &
82                   nrej_limb = nrej_limb + 1
83          end if
85          ! satzen  = rad%satzen
86          ! if (abs(satzen) > 45.0) iv%instid(i)%tb_qc(:,n)  =  qc_bad
88          !  d. check precipitation 
89          !-----------------------------------------------------------
90          if (ob%instid(i)%tb(1,n) > 0.0 .and. &
91              ob%instid(i)%tb(15,n) > 0.0) then
92             si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(15,n)
93             if (si >= 3.0) then
94                iv%instid(i)%tb_qc(:,n) = qc_bad
95                iv%instid(i)%cloud_flag(:,n) = qc_bad
96                if (iv%instid(i)%info%proc_domain(1,n)) &
97                   nrej_si = nrej_si + 1
98             end if
99          end if
101          if (iv%instid(i)%clwp(n) >= 0.2) then
102             iv%instid(i)%tb_qc(:,n) = qc_bad
103             iv%instid(i)%cloud_flag(:,n) = qc_bad
104             if (iv%instid(i)%info%proc_domain(1,n)) &
105                nrej_clw = nrej_clw + 1
106          end if
108          !   3.1 Estimate Cloud Liquid Water (CLW) in mm over sea
109          !       (Grody etal. 2001, JGR, Equation 5b,7c,7d,9)
110          !---------------------------------------------------------
111          ! if (isflg == 0) then
112          !    coszen =  cos(iv%instid(i)%satzen(n))
113          !    d0     =  8.24-(2.622-1.846*coszen)*coszen
114          !    d1     =  0.754
115          !    d2     =  -2.265
116          !    ts     =  iv%instid(i)%ts(n)
117          !    tb1    =  ob%instid(i)%tb(1,n)
118          !    tb2    =  ob%instid(i)%tb(2,n)
119          !    clw    =  coszen*(d0+d1*log(ts-tb1)+d2*log(ts-tb2))
120          !    clw    =  clw - 0.03
121          ! end if
124          !  e. check surface height/pressure
125          !-----------------------------------------------------------
126          ! sfchgt = ivrad%info(n)%elv
127          ! if (sfchgt >=) then
128          ! else 
129          ! end if
131          if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.0)) then
132             iv%instid(i)%tb_qc(5,n)  = qc_bad
133             if (iv%instid(i)%info%proc_domain(1,n)) &
134                nrej_topo = nrej_topo + 1
135          end if
137          !  g. check iuse
138          !-----------------------------------------------------------
139          do k = 1, nchan
140             if (satinfo(i)%iuse(k) .eq. -1) &
141                iv%instid(i)%tb_qc(k,n)  = qc_bad
142          end do
144          !  f. check innovation
145          !-----------------------------------------------------------
146          do k = 1, nchan
148          ! absolute departure check
149             if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
150                iv%instid(i)%tb_qc(k,n)  = qc_bad
151                if (iv%instid(i)%info%proc_domain(1,n)) &
152                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
153             end if
155          ! relative departure check
156             if (use_error_factor_rad) then
157                iv%instid(i)%tb_error(k,n) = &
158                    satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
159             else
160                iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
161             end if
163             if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
164                 iv%instid(i)%tb_qc(k,n)  = qc_bad
165                 if (iv%instid(i)%info%proc_domain(1,n)) &
166                    nrej_omb_std(k) = nrej_omb_std(k) + 1
167             end if
169          ! final QC decsion
170             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
171                iv%instid(i)%tb_error(k,n) = 500.0
172                if (iv%instid(i)%info%proc_domain(1,n)) &
173                   nrej(k) = nrej(k) + 1
174             else
175                if (iv%instid(i)%info%proc_domain(1,n)) &
176                   ngood(k) = ngood(k) + 1
177             end if
179          end do ! chan
180       end do ! end loop pixel
182    ! Do inter-processor communication to gather statistics.
183    call da_proc_sum_int (num_proc_domain)
184    call da_proc_sum_int (nrej_mixsurface)
185    call da_proc_sum_int (nrej_windowchanl)
186    call da_proc_sum_int (nrej_si )
187    call da_proc_sum_int (nrej_clw)
188    call da_proc_sum_int (nrej_topo)
189    call da_proc_sum_int (nrej_limb)
190    call da_proc_sum_ints (nrej_omb_abs(:))
191    call da_proc_sum_ints (nrej_omb_std(:))
192    call da_proc_sum_ints (nrej(:))
193    call da_proc_sum_ints (ngood(:))
195    if (rootproc) then
196       if (num_fgat_time > 1) then
197          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
198       else
199          write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
200       end if
202       call da_get_unit(fgat_rad_unit)
203       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
204       if (ios /= 0) then
205          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
206          call da_error(__FILE__,__LINE__,message(1:1))
207       end if
209       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
210       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
211       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
212       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
213       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
214       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
215       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
216       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
217       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
218       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
219       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
220       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
221       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
222       write(fgat_rad_unit,'(10i7)')     nrej(:)
223       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
224       write(fgat_rad_unit,'(10i7)')     ngood(:)
226       close(fgat_rad_unit)
227       call da_free_unit(fgat_rad_unit)
228    end if
230    if (trace_use) call da_trace_exit("da_qc_amsua")
232 end subroutine da_qc_amsua