LCOV - code coverage report
Current view: | top level - utils - air_composition.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 235 | 467 | 50.3 %
Date: | 2025-01-13 21:54:50 |  | 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        4608 :             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           8 :          do indx = 1, water_species_in_air_num
         221           8 :             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        6144 :       thermodynamic_active_species_idx        = -HUGE(1)
         343        6144 :       thermodynamic_active_species_idx_dycore = -HUGE(1)
         344        7680 :       thermodynamic_active_species_cp         = 0.0_r8
         345        7680 :       thermodynamic_active_species_cv         = 0.0_r8
         346        7680 :       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     2744592 :       cpairv(:pcols, :pver, begchunk:endchunk) = cpair
         354     2744592 :       rairv(:pcols,  :pver, begchunk:endchunk) = rair
         355     2744592 :       cappav(:pcols, :pver, begchunk:endchunk) = rair / cpair
         356     2744592 :       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        6144 :       do idx = 1, water_species_in_air_num
         478        9216 :          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           0 :             call air_species_info('RAINQM', ix, mw)
         519           0 :             thermodynamic_active_species_idx(icnst) = ix
         520           0 :             thermodynamic_active_species_cp (icnst) = cpliq
         521           0 :             thermodynamic_active_species_cv (icnst) = cpliq
         522           0 :             liq_num           = liq_num+1
         523           0 :             liq_idx(liq_num)  = ix
         524           0 :             icnst = icnst + 1
         525           0 :             has_liq = .true.
         526             :             !
         527             :             ! SNOWQM
         528             :             !
         529             :          case('SNOWQM')
         530           0 :             call air_species_info('SNOWQM', ix, mw)
         531           0 :             thermodynamic_active_species_idx(icnst) = ix
         532           0 :             thermodynamic_active_species_cp (icnst) = cpice
         533           0 :             thermodynamic_active_species_cv (icnst) = cpice
         534           0 :             ice_num           = ice_num+1
         535           0 :             ice_idx(ice_num)  = ix
         536           0 :             icnst = icnst + 1
         537           0 :             has_ice = .true.
         538             :             !
         539             :             ! GRAUQM
         540             :             !
         541             :          case('GRAUQM')
         542           0 :             call air_species_info('GRAUQM', ix, mw)
         543           0 :             thermodynamic_active_species_idx(icnst) = ix
         544           0 :             thermodynamic_active_species_cp (icnst) = cpice
         545           0 :             thermodynamic_active_species_cv (icnst) = cpice
         546           0 :             ice_num           = ice_num+1
         547           0 :             ice_idx(ice_num)  = ix
         548           0 :             icnst = icnst + 1
         549           0 :             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        9216 :             call endrun(subname//': moist air component not found')
         557             :          end select
         558             :          !
         559             :          !
         560             :          !
         561        4608 :          if (masterproc) then
         562           6 :             write(iulog, *) "Thermodynamic active species ",                  &
         563          12 :                  TRIM(water_species_in_air(idx))
         564           6 :             write(iulog, *) "   global index                           : ",   &
         565          12 :                  icnst-1
         566           6 :             write(iulog, *) "   thermodynamic_active_species_idx       : ",   &
         567          12 :                  thermodynamic_active_species_idx(icnst-1)
         568           6 :             write(iulog, *) "   cp                                     : ",   &
         569          12 :                  thermodynamic_active_species_cp(icnst-1)
         570           6 :             write(iulog, *) "   cv                                     : ",   &
         571          12 :                  thermodynamic_active_species_cv(icnst-1)
         572           6 :             if (has_liq) then
         573           2 :                write(iulog, *) "   register phase (liquid or ice)         :", &
         574           4 :                     " liquid"
         575             :             end if
         576           6 :             if (has_ice) then
         577           2 :                write(iulog, *) "   register phase (liquid or ice)         :", &
         578           4 :                     " ice"
         579             :             end if
         580           6 :             write(iulog, *) "  "
         581             :          end if
         582        4608 :          has_liq = .false.
         583        6144 :          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        4608 :       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        3072 :       thermodynamic_active_species_liq_idx_dycore = -99
         608             :
         609        4608 :       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        3072 :       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 40831128000 :          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   326610720 :          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   649488528 :          factor_local = factor
         841             :       else
         842  8816919336 :          factor_local = 1.0_r8
         843             :       end if
         844             :
         845  9464918688 :       sum_species = 1.0_r8 ! all dry air species sum to 1
         846   261288576 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         847   195966432 :         itrac = idx_local(qdx)
         848 28460078208 :         sum_species(:,:) = sum_species(:,:) + (tracer(:,:,itrac) * factor_local(:,:))
         849             :       end do
         850             :
         851    65322144 :       if (dry_air_species_num == 0) then
         852  9464918688 :          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   261288576 :       do qdx = dry_air_species_num + 1, thermodynamic_active_species_num
         862   195966432 :         itrac = idx_local(qdx)
         863             :         sum_cp(:,:) = sum_cp(:,:)+ &
         864 28460078208 :              thermodynamic_active_species_cp(qdx) * tracer(:,:,itrac)* factor_local(:,:)
         865             :       end do
         866    65322144 :       if (inv_cp) then
         867  8228563200 :         cp = sum_species / sum_cp
         868             :       else
         869  1301677632 :         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 81662256000 :          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        4608 :    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        4608 :       call cnst_get_ind(trim(name), index, abort=.false.)
        1183        4608 :       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        4608 :          molec_weight = cnst_mw(index)
        1197             :       end if
        1198             :
        1199        4608 :    end subroutine air_species_info
        1200             :
        1201             :
        1202             : end module air_composition

Generated by: LCOV version 1.14