init_reference_pressure Subroutine

    subroutine init_reference_pressure()
        ! Module(s) from CAM-SIMA.
        use cam_abortutils, only: check_allocate
        use cam_history_support, only: add_vert_coord
        use cam_logfile, only: debugout_debug, debugout_verbose
        use dyn_comp, only: dyn_debug_print, mpas_dynamical_core
        use dynconst, only: constant_p0 => pref
        use ref_pres, only: ref_pres_init
        use std_atm_profile, only: std_atm_pres
        use string_utils, only: stringify
        use vert_coord, only: pver, pverp
        ! Module(s) from CESM Share.
        use shr_kind_mod, only: kind_r8 => shr_kind_r8
        ! Module(s) from MPAS.
        use dyn_mpas_subdriver, only: kind_dyn_mpas => mpas_dynamical_core_real_kind

        character(*), parameter :: subname = 'dyn_grid::init_reference_pressure'
        ! Number of pure pressure levels at model top.
        integer, parameter :: num_pure_p_lev = 0
        integer :: ierr
        integer :: k, l
        ! Reference pressure (Pa) at layer interfaces and midpoints.
        real(kind_r8), allocatable :: p_ref_int(:)
        real(kind_r8), allocatable :: p_ref_mid(:)
        ! In MPAS, the vertical coordinate is denoted as lowercase "zeta", a Greek alphabet.
        ! `zu` denotes zeta at u-wind levels (i.e., at layer midpoints).
        ! `zw` denotes zeta at w-wind levels (i.e., at layer interfaces).
        ! `dzw` denotes the delta/difference between `zw`.
        ! `rdzw` denotes the reciprocal of `dzw`.
        real(kind_r8), allocatable :: dzw(:)
        real(kind_dyn_mpas), pointer :: rdzw(:)
        real(kind_r8), pointer :: zu(:) ! CANNOT be safely deallocated because `add_vert_coord`
                                        ! just uses pointers to point at it internally.
        real(kind_r8), pointer :: zw(:) ! CANNOT be safely deallocated because `add_vert_coord`
                                        ! just uses pointers to point at it internally.

        call dyn_debug_print(debugout_debug, subname // ' entered')

        nullify(rdzw)
        nullify(zu)
        nullify(zw)

        ! Compute reference height.
        call mpas_dynamical_core % get_variable_pointer(rdzw, 'mesh', 'rdzw')

        allocate(dzw(pver), stat=ierr)
        call check_allocate(ierr, subname, 'dzw(pver)', 'dyn_grid', __LINE__)

        dzw(:) = 1.0_kind_r8 / real(rdzw(:), kind_r8)

        nullify(rdzw)

        allocate(zw(pverp), stat=ierr)
        call check_allocate(ierr, subname, 'zw(pverp)', 'dyn_grid', __LINE__)
        allocate(zu(pver), stat=ierr)
        call check_allocate(ierr, subname, 'zu(pver)', 'dyn_grid', __LINE__)

        ! In MPAS, zeta coordinates are stored in increasing order (i.e., bottom to top of atmosphere).
        ! In CAM-SIMA, however, index order is reversed (i.e., top to bottom of atmosphere).
        ! Compute in reverse below.
        zw(pverp) = 0.0_kind_r8

        do k = pver, 1, -1
            zw(k) = zw(k + 1) + dzw(pver - k + 1)
            zu(k) = 0.5_kind_r8 * (zw(k + 1) + zw(k))
        end do

        deallocate(dzw)

        ! Register zeta coordinates with history.
        call add_vert_coord('ilev', pverp, 'Height (zeta) level at layer interfaces', 'm', zw, &
            positive='up')
        call add_vert_coord('lev', pver, 'Height (zeta) level at layer midpoints', 'm', zu, &
            positive='up')

        ! Compute reference pressure from reference height.
        allocate(p_ref_int(pverp), stat=ierr)
        call check_allocate(ierr, subname, 'p_ref_int(pverp)', 'dyn_grid', __LINE__)

        call std_atm_pres(zw, p_ref_int, user_specified_ps=constant_p0)

        allocate(p_ref_mid(pver), stat=ierr)
        call check_allocate(ierr, subname, 'p_ref_mid(pver)', 'dyn_grid', __LINE__)

        p_ref_mid(:) = 0.5_kind_r8 * (p_ref_int(1:pver) + p_ref_int(2:pverp))

        ! Print a nice table of reference height and pressure.
        call dyn_debug_print(debugout_verbose, 'Reference layer information:')
        call dyn_debug_print(debugout_verbose, '----- | -------------- | --------------')
        call dyn_debug_print(debugout_verbose, 'Index |     Height (m) | Pressure (hPa)')

        do k = 1, pver
            call dyn_debug_print(debugout_verbose, repeat('-', 5) // ' |  ' // stringify([zw(k)]) // &
                ' |  ' // stringify([p_ref_int(k) / 100.0_kind_r8]))

            l = len(stringify([k]))

            call dyn_debug_print(debugout_verbose, repeat(' ', 5 - l) // stringify([k]) // ' |  ' // stringify([zu(k)]) // &
                ' |  ' // stringify([p_ref_mid(k) / 100.0_kind_r8]))
        end do

        k = pverp

        call dyn_debug_print(debugout_verbose, repeat('-', 5) // ' |  ' // stringify([zw(k)]) // &
            ' |  ' // stringify([p_ref_int(k) / 100.0_kind_r8]))

        call ref_pres_init(p_ref_int, p_ref_mid, num_pure_p_lev)

        deallocate(p_ref_int)
        deallocate(p_ref_mid)

        nullify(zu)
        nullify(zw)

        call dyn_debug_print(debugout_debug, subname // ' completed')
    end subroutine init_reference_pressure