LCOV - code coverage report
Current view: | top level - utils - cam_thermo.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 209 | 555 | 37.7 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 17 | 32 | 53.1 %

              Line data    Source code
           1             : ! cam_thermo module provides interfaces to compute thermodynamic quantities
           2             : module cam_thermo
           3             :
           4             :    use shr_kind_mod,    only: r8 => shr_kind_r8
           5             :    use cam_abortutils,  only: endrun
           6             :    use air_composition, only: thermodynamic_active_species_num
           7             :    use air_composition, only: thermodynamic_active_species_idx
           8             :    use air_composition, only: thermodynamic_active_species_idx_dycore
           9             :    use air_composition, only: thermodynamic_active_species_cp
          10             :    use air_composition, only: thermodynamic_active_species_R
          11             :    use air_composition, only: thermodynamic_active_species_mwi
          12             :    use air_composition, only: thermodynamic_active_species_kv
          13             :    use air_composition, only: thermodynamic_active_species_kc
          14             :    use air_composition, only: thermodynamic_active_species_liq_num
          15             :    use air_composition, only: thermodynamic_active_species_ice_num
          16             :    use air_composition, only: thermodynamic_active_species_liq_idx
          17             :    use air_composition, only: thermodynamic_active_species_liq_idx_dycore
          18             :    use air_composition, only: thermodynamic_active_species_ice_idx
          19             :    use air_composition, only: thermodynamic_active_species_ice_idx_dycore
          20             :    use air_composition, only: dry_air_species_num
          21             :    use air_composition, only: enthalpy_reference_state
          22             :    use air_composition, only: mmro2, mmrn2, o2_mwi, n2_mwi, mbar
          23             :
          24             :    implicit none
          25             :    private
          26             :    save
          27             :
          28             :    ! subroutines to compute thermodynamic quantities
          29             :    !
          30             :    ! See Lauritzen et al. (2018) for formulae
          31             :    !     DOI: 10.1029/2017MS001257
          32             :    !     https://opensky.ucar.edu/islandora/object/articles:21929
          33             :
          34             :    ! cam_thermo_init: Initialize constituent dependent properties
          35             :    public :: cam_thermo_init
          36             :    ! cam_thermo_dry_air_update: Update dry air composition dependent properties
          37             :    public :: cam_thermo_dry_air_update
          38             :    ! cam_thermo_water_update: Update water dependent properties
          39             :    public :: cam_thermo_water_update
          40             :    ! get_enthalpy: enthalpy quantity = dp*cp*T
          41             :    public :: get_enthalpy
          42             :    ! get_virtual_temp: virtual temperature
          43             :    public :: get_virtual_temp
          44             :    ! get_sum_species: sum of thermodynamically active species:
          45             :    !                  Note: dp = dp_dry * sum_species
          46             :    public :: get_sum_species
          47             :    ! get_virtual_theta: virtual potential temperature
          48             :    public :: get_virtual_theta
          49             :    ! cam_thermo_calc_kappav: update species dependent kappa for FV dycore
          50             :    public :: cam_thermo_calc_kappav
          51             :    ! get_dp: pressure level thickness from dry dp and dry mixing ratios
          52             :    public :: get_dp
          53             :    ! get_pmid_from_dp: full level pressure from dp (approximation depends on dycore)
          54             :    public :: get_pmid_from_dp
          55             :    ! get_ps: surface pressure
          56             :    public :: get_ps
          57             :    ! get_gz: geopotential
          58             :    public :: get_gz
          59             :    ! get_Richardson_number: Richardson number at layer interfaces
          60             :    public :: get_Richardson_number
          61             :    ! get_kappa_dry: (generalized) dry kappa = R_dry/cp_dry
          62             :    public :: get_kappa_dry
          63             :    ! get_dp_ref: reference pressure layer thickness (include topography)
          64             :    public :: get_dp_ref
          65             :    ! get_molecular_diff_coef: molecular diffusion and thermal conductivity
          66             :    public :: get_molecular_diff_coef
          67             :    ! get_molecular_diff_coef_reference: reference vertical profile of density,
          68             :    !                             molecular diffusion and thermal conductivity
          69             :    public :: get_molecular_diff_coef_reference
          70             :    ! get_rho_dry: dry density from temperature (temp) and
          71             :    !              pressure (dp_dry and tracer)
          72             :    public :: get_rho_dry
          73             :    ! get_exner: Exner pressure
          74             :    public :: get_exner
          75             :    ! get_hydrostatic_energy: Vertically integrated total energy
          76             :    public :: get_hydrostatic_energy
          77             :
          78             :    ! Public variables
          79             :    ! mixing_ratio options
          80             :    integer, public, parameter :: DRY_MIXING_RATIO = 1
          81             :    integer, public, parameter :: MASS_MIXING_RATIO = 2
          82             :    !---------------  Variables below here are for WACCM-X ---------------------
          83             :    ! kmvis: molecular viscosity      kg/m/s
          84             :    real(r8), public, protected, allocatable :: kmvis(:,:,:)
          85             :    ! kmcnd: molecular conductivity   J/m/s/K
          86             :    real(r8), public, protected, allocatable :: kmcnd(:,:,:)
          87             :
          88             :    !------------- Variables for consistent themodynamics --------------------
          89             :    !
          90             :
          91             :    !
          92             :    ! Interfaces for public routines
          93             :    interface get_gz
          94             :       ! get_gz_geopotential (with dp_dry, ptop, temp, and phis as input)
          95             :       module procedure get_gz_from_dp_dry_ptop_temp_1hd
          96             :       ! get_gz_given_dp_Tv_Rdry: geopotential (with dp,dry R and Tv as input)
          97             :       module procedure get_gz_given_dp_Tv_Rdry_1hd
          98             :       module procedure get_gz_given_dp_Tv_Rdry_2hd
          99             :    end interface get_gz
         100             :
         101             :    interface get_enthalpy
         102             :       module procedure get_enthalpy_1hd
         103             :       module procedure get_enthalpy_2hd
         104             :    end interface get_enthalpy
         105             :
         106             :    interface get_virtual_temp
         107             :       module procedure get_virtual_temp_1hd
         108             :       module procedure get_virtual_temp_2hd
         109             :    end interface get_virtual_temp
         110             :
         111             :    interface get_sum_species
         112             :       module procedure get_sum_species_1hd
         113             :       module procedure get_sum_species_2hd
         114             :    end interface get_sum_species
         115             :
         116             :    interface get_dp
         117             :       module procedure get_dp_1hd
         118             :       module procedure get_dp_2hd
         119             :    end interface get_dp
         120             :
         121             :    interface get_pmid_from_dp
         122             :       module procedure get_pmid_from_dpdry_1hd
         123             :       module procedure get_pmid_from_dp_1hd
         124             :    end interface get_pmid_from_dp
         125             :
         126             :    interface get_exner
         127             :       module procedure get_exner_1hd
         128             :    end interface get_exner
         129             :
         130             :    interface get_virtual_theta
         131             :       module procedure get_virtual_theta_1hd
         132             :    end interface get_virtual_theta
         133             :
         134             :    interface get_Richardson_number
         135             :       module procedure get_Richardson_number_1hd
         136             :    end interface get_Richardson_number
         137             :
         138             :    interface get_ps
         139             :       module procedure get_ps_1hd
         140             :       module procedure get_ps_2hd
         141             :    end interface get_ps
         142             :
         143             :    interface get_kappa_dry
         144             :       module procedure get_kappa_dry_1hd
         145             :       module procedure get_kappa_dry_2hd
         146             :    end interface get_kappa_dry
         147             :
         148             :    interface get_dp_ref
         149             :       module procedure get_dp_ref_1hd
         150             :       module procedure get_dp_ref_2hd
         151             :    end interface get_dp_ref
         152             :
         153             :    interface get_rho_dry
         154             :       module procedure get_rho_dry_1hd
         155             :       module procedure get_rho_dry_2hd
         156             :    end interface get_rho_dry
         157             :
         158             :    interface get_molecular_diff_coef
         159             :       module procedure get_molecular_diff_coef_1hd
         160             :       module procedure get_molecular_diff_coef_2hd
         161             :    end interface get_molecular_diff_coef
         162             :
         163             :    interface cam_thermo_calc_kappav
         164             :       ! Since this routine is currently only used by the FV dycore,
         165             :       !    a 1-d interface is not needed (but can easily be added)
         166             :       module procedure cam_thermo_calc_kappav_2hd
         167             :    end interface cam_thermo_calc_kappav
         168             :
         169             :    interface get_hydrostatic_energy
         170             :       module procedure get_hydrostatic_energy_1hd
         171             :       ! This routine is currently only called from the physics so a
         172             :       !    2-d interface is not needed (but can easily be added)
         173             :    end interface get_hydrostatic_energy
         174             :
         175             :    integer, public, parameter :: thermo_budget_num_vars = 10
         176             :    integer, public, parameter :: wvidx = 1
         177             :    integer, public, parameter :: wlidx = 2
         178             :    integer, public, parameter :: wiidx = 3
         179             :    integer, public, parameter :: seidx = 4 ! enthalpy or internal energy (W/m2) index
         180             :    integer, public, parameter :: poidx = 5 ! surface potential or potential energy index
         181             :    integer, public, parameter :: keidx = 6 ! kinetic energy index
         182             :    integer, public, parameter :: mridx = 7
         183             :    integer, public, parameter :: moidx = 8
         184             :    integer, public, parameter :: ttidx = 9
         185             :    integer, public, parameter :: teidx = 10
         186             :    character (len = 2)  ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars  = &
         187             :         (/"WV"  ,"WL"  ,"WI"  ,"SE"   ,"PO"   ,"KE"   ,"MR"   ,"MO"   ,"TT"   ,"TE"   /)
         188             :    character (len = 46) ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars_descriptor = (/&
         189             :         "Total column water vapor                      ",&
         190             :         "Total column liquid water                     ",&
         191             :         "Total column frozen water                     ",&
         192             :         "Total column enthalpy or internal energy      ",&
         193             :         "Total column srf potential or potential energy",&
         194             :         "Total column kinetic energy                   ",&
         195             :         "Total column wind axial angular momentum      ",&
         196             :         "Total column mass axial angular momentum      ",&
         197             :         "Total column test_tracer                      ",&
         198             :         "Total column energy (ke + se + po)            "/)
         199             :
         200             :    character (len = 14), public, dimension(thermo_budget_num_vars)  :: &
         201             :         thermo_budget_vars_unit = (/&
         202             :         "kg/m2        ","kg/m2        ","kg/m2        ","J/m2         ",&
         203             :         "J/m2         ","J/m2         ","kg*m2/s*rad2 ","kg*m2/s*rad2 ",&
         204             :         "kg/m2        ","J/m2         "/)
         205             :    logical              ,public, dimension(thermo_budget_num_vars) :: thermo_budget_vars_massv = (/&
         206             :         .true.,.true.,.true.,.false.,.false.,.false.,.false.,.false.,.true.,.false./)
         207             : CONTAINS
         208             :
         209             :    !===========================================================================
         210             :
         211        1536 :    subroutine cam_thermo_init()
         212             :       use shr_infnan_mod,  only: assignment(=), shr_infnan_qnan
         213             :       use ppgrid,          only: pcols, pver, pverp, begchunk, endchunk
         214             :
         215             :       integer                     :: ierr
         216             :       character(len=*), parameter :: subname = "cam_thermo_init"
         217             :       character(len=*), parameter :: errstr = subname//": failed to allocate "
         218             :
         219             :       !------------------------------------------------------------------------
         220             :       !  Allocate constituent dependent properties
         221             :       !------------------------------------------------------------------------
         222        4608 :       allocate(kmvis(pcols,pverp,begchunk:endchunk), stat=ierr)
         223        1536 :       if (ierr /= 0) then
         224           0 :          call endrun(errstr//"kmvis")
         225             :       end if
         226        4608 :       allocate(kmcnd(pcols,pverp,begchunk:endchunk), stat=ierr)
         227        1536 :       if (ierr /= 0) then
         228           0 :          call endrun(errstr//"kmcnd")
         229             :       end if
         230             :
         231             :       !------------------------------------------------------------------------
         232             :       !  Initialize constituent dependent properties
         233             :       !------------------------------------------------------------------------
         234        1536 :       kmvis(:pcols,  :pver, begchunk:endchunk) = shr_infnan_qnan
         235        1536 :       kmcnd(:pcols,  :pver, begchunk:endchunk) = shr_infnan_qnan
         236             :
         237        1536 :    end subroutine cam_thermo_init
         238             :    !
         239             :    !***************************************************************************
         240             :    !
         241             :    ! cam_thermo_dry_air_update: update dry air species dependent constants for physics
         242             :    !
         243             :    !***************************************************************************
         244             :    !
         245           0 :    subroutine cam_thermo_dry_air_update(mmr, T, lchnk, ncol, to_dry_factor)
         246             :       use air_composition, only: dry_air_composition_update
         247             :       use string_utils,    only: int2str
         248             :       !------------------------------Arguments----------------------------------
         249             :       !(mmr = dry mixing ratio, if not use to_dry_factor to convert)
         250             :       real(r8),           intent(in) :: mmr(:,:,:) ! constituents array
         251             :       real(r8),           intent(in) :: T(:,:)     ! temperature
         252             :       integer,            intent(in) :: lchnk      ! Chunk number
         253             :       integer,            intent(in) :: ncol       ! number of columns
         254             :       real(r8), optional, intent(in) :: to_dry_factor(:,:)!if mmr moist convert
         255             :       !
         256             :       !---------------------------Local storage-------------------------------
         257           0 :       real(r8):: sponge_factor(SIZE(mmr, 2))
         258             :       character(len=*), parameter :: subname = 'cam_thermo_update: '
         259             :
         260           0 :       if (present(to_dry_factor)) then
         261           0 :         if (SIZE(to_dry_factor, 1) /= ncol) then
         262           0 :           call endrun(subname//'DIM 1 of to_dry_factor is'//int2str(SIZE(to_dry_factor,1))//'but should be'//int2str(ncol))
         263             :         end if
         264             :       end if
         265             :
         266           0 :       sponge_factor = 1.0_r8
         267           0 :       call dry_air_composition_update(mmr, lchnk, ncol, to_dry_factor=to_dry_factor)
         268           0 :       call get_molecular_diff_coef(T(:ncol,:), .true., sponge_factor, kmvis(:ncol,:,lchnk), &
         269           0 :            kmcnd(:ncol,:,lchnk), tracer=mmr(:ncol,:,:), fact=to_dry_factor,  &
         270           0 :            active_species_idx_dycore=thermodynamic_active_species_idx)
         271           0 :     end subroutine cam_thermo_dry_air_update
         272             :     !
         273             :     !***************************************************************************
         274             :     !
         275             :     ! cam_thermo_water+update: update water species dependent constants for physics
         276             :     !
         277             :     !***************************************************************************
         278             :     !
         279     2984544 :     subroutine cam_thermo_water_update(mmr, lchnk, ncol, vcoord, to_dry_factor)
         280             :       use air_composition, only: water_composition_update
         281             :       !-----------------------------------------------------------------------
         282             :       ! Update the physics "constants" that vary
         283             :       !-------------------------------------------------------------------------
         284             :
         285             :       !------------------------------Arguments----------------------------------
         286             :
         287             :       real(r8),           intent(in) :: mmr(:,:,:) ! constituents array
         288             :       integer,            intent(in) :: lchnk      ! Chunk number
         289             :       integer,            intent(in) :: ncol       ! number of columns
         290             :       integer,            intent(in) :: vcoord
         291             :       real(r8), optional, intent(in) :: to_dry_factor(:,:)
         292             :       !
         293             :       logical :: lcp
         294             :
         295     4479912 :       call water_composition_update(mmr, lchnk, ncol, vcoord, to_dry_factor=to_dry_factor)
         296     2984544 :     end subroutine cam_thermo_water_update
         297             :
         298             :    !===========================================================================
         299             :
         300             :    !
         301             :    !***********************************************************************
         302             :    !
         303             :    ! Compute enthalpy = cp*T*dp, where dp is pressure level thickness,
         304             :    !    cp is generalized cp and T temperature
         305             :    !
         306             :    ! Note: tracer is in units of m*dp_dry ("mass")
         307             :    !
         308             :    !***********************************************************************
         309             :    !
         310    41558400 :    subroutine get_enthalpy_1hd(tracer_mass, temp, dp_dry,               &
         311    41558400 :         enthalpy, active_species_idx_dycore)
         312             :       use air_composition, only: dry_air_species_num, get_cp_dry
         313             :       ! Dummy arguments
         314             :       ! tracer_mass: tracer array (mass weighted)
         315             :       real(r8),          intent(in)  :: tracer_mass(:,:,:)
         316             :       ! temp: temperature
         317             :       real(r8),          intent(in)  :: temp(:,:)
         318             :       ! dp_dry: dry presure level thickness
         319             :       real(r8),          intent(in)  :: dp_dry(:,:)
         320             :       ! enthalpy: enthalpy in each column: sum cp*T*dp
         321             :       real(r8),          intent(out) :: enthalpy(:,:)
         322             :       !
         323             :       ! active_species_idx_dycore:
         324             :       !    array of indicies for index of thermodynamic active species in
         325             :       !    dycore tracer array (if different from physics index)
         326             :       !
         327             :       integer, optional, intent(in)  :: active_species_idx_dycore(:)
         328             :
         329             :       ! Local vars
         330             :       integer                     :: qdx, itrac
         331             :       character(len=*), parameter :: subname = 'get_enthalpy: '
         332             :
         333             :       !
         334             :       ! "mass-weighted" cp (dp must be dry)
         335             :       !
         336    41558400 :       if (dry_air_species_num == 0) then
         337           0 :          enthalpy(:,:) = thermodynamic_active_species_cp(0) *         &
         338  5444150400 :               dp_dry(:,:)
         339             :       else
         340           0 :          if (present(active_species_idx_dycore)) then
         341             :             call get_cp_dry(tracer_mass, active_species_idx_dycore,           &
         342           0 :                  enthalpy, fact=1.0_r8/dp_dry(:,:))
         343             :          else
         344             :             call get_cp_dry(tracer_mass, thermodynamic_active_species_idx,    &
         345           0 :                  enthalpy, fact=1.0_r8/dp_dry(:,:))
         346             :          end if
         347           0 :          enthalpy(:,:) = enthalpy(:,:) * dp_dry(:,:)
         348             :       end if
         349             :       !
         350             :       ! tracer is in units of m*dp ("mass"), where:
         351             :       !    m is the dry mixing ratio
         352             :       !    dp is the dry pressure level thickness
         353             :       !
         354   166233600 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         355   124675200 :          if (present(active_species_idx_dycore)) then
         356   124675200 :             itrac = active_species_idx_dycore(qdx)
         357             :          else
         358           0 :             itrac = thermodynamic_active_species_idx(qdx)
         359             :          end if
         360             :          enthalpy(:,:) = enthalpy(:,:) +                      &
         361 16374009600 :               (thermodynamic_active_species_cp(qdx) * tracer_mass(:,:,itrac))
         362             :       end do
         363  5485708800 :       enthalpy(:,:) = enthalpy(:,:) * temp(:,:)
         364             :
         365    41558400 :    end subroutine get_enthalpy_1hd
         366             :
         367             :    !===========================================================================
         368             :
         369    10389600 :    subroutine get_enthalpy_2hd(tracer_mass, temp, dp_dry,               &
         370    10389600 :         enthalpy, active_species_idx_dycore)
         371             :       ! Dummy arguments
         372             :       ! tracer_mass: tracer array (mass weighted)
         373             :       real(r8),          intent(in)  :: tracer_mass(:,:,:,:)
         374             :       ! temp: temperature
         375             :       real(r8),          intent(in)  :: temp(:,:,:)
         376             :       ! dp_dry: dry presure level thickness
         377             :       real(r8),          intent(in)  :: dp_dry(:,:,:)
         378             :       ! enthalpy: enthalpy in each column: sum cp*T*dp
         379             :       real(r8),          intent(out) :: enthalpy(:,:,:)
         380             :       !
         381             :       ! active_species_idx_dycore:
         382             :       !    array of indicies for index of thermodynamic active species in
         383             :       !    dycore tracer array (if different from physics index)
         384             :       !
         385             :       integer, optional, intent(in)  :: active_species_idx_dycore(:)
         386             :
         387             :       ! Local variables
         388             :       integer                     :: jdx
         389             :       character(len=*), parameter :: subname = 'get_enthalpy_2hd: '
         390             :
         391    51948000 :       do jdx = 1, SIZE(tracer_mass, 2)
         392             :          call get_enthalpy(tracer_mass(:, jdx, :, :), temp(:, jdx, :),     &
         393             :               dp_dry(:, jdx, :), enthalpy(:, jdx, :),                   &
         394    51948000 :               active_species_idx_dycore=active_species_idx_dycore)
         395             :       end do
         396             :
         397    10389600 :    end subroutine get_enthalpy_2hd
         398             :
         399             :    !===========================================================================
         400             :
         401             :    !**************************************************************************
         402             :    !
         403             :    ! get_virtual_temp: Compute virtual temperature T_v
         404             :    !
         405             :    ! tracer is in units of dry mixing ratio unless optional argument
         406             :    !    dp_dry is present in which case tracer is in units of "mass" (=m*dp)
         407             :    !
         408             :    ! If temperature is not supplied then just return factor that T
         409             :    !    needs to be multiplied by to get T_v
         410             :    !
         411             :    !**************************************************************************
         412             :    !
         413   623376000 :    subroutine get_virtual_temp_1hd(tracer, T_v, temp, dp_dry, sum_q,          &
         414   311688000 :         active_species_idx_dycore)
         415             :       use cam_abortutils,  only: endrun
         416             :       use string_utils,    only: int2str
         417             :       use air_composition, only: dry_air_species_num, get_R_dry
         418             :
         419             :       ! Dummy Arguments
         420             :       ! tracer: tracer array
         421             :       real(r8),           intent(in)  :: tracer(:, :, :)
         422             :       ! T_v: virtual temperature
         423             :       real(r8),           intent(out) :: T_v(:, :)
         424             :       ! temp: temperature
         425             :       real(r8), optional, intent(in)  :: temp(:, :)
         426             :       ! dp_dry: dry pressure level thickness
         427             :       real(r8), optional, intent(in)  :: dp_dry(:, :)
         428             :       ! sum_q: sum tracer
         429             :       real(r8), optional, intent(out) :: sum_q(:, :)
         430             :       !
         431             :       ! array of indicies for index of thermodynamic active species in
         432             :       !    dycore tracer array (if different from physics index)
         433             :       !
         434             :       integer, optional,  intent(in) :: active_species_idx_dycore(:)
         435             :
         436             :       ! Local Variables
         437             :       integer                     :: itrac, qdx
         438   623376000 :       real(r8)                    :: sum_species(SIZE(tracer, 1), SIZE(tracer, 2))
         439   623376000 :       real(r8)                    :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
         440   623376000 :       real(r8)                    :: Rd(SIZE(tracer, 1), SIZE(tracer, 2))
         441   623376000 :       integer                     :: idx_local(thermodynamic_active_species_num)
         442             :       character(len=*), parameter :: subname = 'get_virtual_temp_1hd: '
         443             :
         444   311688000 :       if (present(active_species_idx_dycore)) then
         445   311688000 :          if (SIZE(active_species_idx_dycore) /=                               &
         446             :               thermodynamic_active_species_num) then
         447             :             call endrun(subname//"SIZE mismatch "//                           &
         448             :                  int2str(SIZE(active_species_idx_dycore))//' /= '//           &
         449           0 :                  int2str(thermodynamic_active_species_num))
         450             :          end if
         451  1558440000 :          idx_local = active_species_idx_dycore
         452             :       else
         453           0 :          idx_local = thermodynamic_active_species_idx
         454             :       end if
         455             :
         456   623376000 :       call get_sum_species(tracer, idx_local, sum_species, dp_dry=dp_dry, factor=factor)
         457             :
         458   311688000 :       call get_R_dry(tracer, idx_local, Rd, fact=factor)
         459 41142816000 :       t_v(:, :) = Rd(:, :)
         460  1246752000 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         461   935064000 :          itrac = idx_local(qdx)
         462           0 :          t_v(:, :) = t_v(:, :) + (thermodynamic_active_species_R(qdx) *       &
         463 >12280*10^7 :               tracer(:, :, itrac) * factor(:, :))
         464             :       end do
         465   311688000 :       if (present(temp)) then
         466 41142816000 :          t_v(:, :) = t_v(:, :) * temp(:, :) / (Rd(:, :) * sum_species)
         467             :       else
         468           0 :          t_v(:, :) = t_v(:, :) / (Rd(:, :) * sum_species)
         469             :       end if
         470   311688000 :       if (present(sum_q)) then
         471 41142816000 :          sum_q = sum_species
         472             :       end if
         473             :
         474   311688000 :    end subroutine get_virtual_temp_1hd
         475             :
         476             :    !===========================================================================
         477             :
         478    77922000 :    subroutine get_virtual_temp_2hd(tracer, T_v, temp, dp_dry, sum_q,          &
         479    77922000 :         active_species_idx_dycore)
         480             :
         481             :       ! Dummy Arguments
         482             :       ! tracer: tracer array
         483             :       real(r8),           intent(in)  :: tracer(:, :, :, :)
         484             :       ! T_v: virtual temperature
         485             :       real(r8),           intent(out) :: T_v(:, :, :)
         486             :       ! temp: temperature
         487             :       real(r8), optional, intent(in)  :: temp(:, :, :)
         488             :       ! dp_dry: dry pressure level thickness
         489             :       real(r8), optional, intent(in)  :: dp_dry(:, :, :)
         490             :       ! sum_q: sum tracer
         491             :       real(r8), optional, intent(out) :: sum_q(:, :, :)
         492             :       !
         493             :       ! array of indicies for index of thermodynamic active species in
         494             :       !    dycore tracer array (if different from physics index)
         495             :       !
         496             :       integer, optional,  intent(in) :: active_species_idx_dycore(:)
         497             :
         498             :       ! Local vars
         499             :       integer                     :: jdx
         500             :       character(len=*), parameter :: subname = 'get_virtual_temp_2hd: '
         501             :
         502             :       ! Rather than do a bunch of copying into temp variables, do the
         503             :       !    combinatorics
         504   389610000 :       do jdx = 1, SIZE(tracer, 2)
         505   389610000 :          if (present(temp) .and. present(dp_dry) .and. present(sum_q)) then
         506             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         507             :                  temp=temp(:, jdx, :), dp_dry=dp_dry(:, jdx, :),              &
         508             :                  sum_q=sum_q(:, jdx, :),                                      &
         509           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         510   311688000 :          else if (present(temp) .and. present(dp_dry)) then
         511             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         512             :                  temp=temp(:, jdx, :), dp_dry=dp_dry(:, jdx, :),              &
         513           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         514   311688000 :          else if (present(temp) .and. present(sum_q)) then
         515             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         516             :                  temp=temp(:, jdx, :), sum_q=sum_q(:, jdx, :),                &
         517   311688000 :                  active_species_idx_dycore=active_species_idx_dycore)
         518           0 :          else if (present(dp_dry) .and. present(sum_q)) then
         519             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         520             :                  dp_dry=dp_dry(:, jdx, :), sum_q=sum_q(:, jdx, :),            &
         521           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         522           0 :          else if (present(temp)) then
         523             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         524             :                  temp=temp(:, jdx, :),                                        &
         525           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         526           0 :          else if (present(dp_dry)) then
         527             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         528             :                  dp_dry=dp_dry(:, jdx, :),                                    &
         529           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         530           0 :          else if (present(sum_q)) then
         531             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         532             :                  sum_q=sum_q(:, jdx, :),                                      &
         533           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         534             :          else
         535             :             call get_virtual_temp(tracer(:, jdx, :, :), T_v(:, jdx, :),       &
         536           0 :                  active_species_idx_dycore=active_species_idx_dycore)
         537             :          end if
         538             :       end do
         539             :
         540    77922000 :    end subroutine get_virtual_temp_2hd
         541             :
         542             :    !===========================================================================
         543             :
         544             :    !
         545             :    !***************************************************************************
         546             :    !
         547             :    ! get_sum_species:
         548             :    !
         549             :    ! Compute sum of thermodynamically active species
         550             :    !
         551             :    ! tracer is in units of dry mixing ratio unless optional argument
         552             :    !    dp_dry is present in which case tracer is in units of "mass" (=m*dp)
         553             :    !
         554             :    !***************************************************************************
         555             :    !
         556   311688000 :    subroutine get_sum_species_1hd(tracer, active_species_idx,                 &
         557   311688000 :         sum_species, dp_dry, factor)
         558             :       use air_composition, only: dry_air_species_num
         559             :
         560             :       ! Dummy arguments
         561             :       ! tracer: Tracer array
         562             :       real(r8),           intent(in)  :: tracer(:, :, :)
         563             :       ! active_species_idx: Index for thermodynamic active tracers
         564             :       integer,            intent(in)  :: active_species_idx(:)
         565             :       ! dp_dry: Dry pressure level thickness.
         566             :       !         If present, then tracer is in units of mass
         567             :       real(r8), optional, intent(in)  :: dp_dry(:, :)
         568             :       ! sum_species: sum species
         569             :       real(r8),           intent(out) :: sum_species(:, :)
         570             :       ! factor: to moist factor
         571             :       real(r8), optional, intent(out) :: factor(:, :)
         572             :       ! Local variables
         573   623376000 :       real(r8) :: factor_loc(SIZE(tracer, 1), SIZE(tracer, 2))
         574             :       integer  :: qdx, itrac
         575   311688000 :       if (present(dp_dry)) then
         576           0 :          factor_loc = 1.0_r8 / dp_dry(:,:)
         577             :       else
         578 40831128000 :          factor_loc = 1.0_r8
         579             :       end if
         580 40831128000 :       sum_species = 1.0_r8 ! all dry air species sum to 1
         581  1246752000 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         582   935064000 :          itrac = active_species_idx(qdx)
         583 >12280*10^7 :          sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_loc(:,:))
         584             :       end do
         585   311688000 :       if (present(factor)) then
         586 41142816000 :          factor = factor_loc
         587             :       end if
         588   311688000 :    end subroutine get_sum_species_1hd
         589             :
         590             :    !===========================================================================
         591             :
         592           0 :    subroutine get_sum_species_2hd(tracer, active_species_idx,                 &
         593           0 :         sum_species,dp_dry, factor)
         594             :
         595             :       ! Dummy arguments
         596             :       ! tracer: Tracer array
         597             :       real(r8),           intent(in)  :: tracer(:, :, :, :)
         598             :       ! active_species_idx: Index for thermodynamic active tracers
         599             :       integer,            intent(in)  :: active_species_idx(:)
         600             :       ! dp_dry: Dry pressure level thickness.
         601             :       !         If present, then tracer is in units of mass
         602             :       real(r8), optional, intent(in)  :: dp_dry(:, :, :)
         603             :       ! sum_species: sum species
         604             :       real(r8),           intent(out) :: sum_species(:, :, :)
         605             :       ! factor: to moist factor
         606             :       real(r8), optional, intent(out) :: factor(:, :, :)
         607             :       ! Local variable
         608             :       integer                     :: jdx
         609             :
         610           0 :       do jdx = 1, SIZE(tracer, 2)
         611           0 :          if (present(dp_dry) .and. present(factor)) then
         612             :             call get_sum_species(tracer(:, jdx, :, :), active_species_idx,    &
         613           0 :                  sum_species(:, jdx, :), dp_dry=dp_dry(:, jdx, :), factor=factor(:, jdx, :))
         614           0 :          else if (present(dp_dry)) then
         615             :             call get_sum_species(tracer(:, jdx, :, :), active_species_idx,    &
         616           0 :                  sum_species(:, jdx, :), dp_dry=dp_dry(:, jdx, :))
         617           0 :          else if (present(factor)) then
         618             :             call get_sum_species(tracer(:, jdx, :, :), active_species_idx,    &
         619           0 :                  sum_species(:, jdx, :), factor=factor(:, jdx, :))
         620             :          else
         621             :             call get_sum_species(tracer(:, jdx, :, :), active_species_idx,    &
         622           0 :                  sum_species(:, jdx, :))
         623             :          end if
         624             :       end do
         625             :
         626           0 :    end subroutine get_sum_species_2hd
         627             :
         628             :    !===========================================================================
         629             :
         630             :    !***************************************************************************
         631             :    !
         632             :    ! get_dp: Compute pressure level thickness from dry pressure and
         633             :    !         thermodynamic active species mixing ratios
         634             :    !
         635             :    ! Tracer can either be in units of dry mixing ratio (mixing_ratio=1) or
         636             :    !    "mass" (=m*dp_dry) (mixing_ratio=2)
         637             :    !
         638             :    !***************************************************************************
         639             :    !
         640    51991200 :    subroutine get_dp_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, dp, ps, ptop)
         641             :      use air_composition,  only: dry_air_species_num
         642             :      use string_utils,    only: int2str
         643             :
         644             :      real(r8), intent(in)  :: tracer(:, :, :)                    ! tracers; quantity specified by mixing_ratio arg
         645             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is dry mixing ratio
         646             :                                                                  ! 2 => tracer is mass (q*dp)
         647             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         648             :      real(r8), intent(in)  :: dp_dry(:, :)                       ! dry pressure level thickness
         649             :      real(r8), intent(out) :: dp(:, :)                           ! pressure level thickness
         650             :      real(r8), optional,intent(out) :: ps(:)                     ! surface pressure (if ps present then ptop
         651             :                                                                  !                   must be present)
         652             :      real(r8), optional,intent(in)  :: ptop                      ! pressure at model top
         653             :
         654             :      integer :: idx, kdx, m_cnst, qdx
         655             :
         656             :      character(len=*), parameter :: subname = 'get_dp_1hd: '
         657             :
         658  6862838400 :      dp = dp_dry
         659    51991200 :      if (mixing_ratio == DRY_MIXING_RATIO) then
         660           0 :        do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         661           0 :          m_cnst = active_species_idx(qdx)
         662           0 :          do kdx = 1, SIZE(tracer, 2)
         663           0 :            do idx = 1, SIZE(tracer, 1)
         664           0 :              dp(idx, kdx) = dp(idx, kdx) + dp_dry(idx, kdx)*tracer(idx, kdx, m_cnst)
         665             :            end do
         666             :          end do
         667             :        end do
         668    51991200 :      else if (mixing_ratio == MASS_MIXING_RATIO) then
         669   207964800 :        do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         670   155973600 :          m_cnst = active_species_idx(qdx)
         671  4263278400 :          do kdx = 1, SIZE(tracer, 2)
         672 20432541600 :            do idx = 1, SIZE(tracer, 1)
         673 20276568000 :              dp(idx, kdx) = dp(idx, kdx) + tracer(idx, kdx, m_cnst)
         674             :            end do
         675             :          end do
         676             :        end do
         677             :      else
         678           0 :         call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
         679             :      end if
         680    51991200 :      if (present(ps)) then
         681           0 :        if (present(ptop)) then
         682           0 :          ps = ptop
         683           0 :          do kdx = 1, SIZE(tracer, 2)
         684           0 :            do idx = 1, SIZE(tracer, 1)
         685           0 :              ps(idx) = ps(idx) + dp(idx, kdx)
         686             :            end do
         687             :          end do
         688             :        else
         689           0 :          call endrun(subname//'if ps is present ptop must be present')
         690             :        end if
         691             :      end if
         692    51991200 :    end subroutine get_dp_1hd
         693             :
         694    12997800 :    subroutine get_dp_2hd(tracer, mixing_ratio, active_species_idx, dp_dry, dp, ps, ptop)
         695             :      ! Version of get_dp for arrays that have a second horizontal index
         696             :      real(r8), intent(in)  :: tracer(:,:,:,:)                    ! tracers; quantity specified by mixing_ratio arg
         697             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is dry mixing ratio
         698             :                                                                  ! 2 => tracer is mass (q*dp)
         699             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         700             :      real(r8), intent(in)  :: dp_dry(:,:,:)                      ! dry pressure level thickness
         701             :      real(r8), intent(out) :: dp(:,:,:)                          ! pressure level thickness
         702             :      real(r8), optional,intent(out) :: ps(:,:)                   ! surface pressure
         703             :      real(r8), optional,intent(in)  :: ptop                      ! pressure at model top
         704             :
         705             :      integer :: jdx
         706             :
         707    64989000 :      do jdx = 1, SIZE(tracer, 2)
         708    64989000 :        if (present(ps)) then
         709             :           call get_dp(tracer(:, jdx, :, :), mixing_ratio, active_species_idx, &
         710           0 :                   dp_dry(:, jdx, :), dp(:, jdx, :), ps=ps(:,jdx), ptop=ptop)
         711             :        else
         712             :           call get_dp(tracer(:, jdx, :, :), mixing_ratio, active_species_idx,  &
         713    51991200 :                   dp_dry(:, jdx, :), dp(:, jdx, :), ptop=ptop)
         714             :        end if
         715             :      end do
         716             :
         717    12997800 :    end subroutine get_dp_2hd
         718             :    !===========================================================================
         719             :
         720             :    !*************************************************************************************************************************
         721             :    !
         722             :    ! compute mid-level (full level) pressure from dry pressure and water tracers
         723             :    !
         724             :    !*************************************************************************************************************************
         725             :    !
         726           0 :    subroutine get_pmid_from_dpdry_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, pmid, pint, dp)
         727             :
         728             :      real(r8), intent(in)  :: tracer(:,:,:)                      ! tracers; quantity specified by mixing_ratio arg
         729             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is mixing ratio
         730             :                                                                  ! 2 => tracer is mass (q*dp)
         731             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         732             :      real(r8), intent(in)  :: dp_dry(:,:)                        ! dry pressure level thickness
         733             :      real(r8), intent(in)  :: ptop                               ! model top pressure
         734             :      real(r8), intent(out) :: pmid(:,:)                          ! mid-level pressure
         735             :      real(r8), optional, intent(out) :: pint(:,:)                ! half-level pressure
         736             :      real(r8), optional, intent(out) :: dp(:,:)                  ! presure level thickness
         737             :
         738           0 :      real(r8) :: dp_local(SIZE(tracer, 1), SIZE(tracer, 2))      ! local pressure level thickness
         739           0 :      real(r8) :: pint_local(SIZE(tracer, 1), SIZE(tracer, 2) + 1)! local interface pressure
         740             :
         741           0 :      call get_dp(tracer, mixing_ratio, active_species_idx, dp_dry, dp_local)
         742             :
         743           0 :      call get_pmid_from_dp(dp_local, ptop, pmid, pint_local)
         744             :
         745           0 :      if (present(pint)) pint=pint_local
         746           0 :      if (present(dp)) dp=dp_local
         747           0 :    end subroutine get_pmid_from_dpdry_1hd
         748             :
         749             :    !===========================================================================
         750             :
         751             :    !*************************************************************************************************************************
         752             :    !
         753             :    ! compute mid-level (full level) pressure
         754             :    !
         755             :    !*************************************************************************************************************************
         756             :    !
         757   311688000 :    subroutine get_pmid_from_dp_1hd(dp, ptop, pmid, pint)
         758             :      use dycore, only: dycore_is
         759             :      real(r8), intent(in)            :: dp(:,:)     ! pressure level thickness
         760             :      real(r8), intent(in)            :: ptop        ! pressure at model top
         761             :      real(r8), intent(out)           :: pmid(:,:)   ! mid (full) level pressure
         762             :      real(r8), optional, intent(out) :: pint(:,:)   ! pressure at interfaces (half levels)
         763             :
         764   623376000 :      real(r8) :: pint_local(SIZE(dp, 1), SIZE(dp,2) + 1)
         765             :      integer  :: kdx
         766             :
         767  1558440000 :      pint_local(:, 1) = ptop
         768  8415576000 :      do kdx = 2, SIZE(dp, 2) + 1
         769 40831128000 :        pint_local(:, kdx) = dp(:, kdx - 1) + pint_local(:, kdx - 1)
         770             :      end do
         771             :
         772   311688000 :      if (dycore_is('LR') .or. dycore_is('FV3')) then
         773           0 :        do kdx = 1, SIZE(dp, 2)
         774           0 :          pmid(:, kdx) = dp(:, kdx) / (log(pint_local(:, kdx + 1)) - log(pint_local(:, kdx)))
         775             :        end do
         776             :      else
         777  8415576000 :        do kdx = 1, SIZE(dp, 2)
         778 40831128000 :          pmid(:, kdx) = 0.5_r8 * (pint_local(:, kdx) + pint_local(:, kdx + 1))
         779             :        end do
         780             :      end if
         781 42389568000 :      if (present(pint)) pint=pint_local
         782   311688000 :    end subroutine get_pmid_from_dp_1hd
         783             :
         784             :    !===========================================================================
         785             :
         786             :    !****************************************************************************************************************
         787             :    !
         788             :    ! Compute Exner pressure
         789             :    !
         790             :    !****************************************************************************************************************
         791             :    !
         792           0 :    subroutine get_exner_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, inv_exner, exner, poverp0)
         793             :      use string_utils,    only: int2str
         794             :      real(r8), intent(in)  :: tracer(:,:,:)                      ! tracers; quantity specified by mixing_ratio arg
         795             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is mixing ratio
         796             :                                                                  ! 2 => tracer is mass (q*dp)
         797             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         798             :      real(r8), intent(in)  :: dp_dry(:,:)                        ! dry pressure level thickness
         799             :      real(r8), intent(in)  :: ptop                               ! pressure at model top
         800             :      real(r8), intent(in)  :: p00                                ! reference pressure for Exner pressure (usually 1000hPa)
         801             :      logical , intent(in)  :: inv_exner                          ! logical for outputting inverse Exner or Exner pressure
         802             :      real(r8), intent(out) :: exner(:,:)
         803             :      real(r8), optional, intent(out) :: poverp0(:,:)             ! for efficiency when a routine needs this variable
         804             :
         805           0 :      real(r8) :: pmid(SIZE(tracer, 1), SIZE(tracer, 2))
         806           0 :      real(r8) :: kappa_dry(SIZE(tracer, 1), SIZE(tracer, 2))
         807             :      character(len=*), parameter :: subname = 'get_exner_1hd: '
         808             :      !
         809             :      ! compute mid level pressure
         810             :      !
         811           0 :      call get_pmid_from_dp(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, pmid)
         812             :      !
         813             :      ! compute kappa = Rd / cpd
         814             :      !
         815           0 :      if (mixing_ratio == DRY_MIXING_RATIO) then
         816           0 :        call get_kappa_dry(tracer, active_species_idx, kappa_dry)
         817           0 :      else if (mixing_ratio == MASS_MIXING_RATIO) then
         818           0 :        call get_kappa_dry(tracer, active_species_idx, kappa_dry, 1.0_r8 / dp_dry)
         819             :      else
         820           0 :        call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
         821             :      end if
         822           0 :      if (inv_exner) then
         823           0 :        exner(:,:) = (p00 / pmid(:,:)) ** kappa_dry(:,:)
         824             :      else
         825           0 :        exner(:,:) = (pmid(:,:) / p00) ** kappa_dry(:,:)
         826             :      end if
         827           0 :      if (present(poverp0)) poverp0 = pmid(:,:) / p00
         828           0 :    end subroutine get_exner_1hd
         829             :
         830             :    !===========================================================================
         831             :
         832             :    !****************************************************************************************************************
         833             :    !
         834             :    ! Compute virtual potential temperature from dp_dry, m, T and ptop.
         835             :    !
         836             :    !****************************************************************************************************************
         837             :    !
         838           0 :    subroutine get_virtual_theta_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, temp, theta_v)
         839             :      real(r8), intent(in)  :: tracer(:,:,:)                      ! tracers; quantity specified by mixing_ratio arg
         840             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is dry mixing ratio
         841             :                                                                  ! 2 => tracer is mass (q*dp)
         842             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         843             :      real(r8), intent(in)  :: dp_dry(:,:)                        ! dry pressure level thickness
         844             :      real(r8), intent(in)  :: ptop                               ! pressure at model top
         845             :      real(r8), intent(in)  :: p00                                ! reference pressure for Exner pressure (usually 1000hPa)
         846             :      real(r8), intent(in)  :: temp(:,:)                          ! temperature
         847             :      real(r8), intent(out) :: theta_v(:,:)                       ! virtual potential temperature
         848             :
         849           0 :      real(r8) :: iexner(SIZE(tracer, 1), SIZE(tracer, 2))
         850             :
         851           0 :      call get_exner(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, .true., iexner)
         852             :
         853           0 :      theta_v(:,:) = temp(:,:) * iexner(:,:)
         854             :
         855           0 :    end subroutine get_virtual_theta_1hd
         856             :
         857             :    !===========================================================================
         858             :
         859             :    !****************************************************************************************************************
         860             :    !
         861             :    ! Compute geopotential from dry pressure level thichkness, water tracers, model top pressure and temperature
         862             :    !
         863             :    !****************************************************************************************************************
         864             :    !
         865           0 :    subroutine get_gz_from_dp_dry_ptop_temp_1hd(tracer, mixing_ratio, active_species_idx, &
         866           0 :                      dp_dry, ptop, temp, phis, gz, pmid, dp, T_v)
         867             :      use air_composition, only: get_R_dry
         868             :      use string_utils,    only: int2str
         869             :      real(r8), intent(in)  :: tracer(:,:,:)                      ! tracer; quantity specified by mixing_ratio arg
         870             :      integer,  intent(in)  :: mixing_ratio                       ! 1 => tracer is dry mixing ratio
         871             :                                                                  ! 2 => tracer is mass (q*dp)
         872             :      integer,  intent(in)  :: active_species_idx(:)              ! index for thermodynamic species in tracer array
         873             :      real(r8), intent(in)  :: dp_dry(:,:)                        ! dry pressure level thickness
         874             :      real(r8), intent(in)  :: ptop                               ! pressure at model top
         875             :      real(r8), intent(in)  :: temp(:,:)                          ! temperature
         876             :      real(r8), intent(in)  :: phis(:)                            ! surface geopotential
         877             :      real(r8), intent(out) :: gz(:,:)                            ! geopotential
         878             :      real(r8), optional, intent(out) :: pmid(:,:)                ! mid-level pressure
         879             :      real(r8), optional, intent(out) :: dp(:,:)                  ! pressure level thickness
         880             :      real(r8), optional, intent(out) :: t_v(:,:)                 ! virtual temperature
         881             :
         882             :
         883           0 :      real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2))     :: pmid_local, t_v_local, dp_local, R_dry
         884           0 :      real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2) + 1) :: pint
         885             :      character(len=*), parameter                               :: subname = 'get_gz_from_dp_dry_ptop_temp_1hd: '
         886             :
         887             :
         888             :      call get_pmid_from_dp(tracer, mixing_ratio, active_species_idx, &
         889           0 :                               dp_dry, ptop, pmid_local, pint=pint, dp=dp_local)
         890           0 :      if (mixing_ratio == DRY_MIXING_RATIO) then
         891           0 :        call get_virtual_temp(tracer, t_v_local, temp=temp, active_species_idx_dycore=active_species_idx)
         892           0 :        call get_R_dry(tracer, active_species_idx, R_dry)
         893           0 :      else if (mixing_ratio == MASS_MIXING_RATIO) then
         894           0 :        call get_virtual_temp(tracer, t_v_local, temp=temp, dp_dry=dp_dry, active_species_idx_dycore=active_species_idx)
         895           0 :        call get_R_dry(tracer,active_species_idx, R_dry, fact=1.0_r8 / dp_dry)
         896             :      else
         897           0 :        call endrun(subname//'unrecognized input ('//int2str(mixing_ratio)//') for mixing_ratio')
         898             :      end if
         899           0 :      call get_gz(dp_local, T_v_local, R_dry, phis, ptop, gz, pmid_local)
         900             :
         901           0 :      if (present(pmid)) pmid=pmid_local
         902           0 :      if (present(T_v))  T_v=T_v_local
         903           0 :      if (present(dp))   dp=dp_local
         904           0 :   end subroutine get_gz_from_dp_dry_ptop_temp_1hd
         905             :
         906             :    !===========================================================================
         907             :
         908             :    !***************************************************************************
         909             :    !
         910             :    ! Compute geopotential from pressure level thickness and virtual temperature
         911             :    !
         912             :    !***************************************************************************
         913             :    !
         914   311688000 :    subroutine get_gz_given_dp_Tv_Rdry_1hd(dp, T_v, R_dry, phis, ptop, gz, pmid)
         915             :      use dycore, only: dycore_is
         916             :      real(r8), intent(in)  :: dp   (:,:)                       ! pressure level thickness
         917             :      real(r8), intent(in)  :: T_v  (:,:)                       ! virtual temperature
         918             :      real(r8), intent(in)  :: R_dry(:,:)                       ! R dry
         919             :      real(r8), intent(in)  :: phis (:)                         ! surface geopotential
         920             :      real(r8), intent(in)  :: ptop                             ! model top presure
         921             :      real(r8), intent(out) :: gz(:,:)                          ! geopotential
         922             :      real(r8), optional, intent(out) :: pmid(:,:)              ! mid-level pressure
         923             :
         924             :
         925   623376000 :      real(r8), dimension(SIZE(dp, 1), SIZE(dp, 2))      :: pmid_local
         926   623376000 :      real(r8), dimension(SIZE(dp, 1), SIZE(dp, 2) + 1)  :: pint
         927   623376000 :      real(r8), dimension(SIZE(dp, 1))                   :: gzh, Rdry_tv
         928             :      integer :: kdx
         929             :
         930   311688000 :      call get_pmid_from_dp(dp, ptop, pmid_local, pint)
         931             :
         932             :      !
         933             :      ! integrate hydrostatic eqn
         934             :      !
         935  1558440000 :      gzh = phis
         936   311688000 :      if (dycore_is('LR') .or. dycore_is('FV3')) then
         937           0 :        do kdx = SIZE(dp, 2), 1, -1
         938           0 :          Rdry_tv(:) = R_dry(:, kdx) * T_v(:, kdx)
         939           0 :          gz(:, kdx) = gzh(:) + Rdry_tv(:) * (1.0_r8 - pint(:, kdx) / pmid_local(:, kdx))
         940           0 :          gzh(:)  = gzh(:) + Rdry_tv(:) * (log(pint(:, kdx + 1)) - log(pint(:, kdx)))
         941             :        end do
         942             :      else
         943  8415576000 :        do kdx = SIZE(dp,2), 1, -1
         944 40519440000 :          Rdry_tv(:) = R_dry(:,kdx) * T_v(:, kdx)
         945 40519440000 :          gz(:,kdx) = gzh(:) + Rdry_tv(:) * 0.5_r8 * dp(:, kdx) / pmid_local(:, kdx)
         946 40831128000 :          gzh(:)  = gzh(:) + Rdry_tv(:) * dp(:, kdx) / pmid_local(:, kdx)
         947             :        end do
         948             :      end if
         949 40831128000 :      if (present(pmid)) pmid=pmid_local
         950   311688000 :    end subroutine get_gz_given_dp_Tv_Rdry_1hd
         951             :
         952    77922000 :    subroutine get_gz_given_dp_Tv_Rdry_2hd(dp, T_v, R_dry, phis, ptop, gz, pmid)
         953             :      ! Version of get_gz_given_dp_Tv_Rdry for arrays that have a second horizontal index
         954             :      real(r8), intent(in)  :: dp   (:,:,:)                     ! pressure level thickness
         955             :      real(r8), intent(in)  :: T_v  (:,:,:)                     ! virtual temperature
         956             :      real(r8), intent(in)  :: R_dry(:,:,:)                     ! R dry
         957             :      real(r8), intent(in)  :: phis (:,:)                       ! surface geopotential
         958             :      real(r8), intent(in)  :: ptop                             ! model top presure
         959             :      real(r8), intent(out) :: gz(:,:,:)                        ! geopotential
         960             :      real(r8), optional, intent(out) :: pmid(:,:,:)            ! mid-level pressure
         961             :
         962             :      integer :: jdx
         963             :
         964   389610000 :      do jdx = 1, SIZE(dp, 2)
         965   389610000 :        if (present(pmid)) then
         966             :          call get_gz(dp(:, jdx, :), T_v(:, jdx, :), R_dry(:, jdx, :), phis(:, jdx), &
         967   311688000 :               ptop, gz(:, jdx, :), pmid=pmid(:, jdx, :))
         968             :        else
         969           0 :          call get_gz(dp(:, jdx, :), T_v(:, jdx, :), R_dry(:, jdx, :), phis(:, jdx), ptop, gz(:, jdx, :))
         970             :        end if
         971             :      end do
         972             :
         973             :
         974    77922000 :    end subroutine get_gz_given_dp_Tv_Rdry_2hd
         975             :
         976             :    !===========================================================================
         977             :
         978             :    !***************************************************************************
         979             :    !
         980             :    ! Compute Richardson number at cell interfaces (half levels)
         981             :    !
         982             :    !***************************************************************************
         983             :    !
         984           0 :    subroutine get_Richardson_number_1hd(tracer,mixing_ratio, active_species_idx, dp_dry, ptop, &
         985           0 :         p00, temp, v, Richardson_number, pmid, dp)
         986             :      real(r8), intent(in)  :: tracer(:,:,:)                        ! tracer; quantity specified by mixing_ratio arg
         987             :      integer,  intent(in)  :: mixing_ratio                         ! 1 => tracer is dry mixing ratio
         988             :                                                                    ! 2 => tracer is mass (q*dp)
         989             :      integer,  intent(in)  :: active_species_idx(:)                ! index for thermodynamic species in tracer array
         990             :      real(r8), intent(in)  :: dp_dry(:,:)                          ! dry pressure level thickness
         991             :      real(r8), intent(in)  :: ptop                                 ! pressure at model top
         992             :      real(r8), intent(in)  :: p00                                  ! reference pressure for Exner pressure (usually 1000hPa)
         993             :      real(r8), intent(in)  :: temp(:,:)                            ! temperature
         994             :      real(r8), intent(in)  :: v(:,:,:)                             ! velocity components
         995             :      real(r8), intent(out) :: Richardson_number(:,:)
         996             :      real(r8), optional, intent(out) :: pmid(:,:)
         997             :      real(r8), optional, intent(out) :: dp(:,:)
         998             :
         999           0 :      real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2)) :: gz, theta_v
        1000           0 :      real(r8), dimension(SIZE(tracer, 1))                  :: pt1, pt2, phis
        1001             :      integer :: kdx, kdxm1
        1002             :      real(r8), parameter:: ustar2 = 1.E-4_r8
        1003             :
        1004           0 :      phis = 0.0_r8
        1005           0 :      call get_gz(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, temp, phis, gz, pmid=pmid, dp=dp)
        1006           0 :      call get_virtual_theta(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, p00, temp, theta_v)
        1007           0 :      Richardson_number(:, 1)                   = 0.0_r8
        1008           0 :      Richardson_number(:, SIZE(tracer, 2) + 1) = 0.0_r8
        1009           0 :      do kdx = SIZE(tracer, 2), 2, -1
        1010           0 :        kdxm1 = kdx - 1
        1011           0 :        pt1(:) = theta_v(:, kdxm1)
        1012           0 :        pt2(:) = theta_v(:, kdx)
        1013           0 :        Richardson_number(:, kdx) = (gz(:, kdxm1) - gz(:, kdx)) * (pt1 - pt2) / ( 0.5_r8*(pt1 + pt2) *        &
        1014           0 :             ((v(:, 1, kdxm1) - v(:, 1, kdx)) ** 2 + (v(:, 2, kdxm1) - v(:, 2, kdx)) ** 2 + ustar2) )
        1015             :      end do
        1016           0 :    end subroutine get_Richardson_number_1hd
        1017             :
        1018             :    !
        1019             :    !****************************************************************************************************************
        1020             :    !
        1021             :    ! get surface pressure from dry pressure and thermodynamic active species (e.g., forms of water: water vapor, cldliq, etc.)
        1022             :    !
        1023             :    !****************************************************************************************************************
        1024             :    !
        1025           0 :    subroutine get_ps_1hd(tracer_mass, active_species_idx, dp_dry, ps, ptop)
        1026             :      use air_composition,  only: dry_air_species_num
        1027             :
        1028             :      real(r8), intent(in)   :: tracer_mass(:,:,:)                      ! Tracer array (q*dp)
        1029             :      real(r8), intent(in)   :: dp_dry(:,:)                             ! dry pressure level thickness
        1030             :      real(r8), intent(out)  :: ps(:)                                   ! surface pressure
        1031             :      real(r8), intent(in)   :: ptop
        1032             :      integer,  intent(in)   :: active_species_idx(:)
        1033             :
        1034             :      integer                    :: idx, kdx, m_cnst, qdx
        1035           0 :      real(r8)                   :: dp(SIZE(tracer_mass, 1), SIZE(tracer_mass, 2))  ! dry pressure level thickness
        1036             :
        1037           0 :      dp = dp_dry
        1038           0 :      do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
        1039           0 :        m_cnst = active_species_idx(qdx)
        1040           0 :        do kdx = 1, SIZE(tracer_mass, 2)
        1041           0 :          do idx = 1, SIZE(tracer_mass, 1)
        1042           0 :            dp(idx, kdx) = dp(idx, kdx) + tracer_mass(idx, kdx, m_cnst)
        1043             :          end do
        1044             :        end do
        1045             :      end do
        1046           0 :      ps = ptop
        1047           0 :      do kdx = 1, SIZE(tracer_mass, 2)
        1048           0 :        do idx = 1, SIZE(tracer_mass, 1)
        1049           0 :          ps(idx) = ps(idx) + dp(idx, kdx)
        1050             :        end do
        1051             :      end do
        1052           0 :    end subroutine get_ps_1hd
        1053             :
        1054           0 :    subroutine get_ps_2hd(tracer_mass, active_species_idx, dp_dry, ps, ptop)
        1055             :      ! Version of get_ps for arrays that have a second horizontal index
        1056             :      real(r8), intent(in)   :: tracer_mass(:,:,:,:)                      ! Tracer array (q*dp)
        1057             :      real(r8), intent(in)   :: dp_dry(:,:,:)                             ! dry pressure level thickness
        1058             :      real(r8), intent(out)  :: ps(:,:)                                   ! surface pressure
        1059             :      real(r8), intent(in)   :: ptop
        1060             :      integer,  intent(in)   :: active_species_idx(:)
        1061             :
        1062             :      integer :: jdx
        1063             :
        1064           0 :      do jdx = 1, SIZE(tracer_mass, 2)
        1065           0 :        call get_ps(tracer_mass(:, jdx, :, :), active_species_idx, dp_dry(:, jdx, :), ps(:, jdx), ptop)
        1066             :      end do
        1067             :
        1068           0 :    end subroutine get_ps_2hd
        1069             :
        1070             :    !===========================================================================
        1071             :
        1072             :    !*************************************************************************************************************************
        1073             :    !
        1074             :    ! compute generalized kappa =Rdry/cpdry
        1075             :    !
        1076             :    !*************************************************************************************************************************
        1077             :    !
        1078    62337600 :    subroutine get_kappa_dry_1hd(tracer, active_species_idx, kappa_dry, fact)
        1079             :      use air_composition,  only: dry_air_species_num, get_R_dry, get_cp_dry
        1080             :      use physconst,        only: rair, cpair
        1081             :
        1082             :      real(r8), intent(in)  :: tracer(:,:,:)              !tracer array
        1083             :      integer,  intent(in)  :: active_species_idx(:)      !index of thermodynamic active tracers
        1084             :      real(r8), intent(out) :: kappa_dry(:,:)             !kappa dry
        1085             :      real(r8), optional, intent(in) :: fact(:,:)         !factor for converting tracer to dry mixing ratio
        1086             :      !
        1087    62337600 :      real(r8), allocatable, dimension(:,:) :: cp_dry,R_dry
        1088             :      integer                     :: ierr
        1089             :      character(len=*), parameter :: subname = "get_kappa_dry_1hd"
        1090             :      character(len=*), parameter :: errstr = subname//": failed to allocate "
        1091             :      !
        1092             :      ! dry air not species dependent
        1093    62337600 :      if (dry_air_species_num==0) then
        1094  8166225600 :        kappa_dry = rair / cpair
        1095             :      else
        1096           0 :        allocate(R_dry(SIZE(kappa_dry, 1), SIZE(kappa_dry, 2)), stat=ierr)
        1097           0 :        if (ierr /= 0) then
        1098           0 :          call endrun(errstr//"R_dry")
        1099             :        end if
        1100           0 :        allocate(cp_dry(SIZE(kappa_dry, 1), SIZE(kappa_dry, 2)), stat=ierr)
        1101           0 :        if (ierr /= 0) then
        1102           0 :          call endrun(errstr//"cp_dry")
        1103             :        end if
        1104           0 :        call get_cp_dry(tracer, active_species_idx, cp_dry, fact=fact)
        1105           0 :        call get_R_dry( tracer, active_species_idx, R_dry,  fact=fact)
        1106           0 :        kappa_dry = R_dry / cp_dry
        1107           0 :        deallocate(R_dry, cp_dry)
        1108             :      end if
        1109    62337600 :    end subroutine get_kappa_dry_1hd
        1110             :
        1111    15584400 :    subroutine get_kappa_dry_2hd(tracer, active_species_idx, kappa_dry, fact)
        1112             :      ! Version of get_kappa_dry for arrays that have a second horizontal index
        1113             :      real(r8), intent(in)  :: tracer(:,:,:,:)              !tracer array
        1114             :      integer,  intent(in)  :: active_species_idx(:)        !index of thermodynamic active tracers
        1115             :      real(r8), intent(out) :: kappa_dry(:,:,:)             !kappa dry
        1116             :      real(r8), optional, intent(in) :: fact(:,:,:)         !factor for converting tracer to dry mixing ratio
        1117             :
        1118             :      integer :: jdx
        1119             :
        1120    77922000 :      do jdx = 1, SIZE(tracer, 2)
        1121    77922000 :        if (present(fact)) then
        1122           0 :          call get_kappa_dry(tracer(:, jdx, :, :), active_species_idx, kappa_dry(:, jdx, :), fact=fact(:, jdx, :))
        1123             :        else
        1124    62337600 :          call get_kappa_dry(tracer(:, jdx, :, :), active_species_idx, kappa_dry(:, jdx, :))
        1125             :        end if
        1126             :      end do
        1127             :
        1128    15584400 :    end subroutine get_kappa_dry_2hd
        1129             :
        1130             :    !===========================================================================
        1131             :
        1132             :    !*************************************************************************************************************************
        1133             :    !
        1134             :    ! compute reference pressure levels
        1135             :    !
        1136             :    !*************************************************************************************************************************
        1137             :    !
        1138       43200 :    subroutine get_dp_ref_1hd(hyai, hybi, ps0, phis, dp_ref, ps_ref)
        1139             :      use physconst,  only: tref, rair
        1140             :      real(r8), intent(in)           :: hyai(:)
        1141             :      real(r8), intent(in)           :: hybi(:)
        1142             :      real(r8), intent(in)           :: ps0
        1143             :      real(r8), intent(in)           :: phis(:)
        1144             :      real(r8), intent(out)          :: dp_ref(:,:)
        1145             :      real(r8), intent(out)          :: ps_ref(:)
        1146             :      integer :: kdx
        1147             :      !
        1148             :      ! use static reference pressure (hydrostatic balance incl. effect of topography)
        1149             :      !
        1150      216000 :      ps_ref(:) = ps0 * exp(-phis(:) / (rair * tref))
        1151     1166400 :      do kdx = 1, SIZE(dp_ref, 2)
        1152     5659200 :        dp_ref(:,kdx) = ((hyai(kdx + 1) - hyai(kdx)) * ps0 + (hybi(kdx + 1) - hybi(kdx)) * ps_ref(:))
        1153             :      end do
        1154       43200 :    end subroutine get_dp_ref_1hd
        1155             :
        1156       10800 :    subroutine get_dp_ref_2hd(hyai, hybi, ps0, phis, dp_ref, ps_ref)
        1157             :      ! Version of get_dp_ref for arrays that have a second horizontal index
        1158             :      real(r8), intent(in)           :: hyai(:)
        1159             :      real(r8), intent(in)           :: hybi(:)
        1160             :      real(r8), intent(in)           :: ps0
        1161             :      real(r8), intent(in)           :: phis(:,:)
        1162             :      real(r8), intent(out)          :: dp_ref(:,:,:)
        1163             :      real(r8), intent(out)          :: ps_ref(:,:)
        1164             :      integer :: jdx
        1165             :
        1166       54000 :      do jdx = 1, SIZE(dp_ref, 2)
        1167       54000 :        call get_dp_ref(hyai, hybi, ps0, phis(:, jdx), dp_ref(:, jdx, :), ps_ref(:, jdx))
        1168             :      end do
        1169             :
        1170       10800 :    end subroutine get_dp_ref_2hd
        1171             :
        1172             :    !===========================================================================
        1173             :
        1174             :    !*************************************************************************************************************************
        1175             :    !
        1176             :    ! compute dry densisty from temperature (temp) and pressure (dp_dry and tracer)
        1177             :    !
        1178             :    !*************************************************************************************************************************
        1179             :    !
        1180           0 :    subroutine get_rho_dry_1hd(tracer, temp, ptop, dp_dry, tracer_mass, rho_dry, rhoi_dry, &
        1181           0 :               active_species_idx_dycore)
        1182             :      use air_composition, only: get_R_dry
        1183             :      ! args
        1184             :      real(r8), intent(in)           :: tracer(:,:,:)      ! Tracer array
        1185             :      real(r8), intent(in)           :: temp(:,:)          ! Temperature
        1186             :      real(r8), intent(in)           :: ptop
        1187             :      real(r8), intent(in)           :: dp_dry(:,:)
        1188             :      logical,  intent(in)           :: tracer_mass
        1189             :      real(r8), optional,intent(out) :: rho_dry(:,:)
        1190             :      real(r8), optional,intent(out) :: rhoi_dry(:,:)
        1191             :      !
        1192             :      ! array of indicies for index of thermodynamic active species in dycore tracer array
        1193             :      ! (if different from physics index)
        1194             :      !
        1195             :      integer, optional, intent(in)   :: active_species_idx_dycore(:)
        1196             :
        1197             :      ! local vars
        1198             :      integer :: idx, kdx
        1199           0 :      real(r8),  dimension(SIZE(tracer, 1), SIZE(tracer, 2))        :: pmid
        1200           0 :      real(r8),  dimension(SIZE(tracer, 1), SIZE(tracer, 2) + 1)    :: pint
        1201           0 :      real(r8),  allocatable                                        :: R_dry(:,:)
        1202           0 :      integer,  dimension(thermodynamic_active_species_num)         :: idx_local
        1203             :      integer                     :: ierr
        1204             :      character(len=*), parameter :: subname = "get_rho_dry_1hd"
        1205             :      character(len=*), parameter :: errstr = subname//": failed to allocate "
        1206             :
        1207           0 :      if (present(active_species_idx_dycore)) then
        1208           0 :        idx_local = active_species_idx_dycore
        1209             :      else
        1210           0 :        idx_local = thermodynamic_active_species_idx
        1211             :      end if
        1212             :      !
        1213             :      ! we assume that air is dry where molecular viscosity may be significant
        1214             :      !
        1215           0 :      call get_pmid_from_dp(dp_dry, ptop, pmid, pint=pint)
        1216           0 :      if (present(rhoi_dry)) then
        1217           0 :        allocate(R_dry(SIZE(tracer, 1), SIZE(tracer, 2) + 1), stat=ierr)
        1218           0 :        if (ierr /= 0) then
        1219           0 :          call endrun(errstr//"R_dry")
        1220             :        end if
        1221           0 :        if (tracer_mass) then
        1222           0 :          call get_R_dry(tracer, idx_local, R_dry, fact=1.0_r8 / dp_dry)
        1223             :        else
        1224           0 :          call get_R_dry(tracer, idx_local, R_dry)
        1225             :        end if
        1226           0 :        do kdx = 2, SIZE(tracer, 2) + 1
        1227           0 :          rhoi_dry(:, kdx) = 0.5_r8 * (temp(:, kdx) + temp(:, kdx - 1))!could be more accurate!
        1228           0 :          rhoi_dry(:, kdx) = pint(:,kdx) / (rhoi_dry(:, kdx) * R_dry(:, kdx)) !ideal gas law for dry air
        1229             :        end do
        1230             :        !
        1231             :        ! extrapolate top level value
        1232             :        !
        1233           0 :        kdx=1
        1234           0 :        rhoi_dry(:, kdx) = 1.5_r8 * (temp(:, kdx) - 0.5_r8 * temp(:, kdx + 1))
        1235           0 :        rhoi_dry(:, kdx) = pint(:, kdx) / (rhoi_dry(:, kdx) * R_dry(:, kdx)) !ideal gas law for dry air
        1236           0 :        deallocate(R_dry)
        1237             :      end if
        1238           0 :      if (present(rho_dry)) then
        1239           0 :        allocate(R_dry(SIZE(tracer, 1), size(rho_dry, 2)), stat=ierr)
        1240           0 :        if (ierr /= 0) then
        1241           0 :          call endrun(errstr//"R_dry")
        1242             :        end if
        1243           0 :        if (tracer_mass) then
        1244           0 :          call get_R_dry(tracer, idx_local, R_dry, fact=1.0_r8 / dp_dry)
        1245             :        else
        1246           0 :          call get_R_dry(tracer, idx_local, R_dry)
        1247             :        end if
        1248           0 :        do kdx = 1, SIZE(rho_dry, 2)
        1249           0 :          do idx = 1, SIZE(rho_dry, 1)
        1250           0 :            rho_dry(idx, kdx) = pmid(idx, kdx) / (temp(idx, kdx) * R_dry(idx, kdx)) !ideal gas law for dry air
        1251             :          end do
        1252             :        end do
        1253           0 :        deallocate(R_dry)
        1254             :      end if
        1255           0 :    end subroutine get_rho_dry_1hd
        1256             :
        1257           0 :    subroutine get_rho_dry_2hd(tracer, temp, ptop, dp_dry, tracer_mass, rho_dry, rhoi_dry, &
        1258           0 :               active_species_idx_dycore)
        1259             :      ! Version of get_rho_dry for arrays that have a second horizontal index
        1260             :      real(r8), intent(in)           :: tracer(:,:,:,:)      ! Tracer array
        1261             :      real(r8), intent(in)           :: temp(:,:,:)          ! Temperature
        1262             :      real(r8), intent(in)           :: ptop
        1263             :      real(r8), intent(in)           :: dp_dry(:,:,:)
        1264             :      logical,  intent(in)           :: tracer_mass
        1265             :      real(r8), optional,intent(out) :: rho_dry(:,:,:)
        1266             :      real(r8), optional,intent(out) :: rhoi_dry(:,:,:)
        1267             :      !
        1268             :      ! array of indicies for index of thermodynamic active species in dycore tracer array
        1269             :      ! (if different from physics index)
        1270             :      !
        1271             :      integer, optional, intent(in)   :: active_species_idx_dycore(:)
        1272             :
        1273             :      integer :: jdx
        1274             :
        1275           0 :      do jdx = 1, SIZE(tracer, 2)
        1276           0 :        if (present(rho_dry) .and. present(rhoi_dry)) then
        1277             :           call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
        1278             :               tracer_mass, rho_dry=rho_dry(:, jdx, :), rhoi_dry=rhoi_dry(:, jdx, :), &
        1279           0 :               active_species_idx_dycore=active_species_idx_dycore)
        1280           0 :        else if (present(rho_dry)) then
        1281             :           call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
        1282           0 :               tracer_mass, rho_dry=rho_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
        1283           0 :        else if (present(rhoi_dry)) then
        1284             :           call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), &
        1285           0 :               tracer_mass, rhoi_dry=rhoi_dry(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
        1286             :        else
        1287             :           call get_rho_dry(tracer(:, jdx, :, :), temp(:, jdx, :), ptop, dp_dry(:, jdx, :), tracer_mass, &
        1288           0 :               active_species_idx_dycore=active_species_idx_dycore)
        1289             :        end if
        1290             :      end do
        1291             :
        1292           0 :    end subroutine get_rho_dry_2hd
        1293             :    !===========================================================================
        1294             :
        1295             :    !*************************************************************************************************************************
        1296             :    !
        1297             :    ! compute 3D molecular diffusion and thermal conductivity
        1298             :    !
        1299             :    !*************************************************************************************************************************
        1300             :    !
        1301           0 :    subroutine get_molecular_diff_coef_1hd(temp, get_at_interfaces, sponge_factor, kmvis, kmcnd, &
        1302           0 :         tracer, fact, active_species_idx_dycore, mbarv_in)
        1303             :      use air_composition,  only: dry_air_species_num, get_mbarv
        1304             :      use air_composition,  only: kv1, kc1, kv2, kc2, kv_temp_exp, kc_temp_exp
        1305             :
        1306             :      ! args
        1307             :      real(r8), intent(in)           :: temp(:,:)                     ! temperature
        1308             :      logical,  intent(in)           :: get_at_interfaces             ! true: compute kmvis and kmcnd at interfaces
        1309             :                                                                      ! false: compute kmvis and kmcnd at mid-levels
        1310             :      real(r8), intent(in)           :: sponge_factor(:)              ! multiply kmvis and kmcnd with sponge_factor
        1311             :                                                                      ! (for sponge layer)
        1312             :      real(r8), intent(out)          :: kmvis(:,:)
        1313             :      real(r8), intent(out)          :: kmcnd(:,:)
        1314             :      real(r8), intent(in)           :: tracer(:,:,:)                 ! tracer array
        1315             :      integer,  intent(in), optional :: active_species_idx_dycore(:)  ! index of active species in tracer
        1316             :      real(r8), intent(in), optional :: fact(:,:)                     ! if tracer is in units of mass or moist
        1317             :                                                                      ! fact converts to dry mixing ratio: tracer/fact
        1318             :      real(r8), intent(in), optional :: mbarv_in(:,:)                 ! composition dependent atmosphere mean mass
        1319             :      !
        1320             :      ! local vars
        1321             :      !
        1322             :      integer :: idx, kdx, icnst, ispecies
        1323             :      real(r8):: mbarvi, mm, residual             ! Mean mass at mid level
        1324             :      real(r8):: cnst_vis, cnst_cnd, temp_local
        1325           0 :      real(r8), dimension(SIZE(tracer,1), SIZE(sponge_factor, 1)) :: factor, mbarv
        1326           0 :      integer,  dimension(thermodynamic_active_species_num)       :: idx_local
        1327             :      character(len=*), parameter :: subname = 'get_molecular_diff_coef_1hd: '
        1328             :
        1329             :      !--------------------------------------------
        1330             :      ! Set constants needed for updates
        1331             :      !--------------------------------------------
        1332             :
        1333           0 :      if (dry_air_species_num==0) then
        1334             :
        1335           0 :        cnst_vis = (kv1 * mmro2 * o2_mwi + kv2 * mmrn2 * n2_mwi) * mbar
        1336           0 :        cnst_cnd = (kc1 * mmro2 * o2_mwi + kc2 * mmrn2 * n2_mwi) * mbar
        1337           0 :        if (get_at_interfaces) then
        1338           0 :            do kdx = 2, SIZE(sponge_factor, 1)
        1339           0 :              do idx = 1, SIZE(tracer, 1)
        1340           0 :                temp_local   = 0.5_r8 * (temp(idx, kdx) + temp(idx, kdx - 1))
        1341           0 :                kmvis(idx, kdx) = sponge_factor(kdx) * cnst_vis * temp_local ** kv_temp_exp
        1342           0 :                kmcnd(idx, kdx) = sponge_factor(kdx) * cnst_cnd * temp_local ** kc_temp_exp
        1343             :              end do
        1344             :            end do
        1345             :            !
        1346             :            ! extrapolate top level value
        1347             :            !
        1348           0 :            kmvis(1:SIZE(tracer, 1), 1) = 1.5_r8 * kmvis(1:SIZE(tracer, 1), 2) - 0.5_r8 * kmvis(1:SIZE(tracer, 1), 3)
        1349           0 :            kmcnd(1:SIZE(tracer, 1), 1) = 1.5_r8 * kmcnd(1:SIZE(tracer, 1), 2) - 0.5_r8 * kmcnd(1:SIZE(tracer, 1), 3)
        1350             :        else if (.not. get_at_interfaces) then
        1351           0 :          do kdx = 1, SIZE(sponge_factor, 1)
        1352           0 :            do idx = 1, SIZE(tracer, 1)
        1353           0 :              kmvis(idx, kdx) = sponge_factor(kdx) * cnst_vis * temp(idx, kdx) ** kv_temp_exp
        1354           0 :              kmcnd(idx, kdx) = sponge_factor(kdx) * cnst_cnd * temp(idx, kdx) ** kc_temp_exp
        1355             :            end do
        1356             :          end do
        1357             :        else
        1358             :          call endrun(subname//'get_at_interfaces must be .true. or .false.')
        1359             :        end if
        1360             :      else
        1361           0 :        if (present(active_species_idx_dycore)) then
        1362           0 :          idx_local = active_species_idx_dycore
        1363             :        else
        1364           0 :          idx_local = thermodynamic_active_species_idx
        1365             :        end if
        1366           0 :        if (present(fact)) then
        1367           0 :          factor = fact(:,:)
        1368             :        else
        1369           0 :          factor = 1.0_r8
        1370             :        endif
        1371           0 :        if (present(mbarv_in)) then
        1372           0 :          mbarv = mbarv_in
        1373             :        else
        1374           0 :          call get_mbarv(tracer, idx_local, mbarv, fact=factor)
        1375             :        end if
        1376             :        !
        1377             :        ! major species dependent code
        1378             :        !
        1379           0 :        if (get_at_interfaces) then
        1380           0 :          do kdx = 2, SIZE(sponge_factor, 1)
        1381           0 :            do idx = 1, SIZE(tracer, 1)
        1382           0 :              kmvis(idx, kdx) = 0.0_r8
        1383           0 :              kmcnd(idx, kdx) = 0.0_r8
        1384           0 :              residual = 1.0_r8
        1385           0 :              do icnst = 1, dry_air_species_num
        1386           0 :                ispecies = idx_local(icnst)
        1387           0 :                mm       = 0.5_r8 * (tracer(idx, kdx, ispecies) * factor(idx, kdx) + &
        1388           0 :                           tracer(idx, kdx - 1, ispecies) * factor(idx, kdx-1))
        1389           0 :                kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
        1390           0 :                                            thermodynamic_active_species_mwi(icnst) * mm
        1391           0 :                kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
        1392           0 :                                            thermodynamic_active_species_mwi(icnst) * mm
        1393           0 :                residual        = residual - mm
        1394             :              end do
        1395           0 :              icnst = 0 ! N2
        1396           0 :              kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
        1397           0 :                                          thermodynamic_active_species_mwi(icnst) * residual
        1398           0 :              kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
        1399           0 :                                          thermodynamic_active_species_mwi(icnst) * residual
        1400             :
        1401           0 :              temp_local = 0.5_r8 * (temp(idx, kdx - 1) + temp(idx, kdx))
        1402           0 :              mbarvi = 0.5_r8 * (mbarv(idx, kdx - 1) + mbarv(idx, kdx))
        1403           0 :              kmvis(idx, kdx) = kmvis(idx, kdx) * mbarvi * temp_local ** kv_temp_exp
        1404           0 :              kmcnd(idx, kdx) = kmcnd(idx, kdx) * mbarvi * temp_local ** kc_temp_exp
        1405             :            enddo
        1406             :          end do
        1407           0 :          do idx = 1, SIZE(tracer, 1)
        1408           0 :            kmvis(idx, 1) = 1.5_r8 * kmvis(idx, 2) - .5_r8 * kmvis(idx, 3)
        1409           0 :            kmcnd(idx, 1) = 1.5_r8 * kmcnd(idx, 2) - .5_r8 * kmcnd(idx, 3)
        1410           0 :            kmvis(idx, SIZE(sponge_factor, 1) + 1) = kmvis(idx, SIZE(sponge_factor, 1))
        1411           0 :            kmcnd(idx, SIZE(sponge_factor, 1) + 1) = kmcnd(idx, SIZE(sponge_factor, 1))
        1412             :          end do
        1413             :        else if (.not. get_at_interfaces) then
        1414           0 :          do kdx = 1, SIZE(sponge_factor, 1)
        1415           0 :            do idx = 1, SIZE(tracer, 1)
        1416           0 :              kmvis(idx, kdx) = 0.0_r8
        1417           0 :              kmcnd(idx, kdx) = 0.0_r8
        1418           0 :              residual = 1.0_r8
        1419           0 :              do icnst = 1, dry_air_species_num - 1
        1420           0 :                ispecies = idx_local(icnst)
        1421           0 :                mm = tracer(idx, kdx, ispecies) * factor(idx, kdx)
        1422           0 :                kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
        1423           0 :                                  thermodynamic_active_species_mwi(icnst) * mm
        1424           0 :                kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
        1425           0 :                                  thermodynamic_active_species_mwi(icnst) * mm
        1426           0 :                residual        = residual - mm
        1427             :              end do
        1428           0 :              icnst = dry_air_species_num
        1429           0 :              kmvis(idx, kdx) = kmvis(idx, kdx) + thermodynamic_active_species_kv(icnst) * &
        1430           0 :                                thermodynamic_active_species_mwi(icnst) * residual
        1431           0 :              kmcnd(idx, kdx) = kmcnd(idx, kdx) + thermodynamic_active_species_kc(icnst) * &
        1432           0 :                                thermodynamic_active_species_mwi(icnst) * residual
        1433             :
        1434           0 :              kmvis(idx, kdx) = kmvis(idx, kdx) * mbarv(idx, kdx) * temp(idx, kdx) ** kv_temp_exp
        1435           0 :              kmcnd(idx, kdx) = kmcnd(idx, kdx) * mbarv(idx, kdx) * temp(idx, kdx) ** kc_temp_exp
        1436             :            end do
        1437             :          end do
        1438             :        else
        1439             :          call endrun(subname//'get_at_interfaces must be .true. or .false.')
        1440             :        end if
        1441             :      end if
        1442           0 :    end subroutine get_molecular_diff_coef_1hd
        1443             :
        1444           0 :    subroutine get_molecular_diff_coef_2hd(temp, get_at_interfaces, sponge_factor, kmvis, kmcnd, &
        1445           0 :         tracer, fact, active_species_idx_dycore, mbarv_in)
        1446             :      ! Version of get_molecular_diff_coef for arrays that have a second horizontal index
        1447             :      real(r8), intent(in)           :: temp(:,:,:)                     ! temperature
        1448             :      logical,  intent(in)           :: get_at_interfaces               ! true: compute kmvis and kmcnd at interfaces
        1449             :                                                                        ! false: compute kmvis and kmcnd at mid-levels
        1450             :      real(r8), intent(in)           :: sponge_factor(:)                ! multiply kmvis and kmcnd with sponge_factor
        1451             :                                                                        ! (for sponge layer)
        1452             :      real(r8), intent(out)          :: kmvis(:,:,:)
        1453             :      real(r8), intent(out)          :: kmcnd(:,:,:)
        1454             :      real(r8), intent(in)           :: tracer(:,:,:,:)                 ! tracer array
        1455             :      integer,  intent(in), optional :: active_species_idx_dycore(:)    ! index of active species in tracer
        1456             :      real(r8), intent(in), optional :: fact(:,:,:)                     ! if tracer is in units of mass or moist
        1457             :                                                                        ! fact converts to dry mixing ratio: tracer/fact
        1458             :      real(r8), intent(in), optional :: mbarv_in(:,:,:)                 ! composition dependent atmosphere mean mass
        1459             :      integer :: jdx
        1460             :
        1461           0 :      do jdx = 1, SIZE(tracer, 2)
        1462           0 :        if (present(fact) .and. present(mbarv_in)) then
        1463             :          call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
        1464             :               kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), fact=fact(:, jdx, :), &
        1465           0 :               active_species_idx_dycore=active_species_idx_dycore, mbarv_in=mbarv_in(:, jdx, :))
        1466           0 :        else if (present(fact)) then
        1467             :          call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
        1468             :               kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), fact=fact(:, jdx, :), &
        1469           0 :               active_species_idx_dycore=active_species_idx_dycore)
        1470           0 :        else if (present(mbarv_in)) then
        1471             :          call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
        1472             :               kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), &
        1473           0 :               active_species_idx_dycore=active_species_idx_dycore, mbarv_in=mbarv_in(:, jdx, :))
        1474             :        else
        1475             :          call get_molecular_diff_coef(temp(:, jdx, :), get_at_interfaces, sponge_factor, &
        1476             :               kmvis(:, jdx, :), kmcnd(:, jdx, :), tracer(:, jdx, :, :), &
        1477           0 :               active_species_idx_dycore=active_species_idx_dycore)
        1478             :        end if
        1479             :      end do
        1480             :
        1481           0 :    end subroutine get_molecular_diff_coef_2hd
        1482             :    !===========================================================================
        1483             :
        1484             :    !***************************************************************************
        1485             :    !
        1486             :    ! compute reference vertical profile of density, molecular diffusion and thermal conductivity
        1487             :    !
        1488             :    !***************************************************************************
        1489             :    !
        1490           0 :    subroutine get_molecular_diff_coef_reference(tref,press,sponge_factor,kmvis_ref,kmcnd_ref,rho_ref)
        1491             :      use physconst,       only: rair
        1492             :      use air_composition, only: kv1, kv2, kc1, kc2, kv_temp_exp, kc_temp_exp
        1493             :      ! args
        1494             :      real(r8), intent(in)           :: tref                 !reference temperature
        1495             :      real(r8), intent(in)           :: press(:)             !pressure
        1496             :      real(r8), intent(in)           :: sponge_factor(:)     !multiply kmvis and kmcnd with sponge_factor (for sponge layer)
        1497             :      real(r8), intent(out)          :: kmvis_ref(:)         !reference molecular diffusion coefficient
        1498             :      real(r8), intent(out)          :: kmcnd_ref(:)         !reference thermal conductivity coefficient
        1499             :      real(r8), intent(out)          :: rho_ref(:)           !reference density
        1500             :
        1501             :      ! local vars
        1502             :      integer :: kdx
        1503             :
        1504             :      !--------------------------------------------
        1505             :      ! Set constants needed for updates
        1506             :      !--------------------------------------------
        1507             :
        1508           0 :      do kdx = 1, SIZE(press, 1)
        1509           0 :        rho_ref(kdx) = press(kdx) / (tref * rair) !ideal gas law for dry air
        1510           0 :        kmvis_ref(kdx) = sponge_factor(kdx) * &
        1511             :             (kv1 * mmro2 * o2_mwi +         &
        1512             :              kv2 * mmrn2 * n2_mwi) * mbar *    &
        1513           0 :              tref ** kv_temp_exp
        1514           0 :        kmcnd_ref(kdx) = sponge_factor(kdx) * &
        1515             :             (kc1 * mmro2 * o2_mwi +         &
        1516             :              kc2 * mmrn2 * n2_mwi) * mbar *    &
        1517           0 :              tref ** kc_temp_exp
        1518             :      end do
        1519           0 :    end subroutine get_molecular_diff_coef_reference
        1520             :
        1521             :    !==========================================================================
        1522             :
        1523             :    !
        1524             :    !***************************************************************************
        1525             :    !
        1526             :    ! cam_thermo_calc_kappav: update species dependent kappa for FV dycore
        1527             :    !
        1528             :    !***************************************************************************
        1529             :    !
        1530           0 :    subroutine cam_thermo_calc_kappav_2hd(tracer, kappav, cpv)
        1531             :       use air_composition, only: get_R_dry, get_cp_dry
        1532             :       ! assumes moist MMRs
        1533             :
        1534             :       ! Dummy arguments
        1535             :       real(r8),           intent(in)  :: tracer(:, :, :, :)
        1536             :       real(r8),           intent(out) :: kappav(:, :, :)
        1537             :       real(r8), optional, intent(out) :: cpv(:, :, :)
        1538             :
        1539             :       ! Local variables
        1540           0 :       real(r8) :: rgas_var(SIZE(tracer, 1), SIZE(tracer, 2), SIZE(tracer, 3))
        1541           0 :       real(r8) :: cp_var(SIZE(tracer, 1), SIZE(tracer, 2), SIZE(tracer, 3))
        1542             :       integer :: ind, jnd, knd
        1543             :
        1544             :       !-----------------------------------------------------------------------
        1545             :       !  Calculate constituent dependent specific heat, gas constant and cappa
        1546             :       !-----------------------------------------------------------------------
        1547           0 :       call get_R_dry(tracer, thermodynamic_active_species_idx, rgas_var)
        1548           0 :       call get_cp_dry(tracer, thermodynamic_active_species_idx, cp_var)
        1549             :       !$omp parallel do private(ind,jnd,knd)
        1550           0 :       do knd = 1, SIZE(tracer, 3)
        1551           0 :          do jnd = 1, SIZE(tracer, 2)
        1552           0 :             do ind = 1, SIZE(tracer, 1)
        1553           0 :                kappav(ind,jnd,knd) = rgas_var(ind,jnd,knd) / cp_var(ind,jnd,knd)
        1554             :             end do
        1555             :          end do
        1556             :       end do
        1557             :
        1558           0 :       if (present(cpv)) then
        1559           0 :          cpv(:,:,:) = cp_var(:,:,:)
        1560             :       end if
        1561             :
        1562           0 :    end subroutine cam_thermo_calc_kappav_2hd
        1563             :
        1564             :    !===========================================================================
        1565             :    !
        1566             :    !***************************************************************************
        1567             :    !
        1568             :    ! compute column integrated total energy consistent with vertical
        1569             :    !    coordinate as well as vertical integrals of water mass (H2O,wv,liq,ice)
        1570             :    !
        1571             :    ! if subroutine is asked to compute "te" then the latent heat terms are
        1572             :    !    added to the kinetic (ke), internal + geopotential (se)  energy terms
        1573             :    !
        1574             :    ! subroutine assumes that enthalpy term (rho*cp*T) uses dry air heat capacity
        1575             :    !
        1576             :    !***************************************************************************
        1577             :    !
        1578    35814528 :    subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, &
        1579    35814528 :         cp_or_cv, U, V, T, vcoord, ptop, phis, z_mid, dycore_idx, qidx,       &
        1580    35814528 :         te, se, po, ke, wv, H2O, liq, ice)
        1581             :
        1582             :       use cam_logfile,     only: iulog
        1583             :       use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure
        1584             :       use air_composition, only: wv_idx
        1585             :       use physconst,       only: rga, latvap, latice
        1586             :
        1587             :       ! Dummy arguments
        1588             :       ! tracer: tracer mixing ratio
        1589             :       !
        1590             :       ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry
        1591             :       real(r8), intent(in)            :: tracer(:,:,:)
        1592             :       logical, intent(in)             :: moist_mixing_ratio
        1593             :       ! pdel: pressure level thickness
        1594             :       real(r8), intent(in)            :: pdel_in(:,:)
        1595             :       ! cp_or_cv: dry air heat capacity under constant pressure or
        1596             :       !           constant volume (depends on vcoord)
        1597             :       real(r8), intent(in)            :: cp_or_cv(:,:)
        1598             :       real(r8), intent(in)            :: U(:,:)
        1599             :       real(r8), intent(in)            :: V(:,:)
        1600             :       real(r8), intent(in)            :: T(:,:)
        1601             :       integer,  intent(in)            :: vcoord ! vertical coordinate
        1602             :       real(r8), intent(in),  optional :: ptop(:)
        1603             :       real(r8), intent(in),  optional :: phis(:)
        1604             :       real(r8), intent(in),  optional :: z_mid(:,:)
        1605             :       ! dycore_idx: use dycore index for thermodynamic active species
        1606             :       logical,  intent(in),  optional :: dycore_idx
        1607             :       ! qidx: Index of water vapor
        1608             :       integer,  intent(in),  optional :: qidx
        1609             :       ! H2O: vertically integrated total water
        1610             :       real(r8), intent(out), optional :: H2O(:)
        1611             :       ! TE: vertically integrated total energy
        1612             :       real(r8), intent(out), optional :: te (:)
        1613             :       ! KE: vertically integrated kinetic energy
        1614             :       real(r8), intent(out), optional :: ke (:)
        1615             :       ! SE: vertically integrated enthalpy (pressure coordinate)
        1616             :       !     or internal energy (z coordinate)
        1617             :       real(r8), intent(out), optional :: se (:)
        1618             :       ! PO: vertically integrated PHIS term (pressure coordinate)
        1619             :       !     or potential energy (z coordinate)
        1620             :       real(r8), intent(out), optional :: po (:)
        1621             :       ! WV: vertically integrated water vapor
        1622             :       real(r8), intent(out), optional :: wv (:)
        1623             :       ! liq: vertically integrated liquid
        1624             :       real(r8), intent(out), optional :: liq(:)
        1625             :       ! ice: vertically integrated ice
        1626             :       real(r8), intent(out), optional :: ice(:)
        1627             :
        1628             :       ! Local variables
        1629    71629056 :       real(r8) :: ke_vint(SIZE(tracer, 1))  ! Vertical integral of KE
        1630    71629056 :       real(r8) :: se_vint(SIZE(tracer, 1))  ! Vertical integral of enthalpy or internal energy
        1631    71629056 :       real(r8) :: po_vint(SIZE(tracer, 1))  ! Vertical integral of PHIS or potential energy
        1632    71629056 :       real(r8) :: wv_vint(SIZE(tracer, 1))  ! Vertical integral of wv
        1633    71629056 :       real(r8) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq
        1634    71629056 :       real(r8) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice
        1635    71629056 :       real(r8) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) !moist pressure level thickness
        1636             :       real(r8)                      :: latsub ! latent heat of sublimation
        1637             :
        1638             :       integer                       :: ierr
        1639             :       integer                       :: kdx, idx ! coord indices
        1640             :       integer                       :: qdx      ! tracer index
        1641             :       integer                       :: wvidx    ! water vapor index
        1642    35814528 :       integer,          allocatable :: species_idx(:)
        1643    35814528 :       integer,          allocatable :: species_liq_idx(:)
        1644    35814528 :       integer,          allocatable :: species_ice_idx(:)
        1645             :       character(len=*), parameter   :: subname = 'get_hydrostatic_energy'
        1646             :
        1647   107443584 :       allocate(species_idx(thermodynamic_active_species_num), stat=ierr)
        1648    35814528 :       if ( ierr /= 0 ) then
        1649           0 :          call endrun(subname//': allocation error for species_idx array')
        1650             :       end if
        1651   107443584 :       allocate(species_liq_idx(thermodynamic_active_species_liq_num), stat=ierr)
        1652    35814528 :       if ( ierr /= 0 ) then
        1653           0 :          call endrun(subname//': allocation error for species_liq_idx array')
        1654             :       end if
        1655   107443584 :       allocate(species_ice_idx(thermodynamic_active_species_ice_num), stat=ierr)
        1656    35814528 :       if ( ierr /= 0 ) then
        1657           0 :          call endrun(subname//': allocation error for species_ice_idx array')
        1658             :       end if
        1659             :
        1660    35814528 :       if (present(dycore_idx))then
        1661           0 :          if (dycore_idx) then
        1662           0 :             species_idx(:) = thermodynamic_active_species_idx_dycore(:)
        1663           0 :             species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:)
        1664           0 :             species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:)
        1665             :          else
        1666           0 :             species_idx(:) = thermodynamic_active_species_idx(:)
        1667           0 :             species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
        1668           0 :             species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
        1669             :          end if
        1670             :       else
        1671   143258112 :          species_idx(:) = thermodynamic_active_species_idx(:)
        1672    71629056 :          species_liq_idx(:) = thermodynamic_active_species_liq_idx(:)
        1673    71629056 :          species_ice_idx(:) = thermodynamic_active_species_ice_idx(:)
        1674             :       end if
        1675             :
        1676    35814528 :       if (present(qidx)) then
        1677           0 :          wvidx = qidx
        1678             :       else
        1679    35814528 :          wvidx = wv_idx
        1680             :       end if
        1681             :
        1682    35814528 :       if (moist_mixing_ratio) then
        1683 15584317056 :         pdel     = pdel_in
        1684             :       else
        1685           0 :         pdel     = pdel_in
        1686           0 :         do qdx = dry_air_species_num+1, thermodynamic_active_species_num
        1687           0 :           pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx))
        1688             :         end do
        1689             :       end if
        1690             :
        1691   598019328 :       ke_vint = 0._r8
        1692   598019328 :       se_vint = 0._r8
        1693    71629056 :       select case (vcoord)
        1694             :       case(vc_moist_pressure, vc_dry_pressure)
        1695    35814528 :          if (.not. present(ptop).or. (.not. present(phis))) then
        1696           0 :             write(iulog, *) subname, ' ptop and phis must be present for ',     &
        1697           0 :                  'moist/dry pressure vertical coordinate'
        1698             :             call endrun(subname//':  ptop and phis must be present for '//      &
        1699           0 :                  'moist/dry pressure vertical coordinate')
        1700             :          end if
        1701   598019328 :          po_vint = ptop
        1702   966992256 :          do kdx = 1, SIZE(tracer, 2)
        1703 15584317056 :             do idx = 1, SIZE(tracer, 1)
        1704 29234649600 :                ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) *                &
        1705 43851974400 :                     0.5_r8 * (U(idx, kdx)**2 + V(idx, kdx)**2)) * rga
        1706 14617324800 :                se_vint(idx) = se_vint(idx) + (T(idx, kdx) *                   &
        1707 14617324800 :                     cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
        1708 15548502528 :                po_vint(idx) =  po_vint(idx)+pdel(idx, kdx)
        1709             :
        1710             :             end do
        1711             :          end do
        1712   633833856 :          do idx = 1, SIZE(tracer, 1)
        1713   598019328 :             po_vint(idx) =  (phis(idx) * po_vint(idx) * rga)
        1714             :          end do
        1715             :       case(vc_height)
        1716           0 :          if (.not. present(phis)) then
        1717           0 :             write(iulog, *) subname, ' phis must be present for ',     &
        1718           0 :                  'heigt-based vertical coordinate'
        1719             :             call endrun(subname//':  phis must be present for '//      &
        1720           0 :                  'height-based vertical coordinate')
        1721             :          end if
        1722           0 :          po_vint = 0._r8
        1723           0 :          do kdx = 1, SIZE(tracer, 2)
        1724           0 :             do idx = 1, SIZE(tracer, 1)
        1725           0 :                ke_vint(idx) = ke_vint(idx) + (pdel(idx, kdx) *                &
        1726           0 :                     0.5_r8 * (U(idx, kdx)**2 + V(idx, kdx)**2) * rga)
        1727           0 :                se_vint(idx) = se_vint(idx) + (T(idx, kdx) *                   &
        1728           0 :                     cp_or_cv(idx, kdx) * pdel(idx, kdx) * rga)
        1729             :                ! z_mid is height above ground
        1730           0 :                po_vint(idx) = po_vint(idx) + (z_mid(idx, kdx) +               &
        1731           0 :                     phis(idx) * rga) * pdel(idx, kdx)
        1732             :             end do
        1733             :          end do
        1734             :       case default
        1735           0 :          write(iulog, *) subname, ' vertical coordinate not supported: ', vcoord
        1736    35814528 :          call endrun(subname//': vertical coordinate not supported')
        1737             :       end select
        1738    35814528 :       if (present(te)) then
        1739   598019328 :          te  = se_vint + po_vint+ ke_vint
        1740             :       end if
        1741    35814528 :       if (present(se)) then
        1742   274040496 :          se = se_vint
        1743             :       end if
        1744    35814528 :       if (present(po)) then
        1745   274040496 :          po = po_vint
        1746             :       end if
        1747    35814528 :       if (present(ke)) then
        1748   274040496 :          ke = ke_vint
        1749             :       end if
        1750             :       !
        1751             :       ! vertical integral of total liquid water
        1752             :       !
        1753    35814528 :       if (.not.moist_mixing_ratio) then
        1754           0 :         pdel = pdel_in! set pseudo density to dry
        1755             :       end if
        1756             :
        1757   598019328 :       wv_vint = 0._r8
        1758   966992256 :       do kdx = 1, SIZE(tracer, 2)
        1759 15584317056 :         do idx = 1, SIZE(tracer, 1)
        1760 43851974400 :           wv_vint(idx) = wv_vint(idx) + (tracer(idx, kdx, wvidx) *       &
        1761 59400476928 :                pdel(idx, kdx) * rga)
        1762             :         end do
        1763             :       end do
        1764   293443128 :       if (present(wv)) wv = wv_vint
        1765             :
        1766   598019328 :       liq_vint = 0._r8
        1767    71629056 :       do qdx = 1, thermodynamic_active_species_liq_num
        1768  1002806784 :          do kdx = 1, SIZE(tracer, 2)
        1769 15584317056 :             do idx = 1, SIZE(tracer, 1)
        1770 29234649600 :                liq_vint(idx) = liq_vint(idx) + (pdel(idx, kdx) *         &
        1771 44783152128 :                     tracer(idx, kdx, species_liq_idx(qdx)) * rga)
        1772             :             end do
        1773             :          end do
        1774             :       end do
        1775   293443128 :       if (present(liq)) liq = liq_vint
        1776             :
        1777             :       !
        1778             :       ! vertical integral of total frozen (ice) water
        1779             :       !
        1780   598019328 :       ice_vint = 0._r8
        1781    71629056 :       do qdx = 1, thermodynamic_active_species_ice_num
        1782  1002806784 :          do kdx = 1, SIZE(tracer, 2)
        1783 15584317056 :             do idx = 1, SIZE(tracer, 1)
        1784 29234649600 :                ice_vint(idx) = ice_vint(idx) + (pdel(idx, kdx) *              &
        1785 44783152128 :                     tracer(idx, kdx, species_ice_idx(qdx))  * rga)
        1786             :             end do
        1787             :          end do
        1788             :       end do
        1789   293443128 :       if (present(ice)) ice = ice_vint
        1790             :       ! Compute vertical integrals of total water.
        1791    35814528 :       if (present(H2O)) then
        1792   598019328 :          H2O = wv_vint + liq_vint + ice_vint
        1793             :       end if
        1794             :       !
        1795             :       ! latent heat terms depend on enthalpy reference state
        1796             :       !
        1797    35814528 :       latsub = latvap + latice
        1798    35814528 :       if (present(te)) then
        1799    71629056 :          select case (TRIM(enthalpy_reference_state))
        1800             :          case('ice')
        1801   598019328 :             te = te + (latsub * wv_vint) + (latice * liq_vint)
        1802             :          case('liq')
        1803           0 :             te = te + (latvap * wv_vint) - (latice * ice_vint)
        1804             :          case('wv')
        1805           0 :             te = te - (latvap * liq_vint) - (latsub * ice_vint)
        1806             :          case default
        1807           0 :             write(iulog, *) subname, ' enthalpy reference state not ',        &
        1808           0 :                  'supported: ', TRIM(enthalpy_reference_state)
        1809    71629056 :             call endrun(subname//': enthalpy reference state not supported')
        1810             :          end select
        1811             :       end if
        1812    35814528 :       deallocate(species_idx, species_liq_idx, species_ice_idx)
        1813    35814528 :     end subroutine get_hydrostatic_energy_1hd
        1814             :
        1815             : end module cam_thermo

Generated by: LCOV version 1.14