Fortran 90 SUBROUTINE CCM2_Z2_sigma
!>============================================================================<
! FILE: SUB_CCM2_Z2_sigma.f90
! AUTHOR: David Stepaniak, NCAR/CGD/CAS
! DATE INITIATED: 9 September 1998
! LAST MODIFIED: Mon Sep 14 13:41:39 MDT 1998
!
! DESCRIPTION: Computes Z2 on sigma levels. Z2 (a CCM Processor Code-Defined
! Derived Field) is the geopotential height based on the CCM2
! hydrostatic formulation. For a fixed time, latitude, and
! longitude, the geopotential height at model mid-levels may
! be written
!
! Z2 = PHIS/g + (R/g) * H * TV
!
! where PHIS is the surface geopotential, g the gravitational
! acceleration, R the dry gas constant, H the CCM2 hydrostatic
! matrix (which is a function of surface pressure and pressure
! at each model mid-level), and TV the virtual temperature.
! Subroutine CCM2_Z2_sigma contains ccm2_hydro_mat, a subrou-
! tine which computes H.
!
! The ordering of the vertical profiles of T, Q, and sigma
! must be TOP to BOTTOM. The vertical profile of Z2 on exit is
! also TOP to BOTTOM.
!
!>============================================================================<
SUBROUTINE CCM2_Z2_sigma( nlvl, T, Q, PS, PHIS, Z2, g, Rd, Qmf, sigma )
IMPLICIT NONE
INTEGER, INTENT(IN) :: nlvl
! Number of levels.
REAL, DIMENSION(1:nlvl), INTENT(IN) :: T
! Temperature, Kelvin.
REAL, DIMENSION(1:nlvl), INTENT(IN) :: Q
! Specific humidity,
! kg/kg.
REAL, INTENT(IN) :: PS
! Surface pressure, Pa.
REAL, INTENT(IN) :: PHIS
! Surface geopotential,
! m^2/s^2.
REAL, DIMENSION(1:nlvl), INTENT(OUT) :: Z2
! Geopotential height
! based on CCM2 hydro-
! static formulation.
REAL, INTENT(IN) :: g
! Gravitational ac-
! celeration, m/s^2.
REAL, INTENT(IN) :: Rd
! Dry gas constant,
! J / (kg K)
REAL, INTENT(IN) :: Qmf
! Multiplicative factor
! of Q in formula for
! virtual temperature,
! i.e.:
! TV = T(1 + Qmf*Q)
REAL, DIMENSION(1:nlvl), INTENT(IN) :: sigma
! Sigma levels top to
! bottom
!>----------------------------------------------------------------------------<
! Local variables:
REAL :: Rdog
! Dry gas constant di-
! vided by gravitational
! acceleration.
REAL, DIMENSION(:), ALLOCATABLE :: TV
! Virtual temperature,
! K.
REAL :: ZS
! Geopotential height
! at surface, m.
REAL, DIMENSION(:,:), ALLOCATABLE :: H
! CCM2 hydrostatic
! matrix.
REAL, DIMENSION(:), ALLOCATABLE :: pmid
! Mid-level pressures
! computed from sigma,
! Pa.
!>----------------------------------------------------------------------------<
! Assign dry gas constant divided by gravitational acceleration:
Rdog = Rd / g
!>----------------------------------------------------------------------------<
! Allocate and compute virtual temperature and geopotential height of surface:
ALLOCATE( TV(1:nlvl) )
TV = T * ( 1. + Qmf*Q )
ZS = PHIS / g
!>----------------------------------------------------------------------------<
! Compute mid-level pressures from sigma and surface pressure, compute CCM2
! hydrostatic matrix H, and compute Z2. The subroutine ccm2_hydro_mat is called
! to compute H.
ALLOCATE( pmid(1:nlvl ) )
ALLOCATE( H(1:nlvl,1:nlvl) )
pmid = sigma * PS
CALL ccm2_hydro_mat(nlvl, PS, pmid, H)
Z2(:) = ZS + Rdog * MATMUL( H, TV(:) )
!>----------------------------------------------------------------------------<
! Deallocate virtual temperature, surface geopotential height, mid-level
! pressures, and CCM2 hydrostatic matrix.
DEALLOCATE( TV )
DEALLOCATE( pmid )
DEALLOCATE( H )
!>----------------------------------------------------------------------------<
CONTAINS
SUBROUTINE ccm2_hydro_mat(ucK, pi, p, H)
! DESCRIPTION: Denoting the number of midpoint model levels by K,
! SUBROUTINE ccm2_hydro_mat returns the [K x K] CCM2
! hydostatic matrix H defined by Equation (3.a.109)
! in NCAR Technical Note 382 (Hack et al., 1993). The
! notation in this code follows the notation found in
! the Technical Note.
!
! Let k be the row index of H, l the column index of
! H, p_k the midpoint pressure where k = 1,2...,K arranged
! from top to bottom, and pi the surface pressure.
! (The midpoint pressure values are derived, for example,
! from purely sigma or hybrid coordinate midpoint
! values.) Then, from Hack et al., p. 27, H(k,l) is
! defined as:
!
! / 0, l < k
! |
! | .5 ln( p_(k+1)/p_k ), l = k, k < K
! |
! H(k,l) = < .5 ln( p_(l+1)/p_(l-1) ), l > k, k < K
! |
! | .5 ln( pi^2/(p_(K-1)p_K) ), l = K, k < K
! |
! \ ln( pi/p_K ) l = K, k = K
!
!
!
! A prime example of the utility of H is in the calculation
! of Z2 (a code-defined derived field) from the CCM processor.
! Z2 is the geopotential height based on the CCM2 hydrostatic
! formulation. For a given vertical profile of virtual tempera-
! ture, say TV_k, Z2 in Fortran 90 is given by the matrix
! equation
!
! Z2 = PHIS/g + (R/g) * H * TV
!
! where PHIS is the surface geopotential and R is the gas
! constant for dry air. In this case the column vector Z2
! represents the geopoential height arranged from top to
! bottom.
!
! REFERENCE: Hack, J.J., B.A. Boville, B.P. Briegleb, J.T. Kiehl, P.J.
! Rasch, D.L. Williamson, 1993: Description of the NCAR
! Community Climate Model (CCM2). NCAR Technical Note
! NCAR/TN-382+STR, 108 pp.
!
! NOTE: K and ucK are used interchangeably in comments.
IMPLICIT NONE
INTEGER, INTENT(IN) :: ucK ! Number of midpoint model levels
! (uc refers to upper case).
REAL, INTENT(IN) :: pi ! Surface pressure, Pa or mb.
REAL, DIMENSION(1:ucK), INTENT(IN) :: p ! Midpoint pressure values,
! ARRANGED FROM TOP TO BOTTOM.
! Units Pa or mb, but same as
! units used for pi. The midpoint
! pressure values are derived,
! for example, from purely sigma
! or hybrid coordinate midpoint
! values.
REAL, DIMENSION(1:ucK,1:ucK), INTENT(OUT) :: H
! CCM2 hydrostatic matrix defined
! by Equation (3.a.109) in Hack
! et al. (1993). H is indexed
! as H(k,l) where k and l are
! defined below.
! Local variables:
INTEGER :: k ! Lower case k representing the
! row index of H.
INTEGER :: l ! Lower case l representing the
! column index of H.
! Initialize all elements of H as 0.:
H = 0.
! Compute all diagonal elements of H except H(K,K):
DO k = 1, ucK - 1
H(k,k) = .5 * LOG( p(k+1) / p(k) )
END DO
! Compute all off-diagonal elements of H, except for last two rows, and last
! column:
DO k = 1, ucK - 2
DO l = k + 1, ucK - 1
H(k,l) = .5 * LOG( p(l+1) / p(l-1) )
END DO
END DO
! Compute last column of H, except k = K:
DO k = 1, ucK - 1
H(k,ucK) = .5 * LOG( pi*pi / (p(ucK-1)*p(ucK)) )
END DO
! Compute H(K,K):
H(ucK,ucK) = LOG( pi / p(ucK) )
END SUBROUTINE ccm2_hydro_mat
!>----------------------------------------------------------------------------<
END SUBROUTINE CCM2_Z2_sigma
!>----------------------------------------------------------------------------<