var/da/da_radiance/da_qc_goesimg.inc

   1 subroutine da_qc_goesimg(it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for GOES-image radiance data.
   5    !
   6    ! Method: Yang et al., 2017: Impact of assimilating GOES imager
   7    !          clear-sky radiance with a rapid refresh assimilation
   8    !          system for convection-permitting forecast over Mexico.
   9    !          J. Geophys. Res. Atmos., 122, 5472–5490
  10    !---------------------------------------------------------------------------
  11
  12    implicit none
  13
  14    integer, intent(in)             :: it         ! outer loop count
  15    integer, intent(in)             :: i          ! sensor index.
  16    integer, intent(in)             :: nchan      ! number of channel
  17    type (y_type),  intent(in)      :: ob         ! Observation structure.
  18    type (iv_type), intent(inout)   :: iv         ! O-B structure.
  19
  20
  21    ! local variables
  22    logical   :: lmix,lcould_read
  23    real      :: satzen
  24    integer   :: n,k,isflg,ios,fgat_rad_unit,sensor_id
  25    integer   :: scanpos
  26    integer   :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
  27                 nrej_omb_std(nchan),     &
  28                 nrej_clw,nrej_eccloud, num_proc_domain, nrej_mixsurface
  29
  30    real      :: inv_grosscheck
  31
  32    character(len=30)  :: filename
  33
  34    ! mmr or pf Cloud Detection Variables
  35    integer              :: kts_100hPa(1), kte_surf, ndim
  36    integer              :: numrad_local(nchan), numrad_global(nchan)
  37    real                 :: tstore
  38    real                 :: bias_local(nchan), bias_global(nchan)
  39    integer              :: kmin, kmax
  40    integer, allocatable  :: k_cloud_flag(:) ! cloud flags
  41    if (trace_use_dull) call da_trace_entry("da_qc_goesimg.inc")
  42
  43    ngood(:)        = 0
  44    nrej(:)         = 0
  45    nrej(:)         = 0
  46    nrej_omb_abs(:) = 0
  47    nrej_omb_std(:) = 0
  48    nrej_clw        = 0
  49    nrej_eccloud    = 0
  50    nrej_mixsurface = 0
  51    num_proc_domain = 0
  52    sensor_id = 22
  53
  54    do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
  55
  56       if (iv%instid(i)%info%proc_domain(1,n)) &
  57          num_proc_domain = num_proc_domain + 1
  58
  59
  60       !  0.0  initialise QC by flags assuming good obs
  61       !---------------------------------------------
  62       iv%instid(i)%tb_qc(:,n) = qc_good
  63
  64       !  a.  reject all channels over mixture surface type
  65       !------------------------------------------------------
  66       isflg = iv%instid(i)%isflg(n)
  67       lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
  68       if (lmix) then
  69          iv%instid(i)%tb_qc(:,n)  =  qc_bad
  70          if (iv%instid(i)%info%proc_domain(1,n)) &
  71             nrej_mixsurface = nrej_mixsurface + 1
  72       end if
  73
  74       if (isflg > 0) then    ! if not over water
  75          do k = 1, nchan     ! IR window channel only used over water
  76             if ( k .ne. 2 ) then
  77                if (only_sea_rad) iv%instid(i)%tb_qc(k,n)  = qc_bad
  78             end if
  79          end do
  80       end if
  81
  82       ! b. cloud detection
  83       !-----------------------------------------------------------
  84       if (.not.crtm_cloud) then
  85          if (iv%instid(i)%clwp(n) >= 0.2) then
  86             iv%instid(i)%tb_qc(:,n) = qc_bad
  87             if (iv%instid(i)%info%proc_domain(1,n)) &
  88                nrej_clw = nrej_clw + 1
  89          end if
  90          !if (imager_format.eq.2) then  ! if CLASS NC GVAR data
  91             if (iv%instid(i)%landsea_mask(n) == 0 ) then
  92                if (iv%instid(i)%tb_xb(3,n)-ob%instid(i)%tb(3,n)>3.5) then
  93                   iv%instid(i)%tb_qc(:,n) = qc_bad
  94                   if (iv%instid(i)%info%proc_domain(1,n)) &
  95                      nrej_eccloud = nrej_eccloud + 1
  96                end if
  97             else
  98                if (iv%instid(i)%tb_xb(3,n)-ob%instid(i)%tb(3,n)>2.5) then
  99                   iv%instid(i)%tb_qc(:,n) = qc_bad
 100                   if (iv%instid(i)%info%proc_domain(1,n)) &
 101                      nrej_eccloud = nrej_eccloud + 1
 102                end if
 103             end if
 104          !else                          ! if CIMSS HDF data
 105          !   if (iv%instid(i)%cloud_flag(1,n) >= 1)then ! only use abs clear pixel
 106          !       iv%instid(i)%tb_qc(:,n) = qc_bad
 107          !      if (iv%instid(i)%info%proc_domain(1,n)) &
 108          !           nrej_eccloud = nrej_eccloud + 1
 109          !   end if
 110          !end if
 111
 112          ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016)
 113          !---------------------------------------------
 114          if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
 115             iv%instid(i)%cloud_flag(:,n) = qc_good
 116
 117             if (rtm_option == rtm_option_rttov) then
 118 #ifdef RTTOV
 119                kte_surf   = iv%instid(i)%nlevels
 120                kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), &
 121                             MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0)
 122                do k=1,nchan
 123                   tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * &
 124                            (ob%instid(i)%tb(k,n) - bias_global(k))
 125                   iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / &
 126                            (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0)
 127                end do
 128 #endif
 129             else if (rtm_option == rtm_option_crtm) then
 130                kte_surf   = kte
 131                kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), &
 132                             MASK = iv%instid(i)%pm(kts:kte,n) < 100.0)
 133
 134                do k = 1, nchan
 135                   CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), &
 136                                             iv%instid(i)%rad_obs(k,n))
 137                end do
 138
 139             end if
 140
 141             ndim = kte_surf - kts_100hPa(1) + 1
 142
 143             call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv)
 144          end if
 145
 146          do k = 1, nchan
 147             if (iv%instid(i)%cloud_flag(k,n) == qc_bad) iv%instid(i)%tb_qc(k,n) = qc_bad
 148          end do
 149       end if   !not crtm
 150
 151       !  c. check innovation
 152       !-----------------------------------------------------------
 153       do k = 1, nchan
 154
 155          !  c.1. check absolute value of innovation
 156          !------------------------------------------------
 157          if (.not.crtm_cloud) then
 158             inv_grosscheck = 15.0
 159             if (use_satcv(2)) inv_grosscheck = 100.0
 160             if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then
 161                iv%instid(i)%tb_qc(k,n)  = qc_bad
 162                if (iv%instid(i)%info%proc_domain(1,n)) &
 163                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
 164             end if
 165          end if
 166
 167          !  c.2. check relative value of innovation
 168          !      and assign of the observation error (standard deviation)
 169          !------------------------------------------------------------------------
 170          if (use_error_factor_rad) then         ! if use error tuning factor
 171             iv%instid(i)%tb_error(k,n) = &
 172                satinfo(i)%error(k)*satinfo(i)%error_factor(k)
 173          else
 174             iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
 175          end if
 176
 177          if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
 178             iv%instid(i)%tb_qc(k,n)  = qc_bad
 179             if (iv%instid(i)%info%proc_domain(1,n)) &
 180                nrej_omb_std(k) = nrej_omb_std(k) + 1
 181          end if
 182
 183       end do ! chan
 184
 185
 186       !  2. Check iuse from information file (channel selection)
 187       !-----------------------------------------------------------
 188       do k = 1, nchan
 189          if (satinfo(i)%iuse(k) .eq. -1) &
 190             iv%instid(i)%tb_qc(k,n)  = qc_bad
 191       end do
 192
 193       ! 3. Final QC decision
 194       !---------------------------------------------
 195       do k = 1, nchan
 196          if (iv%instid(i)%tb_qc(k,n) == qc_bad) then  ! bad obs
 197             iv%instid(i)%tb_error(k,n) = 500.0
 198             if (iv%instid(i)%info%proc_domain(1,n)) &
 199                nrej(k) = nrej(k) + 1
 200          else                                         ! good obs
 201             if (iv%instid(i)%info%proc_domain(1,n)) &
 202                ngood(k) = ngood(k) + 1
 203          end if
 204       end do ! chan
 205    end do ! end loop pixel
 206
 207    ! Do inter-processor communication to gather statistics.
 208    call da_proc_sum_int  (num_proc_domain)
 209    call da_proc_sum_int  (nrej_mixsurface)
 210    call da_proc_sum_int  (nrej_clw)
 211    call da_proc_sum_int  (nrej_eccloud)
 212    call da_proc_sum_ints (nrej_omb_abs(:))
 213    call da_proc_sum_ints (nrej_omb_std(:))
 214    call da_proc_sum_ints (nrej(:))
 215    call da_proc_sum_ints (ngood(:))
 216
 217    if (rootproc) then
 218       if (num_fgat_time > 1) then
 219          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 220       else
 221          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
 222       end if
 223
 224       call da_get_unit(fgat_rad_unit)
 225       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 226       if (ios /= 0) then
 227          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 228          call da_error(__FILE__,__LINE__,message(1:1))
 229       end if
 230
 231       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 232       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 233       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 234       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 235       write(fgat_rad_unit,'(a20,i7)') ' nrej_eccloud     = ', nrej_eccloud
 236       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 237       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 238       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 239       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 240       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 241       write(fgat_rad_unit,'(10i7)')     nrej(:)
 242       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 243       write(fgat_rad_unit,'(10i7)')     ngood(:)
 244
 245       close(fgat_rad_unit)
 246       call da_free_unit(fgat_rad_unit)
 247    end if
 248    if (trace_use_dull) call da_trace_exit("da_qc_goesimg.inc")
 249
 250 end subroutine da_qc_goesimg
 251