var/convertor/average_be/da_interpolate_regcoeff.F

   1 subroutine da_interpolate_regcoeff( kzs, kz, meanp_stats, meanp_xb, &
   2                                     stats_regcoeff1, stats_regcoeff2, stats_regcoeff3, &
   3                                     xb_regcoeff1, xb_regcoeff2, xb_regcoeff3 )
   4
   5 !------------------------------------------------------------------------------
   6 !  PURPOSE: Interpolate statistical regression coefficient to new domain.
   7 !
   8 !------------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    integer, intent(in)                       :: kzs               ! Number of levels in stats.
  13    integer, intent(in)                       :: kz                ! Number of levels in xb.
  14    real, dimension(kzs),         intent(in)  :: meanp_stats       ! Mean pressure on stats levs.
  15    real, dimension(kz),          intent(in)  :: meanp_xb          ! Mean pressure on xb levs.
  16    real, dimension(kzs),         intent(in)  :: stats_regcoeff1, stats_regcoeff2
  17    real, dimension(kzs,kzs),     intent(in)  :: stats_regcoeff3
  18    real, dimension(kz),          intent(out) :: xb_regcoeff1, xb_regcoeff2
  19    real, dimension(kz,kz),       intent(out) :: xb_regcoeff3
  20
  21    integer                :: k1, k2, k, ks            ! Loop counters.
  22    integer                :: ktrap_min, ktrap_max     ! Max/min limits of xb rows.
  23    integer                :: k1s, k2s
  24    integer                :: k_above(1:kz)
  25    real                   :: weight(1:kz)
  26
  27 !---------------------------------------------------------------------------
  28 !  [1.0] Compare stats/xb levels and compute interpolation weights:
  29 !---------------------------------------------------------------------------
  30 !
  31    k_above(1:kz) = 0
  32    weight(1:kz) = 0.0
  33
  34    do k = 1, kz
  35       if ( meanp_xb(k) <= meanp_stats(1) ) then
  36       ktrap_min = k
  37       go to 10
  38       end if
  39    end do
  40    print*,' problem in determining ktrap_min'
  41    stop
  42 10  continue
  43
  44    do k = kz, 1, -1
  45       if ( meanp_xb(k) >= meanp_stats(kzs) ) then
  46       ktrap_max = k
  47       go to 20
  48       end if
  49    end do
  50    print*,' problem in determining ktrap_max'
  51    stop
  52 20  continue
  53
  54    do k = ktrap_min, ktrap_max
  55       do ks = 1, kzs-1
  56          if ( meanp_xb(k) > meanp_stats(ks+1) .AND. &
  57               meanp_xb(k) <= meanp_stats(ks) ) then
  58
  59             weight(k) = ( meanp_xb(k) - meanp_stats(ks+1) ) / &
  60                         ( meanp_stats(ks) - meanp_stats(ks+1) )
  61             k_above(k) = ks+1
  62             exit
  63          end if
  64       end do
  65    end do
  66
  67    do k = 1, ktrap_min - 1
  68       k_above(k) = 2
  69       weight(k) = 1.0
  70    enddo
  71    do k = ktrap_max + 1, kz
  72       k_above(k) = kzs - 1
  73       weight(k) = 0.0
  74    enddo
  75
  76 !---------------------------------------------------------------------------
  77 !  [3.0] Interpolate regression coefficient from stats to xb levels:
  78 !---------------------------------------------------------------------------
  79       do k1 = 1, kz
  80          k1s = k_above(k1)
  81 !print*,k1,' interpolating between ',k1s-1,k1s,' weights ',weight(k1),' stats ',stats_regcoeff1(k1s-1),stats_regcoeff1(k1s)
  82 !            xb_regcoeff1(kz-k1+1) = &
  83             xb_regcoeff1(k1) = &
  84                                (1.0-weight(k1)) * stats_regcoeff1(k1s) + &
  85                                     weight(k1)  * stats_regcoeff1(k1s-1)
  86 !            xb_regcoeff2(kz-k1+1) = &
  87             xb_regcoeff2(k1) = &
  88                                (1.0-weight(k1)) * stats_regcoeff2(k1s) + &
  89                                     weight(k1)  * stats_regcoeff2(k1s-1)
  90
  91          do k2 = 1, kz
  92            k2s = k_above(k2)
  93
  94 !           xb_regcoeff3(kz-k1+1,kz-k2+1) = &
  95            xb_regcoeff3(k1,k2) = &
  96                               (1.0-weight(k1)) * (1.0-weight(k2)) * stats_regcoeff3(k1s,k2s) + &
  97                                    weight(k1)  * (1.0-weight(k2)) * stats_regcoeff3(k1s-1,k2s  ) + &
  98                               (1.0-weight(k1)) *      weight(k2)  * stats_regcoeff3(k1s  ,k2s-1) + &
  99                                    weight(k1)  *      weight(k2)  * stats_regcoeff3(k1s-1,k2s-1)
 100
 101          end do
 102       end do
 103
 104 end subroutine da_interpolate_regcoeff
 105