Fortran 90 FUNCTION V_INT_PRESS
!>============================================================================<
! FILE: FUNC_V_INT_PRESS.f90
! DATE: Wed Oct 29 13:24:08 MST 1997
! MODIFIED: Tue Jan 19 11:19:27 MST 1999
!
! NOTE: Program or subprogram calling this function MUST have an explicit
! interface block as two assumed-shape arrays are utilized. If this
! function is called from a Fortran 77 program or subprogram then
! the number of archived pressure levels must be passed explicitly.
!
! DESCRIPTION:
!
! Given data on J archived pressure levels, this function computes the
! vertical integral for a fixed time, field, and grid point. The pressure
! at the bottom and top of the atmospheric column must be specified so as
! to provide fixed integration limits, and the surface pressure of a real
! or envelope orography must be specified as well. The formulation follows
! Trenberth and Solomon (1994) who define the following pressure levels
! and indexing:
!
! ============================== p_(j=2J) = p_t Prescribed p at top of column.
!
! ------------------- J p_(j=2J-1) Archived p level.
!
! ======================== p_(j=2J-2) Intermediate half-level.
! . . .
! . . .
! . . .
! ======================== p_(j=4) Intermediate half-level.
!
! ------------------- 2 p_(j=3) Archived p level.
!
! ======================== p_(j=2) Intermediate half-level.
!
! ------------------- 1 p_(j=1) Archived p level.
!
! ============================== p_(j=0) = p_o Prescribed p at bottom of
! column.
!
! The complete complement of pressure levels are defined as follows:
!
! Let J = number of archived pressure levels (to be written as nalvl in
! actual code).
!
! Archived pressure levels p_j for j = 1, 3, ..., 2J-1.
!
! Intermediate half-levels p_j = (p_(j+1) + p(j-1))/2, j = 2, 4, ..., 2J-2.
!
! Prescribed limits of in-
! tegration are given by p_2J = p_t < p_(2J-1) everywhere,
! p_0 = p_o > p_s everywhere, where p_s
! is the surface pressure.
!
! Layer thicknesses dp_j = p_(j-1) - p_(j+1), j = 1, 3, ..., 2J-1.
!
! The vertical integral of a quantity M is then given by
!
! --
! I = > (beta_j * M_j * dp_j) / g
! --
! j = 1, 2J-1, 2
!
! where g is the gravitational acceleration and, recalling that p_s is
! the surface pressure
!
! beta_j = 1 if p_(j-1) < p_s,
!
! beta_j = 0 if p_(j+1) > p_s,
!
! beta_j = ( p_s - p_(j+1) ) / ( p_(j-1) - p_(j+1) )
! if p_(j-1) > p_s > p_(j+1).
!
!
! REFERENCE:
!
! Trenberth, K.E., and A. Solomon, 1994: The global heat balance: heat
! transports in the atmosphere and ocean. Clim. Dyn., 10, 107-134. (See
! pages 111-112.)
!
!>============================================================================<
FUNCTION V_INT_PRESS( integrand, msg_val, arch_p_lvl, p_o, p_s, p_t, T_2_B )
IMPLICIT NONE
REAL :: V_INT_PRESS ! Mass weighted vertical integral of specified
! integrand. If, for example, the units of the
! integrand are J kg^(-1), the units of
! V_INT_PRESS are J m^(-2); i.e.:
!
! / p_o
! |
! 1 |
! V_INT_PRESS = - | integrand dp
! g |
! |
! / p_t
REAL, DIMENSION(:), INTENT(IN) :: integrand
! Quantity to be vertically integrated in
! pressure coordinates. Required to have SI
! units such as J kg^(-1). For a fixed time,
! field, and grid point, this will be a
! one dimensional array with nalvl elements
! where nalvl is the number of archived
! pressure levels (more about nalvl below).
REAL, INTENT(IN) :: msg_val ! Missing value of integrand.
REAL, DIMENSION(:), INTENT(IN) :: arch_p_lvl
! Pressure value for each archived pressure
! level, Pascals.
REAL, INTENT(IN) :: p_o ! Fixed pressure at bottom of column.
! Pascals.
REAL, INTENT(IN) :: p_s ! Scalar representing surface pressure for
! fixed time, field, and grid point. Pascals.
REAL, INTENT(IN) :: p_t ! Fixed pressure at top of column. Pascals.
LOGICAL, INTENT(IN) :: T_2_B ! If integrand and arch_p_lvl are ordered top
! to bottom then T_2_B must be assigned the
! logical value .TRUE.
!>----------------------------------------------------------------------------<
! Declaration of allocatable arrays:
REAL, DIMENSION(:), ALLOCATABLE :: integrand_B_2_T
! integrand ordered bottom to top.
REAL, DIMENSION(:), ALLOCATABLE :: arch_p_lvl_B_2_T
! arch_p_lvl ordered bottom to top.
REAL, DIMENSION(:), ALLOCATABLE :: p
! Full complement of archived pressure levels,
! intermediate half levels, and fixed pressures
! at bottom and top of atmospheric column, in
! Pascals.
REAL, DIMENSION(:), ALLOCATABLE :: dp
! Layer thicknesses, where a layer is defined
! as the positive difference between intermedi-
! ate half levels. Pascals.
REAL, DIMENSION(:), ALLOCATABLE :: beta
! Beta multiplier as discusssed in, for example,
! Trenberth and Solomon (1994). Dimensionless.
REAL, DIMENSION(:), ALLOCATABLE :: m
! Array to store nalvl elements of integrand
! mapped into 2*nalvl+1 elements of m, where
! nalvl is the number of archived pressure
! levels.
!>----------------------------------------------------------------------------<
! Declaration of local variables:
INTEGER :: nalvl ! Number of archived pressure levels, ie the
! size of the assumed-shape arrays arch_p_lvl
! and integrand. Minimum number of archived
! pressure levels is 3.
INTEGER :: j ! DO loop control variable and discrete inte-
! gration index employed by Trenberth and
! Solomon (1994). [This is not J (ie nalvl), the
! number of archived pressure levels.]
REAL :: p_sum ! Variable for storing partial sum of discrete
! integration summation.
REAL, PARAMETER :: g = 9.80616
! Gravitational acceleration, m s^(-2). Default.
! (Be sure to check the source of the data with
! regard to the most appropriate value of g.)
INTEGER :: alloc_err ! Memory allocation error flag.
INTEGER :: dealloc_err ! Memory deallocation error flag.
!>----------------------------------------------------------------------------<
! Assign nalvl, the number of archived pressure levels:
nalvl = SIZE(arch_p_lvl) ! Assign number of elements of arch_p_lvl to
! nalvl.
IF ( nalvl < 3 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "The number of archived pressure levels "
WRITE (*,*) "nalvl must be greater than or equal to 3:"
WRITE (*,*) nalvl
WRITE (*,*) "Is there an INTERFACE, END INTERFACE "
WRITE (*,*) "block for this subroutine in calling "
WRITE (*,*) "program? "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
END IF
!>----------------------------------------------------------------------------<
! IF integrand and arch_p_lvl are ordered top to bottom then reorder
! bottom to top in allocated memory to be consistent with algorithm:
ALLOCATE ( integrand_B_2_T(1:nalvl) )
ALLOCATE ( arch_p_lvl_B_2_T(1:nalvl) )
IF ( T_2_B ) THEN
arch_p_lvl_B_2_T(1:nalvl) = arch_p_lvl(nalvl:1:-1)
integrand_B_2_T(1:nalvl) = integrand(nalvl:1:-1)
ELSE IF ( .NOT. ( T_2_B ) ) THEN
arch_p_lvl_B_2_T(:) = arch_p_lvl(:)
integrand_B_2_T(:) = integrand(:)
END IF
!>----------------------------------------------------------------------------<
! Allocate p and compute and assign full complement of p:
ALLOCATE( p(0:2*nalvl), STAT = alloc_err )
IF ( alloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Insufficient memory to allocate p "
WRITE (*,*) "or p has already been allocated. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
p(1:2*nalvl-1:2) = arch_p_lvl_B_2_T(1:nalvl)
! Assign archived pressure levels,
! j = 1, 3, ..., 2*nalvl - 1.
DO j = 3, 2*nalvl-3, 2
IF ( .NOT. ( ( p(j-2) > p(j) ) .AND. ( p(j) > p(j+2) ) ) ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Archived pressure levels must decrease "
WRITE (*,*) "monotonically. "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
END IF
END DO
p(2:2*nalvl-2:2) = ( p(1:2*nalvl-3:2) + p(3:2*nalvl-1:2) ) / 2.
! Compute intermediate half-levels from archived
! pressure levels,
! j = 2, 4, ..., 2*nalvl - 2.
p(0) = p_o ! Assign fixed pressure at bottom of column,
! j = 0.
p(2*nalvl) = p_t ! Assign fixed pressure at top of column,
! j = 2*nalvl.
IF ( ( p(0) < p(1) ) .OR. ( p(2*nalvl) > p(2*nalvl-1) ) ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Either p_o is not greater than the bottom"
WRITE (*,*) "archived pressure, or p_t is not less "
WRITE (*,*) "than the top archived pressure. "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
ELSE IF ( .NOT. ( ( p(0) > p_s ) .AND. ( p_s > p(2*nalvl) ) ) ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "The surface pressure p_s does not lie "
WRITE (*,*) "between p_o and p_t. "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
END IF
!>----------------------------------------------------------------------------<
! Compute and assign layer thicknesses:
ALLOCATE( dp(0:2*nalvl), STAT = alloc_err )
IF ( alloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Insufficient memory to allocate dp "
WRITE (*,*) "or dp has already been allocated. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
dp(1:2*nalvl-1:2) = p(0:2*nalvl-2:2) - p(2:2*nalvl:2)
! Compute layer thicknesses,
! j = 1, 3, ..., 2*nalvl - 1.
dp(0:2*nalvl:2) = 0. ! Assign even indexed dp as 0.,
! j = 0, 2, ..., 2*nalvl.
IF ( ANY( dp < 0. ) ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Layer thicknesses must be nonnegative. "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
END IF
!>----------------------------------------------------------------------------<
! Compute and assign beta:
ALLOCATE( beta(0:2*nalvl), STAT = alloc_err )
IF ( alloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Insufficient memory to allocate beta "
WRITE (*,*) "or beta has already been allocated. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
DO j = 1, 2*nalvl - 1, 2 ! Computing beta for j = 1, 3, ..., 2*nalvl - 1.
IF ( ( p_s < p(j-1) ) .AND. ( p_s > p(j+1) ) ) THEN
beta(j) = ( p_s - p(j+1) ) / ( p(j-1) - p(j+1) )
ELSE IF ( p(j-1) < p_s ) THEN
beta(j) = 1.
ELSE IF ( p(j+1) > p_s ) THEN
beta(j) = 0.
END IF
END DO
beta(0:2*nalvl:2) = 0. ! Assign even indexed beta as 0.,
! j = 0, 2, ..., 2*nalvl.
IF ( ANY( beta < 0. ) ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Beta factors must be nonnegative. "
WRITE (*,*) "Execution halted in function V_INT_PRESS."
STOP
END IF
!>----------------------------------------------------------------------------<
! Map nalvl elements of integrand_B_2_T into 2*nalvl+1 elements of m, accounting
! for missing data after the mapping.
ALLOCATE( m(0:2*nalvl), STAT = alloc_err )
IF ( alloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Insufficient memory to allocate m "
WRITE (*,*) "or m has already been allocated. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
m(1:2*nalvl-1:2) = integrand_B_2_T(1:nalvl)
! Assign elements of m,
! j = 1, 3, ..., 2*nalvl - 1.
m(0:2*nalvl:2) = 0. ! Assign even indexed m as 0.,
! j = 0, 2, ..., 2*nalvl.
WHERE ( m == msg_val ) m = 0. ! Account for missing values so that contribu-
! tion to discrete integration will be zero.
!>----------------------------------------------------------------------------<
! Compute vertical integral.
p_sum = 0. ! Initialize p_sum.
DO j = 1, 2*nalvl - 1, 2
p_sum = p_sum + beta(j) * m(j) * dp(j)
! Cumulative sum of V_INT_PRESS.
END DO
V_INT_PRESS = (1./g) * p_sum ! Assign V_INT_PRESS.
!>----------------------------------------------------------------------------<
! Deallocate allocated arrays:
DEALLOCATE( p, STAT = dealloc_err )
IF ( dealloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Deallocation of memory allocated "
WRITE (*,*) "to p unsuccessful. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
DEALLOCATE( dp, STAT = dealloc_err )
IF ( dealloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Deallocation of memory allocated "
WRITE (*,*) "to dp unsuccessful. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
DEALLOCATE( beta, STAT = dealloc_err )
IF ( dealloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Deallocation of memory allocated "
WRITE (*,*) "to beta unsuccessful. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
DEALLOCATE( m, STAT = dealloc_err )
IF ( dealloc_err /= 0 ) THEN
WRITE (*,*) "ERROR CONDITION: "
WRITE (*,*) "Deallocation of memory allocated "
WRITE (*,*) "to m unsuccessful. "
WRITE (*,*) "Execution halted in function "
WRITE (*,*) "V_INT_PRESS. "
STOP
END IF
DEALLOCATE( integrand_B_2_T )
DEALLOCATE( arch_p_lvl_B_2_T )
!>----------------------------------------------------------------------------<
END FUNCTION V_INT_PRESS
!>----------------------------------------------------------------------------<