var/da/da_varbc/da_varbc_coldstart.inc

   1   subroutine da_varbc_coldstart (iv)
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE: [1]: If a cold-start is needed, calculate mode of histogram
   5    !                of (uncorrected) innovations for each channel
   6    !
   7    !           [2]: Calculate statistics (mean/std) for cold-started predictors
   8    !
   9    !           [3]: Normalize predictors
  10    !
  11    !  Called from da_get_innov_vector_radiance
  12    !
  13    !  HISTORY: 10/26/2007 - Creation                     Tom Auligne
  14    !           11/03/2008 - Bug-fix: exclude HALO data for
  15    !                        computation of predictor statistics    Tom/Hui-Chuan
  16    !---------------------------------------------------------------------------
  17
  18    implicit none
  19
  20    type (iv_type), intent (inout)   :: iv             ! Innovation
  21
  22    integer                          :: inst, k, n, i
  23    integer                          :: npredmax, npred, num_rad, num_rad_domain, num_rad_tot
  24
  25    integer, parameter               :: nbins   = 200                       ! Number of Hist bins.
  26    real,    parameter               :: maxhist = 10.                       ! Maximum bin value.
  27    real,    parameter               :: zbin    = 2 * maxhist / real(nbins) ! Hist bin width.
  28    integer                          :: ibin                                ! Hist bin number.
  29    real                             :: modetmp(1), mode
  30    integer, allocatable             :: hist(:), hist_tot(:)
  31
  32    real, allocatable                :: mean(:), rms(:), mean_tot(:), rms_tot(:)
  33    integer, allocatable             :: ipred(:)
  34    logical                          :: global_cs
  35
  36    if ( iv%num_inst < 1 ) RETURN
  37
  38    if (trace_use) call da_trace_entry("da_varbc_coldstart")
  39
  40    do inst = 1, iv%num_inst                                            !! loop for sensors
  41
  42       npredmax       = iv%instid(inst)%varbc_info%npredmax
  43       num_rad        = iv%instid(inst)%num_rad
  44       if (num_rad > 0) then
  45          num_rad_domain = COUNT(iv%instid(inst)%info%proc_domain(1,1:num_rad)) !do not count HALO
  46       else
  47          num_rad_domain = 0
  48       end if
  49       num_rad_tot    = wrf_dm_sum_integer(num_rad_domain)
  50
  51       if (npredmax <= 0) cycle                                         !! VarBC instr only
  52
  53       allocate( ipred(npredmax) )
  54
  55     !---------------------------------------------------------------------------
  56     ! [1]: Calculate mode of histogram of (uncorrected) innovations
  57     !---------------------------------------------------------------------------
  58       do k = 1, iv%instid(inst)%nchan                                  !! loop for channels
  59          npred          = iv%instid(inst)%varbc(k)%npred
  60          ipred(1:npred) = iv%instid(inst)%varbc(k)%ipred(1:npred)
  61          if (npred <= 0) cycle                                         !! VarBC channels only
  62          if (ALL(iv%instid(inst)%varbc(k)%pred_use /= 0)) cycle        !! Coldstart channels only
  63
  64          where (iv%instid(inst)%varbc(k)%pred_use(ipred(1:npred)) == 0) &
  65                 iv%instid(inst)%varbc(k)%param(1:npred) = 0.0
  66
  67          if (iv%instid(inst)%varbc(k)%pred_use(1) == 0) then
  68             Allocate ( hist(nbins),  hist_tot(nbins))
  69             hist(:) = 0
  70             mode    = 0.0
  71
  72            ! Accumulate statistics for histogram
  73            ! -----------------------------------
  74             do n = 1, num_rad      !! loop for pixel
  75                if (iv%instid(inst)%info%proc_domain(1,n)) then ! do not count HALO
  76                   ibin = NINT( (iv%instid(inst)%tb_inv(k,n)+maxhist)/zbin )
  77                   if ((ibin>0).AND.(ibin<=nbins)) &
  78                      hist(ibin) = hist(ibin) + 1
  79                end if
  80             end do             ! end loop for pixels
  81
  82            ! Do inter-processor communication to gather statistics
  83            ! ------------------------------------------------------
  84             do ibin = 1, nbins
  85                hist_tot(ibin) = wrf_dm_sum_integer(hist(ibin))
  86             end do
  87
  88            ! Determine mode of Histogram
  89            !----------------------------
  90             if ( SUM(hist_tot(:)) > 0 ) then
  91               modetmp(1:1) = MAXLOC(hist_tot(:))*zbin - maxhist
  92               mode = modetmp(1)
  93             end if
  94
  95            ! Use mode to initialize VarBC
  96            !-----------------------------
  97             if (iv%instid(inst)%varbc(k)%ipred(1) == 1) &
  98                 iv%instid(inst)%varbc(k)%param(1) = mode
  99
 100             Deallocate ( hist, hist_tot )
 101             if ( satinfo(inst)%iuse(k) == 1 ) &
 102                write(unit=stdout,fmt='(A,A,I5,A,F5.2)') 'VARBC: Cold-starting ', &
 103                   trim(adjustl(iv%instid(inst)%rttovid_string)),iv%instid(inst)%ichan(k),&
 104                   ' --> ',mode
 105          end if
 106       end do                                                              !  end loop for channels
 107
 108     !---------------------------------------------------------------------------
 109     !  [2]: Calculate statistics for cold-started predictors
 110     !---------------------------------------------------------------------------
 111       global_cs = .true.
 112       do k = 1, iv%instid(inst)%nchan
 113          if (iv%instid(inst)%varbc(k)%npred <= 0) cycle                   !! VarBC channels only
 114          if (ANY(iv%instid(inst)%varbc(k)%pred_use > 0)) global_cs = .false.
 115       end do
 116
 117       if (global_cs) then                                                 !! Instrument coldstart only
 118
 119          allocate (mean(npredmax), rms(npredmax), mean_tot(npredmax), rms_tot(npredmax))
 120
 121         ! Accumulate statistics for predictor mean/std
 122         ! ---------------------------------------------
 123          if (num_rad > 0) then
 124             do i = 1, npredmax
 125                mean(i) = SUM( iv%instid(inst)%varbc_info%pred(i,1:num_rad),    &
 126                          MASK=iv%instid(inst)%info%proc_domain(1,1:num_rad))   ! do not count HALO
 127                rms(i)  = SUM( iv%instid(inst)%varbc_info%pred(i,1:num_rad)**2, &
 128                          MASK=iv%instid(inst)%info%proc_domain(1,1:num_rad))   ! do not count HALO
 129             end do
 130          else
 131             mean = 0.0
 132             rms  = 0.0
 133          end if
 134
 135         ! Do inter-processor communication to gather statistics
 136         ! ------------------------------------------------------
 137          call wrf_dm_sum_reals(mean, mean_tot)
 138          call wrf_dm_sum_reals(rms,  rms_tot)
 139          if (num_rad_tot >= varbc_nobsmin) then
 140             mean_tot = mean_tot / num_rad_tot
 141             rms_tot  = rms_tot  / num_rad_tot
 142          else
 143             mean_tot = 0.0
 144             rms_tot  = 1.0
 145          end if
 146
 147         ! Store statistics
 148         !------------------
 149          iv%instid(inst)%varbc_info%pred_mean = mean_tot
 150          iv%instid(inst)%varbc_info%pred_std  = sqrt(rms_tot - mean_tot**2)
 151
 152          deallocate(mean, rms, mean_tot, rms_tot)
 153
 154       end if
 155       deallocate(ipred)
 156
 157     !---------------------------------------------------------------------------
 158     !  [3]: Normalize predictors
 159     !---------------------------------------------------------------------------
 160       do i = 1,  npredmax
 161          if ( iv%instid(inst)%varbc_info%pred_std(i) <= 0.0 ) cycle
 162          do n = 1, num_rad
 163             iv%instid(inst)%varbc_info%pred(i,n) = &
 164           ( iv%instid(inst)%varbc_info%pred(i,n) - &
 165             iv%instid(inst)%varbc_info%pred_mean(i) ) / &
 166             iv%instid(inst)%varbc_info%pred_std(i)
 167          end do
 168       end do
 169
 170    end do                           !  end loop for sensor
 171
 172    if (trace_use) call da_trace_exit("da_varbc_coldstart")
 173
 174  end subroutine da_varbc_coldstart