subroutine geopotential  ( pcols  , pver   ,                   & 1
           piln   , pmln   , pint   , pmid   , pdel   , rpdel  , &
           t      , q      , rair   , gravit , zvir   ,          &
           zi     , zm     )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute the geopotential height (above the surface) at the midpoints and 
! interfaces using the input temperatures and pressures.
! Author: B.Boville, Feb 2001 from earlier code by Boville and S.J. Lin
!
!-----------------------------------------------------------------------

    implicit none

!------------------------------Arguments--------------------------------
!
! Input arguments
!
    integer, intent(in) :: pcols                 ! Number of longitudes
    integer, intent(in) :: pver                  ! Number of vertical layers

    real,    intent(in) :: piln (pcols,pver+1)   ! Log interface pressures
    real,    intent(in) :: pmln (pcols,pver)     ! Log midpoint pressures
    real,    intent(in) :: pint (pcols,pver+1)   ! Interface pressures
    real,    intent(in) :: pmid (pcols,pver)    ! Midpoint pressures
    real,    intent(in) :: pdel (pcols,pver)    ! layer thickness
    real,    intent(in) :: rpdel(pcols,pver)    ! inverse of layer thickness
    real,    intent(in) :: t    (pcols,pver)    ! temperature
    real,    intent(in) :: q    (pcols,pver)    ! specific humidity
    real,    intent(in) :: rair                 ! Gas constant for dry air
    real,    intent(in) :: gravit               ! Acceleration of gravity
    real,    intent(in) :: zvir                 ! rh2o/rair - 1

! Output arguments

    real,    intent(out) :: zi(pcols,pver+1)    ! Height above surface at interfaces
    real,    intent(out) :: zm(pcols,pver)      ! Geopotential height at mid level
!
!---------------------------Local variables-----------------------------
!
    logical  :: fvdyn              ! finite volume dynamics
    integer  :: i,k                ! Lon, level indices
    real     :: hkk(pcols)         ! diagonal element of hydrostatic matrix
    real     :: hkl(pcols)         ! off-diagonal element
    real     :: rog                ! Rair / gravit
    real     :: tv                 ! virtual temperature
    real     :: tvfac              ! Tv/T
    integer  :: ncol               ! Number of longitudes
!
!-----------------------------------------------------------------------
!
    ncol = pcols
    rog = rair/gravit

! Set dynamics flag

    fvdyn = .true.

! The surface height is zero by definition.

    do i = 1,ncol
       zi(i,pver+1) = 0.0
    end do

! Compute zi, zm from bottom up. 
! Note, zi(i,k) is the interface above zm(i,k)

    do k = pver, 1, -1

! First set hydrostatic elements consistent with dynamics

       if (fvdyn) then
          do i = 1,ncol
             hkl(i) = piln(i,k+1) - piln(i,k)
             hkk(i) = 1. - pint(i,k) * hkl(i) * rpdel(i,k)
          end do
       else
          do i = 1,ncol
             hkl(i) = pdel(i,k) / pmid(i,k)
             hkk(i) = 0.5 * hkl(i)
          end do
       end if

! Now compute tv, zm, zi

       do i = 1,ncol
          tvfac   = 1. + zvir * q(i,k)
          tv      = t(i,k) * tvfac

          zm(i,k) = zi(i,k+1) + rog * tv * hkk(i)
          zi(i,k) = zi(i,k+1) + rog * tv * hkl(i)
       end do
    end do

    return
  end subroutine geopotential