var/da/da_radiance/da_qc_atms.inc

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