var/da/da_radiance/da_qc_hirs.inc

   1 subroutine da_qc_hirs (it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for HIRS 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    :: 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
  25
  26    character(len=30)  :: filename
  27
  28    if (trace_use) call da_trace_entry("da_qc_hirs")
  29
  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
  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)) num_proc_domain = num_proc_domain + 1
  45
  46          !  0.0  initialise QC by flags assuming good obs
  47          !---------------------------------------------
  48          iv%instid(i)%tb_qc(:,n) = qc_good
  49
  50          !  a.  reject all channels over mixture surface type
  51          !------------------------------------------------------
  52          isflg = iv%instid(i)%isflg(n)
  53          lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
  54          if (lmix) then
  55             iv%instid(i)%tb_qc(:,n)  =  qc_bad
  56             if (iv%instid(i)%info%proc_domain(1,n)) &
  57                nrej_mixsurface = nrej_mixsurface + 1
  58          end if
  59
  60          !  b.  reject all channels over land/sea-ice/snow
  61          !------------------------------------------------------
  62          if (isflg > 0) then
  63             iv%instid(i)%tb_qc(:,n)  = qc_bad
  64             if (iv%instid(i)%info%proc_domain(1,n)) &
  65                nrej_windowchanl = nrej_windowchanl + 1
  66             if (only_sea_rad) iv%instid(i)%tb_qc(:,n)    = qc_bad
  67          end if
  68
  69          !  c.  reject channels 13,14(above top model 10mb),15
  70          !------------------------------------------------------
  71          !iv%instid(i)%tb_qc(13:15,n)  = qc_bad
  72
  73          !    reject limb obs
  74          !------------------------------------------------------
  75          scanpos = iv%instid(i)%scanpos(n)
  76          if (scanpos <= 3 .or. scanpos >= 54) 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          !  d. cloud detection to be added
  83          !-----------------------------------------------------------
  84          if (iv%instid(i)%clwp(n) >= 0.2) then
  85             iv%instid(i)%tb_qc(:,n) = qc_bad
  86             iv%instid(i)%cloud_flag(:,n) = qc_bad
  87             if (iv%instid(i)%info%proc_domain(1,n)) &
  88                nrej_clw = nrej_clw + 1
  89          end if
  90
  91          !  e. check surface height/pressure
  92          !-----------------------------------------------------------
  93          ! sfchgt = ivrad%info(n)%elv
  94          ! if (sfchgt >=) then
  95          ! else
  96          ! end if
  97
  98          !if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.0)) then
  99          !   iv%instid(i)%tb_qc(5,n)  = qc_bad
 100          !   if (iv%instid(i)%info%proc_domain(1,n)) &
 101          !      nrej_topo = nrej_topo + 1
 102          !end if
 103
 104          !  g. check iuse from information file (channel selection)
 105          !-----------------------------------------------------------
 106          do k = 1, nchan
 107             if (satinfo(i)%iuse(k) .eq. -1) &
 108                iv%instid(i)%tb_qc(k,n)  = qc_bad
 109          end do
 110
 111          !  f. check innovation
 112          !-----------------------------------------------------------
 113          do k = 1, nchan
 114
 115          !  1. check absolute value of innovation
 116          !------------------------------------------------
 117             if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
 118                iv%instid(i)%tb_qc(k,n)  = qc_bad
 119                if (iv%instid(i)%info%proc_domain(1,n)) &
 120                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
 121             end if
 122
 123          !  2. check relative value of innovation
 124          !      and assign of the observation error (standard deviation)
 125          !------------------------------------------------------------------------
 126             if (use_error_factor_rad) then         ! if use error tuning factor
 127                iv%instid(i)%tb_error(k,n) = &
 128                   satinfo(i)%error(k)*satinfo(i)%error_factor(k)
 129             else
 130                 iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
 131             end if
 132
 133             if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
 134                iv%instid(i)%tb_qc(k,n)  = qc_bad
 135                if (iv%instid(i)%info%proc_domain(1,n)) &
 136                   nrej_omb_std(k) = nrej_omb_std(k) + 1
 137             end if
 138
 139            ! 3. Final QC decision
 140            !---------------------------------------------
 141             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then  ! bad obs
 142                iv%instid(i)%tb_error(k,n) = 500.0
 143                if (iv%instid(i)%info%proc_domain(1,n)) &
 144                   nrej(k) = nrej(k) + 1
 145             else                                         ! good obs
 146                if (iv%instid(i)%info%proc_domain(1,n)) &
 147                   ngood(k) = ngood(k) + 1
 148             end if
 149          end do ! chan
 150       end do ! end loop pixel
 151
 152    ! Do inter-processor communication to gather statistics.
 153    call da_proc_sum_int (num_proc_domain)
 154    call da_proc_sum_int (nrej_mixsurface)
 155    call da_proc_sum_int (nrej_windowchanl)
 156    call da_proc_sum_int (nrej_si )
 157    call da_proc_sum_int (nrej_clw)
 158    call da_proc_sum_int (nrej_topo)
 159    call da_proc_sum_int (nrej_limb)
 160    call da_proc_sum_ints (nrej_omb_abs(:))
 161    call da_proc_sum_ints (nrej_omb_std(:))
 162    call da_proc_sum_ints (nrej(:))
 163    call da_proc_sum_ints (ngood(:))
 164
 165    if (rootproc) then
 166       if (num_fgat_time > 1) then
 167          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 168       else
 169          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
 170       end if
 171
 172       call da_get_unit(fgat_rad_unit)
 173       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 174       if (ios /= 0) then
 175          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 176          call da_error(__FILE__,__LINE__,message(1:1))
 177       end if
 178
 179       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 180       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 181       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 182       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
 183       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
 184       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 185       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
 186       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
 187       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 188       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 189       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 190       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 191       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 192       write(fgat_rad_unit,'(10i7)')     nrej(:)
 193       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 194       write(fgat_rad_unit,'(10i7)')     ngood(:)
 195
 196       close(fgat_rad_unit)
 197       call da_free_unit(fgat_rad_unit)
 198    end if
 199
 200    if (trace_use) call da_trace_exit("da_qc_hirs")
 201
 202 end subroutine da_qc_hirs
 203