From b47246dedfa7c377ac2bfb060e2d465beb931f2a Mon Sep 17 00:00:00 2001 From: Jesse Nusbaumer Date: Thu, 22 Jan 2026 09:10:38 -0700 Subject: [PATCH 1/2] Update CAM4 physics suite for SIMA aquaplanet configurations (#349) Originator(s): nusbaume Description (include issue title and the keyword ['closes', 'fixes', 'resolves'] and issue number): Add namelist modifications and additional helper schemes in order to run the full CAM4 physics suite in CAM-SIMA in aquaplanet mode. List all namelist files that were added or changed: A schemes/radiation_utils/radiative_gas_concentrations_namelist.xml - New namelist variable to set the concentrations of radiatively active gases. M schemes/cloud_fraction/compute_cloud_fraction_namelist.xml M schemes/gravity_wave_drag/gw_common_namelist.xml M schemes/radiation_utils/solar_irradiance_data_namelist.xml M schemes/rrtmgp/rrtmgp_constituents_namelist.xml M schemes/zhang_mcfarlane/zm_conv_options_namelist.xml - Add new namelist values for CAM4 and/or aquaplanet CAM-SIMA configurations. List all files eliminated and why: N/A List all files added and what they do: A schemes/radiation_utils/radiative_gas_concentrations.F90 A schemes/radiation_utils/radiative_gas_concentrations.meta - New scheme which sets the global concentrations of radiatively active gases using hardcoded namelist values. A schemes/utilities/set_surface_coupling_vars.F90 A schemes/utilities/set_surface_coupling_vars.meta - New scheme which sets the outgoing (i.e. `cam_out`) coupler variables. List all existing files that have been modified, and describe the changes: (Helpful git command: `git diff --name-status development...`) M schemes/rasch_kristjansson/prognostic_cloud_water.F90 - Move subroutine call outside of loop to improve performance (found during CAM4 aquaplanet debugging). M suites/suite_cam4.xml - Update CAM4 suite to allow for a complete CAM4 physics run in CAM-SIMA when in aquaplanet mode. M suites/suite_cam7.xml - Add commented-out 'set_surface_coupling_vars' scheme, as it will be needed for an eventual CAM7 simulation in CAM-SIMA. List all automated tests that failed, as well as an explanation for why they weren't fixed: Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc? Yes, at least for full physics CAM4 aquaplanet simulations. This PR introduces new schemes and namelist values to the CAM4 physics suite, which is what allows this configuration to run in CAM-SIMA. If yes to the above question, describe how this code was validated with the new/modified features: Performed a short, three month run and examined the model output, at least for some generic state variables (U,V,T,Q,PS). All of the results seemed reasonable, and the model successfully ran to completion. --------- Co-authored-by: peverwhee Co-authored-by: Cheryl Craig Co-authored-by: Courtney Peverley Co-authored-by: John Truesdale Co-authored-by: Haipeng Lin Co-authored-by: Jian Sun Co-authored-by: Haipeng Lin --- .../compute_cloud_fraction_namelist.xml | 1 + .../gravity_wave_drag/gw_common_namelist.xml | 4 + ...prescribe_radiative_gas_concentrations.F90 | 184 ++++++++++++++ ...rescribe_radiative_gas_concentrations.meta | 56 ++++ ..._radiative_gas_concentrations_namelist.xml | 157 ++++++++++++ .../solar_irradiance_data_namelist.xml | 7 +- .../prognostic_cloud_water.F90 | 6 +- .../rrtmgp/rrtmgp_constituents_namelist.xml | 19 +- .../utilities/set_surface_coupling_vars.F90 | 201 +++++++++++++++ .../utilities/set_surface_coupling_vars.meta | 240 ++++++++++++++++++ .../vertical_diffusion_sponge_layer.meta | 12 + .../zm_conv_options_namelist.xml | 4 + suites/suite_cam4.xml | 34 ++- suites/suite_cam7.xml | 8 +- 14 files changed, 916 insertions(+), 17 deletions(-) create mode 100644 schemes/radiation_utils/prescribe_radiative_gas_concentrations.F90 create mode 100644 schemes/radiation_utils/prescribe_radiative_gas_concentrations.meta create mode 100644 schemes/radiation_utils/prescribe_radiative_gas_concentrations_namelist.xml create mode 100644 schemes/utilities/set_surface_coupling_vars.F90 create mode 100644 schemes/utilities/set_surface_coupling_vars.meta diff --git a/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml index 18c511bd..f91074fe 100644 --- a/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml +++ b/schemes/cloud_fraction/compute_cloud_fraction_namelist.xml @@ -133,6 +133,7 @@ 0.900D0 0.8975D0 0.8975D0 + 0.910D0 diff --git a/schemes/gravity_wave_drag/gw_common_namelist.xml b/schemes/gravity_wave_drag/gw_common_namelist.xml index 167d9543..1cc191ab 100644 --- a/schemes/gravity_wave_drag/gw_common_namelist.xml +++ b/schemes/gravity_wave_drag/gw_common_namelist.xml @@ -32,6 +32,7 @@ 2.5D0 + 0.D0 @@ -47,6 +48,7 @@ 10 + 0 @@ -63,6 +65,7 @@ 2.0D0 + 0.D0 @@ -94,6 +97,7 @@ .true. + .false. diff --git a/schemes/radiation_utils/prescribe_radiative_gas_concentrations.F90 b/schemes/radiation_utils/prescribe_radiative_gas_concentrations.F90 new file mode 100644 index 00000000..3fe7411a --- /dev/null +++ b/schemes/radiation_utils/prescribe_radiative_gas_concentrations.F90 @@ -0,0 +1,184 @@ +!------------------------------------------------------------------------------- +! This module uses namelist variable to set prescribed concentrations +! for radiatively-active gases. + +! Eventually this module should be replaced with a more comprehensive atmospheric +! composition/chemistry system, but is fine to use for now when running low-top, +! non-exoplanet CAM-SIMA configurations with minimal chemistry. +!------------------------------------------------------------------------------- +module prescribe_radiative_gas_concentrations + + implicit none + + private + public :: prescribe_radiative_gas_concentrations_init + + +!------------------------------------------------------------------------------- +contains +!------------------------------------------------------------------------------- + +!> \section arg_table_prescribe_radiative_gas_concentrations_init Argument Table +!! \htmlinclude prescribe_radiative_gas_concentrations_init.html +!! + subroutine prescribe_radiative_gas_concentrations_init(ch4_vmr, co2_vmr, cfc11_vmr, & + cfc12_vmr, n2o_vmr, const_array, & + errmsg, errcode) + + ! Use statements + use ccpp_constituent_prop_mod, only: int_unassigned + use ccpp_scheme_utils, only: ccpp_constituent_index + use ccpp_kinds, only: kind_phys + + ! Use statement from RRTMGP, + ! which should hopefully be replaced once + ! molar masses for these species are included + ! in the constituents properties themselves: + use radiation_utils, only: get_molar_mass_ratio + + + ! Input arguments + real(kind_phys), intent(in) :: ch4_vmr + real(kind_phys), intent(in) :: co2_vmr + real(kind_phys), intent(in) :: cfc11_vmr + real(kind_phys), intent(in) :: cfc12_vmr + real(kind_phys), intent(in) :: n2o_vmr + + ! Input/output arguments + real(kind_phys), intent(inout) :: const_array(:,:,:) ! Constituents array + + ! Output arguments + character(len=512), intent(out) :: errmsg + integer, intent(out) :: errcode + + ! Local variables + integer :: const_idx !Constituents object index + real(kind_phys) :: dry_air_to_const_molar_mass_ratio !Ratio of dry air molar mass to constituent molar mass + + !+++++++++++++++++++++++++++++++++++ + ! Convert number/mole fraction into + ! mass mixing ratio w.r.t dry air + !+++++++++++++++++++++++++++++++++++ + + !---- + ! CH4: + !---- + + ! Check if CH4 is present in constituents object: + call ccpp_constituent_index('CH4', const_idx, errcode, errmsg) + if (errcode /= 0) then + return + else if (const_idx /= int_unassigned) then + + ! Get ratio of molar mass of dry air / constituent molar mass + call get_molar_mass_ratio('CH4', dry_air_to_const_molar_mass_ratio, errmsg, errcode) + if (errcode /= 0) then + return + end if + + ! Convert namelist-provided number/mole fraction to + ! mass mixing ratio w.r.t. dry air, and set constituents + ! array to new coonverted value: + const_array(:,:,const_idx) = ch4_vmr*dry_air_to_const_molar_mass_ratio + + end if + + !---- + ! CO2: + !---- + + ! Check if CO2 is present in constituents object: + call ccpp_constituent_index('CO2', const_idx, errcode, errmsg) + if (errcode /= 0) then + return + else if (const_idx /= int_unassigned) then + + ! Get ratio of molar mass of dry air / constituent molar mass + call get_molar_mass_ratio('CO2', dry_air_to_const_molar_mass_ratio, errmsg, errcode) + if (errcode /= 0) then + return + end if + + ! Convert namelist-provided number/mole fraction to + ! mass mixing ratio w.r.t. dry air, and set constituents + ! array to new coonverted value: + const_array(:,:,const_idx) = co2_vmr*dry_air_to_const_molar_mass_ratio + + end if + + !------ + ! CFC11: + !------ + + ! Check if CFC-11 is present in constituents object: + call ccpp_constituent_index('CFC11', const_idx, errcode, errmsg) + if (errcode /= 0) then + return + else if (const_idx /= int_unassigned) then + + ! Get ratio of molar mass of dry air / constituent molar mass + call get_molar_mass_ratio('CFC11', dry_air_to_const_molar_mass_ratio, errmsg, errcode) + if (errcode /= 0) then + return + end if + + ! Convert namelist-provided number/mole fraction to + ! mass mixing ratio w.r.t. dry air, and set constituents + ! array to new coonverted value: + const_array(:,:,const_idx) = cfc11_vmr*dry_air_to_const_molar_mass_ratio + + end if + + !------ + ! CFC12: + !------ + + ! Check if CFC-12 is present in constituents object: + call ccpp_constituent_index('CFC12', const_idx, errcode, errmsg) + if (errcode /= 0) then + return + else if (const_idx /= int_unassigned) then + + ! Get ratio of molar mass of dry air / constituent molar mass + call get_molar_mass_ratio('CFC12', dry_air_to_const_molar_mass_ratio, errmsg, errcode) + if (errcode /= 0) then + return + end if + + ! Convert namelist-provided number/mole fraction to + ! mass mixing ratio w.r.t. dry air, and set constituents + ! array to new coonverted value: + const_array(:,:,const_idx) = cfc12_vmr*dry_air_to_const_molar_mass_ratio + + end if + + !---- + ! N2O: + !---- + + ! Check if N2O is present in constituents object: + call ccpp_constituent_index('N2O', const_idx, errcode, errmsg) + if (errcode /= 0) then + return + else if (const_idx /= int_unassigned) then + + ! Get ratio of molar mass of dry air / constituent molar mass + call get_molar_mass_ratio('N2O', dry_air_to_const_molar_mass_ratio, errmsg, errcode) + if (errcode /= 0) then + return + end if + + ! Convert namelist-provided number/mole fraction to + ! mass mixing ratio w.r.t. dry air, and set constituents + ! array to new coonverted value: + const_array(:,:,const_idx) = n2o_vmr*dry_air_to_const_molar_mass_ratio + + end if + + ! Set error variables + errmsg = '' + errcode = 0 + + end subroutine prescribe_radiative_gas_concentrations_init + +end module prescribe_radiative_gas_concentrations \ No newline at end of file diff --git a/schemes/radiation_utils/prescribe_radiative_gas_concentrations.meta b/schemes/radiation_utils/prescribe_radiative_gas_concentrations.meta new file mode 100644 index 00000000..55be49a4 --- /dev/null +++ b/schemes/radiation_utils/prescribe_radiative_gas_concentrations.meta @@ -0,0 +1,56 @@ +[ccpp-table-properties] + name = prescribe_radiative_gas_concentrations + type = scheme + dependencies = ../rrtmgp/utils/radiation_utils.F90 + +[ccpp-arg-table] + name = prescribe_radiative_gas_concentrations_init + type = scheme +[ ch4_vmr ] + standard_name = prescribed_volume_mixing_ratio_of_ch4 + units = mol mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ co2_vmr ] + standard_name = prescribed_volume_mixing_ratio_of_co2 + units = mol mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ cfc11_vmr ] + standard_name = prescribed_volume_mixing_ratio_of_cfc11 + units = mol mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ cfc12_vmr ] + standard_name = prescribed_volume_mixing_ratio_of_cfc12 + units = mol mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ n2o_vmr ] + standard_name = prescribed_volume_mixing_ratio_of_n2o + units = mol mol-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ const_array ] + standard_name = ccpp_constituents + units = none + type = real | kind = kind_phys + dimensions = (horizontal_dimension,vertical_layer_dimension,number_of_ccpp_constituents) + intent = inout +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errcode ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/radiation_utils/prescribe_radiative_gas_concentrations_namelist.xml b/schemes/radiation_utils/prescribe_radiative_gas_concentrations_namelist.xml new file mode 100644 index 00000000..7134e532 --- /dev/null +++ b/schemes/radiation_utils/prescribe_radiative_gas_concentrations_namelist.xml @@ -0,0 +1,157 @@ + + + + + + + + + real + kind_phys + ghg_cam + chem_surfvals + prescribed_volume_mixing_ratio_of_ch4 + mol mol-1 + + CH4 volume mixing ratio. This is used as the time invariant value + of CH4 if no time varying values are specified. + + + 1760.0e-9 + 1650.0e-9 + + + + real + kind_phys + ghg_cam + chem_surfvals + prescribed_volume_mixing_ratio_of_co2 + mol mol-1 + + CO2 volume mixing ratio. This is used as the time invariant value + of CO2 if no time varying values are specified. + + + 367.0e-6 + 348.0e-6 + + + + real + kind_phys + ghg_cam + chem_surfvals + prescribed_volume_mixing_ratio_of_cfc11 + mol mol-1 + + CFC-11 volume mixing ratio. This is used as the time invariant value + of CFC-11 if no time varying values are specified. + + + 653.45e-12 + + + + real + kind_phys + ghg_cam + chem_surfvals + prescribed_volume_mixing_ratio_of_cfc12 + mol mol-1 + + CFC-12 volume mixing ratio. This is used as the time invariant value + of CFC-12 if no time varying values are specified. + + + 535.0e-12 + + + + real + kind_phys + ghg_cam + chem_surfvals + prescribed_volume_mixing_ratio_of_n2o + mol mol-1 + + N2O volume mixing ratio. This is used as the time invariant value + of N2O if no time varying values are specified. + + + 316.0e-9 + 306.0e-9 + + + + diff --git a/schemes/radiation_utils/solar_irradiance_data_namelist.xml b/schemes/radiation_utils/solar_irradiance_data_namelist.xml index f1bfce2c..680eaeef 100644 --- a/schemes/radiation_utils/solar_irradiance_data_namelist.xml +++ b/schemes/radiation_utils/solar_irradiance_data_namelist.xml @@ -82,10 +82,13 @@ filename_of_solar_irradiance_data none - The filename of the solar irradiance data. + The filename of the solar irradiance data. If the file + path is set to "NONE" then it is assumed that a + constant solar irradiance is being used. ${DIN_LOC_ROOT}/atm/cam/solar/SolarForcingCMIP7piControl_c20250103.nc + NONE @@ -112,6 +115,7 @@ .true. + .false. @@ -152,6 +156,7 @@ -9999D0 + 1365.0 diff --git a/schemes/rasch_kristjansson/prognostic_cloud_water.F90 b/schemes/rasch_kristjansson/prognostic_cloud_water.F90 index d2f4d117..b239b42d 100644 --- a/schemes/rasch_kristjansson/prognostic_cloud_water.F90 +++ b/schemes/rasch_kristjansson/prognostic_cloud_water.F90 @@ -513,9 +513,13 @@ subroutine prognostic_cloud_water_run( & ! The cloud microphysics is highly nonlinear and coupled with qme ! Both rain processes and qme are calculated iteratively. qme_iter_loop: do l = 1, iter + + ! Adjust relative humidity above the tropopause in higher + ! latitude regions: + call relhum_min_adj(ncol, pver, tropLev, dlat, rhu00, rhu_adj) + qme_update_loop: do i = 1, ncol ! calculation of qme has 4 scenarios - call relhum_min_adj(ncol, pver, tropLev, dlat, rhu00, rhu_adj) if(relhum(i) > rhu_adj(i,k)) then ! 1. whole grid saturation diff --git a/schemes/rrtmgp/rrtmgp_constituents_namelist.xml b/schemes/rrtmgp/rrtmgp_constituents_namelist.xml index 670c2c82..8720ed53 100644 --- a/schemes/rrtmgp/rrtmgp_constituents_namelist.xml +++ b/schemes/rrtmgp/rrtmgp_constituents_namelist.xml @@ -86,7 +86,24 @@ calculation in RRTMGP (H2O not included because it is assumed to be advected if present). - 'N:O2:O2', 'A:CO2:CO2', 'N:ozone:O3', 'A:N2O:N2O', 'A:CH4:CH4', 'N:CFC11STAR:CFC11', 'A:CFC12:CFC12' + + 'N:O2:O2', + 'A:CO2:CO2', + 'N:ozone:O3', + 'A:N2O:N2O', + 'A:CH4:CH4', + 'N:CFC11STAR:CFC11', + 'A:CFC12:CFC12' + + + 'N:O2:O2', + 'N:CO2:CO2', + 'N:ozone:O3', + 'N:N2O:N2O', + 'N:CH4:CH4', + 'N:CFC11:CFC11', + 'N:CFC12:CFC12' + diff --git a/schemes/utilities/set_surface_coupling_vars.F90 b/schemes/utilities/set_surface_coupling_vars.F90 new file mode 100644 index 00000000..6e2eed6f --- /dev/null +++ b/schemes/utilities/set_surface_coupling_vars.F90 @@ -0,0 +1,201 @@ +!----------------------------------------------------------------------- +! Module to handle data that is exchanged between the atmosphere +! model and the surface models (land, sea-ice, ocean, etc.). + +! Please note that currently this is a SIMA-specific module, +! but could be made more host model-independent in the future. +!----------------------------------------------------------------------- +module set_surface_coupling_vars + + implicit none + private + + ! Public interfaces + public :: set_surface_coupling_vars_run + +!=============================================================================== +CONTAINS +!=============================================================================== + +!> \section arg_table_set_surface_coupling_vars_run Argument Table +!! \htmlinclude set_surface_coupling_vars_run.html +!! +subroutine set_surface_coupling_vars_run(ncol, pver, ncnst, gravit, rair, phis, & + surf_pres, air_temp, inv_ps_exner, zm, & + rho, uwnd, vwnd, air_pres, cnst_array, & + prec_dp, snow_dp, prec_sh, snow_sh, & + prec_str, snow_str, air_temp_bot, & + pot_temp_bot, zm_bot, rho_bot, & + uwnd_bot, vwnd_bot, air_pres_bot, & + topo_height, sea_lev_pres, cnst_bot, & + conv_prec, conv_snow, strat_prec, & + strat_snow, errmsg, errcode) + + ! Set surface variables needed for atmosphere-surface coupling. + + ! Use statements + use ccpp_kinds, only: kind_phys + + ! Input arguments + integer, intent(in) :: ncol + integer, intent(in) :: pver + integer, intent(in) :: ncnst + + real(kind_phys), intent(in) :: gravit ! Gravitational acceleration [m s-2] + real(kind_phys), intent(in) :: rair ! Dry air gas constant [J kg-1 K-1] + real(kind_phys), intent(in) :: phis(:) ! Surface geopotential [m2 s-2] + real(kind_phys), intent(in) :: surf_pres(:) ! Surface pressure [Pa] + real(kind_phys), intent(in) :: air_temp(:,:) ! Air temperature [K] + real(kind_phys), intent(in) :: inv_ps_exner(:,:) ! Inverse Exner function w.r.t. surface pressure [1] + real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height wrt surface [m] + real(kind_phys), intent(in) :: rho(:,:) ! Dry air density [kg m-3] + real(kind_phys), intent(in) :: uwnd(:,:) ! Eastward wind [m s-1] + real(kind_phys), intent(in) :: vwnd(:,:) ! Northward wind [m s-1] + real(kind_phys), intent(in) :: air_pres(:,:) ! Air pressure [Pa] + real(kind_phys), intent(in) :: cnst_array(:,:,:) ! CCPP constituents array [none] + + real(kind_phys), intent(in) :: prec_dp(:) ! LWE deep convective precipitation (all phases) [m s-1] + real(kind_phys), intent(in) :: snow_dp(:) ! LWE deep convective frozen precipitation (e.g. snow) [m s-1] + real(kind_phys), intent(in) :: prec_sh(:) ! LWE shallow convective precipitation (all phases) [m s-1] + real(kind_phys), intent(in) :: snow_sh(:) ! LWE shallow convective frozen precipitation (all phases) [m s-1] + real(kind_phys), intent(in) :: prec_str(:) ! LWE stratiform precipitation (all phases) [m s-1] + real(kind_phys), intent(in) :: snow_str(:) ! LWE stratiform frozen precipitation (e.g. snow) [m s-1] + + ! Output arguments + real(kind_phys), intent(out) :: air_temp_bot(:) ! Air temperature at bottom of atmosphere for coupling [K] + real(kind_phys), intent(out) :: pot_temp_bot(:) ! Potential air temprature at bottom of atmosphere for coupling [K] + real(kind_phys), intent(out) :: zm_bot(:) ! Geopotential height wrt surface at bottom of atmosphere for coupling [m] + real(kind_phys), intent(out) :: rho_bot(:) ! Dry air density at bottom of atmosphere for coupling [kg m-3] + real(kind_phys), intent(out) :: uwnd_bot(:) ! Eastward wind at bottom of atmosphere for coupling [m s-1] + real(kind_phys), intent(out) :: vwnd_bot(:) ! Northward wind at bottom of atmosphere for coupling [m s-1] + real(kind_phys), intent(out) :: air_pres_bot(:) ! Air pressure at bottom of atmosphere for coupling [Pa] + real(kind_phys), intent(out) :: topo_height(:) ! Geopotential height of surface topography [m] + real(kind_phys), intent(out) :: sea_lev_pres(:) ! Sea Level Pressure [Pa] + real(kind_phys), intent(out) :: cnst_bot(:,:) ! Constituent mixing ratios at bottom of atmosphere for coupling [kg kg-1] + + real(kind_phys), intent(out) :: conv_prec(:) ! LWE Convective precipitation (all phases) [m s-1] + real(kind_phys), intent(out) :: conv_snow(:) ! LWE Frozen convective precipitation (e.g. snow) [m s-1] + real(kind_phys), intent(out) :: strat_prec(:) ! LWE Stratiform (large-scale) precipitation (all phases) [m s-1] + real(kind_phys), intent(out) :: strat_snow(:) ! LWE Frozen stratiform (large-scale)precipitation (e.g. snow) [m s-1] + + character(len=512), intent(out) :: errmsg ! CCPP error message + integer, intent(out) :: errcode ! CCPP error code + + ! Local variables + + integer :: i ! column index + integer :: m ! constituent index + + !----------------------------------------------------------------------- + + errmsg = '' + errcode = 0 + + ! Copy over atmospheric 3-D variables to surface fields: + do i=1,ncol + air_temp_bot(i) = air_temp(i,pver) + pot_temp_bot(i) = air_temp(i,pver) * inv_ps_exner(i,pver) + zm_bot(i) = zm(i,pver) + rho_bot(i) = rho(i,pver) + uwnd_bot(i) = uwnd(i,pver) + vwnd_bot(i) = vwnd(i,pver) + air_pres_bot(i) = air_pres(i,pver) + topo_height(i) = phis(i)/gravit + end do + + ! Calculate Sea Level Pressure (PSL): + call calc_sea_level_pressure(ncol, gravit, rair, air_pres_bot, phis, & + surf_pres, air_temp_bot, sea_lev_pres) + + ! Set surface constituent values: + do m = 1, ncnst + do i = 1, ncol + cnst_bot(i,m) = cnst_array(i,pver,m) + end do + end do + + ! + ! Precipation and snow rates from shallow convection, deep convection and stratiform processes. + ! Compute total convective and stratiform precipitation and snow rates + ! + do i=1,ncol + conv_prec(i) = prec_dp(i) + prec_sh(i) + conv_snow(i) = snow_dp(i) + snow_sh(i) + + strat_prec(i) = prec_str(i) !Might be better to have microphysics set these directly + strat_snow(i) = snow_str(i) + end do + +end subroutine set_surface_coupling_vars_run + +!*************************************************** +! private subroutine to calculate sea level pressure +!*************************************************** + +subroutine calc_sea_level_pressure(ncol, gravit, rair, pbot, phis, ps, tbot, psl) + +!----------------------------------------------------------------------- +! +! Compute sea level pressure. +! +! Uses ECMWF formulation Algorithm: See section 3.1.b in NCAR NT-396 "Vertical +! Interpolation and Truncation of Model-Coordinate Data +! +!----------------------------------------------------------------------- + + ! Use statements + use ccpp_kinds, only: kind_phys + + !-----------------------------Arguments--------------------------------- + integer , intent(in) :: ncol ! longitude dimension + + real(kind_phys), intent(in) :: gravit ! Gravitational acceleration [m s-2] + real(kind_phys), intent(in) :: rair ! gas constant for dry air [J kg-1 K-1] + real(kind_phys), intent(in) :: pbot(:) ! Bottom layer atmospheric pressure [Pa] + real(kind_phys), intent(in) :: phis(:) ! Surface geopotential [m2 s-2] + real(kind_phys), intent(in) :: ps(:) ! Surface air pressure [Pa] + real(kind_phys), intent(in) :: Tbot(:) ! Bottom layer Atmospheric temperature [K] + + real(kind_phys), intent(out):: psl(:) ! Sea level pressures [Pa] + + !-----------------------------Parameters-------------------------------- + real(kind_phys), parameter :: xlapse = 6.5e-3_kind_phys ! Temperature lapse rate [K m-1] + + !-----------------------------Local Variables--------------------------- + integer :: i ! Loop index + real(kind_phys) :: alpha ! Temperature lapse rate in terms of pressure ratio (unitless) + real(kind_phys) :: Tstar ! Computed surface temperature + real(kind_phys) :: TT0 ! Computed temperature at sea-level + real(kind_phys) :: alph ! Power to raise P/Ps to get rate of increase of T with pressure + real(kind_phys) :: beta ! alpha*phis/(R*T) term used in approximation of PSL + !----------------------------------------------------------------------- + + alpha = rair*xlapse/gravit + do i=1,ncol + if ( abs(phis(i)/gravit) < 1.e-4_kind_phys )then + psl(i)=ps(i) + else + Tstar=Tbot(i)*(1._kind_phys+alpha*(ps(i)/pbot(i)-1._kind_phys)) ! pg 7 eq 5 + + TT0=Tstar + xlapse*phis(i)/gravit ! pg 8 eq 13 + + if ( Tstar<=290.5_kind_phys .and. TT0>290.5_kind_phys ) then ! pg 8 eq 14.1 + alph=rair/phis(i)*(290.5_kind_phys-Tstar) + else if (Tstar>290.5_kind_phys .and. TT0>290.5_kind_phys) then ! pg 8 eq 14.2 + alph=0._kind_phys + Tstar= 0.5_kind_phys * (290.5_kind_phys + Tstar) + else + alph=alpha + if (Tstar<255._kind_phys) then + Tstar= 0.5_kind_phys * (255._kind_phys + Tstar) ! pg 8 eq 14.3 + end if + end if + + beta = phis(i)/(rair*Tstar) + psl(i)=ps(i)*exp( beta*(1._kind_phys-alph*beta/2._kind_phys+((alph*beta)**2)/3._kind_phys)) + end if + end do + +end subroutine calc_sea_level_pressure + +end module set_surface_coupling_vars diff --git a/schemes/utilities/set_surface_coupling_vars.meta b/schemes/utilities/set_surface_coupling_vars.meta new file mode 100644 index 00000000..73f32544 --- /dev/null +++ b/schemes/utilities/set_surface_coupling_vars.meta @@ -0,0 +1,240 @@ +[ccpp-table-properties] + name = set_surface_coupling_vars + type = scheme + +[ccpp-arg-table] + name = set_surface_coupling_vars_run + type = scheme + +#---------------- +#Input Arguments +#---------------- + +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ pver ] + standard_name = vertical_layer_dimension + units = count + type = integer + dimensions = () + intent = in +[ ncnst ] + standard_name = number_of_ccpp_constituents + units = count + type = integer + dimensions = () + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rair ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ phis ] + standard_name = surface_geopotential + type = real | kind = kind_phys + units = m2 s-2 + dimensions = (horizontal_loop_extent) + intent = in +[ surf_pres ] + standard_name = surface_air_pressure + type = real | kind = kind_phys + units = Pa + dimensions = (horizontal_loop_extent) + intent = in +[ air_temp ] + standard_name = air_temperature + type = real | kind = kind_phys + units = K + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ inv_ps_exner ] + standard_name = reciprocal_of_dimensionless_exner_function_wrt_surface_air_pressure + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ zm ] + standard_name = geopotential_height_wrt_surface + type = real | kind = kind_phys + units = m + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rho ] + standard_name = dry_air_density + units = kg m-3 + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + type = real | kind = kind_phys + intent = in +[ uwnd ] + standard_name = eastward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ vwnd ] + standard_name = northward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ air_pres ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ cnst_array ] + standard_name = ccpp_constituents + units = none + dimensions = (horizontal_loop_extent, vertical_layer_dimension, number_of_ccpp_constituents) + type = real | kind = kind_phys + intent = in +[ prec_dp ] + standard_name = lwe_precipitation_rate_at_surface_due_to_deep_convection + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snow_dp ] + standard_name = lwe_frozen_precipitation_rate_at_surface_due_to_deep_convection + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ prec_sh ] + standard_name = lwe_precipitation_rate_at_surface_due_to_shallow_convection + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snow_sh ] + standard_name = lwe_frozen_precipitation_rate_at_surface_due_to_shallow_convection + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ prec_str ] + standard_name = lwe_large_scale_precipitation_rate_at_surface + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snow_str ] + standard_name = lwe_snow_and_cloud_ice_precipitation_rate_at_surface_due_to_microphysics + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in + +#---------------- +#Output Arguments +#---------------- + +[ air_temp_bot ] + standard_name = air_temperature_at_bottom_layer_to_coupler + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ pot_temp_bot ] + standard_name = potential_temperature_wrt_surface_pressure_at_bottom_layer_to_coupler + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ zm_bot ] + standard_name = geopotential_height_wrt_surface_at_bottom_layer_to_coupler + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ rho_bot ] + standard_name = dry_air_density_at_bottom_layer_to_coupler + units = kg m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ uwnd_bot ] + standard_name = eastward_wind_at_bottom_layer_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ vwnd_bot ] + standard_name = northward_wind_at_bottom_layer_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ air_pres_bot ] + standard_name = air_pressure_at_bottom_layer_to_coupler + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ topo_height ] + standard_name = surface_topographic_height_to_coupler + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ sea_lev_pres ] + standard_name = air_pressure_at_sea_level_to_coupler + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cnst_bot ] + standard_name = constituent_mixing_ratio_at_bottom_layer_to_coupler + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, number_of_ccpp_constituents) + intent = out +[ conv_prec ] + standard_name = lwe_convective_precipitation_rate_at_surface_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ conv_snow ] + standard_name = lwe_convective_snowfall_rate_at_surface_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ strat_prec ] + standard_name = lwe_large_scale_precipitation_rate_at_surface_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ strat_snow ] + standard_name = lwe_large_scale_snowfall_rate_at_surface_to_coupler + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + units = none + type = character | kind = len=512 + dimensions = () + intent = out +[ errcode ] + standard_name = ccpp_error_code + units = 1 + type = integer + dimensions = () + intent = out + diff --git a/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta index d091f95c..9ee426cf 100644 --- a/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta +++ b/schemes/vertical_diffusion/vertical_diffusion_sponge_layer.meta @@ -23,6 +23,12 @@ dimensions = () type = real | kind = kind_phys intent = in +[ diff_sponge_fac ] + standard_name = scale_factor_for_sponge_layer_vertical_diffusion + units = 1 + dimensions = () + type = real | kind = kind_phys + intent = in [ errmsg ] standard_name = ccpp_error_message units = none @@ -51,6 +57,12 @@ dimensions = () type = integer intent = in +[ diff_sponge_fac ] + standard_name = scale_factor_for_sponge_layer_vertical_diffusion + units = 1 + dimensions = () + type = real | kind = kind_phys + intent = in [ kvm ] standard_name = eddy_momentum_diffusivity_at_interfaces units = m2 s-1 diff --git a/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml b/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml index 001fde77..2ed47487 100644 --- a/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml +++ b/schemes/zhang_mcfarlane/zm_conv_options_namelist.xml @@ -14,6 +14,7 @@ Autoconversion coefficient over land in ZM deep convection scheme. 0.0075D0 + 0.0035D0 @@ -27,6 +28,7 @@ Autoconversion coefficient over ocean in ZM deep convection scheme. 0.0300D0 + 0.0035D0 @@ -66,6 +68,8 @@ Tunable evaporation efficiency for land in ZM deep convection scheme. 1.0E-5 + 3.0E-6 + 3.0E-6 diff --git a/suites/suite_cam4.xml b/suites/suite_cam4.xml index 3a3f41c5..de4f1484 100644 --- a/suites/suite_cam4.xml +++ b/suites/suite_cam4.xml @@ -7,15 +7,15 @@ Shallow convection Hack Macrophysics RK Microphysics RK - Radiation RRTMGP (commented out) + Radiation RRTMGP Chemistry None (not implemented) Vertical Diffusion HB Gravity Wave Drag Orographic --> - - initialize_constituents + to_be_ccppized_temporary + prescribe_radiative_gas_concentrations zm_conv_options @@ -215,16 +215,16 @@ sima_state_diagnostics - + rrtmgp_cloud_diagnostics - + rrtmgp_lw_mcica_subcol_gen - + - + rrtmgp_diagnostics - + apply_heating_rate + geopotential_temp @@ -269,6 +269,12 @@ tropopause_find tropopause_diagnostics + + + calc_dry_air_ideal_gas_density + set_surface_coupling_vars + @@ -359,6 +365,7 @@ + gravity_wave_drag_common @@ -394,6 +401,7 @@ + check_energy_save_teout diff --git a/suites/suite_cam7.xml b/suites/suite_cam7.xml index b39cebe4..e7508b7f 100644 --- a/suites/suite_cam7.xml +++ b/suites/suite_cam7.xml @@ -5,7 +5,6 @@ - initialize_constituents to_be_ccppized_temporary @@ -75,6 +74,13 @@ check_energy_scaling check_energy_chng + + + + + From d8e279ec5d2666e38140549a52b98a67258ca02c Mon Sep 17 00:00:00 2001 From: Kuan-Chih Wang Date: Tue, 10 Feb 2026 13:48:46 -0700 Subject: [PATCH 2/2] Add MYNN surface layer scheme (#340) ### Originator(s): kuanchihwang ### Descriptions (include the issue title, and the keyword ['closes', 'fixes', 'resolves'] followed by the issue number): This PR adds the Mellor-Yamada-Nakanishi-Niino (MYNN) surface layer scheme to the experimental convection-permitting physics suite. ### List all namelist files that were added or changed: None ### List all files eliminated and why: None ### List all files added and what they do: ``` A .gitattributes * Make GitHub Linguist recognize `*.pf` files as Fortran source code A schemes/mmm/sf_mynn_compat.F90 A schemes/mmm/sf_mynn_compat.meta * Implement MYNN surface layer scheme ``` ### List all existing files that have been modified, and describe the changes: ``` M .github/workflows/unit-tests.yaml * Do not hard code path to working directory in GitHub Actions M .gitmodules M schemes/mmm/mmm_physics * Update git submodule of MMM physics M schemes/mmm/CMakeLists.txt M test/unit-test/tests/mmm/CMakeLists.txt * Synchronize `CMakeLists.txt` of MMM physics with MPAS dycore in CAM-SIMA M schemes/mmm/cu_ntiedtke_compat.F90 M schemes/mmm/cu_ntiedtke_compat.meta * Fix incorrect standard name in new Tiedtke convection scheme M schemes/mmm/mmm_physics_compat.F90 M schemes/mmm/mmm_physics_compat.meta * Implement interstitial schemes for MYNN surface layer scheme M test/test_suites/suite_convection_permitting.xml * Add MYNN surface layer scheme to convection-permitting suite M test/unit-test/tests/mmm/mmm_physics_compat_tests.pf * Implement unit tests ``` ### List all automated tests that failed, as well as an explanation for why they were not fixed: None ### Is this an answer-changing PR? If so, is it a new physics package, algorithm change, tuning change, etc? Answer-changing for the convection-permitting physics suite due to a newly added physics scheme. Nothing is changed for the rest. ### If yes to the above question, describe how this code was validated with the new/modified features: The convection-permitting physics suite is considered an experimental feature. There is no baseline available to validate against because it has never been implemented in CAM-SIMA as well as CAM before. --- .gitattributes | 1 + .github/workflows/unit-tests.yaml | 2 +- .gitmodules | 2 +- schemes/mmm/CMakeLists.txt | 14 +- schemes/mmm/cu_ntiedtke_compat.F90 | 27 +- schemes/mmm/cu_ntiedtke_compat.meta | 20 +- schemes/mmm/mmm_physics | 2 +- schemes/mmm/mmm_physics_compat.F90 | 49 +- schemes/mmm/mmm_physics_compat.meta | 72 ++ schemes/mmm/sf_mynn_compat.F90 | 541 ++++++++++ schemes/mmm/sf_mynn_compat.meta | 982 ++++++++++++++++++ schemes/utilities/state_converters.F90 | 57 +- schemes/utilities/state_converters.meta | 88 ++ .../suite_convection_permitting.xml | 11 + test/unit-test/tests/mmm/CMakeLists.txt | 6 + .../tests/mmm/mmm_physics_compat_tests.pf | 67 +- .../tests/utilities/test_state_converters.pf | 60 ++ 17 files changed, 1960 insertions(+), 41 deletions(-) create mode 100644 .gitattributes create mode 100644 schemes/mmm/sf_mynn_compat.F90 create mode 100644 schemes/mmm/sf_mynn_compat.meta diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 00000000..846e6dbf --- /dev/null +++ b/.gitattributes @@ -0,0 +1 @@ +*.pf linguist-language=Fortran-Free-Form diff --git a/.github/workflows/unit-tests.yaml b/.github/workflows/unit-tests.yaml index 45ef3a57..bb54666f 100644 --- a/.github/workflows/unit-tests.yaml +++ b/.github/workflows/unit-tests.yaml @@ -35,7 +35,7 @@ jobs: - name: Build atmospheric_physics run: | cmake \ - -DCMAKE_PREFIX_PATH=/home/runner/work/atmospheric_physics/atmospheric_physics/pFUnit/build/installed \ + -DCMAKE_PREFIX_PATH=$GITHUB_WORKSPACE/pFUnit/build/installed \ -DATMOSPHERIC_PHYSICS_ENABLE_CODE_COVERAGE=ON \ -B./build \ -S./test/unit-test diff --git a/.gitmodules b/.gitmodules index 6606b271..10431cb9 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,7 +1,7 @@ [submodule "mmm-physics"] path = schemes/mmm/mmm_physics url = https://github.com/NCAR/MMM-physics.git - fxtag = 20250616-MPASv8.3 + fxtag = 20251208-CAM-SIMA fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/NCAR/MMM-physics.git [submodule "pumas"] diff --git a/schemes/mmm/CMakeLists.txt b/schemes/mmm/CMakeLists.txt index e52c1192..82f278de 100644 --- a/schemes/mmm/CMakeLists.txt +++ b/schemes/mmm/CMakeLists.txt @@ -1,8 +1,8 @@ -cmake_minimum_required(VERSION 3.20) +cmake_minimum_required(VERSION 3.21) # `mmm_physics_compat` has not been integrated into the CMake build of any top level projects yet, # and this CMakeLists.txt file is currently for unit testing purposes only. -# Making a change to this CMakeLists.txt file will not impact the build of a parent project at this time. +# Therefore, making a change to this CMakeLists.txt file will not impact the build of a parent project at this time. project(mmm_physics_compat VERSION 0.1.0 @@ -19,11 +19,13 @@ target_sources(mmm_physics_compat ccpp_kind_types.F90 mmm_physics_compat.F90 ) -target_compile_options(mmm_physics_compat - PRIVATE - $<$,$>:-fbacktrace -fcheck=all -std=f2018 -Wall -Wextra -Wpedantic> -) target_include_directories(mmm_physics_compat INTERFACE ${CMAKE_CURRENT_BINARY_DIR} ) +target_compile_options(mmm_physics_compat + PRIVATE + $<$,$>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic> + $<$,$>:-check all -fp-model=precise -stand f18 -traceback -warn all> + $<$,$>:-check all -fp-model=precise -stand f18 -traceback -warn all> +) diff --git a/schemes/mmm/cu_ntiedtke_compat.F90 b/schemes/mmm/cu_ntiedtke_compat.F90 index 19717226..e98a8e25 100644 --- a/schemes/mmm/cu_ntiedtke_compat.F90 +++ b/schemes/mmm/cu_ntiedtke_compat.F90 @@ -14,21 +14,19 @@ module cu_ntiedtke_compat !! \htmlinclude cu_ntiedtke_compat_pre_run.html subroutine cu_ntiedtke_compat_pre_run( & cflx, exner, landfrac, & - lhf, shf, & rthdynten, rthblten, rthratenlw, rthratensw, & rqvdynten, rqvblten, & lndj, & - ptf, pqvf, hfx, evap, & + ptf, pqvf, evap, & errmsg, errflg) use ccpp_kinds, only: kind_phys use ccpp_scheme_utils, only: ccpp_constituent_index real(kind_phys), intent(in) :: cflx(:, :), exner(:, :), landfrac(:), & - lhf(:), shf(:), & rthdynten(:, :), rthblten(:, :), rthratenlw(:, :), rthratensw(:, :), & rqvdynten(:, :), rqvblten(:, :) integer, intent(out) :: lndj(:) - real(kind_phys), intent(out) :: ptf(:, :), pqvf(:, :), hfx(:), evap(:) + real(kind_phys), intent(out) :: ptf(:, :), pqvf(:, :), evap(:) character(*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -42,7 +40,6 @@ subroutine cu_ntiedtke_compat_pre_run( & ptf(:, :) = (rthdynten(:, :) + rthblten(:, :) + rthratenlw(:, :) + rthratensw(:, :)) * exner(:, :) pqvf(:, :) = rqvdynten(:, :) + rqvblten(:, :) - hfx(:) = lhf(:) + shf(:) evap(:) = 0.0_kind_phys call ccpp_constituent_index( & @@ -50,7 +47,7 @@ subroutine cu_ntiedtke_compat_pre_run( & if (errflg /= 0 .or. & water_vapor_mixing_ratio_index < lbound(cflx, 2) .or. water_vapor_mixing_ratio_index > ubound(cflx, 2)) then - errmsg = 'Failed to find desired constituent flux from cflx' + errmsg = 'cu_ntiedtke_compat_pre_run: Failed to find desired constituent flux from cflx' errflg = 1 return @@ -126,36 +123,54 @@ subroutine cu_ntiedtke_compat_run( & allocate(pu_local, source=pu, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pu_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if allocate(pv_local, source=pv, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pv_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if allocate(pt_local, source=pt, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pt_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if allocate(pqv_local, source=pqv, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqv_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if allocate(pqc_local, source=pqc, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqc_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if allocate(pqi_local, source=pqi, errmsg=errmsg, stat=errflg) if (errflg /= 0) then + errmsg = 'cu_ntiedtke_compat_run: Failed to allocate "pqi_local"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + return end if diff --git a/schemes/mmm/cu_ntiedtke_compat.meta b/schemes/mmm/cu_ntiedtke_compat.meta index 6a09d240..5e7d5ed1 100644 --- a/schemes/mmm/cu_ntiedtke_compat.meta +++ b/schemes/mmm/cu_ntiedtke_compat.meta @@ -23,18 +23,6 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent) intent = in -[ lhf ] - standard_name = surface_upward_latent_heat_flux_from_coupler - units = W m-2 - type = real | kind = kind_phys - dimensions = (horizontal_loop_extent) - intent = in -[ shf ] - standard_name = surface_upward_sensible_heat_flux_from_coupler - units = W m-2 - type = real | kind = kind_phys - dimensions = (horizontal_loop_extent) - intent = in [ rthdynten ] standard_name = tendency_of_air_potential_temperature_due_to_dynamics units = K s-1 @@ -89,12 +77,6 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent, vertical_layer_dimension) intent = out -[ hfx ] - standard_name = surface_upward_heat_flux - units = W m-2 - type = real | kind = kind_phys - dimensions = (horizontal_loop_extent) - intent = out [ evap ] standard_name = surface_upward_water_vapor_flux units = kg m-2 s-1 @@ -271,7 +253,7 @@ dimensions = (horizontal_loop_extent) intent = in [ hfx ] - standard_name = surface_upward_heat_flux + standard_name = surface_upward_sensible_heat_flux_from_coupler units = W m-2 type = real | kind = kind_phys dimensions = (horizontal_loop_extent) diff --git a/schemes/mmm/mmm_physics b/schemes/mmm/mmm_physics index a4baf7f3..6cb29bc8 160000 --- a/schemes/mmm/mmm_physics +++ b/schemes/mmm/mmm_physics @@ -1 +1 @@ -Subproject commit a4baf7f3243d1db0dbc5f63473f895bdbdc05c30 +Subproject commit 6cb29bc8f17bc3efdfaaa8a268ca53e2504aadbd diff --git a/schemes/mmm/mmm_physics_compat.F90 b/schemes/mmm/mmm_physics_compat.F90 index 1bd0dfcc..8ab26573 100644 --- a/schemes/mmm/mmm_physics_compat.F90 +++ b/schemes/mmm/mmm_physics_compat.F90 @@ -3,6 +3,7 @@ module mmm_physics_compat implicit none private + public :: mmm_physics_compat_init public :: mmm_physics_compat_run public :: mmm_physics_accumulate_tendencies_timestep_init public :: mmm_physics_accumulate_tendencies_run @@ -12,22 +13,53 @@ module mmm_physics_compat public :: geopotential_height_wrt_sfc_at_if_to_msl_run public :: geopotential_height_wrt_sfc_to_msl_run contains + !> \section arg_table_mmm_physics_compat_init Argument Table + !! \htmlinclude mmm_physics_compat_init.html + pure subroutine mmm_physics_compat_init( & + isfflx, isftcflx, iz0tlnd, & + spp_pbl, & + xice_threshold, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + + integer, intent(out) :: isfflx, isftcflx, iz0tlnd + logical, intent(out) :: spp_pbl + real(kind_phys), intent(out) :: xice_threshold + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + errmsg = '' + errflg = 0 + + ! These options are hardcoded to the same values as in the MPAS physics driver. + ! There are other possible values for them in WRF, but the following combination is the only supported one in MPAS. + isfflx = 1 + isftcflx = 0 + iz0tlnd = 0 + spp_pbl = .false. + xice_threshold = 0.02_kind_phys + end subroutine mmm_physics_compat_init + !> \section arg_table_mmm_physics_compat_run Argument Table !! \htmlinclude mmm_physics_compat_run.html pure subroutine mmm_physics_compat_run( & nstep, & dt, & theta_curr, theta_prev, qv_curr, qv_prev, & + icefrac, xice_threshold, landfrac, & scheme_name, & rthdynten, rqvdynten, & + xland, & errmsg, errflg) use ccpp_kinds, only: kind_phys integer, intent(in) :: nstep real(kind_phys), intent(in) :: dt, & - theta_curr(:, :), theta_prev(:, :), qv_curr(:, :), qv_prev(:, :) + theta_curr(:, :), theta_prev(:, :), qv_curr(:, :), qv_prev(:, :), & + icefrac(:), xice_threshold, landfrac(:) character(256), intent(out) :: scheme_name - real(kind_phys), intent(out) :: rthdynten(:, :), rqvdynten(:, :) + real(kind_phys), intent(out) :: rthdynten(:, :), rqvdynten(:, :), & + xland(:) character(*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -43,6 +75,16 @@ pure subroutine mmm_physics_compat_run( & rthdynten(:, :) = (theta_curr(:, :) - theta_prev(:, :)) / dt rqvdynten(:, :) = (qv_curr(:, :) - qv_prev(:, :)) / dt end if + + ! For MMM physics, land mask (`xland`) is defined as + ! * xland = 1.0 for land cells, including sea ice cells. + ! * xland = 2.0 for water cells. + where (landfrac >= 0.5_kind_phys .or. & + icefrac >= xice_threshold) + xland = 1.0_kind_phys + elsewhere + xland = 2.0_kind_phys + end where end subroutine mmm_physics_compat_run !> \section arg_table_mmm_physics_accumulate_tendencies_timestep_init Argument Table @@ -196,7 +238,8 @@ end subroutine mmm_physics_persist_states_timestep_final !> \section arg_table_compute_characteristic_grid_length_scale_init Argument Table !! \htmlinclude compute_characteristic_grid_length_scale_init.html pure subroutine compute_characteristic_grid_length_scale_init( & - omega, rearth, dx, & + omega, rearth, & + dx, & errmsg, errflg) use ccpp_kinds, only: kind_phys diff --git a/schemes/mmm/mmm_physics_compat.meta b/schemes/mmm/mmm_physics_compat.meta index 5cd2ae93..f6c74782 100644 --- a/schemes/mmm/mmm_physics_compat.meta +++ b/schemes/mmm/mmm_physics_compat.meta @@ -2,6 +2,54 @@ name = mmm_physics_compat type = scheme +[ccpp-arg-table] + name = mmm_physics_compat_init + type = scheme +[ isfflx ] + standard_name = control_for_surface_flux_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = out +[ isftcflx ] + standard_name = control_for_surface_roughness_length_over_water_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = out +[ iz0tlnd ] + standard_name = control_for_surface_roughness_length_over_land_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = out +[ spp_pbl ] + standard_name = flag_for_stochastically_perturbed_parameterization_for_mynn_scheme + units = flag + type = logical + dimensions = () + intent = out +[ xice_threshold ] + standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme + units = fraction + type = real | kind = kind_phys + dimensions = () + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + [ccpp-arg-table] name = mmm_physics_compat_run type = scheme @@ -41,6 +89,24 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent, vertical_layer_dimension) intent = in +[ icefrac ] + standard_name = sea_ice_area_fraction_from_coupler + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ xice_threshold ] + standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme + units = fraction + type = real | kind = kind_phys + dimensions = () + intent = in +[ landfrac ] + standard_name = land_area_fraction_from_coupler + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in [ scheme_name ] standard_name = scheme_name units = none @@ -59,6 +125,12 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent, vertical_layer_dimension) intent = out +[ xland ] + standard_name = land_binary_mask_for_mmm_scheme + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out [ errmsg ] standard_name = ccpp_error_message long_name = Error message for error handling in CCPP diff --git a/schemes/mmm/sf_mynn_compat.F90 b/schemes/mmm/sf_mynn_compat.F90 new file mode 100644 index 00000000..80440a42 --- /dev/null +++ b/schemes/mmm/sf_mynn_compat.F90 @@ -0,0 +1,541 @@ +!> This module contains interstitial schemes that are specific to MYNN surface layer scheme, +!> which is part of MMM physics. +module sf_mynn_compat + use ccpp_kinds, only: kind_phys + + implicit none + + private + public :: sf_mynn_compat_pre_run + public :: sf_mynn_compat_init + public :: sf_mynn_compat_run + public :: sf_mynn_diagnostics_init + public :: sf_mynn_diagnostics_run +contains + !> \section arg_table_sf_mynn_compat_pre_run Argument Table + !! \htmlinclude sf_mynn_compat_pre_run.html + subroutine sf_mynn_compat_pre_run( & + itimestep, & + spp_pbl, & + u, v, t, qv, p, dz, rho, & + icefrac, xice_threshold, landfrac, snowhice, snowhland, & + u1d, v1d, t1d, qv1d, p1d, dz8w1d, rho1d, & + u1d2, v1d2, dz2w1d, & + chs, chs2, cqs2, cpm, rmol, & + znt, ust, zol, mol, regime, psim, & + psih, qfx, & + flhc, flqc, snowh, qgh, qsfc, & + gz1oz0, wspd, br, svp1, svp2, & + svp3, svpt0, qcg, & + rstoch1d, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + + integer, intent(in) :: itimestep + logical, intent(in) :: spp_pbl + real(kind_phys), intent(in) :: u(:, :), v(:, :), t(:, :), qv(:, :), p(:, :), dz(:, :), rho(:, :), & + icefrac(:), xice_threshold, landfrac(:), snowhice(:), snowhland(:) + real(kind_phys), intent(out) :: u1d(:), v1d(:), t1d(:), qv1d(:), p1d(:), dz8w1d(:), rho1d(:), & + u1d2(:), v1d2(:), dz2w1d(:), & + chs(:), chs2(:), cqs2(:), cpm(:), rmol(:), & + znt(:), ust(:), zol(:), mol(:), regime(:), psim(:), & + psih(:), qfx(:), & + flhc(:), flqc(:), snowh(:), qgh(:), qsfc(:), & + gz1oz0(:), wspd(:), br(:), svp1, svp2, & + svp3, svpt0, qcg(:), & + rstoch1d(:) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Provide first guess at the first time step. + if (itimestep == 1) then + ust(:) = max(0.04_kind_phys * sqrt(u(:, 1) ** 2 + v(:, 1) ** 2), 0.001_kind_phys) + mol(:) = 0.0_kind_phys + qsfc(:) = qv(:, 1) ! Note that this `qv` is wet. + end if + + ! In CAM-SIMA, the first vertical index is at top of atmosphere. + ! The last one is at bottom of atmosphere, which is what we want here. + + u1d(:) = u(:, size(u, 2)) + v1d(:) = v(:, size(v, 2)) + t1d(:) = t(:, size(t, 2)) + qv1d(:) = qv(:, size(qv, 2)) + p1d(:) = p(:, size(p, 2)) + dz8w1d(:) = dz(:, size(dz, 2)) + rho1d(:) = rho(:, size(rho, 2)) + + u1d2(:) = u(:, size(u, 2) - 1) + v1d2(:) = v(:, size(v, 2) - 1) + dz2w1d(:) = dz(:, size(dz, 2) - 1) + + chs(:) = 0.0_kind_phys + chs2(:) = 0.0_kind_phys + cqs2(:) = 0.0_kind_phys + cpm(:) = 0.0_kind_phys + rmol(:) = 0.0_kind_phys + + znt(:) = 0.0_kind_phys + zol(:) = 0.0_kind_phys + regime(:) = 0.0_kind_phys + psim(:) = 0.0_kind_phys + + psih(:) = 0.0_kind_phys + qfx(:) = 0.0_kind_phys + + flhc(:) = 0.0_kind_phys + flqc(:) = 0.0_kind_phys + snowh(:) = 0.0_kind_phys + + where (landfrac >= 0.5_kind_phys) + snowh = snowhland + end where + + where (icefrac >= xice_threshold) + snowh = snowhice + end where + + qgh(:) = 0.0_kind_phys + + gz1oz0(:) = 0.0_kind_phys + wspd(:) = 0.0_kind_phys + br(:) = 0.0_kind_phys + + ! Constants in equation 10 from Bolton (1980). See + ! doi:10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2. + svp1 = 6.112_kind_phys + svp2 = 17.67_kind_phys + svp3 = 29.65_kind_phys + svpt0 = 273.15_kind_phys + + qcg(:) = 0.0_kind_phys ! Not used but still appear in the argument list... + + if (spp_pbl) then + errmsg = 'sf_mynn_compat_pre_run: Stochastically perturbed parameterization is not supported' + errflg = 1 + + return + else + rstoch1d(:) = 0.0_kind_phys + end if + + errmsg = '' + errflg = 0 + end subroutine sf_mynn_compat_pre_run + + !> \section arg_table_sf_mynn_compat_init Argument Table + !! \htmlinclude sf_mynn_compat_init.html + subroutine sf_mynn_compat_init( & + ust, mol, qsfc, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + use sf_mynn, only: sf_mynn_init + + real(kind_phys), intent(out) :: ust(:), mol(:), qsfc(:) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Precompute lookup tables in MYNN surface layer scheme. + call sf_mynn_init(errmsg, errflg) + + if (errflg /= 0) then + return + end if + + ! MYNN surface layer scheme takes time averages of these variables internally. + ! As a result, they must be able to persist across time steps. + ust(:) = 0.0_kind_phys + mol(:) = 0.0_kind_phys + qsfc(:) = 0.0_kind_phys + + errmsg = '' + errflg = 0 + end subroutine sf_mynn_compat_init + + !> \section arg_table_sf_mynn_compat_run Argument Table + !! \htmlinclude sf_mynn_compat_run.html + subroutine sf_mynn_compat_run( & + ncol, cflx, icefrac, xice_threshold, sst, & + u1d, v1d, t1d, qv1d, p1d, dz8w1d, rho1d, & + u1d2, v1d2, dz2w1d, cp, g, rovcp, r, xlv, & + psfcpa, chs, chs2, cqs2, cpm, pblh, rmol, & + znt, ust, mavail, zol, mol, regime, psim, & + psih, xland, hfx, qfx, tsk, u10, v10, th2, & + t2, q2, flhc, flqc, snowh, qgh, qsfc, lh, & + gz1oz0, wspd, br, isfflx, dx, svp1, svp2, & + svp3, svpt0, ep1, ep2, karman, qcg, & + itimestep, wstar, qstar, ustm, ck, cka, & + cd, cda, spp_pbl, rstoch1d, isftcflx, & + iz0tlnd, its, ite, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + use ccpp_scheme_utils, only: ccpp_constituent_index + use sf_mynn, only: sf_mynn_run + + ! Typical threshold between sea surface temperature and ice surface temperature. See + ! Algorithm Theoretical Basis Document (ATBD) for the MODIS Snow and Sea Ice-Mapping Algorithms, + ! Section 4.4.3 Ice Surface Temperature (IST) Algorithm. + real(kind_phys), parameter :: sea_ice_temperature_threshold = 271.4_kind_phys + + integer, intent(in) :: ncol, & + isfflx, & + itimestep, & + isftcflx, & + iz0tlnd, its, ite + logical, intent(in) :: spp_pbl + real(kind_phys), intent(in) :: icefrac(:), xice_threshold, sst(:), & + u1d(:), v1d(:), t1d(:), qv1d(:), p1d(:), dz8w1d(:), rho1d(:), & + u1d2(:), v1d2(:), dz2w1d(:), cp, g, rovcp, r, xlv, & + psfcpa(:), pblh(:), & + mavail(:), & + xland(:), tsk(:), & + snowh(:), & + dx(:), svp1, svp2, & + svp3, svpt0, ep1, ep2, karman, qcg(:), & + rstoch1d(:) + real(kind_phys), intent(inout) :: cflx(:, :), & + chs(:), chs2(:), cqs2(:), cpm(:), rmol(:), & + znt(:), ust(:), zol(:), mol(:), regime(:), psim(:), & + psih(:), hfx(:), qfx(:), & + flhc(:), flqc(:), qgh(:), qsfc(:), lh(:), & + gz1oz0(:), wspd(:), br(:) + real(kind_phys), intent(out) :: u10(:), v10(:), th2(:), & + t2(:), q2(:), & + wstar(:), qstar(:), ustm(:), ck(:), cka(:), & + cd(:), cda(:) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: water_vapor_mixing_ratio_index + real(kind_phys), allocatable :: ch(:) + + ! Special handling for sea ice cells like in MMM physics. + logical, allocatable :: mask_sea_ice_cell(:) + real(kind_phys), allocatable :: mavail_sea(:), & + xland_sea(:), tsk_sea(:) + real(kind_phys), allocatable :: chs_sea(:), chs2_sea(:), cqs2_sea(:), cpm_sea(:), rmol_sea(:), & + znt_sea(:), ust_sea(:), zol_sea(:), mol_sea(:), regime_sea(:), psim_sea(:), & + psih_sea(:), hfx_sea(:), qfx_sea(:), & + flhc_sea(:), flqc_sea(:), qgh_sea(:), qsfc_sea(:), lh_sea(:), & + gz1oz0_sea(:), wspd_sea(:), br_sea(:), & + ch_sea(:) + real(kind_phys), allocatable :: u10_sea(:), v10_sea(:), th2_sea(:), & + t2_sea(:), q2_sea(:), & + wstar_sea(:), qstar_sea(:), ustm_sea(:), ck_sea(:), cka_sea(:), & + cd_sea(:), cda_sea(:) + + ! `ch` is duplicate of `chs`, but for unknown reasons it is passed separately in the argument list + ! to MYNN surface layer scheme. + allocate(ch(ncol), errmsg=errmsg, stat=errflg) + + if (errflg /= 0) then + errmsg = 'sf_mynn_compat_run: Failed to allocate "ch"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + + return + end if + + ch(:) = chs(:) + + qfx(:) = 0.0_kind_phys + + call ccpp_constituent_index( & + 'water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water', water_vapor_mixing_ratio_index, errflg, errmsg) + + if (errflg /= 0 .or. & + water_vapor_mixing_ratio_index < lbound(cflx, 2) .or. water_vapor_mixing_ratio_index > ubound(cflx, 2)) then + errmsg = 'sf_mynn_compat_run: Failed to find desired constituent flux from cflx' + errflg = 1 + + return + end if + + qfx(:) = cflx(:, water_vapor_mixing_ratio_index) + + allocate(mask_sea_ice_cell(ncol), errmsg=errmsg, stat=errflg) + + if (errflg /= 0) then + errmsg = 'sf_mynn_compat_run: Failed to allocate "mask_sea_ice_cell"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + + return + end if + + mask_sea_ice_cell(:) = (icefrac >= xice_threshold) + + ! If there are sea ice cells, make local copies of the variables in preparation for the second pass. + if (any(mask_sea_ice_cell)) then + allocate(mavail_sea(ncol), xland_sea(ncol), tsk_sea(ncol), & + errmsg=errmsg, stat=errflg) + + if (errflg /= 0) then + errmsg = 'sf_mynn_compat_run: Failed to allocate "mavail_sea", "xland_sea", "tsk_sea"' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + + return + end if + + mavail_sea(:) = mavail(:) + xland_sea(:) = xland(:) + tsk_sea(:) = tsk(:) + + allocate(chs_sea(ncol), chs2_sea(ncol), cqs2_sea(ncol), cpm_sea(ncol), rmol_sea(ncol), & + znt_sea(ncol), ust_sea(ncol), zol_sea(ncol), mol_sea(ncol), regime_sea(ncol), psim_sea(ncol), & + psih_sea(ncol), hfx_sea(ncol), qfx_sea(ncol), & + flhc_sea(ncol), flqc_sea(ncol), qgh_sea(ncol), qsfc_sea(ncol), lh_sea(ncol), & + gz1oz0_sea(ncol), wspd_sea(ncol), br_sea(ncol), & + ch_sea(ncol), & + errmsg=errmsg, stat=errflg) + + if (errflg /= 0) then + errmsg = 'sf_mynn_compat_run: Failed to allocate "chs_sea", "chs2_sea", "cqs2_sea", ...' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + + return + end if + + chs_sea(:) = chs(:) + chs2_sea(:) = chs2(:) + cqs2_sea(:) = cqs2(:) + cpm_sea(:) = cpm(:) + rmol_sea(:) = rmol(:) + znt_sea(:) = znt(:) + ust_sea(:) = ust(:) + zol_sea(:) = zol(:) + mol_sea(:) = mol(:) + regime_sea(:) = regime(:) + psim_sea(:) = psim(:) + psih_sea(:) = psih(:) + hfx_sea(:) = hfx(:) + qfx_sea(:) = qfx(:) + flhc_sea(:) = flhc(:) + flqc_sea(:) = flqc(:) + qgh_sea(:) = qgh(:) + qsfc_sea(:) = qsfc(:) + lh_sea(:) = lh(:) + gz1oz0_sea(:) = gz1oz0(:) + wspd_sea(:) = wspd(:) + br_sea(:) = br(:) + ch_sea(:) = ch(:) + + allocate(u10_sea(ncol), v10_sea(ncol), th2_sea(ncol), & + t2_sea(ncol), q2_sea(ncol), & + wstar_sea(ncol), qstar_sea(ncol), ustm_sea(ncol), ck_sea(ncol), cka_sea(ncol), & + cd_sea(ncol), cda_sea(ncol), & + errmsg=errmsg, stat=errflg) + + if (errflg /= 0) then + errmsg = 'sf_mynn_compat_run: Failed to allocate "u10_sea", "v10_sea", "th2_sea", ...' // new_line('') // & + 'Allocation returned with error: ' // trim(adjustl(errmsg)) + + return + end if + + u10_sea(:) = u10(:) + v10_sea(:) = v10(:) + th2_sea(:) = th2(:) + t2_sea(:) = t2(:) + q2_sea(:) = q2(:) + wstar_sea(:) = wstar(:) + qstar_sea(:) = qstar(:) + ustm_sea(:) = ustm(:) + ck_sea(:) = ck(:) + cka_sea(:) = cka(:) + cd_sea(:) = cd(:) + cda_sea(:) = cda(:) + + where (mask_sea_ice_cell) + ! Set surface moisture availability to maximum. + mavail_sea = 1.0_kind_phys + + ! Impose minimum skin temperature. + tsk_sea = max(sea_ice_temperature_threshold, sst) + + ! Treat as water cells. + xland_sea = 2.0_kind_phys + + ! Set surface roughness length to calm water. + znt_sea = 0.0001_kind_phys + end where + end if + + ! First pass for all cells. + call sf_mynn_run( & + u1d, v1d, t1d, qv1d, p1d, dz8w1d, rho1d, & + u1d2, v1d2, dz2w1d, cp, g, rovcp, r, xlv, & + psfcpa, chs, chs2, cqs2, cpm, pblh, rmol, & + znt, ust, mavail, zol, mol, regime, psim, & + psih, xland, hfx, qfx, tsk, u10, v10, th2, & + t2, q2, flhc, flqc, snowh, qgh, qsfc, lh, & + gz1oz0, wspd, br, isfflx, dx, svp1, svp2, & + svp3, svpt0, ep1, ep2, karman, ch, qcg, & + itimestep, wstar, qstar, ustm, ck, cka, & + cd, cda, spp_pbl, rstoch1d, isftcflx, & + iz0tlnd, its, ite, & + errmsg, errflg) + + if (errflg /= 0) then + return + end if + + ! The first pass treated sea ice cells as land cells due to how `xland` is defined. + ! The second pass treats sea ice cells as water cells instead. + ! Then, for sea ice cells only, the final results are weighted averages according to their area fractions. + if (any(mask_sea_ice_cell)) then + call sf_mynn_run( & + u1d, v1d, t1d, qv1d, p1d, dz8w1d, rho1d, & + u1d2, v1d2, dz2w1d, cp, g, rovcp, r, xlv, & + psfcpa, chs_sea, chs2_sea, cqs2_sea, cpm_sea, pblh, rmol_sea, & + znt_sea, ust_sea, mavail_sea, zol_sea, mol_sea, regime_sea, psim_sea, & + psih_sea, xland_sea, hfx_sea, qfx_sea, tsk_sea, u10_sea, v10_sea, th2_sea, & + t2_sea, q2_sea, flhc_sea, flqc_sea, snowh, qgh_sea, qsfc_sea, lh_sea, & + gz1oz0_sea, wspd_sea, br_sea, isfflx, dx, svp1, svp2, & + svp3, svpt0, ep1, ep2, karman, ch_sea, qcg, & + itimestep, wstar_sea, qstar_sea, ustm_sea, ck_sea, cka_sea, & + cd_sea, cda_sea, spp_pbl, rstoch1d, isftcflx, & + iz0tlnd, its, ite, & + errmsg, errflg) + + if (errflg /= 0) then + return + end if + + ! Blend the final results between sea ice and sea water. + where (mask_sea_ice_cell) + chs = chs * icefrac + (1.0_kind_phys - icefrac) * chs_sea + chs2 = chs2 * icefrac + (1.0_kind_phys - icefrac) * chs2_sea + cqs2 = cqs2 * icefrac + (1.0_kind_phys - icefrac) * cqs2_sea + cpm = cpm * icefrac + (1.0_kind_phys - icefrac) * cpm_sea + rmol = rmol * icefrac + (1.0_kind_phys - icefrac) * rmol_sea + znt = znt * icefrac + (1.0_kind_phys - icefrac) * znt_sea + ust = ust * icefrac + (1.0_kind_phys - icefrac) * ust_sea + zol = zol * icefrac + (1.0_kind_phys - icefrac) * zol_sea + mol = mol * icefrac + (1.0_kind_phys - icefrac) * mol_sea + psim = psim * icefrac + (1.0_kind_phys - icefrac) * psim_sea + psih = psih * icefrac + (1.0_kind_phys - icefrac) * psih_sea + hfx = hfx * icefrac + (1.0_kind_phys - icefrac) * hfx_sea + qfx = qfx * icefrac + (1.0_kind_phys - icefrac) * qfx_sea + flhc = flhc * icefrac + (1.0_kind_phys - icefrac) * flhc_sea + flqc = flqc * icefrac + (1.0_kind_phys - icefrac) * flqc_sea + qgh = qgh * icefrac + (1.0_kind_phys - icefrac) * qgh_sea + qsfc = qsfc * icefrac + (1.0_kind_phys - icefrac) * qsfc_sea + lh = lh * icefrac + (1.0_kind_phys - icefrac) * lh_sea + gz1oz0 = gz1oz0 * icefrac + (1.0_kind_phys - icefrac) * gz1oz0_sea + wspd = wspd * icefrac + (1.0_kind_phys - icefrac) * wspd_sea + br = br * icefrac + (1.0_kind_phys - icefrac) * br_sea + + u10 = u10 * icefrac + (1.0_kind_phys - icefrac) * u10_sea + v10 = v10 * icefrac + (1.0_kind_phys - icefrac) * v10_sea + th2 = th2 * icefrac + (1.0_kind_phys - icefrac) * th2_sea + t2 = t2 * icefrac + (1.0_kind_phys - icefrac) * t2_sea + q2 = q2 * icefrac + (1.0_kind_phys - icefrac) * q2_sea + wstar = wstar * icefrac + (1.0_kind_phys - icefrac) * wstar_sea + qstar = qstar * icefrac + (1.0_kind_phys - icefrac) * qstar_sea + ustm = ustm * icefrac + (1.0_kind_phys - icefrac) * ustm_sea + ck = ck * icefrac + (1.0_kind_phys - icefrac) * ck_sea + cka = cka * icefrac + (1.0_kind_phys - icefrac) * cka_sea + cd = cd * icefrac + (1.0_kind_phys - icefrac) * cd_sea + cda = cda * icefrac + (1.0_kind_phys - icefrac) * cda_sea + end where + + ! `regime` is a categorical variable. Assign according to whether sea ice or sea water is dominant. + where (mask_sea_ice_cell .and. icefrac < 0.5_kind_phys) + regime = regime_sea + end where + end if + + cflx(:, water_vapor_mixing_ratio_index) = qfx(:) + + errmsg = '' + errflg = 0 + end subroutine sf_mynn_compat_run + + !> \section arg_table_sf_mynn_diagnostics_init Argument Table + !! \htmlinclude sf_mynn_diagnostics_init.html + subroutine sf_mynn_diagnostics_init( & + errmsg, errflg) + use cam_history, only: history_add_field + use cam_history_support, only: horiz_only + + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + call history_add_field('sf_mynn_lh', & + 'surface_upward_latent_heat_flux_from_coupler', horiz_only, 'avg', 'W m-2') + call history_add_field('sf_mynn_hfx', & + 'surface_upward_sensible_heat_flux_from_coupler', horiz_only, 'avg', 'W m-2') + call history_add_field('sf_mynn_qfx', & + 'surface_upward_water_vapor_flux', horiz_only, 'avg', 'kg m-2 s-1') + + call history_add_field('sf_mynn_cd', & + 'drag_coefficient_for_momentum_at_10m', horiz_only, 'avg', '1') + call history_add_field('sf_mynn_cda', & + 'drag_coefficient_for_momentum_at_surface_adjacent_layer', horiz_only, 'avg', '1') + call history_add_field('sf_mynn_ck', & + 'bulk_exchange_coefficient_for_enthalpy_at_2m', horiz_only, 'avg', '1') + call history_add_field('sf_mynn_cka', & + 'bulk_exchange_coefficient_for_enthalpy_at_surface_adjacent_layer', horiz_only, 'avg', '1') + + call history_add_field('sf_mynn_q2', & + 'water_vapor_mixing_ratio_wrt_dry_air_at_2m', horiz_only, 'avg', 'kg kg-1') + call history_add_field('sf_mynn_t2', & + 'air_temperature_at_2m', horiz_only, 'avg', 'K') + call history_add_field('sf_mynn_th2', & + 'air_potential_temperature_at_2m', horiz_only, 'avg', 'K') + call history_add_field('sf_mynn_u10', & + 'eastward_wind_at_10m', horiz_only, 'avg', 'm s-1') + call history_add_field('sf_mynn_v10', & + 'northward_wind_at_10m', horiz_only, 'avg', 'm s-1') + + call history_add_field('sf_mynn_ustm', & + 'surface_friction_velocity_assuming_no_correction_for_surface_convective_velocity_scale', horiz_only, 'avg', 'm s-1') + call history_add_field('sf_mynn_wstar', & + 'surface_convective_velocity_scale', horiz_only, 'avg', 'm s-1') + call history_add_field('sf_mynn_qstar', & + 'surface_moisture_scale', horiz_only, 'avg', 'g kg-1') + + errmsg = '' + errflg = 0 + end subroutine sf_mynn_diagnostics_init + + !> \section arg_table_sf_mynn_diagnostics_run Argument Table + !! \htmlinclude sf_mynn_diagnostics_run.html + subroutine sf_mynn_diagnostics_run( & + lh, hfx, qfx, & + cd, cda, ck, cka, & + q2, t2, th2, u10, v10, & + ustm, wstar, qstar, & + errmsg, errflg) + use cam_history, only: history_out_field + use ccpp_kinds, only: kind_phys + + real(kind_phys), intent(in) :: lh(:), hfx(:), qfx(:), & + cd(:), cda(:), ck(:), cka(:), & + q2(:), t2(:), th2(:), u10(:), v10(:), & + ustm(:), wstar(:), qstar(:) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + call history_out_field('sf_mynn_lh', lh) + call history_out_field('sf_mynn_hfx', hfx) + call history_out_field('sf_mynn_qfx', qfx) + + call history_out_field('sf_mynn_cd', cd) + call history_out_field('sf_mynn_cda', cda) + call history_out_field('sf_mynn_ck', ck) + call history_out_field('sf_mynn_cka', cka) + + call history_out_field('sf_mynn_q2', q2) + call history_out_field('sf_mynn_t2', t2) + call history_out_field('sf_mynn_th2', th2) + call history_out_field('sf_mynn_u10', u10) + call history_out_field('sf_mynn_v10', v10) + + call history_out_field('sf_mynn_ustm', ustm) + call history_out_field('sf_mynn_wstar', wstar) + call history_out_field('sf_mynn_qstar', qstar) + + errmsg = '' + errflg = 0 + end subroutine sf_mynn_diagnostics_run +end module sf_mynn_compat diff --git a/schemes/mmm/sf_mynn_compat.meta b/schemes/mmm/sf_mynn_compat.meta new file mode 100644 index 00000000..53463d22 --- /dev/null +++ b/schemes/mmm/sf_mynn_compat.meta @@ -0,0 +1,982 @@ +[ccpp-table-properties] + name = sf_mynn_compat_pre + type = scheme + +[ccpp-arg-table] + name = sf_mynn_compat_pre_run + type = scheme +[ itimestep ] + standard_name = current_timestep_number + units = count + type = integer + dimensions = () + intent = in +[ spp_pbl ] + standard_name = flag_for_stochastically_perturbed_parameterization_for_mynn_scheme + units = flag + type = logical + dimensions = () + intent = in +[ u ] + standard_name = eastward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ v ] + standard_name = northward_wind + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ t ] + standard_name = air_temperature + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ qv ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ p ] + standard_name = air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ dz ] + standard_name = atmosphere_layer_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rho ] + standard_name = air_density + units = kg m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ icefrac ] + standard_name = sea_ice_area_fraction_from_coupler + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ xice_threshold ] + standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme + units = fraction + type = real | kind = kind_phys + dimensions = () + intent = in +[ landfrac ] + standard_name = land_area_fraction_from_coupler + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowhice ] + standard_name = lwe_surface_snow_depth_over_ice_from_coupler + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ snowhland ] + standard_name = lwe_surface_snow_depth_over_land_from_coupler + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ u1d ] + standard_name = eastward_wind_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ v1d ] + standard_name = northward_wind_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ t1d ] + standard_name = air_temperature_at_surface_adjacent_layer + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ qv1d ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_surface_adjacent_layer + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ p1d ] + standard_name = air_pressure_at_surface_adjacent_layer + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ dz8w1d ] + standard_name = atmosphere_layer_thickness_at_surface_adjacent_layer + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ rho1d ] + standard_name = air_density_at_surface_adjacent_layer + units = kg m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ u1d2 ] + standard_name = eastward_wind_at_second_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ v1d2 ] + standard_name = northward_wind_at_second_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ dz2w1d ] + standard_name = atmosphere_layer_thickness_at_second_surface_adjacent_layer + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ chs ] + standard_name = exchange_coefficient_for_heat_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ chs2 ] + standard_name = exchange_coefficient_for_heat_at_2m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cqs2 ] + standard_name = exchange_coefficient_for_moisture_at_2m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cpm ] + standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure_at_surface_adjacent_layer + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ rmol ] + standard_name = reciprocal_of_monin_obukhov_length + units = m-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ znt ] + standard_name = surface_roughness_length + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ust ] + standard_name = surface_friction_velocity + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ zol ] + standard_name = ratio_of_height_at_surface_adjacent_layer_to_monin_obukhov_length + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ mol ] + standard_name = surface_temperature_scale + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ regime ] + standard_name = control_for_pbl_stability_regime + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ psim ] + standard_name = monin_obukhov_similarity_function_for_momentum + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ psih ] + standard_name = monin_obukhov_similarity_function_for_heat + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ qfx ] + standard_name = surface_upward_water_vapor_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flhc ] + standard_name = surface_kinematic_exchange_coefficient_for_heat + units = W m-2 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flqc ] + standard_name = surface_kinematic_exchange_coefficient_for_moisture + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ snowh ] + standard_name = surface_snow_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ qgh ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_surface_adjacent_layer_assuming_saturation + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ qsfc ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_surface + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ gz1oz0 ] + standard_name = ln_ratio_of_height_at_surface_adjacent_layer_to_surface_roughness_length + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ wspd ] + standard_name = wind_speed_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ br ] + standard_name = bulk_richardson_number_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ svp1 ] + standard_name = constant_1_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = out +[ svp2 ] + standard_name = constant_2_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = out +[ svp3 ] + standard_name = constant_3_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = out +[ svpt0 ] + standard_name = constant_4_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = out +[ qcg ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_dry_air_at_surface + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ rstoch1d ] + standard_name = stochastically_perturbed_weights_for_mynn_scheme + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +# ---------- + +[ccpp-table-properties] + name = sf_mynn_compat + type = scheme + dependencies = ccpp_kind_types.F90, mmm_physics/mynn_shared.F90, mmm_physics/sf_mynn.F90 + +[ccpp-arg-table] + name = sf_mynn_compat_init + type = scheme +[ ust ] + standard_name = surface_friction_velocity + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ mol ] + standard_name = surface_temperature_scale + units = K + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ qsfc ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_surface + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = sf_mynn_compat_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ cflx ] + standard_name = surface_upward_ccpp_constituent_fluxes_from_coupler + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, number_of_ccpp_constituents) + intent = inout +[ icefrac ] + standard_name = sea_ice_area_fraction_from_coupler + units = fraction + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ xice_threshold ] + standard_name = sea_ice_area_fraction_threshold_for_mmm_scheme + units = fraction + type = real | kind = kind_phys + dimensions = () + intent = in +[ sst ] + standard_name = sea_surface_temperature_from_coupler + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ u1d ] + standard_name = eastward_wind_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ v1d ] + standard_name = northward_wind_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ t1d ] + standard_name = air_temperature_at_surface_adjacent_layer + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ qv1d ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_surface_adjacent_layer + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ p1d ] + standard_name = air_pressure_at_surface_adjacent_layer + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ dz8w1d ] + standard_name = atmosphere_layer_thickness_at_surface_adjacent_layer + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ rho1d ] + standard_name = air_density_at_surface_adjacent_layer + units = kg m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ u1d2 ] + standard_name = eastward_wind_at_second_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ v1d2 ] + standard_name = northward_wind_at_second_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ dz2w1d ] + standard_name = atmosphere_layer_thickness_at_second_surface_adjacent_layer + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cp ] + standard_name = specific_heat_of_dry_air_at_constant_pressure + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ g ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ rovcp ] + standard_name = ratio_of_dry_air_gas_constant_to_specific_heat_of_dry_air_at_constant_pressure + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ r ] + standard_name = gas_constant_of_dry_air + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ xlv ] + standard_name = latent_heat_of_vaporization_of_water_at_0c + units = J kg-1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ psfcpa ] + standard_name = surface_air_pressure + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ chs ] + standard_name = exchange_coefficient_for_heat_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ chs2 ] + standard_name = exchange_coefficient_for_heat_at_2m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ cqs2 ] + standard_name = exchange_coefficient_for_moisture_at_2m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ cpm ] + standard_name = composition_dependent_specific_heat_of_dry_air_at_constant_pressure_at_surface_adjacent_layer + units = J kg-1 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ pblh ] + standard_name = atmosphere_boundary_layer_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ rmol ] + standard_name = reciprocal_of_monin_obukhov_length + units = m-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ znt ] + standard_name = surface_roughness_length + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ ust ] + standard_name = surface_friction_velocity + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ mavail ] + standard_name = surface_moisture_availability + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ zol ] + standard_name = ratio_of_height_at_surface_adjacent_layer_to_monin_obukhov_length + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ mol ] + standard_name = surface_temperature_scale + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ regime ] + standard_name = control_for_pbl_stability_regime + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ psim ] + standard_name = monin_obukhov_similarity_function_for_momentum + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ psih ] + standard_name = monin_obukhov_similarity_function_for_heat + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ xland ] + standard_name = land_binary_mask_for_mmm_scheme + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ hfx ] + standard_name = surface_upward_sensible_heat_flux_from_coupler + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ qfx ] + standard_name = surface_upward_water_vapor_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ tsk ] + standard_name = skin_temperature_at_surface + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ u10 ] + standard_name = eastward_wind_at_10m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ v10 ] + standard_name = northward_wind_at_10m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ th2 ] + standard_name = air_potential_temperature_at_2m + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ t2 ] + standard_name = air_temperature_at_2m + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ q2 ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_2m + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ flhc ] + standard_name = surface_kinematic_exchange_coefficient_for_heat + units = W m-2 K-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ flqc ] + standard_name = surface_kinematic_exchange_coefficient_for_moisture + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ snowh ] + standard_name = surface_snow_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ qgh ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_surface_adjacent_layer_assuming_saturation + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ qsfc ] + standard_name = water_vapor_mixing_ratio_wrt_moist_air_and_condensed_water_at_surface + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ lh ] + standard_name = surface_upward_latent_heat_flux_from_coupler + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ gz1oz0 ] + standard_name = ln_ratio_of_height_at_surface_adjacent_layer_to_surface_roughness_length + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ wspd ] + standard_name = wind_speed_at_surface_adjacent_layer + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ br ] + standard_name = bulk_richardson_number_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = inout +[ isfflx ] + standard_name = control_for_surface_flux_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = in +[ dx ] + standard_name = characteristic_grid_lengthscale + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ svp1 ] + standard_name = constant_1_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ svp2 ] + standard_name = constant_2_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ svp3 ] + standard_name = constant_3_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ svpt0 ] + standard_name = constant_4_in_formula_for_water_vapor_partial_pressure_assuming_saturation + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ ep1 ] + standard_name = ratio_of_water_vapor_to_dry_air_gas_constants_minus_one + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ ep2 ] + standard_name = ratio_of_water_vapor_to_dry_air_molecular_weights + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ karman ] + standard_name = von_karman_constant + units = 1 + type = real | kind = kind_phys + dimensions = () + intent = in +[ qcg ] + standard_name = cloud_liquid_water_mixing_ratio_wrt_dry_air_at_surface + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ itimestep ] + standard_name = current_timestep_number + units = count + type = integer + dimensions = () + intent = in +[ wstar ] + standard_name = surface_convective_velocity_scale + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ qstar ] + standard_name = surface_moisture_scale + units = g kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ustm ] + standard_name = surface_friction_velocity_assuming_no_correction_for_surface_convective_velocity_scale + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ ck ] + standard_name = bulk_exchange_coefficient_for_enthalpy_at_2m + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cka ] + standard_name = bulk_exchange_coefficient_for_enthalpy_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cd ] + standard_name = drag_coefficient_for_momentum_at_10m + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ cda ] + standard_name = drag_coefficient_for_momentum_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out +[ spp_pbl ] + standard_name = flag_for_stochastically_perturbed_parameterization_for_mynn_scheme + units = flag + type = logical + dimensions = () + intent = in +[ rstoch1d ] + standard_name = stochastically_perturbed_weights_for_mynn_scheme + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ isftcflx ] + standard_name = control_for_surface_roughness_length_over_water_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = in +[ iz0tlnd ] + standard_name = control_for_surface_roughness_length_over_land_for_mmm_scheme + units = 1 + type = integer + dimensions = () + intent = in +[ its ] + standard_name = horizontal_loop_begin + units = count + type = integer + dimensions = () + intent = in +[ ite ] + standard_name = horizontal_loop_end + units = count + type = integer + dimensions = () + intent = in +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +# ---------- + +[ccpp-table-properties] + name = sf_mynn_diagnostics + type = scheme + +[ccpp-arg-table] + name = sf_mynn_diagnostics_init + type = scheme +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +[ccpp-arg-table] + name = sf_mynn_diagnostics_run + type = scheme +[ lh ] + standard_name = surface_upward_latent_heat_flux_from_coupler + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ hfx ] + standard_name = surface_upward_sensible_heat_flux_from_coupler + units = W m-2 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ qfx ] + standard_name = surface_upward_water_vapor_flux + units = kg m-2 s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cd ] + standard_name = drag_coefficient_for_momentum_at_10m + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cda ] + standard_name = drag_coefficient_for_momentum_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ck ] + standard_name = bulk_exchange_coefficient_for_enthalpy_at_2m + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ cka ] + standard_name = bulk_exchange_coefficient_for_enthalpy_at_surface_adjacent_layer + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ q2 ] + standard_name = water_vapor_mixing_ratio_wrt_dry_air_at_2m + units = kg kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ t2 ] + standard_name = air_temperature_at_2m + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ th2 ] + standard_name = air_potential_temperature_at_2m + units = K + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ u10 ] + standard_name = eastward_wind_at_10m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ v10 ] + standard_name = northward_wind_at_10m + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ ustm ] + standard_name = surface_friction_velocity_assuming_no_correction_for_surface_convective_velocity_scale + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ wstar ] + standard_name = surface_convective_velocity_scale + units = m s-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ qstar ] + standard_name = surface_moisture_scale + units = g kg-1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out diff --git a/schemes/utilities/state_converters.F90 b/schemes/utilities/state_converters.F90 index 5d49ec97..a4b5c578 100644 --- a/schemes/utilities/state_converters.F90 +++ b/schemes/utilities/state_converters.F90 @@ -10,9 +10,15 @@ module state_converters public :: temp_to_potential_temp_run public :: potential_temp_to_temp_run - ! Calculate density from equation of state/ideal gas law + ! Calculate dry air density by equation of state/ideal gas law public :: calc_dry_air_ideal_gas_density_run + ! Calculate air density by hydrostatic equation + public :: calc_hydrostatic_air_density_run + + ! Calculate atmosphere layer thickness + public :: calc_atmosphere_layer_thickness_run + ! Calculate exner public :: calc_exner_run @@ -75,7 +81,7 @@ end subroutine potential_temp_to_temp_run subroutine calc_dry_air_ideal_gas_density_run(ncol, nz, rair, pmiddry, temp, rho, errmsg, errflg) integer, intent(in) :: ncol ! Number of columns integer, intent(in) :: nz ! Number of vertical levels - real(kind_phys), intent(in) :: rair(:,:) ! gas constant for dry air (J kg-1) + real(kind_phys), intent(in) :: rair(:,:) ! Gas constant of dry air (J kg-1 K-1) real(kind_phys), intent(in) :: pmiddry(:,:) ! Air pressure of dry air (Pa) real(kind_phys), intent(in) :: temp(:,:) ! Air temperature (K) real(kind_phys), intent(out) :: rho(:,:) ! Dry air density (kg m-3) @@ -93,6 +99,53 @@ subroutine calc_dry_air_ideal_gas_density_run(ncol, nz, rair, pmiddry, temp, rho end subroutine calc_dry_air_ideal_gas_density_run + !> \section arg_table_calc_hydrostatic_air_density_run Argument Table + !! \htmlinclude calc_hydrostatic_air_density_run.html + pure subroutine calc_hydrostatic_air_density_run( & + pdel, gravit, dz, & + rho, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + + real(kind_phys), intent(in) :: pdel(:, :), gravit, dz(:, :) + real(kind_phys), intent(out) :: rho(:, :) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Calculate air density by hydrostatic equation. + rho(:, :) = pdel(:, :) / (gravit * dz(:, :)) + + errmsg = '' + errflg = 0 + end subroutine calc_hydrostatic_air_density_run + + !> \section arg_table_calc_atmosphere_layer_thickness_run Argument Table + !! \htmlinclude calc_atmosphere_layer_thickness_run.html + pure subroutine calc_atmosphere_layer_thickness_run( & + ncol, & + zisfc, & + dz, & + errmsg, errflg) + use ccpp_kinds, only: kind_phys + + integer, intent(in) :: ncol + real(kind_phys), intent(in) :: zisfc(:, :) + real(kind_phys), intent(out) :: dz(:, :) + character(*), intent(out) :: errmsg + integer, intent(out) :: errflg + + integer :: i + + ! In CAM-SIMA, the first vertical index is at top of atmosphere. + ! The last one is at bottom of atmosphere. The resulting `dz` is positive. + do i = 1, ncol + dz(i, :) = zisfc(i, 1:size(zisfc, 2) - 1) - zisfc(i, 2:size(zisfc, 2)) + end do + + errmsg = '' + errflg = 0 + end subroutine calc_atmosphere_layer_thickness_run + !> \section arg_table_calc_exner_run Argument Table !! \htmlinclude calc_exner_run.html subroutine calc_exner_run(ncol, nz, cpair, rair, ref_pres, pmid, exner, & diff --git a/schemes/utilities/state_converters.meta b/schemes/utilities/state_converters.meta index 1bbfc5bb..9389d8f5 100644 --- a/schemes/utilities/state_converters.meta +++ b/schemes/utilities/state_converters.meta @@ -173,6 +173,94 @@ type = integer intent = out +######################################################### +[ccpp-table-properties] + name = calc_hydrostatic_air_density + type = scheme + +[ccpp-arg-table] + name = calc_hydrostatic_air_density_run + type = scheme +[ pdel ] + standard_name = air_pressure_thickness + units = Pa + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ gravit ] + standard_name = standard_gravitational_acceleration + units = m s-2 + type = real | kind = kind_phys + dimensions = () + intent = in +[ dz ] + standard_name = atmosphere_layer_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = in +[ rho ] + standard_name = air_density + units = kg m-3 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + +######################################################### +[ccpp-table-properties] + name = calc_atmosphere_layer_thickness + type = scheme + +[ccpp-arg-table] + name = calc_atmosphere_layer_thickness_run + type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in +[ zisfc ] + standard_name = geopotential_height_wrt_surface_at_interface + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_interface_dimension) + intent = in +[ dz ] + standard_name = atmosphere_layer_thickness + units = m + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent, vertical_layer_dimension) + intent = out +[ errmsg ] + standard_name = ccpp_error_message + long_name = Error message for error handling in CCPP + units = none + type = character | kind = len=* + dimensions = () + intent = out +[ errflg ] + standard_name = ccpp_error_code + long_name = Error flag for error handling in CCPP + units = 1 + type = integer + dimensions = () + intent = out + ######################################################### [ccpp-table-properties] name = calc_exner diff --git a/test/test_suites/suite_convection_permitting.xml b/test/test_suites/suite_convection_permitting.xml index 4c91fda1..884a0c1f 100644 --- a/test/test_suites/suite_convection_permitting.xml +++ b/test/test_suites/suite_convection_permitting.xml @@ -2,7 +2,10 @@ + + calc_atmosphere_layer_thickness calc_exner + calc_hydrostatic_air_density compute_characteristic_grid_length_scale geopotential_height_wrt_sfc_at_if_to_msl geopotential_height_wrt_sfc_to_msl @@ -13,14 +16,22 @@ mmm_physics_compat + + sf_mynn_compat_pre + sf_mynn_compat + sf_mynn_diagnostics + + bl_gwdo_compat_pre bl_gwdo_compat bl_gwdo_diagnostics + cu_ntiedtke_compat_pre cu_ntiedtke_compat cu_ntiedtke_diagnostics + mmm_physics_accumulate_tendencies apply_tendency_of_eastward_wind apply_tendency_of_northward_wind diff --git a/test/unit-test/tests/mmm/CMakeLists.txt b/test/unit-test/tests/mmm/CMakeLists.txt index 69528847..269c7b6e 100644 --- a/test/unit-test/tests/mmm/CMakeLists.txt +++ b/test/unit-test/tests/mmm/CMakeLists.txt @@ -6,3 +6,9 @@ add_pfunit_ctest(mmm_physics_compat_tests LINK_LIBRARIES mmm_physics_compat ) +target_compile_options(mmm_physics_compat_tests + PRIVATE + $<$,$>:-fbacktrace -fcheck=all -ffpe-trap=invalid,overflow,zero -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans -std=f2018 -Wall -Wextra -Wpedantic> + $<$,$>:-check all -fp-model=precise -stand f18 -traceback -warn all> + $<$,$>:-check all -fp-model=precise -stand f18 -traceback -warn all> +) diff --git a/test/unit-test/tests/mmm/mmm_physics_compat_tests.pf b/test/unit-test/tests/mmm/mmm_physics_compat_tests.pf index 0160fd60..d1e3e0af 100644 --- a/test/unit-test/tests/mmm/mmm_physics_compat_tests.pf +++ b/test/unit-test/tests/mmm/mmm_physics_compat_tests.pf @@ -1,3 +1,37 @@ +@test +subroutine test_mmm_physics_compat_init() + use ccpp_kinds, only: kind_phys + use funit + use mmm_physics_compat, only: mmm_physics_compat_init + + integer :: isfflx, isftcflx, iz0tlnd + logical :: spp_pbl + real(kind_phys) :: xice_threshold + character(100) :: errmsg + integer :: errflg + + isfflx = huge(0) + isftcflx = huge(0) + iz0tlnd = huge(0) + spp_pbl = .false. + xice_threshold = huge(0.0_kind_phys) + + call mmm_physics_compat_init( & + isfflx, isftcflx, iz0tlnd, & + spp_pbl, & + xice_threshold, & + errmsg, errflg) + + ! Should set options to the same values as in the MPAS physics driver. + @assertEqual(1, isfflx) + @assertEqual(0, isftcflx) + @assertEqual(0, iz0tlnd) + @assertEqual(.false., spp_pbl) + @assertEqual(0.02_kind_phys, xice_threshold) + @assertEqual('', errmsg) + @assertEqual(0, errflg) +end subroutine test_mmm_physics_compat_init + @test subroutine test_mmm_physics_compat_run() use ccpp_kinds, only: kind_phys @@ -8,11 +42,15 @@ subroutine test_mmm_physics_compat_run() integer :: nstep real(kind_phys) :: dt real(kind_phys) :: theta_curr(ncol, pver), theta_prev(ncol, pver), qv_curr(ncol, pver), qv_prev(ncol, pver) + real(kind_phys) :: icefrac(ncol), xice_threshold, landfrac(ncol) character(256) :: scheme_name real(kind_phys) :: rthdynten(ncol, pver), rqvdynten(ncol, pver) + real(kind_phys) :: xland(ncol) character(100) :: errmsg integer :: errflg + integer :: i + nstep = 0 dt = 20.0_kind_phys @@ -20,41 +58,66 @@ subroutine test_mmm_physics_compat_run() theta_prev(:, :) = 273.15_kind_phys qv_curr(:, :) = 0.03_kind_phys qv_prev(:, :) = 0.01_kind_phys + icefrac(1:50) = 0.0_kind_phys + icefrac(51:100) = [(0.01_kind_phys + real(i - 1, kind_phys) * 0.02_kind_phys, i = 1, 50)] + xice_threshold = 0.02_kind_phys + landfrac(1:50) = [(0.01_kind_phys + real(i - 1, kind_phys) * 0.02_kind_phys, i = 1, 50)] + landfrac(51:100) = 0.0_kind_phys + scheme_name = '' rthdynten(:, :) = huge(0.0_kind_phys) rqvdynten(:, :) = huge(0.0_kind_phys) + xland(:) = huge(0.0_kind_phys) call mmm_physics_compat_run( & nstep, & dt, & theta_curr, theta_prev, qv_curr, qv_prev, & + icefrac, xice_threshold, landfrac, & scheme_name, & rthdynten, rqvdynten, & + xland, & errmsg, errflg) - ! Tendencies should be zero at start (nstep = 0). + ! Should set scheme name. @assertEqual('mmm_physics_compat_run', scheme_name) + ! Tendencies should be zero at start (nstep = 0). @assertEqual(0.0_kind_phys, rthdynten) @assertEqual(0.0_kind_phys, rqvdynten) + ! Should set land mask according to ice and land fractions. + @assertEqual(2.0_kind_phys, xland(1:25)) + @assertEqual(1.0_kind_phys, xland(26:50)) + @assertEqual(2.0_kind_phys, xland(51)) + @assertEqual(1.0_kind_phys, xland(52:100)) @assertEqual('', errmsg) @assertEqual(0, errflg) nstep = 1 + scheme_name = '' rthdynten(:, :) = huge(0.0_kind_phys) rqvdynten(:, :) = huge(0.0_kind_phys) + xland(:) = huge(0.0_kind_phys) call mmm_physics_compat_run( & nstep, & dt, & theta_curr, theta_prev, qv_curr, qv_prev, & + icefrac, xice_threshold, landfrac, & scheme_name, & rthdynten, rqvdynten, & + xland, & errmsg, errflg) - ! Should compute tendencies correctly afterwards (nstep > 0). + ! Should set scheme name. @assertEqual('mmm_physics_compat_run', scheme_name) + ! Should compute tendencies correctly afterwards (nstep > 0). @assertEqual(1.0_kind_phys, rthdynten, epsilon(0.0_kind_phys)) @assertEqual(0.001_kind_phys, rqvdynten, epsilon(0.0_kind_phys)) + ! Should set land mask according to ice and land fractions. + @assertEqual(2.0_kind_phys, xland(1:25)) + @assertEqual(1.0_kind_phys, xland(26:50)) + @assertEqual(2.0_kind_phys, xland(51)) + @assertEqual(1.0_kind_phys, xland(52:100)) @assertEqual('', errmsg) @assertEqual(0, errflg) end subroutine test_mmm_physics_compat_run diff --git a/test/unit-test/tests/utilities/test_state_converters.pf b/test/unit-test/tests/utilities/test_state_converters.pf index bc56f038..dcd7bc8e 100644 --- a/test/unit-test/tests/utilities/test_state_converters.pf +++ b/test/unit-test/tests/utilities/test_state_converters.pf @@ -25,3 +25,63 @@ subroutine test_temp_to_potential_temp() @assertEqual(0, errflg) end subroutine test_temp_to_potential_temp + +@test +subroutine test_calc_hydrostatic_air_density_run() + use ccpp_kinds, only: kind_phys + use funit + use state_converters, only: calc_hydrostatic_air_density_run + + integer, parameter :: ncol = 100, pver = 10 + real(kind_phys) :: pdel(ncol, pver), gravit, dz(ncol, pver) + real(kind_phys) :: rho(ncol, pver) + character(100) :: errmsg + integer :: errflg + + pdel(:, :) = 98.0_kind_phys + gravit = 9.8_kind_phys + dz(:, :) = 10.0_kind_phys + rho(:, :) = huge(0.0_kind_phys) + + call calc_hydrostatic_air_density_run( & + pdel, gravit, dz, & + rho, & + errmsg, errflg) + + ! Should compute air density correctly. + @assertEqual(1.0_kind_phys, rho, spacing(1.0_kind_phys)) + @assertEqual('', errmsg) + @assertEqual(0, errflg) +end subroutine test_calc_hydrostatic_air_density_run + +@test +subroutine test_calc_atmosphere_layer_thickness_run() + use ccpp_kinds, only: kind_phys + use funit + use state_converters, only: calc_atmosphere_layer_thickness_run + + integer, parameter :: ncol = 100, pver = 10, pverp = 11 + real(kind_phys) :: zisfc(ncol, pverp) + real(kind_phys) :: dz(ncol, pver) + character(100) :: errmsg + integer :: errflg + + integer :: i + + do i = 1, pverp + zisfc(:, i) = 10000.0_kind_phys - real(i - 1, kind_phys) * 1000.0_kind_phys + end do + + dz(:, :) = huge(0.0_kind_phys) + + call calc_atmosphere_layer_thickness_run( & + ncol, & + zisfc, & + dz, & + errmsg, errflg) + + ! Should compute atmosphere layer thickness correctly. + @assertEqual(1000.0_kind_phys, dz, spacing(1000.0_kind_phys)) + @assertEqual('', errmsg) + @assertEqual(0, errflg) +end subroutine test_calc_atmosphere_layer_thickness_run