var/da/da_radiance/da_qc_mwhs2.inc

   1 subroutine da_qc_mwhs2 (it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for fy3c mwhs2 data.
   5    ! Wei Sun modified from da_qc_mwhs 2017/10/15
   6    !---------------------------------------------------------------------------
   7
   8    implicit none
   9
  10    integer, intent(in)             :: it         ! outer loop count
  11    integer, intent(in)             :: i          ! sensor index.
  12    integer, intent(in)             :: nchan      ! number of channel
  13    type (y_type),  intent(in)      :: ob         ! Observation structure.
  14    type (iv_type), intent(inout)   :: iv         ! O-B structure.
  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_mwhs2")
  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
  42    num_proc_domain = 0
  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 flags by 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
  63          !  b.  reject channels 1 and 10 over land/sea-ice/snow
  64          !------------------------------------------------------
  65          if (isflg > 0) then
  66             iv%instid(i)%tb_qc(1,n)     = qc_bad
  67             iv%instid(i)%tb_qc(10,n)    = qc_bad
  68             if (iv%instid(i)%info%proc_domain(1,n)) &
  69                nrej_windowchanl = nrej_windowchanl + 1
  70             if (only_sea_rad) iv%instid(i)%tb_qc(:,n)    = qc_bad
  71          end if
  72
  73          ! c. reject limb obs
  74          !------------------------------------------------------
  75          scanpos = iv%instid(i)%scanpos(n)
  76          if (scanpos <= 8 .or. scanpos >= 90) then
  77              iv%instid(i)%tb_qc(:,n)  =  qc_bad
  78              if (iv%instid(i)%info%proc_domain(1,n)) &
  79                  nrej_limb = nrej_limb + 1
  80          end if
  81
  82          ! satzen  = rad%satzen
  83          ! if (abs(satzen) > 45.0) &
  84          !    iv%instid(i)%tb_qc(:,n)  =  qc_bad
  85
  86          !  d. check cloud/precipitation
  87          !-----------------------------------------------------------
  88          if (ob%instid(i)%tb(1,n) > 0.0 .and. ob%instid(i)%tb(10,n) > 0.0 ) then
  89             si = ob%instid(i)%tb(1,n)-ob%instid(i)%tb(10,n)
  90             if (si >= 4.0) then
  91                iv%instid(i)%tb_qc(:,n) = qc_bad
  92                iv%instid(i)%cloud_flag(:,n) = qc_bad
  93                if (iv%instid(i)%info%proc_domain(1,n)) &
  94                   nrej_si = nrej_si + 1
  95             end if
  96          end if
  97
  98          if (iv%instid(i)%clwp(n) >= 0.2) then
  99             iv%instid(i)%tb_qc(:,n) = qc_bad
 100             iv%instid(i)%cloud_flag(:,n) = qc_bad
 101             if (iv%instid(i)%info%proc_domain(1,n)) &
 102                nrej_clw = nrej_clw + 1
 103          end if
 104
 105          !  e. check surface height/pressure
 106          !------------------------------------------------------
 107          ! sfchgt = ivrad%info(n)%elv
 108          ! if (sfchgt >=) then
 109          !
 110          ! else
 111          ! end if
 112
 113         !!!!! chanel 5 when ps ???? !!!!!
 114          if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 800.0)) then
 115             iv%instid(i)%tb_qc(15,n)  = qc_bad
 116             if (iv%instid(i)%info%proc_domain(1,n)) &
 117                nrej_topo = nrej_topo + 1
 118          end if
 119
 120          !  g. check iuse (pre-rejected channels by .info files)
 121          !------------------------------------------------------
 122          do k = 1, nchan
 123            if (satinfo(i)%iuse(k) .eq. -1) &
 124                 iv%instid(i)%tb_qc(k,n) = qc_bad
 125          end do
 126
 127          !  f. check innovation
 128          !-----------------------------------------------------------
 129          do k = 1, nchan
 130           ! absolute departure check
 131             if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
 132                iv%instid(i)%tb_qc(k,n)  = qc_bad
 133                if (iv%instid(i)%info%proc_domain(1,n)) &
 134                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
 135             end if
 136
 137           ! relative departure check
 138             if (use_error_factor_rad) then
 139                  iv%instid(i)%tb_error(k,n) =  &
 140                     satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
 141             else
 142                  iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
 143             end if
 144
 145             if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
 146                iv%instid(i)%tb_qc(k,n)  = qc_bad
 147                if (iv%instid(i)%info%proc_domain(1,n)) &
 148                      nrej_omb_std(k) = nrej_omb_std(k) + 1
 149             end if
 150
 151            ! final QC decision
 152             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
 153                iv%instid(i)%tb_error(k,n) = 500.0
 154                if (iv%instid(i)%info%proc_domain(1,n)) &
 155                   nrej(k) = nrej(k) + 1
 156             else
 157                if (iv%instid(i)%info%proc_domain(1,n)) &
 158                   ngood(k) = ngood(k) + 1
 159             end if
 160          end do ! chan
 161       end do ! end loop pixel
 162
 163    ! Do inter-processor communication to gather statistics.
 164
 165    call da_proc_sum_int (num_proc_domain)
 166    call da_proc_sum_int (nrej_mixsurface)
 167    call da_proc_sum_int (nrej_windowchanl)
 168    call da_proc_sum_int (nrej_si)
 169    call da_proc_sum_int (nrej_clw)
 170    call da_proc_sum_int (nrej_topo)
 171    call da_proc_sum_int (nrej_limb)
 172    call da_proc_sum_ints (nrej_omb_abs(:))
 173    call da_proc_sum_ints (nrej_omb_std(:))
 174    call da_proc_sum_ints (nrej(:))
 175    call da_proc_sum_ints (ngood(:))
 176
 177    if (rootproc) then
 178       if (num_fgat_time > 1) then
 179          write(filename,'(i2.2,a,i2.2)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 180       else
 181          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
 182       end if
 183       call da_get_unit(fgat_rad_unit)
 184       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 185       if (ios /= 0) then
 186          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 187          call da_error(__FILE__,__LINE__,message(1:1))
 188       end if
 189
 190       write(fgat_rad_unit, fmt='(/a/)') &
 191          'Quality Control Statistics for '//iv%instid(i)%rttovid_string
 192       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 193       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 194       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
 195       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
 196       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 197       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
 198       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
 199       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 200       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 201       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 202       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 203       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 204       write(fgat_rad_unit,'(10i7)')     nrej(:)
 205       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 206       write(fgat_rad_unit,'(10i7)')     ngood(:)
 207
 208       close(fgat_rad_unit)
 209       call da_free_unit(fgat_rad_unit)
 210    end if
 211
 212    if (trace_use) call da_trace_exit("da_qc_mwhs2")
 213
 214 end subroutine da_qc_mwhs2