LCOV - code coverage report
Current view: top level - utils - air_composition.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 259 467 55.5 %
Date: 2024-12-17 17:57:11 Functions: 10 14 71.4 %

          Line data    Source code
       1             : ! air_composition module defines major species of the atmosphere and manages
       2             : ! the physical properties that are dependent on the composition of air
       3             : module air_composition
       4             : 
       5             :    use shr_kind_mod,   only: r8 => shr_kind_r8
       6             :    use cam_abortutils, only: endrun
       7             : 
       8             :    implicit none
       9             :    private
      10             :    save
      11             : 
      12             :    public  :: air_composition_readnl
      13             :    public  :: air_composition_init
      14             :    public  :: dry_air_composition_update
      15             :    public  :: water_composition_update
      16             : 
      17             :    ! get_cp_dry: (generalized) heat capacity for dry air
      18             :    public :: get_cp_dry
      19             :    ! get_cp: (generalized) heat capacity
      20             :    public :: get_cp
      21             :    ! get_R_dry: (generalized) dry air gas constant
      22             :    public :: get_R_dry
      23             :    ! get_R: Compute generalized R
      24             :    public :: get_R
      25             :    ! get_mbarv: molecular weight of dry air
      26             :    public :: get_mbarv
      27             : 
      28             :    private :: air_species_info
      29             : 
      30             :    integer,  parameter :: unseti = -HUGE(1)
      31             :    real(r8), parameter :: unsetr = HUGE(1.0_r8)
      32             : 
      33             :    ! composition of air
      34             :    !
      35             :    integer, parameter :: num_names_max = 20 ! Should match namelist definition
      36             :    character(len=6)   :: dry_air_species(num_names_max)
      37             :    character(len=6)   :: water_species_in_air(num_names_max)
      38             : 
      39             :    integer, protected, public :: dry_air_species_num
      40             :    integer, protected, public :: water_species_in_air_num
      41             : 
      42             :    ! Thermodynamic variables
      43             :    integer,               protected, public :: thermodynamic_active_species_num = unseti
      44             :    integer,  allocatable, protected, public :: thermodynamic_active_species_idx(:)
      45             :    integer,  allocatable,            public :: thermodynamic_active_species_idx_dycore(:)
      46             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_cp(:)
      47             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_cv(:)
      48             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_R(:)
      49             :    ! thermodynamic_active_species_mwi: inverse molecular weights dry air
      50             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_mwi(:)
      51             :    ! thermodynamic_active_species_kv: molecular diffusion
      52             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_kv(:)
      53             :    ! thermodynamic_active_species_kc: thermal conductivity
      54             :    real(r8), allocatable, protected, public :: thermodynamic_active_species_kc(:)
      55             :    !
      56             :    ! for energy computations liquid and ice species need to be identified
      57             :    !
      58             :    ! thermodynamic_active_species_liq_num: number of liquid water species
      59             :    integer,               protected, public :: thermodynamic_active_species_liq_num = unseti
      60             :    ! thermodynamic_active_species_ice_num: number of frozen water species
      61             :    integer,               protected, public :: thermodynamic_active_species_ice_num = unseti
      62             :    ! thermodynamic_active_species_liq_idx: index of liquid water species
      63             :    integer,  allocatable, protected, public :: thermodynamic_active_species_liq_idx(:)
      64             :    ! thermodynamic_active_species_liq_idx_dycore: index of liquid water species
      65             :    integer,  allocatable,            public :: thermodynamic_active_species_liq_idx_dycore(:)
      66             :    ! thermodynamic_active_species_ice_idx: index of ice water species
      67             :    integer,  allocatable, protected, public :: thermodynamic_active_species_ice_idx(:)
      68             :    ! thermodynamic_active_species_ice_idx_dycore: index of ice water species
      69             :    integer,  allocatable,            public :: thermodynamic_active_species_ice_idx_dycore(:)
      70             :    ! enthalpy_reference_state: choices: 'ice', 'liq', 'wv'
      71             :    character(len=3), public, protected :: enthalpy_reference_state = 'xxx'
      72             : 
      73             :    integer, protected, public :: wv_idx = -1 ! Water vapor index
      74             : 
      75             :    !------------- Variables for consistent themodynamics --------------------
      76             :    !
      77             : 
      78             :    ! standard dry air (constant composition)
      79             :    real(r8), public, protected :: mmro2 = unsetr  ! Mass mixing ratio of O2
      80             :    real(r8), public, protected :: mmrn2 = unsetr  ! Mass mixing ratio of N2
      81             :    real(r8), public, protected :: o2_mwi = unsetr ! Inverse mol. weight of O2
      82             :    real(r8), public, protected :: n2_mwi = unsetr ! Inverse mol. weight of N2
      83             :    real(r8), public, protected :: mbar = unsetr   ! Mean mass at mid level
      84             : 
      85             :    ! coefficients in expressions for molecular diffusion coefficients
      86             :    ! kv1,..,kv3 are coefficients for kmvis calculation
      87             :    ! kc1,..,kc3 are coefficients for kmcnd calculation
      88             :    ! Liu, H.-L., et al. (2010), Thermosphere extension of the Whole Atmosphere Community Climate Model,
      89             :    !    J. Geophys. Res., 115, A12302, doi:10.1029/2010JA015586.
      90             :    real(r8), public, parameter :: kv1 = 4.03_r8 * 1.e-7_r8
      91             :    real(r8), public, parameter :: kv2 = 3.42_r8 * 1.e-7_r8
      92             :    real(r8), public, parameter :: kv3 = 3.9_r8 * 1.e-7_r8
      93             :    real(r8), public, parameter :: kc1 = 56._r8 * 1.e-5_r8
      94             :    real(r8), public, parameter :: kc2 = 56._r8 * 1.e-5_r8
      95             :    real(r8), public, parameter :: kc3 = 75.9_r8 * 1.e-5_r8
      96             : 
      97             :    real(r8), public, parameter :: kv_temp_exp = 0.69_r8
      98             :    real(r8), public, parameter :: kc_temp_exp = 0.69_r8
      99             : 
     100             :    ! cpairv:  composition dependent specific heat at constant pressure
     101             :    real(r8), public, protected, allocatable :: cpairv(:,:,:)
     102             :    ! rairv: composition dependent gas "constant"
     103             :    real(r8), public, protected, allocatable :: rairv(:,:,:)
     104             :    ! cappav: rairv / cpairv
     105             :    real(r8), public, protected, allocatable :: cappav(:,:,:)
     106             :    ! mbarv: composition dependent atmosphere mean mass
     107             :    real(r8), public, protected, allocatable :: mbarv(:,:,:)
     108             :    ! cp_or_cv_dycore:  enthalpy or internal energy scaling factor for 
     109             :    !                   energy consistency
     110             :    real(r8), public, protected, allocatable :: cp_or_cv_dycore(:,:,:)
     111             :    !
     112             :    ! Interfaces for public routines
     113             :    interface get_cp_dry
     114             :       module procedure get_cp_dry_1hd
     115             :       module procedure get_cp_dry_2hd
     116             :    end interface get_cp_dry
     117             : 
     118             :    interface get_cp
     119             :       module procedure get_cp_1hd
     120             :       module procedure get_cp_2hd
     121             :    end interface get_cp
     122             : 
     123             :    interface get_R_dry
     124             :       module procedure get_R_dry_1hd
     125             :       module procedure get_R_dry_2hd
     126             :    end interface get_R_dry
     127             : 
     128             :    interface get_R
     129             :       module procedure get_R_1hd
     130             :       module procedure get_R_2hd
     131             :    end interface get_R
     132             : 
     133             :    interface get_mbarv
     134             :       module procedure get_mbarv_1hd
     135             :    end interface get_mbarv
     136             : 
     137             : CONTAINS
     138             : 
     139             :    ! Read namelist variables.
     140        1536 :    subroutine air_composition_readnl(nlfile)
     141             :       use namelist_utils, only: find_group_name
     142             :       use spmd_utils,     only: masterproc, mpicom, masterprocid
     143             :       use spmd_utils,     only: mpi_character
     144             :       use cam_logfile,    only: iulog
     145             : 
     146             :       ! Dummy argument: filepath for file containing namelist input
     147             :       character(len=*), intent(in) :: nlfile
     148             : 
     149             :       ! Local variables
     150             :       integer                     :: unitn, ierr, indx
     151             :       integer,          parameter :: lsize = 76
     152             :       character(len=*), parameter :: subname = 'air_composition_readnl :: '
     153             :       character(len=lsize)        :: banner
     154             :       character(len=lsize)        :: bline
     155             : 
     156             :       ! Variable components of dry air and water species in air
     157             :       namelist /air_composition_nl/ dry_air_species, water_species_in_air
     158             :       !-----------------------------------------------------------------------
     159             : 
     160        1536 :       banner = repeat('*', lsize)
     161        1536 :       bline = "***"//repeat(' ', lsize - 6)//"***"
     162             : 
     163             :       ! Read variable components of dry air and water species in air
     164       32256 :       dry_air_species = (/ (' ', indx = 1, num_names_max) /)
     165       32256 :       water_species_in_air = (/ (' ', indx = 1, num_names_max) /)
     166             : 
     167        1536 :       if (masterproc) then
     168           2 :          open(newunit=unitn, file=trim(nlfile), status='old')
     169           2 :          call find_group_name(unitn, 'air_composition_nl', status=ierr)
     170           2 :          if (ierr == 0) then
     171           2 :             read(unitn, air_composition_nl, iostat=ierr)
     172           2 :             if (ierr /= 0) then
     173           0 :                call endrun(subname//'ERROR reading namelist, air_composition_nl')
     174             :             end if
     175             :          end if
     176           2 :          close(unitn)
     177             :       end if
     178             : 
     179             :       call mpi_bcast(dry_air_species, len(dry_air_species)*num_names_max,     &
     180        1536 :            mpi_character, masterprocid, mpicom, ierr)
     181        1536 :       if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: dry_air_species")
     182             :       call mpi_bcast(water_species_in_air,                                    &
     183             :            len(water_species_in_air)*num_names_max, mpi_character,            &
     184        1536 :            masterprocid, mpicom, ierr)
     185        1536 :       if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: water_species_in_air")
     186             : 
     187        1536 :       dry_air_species_num = 0
     188        1536 :       water_species_in_air_num = 0
     189       32256 :       do indx = 1, num_names_max
     190       30720 :          if ( (LEN_TRIM(dry_air_species(indx)) > 0) .and.                     &
     191             :               (TRIM(dry_air_species(indx)) /= 'N2')) then
     192           0 :             dry_air_species_num = dry_air_species_num + 1
     193             :          end if
     194       32256 :          if (LEN_TRIM(water_species_in_air(indx)) > 0) then
     195        9216 :             water_species_in_air_num = water_species_in_air_num + 1
     196             :          end if
     197             :       end do
     198             : 
     199             :       ! Initialize number of thermodynamically active species
     200             :       thermodynamic_active_species_num =                                      &
     201        1536 :            dry_air_species_num + water_species_in_air_num
     202             : 
     203        1536 :       if (masterproc) then
     204           2 :          write(iulog, *) banner
     205           2 :          write(iulog, *) bline
     206             : 
     207           2 :          if (dry_air_species_num == 0) then
     208           2 :             write(iulog, *) " Thermodynamic properties of dry air are ",      &
     209           4 :                  "fixed at troposphere values"
     210             :          else
     211           0 :             write(iulog, *) " Thermodynamic properties of dry air are ",      &
     212           0 :                  "based on variable composition of the following species:"
     213           0 :             do indx = 1, dry_air_species_num
     214           0 :                write(iulog, *) '   ',  trim(dry_air_species(indx))
     215             :             end do
     216           0 :             write(iulog,*) ' '
     217             :          end if
     218           2 :          write(iulog,*) " Thermodynamic properties of moist air are ",        &
     219           4 :               "based on variable composition of the following water species:"
     220          14 :          do indx = 1, water_species_in_air_num
     221          14 :             write(iulog, *) '   ', trim(water_species_in_air(indx))
     222             :          end do
     223           2 :          write(iulog, *) bline
     224           2 :          write(iulog, *) banner
     225             :       end if
     226             : 
     227        1536 :    end subroutine air_composition_readnl
     228             : 
     229             :    !===========================================================================
     230             : 
     231        1536 :    subroutine air_composition_init()
     232             :       use string_utils, only: int2str
     233             :       use spmd_utils,   only: masterproc
     234             :       use cam_logfile,  only: iulog
     235             :       use physconst,    only: r_universal, cpair, rair, cpwv, rh2o, cpliq, cpice, mwdry
     236             :       use constituents, only: cnst_get_ind, cnst_mw
     237             :       use ppgrid,       only: pcols, pver, begchunk, endchunk
     238             :       integer  :: icnst, ix, isize, ierr, idx
     239             :       integer  :: liq_num, ice_num
     240        3072 :       integer  :: liq_idx(water_species_in_air_num)
     241        3072 :       integer  :: ice_idx(water_species_in_air_num)
     242             :       logical  :: has_liq, has_ice
     243             :       real(r8) :: mw
     244             : 
     245             :       character(len=*), parameter :: subname = 'composition_init'
     246             :       character(len=*), parameter :: errstr = subname//": failed to allocate "
     247             : 
     248             :       !
     249             :       ! define cp and R for species in species_name
     250             :       !
     251             :       ! Last major species in namelist dry_air_species is derived from the
     252             :       !    other major species (since the sum of dry mixing ratios for
     253             :       !    major species of dry air add must add to one)
     254             :       !
     255             :       ! cv = R * dofx / 2;   cp = R * (1 + (dofx / 2))
     256             :       ! DOF == Degrees of Freedom
     257             :       ! dof1 = monatomic ideal gas, 3 translational DOF
     258             :       real(r8), parameter :: dof1 = 3._r8
     259             :       real(r8), parameter :: cv1 = 0.5_r8 * r_universal * dof1
     260             :       real(r8), parameter :: cp1 = 0.5_r8 * r_universal * (2._r8 + dof1)
     261             :       ! dof2 = diatomic ideal gas, 3 translational + 2 rotational = 5 DOF
     262             :       real(r8), parameter :: dof2 = 5._r8
     263             :       real(r8), parameter :: cv2 = 0.5_r8 * r_universal * dof2
     264             :       real(r8), parameter :: cp2 = 0.5_r8 * r_universal * (2._r8 + dof2)
     265             :       ! dof3 = polyatomic ideal gas, 3 translational + 3 rotational = 6 DOF
     266             :       real(r8), parameter :: dof3 = 6._r8
     267             :       real(r8), parameter :: cv3 = 0.5_r8 * r_universal * dof3
     268             :       real(r8), parameter :: cp3 = 0.5_r8 * r_universal * (2._r8 + dof3)
     269             : 
     270        1536 :       liq_num = 0
     271        1536 :       ice_num = 0
     272        1536 :       has_liq = .false.
     273        1536 :       has_ice = .false.
     274             :       ! standard dry air (constant composition)
     275        1536 :       o2_mwi = 1._r8 / 32._r8
     276        1536 :       n2_mwi = 1._r8 / 28._r8
     277        1536 :       mmro2 = 0.235_r8
     278        1536 :       mmrn2 = 0.765_r8
     279        1536 :       mbar = 1._r8 / ((mmro2 * o2_mwi) + (mmrn2 * n2_mwi))
     280             : 
     281             :       ! init for variable composition dry air
     282             : 
     283        1536 :       isize = dry_air_species_num + water_species_in_air_num
     284        4608 :       allocate(thermodynamic_active_species_idx(isize), stat=ierr)
     285        1536 :       if (ierr /= 0) then
     286           0 :          call endrun(errstr//"thermodynamic_active_species_idx")
     287             :       end if
     288        3072 :       allocate(thermodynamic_active_species_idx_dycore(isize), stat=ierr)
     289        1536 :       if (ierr /= 0) then
     290           0 :          call endrun(errstr//"thermodynamic_active_species_idx_dycore")
     291             :       end if
     292        4608 :       allocate(thermodynamic_active_species_cp(0:isize), stat=ierr)
     293        1536 :       if (ierr /= 0) then
     294           0 :          call endrun(errstr//"thermodynamic_active_species_cp")
     295             :       end if
     296        3072 :       allocate(thermodynamic_active_species_cv(0:isize), stat=ierr)
     297        1536 :       if (ierr /= 0) then
     298           0 :          call endrun(errstr//"thermodynamic_active_species_cv")
     299             :       end if
     300        3072 :       allocate(thermodynamic_active_species_R(0:isize), stat=ierr)
     301        1536 :       if (ierr /= 0) then
     302           0 :          call endrun(errstr//"thermodynamic_active_species_R")
     303             :       end if
     304             : 
     305        1536 :       isize = dry_air_species_num
     306        4608 :       allocate(thermodynamic_active_species_mwi(0:isize), stat=ierr)
     307        1536 :       if (ierr /= 0) then
     308           0 :          call endrun(errstr//"thermodynamic_active_species_mwi")
     309             :       end if
     310        3072 :       allocate(thermodynamic_active_species_kv(0:isize), stat=ierr)
     311        1536 :       if (ierr /= 0) then
     312           0 :          call endrun(errstr//"thermodynamic_active_species_kv")
     313             :       end if
     314        3072 :       allocate(thermodynamic_active_species_kc(0:isize), stat=ierr)
     315        1536 :       if (ierr /= 0) then
     316           0 :          call endrun(errstr//"thermodynamic_active_species_kc")
     317             :       end if
     318             :       !------------------------------------------------------------------------
     319             :       !  Allocate constituent dependent properties
     320             :       !------------------------------------------------------------------------
     321        4608 :       allocate(cpairv(pcols,pver,begchunk:endchunk), stat=ierr)
     322        1536 :       if (ierr /= 0) then
     323           0 :          call endrun(errstr//"cpairv")
     324             :       end if
     325        4608 :       allocate(rairv(pcols,pver,begchunk:endchunk),  stat=ierr)
     326        1536 :       if (ierr /= 0) then
     327           0 :          call endrun(errstr//"rairv")
     328             :       end if
     329        4608 :       allocate(cappav(pcols,pver,begchunk:endchunk), stat=ierr)
     330        1536 :       if (ierr /= 0) then
     331           0 :          call endrun(errstr//"cappav")
     332             :       end if
     333        4608 :       allocate(mbarv(pcols,pver,begchunk:endchunk),  stat=ierr)
     334        1536 :       if (ierr /= 0) then
     335           0 :          call endrun(errstr//"mbarv")
     336             :       end if
     337        4608 :       allocate(cp_or_cv_dycore(pcols,pver,begchunk:endchunk), stat=ierr)
     338        1536 :       if (ierr /= 0) then
     339           0 :          call endrun(errstr//"cp_or_cv_dycore")
     340             :       end if
     341             : 
     342       10752 :       thermodynamic_active_species_idx        = -HUGE(1)
     343       10752 :       thermodynamic_active_species_idx_dycore = -HUGE(1)
     344       12288 :       thermodynamic_active_species_cp         = 0.0_r8
     345       12288 :       thermodynamic_active_species_cv         = 0.0_r8
     346       12288 :       thermodynamic_active_species_R          = 0.0_r8
     347        3072 :       thermodynamic_active_species_mwi        = 0.0_r8
     348        3072 :       thermodynamic_active_species_kv         = 0.0_r8
     349        3072 :       thermodynamic_active_species_kc         = 0.0_r8
     350             :       !------------------------------------------------------------------------
     351             :       !  Initialize constituent dependent properties
     352             :       !------------------------------------------------------------------------
     353     9797280 :       cpairv(:pcols, :pver, begchunk:endchunk) = cpair
     354     9797280 :       rairv(:pcols,  :pver, begchunk:endchunk) = rair
     355     9797280 :       cappav(:pcols, :pver, begchunk:endchunk) = rair / cpair
     356     9797280 :       mbarv(:pcols,  :pver, begchunk:endchunk) = mwdry
     357             :       !
     358        1536 :       if (dry_air_species_num > 0) then
     359             :          !
     360             :          ! The last major species in dry_air_species is derived from the
     361             :          !    others and constants associated with it are initialized here
     362             :          !
     363           0 :          if (TRIM(dry_air_species(dry_air_species_num + 1)) == 'N2') then
     364           0 :             call air_species_info('N', ix, mw)
     365           0 :             mw = 2.0_r8 * mw
     366           0 :             icnst = 0 ! index for the derived tracer N2
     367           0 :             thermodynamic_active_species_cp(icnst) = cp2 / mw
     368           0 :             thermodynamic_active_species_cv(icnst) = cv2 / mw !N2
     369           0 :             thermodynamic_active_species_R  (icnst) = r_universal / mw
     370           0 :             thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
     371           0 :             thermodynamic_active_species_kv(icnst)  = kv2
     372           0 :             thermodynamic_active_species_kc(icnst)  = kc2
     373             :             !
     374             :             ! if last major species is not N2 then add code here
     375             :             !
     376             :          else
     377           0 :             write(iulog, *) subname, ' derived major species not found: ',    &
     378           0 :                  dry_air_species(dry_air_species_num)
     379           0 :             call endrun(subname//': derived major species not found')
     380             :          end if
     381             :       else
     382             :          !
     383             :          ! dry air is not species dependent
     384             :          !
     385        1536 :          icnst = 0
     386        1536 :          thermodynamic_active_species_cp (icnst) = cpair
     387        1536 :          thermodynamic_active_species_cv (icnst) = cpair - rair
     388        1536 :          thermodynamic_active_species_R  (icnst) = rair
     389             :       end if
     390             :       !
     391             :       !************************************************************************
     392             :       !
     393             :       ! add prognostic components of dry air
     394             :       !
     395             :       !************************************************************************
     396             :       !
     397        1536 :       icnst = 1
     398        1536 :       do idx = 1, dry_air_species_num
     399           0 :          select case (TRIM(dry_air_species(idx)))
     400             :             !
     401             :             ! O
     402             :             !
     403             :          case('O')
     404           0 :             call air_species_info('O', ix, mw)
     405           0 :             thermodynamic_active_species_idx(icnst) = ix
     406           0 :             thermodynamic_active_species_cp (icnst) = cp1 / mw
     407           0 :             thermodynamic_active_species_cv (icnst) = cv1 / mw
     408           0 :             thermodynamic_active_species_R  (icnst) = r_universal / mw
     409           0 :             thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
     410           0 :             thermodynamic_active_species_kv(icnst)  = kv3
     411           0 :             thermodynamic_active_species_kc(icnst)  = kc3
     412           0 :             icnst = icnst + 1
     413             :             !
     414             :             ! O2
     415             :             !
     416             :          case('O2')
     417           0 :             call air_species_info('O2', ix, mw)
     418           0 :             thermodynamic_active_species_idx(icnst) = ix
     419           0 :             thermodynamic_active_species_cp (icnst) = cp2 / mw
     420           0 :             thermodynamic_active_species_cv (icnst) = cv2 / mw
     421           0 :             thermodynamic_active_species_R  (icnst) = r_universal / mw
     422           0 :             thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
     423           0 :             thermodynamic_active_species_kv(icnst)  = kv1
     424           0 :             thermodynamic_active_species_kc(icnst)  = kc1
     425           0 :             icnst = icnst + 1
     426             :             !
     427             :             ! H
     428             :             !
     429             :          case('H')
     430           0 :             call air_species_info('H', ix, mw)
     431           0 :             thermodynamic_active_species_idx(icnst) = ix
     432           0 :             thermodynamic_active_species_cp (icnst) = cp1 / mw
     433           0 :             thermodynamic_active_species_cv (icnst) = cv1 / mw
     434           0 :             thermodynamic_active_species_R  (icnst) = r_universal / mw
     435           0 :             thermodynamic_active_species_mwi(icnst) = 1.0_r8 / mw
     436             :             ! Hydrogen not included in calculation of diffusivity and conductivity
     437           0 :             thermodynamic_active_species_kv(icnst)  = 0.0_r8
     438           0 :             thermodynamic_active_species_kc(icnst)  = 0.0_r8
     439           0 :             icnst = icnst + 1
     440             :             !
     441             :             ! If support for more major species is to be included add code here
     442             :             !
     443             :          case default
     444           0 :             write(iulog, *) subname, ' dry air component not found: ',        &
     445           0 :                  dry_air_species(idx)
     446           0 :             call endrun(subname//': dry air component not found')
     447             :          end select
     448             : 
     449        1536 :          if (masterproc) then
     450           0 :             write(iulog, *) "Dry air composition ",                           &
     451           0 :                  TRIM(dry_air_species(idx)),                                  &
     452           0 :                  icnst-1,thermodynamic_active_species_idx(icnst-1),           &
     453           0 :                  thermodynamic_active_species_mwi(icnst-1),                   &
     454           0 :                  thermodynamic_active_species_cp(icnst-1),                    &
     455           0 :                  thermodynamic_active_species_cv(icnst-1)
     456             :          end if
     457             :       end do
     458        1536 :       isize = dry_air_species_num+1
     459        1536 :       icnst = 0 ! N2
     460        1536 :       if(isize > 0) then
     461        1536 :          if(masterproc) then
     462           2 :             write(iulog, *) "Dry air composition ",                           &
     463           2 :                  TRIM(dry_air_species(idx)),                                  &
     464           2 :                  icnst, -1, thermodynamic_active_species_mwi(icnst),          &
     465           2 :                  thermodynamic_active_species_cp(icnst),                      &
     466           4 :                  thermodynamic_active_species_cv(icnst)
     467             :          end if
     468             :       end if
     469             :       !
     470             :       !************************************************************************
     471             :       !
     472             :       ! Add non-dry components of moist air (water vapor and condensates)
     473             :       !
     474             :       !************************************************************************
     475             :       !
     476        1536 :       icnst = dry_air_species_num + 1
     477       10752 :       do idx = 1, water_species_in_air_num
     478       18432 :          select case (TRIM(water_species_in_air(idx)))
     479             :             !
     480             :             ! Q
     481             :             !
     482             :          case('Q')
     483        1536 :             call air_species_info('Q', ix, mw)
     484        1536 :             wv_idx = ix
     485        1536 :             thermodynamic_active_species_idx(icnst) = ix
     486        1536 :             thermodynamic_active_species_cp (icnst) = cpwv
     487        1536 :             thermodynamic_active_species_cv (icnst) = cv3 / mw
     488        1536 :             thermodynamic_active_species_R  (icnst) = rh2o
     489        1536 :             icnst = icnst + 1
     490             :             !
     491             :             ! CLDLIQ
     492             :             !
     493             :          case('CLDLIQ')
     494        1536 :             call air_species_info('CLDLIQ', ix, mw)
     495        1536 :             thermodynamic_active_species_idx(icnst) = ix
     496        1536 :             thermodynamic_active_species_cp (icnst) = cpliq
     497        1536 :             thermodynamic_active_species_cv (icnst) = cpliq
     498        1536 :             liq_num           = liq_num+1
     499        1536 :             liq_idx (liq_num) = ix
     500        1536 :             icnst = icnst + 1
     501        1536 :             has_liq = .true.
     502             :             !
     503             :             ! CLDICE
     504             :             !
     505             :          case('CLDICE')
     506        1536 :             call air_species_info('CLDICE', ix, mw)
     507        1536 :             thermodynamic_active_species_idx(icnst) = ix
     508        1536 :             thermodynamic_active_species_cp (icnst) = cpice
     509        1536 :             thermodynamic_active_species_cv (icnst) = cpice
     510        1536 :             ice_num           = ice_num+1
     511        1536 :             ice_idx(ice_num)  = ix
     512        1536 :             icnst = icnst + 1
     513        1536 :             has_ice = .true.
     514             :             !
     515             :             ! RAINQM
     516             :             !
     517             :          case('RAINQM')
     518        1536 :             call air_species_info('RAINQM', ix, mw)
     519        1536 :             thermodynamic_active_species_idx(icnst) = ix
     520        1536 :             thermodynamic_active_species_cp (icnst) = cpliq
     521        1536 :             thermodynamic_active_species_cv (icnst) = cpliq
     522        1536 :             liq_num           = liq_num+1
     523        1536 :             liq_idx(liq_num)  = ix
     524        1536 :             icnst = icnst + 1
     525        1536 :             has_liq = .true.
     526             :             !
     527             :             ! SNOWQM
     528             :             !
     529             :          case('SNOWQM')
     530        1536 :             call air_species_info('SNOWQM', ix, mw)
     531        1536 :             thermodynamic_active_species_idx(icnst) = ix
     532        1536 :             thermodynamic_active_species_cp (icnst) = cpice
     533        1536 :             thermodynamic_active_species_cv (icnst) = cpice
     534        1536 :             ice_num           = ice_num+1
     535        1536 :             ice_idx(ice_num)  = ix
     536        1536 :             icnst = icnst + 1
     537        1536 :             has_ice = .true.
     538             :             !
     539             :             ! GRAUQM
     540             :             !
     541             :          case('GRAUQM')
     542        1536 :             call air_species_info('GRAUQM', ix, mw)
     543        1536 :             thermodynamic_active_species_idx(icnst) = ix
     544        1536 :             thermodynamic_active_species_cp (icnst) = cpice
     545        1536 :             thermodynamic_active_species_cv (icnst) = cpice
     546        1536 :             ice_num           = ice_num+1
     547        1536 :             ice_idx(ice_num)  = ix
     548        1536 :             icnst = icnst + 1
     549        1536 :             has_ice = .true.
     550             :             !
     551             :             ! If support for more major species is to be included add code here
     552             :             !
     553             :          case default
     554           0 :             write(iulog, *) subname, ' moist air component not found: ',      &
     555           0 :                  water_species_in_air(idx)
     556       18432 :             call endrun(subname//': moist air component not found')
     557             :          end select
     558             :          !
     559             :          !
     560             :          !
     561        9216 :          if (masterproc) then
     562          12 :             write(iulog, *) "Thermodynamic active species ",                  &
     563          24 :                  TRIM(water_species_in_air(idx))
     564          12 :             write(iulog, *) "   global index                           : ",   &
     565          24 :                  icnst-1
     566          12 :             write(iulog, *) "   thermodynamic_active_species_idx       : ",   &
     567          24 :                  thermodynamic_active_species_idx(icnst-1)
     568          12 :             write(iulog, *) "   cp                                     : ",   &
     569          24 :                  thermodynamic_active_species_cp(icnst-1)
     570          12 :             write(iulog, *) "   cv                                     : ",   &
     571          24 :                  thermodynamic_active_species_cv(icnst-1)
     572          12 :             if (has_liq) then
     573           4 :                write(iulog, *) "   register phase (liquid or ice)         :", &
     574           8 :                     " liquid"
     575             :             end if
     576          12 :             if (has_ice) then
     577           6 :                write(iulog, *) "   register phase (liquid or ice)         :", &
     578          12 :                     " ice"
     579             :             end if
     580          12 :             write(iulog, *) "  "
     581             :          end if
     582        9216 :          has_liq = .false.
     583       10752 :          has_ice = .false.
     584             :       end do
     585             : 
     586        4608 :       allocate(thermodynamic_active_species_liq_idx(liq_num), stat=ierr)
     587        1536 :       if (ierr /= 0) then
     588           0 :          call endrun(errstr//"thermodynamic_active_species_liq_idx")
     589             :       end if
     590        3072 :       allocate(thermodynamic_active_species_liq_idx_dycore(liq_num), stat=ierr)
     591        1536 :       if (ierr /= 0) then
     592           0 :          call endrun(errstr//"thermodynamic_active_species_liq_idx_dycore")
     593             :       end if
     594        4608 :       allocate(thermodynamic_active_species_ice_idx(ice_num), stat=ierr)
     595        1536 :       if (ierr /= 0) then
     596           0 :          call endrun(errstr//"thermodynamic_active_species_ice_idx")
     597             :       end if
     598        3072 :       allocate(thermodynamic_active_species_ice_idx_dycore(ice_num), stat=ierr)
     599        1536 :       if (ierr /= 0) then
     600           0 :          call endrun(errstr//"thermodynamic_active_species_ice_idx_dycore")
     601             :       end if
     602             : 
     603        6144 :       thermodynamic_active_species_liq_idx = liq_idx(1:liq_num)
     604        1536 :       thermodynamic_active_species_liq_num = liq_num
     605             : 
     606             :       ! array initialized by the dycore
     607        4608 :       thermodynamic_active_species_liq_idx_dycore = -99
     608             : 
     609        7680 :       thermodynamic_active_species_ice_idx = ice_idx(1:ice_num)
     610        1536 :       thermodynamic_active_species_ice_num = ice_num
     611             : 
     612             :       ! array initialized by the dycore
     613        6144 :       thermodynamic_active_species_ice_idx_dycore = -99
     614             : 
     615        1536 :       if (water_species_in_air_num /= 1 + liq_num+ice_num) then
     616           0 :          write(iulog, '(2a,2(i0,a))') subname,                                &
     617           0 :               "  water_species_in_air_num = ",                                &
     618           0 :               water_species_in_air_num, ", should be ",              &
     619           0 :               (1 + liq_num + ice_num), " (1 + liq_num + ice_num)"
     620           0 :          call endrun(subname//': water_species_in_air_num /= 1+liq_num+ice_num')
     621             :       end if
     622        1536 :       enthalpy_reference_state = 'ice'
     623        1536 :       if (masterproc) then
     624           2 :          write(iulog, *)   'Enthalpy reference state           : ',           &
     625           4 :               TRIM(enthalpy_reference_state)
     626             :       end if
     627        1536 :    end subroutine air_composition_init
     628             : 
     629             :    !===========================================================================
     630             :    !-----------------------------------------------------------------------
     631             :    ! dry_air_composition_update: Update the physics "constants" that vary
     632             :    !-------------------------------------------------------------------------
     633             :    !===========================================================================
     634             : 
     635           0 :    subroutine dry_air_composition_update(mmr, lchnk, ncol, to_dry_factor)
     636             :       use cam_abortutils,  only: endrun
     637             :       !(mmr = dry mixing ratio, if not, use to_dry_factor to convert!)
     638             :       real(r8),           intent(in) :: mmr(:,:,:) ! mixing ratios for species dependent dry air
     639             :       integer,            intent(in) :: lchnk      ! Chunk number
     640             :       integer,            intent(in) :: ncol       ! number of columns
     641             :       real(r8), optional, intent(in) :: to_dry_factor(:,:)
     642             : 
     643           0 :       call get_R_dry(mmr(:ncol, :, :), thermodynamic_active_species_idx, &
     644           0 :            rairv(:ncol, :, lchnk), fact=to_dry_factor)
     645             :       call get_cp_dry(mmr(:ncol,:,:), thermodynamic_active_species_idx,  &
     646           0 :            cpairv(:ncol,:,lchnk), fact=to_dry_factor)
     647             :       call get_mbarv(mmr(:ncol,:,:), thermodynamic_active_species_idx,   &
     648           0 :            mbarv(:ncol,:,lchnk), fact=to_dry_factor)
     649           0 :       cappav(:ncol,:,lchnk) = rairv(:ncol,:,lchnk) / cpairv(:ncol,:,lchnk)
     650           0 :    end subroutine dry_air_composition_update
     651             : 
     652             :    !===========================================================================
     653             :    !---------------------------------------------------------------------------
     654             :    ! water_composition_update: Update generalized cp or cv depending on dycore
     655             :    !---------------------------------------------------------------------------
     656             :    !===========================================================================
     657             : 
     658     2984544 :    subroutine water_composition_update(mmr, lchnk, ncol, vcoord, to_dry_factor)
     659             :       use cam_abortutils,  only: endrun
     660             :       use string_utils,    only: int2str
     661             :       use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure
     662             :       real(r8),           intent(in) :: mmr(:,:,:) ! constituents array
     663             :       integer,            intent(in) :: lchnk      ! Chunk number
     664             :       integer,            intent(in) :: ncol       ! number of columns
     665             :       integer,            intent(in) :: vcoord
     666             :       real(r8), optional, intent(in) :: to_dry_factor(:,:)
     667             : 
     668             :       character(len=*), parameter :: subname = 'water_composition_update'
     669             : 
     670     2984544 :       if (vcoord==vc_dry_pressure) then
     671     2984544 :         call get_cp(mmr(:ncol,:,:),.false.,cp_or_cv_dycore(:ncol,:,lchnk), factor=to_dry_factor,    &
     672     7464456 :              active_species_idx_dycore=thermodynamic_active_species_idx,cpdry=cpairv(:ncol,:,lchnk))
     673           0 :       else if (vcoord==vc_height) then
     674           0 :         call get_R(mmr(:ncol,:,:), thermodynamic_active_species_idx, &
     675           0 :              cp_or_cv_dycore(:ncol,:,lchnk), fact=to_dry_factor, Rdry=rairv(:ncol,:,lchnk))
     676             :         !
     677             :         ! internal energy coefficient for MPAS 
     678             :         ! (equation 92 in Eldred et al. 2023; https://rmets.onlinelibrary.wiley.com/doi/epdf/10.1002/qj.4353)
     679             :         !
     680           0 :         cp_or_cv_dycore(:ncol,:,lchnk)=cp_or_cv_dycore(:ncol,:,lchnk)*&
     681           0 :              (cpairv(:ncol,:,lchnk)-rairv(:ncol,:,lchnk)) /rairv(:ncol,:,lchnk)
     682           0 :       else if (vcoord==vc_moist_pressure) then
     683             :         ! no update needed for moist pressure vcoord
     684             :       else
     685           0 :         call endrun(subname//" vertical coordinate not supported; vcoord="//  int2str(vcoord))
     686             :       end if
     687     2984544 :    end subroutine water_composition_update
     688             : 
     689             :    !===========================================================================
     690             :    !***************************************************************************
     691             :    !
     692             :    ! get_cp_dry: Compute dry air heat capacity under constant pressure
     693             :    !
     694             :    !***************************************************************************
     695             :    !
     696   311688000 :    subroutine get_cp_dry_1hd(tracer, active_species_idx, cp_dry, fact)
     697             :       use cam_abortutils,  only: endrun
     698             :       use string_utils,    only: int2str
     699             :       use physconst,       only: cpair
     700             : 
     701             :       ! Dummy arguments
     702             :       ! tracer: tracer array
     703             :       real(r8),           intent(in)  :: tracer(:,:,:)
     704             :       integer,            intent(in)  :: active_species_idx(:)
     705             :       ! fact: optional dry pressure level thickness
     706             :       real(r8), optional, intent(in)  :: fact(:,:)
     707             :       ! cp_dry: dry air heat capacity under constant pressure
     708             :       real(r8),           intent(out) :: cp_dry(:,:)
     709             : 
     710             :       ! Local variables
     711             :       integer  :: idx, kdx , m_cnst, qdx
     712             :       ! factor: dry pressure level thickness
     713   623376000 :       real(r8) :: factor(SIZE(cp_dry, 1), SIZE(cp_dry, 2))
     714   623376000 :       real(r8) :: residual(SIZE(cp_dry, 1), SIZE(cp_dry, 2))
     715             :       real(r8) :: mmr
     716             :       character(len=*), parameter :: subname = 'get_cp_dry_1hd: '
     717             : 
     718   311688000 :       if (dry_air_species_num == 0) then
     719             :          ! dry air heat capacity not species dependent
     720 >14524*10^7 :          cp_dry = cpair
     721             :       else
     722             :          ! dry air heat capacity is species dependent
     723           0 :          if (present(fact)) then
     724           0 :             if (SIZE(fact, 1) /= SIZE(factor, 1)) then
     725             :                call endrun(subname//"SIZE mismatch in dimension 1 "//         &
     726           0 :                     int2str(SIZE(fact, 1))//' /= '//int2str(SIZE(factor, 1)))
     727             :             end if
     728           0 :             if (SIZE(fact, 2) /= SIZE(factor, 2)) then
     729             :                call endrun(subname//"SIZE mismatch in dimension 2 "//         &
     730           0 :                     int2str(SIZE(fact, 2))//' /= '//int2str(SIZE(factor, 2)))
     731             :             end if
     732           0 :             factor = fact(:,:)
     733             :          else
     734           0 :             factor = 1.0_r8
     735             :          end if
     736             : 
     737           0 :          cp_dry = 0.0_r8
     738           0 :          residual = 1.0_r8
     739           0 :          do qdx = 1, dry_air_species_num
     740           0 :             m_cnst = active_species_idx(qdx)
     741           0 :             do kdx = 1, SIZE(cp_dry, 2)
     742           0 :                do idx = 1, SIZE(cp_dry, 1)
     743           0 :                   mmr = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
     744           0 :                   cp_dry(idx, kdx) = cp_dry(idx, kdx) +             &
     745           0 :                        (thermodynamic_active_species_cp(qdx) * mmr)
     746           0 :                   residual(idx, kdx) = residual(idx, kdx) - mmr
     747             :                end do
     748             :             end do
     749             :          end do
     750           0 :          qdx = 0 ! N2
     751           0 :          do kdx = 1, SIZE(cp_dry, 2)
     752           0 :             do idx = 1, SIZE(cp_dry, 1)
     753           0 :                cp_dry(idx, kdx) = cp_dry(idx, kdx) +                          &
     754           0 :                     (thermodynamic_active_species_cp(qdx) * residual(idx, kdx))
     755             :             end do
     756             :          end do
     757             :       end if
     758   311688000 :    end subroutine get_cp_dry_1hd
     759             : 
     760             :    !===========================================================================
     761             : 
     762    77922000 :    subroutine get_cp_dry_2hd(tracer, active_species_idx, cp_dry, fact)
     763             :       ! Version of get_cp_dry for arrays that have a second horizontal index
     764             : 
     765             :       ! Dummy arguments
     766             :       ! tracer: tracer array
     767             :       real(r8),           intent(in)  :: tracer(:,:,:,:)
     768             :       integer,            intent(in)  :: active_species_idx(:)
     769             :       ! fact:        optional dry pressure level thickness
     770             :       real(r8), optional, intent(in)  :: fact(:,:,:)
     771             :       ! cp_dry: dry air heat capacity under constant pressure
     772             :       real(r8),           intent(out) :: cp_dry(:,:,:)
     773             : 
     774             :       ! Local variable
     775             :       integer  :: jdx
     776             : 
     777   389610000 :       do jdx = 1, SIZE(cp_dry, 2)
     778   389610000 :          if (present(fact)) then
     779             :             call get_cp_dry(tracer(:,jdx,:,:), active_species_idx,            &
     780           0 :                  cp_dry(:,jdx,:), fact=fact(:,jdx,:))
     781             :          else
     782             :             call get_cp_dry(tracer(:,jdx,:,:), active_species_idx,            &
     783   311688000 :                  cp_dry(:,jdx,:))
     784             :          end if
     785             :       end do
     786             : 
     787    77922000 :    end subroutine get_cp_dry_2hd
     788             : 
     789             :    !===========================================================================
     790             :    !
     791             :    !***************************************************************************
     792             :    !
     793             :    ! get_cp: Compute generalized heat capacity at constant pressure
     794             :    !
     795             :    !***************************************************************************
     796             :    !
     797    65322144 :    subroutine get_cp_1hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry)
     798             :       use cam_abortutils,  only: endrun
     799             :       use string_utils,    only: int2str
     800             : 
     801             :       ! Dummy arguments
     802             :       ! tracer: Tracer array
     803             :       !
     804             :       ! factor not present then tracer must be dry mixing ratio
     805             :       ! if factor present tracer*factor must be dry mixing ratio
     806             :       !
     807             :       real(r8),           intent(in)  :: tracer(:,:,:)
     808             :       ! inv_cp: output inverse cp instead of cp
     809             :       logical,            intent(in)  :: inv_cp
     810             :       real(r8),           intent(out) :: cp(:,:)
     811             :       ! dp: if provided then tracer is mass not mixing ratio
     812             :       real(r8), optional, intent(in)  :: factor(:,:)
     813             :       ! active_species_idx_dycore: array of indices for index of
     814             :       !    thermodynamic active species in dycore tracer array
     815             :       !    (if different from physics index)
     816             :       integer, optional,  intent(in)  :: active_species_idx_dycore(:)
     817             :       real(r8),optional,  intent(in)  :: cpdry(:,:)
     818             : 
     819             :       ! LOCAL VARIABLES
     820             :       integer  :: qdx, itrac
     821   130644288 :       real(r8) :: sum_species(SIZE(cp, 1), SIZE(cp, 2))
     822   130644288 :       real(r8) :: sum_cp(SIZE(cp, 1), SIZE(cp, 2))
     823   130644288 :       real(r8) :: factor_local(SIZE(cp, 1), SIZE(cp, 2))
     824   130644288 :       integer  :: idx_local(thermodynamic_active_species_num)
     825             :       character(LEN=*), parameter :: subname = 'get_cp_1hd: '
     826             : 
     827    65322144 :       if (present(active_species_idx_dycore)) then
     828    65322144 :          if (SIZE(active_species_idx_dycore) /=                               &
     829             :               thermodynamic_active_species_num) then
     830             :             call endrun(subname//"SIZE mismatch "//                           &
     831             :                  int2str(SIZE(active_species_idx_dycore))//' /= '//           &
     832           0 :                  int2str(thermodynamic_active_species_num))
     833             :         end if
     834   522577152 :          idx_local = active_species_idx_dycore
     835             :       else
     836           0 :          idx_local = thermodynamic_active_species_idx
     837             :       end if
     838             : 
     839    65322144 :       if (present(factor)) then
     840  2315495520 :          factor_local = factor
     841             :       else
     842 31372949592 :          factor_local = 1.0_r8
     843             :       end if
     844             : 
     845 33686955936 :       sum_species = 1.0_r8 ! all dry air species sum to 1
     846   457255008 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
     847   391932864 :         itrac = idx_local(qdx)
     848 >20218*10^7 :         sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:))
     849             :       end do
     850             : 
     851    65322144 :       if (dry_air_species_num == 0) then
     852 33686955936 :          sum_cp = thermodynamic_active_species_cp(0)
     853           0 :       else if (present(cpdry)) then
     854             :         !
     855             :         ! if cpdry is known don't recompute
     856             :         !
     857           0 :         sum_cp = cpdry
     858             :       else
     859           0 :          call get_cp_dry(tracer, idx_local, sum_cp, fact=factor_local)
     860             :       end if
     861   457255008 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
     862   391932864 :         itrac = idx_local(qdx)
     863             :         sum_cp(:,:) = sum_cp(:,:)+ &
     864 >20218*10^7 :              thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac)* factor_local(:,:)
     865             :       end do
     866    65322144 :       if (inv_cp) then
     867 29111659200 :         cp = sum_species / sum_cp
     868             :       else
     869  4640618880 :         cp = sum_cp / sum_species
     870             :       end if
     871    65322144 :     end subroutine get_cp_1hd
     872             : 
     873             :    !===========================================================================
     874             : 
     875    15584400 :    subroutine get_cp_2hd(tracer, inv_cp, cp, factor, active_species_idx_dycore, cpdry)
     876             :       ! Version of get_cp for arrays that have a second horizontal index
     877             :       use cam_abortutils, only: endrun
     878             :       use string_utils,   only: int2str
     879             : 
     880             :       ! Dummy arguments
     881             :       ! tracer: Tracer array
     882             :       !
     883             :       real(r8),           intent(in)  :: tracer(:,:,:,:)
     884             :       ! inv_cp: output inverse cp instead of cp
     885             :       logical,            intent(in)  :: inv_cp
     886             :       real(r8),           intent(out) :: cp(:,:,:)
     887             :       real(r8), optional, intent(in)  :: factor(:,:,:)
     888             :       real(r8), optional, intent(in)  :: cpdry(:,:,:)
     889             : 
     890             :       ! active_species_idx_dycore: array of indicies for index of
     891             :       !    thermodynamic active species in dycore tracer array
     892             :       !    (if different from physics index)
     893             :       integer, optional,  intent(in)  :: active_species_idx_dycore(:)
     894             : 
     895             :       ! Local variables
     896             :       integer  :: jdx
     897             :       integer  :: idx_local(thermodynamic_active_species_num)
     898             :       character(len=*), parameter :: subname = 'get_cp_2hd: '
     899             : 
     900    77922000 :       do jdx = 1, SIZE(cp, 2)
     901    77922000 :          if (present(factor).and.present(cpdry)) then
     902             :             call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
     903           0 :                  factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:))
     904    62337600 :          else if (present(factor)) then
     905             :             call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
     906           0 :                  factor=factor(:, jdx, :), active_species_idx_dycore=active_species_idx_dycore)
     907    62337600 :          else if (present(cpdry)) then
     908             :             call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
     909           0 :                  active_species_idx_dycore=active_species_idx_dycore, cpdry=cpdry(:,jdx,:))
     910             :          else
     911             :             call get_cp(tracer(:, jdx, :, :), inv_cp, cp(:, jdx, :),&
     912    62337600 :                  active_species_idx_dycore=active_species_idx_dycore)
     913             :          end if
     914             :       end do
     915             : 
     916    15584400 :    end subroutine get_cp_2hd
     917             : 
     918             :    !===========================================================================
     919             : 
     920             :    !***************************************************************************
     921             :    !
     922             :    ! get_R_dry: Compute generalized dry air gas constant R
     923             :    !
     924             :    !***************************************************************************
     925             :    !
     926   623376000 :    subroutine get_R_dry_1hd(tracer, active_species_idx_dycore, R_dry, fact)
     927             :       use physconst,       only: rair
     928             : 
     929             :       ! tracer: tracer array
     930             :       real(r8),           intent(in)  :: tracer(:, :, :)
     931             :       ! active_species_idx_dycore: index of active species in tracer
     932             :       integer,            intent(in)  :: active_species_idx_dycore(:)
     933             :       ! R_dry: dry air R
     934             :       real(r8),           intent(out) :: R_dry(:, :)
     935             :       ! fact:   optional factor for converting tracer to dry mixing ratio
     936             :       real(r8), optional, intent(in)  :: fact(:, :)
     937             : 
     938             :       ! Local variables
     939             :       integer  :: idx, kdx, m_cnst, qdx
     940  1246752000 :       real(r8) :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
     941  1246752000 :       real(r8) :: residual(SIZE(R_dry, 1), SIZE(R_dry, 2))
     942             :       real(r8) :: mmr
     943             : 
     944   623376000 :       if (dry_air_species_num == 0) then
     945             :          !
     946             :          ! dry air not species dependent
     947             :          !
     948 >29049*10^7 :          R_dry = rair
     949             :       else
     950           0 :          if (present(fact)) then
     951           0 :             factor = fact(:,:)
     952             :          else
     953           0 :             factor = 1.0_r8
     954             :          end if
     955             : 
     956           0 :          R_dry = 0.0_r8
     957           0 :          residual = 1.0_r8
     958           0 :          do qdx = 1, dry_air_species_num
     959           0 :             m_cnst = active_species_idx_dycore(qdx)
     960           0 :             do kdx = 1, SIZE(R_dry, 2)
     961           0 :                do idx = 1, SIZE(R_dry, 1)
     962           0 :                   mmr = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
     963           0 :                   R_dry(idx, kdx) = R_dry(idx, kdx) +                         &
     964           0 :                        (thermodynamic_active_species_R(qdx) * mmr)
     965           0 :                   residual(idx, kdx) = residual(idx, kdx) - mmr
     966             :                end do
     967             :             end do
     968             :          end do
     969             :          !
     970             :          ! N2 derived from the others
     971             :          !
     972           0 :          qdx = 0
     973           0 :          do kdx = 1, SIZE(R_dry, 2)
     974           0 :             do idx = 1, SIZE(R_dry, 1)
     975           0 :                R_dry(idx, kdx) = R_dry(idx, kdx) +                            &
     976           0 :                     (thermodynamic_active_species_R(qdx) * residual(idx, kdx))
     977             :             end do
     978             :          end do
     979             :       end if
     980   623376000 :    end subroutine get_R_dry_1hd
     981             : 
     982             :    !===========================================================================
     983             : 
     984    77922000 :    subroutine get_R_dry_2hd(tracer, active_species_idx_dycore, R_dry, fact)
     985             :       ! Version of get_R_dry for arrays that have a second horizontal index
     986             : 
     987             :       ! tracer: tracer array
     988             :       real(r8),           intent(in)  :: tracer(:, :, :, :)
     989             :       ! active_species_idx_dycore: index of active species in tracer
     990             :       integer,            intent(in)  :: active_species_idx_dycore(:)
     991             :       ! R_dry: dry air R
     992             :       real(r8),           intent(out) :: R_dry(:, :, :)
     993             :       ! fact:   optional factor for converting tracer to dry mixing ratio
     994             :       real(r8), optional, intent(in)  :: fact(:, :, :)
     995             : 
     996             :       ! Local variable
     997             :       integer  :: jdx
     998             : 
     999   389610000 :       do jdx = 1, SIZE(tracer, 2)
    1000   389610000 :          if (present(fact)) then
    1001             :             call get_R_dry(tracer(:, jdx, :, :), active_species_idx_dycore,      &
    1002           0 :                  R_dry(:, jdx, :), fact=fact(:, jdx, :))
    1003             :          else
    1004             :             call get_R_dry(tracer(:, jdx, :, :), active_species_idx_dycore,      &
    1005   311688000 :                  R_dry(:, jdx, :))
    1006             :          end if
    1007             :       end do
    1008             : 
    1009    77922000 :    end subroutine get_R_dry_2hd
    1010             : 
    1011             :    !===========================================================================
    1012             :    !
    1013             :    !***************************************************************************
    1014             :    !
    1015             :    ! get_R: Compute generalized R
    1016             :    !        This code (both 1hd and 2hd) is currently unused and untested
    1017             :    !
    1018             :    !***************************************************************************
    1019             :    !
    1020           0 :    subroutine get_R_1hd(tracer, active_species_idx, R, fact, Rdry)
    1021             :       use cam_abortutils,  only: endrun
    1022             :       use string_utils,    only: int2str
    1023             :       use physconst,       only: rair
    1024             : 
    1025             :       ! Dummy arguments
    1026             :       ! tracer: !tracer array
    1027             :       real(r8), intent(in)  :: tracer(:, :, :)
    1028             :       ! active_species_idx: index of active species in tracer
    1029             :       integer,  intent(in)  :: active_species_idx(:)
    1030             :       ! R: generalized gas constant
    1031             :       real(r8), intent(out) :: R(:, :)
    1032             :       ! fact: optional factor for converting tracer to dry mixing ratio
    1033             :       real(r8), optional, intent(in) :: fact(:, :)
    1034             :       real(r8), optional, intent(in) :: Rdry(:, :)
    1035             : 
    1036             :       ! Local variables
    1037             :       integer  :: qdx, itrac
    1038           0 :       real(r8) :: factor(SIZE(tracer, 1), SIZE(tracer, 2))
    1039           0 :       real(r8) :: sum_species(SIZE(R, 1), SIZE(R, 2))
    1040           0 :       integer  :: idx_local(thermodynamic_active_species_num)
    1041             : 
    1042             :       character(len=*), parameter :: subname = 'get_R_1hd: '
    1043             : 
    1044           0 :       if (present(fact)) then
    1045           0 :          if (SIZE(fact, 1) /= SIZE(factor, 1)) then
    1046             :             call endrun(subname//"SIZE mismatch in dimension 1 "//            &
    1047           0 :                  int2str(SIZE(fact, 1))//' /= '//int2str(SIZE(factor, 1)))
    1048             :          end if
    1049           0 :          if (SIZE(fact, 2) /= SIZE(factor, 2)) then
    1050             :             call endrun(subname//"SIZE mismatch in dimension 2 "//            &
    1051           0 :                  int2str(SIZE(fact, 2))//' /= '//int2str(SIZE(factor, 2)))
    1052             :          end if
    1053           0 :          factor = fact(:,:)
    1054             :       else
    1055           0 :          factor = 1.0_r8
    1056             :       end if
    1057             : 
    1058           0 :       if (dry_air_species_num == 0) then
    1059           0 :         R = rair
    1060           0 :       else if (present(Rdry)) then
    1061           0 :         R = Rdry
    1062             :       else
    1063           0 :         call get_R_dry(tracer, active_species_idx, R, fact=factor)
    1064             :       end if
    1065             : 
    1066           0 :       idx_local = active_species_idx
    1067           0 :       sum_species = 1.0_r8 ! all dry air species sum to 1
    1068           0 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
    1069           0 :          itrac = idx_local(qdx)
    1070             :          sum_species(:,:) = sum_species(:,:) +                            &
    1071           0 :               (tracer(:,:,itrac) * factor(:,:))
    1072             :       end do
    1073           0 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
    1074           0 :          itrac = idx_local(qdx)
    1075             :          R(:,:) = R(:,:) +                                                &
    1076           0 :               (thermodynamic_active_species_R(qdx) * tracer(:,:,itrac) *    &
    1077           0 :               factor(:,:))
    1078             :       end do
    1079           0 :       R = R / sum_species
    1080           0 :    end subroutine get_R_1hd
    1081             : 
    1082             :    !===========================================================================
    1083             : 
    1084           0 :    subroutine get_R_2hd(tracer, active_species_idx, R, fact)
    1085             : 
    1086             :       ! Dummy arguments
    1087             :       ! tracer: !tracer array
    1088             :       real(r8), intent(in)  :: tracer(:, :, :, :)
    1089             :       ! active_species_idx: index of active species in tracer
    1090             :       integer,  intent(in)  :: active_species_idx(:)
    1091             :       ! R: generalized gas constant
    1092             :       real(r8), intent(out) :: R(:, :, :)
    1093             :       ! fact: optional factor for converting tracer to dry mixing ratio
    1094             :       real(r8), optional, intent(in) :: fact(:, :, :)
    1095             : 
    1096             :       ! Local variable
    1097             :       integer  :: jdx
    1098             : 
    1099           0 :       do jdx = 1, SIZE(tracer, 2)
    1100           0 :          if (present(fact)) then
    1101             :             call get_R(tracer(:, jdx, :, :), active_species_idx,          &
    1102           0 :                  R(:, jdx, :), fact=fact(:, jdx, :))
    1103             :          else
    1104             :             call get_R(tracer(:, jdx, :, :), active_species_idx,          &
    1105           0 :                  R(:, jdx, :))
    1106             :          end if
    1107             :       end do
    1108             : 
    1109           0 :    end subroutine get_R_2hd
    1110             : 
    1111             :    !===========================================================================
    1112             : 
    1113             :    !*************************************************************************************************************************
    1114             :    !
    1115             :    ! compute molecular weight dry air
    1116             :    !
    1117             :    !*************************************************************************************************************************
    1118             :    !
    1119           0 :    subroutine get_mbarv_1hd(tracer, active_species_idx, mbarv_in, fact)
    1120             :      use physconst,        only: mwdry
    1121             :      real(r8), intent(in)  :: tracer(:,:,:)                      !tracer array
    1122             :      integer,  intent(in)  :: active_species_idx(:)              !index of active species in tracer
    1123             :      real(r8), intent(out) :: mbarv_in(:,:)                        !molecular weight of dry air
    1124             :      real(r8), optional, intent(in) :: fact(:,:)                 !factor for converting tracer to dry mixing ratio
    1125             : 
    1126             :      integer :: idx, kdx, m_cnst, qdx
    1127           0 :      real(r8):: factor(SIZE(mbarv_in, 1), SIZE(mbarv_in, 2))
    1128           0 :      real(r8):: residual(SIZE(tracer, 1), SIZE(mbarv_in, 2))
    1129             :      real(r8):: mm
    1130             :      !
    1131             :      ! dry air not species dependent
    1132             :      !
    1133           0 :      if (dry_air_species_num==0) then
    1134           0 :        mbarv_in = mwdry
    1135             :      else
    1136           0 :        if (present(fact)) then
    1137           0 :          factor(:,:) = fact(:,:)
    1138             :        else
    1139           0 :          factor(:,:) = 1.0_r8
    1140             :        endif
    1141             : 
    1142           0 :        mbarv_in = 0.0_r8
    1143           0 :        residual = 1.0_r8
    1144           0 :        do qdx = 1, dry_air_species_num
    1145           0 :          m_cnst = active_species_idx(qdx)
    1146           0 :          do kdx = 1, SIZE(mbarv_in, 2)
    1147           0 :            do idx = 1, SIZE(mbarv_in, 1)
    1148           0 :              mm = tracer(idx, kdx, m_cnst) * factor(idx, kdx)
    1149           0 :              mbarv_in(idx, kdx) = mbarv_in(idx, kdx) + thermodynamic_active_species_mwi(qdx) * mm
    1150           0 :              residual(idx, kdx) = residual(idx, kdx) - mm
    1151             :            end do
    1152             :          end do
    1153             :        end do
    1154           0 :        qdx = 0 ! N2
    1155           0 :        do kdx = 1, SIZE(mbarv_in, 2)
    1156           0 :          do idx = 1, SIZE(mbarv_in, 1)
    1157           0 :            mbarv_in(idx, kdx) = mbarv_in(idx, kdx) + thermodynamic_active_species_mwi(qdx) * residual(idx, kdx)
    1158             :          end do
    1159             :        end do
    1160           0 :        mbarv_in(:,:) = 1.0_r8 / mbarv_in(:,:)
    1161             :      end if
    1162           0 :    end subroutine get_mbarv_1hd
    1163             : 
    1164             :    !===========================================================================
    1165             : 
    1166        9216 :    subroutine air_species_info(name, index, molec_weight, caller)
    1167             :       use cam_abortutils, only: endrun
    1168             :       use cam_logfile,    only: iulog
    1169             :       use constituents,   only: cnst_get_ind, cnst_mw
    1170             :       ! Find the constituent index of <name> and return it in
    1171             :       !    <index>. Return the constituent molecular weight in
    1172             :       !    <molec_weight>
    1173             : 
    1174             :       ! Dummy arguments
    1175             :       character(len=*),           intent(in)  :: name
    1176             :       integer,                    intent(out) :: index
    1177             :       real(r8),                   intent(out) :: molec_weight
    1178             :       character(len=*), optional, intent(in)  :: caller
    1179             :       ! Local parameter
    1180             :       character(len=*), parameter :: subname = 'air_species_info: '
    1181             : 
    1182        9216 :       call cnst_get_ind(trim(name), index, abort=.false.)
    1183        9216 :       if (index < 1) then
    1184           0 :          if (present(caller)) then
    1185           0 :             write(iulog, *) trim(caller), ": air component not found, '", &
    1186           0 :                  trim(name), "'"
    1187             :             call endrun(trim(caller)//": air component not found, '"//    &
    1188           0 :                  trim(name)//"'")
    1189             :          else
    1190           0 :             write(iulog, *) subname, "air component not found, '",        &
    1191           0 :                  trim(name), "'"
    1192             :             call endrun(subname//"air component not found, '"//           &
    1193           0 :                  trim(name)//"'")
    1194             :          end if
    1195             :       else
    1196        9216 :          molec_weight = cnst_mw(index)
    1197             :       end if
    1198             : 
    1199        9216 :    end subroutine air_species_info
    1200             : 
    1201             : 
    1202             : end module air_composition

Generated by: LCOV version 1.14