Fortran 90 SUBROUTINE ccm2_hydro_mat
!>============================================================================<
! FILE: /fs/cgd/home0/davestep/f90_routines/SUB_ccm2_hydro_mat.f90
! AUTHOR: David Stepaniak, NCAR/CGD/CAS
! DATE INITIATED: 21 August 1998
! LAST MODIFIED: Tue Aug 25 12:03:56 MDT 1998
!
!>----------------------------------------------------------------------------<
!
! 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.
!
!>============================================================================<
SUBROUTINE ccm2_hydro_mat(ucK, pi, p, H)
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 varibles:
INTEGER :: k ! Lower case k representing the
! row index of H.
INTEGER :: l ! Lower case l representing the
! column index of H.
!>----------------------------------------------------------------------------<
! Check that p increases monotonically (i.e. midpoint pressure values arranged
! from top to bottom), and that p(K) < pi. Both conditions must be satisfied.
! (This part may be removed to avoid error-checking every time SUBROUTINE
! ccm2_hydro_mat is invoked.)
DO k = 2, uck - 1
IF ( .NOT. ( ( p(k-1) < p(k) ) .AND. ( p(k) < p(k+1) ) ) ) THEN
WRITE (*,*) "FATAL ERROR CONDITION: "
WRITE (*,*) "Midpoint pressure p must "
WRITE (*,*) "increase monotonically. "
WRITE (*,*) "Execution halted in "
WRITE (*,*) "SUBROUTINE ccm2_hydro_mat."
STOP
END IF
END DO
IF ( .NOT. ( p(ucK) < pi ) ) THEN
WRITE (*,*) "FATAL ERROR CONDITION: "
WRITE (*,*) "Surface pressure pi must be "
WRITE (*,*) "greater then p(K). "
WRITE (*,*) "Execution halted in "
WRITE (*,*) "SUBROUTINE ccm2_hydro_mat. "
STOP
END IF
!>----------------------------------------------------------------------------<
! 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
!>----------------------------------------------------------------------------<