var/da/da_radiance/da_qc_amsua.inc

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