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