wrf svn trunk commit r4103
[wrffire.git] / wrfv2_fire / var / da / da_radiance / da_qc_mhs.inc
blob9ae243c743a492b969e3887bf8a3282a32a4a20c
1 subroutine da_qc_mhs (it, i, nchan, ob, iv)
3    !---------------------------------------------------------------------------
4    ! Purpose: perform quality control for mhs 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.
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_proc_domain,  &
24                 nrej_limb
26    character(len=30)  :: filename
28    if (trace_use) call da_trace_entry("da_qc_mhs")
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
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 flags by 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
62          !  b.  reject channels 1~2 over land/sea-ice/snow
63          !------------------------------------------------------
64          if (isflg > 0) then
65             iv%instid(i)%tb_qc(1:2,n)    = qc_bad
66             if (iv%instid(i)%info%proc_domain(1,n)) &
67                nrej_windowchanl = nrej_windowchanl + 1
68             if (only_sea_rad) iv%instid(i)%tb_qc(:,n)    = qc_bad
69          end if
71          ! reject limb obs 
72          !------------------------------------------------------
73          scanpos = iv%instid(i)%scanpos(n)
74          if (scanpos <= 8 .or. scanpos >= 83) 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
80          ! satzen  = rad%satzen
81          ! if (abs(satzen) > 45.0) &
82          !    iv%instid(i)%tb_qc(:,n)  =  qc_bad
84          !  d. check cloud/precipitation 
85          !-----------------------------------------------------------
86          if (ob%instid(i)%tb(1,n) > 0.0 .and. &
87               ob%instid(i)%tb(2,n) > 0.0) then
88             si = ob%instid(i)%tb(1,n) - ob%instid(i)%tb(2,n)
89             if (si >= 3.0) then
90                iv%instid(i)%tb_qc(:,n) = qc_bad
91                iv%instid(i)%cloud_flag(:,n) = qc_bad
92                if (iv%instid(i)%info%proc_domain(1,n)) &
93                   nrej_si = nrej_si + 1
94             end if
95          end if
97          if (iv%instid(i)%clwp(n) >= 0.2) then
98             iv%instid(i)%tb_qc(:,n) = qc_bad
99             iv%instid(i)%cloud_flag(:,n) = qc_bad
100             if (iv%instid(i)%info%proc_domain(1,n)) &
101                nrej_clw = nrej_clw + 1
102          end if
104          !  e. check surface height/pressure
105          !------------------------------------------------------
106          ! sfchgt = ivrad%info(n)%elv
107          ! if (sfchgt >=) then
108          !
109          ! else 
110          ! end if
112          if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 800.0)) then
113             iv%instid(i)%tb_qc(5,n)  = qc_bad
114             if (iv%instid(i)%info%proc_domain(1,n)) &
115                nrej_topo = nrej_topo + 1
116          end if
118          !  g. check iuse (pre-rejected channels by .info files)
119          !------------------------------------------------------
120          do k = 1, nchan
121            if (satinfo(i)%iuse(k) .eq. -1) &
122                 iv%instid(i)%tb_qc(k,n) = qc_bad
123          end do
125          !  f. check innovation
126          !-----------------------------------------------------------
127          do k = 1, nchan
128           ! absolute departure check
129             if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
130                iv%instid(i)%tb_qc(k,n)  = qc_bad
131                if (iv%instid(i)%info%proc_domain(1,n)) &
132                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
133             end if
135           ! relative departure check
136             if (use_error_factor_rad) then
137                  iv%instid(i)%tb_error(k,n) =  &
138                     satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
139             else
140                  iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
141             end if
143             if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
144                iv%instid(i)%tb_qc(k,n)  = qc_bad
145                if (iv%instid(i)%info%proc_domain(1,n)) &
146                      nrej_omb_std(k) = nrej_omb_std(k) + 1
147             end if
149            ! final QC decision
150             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
151                iv%instid(i)%tb_error(k,n) = 500.0
152                if (iv%instid(i)%info%proc_domain(1,n)) &
153                   nrej(k) = nrej(k) + 1
154             else
155                if (iv%instid(i)%info%proc_domain(1,n)) &
156                   ngood(k) = ngood(k) + 1
157             end if
158          end do ! chan
159       end do ! end loop pixel
161    ! Do inter-processor communication to gather statistics.
163    call da_proc_sum_int (num_proc_domain)
164    call da_proc_sum_int (nrej_mixsurface)
165    call da_proc_sum_int (nrej_windowchanl)
166    call da_proc_sum_int (nrej_si)
167    call da_proc_sum_int (nrej_clw)
168    call da_proc_sum_int (nrej_topo)
169    call da_proc_sum_int (nrej_limb)
170    call da_proc_sum_ints (nrej_omb_abs(:))
171    call da_proc_sum_ints (nrej_omb_std(:))
172    call da_proc_sum_ints (nrej(:))
173    call da_proc_sum_ints (ngood(:))
175    if (rootproc) then
176       if (num_fgat_time > 1) then
177          write(filename,'(i2.2,a,i2.2)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
178       else
179          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
180       end if
181       call da_get_unit(fgat_rad_unit)
182       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
183       if (ios /= 0) then
184          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
185          call da_error(__FILE__,__LINE__,message(1:1))
186       end if
188       write(fgat_rad_unit, fmt='(/a/)') &
189          'Quality Control Statistics for '//iv%instid(i)%rttovid_string
190       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
191       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
192       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
193       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
194       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
195       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
196       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
197       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
198       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
199       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
200       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
201       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
202       write(fgat_rad_unit,'(10i7)')     nrej(:)
203       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
204       write(fgat_rad_unit,'(10i7)')     ngood(:)
206       close(fgat_rad_unit)
207       call da_free_unit(fgat_rad_unit)
208    end if
210    if (trace_use) call da_trace_exit("da_qc_mhs")
212 end subroutine da_qc_mhs