var/da/da_radiance/da_qc_ssmis.inc

   1 subroutine da_qc_ssmis (it,i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for ssmis 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,k,isflg,ios,fgat_rad_unit
  18    logical   :: lmix
  19    real      :: si37, q19, q37, term22v
  20    integer   :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
  21                 nrej_omb_std(nchan),      &
  22                 nrej_mixsurface,nrej_windowchanl, nrej_rain, nrej_si,    &
  23                 nrej_clw,nrej_clw_rv,nrej_topo,num_proc_domain
  24
  25    character(len=30)  :: filename
  26
  27    if (trace_use) call da_trace_entry("da_qc_ssmis")
  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_rain       = 0
  36    nrej_si         = 0
  37    nrej_clw        = 0
  38    nrej_clw_rv     = 0
  39    nrej_topo       = 0
  40    num_proc_domain = 0
  41
  42    do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
  43
  44       if (iv%instid(i)%info%proc_domain(1,n)) &
  45             num_proc_domain = num_proc_domain + 1
  46
  47       !  0.0  initialise QC by flags assuming good obs
  48       !---------------------------------------------
  49       iv%instid(i)%tb_qc(:,n) = qc_good
  50
  51       !  a.  reject all channels over mixture surface type
  52       !------------------------------------------------------
  53       isflg = iv%instid(i)%isflg(n)
  54       lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
  55       if (lmix) then
  56          iv%instid(i)%tb_qc(:,n)  =  qc_bad
  57          if (iv%instid(i)%info%proc_domain(1,n)) &
  58             nrej_mixsurface = nrej_mixsurface + 1
  59       end if
  60       !  b.  reject channels 1~2,8 over land/sea-ice/snow
  61       !------------------------------------------------------
  62       if (isflg > 0) then
  63          iv%instid(i)%tb_qc(1:2,n)  = qc_bad
  64          iv%instid(i)%tb_qc(8,n)    = qc_bad
  65          if (iv%instid(i)%info%proc_domain(1,n)) &
  66             nrej_windowchanl = nrej_windowchanl + 1
  67          if (only_sea_rad) iv%instid(i)%tb_qc(:,n)  = qc_bad
  68       end if
  69
  70       !  c. reject rain_flagged data
  71       !------------------------------------------------------
  72       if (iv%instid(i)%rain_flag(n) == 1) then
  73          iv%instid(i)%tb_qc(:,n) = qc_bad
  74          iv%instid(i)%cloud_flag(:,n) = qc_bad
  75          if (iv%instid(i)%info%proc_domain(1,n)) &
  76             nrej_rain = nrej_rain + 1
  77       end if
  78
  79       !  d. check precipitation
  80       !  Ferraro, 1997: Journal of Geophysical Research Vol 102, 16715-16735
  81       !  SI37 = 62.18 + 0.773 * TB19v - TB37v
  82       !  Q19 = -2.70 * ( ln(290-TB19v) - 2.84 - 0.40 * ln(290-TB22v) )
  83       !  Q37 = -1.15 * ( ln(290-TB37v) - 2.99 - 0.32 * ln(290-TB22v) )
  84       !-----------------------------------------------------------
  85
  86       if ( isflg >= 2 ) then
  87          si37 = 62.18 + 0.773 * ob%instid(i)%tb(13,n) - ob%instid(i)%tb(16,n)
  88          if ( si37 >= 5.0 ) then
  89             if ( ob%instid(i)%tb(14,n) >= 264.0 ) then  ! snow check
  90                if ( ob%instid(i)%tb(13,n)-ob%instid(i)%tb(12,n) <= 20.0 ) then  ! desert check
  91                   if ( ob%instid(i)%tb(16,n) <= 253.0 .and.   &
  92                        ob%instid(i)%tb(13,n)-ob%instid(i)%tb(12,n) <= 7.0 ) then  ! arid soil check
  93                      iv%instid(i)%tb_qc(:,n) = qc_bad
  94                      iv%instid(i)%cloud_flag(:,n) = qc_bad
  95                      if (iv%instid(i)%info%proc_domain(1,n)) &
  96                         nrej_si = nrej_si + 1
  97                    end if
  98                end if
  99             end if
 100           end if
 101       end if
 102
 103       if ( isflg <= 1 ) then
 104          if ( ob%instid(i)%tb(14,n) < 290.0 ) then  ! assure positive for log
 105             term22v = log(290.0-ob%instid(i)%tb(14,n))
 106             if ( ob%instid(i)%tb(13,n) < 290.0 ) then  ! assure positive for log
 107                q19 = -2.70*(log(290.0-ob%instid(i)%tb(13,n))-2.84-0.40*term22v)
 108             end if
 109             if ( ob%instid(i)%tb(16,n) < 290.0 ) then  ! assure positive for log
 110                q37 = -1.15*(log(290.0-ob%instid(i)%tb(16,n))-2.99-0.32*term22v)
 111             end if
 112             if ( q19 >= 0.60 .or. q37 >= 0.20 ) then
 113                if ( ob%instid(i)%tb(14,n) >= 44.0+0.85*ob%instid(i)%tb(13,n) .or.  &
 114                     (ob%instid(i)%tb(14,n) <= 264.0 .and.    &
 115                      (ob%instid(i)%tb(14,n)-ob%instid(i)%tb(13,n)) >= 2.0) ) then
 116                   iv%instid(i)%tb_qc(:,n) = qc_bad
 117                   iv%instid(i)%cloud_flag(:,n) = qc_bad
 118                   if (iv%instid(i)%info%proc_domain(1,n)) &
 119                      nrej_clw_rv = nrej_clw_rv + 1      ! clw retrieval
 120                end if
 121             end if
 122          end if
 123       end if
 124
 125       if (iv%instid(i)%clwp(n) >= 0.2) then
 126          iv%instid(i)%tb_qc(:,n) = qc_bad
 127          iv%instid(i)%cloud_flag(:,n) = qc_bad
 128          if (iv%instid(i)%info%proc_domain(1,n)) &
 129             nrej_clw = nrej_clw + 1               ! clw model
 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_rain )
 188    call da_proc_sum_int (nrej_si )
 189    call da_proc_sum_int (nrej_clw_rv)
 190    call da_proc_sum_int (nrej_clw)
 191 !   call da_proc_sum_int (nrej_topo)
 192    call da_proc_sum_ints (nrej_omb_abs(:))
 193    call da_proc_sum_ints (nrej_omb_std(:))
 194    call da_proc_sum_ints (nrej(:))
 195    call da_proc_sum_ints (ngood(:))
 196
 197    if (rootproc) then
 198       if (num_fgat_time > 1) then
 199          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 200       else
 201          write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
 202       end if
 203
 204       call da_get_unit(fgat_rad_unit)
 205       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 206       if (ios /= 0) then
 207          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 208          call da_error(__FILE__,__LINE__,message(1:1))
 209       end if
 210
 211       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 212       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 213       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 214       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
 215       write(fgat_rad_unit,'(a20,i7)') ' nrej_rain        = ', nrej_rain
 216       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
 217       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw_rv      = ', nrej_clw_rv
 218       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 219 !      write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
 220       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 221       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 222       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 223       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 224       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 225       write(fgat_rad_unit,'(10i7)')     nrej(:)
 226       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 227       write(fgat_rad_unit,'(10i7)')     ngood(:)
 228
 229       close(fgat_rad_unit)
 230       call da_free_unit(fgat_rad_unit)
 231    end if
 232
 233    if (trace_use) call da_trace_exit("da_qc_ssmis")
 234
 235 end subroutine da_qc_ssmis
 236
 237