var/da/da_radiance/da_qc_seviri.inc

   1 subroutine da_qc_seviri (it, i, nchan, ob, iv)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: perform quality control for seviri 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      :: si
  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_seviri")
  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
  43       do n= iv%instid(i)%info%n1,iv%instid(i)%info%n2
  44
  45          if (iv%instid(i)%info%proc_domain(1,n)) &
  46                num_proc_domain = num_proc_domain + 1
  47
  48          !  0.0  initialise QC by flags assuming good obs
  49          !---------------------------------------------
  50          iv%instid(i)%tb_qc(:,n) = qc_good
  51
  52          ! observation errors
  53          ! assigned to satinfo%error_std instead of satinfo%error
  54          ! this is to make the function (read_biascoef) in da_radiance_init.inc work for seviri
  55          do k = 1, nchan
  56             if (use_error_factor_rad) then
  57                iv%instid(i)%tb_error(k,n) = &
  58                    satinfo(i)%error_std(k)*satinfo(i)%error_factor(k)
  59             else
  60                iv%instid(i)%tb_error(k,n) = satinfo(i)%error_std(k)
  61             end if
  62         end do
  63
  64          !  1.0 check iuse
  65          !-----------------------------------------------------------
  66          do k = 1, nchan
  67             if (satinfo(i)%iuse(k) .eq. -1) &
  68                iv%instid(i)%tb_qc(k,n)  = qc_bad
  69          end do
  70
  71          !  2.0 check innovation
  72          !-----------------------------------------------------------
  73          do k = 1, nchan
  74             if ( iv%instid(i)%tb_qc(k,n) .eq. qc_good ) then
  75
  76                ! absolute departure check
  77                if (abs(iv%instid(i)%tb_inv(k,n)) > 15.0) then
  78                   iv%instid(i)%tb_qc(k,n)  = qc_bad
  79                   if (iv%instid(i)%info%proc_domain(1,n)) &
  80                   nrej_omb_abs(k) = nrej_omb_abs(k) + 1
  81                end if
  82
  83                ! relative departure check
  84                if (abs(iv%instid(i)%tb_inv(k,n)) > 3.0*iv%instid(i)%tb_error(k,n)) then
  85                   iv%instid(i)%tb_qc(k,n)  = qc_bad
  86                   !iv%instid(i)%tb_error(k,n) = 500.0
  87                   if (iv%instid(i)%info%proc_domain(1,n)) &
  88                   nrej_omb_std(k) = nrej_omb_std(k) + 1
  89                end if
  90
  91                ! final QC decsion
  92                if (iv%instid(i)%tb_qc(k,n) == qc_bad) then
  93                   iv%instid(i)%tb_error(k,n) = 500.0
  94                   if (iv%instid(i)%info%proc_domain(1,n)) &
  95                      nrej(k) = nrej(k) + 1
  96                else
  97                   if (iv%instid(i)%info%proc_domain(1,n)) &
  98                      ngood(k) = ngood(k) + 1
  99                end if
 100
 101             end if   ! qc_good
 102          end do      ! nchan
 103       end do ! end loop pixel
 104
 105    ! Do inter-processor communication to gather statistics.
 106    call da_proc_sum_int (num_proc_domain)
 107    call da_proc_sum_int (nrej_mixsurface)
 108    call da_proc_sum_int (nrej_windowchanl)
 109    call da_proc_sum_int (nrej_si )
 110    call da_proc_sum_int (nrej_clw)
 111    call da_proc_sum_int (nrej_topo)
 112    call da_proc_sum_int (nrej_limb)
 113    call da_proc_sum_ints (nrej_omb_abs(:))
 114    call da_proc_sum_ints (nrej_omb_std(:))
 115    call da_proc_sum_ints (nrej(:))
 116    call da_proc_sum_ints (ngood(:))
 117
 118    if (rootproc) then
 119       if (num_fgat_time > 1) then
 120          write(filename,'(i2.2,a,i2.2)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)//'_',iv%time
 121       else
 122          write(filename,'(i2.2,a)') it,'_qcstat_'//trim(iv%instid(i)%rttovid_string)
 123       end if
 124
 125       call da_get_unit(fgat_rad_unit)
 126       open(fgat_rad_unit,file=trim(filename),form='formatted',iostat=ios)
 127       if (ios /= 0) then
 128          write(unit=message(1),fmt='(A,A)') 'error opening the output file ', filename
 129          call da_error(__FILE__,__LINE__,message(1:1))
 130       end if
 131
 132       write(fgat_rad_unit, fmt='(/a/)') ' Quality Control Statistics for '//iv%instid(i)%rttovid_string
 133       write(fgat_rad_unit,'(a20,i7)') ' num_proc_domain  = ', num_proc_domain
 134       write(fgat_rad_unit,'(a20,i7)') ' nrej_mixsurface  = ', nrej_mixsurface
 135       write(fgat_rad_unit,'(a20,i7)') ' nrej_windowchanl = ', nrej_windowchanl
 136       write(fgat_rad_unit,'(a20,i7)') ' nrej_si          = ', nrej_si
 137       write(fgat_rad_unit,'(a20,i7)') ' nrej_clw         = ', nrej_clw
 138       write(fgat_rad_unit,'(a20,i7)') ' nrej_topo        = ', nrej_topo
 139       write(fgat_rad_unit,'(a20,i7)') ' nrej_limb        = ', nrej_limb
 140       write(fgat_rad_unit,'(a20)')    ' nrej_omb_abs(:)  = '
 141       write(fgat_rad_unit,'(10i7)')     nrej_omb_abs(:)
 142       write(fgat_rad_unit,'(a20)')    ' nrej_omb_std(:)  = '
 143       write(fgat_rad_unit,'(10i7)')     nrej_omb_std(:)
 144       write(fgat_rad_unit,'(a20)')    ' nrej(:)          = '
 145       write(fgat_rad_unit,'(10i7)')     nrej(:)
 146       write(fgat_rad_unit,'(a20)')    ' ngood(:)         = '
 147       write(fgat_rad_unit,'(10i7)')     ngood(:)
 148
 149       close(fgat_rad_unit)
 150       call da_free_unit(fgat_rad_unit)
 151    end if
 152    if (trace_use) call da_trace_exit("da_qc_seviri")
 153
 154 end subroutine da_qc_seviri
 155
 156