share/module_compute_geop.F

   1 MODULE module_compute_geop
   2
   3 CONTAINS
   4   SUBROUTINE compute_500mb_height  ( ph, phb, p, pb,                  &
   5                                    height, track_level,             &
   6                                    ids, ide, jds, jde, kds, kde,    &
   7                                    ims, ime, jms, jme, kms, kme,    &
   8                                    its, ite, jts, jte, kts, kte    )
   9
  10    IMPLICIT NONE
  11
  12
  13    !  Input data.
  14
  15    INTEGER ,       INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
  16                                     ims, ime, jms, jme, kms, kme, &
  17                                     its, ite, jts, jte, kts, kte
  18
  19    REAL , DIMENSION(  ims:ime , kms:kme, jms:jme ) ,                      &
  20                                                INTENT(IN   ) ::           &
  21                                                                  ph,      &
  22                                                                  phb,     &
  23                                                                  pb,      &
  24                                                                  p
  25
  26    REAL , DIMENSION( ims:ime , jms:jme ) ,    INTENT(  OUT) :: height
  27
  28    INTEGER , INTENT(IN) :: track_level
  29
  30 !  local variables
  31
  32    integer :: i,j,k
  33    real, dimension(kms:kme) :: pressure,geopotential
  34    real :: interp_output
  35    real :: track_level_p
  36
  37 !  slow version of code, we'll call interp routine for each column
  38
  39    track_level_p = float(track_level)
  40
  41    do j = jts, min(jde-1,jte)
  42    do i = its, min(ide-1,ite)
  43
  44       do k=kds,kde-1
  45         pressure(k) = p(i,k,j) + pb(i,k,j)
  46         geopotential(k) = 0.5*( ph(i,k  ,j)+phb(i,k  ,j)  &
  47                                +ph(i,k+1,j)+phb(i,k+1,j) )
  48       enddo
  49
  50 !     call interp_p( geopotential, pressure, 70000., interp_output,  &
  51       call interp_p( geopotential, pressure, track_level_p, interp_output,  &
  52                      kds,kde-1,kms,kme, i,j )
  53
  54       height(i,j) = interp_output/9.81  !  500 mb height in meters
  55
  56    enddo
  57    enddo
  58
  59    end subroutine compute_500mb_height
  60
  61 !--------
  62
  63   subroutine interp_p(a,p,p_loc,a_interp,ks,ke,kms,kme,i,j)
  64   implicit none
  65
  66   integer, intent(in) :: ks,ke,kms,kme,i,j
  67   real, dimension(kms:kme), intent(in) :: a,p
  68   real, intent(in)  :: p_loc
  69   real, intent(out) :: a_interp
  70
  71 !---  local variables
  72
  73   integer :: kp, pk, k
  74   real    :: wght1, wght2, dp, pressure
  75   character*256 mess
  76
  77 !cys_change: set high value at below-ground point
  78     if (p(ks).lt.p_loc) then
  79        a_interp=9.81e5
  80     else
  81
  82     kp = ks+1
  83     pk = p(kp)
  84     pressure = p_loc
  85     do while( pk .gt. pressure )
  86
  87       kp = kp+1
  88
  89       if(kp .gt. ke) then
  90         write(mess,*) ' interp too high: pressure, p(ke), i, j = ',pressure,p(ke),i,j
  91         call wrf_message ( mess )
  92         write(mess,*)'p: ',p
  93         call wrf_error_fatal( mess )
  94       end if
  95
  96       pk = p(kp)
  97
  98     enddo
  99
 100     dp = p(kp-1) - p(kp)
 101     wght2 = (p(kp-1)-pressure)/dp
 102     wght1 = 1.-wght2
 103
 104     a_interp = wght1*a(kp-1) + wght2*a(kp)
 105
 106     endif   !cys_change
 107
 108     end subroutine interp_p
 109
 110 END MODULE module_compute_geop