var/da/da_radiance/da_qc_iasi.inc

   1 subroutine da_qc_iasi (it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for metop-a IASI 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    logical   :: lmix
  18    integer   :: n,k,isflg,ios,fgat_rad_unit,sensor_id
  19    integer :: scanpos
  20    integer   :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
  21                 nrej_omb_std(nchan),nrej_limb,     &
  22                 nrej_landsurface,nrej_windowchshort,nrej_windowchlong,    &
  23                 nrej_clw,nrej_sst,nrej_eccloud, num_proc_domain, nrej_mixsurface
  24
  25    real      :: inv_grosscheck
  26
  27    character(len=30)  :: filename
  28
  29    ! iasi Cloud Detection Variables
  30    integer              :: kmin, kmax, ndim
  31    integer, allocatable  :: k_cloud_flag(:) ! cloud flags
  32    ! mmr Cloud Detection Variables
  33    integer              :: kts_100hPa(1), kte_surf,kts_200hPa(1)
  34    integer              :: numrad_local(nchan), numrad_global(nchan)
  35    real                 :: tstore
  36    real                 :: bias_local(nchan), bias_global(nchan)
  37    if (trace_use_dull) call da_trace_entry("da_qc_iasi")
  38
  39    ngood(:)        = 0
  40    nrej(:)         = 0
  41    nrej_omb_abs(:) = 0
  42    nrej_omb_std(:) = 0
  43    nrej_landsurface = 0
  44    nrej_windowchshort= 0
  45    nrej_windowchlong= 0
  46    nrej_sst= 0
  47    nrej_clw        = 0
  48    nrej_eccloud       = 0
  49
  50 !   nrej_mixsurface = 0
  51
  52    nrej_limb       = 0
  53    num_proc_domain = 0
  54    numrad_local    = 0
  55    bias_local      = 0.0
  56    sensor_id = 16
  57
  58
  59
  60       do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
  61
  62          if (iv%instid(i)%info%proc_domain(1,n)) &
  63                num_proc_domain = num_proc_domain + 1
  64
  65          !  0.0  initialise QC by flags assuming good obs
  66          !---------------------------------------------
  67          iv%instid(i)%tb_qc(:,n) = qc_good
  68
  69          !  a.  reject all channels over mixture surface type
  70          !------------------------------------------------------
  71          isflg = iv%instid(i)%isflg(n)
  72          lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
  73          if (lmix) then
  74             iv%instid(i)%tb_qc(:,n)  =  qc_bad
  75             if (iv%instid(i)%info%proc_domain(1,n)) &
  76                nrej_mixsurface = nrej_mixsurface + 1
  77          end if
  78
  79         !  b.  reject limb obs
  80          !------------------------------------------------------
  81          scanpos = iv%instid(i)%scanpos(n)
  82          if (scanpos <= 5 .or. scanpos >= 56) 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          ! c.  reject channels from 565(Reject wavenumber > 2400 )
  88          !------------------------------------------------------
  89          iv%instid(i)%tb_qc(565:nchan,n)  = qc_bad
  90          ! c. cloud detection
  91          !-----------------------------------------------------------
  92          if (iv%instid(i)%clwp(n) >= 0.2) then
  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_clw = nrej_clw + 1
  97          end if
  98
  99          !  d. check innovation
 100          !-----------------------------------------------------------
 101          do k = 1, nchan
 102
 103            !  d.1. check absolute value of innovation
 104            !------------------------------------------------
 105                inv_grosscheck = 15.0
 106                if (use_satcv(2)) inv_grosscheck = 100.0
 107               if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then
 108                  iv%instid(i)%tb_qc(k,n)  = qc_bad
 109                  if (iv%instid(i)%info%proc_domain(1,n)) &
 110                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
 111               end if
 112
 113
 114            !  d.2. check relative value of innovation
 115            !      and assign of the observation error (standard deviation)
 116            !------------------------------------------------------------------------
 117            if (use_error_factor_rad) then         ! if use error tuning factor
 118                iv%instid(i)%tb_error(k,n) = &
 119                    satinfo(i)%error(k)*satinfo(i)%error_factor(k)
 120            else
 121                iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
 122            end if
 123
 124            if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
 125                iv%instid(i)%tb_qc(k,n)  = qc_bad
 126                if (iv%instid(i)%info%proc_domain(1,n)) &
 127                   nrej_omb_std(k) = nrej_omb_std(k) + 1
 128            end if
 129
 130                if ( (iv%instid(i)%tb_qc     (k,n) == qc_good) .and. &
 131                  (iv%instid(i)%cloud_flag(k,n) == qc_good) ) then
 132                 bias_local(k)   = bias_local(k) + ob%instid(i)%tb(k,n) - iv%instid(i)%tb_xb(k,n)
 133                 numrad_local(k) = numrad_local(k) + 1
 134                end if
 135          end do ! chan
 136
 137          !  k. MW03 ecmwf method in McNally AP and Watts PD (2003)
 138          !-----------------------------------------------------------
 139          if (use_clddet==3) then
 140             iv%instid(i)%cloud_flag(:,n) = qc_good
 141             allocate ( k_cloud_flag(1:nchan) )
 142             call da_error(__FILE__,__LINE__, &
 143               (/"Cloud_Detect_Setup is not implemented here, please contact the author of this subroutine."/))
 144             write(unit=stdout,fmt='(A)') 'conducting ECMWF cloud detection setup'
 145 !            CALL Cloud_Detect_Setup(sensor_id)
 146             write(unit=stdout,fmt='(A)') 'finish cloud detection setup'
 147             kmin = iv%instid(i)%kmin_t(n)
 148             kmax = iv%instid(i)%kmax_p(n)
 149             call da_error(__FILE__,__LINE__, &
 150               (/"Cloud_Detect is not implemented here, please contact the author of this subroutine."/))
 151 !            CALL Cloud_Detect( &
 152 !            sensor_id,     &  ! in
 153 !            nchan,     &  ! in
 154 !            iv%instid(i)%ichan,     &  ! in
 155 !            iv%instid(i)%tb_xb(:,n),   &  ! in
 156 !            iv%instid(i)%tb_inv(:,n),     &  ! in
 157 !            iv%instid(i)%p_chan_level(:,n), &  ! in
 158 !            k_cloud_flag, &  ! out
 159 !            kmin,     &  ! in
 160 !            kmax)        ! in  !dm cloud mod
 161
 162             do k = 1, nchan
 163                ! remove channels above the model top
 164                if (iv%instid(i)%p_chan_level(k,n)==0) then
 165                   iv%instid(i)%tb_qc(k,n)  = qc_bad
 166                   if (iv%instid(i)%info%proc_domain(1,n)) &
 167                      nrej_eccloud = nrej_eccloud + 1
 168                end if
 169                if (k_cloud_flag(k) .eq. 1) then
 170                   iv%instid(i)%tb_qc(k,n)  = qc_bad
 171                   if (iv%instid(i)%info%proc_domain(1,n)) &
 172                      nrej_eccloud = nrej_eccloud + 1
 173                 end if
 174             end do
 175             deallocate ( k_cloud_flag )
 176         end if  ! ecmwf
 177
 178       end do ! end loop pixel
 179
 180
 181       ! Do inter-processor communication to gather statistics.
 182
 183       if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
 184          do k = 1, nchan
 185             bias_global(k)   = wrf_dm_sum_real(bias_local(k))
 186             numrad_global(k) = wrf_dm_sum_integer(numrad_local(k))
 187             if (numrad_global(k) > 0) bias_global(k) = bias_global(k) / numrad_global(k)
 188          end do
 189       end if
 190
 191       do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
 192
 193          ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016)
 194          !---------------------------------------------
 195          if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
 196             iv%instid(i)%cloud_flag(:,n) = qc_good
 197
 198             if (rtm_option == rtm_option_rttov) then
 199 #ifdef RTTOV
 200                kte_surf   = iv%instid(i)%nlevels
 201                kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), &
 202                             MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0)
 203                do k=1,nchan
 204                   tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * &
 205                            (ob%instid(i)%tb(k,n) - bias_global(k))
 206                   iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / &
 207                            (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0)
 208                end do
 209 #endif
 210             else if (rtm_option == rtm_option_crtm) then
 211                kte_surf   = kte
 212                kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), &
 213                             MASK = iv%instid(i)%pm(kts:kte,n) < 100.0)
 214
 215                do k = 1, nchan
 216                   CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), &
 217                                             iv%instid(i)%rad_obs(k,n))
 218                end do
 219
 220             end if
 221
 222             ndim = kte_surf - kts_100hPa(1) + 1
 223
 224             call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv)
 225          end if
 226
 227          do k = 1, nchan
 228             if (iv%instid(i)%cloud_flag(k,n) == qc_bad) iv%instid(i)%tb_qc(k,n) = qc_bad
 229          end do
 230
 231          !  2. Check iuse from information file (channel selection)
 232          !-----------------------------------------------------------
 233          do k = 1, nchan
 234             if (satinfo(i)%iuse(k) .eq. -1) &
 235                iv%instid(i)%tb_qc(k,n)  = qc_bad
 236          end do
 237
 238          ! 3. Final QC decision
 239          !---------------------------------------------
 240          do k = 1, nchan
 241             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then  ! bad obs
 242                iv%instid(i)%tb_error(k,n) = 500.0
 243                if (iv%instid(i)%info%proc_domain(1,n)) &
 244                   nrej(k) = nrej(k) + 1
 245             else                                         ! good obs
 246                if (iv%instid(i)%info%proc_domain(1,n)) &
 247                   ngood(k) = ngood(k) + 1
 248             end if
 249          end do ! chan
 250
 251       end do ! end loop pixel
 252
 253    ! Do inter-processor communication to gather statistics.
 254    call da_proc_sum_int  (num_proc_domain)
 255    call da_proc_sum_int  (nrej_landsurface )
 256    call da_proc_sum_int  (nrej_windowchlong)
 257    call da_proc_sum_int  (nrej_windowchshort)
 258    call da_proc_sum_int  (nrej_sst)
 259    call da_proc_sum_int  (nrej_clw  )
 260    call da_proc_sum_int  (nrej_eccloud )
 261    call da_proc_sum_int  (nrej_limb)
 262    call da_proc_sum_ints (nrej_omb_abs(:))
 263    call da_proc_sum_ints (nrej_omb_std(:))
 264    call da_proc_sum_ints (nrej(:))
 265    call da_proc_sum_ints (ngood(:))
 266
 267    if (rootproc) then
 268       if (num_fgat_time > 1) then
 269          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 270       else
 271          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
 272       end if
 273
 274       call da_get_unit(fgat_rad_unit)
 275       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 276       if (ios /= 0) then
 277          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 278          call da_error(__FILE__,__LINE__,message(1:1))
 279       end if
 280
 281       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 282       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 283       write(fgat_rad_unit,'(a20,i7)') ' nrej_landsurface  = ', nrej_landsurface
 284       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchlong = ', nrej_windowchlong
 285       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchshort = ', nrej_windowchshort
 286       write(fgat_rad_unit,'(a20,i7)') ' nrej_sst = ', nrej_sst
 287       !      write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 288       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 289       write(fgat_rad_unit,'(a20,i7)') ' nrej_eccloud        = ', nrej_eccloud
 290       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
 291       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 292       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 293       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 294       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 295       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 296       write(fgat_rad_unit,'(10i7)')     nrej(:)
 297       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 298       write(fgat_rad_unit,'(10i7)')     ngood(:)
 299
 300       close(fgat_rad_unit)
 301       call da_free_unit(fgat_rad_unit)
 302    end if
 303    if (trace_use_dull) call da_trace_exit("da_qc_iasi")
 304
 305 end subroutine da_qc_iasi
 306