var/da/da_radiance/da_qc_airs.inc

   1 subroutine da_qc_airs (it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for AQUA/EOS-2-AIRS 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    integer :: scanpos
  19    ! logical   :: lmix
  20    ! real    :: satzen
  21
  22    integer   :: ngood(nchan),nrej(nchan),nrej_omb_abs(nchan), &
  23                 nrej_omb_std(nchan),nrej_limb,     &
  24                 nrej_landsurface,nrej_windowchshort,nrej_windowchlong,    &
  25                 nrej_clw,nrej_sst,nrej_topo, num_proc_domain,sensor_id,nrej_eccloud
  26
  27    real      :: SST_model, SST_airs, SST_pred, diffSST, diffSST2
  28    real      :: inv_grosscheck
  29
  30    character(len=30)  :: filename
  31
  32    ! AIRS Cloud Detection Variables
  33    integer              :: kts_100hPa(1), kte_surf, ndim
  34    integer              :: numrad_local(nchan), numrad_global(nchan)
  35    real                 :: tstore
  36    real                 :: bias_local(nchan), bias_global(nchan)
  37    integer              :: kmin, kmax
  38    integer, allocatable  :: k_cloud_flag(:) ! cloud flags
  39    if (trace_use_dull) call da_trace_entry("da_qc_airs")
  40
  41    ngood(:)        = 0
  42    nrej(:)         = 0
  43    nrej_omb_abs(:) = 0
  44    nrej_omb_std(:) = 0
  45    nrej_landsurface = 0
  46    nrej_windowchshort= 0
  47    nrej_windowchlong= 0
  48    nrej_sst= 0
  49    nrej_clw        = 0
  50    nrej_topo       = 0
  51    nrej_eccloud       = 0
  52    sensor_id = 11
  53 !   nrej_mixsurface = 0
  54
  55    nrej_limb       = 0
  56    num_proc_domain = 0
  57    numrad_local    = 0
  58    bias_local      = 0.0
  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 channels over land/sea-ice/snow and mixture
  70          !------------------------------------------------------------
  71          isflg = iv%instid(i)%isflg(n)
  72          if (isflg > 0) then
  73             ! Check surface emissivity Jacobian
  74             !-----------------------------------
  75             if (rtm_option == rtm_option_crtm .and. use_crtm_kmatrix) then
  76                do k = 1, nchan
  77                   if ( abs(iv%instid(i)%emiss_jacobian(k,n)) > 0.1 ) &
  78                   iv%instid(i)%tb_qc(k,n)  =  qc_bad
  79                end do
  80             end if
  81
  82             if (use_rttov_kmatrix) then
  83                do k = 1, nchan
  84                   if ( abs(iv%instid(i)%emiss_jacobian(k,n)) > 0.1 ) &
  85                   iv%instid(i)%tb_qc(k,n)  =  qc_bad
  86                end do
  87             end if
  88
  89
  90             ! reject all channels
  91             !--------------------
  92             if (only_sea_rad) then
  93                iv%instid(i)%tb_qc(:,n) = qc_bad
  94                if (iv%instid(i)%info%proc_domain(1,n)) &
  95                   nrej_landsurface = nrej_landsurface + 1
  96             end if
  97          end if
  98
  99          ! a (bis) Check T and Q Jacobians for sensitivity to model top
 100          !-----------------------------------------------------------
 101          if (rtm_option == rtm_option_crtm .and. use_crtm_kmatrix) then
 102             do k = 1, nchan
 103               if ( abs(iv%instid(i)%t_jacobian(k,1,n)) > 0.1 * SUM( &
 104                    abs(iv%instid(i)%t_jacobian(k,1:41,n))) ) &
 105                  iv%instid(i)%tb_qc(k,n)  =  qc_bad
 106               if ( abs(iv%instid(i)%q_jacobian(k,1,n)) > 0.1 * SUM( &
 107                    abs(iv%instid(i)%q_jacobian(k,1:41,n))) ) &
 108                  iv%instid(i)%tb_qc(k,n)  =  qc_bad
 109             end do
 110          end if
 111
 112          if (use_rttov_kmatrix) then
 113             do k = 1, nchan
 114               if ( abs(iv%instid(i)%t_jacobian(k,1,n)) > 0.1 * SUM( &
 115                    abs(iv%instid(i)%t_jacobian(k,1:41,n))) ) &
 116                  iv%instid(i)%tb_qc(k,n)  =  qc_bad
 117               if ( abs(iv%instid(i)%q_jacobian(k,1,n)) > 0.1 * SUM( &
 118                    abs(iv%instid(i)%q_jacobian(k,1:41,n))) ) &
 119                  iv%instid(i)%tb_qc(k,n)  =  qc_bad
 120             end do
 121          end if
 122
 123 !         !  a.  reject all channels over mixture surface type
 124 !         !------------------------------------------------------
 125 !         isflg = iv%instid(i)%isflg(n)
 126 !         lmix  = (isflg==4) .or. (isflg==5) .or. (isflg==6) .or. (isflg==7)
 127 !         if (lmix) then
 128 !            iv%instid(i)%tb_qc(:,n)  =  qc_bad
 129 !            if (iv%instid(i)%info%proc_domain(1,n)) &
 130 !               nrej_mixsurface = nrej_mixsurface + 1
 131 !         end if
 132
 133          !  b.  reject limb obs
 134          !------------------------------------------------------
 135          scanpos = iv%instid(i)%scanpos(n)
 136          if (scanpos <= 3 .or. scanpos >= 88) then
 137             iv%instid(i)%tb_qc(:,n)  =  qc_bad
 138             if (iv%instid(i)%info%proc_domain(1,n)) &
 139                   nrej_limb = nrej_limb + 1
 140          end if
 141
 142
 143          !  c. Check for model clouds
 144          !-----------------------------------------------------------
 145          if (iv%instid(i)%clwp(n) > 0.05) then
 146             iv%instid(i)%cloud_flag(:,n) = qc_bad
 147             if (iv%instid(i)%info%proc_domain(1,n)) &
 148                nrej_clw = nrej_clw + 1
 149          end if
 150
 151          !  d. Crude check for clouds in obs (assuming obs are used over ocean only)
 152          !   Use long wave window channel #914 - 10.662 nm (965.43 cm^-1)
 153          !   should be warmer than freezing temperature of the sea
 154          !-----------------------------------------------------------
 155          !
 156          if(ob%instid(i)%tb(129,n) < 271.0) then
 157             iv%instid(i)%cloud_flag(:,n) = qc_bad
 158             if (iv%instid(i)%info%proc_domain(1,n)) &
 159                nrej_windowchlong = nrej_windowchlong + 1
 160          end if
 161
 162          !  e. Check for contaminated obs in warmest near-infrared: Sun contamination during day
 163          !-----------------------------------------------------------
 164          !
 165 !         SST_airs=ob%instid(i)%tb(272,n)   !! short wave window channel #2333 - 3.882 nm (2616.38 cm^-1)
 166 !         if(SST_airs > 307.0) then
 167 !           iv%instid(i)%tb_qc(257:281,n)  = qc_bad
 168 !            if (iv%instid(i)%info%proc_domain(1,n)) &
 169 !               nrej_windowchshort = nrej_windowchshort + 1
 170 !         end if
 171
 172
 173          !  f. Check for cloud free in obs (assuming obs are used over ocean only)
 174          !  Criterion: model SST within 2 K of transparent (hottest) short wave window channel
 175          !             includes check for contaminated near-infrared
 176          !-----------------------------------------------------------
 177          !
 178          SST_model=iv%instid(i)%ts(n)       !! SST in the model
 179 !         diffSST=abs(SST_model-SST_airs)
 180 !         if(iv%instid(i)%solzen(n)>85.0 .and. diffSST > 2.0) then !! night-time only
 181 !            iv%instid(i)%tb_qc(:,n)  = qc_bad
 182 !            if (iv%instid(i)%info%proc_domain(1,n)) &
 183 !               nrej_sst = nrej_sst + 1
 184 !         end if
 185
 186          !  g. Test on SST versus AIRS predicted SST from shortwave and longwave
 187          !  Use channels #791, 914, 1285 and 1301.
 188          !----------------------------------------------------------
 189          !
 190          SST_pred=8.28206-0.97957*ob%instid(i)%tb(126,n)+0.60529*ob%instid(i)%tb(129,n) + &
 191                   1.7444*ob%instid(i)%tb(165,n)-0.40379*ob%instid(i)%tb(166,n)
 192          diffSST2=SST_model-SST_pred
 193
 194          if ((diffSST2<-0.6).or.(diffSST2>3.5)) then
 195             iv%instid(i)%cloud_flag(:,n) = qc_bad
 196             if (iv%instid(i)%info%proc_domain(1,n)) &
 197                nrej_sst = nrej_sst + 1
 198          end if
 199
 200
 201          !  h. Test AIRS/VISNIR cloud fraction (in %)
 202          !  Criterion : 5% cloud coverage within AIRS pixel
 203          !----------------------------------------------------------
 204          !
 205          if ((iv%instid(i)%rain_flag(n)>5).and.(iv%instid(i)%rain_flag(n)<100)) then
 206             iv%instid(i)%cloud_flag(:,n) = qc_bad
 207             if (iv%instid(i)%info%proc_domain(1,n)) &
 208                nrej_sst = nrej_sst + 1
 209          end if
 210
 211
 212          !  i. check surface height/pressure
 213          !-----------------------------------------------------------
 214          ! sfchgt = ivrad%info(n)%elv
 215          ! if (sfchgt >=) then
 216          ! else
 217          ! end if
 218
 219          !if ((isflg .ne. 0) .and. (iv%instid(i)%ps(n) < 850.0)) then
 220          !   iv%instid(i)%tb_qc(5,n)  = qc_bad
 221          !   if (iv%instid(i)%info%proc_domain(1,n)) &
 222          !      nrej_topo = nrej_topo + 1
 223          !end if
 224
 225          !  j. check innovation
 226          !-----------------------------------------------------------
 227          do k = 1, nchan
 228
 229            !  j.1. check absolute value of innovation
 230            !------------------------------------------------
 231             inv_grosscheck = 15.0
 232             if (use_satcv(2)) inv_grosscheck = 100.0
 233             if (abs(iv%instid(i)%tb_inv(k,n)) > inv_grosscheck) then
 234                iv%instid(i)%tb_qc(k,n)  = qc_bad
 235                if (iv%instid(i)%info%proc_domain(1,n)) &
 236                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
 237             end if
 238             if (use_satcv(2)) cycle
 239
 240            !  j.2. check relative value of innovation
 241            !      and assign of the observation error (standard deviation)
 242            !------------------------------------------------------------------------
 243             if (use_error_factor_rad) then         ! if use error tuning factor
 244                iv%instid(i)%tb_error(k,n) = &
 245                    satinfo(i)%error(k)*satinfo(i)%error_factor(k)
 246             else
 247                iv%instid(i)%tb_error(k,n) = satinfo(i)%error(k)
 248             end if
 249
 250             ! M-estimator using Huber function (with k=sigmaO)
 251 !           if (abs(iv%instid(i)%tb_inv(k,n)) > iv%instid(i)%tb_error(k,n)) &
 252 !              iv%instid(i)%tb_error(k,n) = sqrt( &
 253 !              iv%instid(i)%tb_error(k,n) * abs(iv%instid(i)%tb_inv(k,n)) )
 254
 255             if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
 256                iv%instid(i)%tb_qc(k,n)  = qc_bad
 257                if (iv%instid(i)%info%proc_domain(1,n)) &
 258                   nrej_omb_std(k) = nrej_omb_std(k) + 1
 259             end if
 260
 261             if ( (iv%instid(i)%tb_qc     (k,n) == qc_good) .and. &
 262                  (iv%instid(i)%cloud_flag(k,n) == qc_good) ) then
 263                bias_local(k)   = bias_local(k) + ob%instid(i)%tb(k,n) - iv%instid(i)%tb_xb(k,n)
 264                numrad_local(k) = numrad_local(k) + 1
 265             end if
 266
 267          end do ! chan
 268
 269          !  k. MW03 ecmwf method in McNally AP and Watts PD (2003)
 270          !-----------------------------------------------------------
 271          if (use_clddet==3) then
 272             iv%instid(i)%cloud_flag(:,n) = qc_good
 273             allocate ( k_cloud_flag(1:nchan) )
 274             call da_error(__FILE__,__LINE__, &
 275               (/"Cloud_Detect_Setup is not implemented here, please contact the author of this subroutine."/))
 276 !            CALL Cloud_Detect_Setup(sensor_id)
 277             write(unit=stdout,fmt='(A)') 'finish cloud detection setup for airs'
 278             kmin = iv%instid(i)%kmin_t(n)
 279             kmax = iv%instid(i)%kmax_p(n)
 280             call da_error(__FILE__,__LINE__, &
 281                (/"Cloud_Detect is not implemented here, please contact the author of this subroutine."/))
 282 !            CALL Cloud_Detect( &
 283 !            sensor_id,     &  ! in
 284 !            nchan,     &  ! in
 285 !            iv%instid(i)%ichan,     &  ! in
 286 !            iv%instid(i)%tb_xb(:,n),   &  ! in
 287 !            iv%instid(i)%tb_inv(:,n),     &  ! in
 288 !            iv%instid(i)%p_chan_level(:,n), &  ! in
 289 !            k_cloud_flag, &  ! out
 290 !            kmin,     &  ! in
 291 !            kmax)        ! in
 292
 293             do k = 1, nchan
 294                ! remove channels above the model top
 295                if(iv%instid(i)%p_chan_level(k,n)==0) then
 296                   iv%instid(i)%tb_qc(k,n)  = qc_bad
 297                   if (iv%instid(i)%info%proc_domain(1,n)) &
 298                      nrej_eccloud = nrej_eccloud + 1
 299                end if
 300                if (k_cloud_flag(k) .eq. 1) then
 301                   iv%instid(i)%tb_qc(k,n)  = qc_bad
 302                   if (iv%instid(i)%info%proc_domain(1,n)) &
 303                      nrej_eccloud = nrej_eccloud + 1
 304                end if
 305             end do !chan
 306             deallocate ( k_cloud_flag )
 307          end if !ecmwf
 308
 309       end do ! end loop pixel
 310
 311       ! Do inter-processor communication to gather statistics.
 312
 313       if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
 314          do k = 1, nchan
 315             bias_global(k)   = wrf_dm_sum_real(bias_local(k))
 316             numrad_global(k) = wrf_dm_sum_integer(numrad_local(k))
 317             if (numrad_global(k) > 0) bias_global(k) = bias_global(k) / numrad_global(k)
 318          end do
 319       end if
 320
 321       do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
 322
 323          ! 1. Cloud detection scheme MMR in Auligné (2014).or. PF in Xu et al., (2016)
 324          !---------------------------------------------
 325          if ((use_clddet==1 .or. use_clddet==2) .and. (.not.use_satcv(2))) then
 326             iv%instid(i)%cloud_flag(:,n) = qc_good
 327
 328             if (rtm_option == rtm_option_rttov) then
 329 #ifdef RTTOV
 330                kte_surf   = iv%instid(i)%nlevels
 331                kts_100hPa = MAXLOC(coefs(i)%coef%ref_prfl_p(1:kte_surf), &
 332                             MASK = coefs(i)%coef%ref_prfl_p(1:kte_surf) < 100.0)
 333                do k=1,nchan
 334                   tstore = coefs(i)%coef%ff_bco(k) + coefs(i)%coef%ff_bcs(k) * &
 335                            (ob%instid(i)%tb(k,n) - bias_global(k))
 336                   iv%instid(i)%rad_obs(k,n) = coefs(i)%coef%planck1(k) / &
 337                            (EXP(coefs(i)%coef%planck2(k)/tstore) - 1.0)
 338                end do
 339 #endif
 340             elseif (rtm_option == rtm_option_crtm) then
 341                kte_surf   = kte
 342                kts_100hPa = MAXLOC(iv%instid(i)%pm(kts:kte,n), &
 343                             MASK = iv%instid(i)%pm(kts:kte,n) < 100.0)
 344
 345                do k = 1, nchan
 346                   CALL CRTM_Planck_Radiance(i,k,ob%instid(i)%tb(k,n) - bias_global(k), &
 347                                             iv%instid(i)%rad_obs(k,n))
 348                end do
 349
 350             end if
 351
 352             ndim = kte_surf - kts_100hPa(1) + 1
 353
 354             call da_cloud_detect(i,nchan,ndim,kts_100hPa(1),kte_surf,n,iv)
 355          end if
 356
 357          do k = 1, nchan
 358             if (iv%instid(i)%cloud_flag(k,n) == qc_bad) iv%instid(i)%tb_qc(k,n) = qc_bad
 359          end do
 360
 361          !  2. Check iuse from information file (channel selection)
 362          !-----------------------------------------------------------
 363          do k = 1, nchan
 364             if (satinfo(i)%iuse(k) .eq. -1) &
 365                iv%instid(i)%tb_qc(k,n)  = qc_bad
 366          end do
 367
 368          ! 3. Final QC decision
 369          !---------------------------------------------
 370          do k = 1, nchan
 371             if (iv%instid(i)%tb_qc(k,n) == qc_bad) then  ! bad obs
 372                iv%instid(i)%tb_error(k,n) = 500.0
 373                if (iv%instid(i)%info%proc_domain(1,n)) &
 374                   nrej(k) = nrej(k) + 1
 375             else                                         ! good obs
 376                if (iv%instid(i)%info%proc_domain(1,n)) &
 377                   ngood(k) = ngood(k) + 1
 378             end if
 379          end do ! chan
 380
 381       end do ! end loop pixel
 382
 383    ! Do inter-processor communication to gather statistics.
 384    call da_proc_sum_int  (num_proc_domain)
 385    call da_proc_sum_int  (nrej_landsurface )
 386    call da_proc_sum_int  (nrej_windowchlong)
 387    call da_proc_sum_int  (nrej_windowchshort)
 388    call da_proc_sum_int  (nrej_sst)
 389    call da_proc_sum_int  (nrej_clw  )
 390    call da_proc_sum_int  (nrej_eccloud )
 391    call da_proc_sum_int  (nrej_topo )
 392    call da_proc_sum_int  (nrej_limb)
 393    call da_proc_sum_ints (nrej_omb_abs(:))
 394    call da_proc_sum_ints (nrej_omb_std(:))
 395    call da_proc_sum_ints (nrej(:))
 396    call da_proc_sum_ints (ngood(:))
 397
 398    if (rootproc) then
 399       if (num_fgat_time > 1) then
 400          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 401       else
 402          write(filename,'(i2.2,a)') it, '_qcstat_'//trim(iv%instid(i)%rttovid_string)
 403       end if
 404
 405       call da_get_unit(fgat_rad_unit)
 406       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 407       if (ios /= 0) then
 408          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 409          call da_error(__FILE__,__LINE__,message(1:1))
 410       end if
 411
 412       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 413       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 414       write(fgat_rad_unit,'(a20,i7)') ' nrej_landsurface  = ', nrej_landsurface
 415       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchlong = ', nrej_windowchlong
 416       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchshort = ', nrej_windowchshort
 417       write(fgat_rad_unit,'(a20,i7)') ' nrej_sst = ', nrej_sst
 418       !      write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 419       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 420       write(fgat_rad_unit,'(a20,i7)') ' nrej_eccloud        = ', nrej_eccloud
 421       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
 422       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
 423       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 424       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 425       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 426       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 427       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 428       write(fgat_rad_unit,'(10i7)')     nrej(:)
 429       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 430       write(fgat_rad_unit,'(10i7)')     ngood(:)
 431
 432       close(fgat_rad_unit)
 433       call da_free_unit(fgat_rad_unit)
 434    end if
 435
 436    if (trace_use_dull) call da_trace_exit("da_qc_airs")
 437
 438 end subroutine da_qc_airs
 439