From 396ad374afda6aa0ab2727b4eeaa12036eacd8c2 Mon Sep 17 00:00:00 2001 From: Charles Kawczynski Date: Mon, 16 Aug 2021 13:35:31 -0700 Subject: [PATCH] Remove const g parameter -> use CLIMAParameters --- integration_tests/utils/Cases.jl | 20 +++++++++++++-- src/EDMF_Environment.jl | 12 ++++++--- src/EDMF_Updrafts.jl | 9 ++++--- src/Surface.jl | 44 +++++++++++++++++++++++++------- src/Turbulence_PrognosticTKE.jl | 7 ++++- src/Variables.jl | 3 ++- src/microphysics_functions.jl | 9 ++++--- src/parameters.jl | 1 - src/surface_functions.jl | 3 ++- src/thermodynamic_functions.jl | 3 ++- src/turbulence_functions.jl | 3 ++- src/types.jl | 10 +++++--- 12 files changed, 93 insertions(+), 31 deletions(-) diff --git a/integration_tests/utils/Cases.jl b/integration_tests/utils/Cases.jl index a2a539903b..b4e9aaf9fe 100644 --- a/integration_tests/utils/Cases.jl +++ b/integration_tests/utils/Cases.jl @@ -2,6 +2,10 @@ module Cases using Random +import CLIMAParameters +const CP = CLIMAParameters +const CPP = CP.Planet + import ..TurbulenceConvection using ..TurbulenceConvection: CasesBase using ..TurbulenceConvection: set_bcs @@ -18,7 +22,6 @@ using ..TurbulenceConvection: dycoms_L using ..TurbulenceConvection: qv_star_c using ..TurbulenceConvection: t_to_thetali_c using ..TurbulenceConvection: rho_c -using ..TurbulenceConvection: g using ..TurbulenceConvection: dycoms_Rd using ..TurbulenceConvection: dycoms_cp using ..TurbulenceConvection: add_ts @@ -180,6 +183,8 @@ function initialize_profiles(self::CasesBase{SoaresCase}, Gr::Grid, GMV::GridMea end function initialize_surface(self::CasesBase{SoaresCase}, Gr::Grid, Ref::ReferenceState) + param_set = TC.parameter_set(Ref) + g = CPP.grav(param_set) self.Sur.zrough = 0.16 #1.0e-4 0.16 is the value specified in the Nieuwstadt paper. self.Sur.Tsurface = 300.0 self.Sur.qsurface = 5.0e-3 @@ -273,6 +278,8 @@ function initialize_profiles(self::CasesBase{Nieuwstadt}, Gr::Grid, GMV::GridMea satadjust(GMV) end function initialize_surface(self::CasesBase{Nieuwstadt}, Gr::Grid, Ref::ReferenceState) + param_set = TC.parameter_set(Ref) + g = CPP.grav(param_set) self.Sur.zrough = 0.16 #1.0e-4 0.16 is the value specified in the Nieuwstadt paper. self.Sur.Tsurface = 300.0 self.Sur.qsurface = 0.0 @@ -528,6 +535,8 @@ function initialize_profiles(self::CasesBase{life_cycle_Tan2018}, Gr::Grid, GMV: end function initialize_surface(self::CasesBase{life_cycle_Tan2018}, Gr::Grid, Ref::ReferenceState) + param_set = TC.parameter_set(Ref) + g = CPP.grav(param_set) self.Sur.zrough = 1.0e-4 # not actually used, but initialized to reasonable value self.Sur.Tsurface = 299.1 * exner_c(Ref.Pg) self.Sur.qsurface = 22.45e-3 # kg/kg @@ -598,6 +607,8 @@ function TC.io(self::CasesBase{life_cycle_Tan2018}, Stats::NetCDFIO_Stats) io(self, Stats, BaseCase()) end function update_surface(self::CasesBase{life_cycle_Tan2018}, GMV::GridMeanVariables, TS::TimeStepping) + param_set = TC.parameter_set(GMV) + g = CPP.grav(param_set) weight = 1.0 weight_factor = 0.01 + 0.99 * (cos(2.0 * π * TS.t / 3600.0) + 1.0) / 2.0 weight = weight * weight_factor @@ -1400,6 +1411,7 @@ function dycoms_sat_adjst(self::CasesBase{DYCOMS_RF01}, p_, thetal_, qt_) end function initialize_profiles(self::CasesBase{DYCOMS_RF01}, Gr::Grid, GMV::GridMeanVariables, Ref::ReferenceState) + param_set = parameter_set(GMV) thetal = TC.center_field(Gr) # helper variable to recalculate temperature ql = TC.center_field(Gr) # DYCOMS case is saturated qi = 0.0 # no ice @@ -1435,7 +1447,7 @@ function initialize_profiles(self::CasesBase{DYCOMS_RF01}, Gr::Grid, GMV::GridMe # buoyancy profile qv = GMV.QT.values[k] - qi - GMV.QL.values[k] rho = rho_c(Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], qv) - GMV.B.values[k] = buoyancy_c(Ref.rho0_half[k], rho) + GMV.B.values[k] = buoyancy_c(param_set, Ref.rho0_half[k], rho) # velocity profile (geostrophic) GMV.U.values[k] = 7.0 @@ -1454,6 +1466,8 @@ function initialize_profiles(self::CasesBase{DYCOMS_RF01}, Gr::Grid, GMV::GridMe end function initialize_surface(self::CasesBase{DYCOMS_RF01}, Gr::Grid, Ref::ReferenceState) + param_set = TC.parameter_set(Ref) + g = CPP.grav(param_set) self.Sur.zrough = 1.0e-4 self.Sur.ustar_fixed = false self.Sur.cm = 0.0011 @@ -1687,6 +1701,8 @@ function initialize_profiles(self::CasesBase{SP}, Gr::Grid, GMV::GridMeanVariabl end function initialize_surface(self::CasesBase{SP}, Gr::Grid, Ref::ReferenceState) + param_set = TC.parameter_set(Ref) + g = CPP.grav(param_set) self.Sur.Gr = Gr self.Sur.Ref = Ref self.Sur.zrough = 0.1 diff --git a/src/EDMF_Environment.jl b/src/EDMF_Environment.jl index 3b440df37d..319730f26b 100644 --- a/src/EDMF_Environment.jl +++ b/src/EDMF_Environment.jl @@ -66,11 +66,12 @@ function env_cloud_diagnostics(self::EnvironmentVariables, Ref::ReferenceState) end function update_EnvVar(self::EnvironmentThermodynamics, k, EnvVar::EnvironmentVariables, T, H, qt, ql, rho) + param_set = parameter_set(EnvVar) EnvVar.T.values[k] = T EnvVar.H.values[k] = H EnvVar.QT.values[k] = qt EnvVar.QL.values[k] = ql - EnvVar.B.values[k] = buoyancy_c(self.Ref.rho0_half[k], rho) + EnvVar.B.values[k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) EnvVar.RH.values[k] = relative_humidity_c(self.Ref.p0_half[k], qt, ql, 0.0, T) return end @@ -99,7 +100,7 @@ function update_cloud_dry(self::EnvironmentThermodynamics, k, EnvVar::Environmen end function saturation_adjustment(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables) - + param_set = parameter_set(EnvVar) sa = eos_struct() mph = mph_struct() @@ -114,7 +115,7 @@ function saturation_adjustment(self::EnvironmentThermodynamics, EnvVar::Environm EnvVar.QT.values[k], EnvVar.QT.values[k] - EnvVar.QL.values[k], ) - EnvVar.B.values[k] = buoyancy_c(self.Ref.rho0_half[k], rho) + EnvVar.B.values[k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) update_cloud_dry( self, @@ -133,6 +134,7 @@ end function sgs_mean(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables, Rain::RainVariables, dt) + param_set = parameter_set(EnvVar) sa = eos_struct() mph = mph_struct() @@ -145,6 +147,7 @@ function sgs_mean(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables, sa = eos(self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k]) # autoconversion and accretion mph = microphysics_rain_src( + param_set, Rain.rain_model, Rain.max_supersaturation, Rain.C_drag, @@ -169,6 +172,7 @@ function sgs_mean(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables, end function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVariables, Rain::RainVariables, dt) + param_set = parameter_set(EnvVar) # TODO: double check this python-> julia translation # a, w = np.polynomial.hermite.hermgauss(self.quadrature_order) a, w = FastGaussQuadrature.gausshermite(self.quadrature_order) @@ -283,6 +287,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari sa = eos(self.Ref.p0_half[k], qt_hat, h_hat) # autoconversion and accretion mph = microphysics_rain_src( + param_set, Rain.rain_model, Rain.max_supersaturation, Rain.C_drag, @@ -368,6 +373,7 @@ function sgs_quadrature(self::EnvironmentThermodynamics, EnvVar::EnvironmentVari # if variance and covariance are zero do the same as in SA_mean sa = eos(self.Ref.p0_half[k], EnvVar.QT.values[k], EnvVar.H.values[k]) mph = microphysics_rain_src( + param_set, Rain.rain_model, Rain.max_supersaturation, Rain.C_drag, diff --git a/src/EDMF_Updrafts.jl b/src/EDMF_Updrafts.jl index 928d60bf61..8d8b40a347 100644 --- a/src/EDMF_Updrafts.jl +++ b/src/EDMF_Updrafts.jl @@ -359,6 +359,7 @@ function buoyancy( grid = self.Gr kc_surf = kc_surface(grid) + param_set = parameter_set(GMV) UpdVar.Area.bulkvalues .= up_sum(UpdVar.Area.values) @@ -368,7 +369,7 @@ function buoyancy( if UpdVar.Area.values[i, k] > 0.0 qv = UpdVar.QT.values[i, k] - UpdVar.QL.values[i, k] rho = rho_c(self.Ref.p0_half[k], UpdVar.T.values[i, k], UpdVar.QT.values[i, k], qv) - UpdVar.B.values[i, k] = buoyancy_c(self.Ref.rho0_half[k], rho) + UpdVar.B.values[i, k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) else UpdVar.B.values[i, k] = EnvVar.B.values[k] end @@ -390,7 +391,7 @@ function buoyancy( h = UpdVar.H.values[i, k] t = UpdVar.T.values[i, k] rho = rho_c(self.Ref.p0_half[k], t, qt, qv) - UpdVar.B.values[i, k] = buoyancy_c(self.Ref.rho0_half[k], rho) + UpdVar.B.values[i, k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) UpdVar.RH.values[i, k] = relative_humidity_c(self.Ref.p0_half[k], qt, qt - qv, 0.0, t) elseif UpdVar.Area.values[i, k - 1] > 0.0 && k > kc_surf qt = UpdVar.QT.values[i, k - 1] @@ -399,7 +400,7 @@ function buoyancy( qv = qt - sa.ql t = sa.T rho = rho_c(self.Ref.p0_half[k], t, qt, qv) - UpdVar.B.values[i, k] = buoyancy_c(self.Ref.rho0_half[k], rho) + UpdVar.B.values[i, k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) UpdVar.RH.values[i, k] = relative_humidity_c(self.Ref.p0_half[k], qt, qt - qv, 0.0, t) else UpdVar.B.values[i, k] = EnvVar.B.values[k] @@ -431,12 +432,14 @@ function microphysics(self::UpdraftThermodynamics, UpdVar::UpdraftVariables, Rai rst = rain_struct() mph = mph_struct() sa = eos_struct() + param_set = parameter_set() @inbounds for i in xrange(self.n_updraft) @inbounds for k in center_indicies(self.Gr) # autoconversion and accretion mph = microphysics_rain_src( + param_set, Rain.rain_model, Rain.max_supersaturation, Rain.C_drag, diff --git a/src/Surface.jl b/src/Surface.jl index adca9886aa..db44eadf54 100644 --- a/src/Surface.jl +++ b/src/Surface.jl @@ -4,13 +4,14 @@ update(self::SurfaceBase, GMV::GridMeanVariables, ::BaseCase) = nothing function free_convection_windspeed(self::SurfaceBase, GMV::GridMeanVariables, ::BaseCase) theta_rho = center_field(self.Gr) + param_set = parameter_set(GMV) # Need to get theta_rho @inbounds for k in center_indicies(self.Gr) qv = GMV.QT.values[k] - GMV.QL.values[k] theta_rho[k] = theta_rho_c(self.Ref.p0_half[k], GMV.T.values[k], GMV.QT.values[k], qv) end - zi = get_inversion(theta_rho, GMV.U.values, GMV.V.values, self.Gr, self.Ri_bulk_crit) + zi = get_inversion(param_set, theta_rho, GMV.U.values, GMV.V.values, self.Gr, self.Ri_bulk_crit) wstar = get_wstar(self.bflux, zi) # yair here zi in TRMM should be adjusted self.windspeed = sqrt(self.windspeed * self.windspeed + (1.2 * wstar) * (1.2 * wstar)) return @@ -32,6 +33,7 @@ free_convection_windspeed(self::SurfaceBase{SurfaceNone}, GMV::GridMeanVariables initialize(self::SurfaceBase{SurfaceFixedFlux}) = nothing function update(self::SurfaceBase{SurfaceFixedFlux}, GMV::GridMeanVariables) + param_set = parameter_set(GMV) kc_surf = kc_surface(self.Gr) kf_surf = kf_surface(self.Gr) rho_tflux = self.shf / (cpm_c(self.qsurface)) @@ -40,8 +42,14 @@ function update(self::SurfaceBase{SurfaceFixedFlux}, GMV::GridMeanVariables) self.rho_qtflux = self.lhf / (latent_heat(self.Tsurface)) self.rho_hflux = rho_tflux / exner_c(self.Ref.Pg) - self.bflux = - buoyancy_flux(self.shf, self.lhf, GMV.T.values[kc_surf], GMV.QT.values[kc_surf], self.Ref.alpha0[kf_surf]) + self.bflux = buoyancy_flux( + param_set, + self.shf, + self.lhf, + GMV.T.values[kc_surf], + GMV.QT.values[kc_surf], + self.Ref.alpha0[kf_surf], + ) if !self.ustar_fixed # Correction to windspeed for free convective cases (Beljaars, QJRMS (1994), 121, pp. 255-270) @@ -85,6 +93,7 @@ function initialize(self::SurfaceBase{SurfaceFixedCoeffs}) end function update(self::SurfaceBase{SurfaceFixedCoeffs}, GMV::GridMeanVariables) + param_set = parameter_set(GMV) kc_surf = kc_surface(self.Gr) kf_surf = kf_surface(self.Gr) windspeed = @@ -98,8 +107,14 @@ function update(self::SurfaceBase{SurfaceFixedCoeffs}, GMV::GridMeanVariables) -self.ch * windspeed * (GMV.H.values[kc_surf] - self.Tsurface / exner_c(self.Ref.Pg)) * self.Ref.rho0[kf_surf] self.shf = cp_ * self.rho_hflux - self.bflux = - buoyancy_flux(self.shf, self.lhf, GMV.T.values[kc_surf], GMV.QT.values[kc_surf], self.Ref.alpha0[kf_surf]) + self.bflux = buoyancy_flux( + param_set, + self.shf, + self.lhf, + GMV.T.values[kc_surf], + GMV.QT.values[kc_surf], + self.Ref.alpha0[kf_surf], + ) self.ustar = sqrt(self.cm) * windspeed # CK--testing this--EDMF scheme checks greater or less than zero, @@ -122,6 +137,8 @@ end initialize(self::SurfaceBase{SurfaceMoninObukhov}) = nothing function update(self::SurfaceBase{SurfaceMoninObukhov}, GMV::GridMeanVariables) + param_set = parameter_set(GMV) + g = CPP.grav(param_set) kc_surf = kc_surface(self.Gr) kf_surf = kf_surface(self.Gr) self.qsurface = qv_star_t(self.Ref.Pg, self.Tsurface) @@ -148,8 +165,14 @@ function update(self::SurfaceBase{SurfaceMoninObukhov}, GMV::GridMeanVariables) self.shf = cpm_c(GMV.QT.values[kc_surf]) * self.rho_hflux - self.bflux = - buoyancy_flux(self.shf, self.lhf, GMV.T.values[kc_surf], GMV.QT.values[kc_surf], self.Ref.alpha0[kf_surf]) + self.bflux = buoyancy_flux( + param_set, + self.shf, + self.lhf, + GMV.T.values[kc_surf], + GMV.QT.values[kc_surf], + self.Ref.alpha0[kf_surf], + ) self.ustar = sqrt(self.cm) * self.windspeed # CK--testing this--EDMF scheme checks greater or less than zero, if abs(self.bflux) < 1e-10 @@ -169,6 +192,8 @@ end initialize(self::SurfaceBase{SurfaceMoninObukhovDry}) = nothing function update(self::SurfaceBase{SurfaceMoninObukhovDry}, GMV::GridMeanVariables) + param_set = parameter_set(GMV) + g = CPP.grav(param_set) self.qsurface = qv_star_t(self.Ref.Pg, self.Tsurface) kc_surf = kc_surface(self.Gr) kf_surf = kf_surface(self.Gr) @@ -195,7 +220,7 @@ function update(self::SurfaceBase{SurfaceMoninObukhovDry}, GMV::GridMeanVariable self.shf = cpm_c(0.0) * self.rho_hflux - self.bflux = buoyancy_flux(self.shf, self.lhf, GMV.T.values[kc_surf], 0.0, self.Ref.alpha0[kf_surf]) + self.bflux = buoyancy_flux(param_set, self.shf, self.lhf, GMV.T.values[kc_surf], 0.0, self.Ref.alpha0[kf_surf]) self.ustar = sqrt(self.cm) * self.windspeed # CK--testing this--EDMF scheme checks greater or less than zero, if abs(self.bflux) < 1e-10 @@ -213,13 +238,14 @@ free_convection_windspeed(self::SurfaceBase{SurfaceMoninObukhovDry}, GMV::GridMe initialize(self::SurfaceBase{SurfaceSullivanPatton}) = nothing function update(self::SurfaceBase{SurfaceSullivanPatton}, GMV::GridMeanVariables) + param_set = parameter_set(GMV) + g = CPP.grav(param_set) kc_surf = kc_surface(self.Gr) kf_surf = kf_surface(self.Gr) zb = self.Gr.z_half[kc_surf] theta_rho_g = theta_rho_c(self.Ref.Pg, self.Tsurface, self.qsurface, self.qsurface) theta_rho_b = theta_rho_c(self.Ref.p0_half[kc_surf], GMV.T.values[kc_surf], self.qsurface, self.qsurface) lv = latent_heat(GMV.T.values[kc_surf]) - g = 9.81 T0 = self.Ref.p0_half[kc_surf] * self.Ref.alpha0_half[kc_surf] / Rd theta_flux = 0.24 diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 6a6be5013d..f41ef50249 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -342,6 +342,8 @@ end function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariables) grid = get_grid(self) + param_set = parameter_set(GMV) + g = CPP.grav(param_set) ref_state = reference_state(self) kc_surf = kc_surface(grid) tau = get_mixing_tau(self.zi, self.wstar) @@ -1443,6 +1445,8 @@ end function compute_tke_buoy(self::EDMF_PrognosticTKE, GMV::GridMeanVariables) grid = get_grid(self) ref_state = reference_state(self) + param_set = parameter_set(GMV) + g = CPP.grav(param_set) grad_thl_minus = 0.0 grad_qt_minus = 0.0 grad_thl_plus = 0 @@ -2279,6 +2283,7 @@ function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, opti theta_rho = center_field(self.Gr) ∇θ_liq_max = 0.0 kc_surf = kc_surface(self.Gr) + param_set = parameter_set(GMV) @inbounds for k in real_center_indicies(self.Gr) qv = GMV.QT.values[k] - GMV.QL.values[k] @@ -2303,7 +2308,7 @@ function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, opti end end elseif option == "critical_Ri" - self.zi = get_inversion(theta_rho, GMV.U.values, GMV.V.values, self.Gr, Ri_bulk_crit(self)) + self.zi = get_inversion(param_set, theta_rho, GMV.U.values, GMV.V.values, self.Gr, Ri_bulk_crit(self)) else error("INVERSION HEIGHT OPTION NOT RECOGNIZED") diff --git a/src/Variables.jl b/src/Variables.jl index 2190f0856d..ea6d44da47 100644 --- a/src/Variables.jl +++ b/src/Variables.jl @@ -113,6 +113,7 @@ end function satadjust(self::GridMeanVariables) sa = eos_struct() + param_set = parameter_set(self) @inbounds for k in center_indicies(self.Gr) h = self.H.values[k] qt = self.QT.values[k] @@ -122,7 +123,7 @@ function satadjust(self::GridMeanVariables) self.T.values[k] = sa.T qv = qt - sa.ql rho = rho_c(p0, sa.T, qt, qv) - self.B.values[k] = buoyancy_c(self.Ref.rho0_half[k], rho) + self.B.values[k] = buoyancy_c(param_set, self.Ref.rho0_half[k], rho) self.RH.values[k] = relative_humidity_c(self.Ref.p0_half[k], qt, qt - qv, 0.0, self.T.values[k]) end return diff --git a/src/microphysics_functions.jl b/src/microphysics_functions.jl index 3ba03e0bc3..778b190e6f 100644 --- a/src/microphysics_functions.jl +++ b/src/microphysics_functions.jl @@ -57,12 +57,12 @@ function conv_q_vap_to_q_liq(param_set::APS, q_sat_liq, q_liq) return (q_sat_liq - q_liq) / CPMP.τ_cond_evap(param_set) end -function conv_q_liq_to_q_rai_acnv(q_liq_threshold, tau_acnv, q_liq) +function conv_q_liq_to_q_rai_acnv(param_set, q_liq_threshold, tau_acnv, q_liq) return max(0.0, q_liq - q_liq_threshold) / tau_acnv end -function conv_q_liq_to_q_rai_accr(C_drag, MP_n_0, E_col, q_liq, q_rai, rho) +function conv_q_liq_to_q_rai_accr(param_set, C_drag, MP_n_0, E_col, q_liq, q_rai, rho) v_c = terminal_velocity_single_drop_coeff(C_drag, rho) gamma_7_2 = 3.3233509704478426 @@ -102,6 +102,7 @@ return rates: qr_src, thl_rain_src """ function microphysics_rain_src( + param_set::APS, rain_model, max_supersaturation, C_drag, @@ -145,8 +146,8 @@ function microphysics_rain_src( _ret.qr_src = min( ql, ( - conv_q_liq_to_q_rai_acnv(q_liq_threshold, tau_acnv, ql) + - conv_q_liq_to_q_rai_accr(C_drag, MP_n_0, E_col, ql, qr, rho) + conv_q_liq_to_q_rai_acnv(param_set, q_liq_threshold, tau_acnv, ql) + + conv_q_liq_to_q_rai_accr(param_set, C_drag, MP_n_0, E_col, ql, qr, rho) ) * dt, ) end diff --git a/src/parameters.jl b/src/parameters.jl index 692eb0ff36..1c2acdb15f 100644 --- a/src/parameters.jl +++ b/src/parameters.jl @@ -1,5 +1,4 @@ #Adapated from PyCLES: https://github.com/pressel/pycles -const g = 9.80665 const Rd = 287.1 const Rv = 461.5 const eps_v = 0.62210184182 # Rd / Rv diff --git a/src/surface_functions.jl b/src/surface_functions.jl index 1e1ea9ca2d..8978fdf54e 100644 --- a/src/surface_functions.jl +++ b/src/surface_functions.jl @@ -1,5 +1,6 @@ -function buoyancy_flux(shf, lhf, T_b, qt_b, alpha0_0) +function buoyancy_flux(param_set, shf, lhf, T_b, qt_b, alpha0_0) + g = CPP.grav(param_set) cp_ = cpm_c(qt_b) lv = latent_heat(T_b) return (g * alpha0_0 / cp_ / T_b * (shf + (eps_vi - 1.0) * cp_ * T_b * lhf / lv)) diff --git a/src/thermodynamic_functions.jl b/src/thermodynamic_functions.jl index 59c55303ba..73c75dd2fd 100644 --- a/src/thermodynamic_functions.jl +++ b/src/thermodynamic_functions.jl @@ -54,7 +54,8 @@ function relative_humidity_c(p0, qt, ql, qi, T) return 100.0 * pv / pv_star_ end -function buoyancy_c(rho0, rho) +function buoyancy_c(param_set, rho0, rho) + g = CPP.grav(param_set) return g * (rho0 - rho) / rho0 end function qv_star_c(p0, qt, pv) diff --git a/src/turbulence_functions.jl b/src/turbulence_functions.jl index 0bc4b05eda..5fc552a52b 100644 --- a/src/turbulence_functions.jl +++ b/src/turbulence_functions.jl @@ -4,7 +4,8 @@ function get_wstar(bflux, zi) end # BL height -function get_inversion(theta_rho, u, v, grid::Grid, Ri_bulk_crit) +function get_inversion(param_set, theta_rho, u, v, grid::Grid, Ri_bulk_crit) + g = CPP.grav(param_set) kc_surf = kc_surface(grid) kmin = kc_surf theta_rho_b = theta_rho[kmin] diff --git a/src/types.jl b/src/types.jl index fe227042f7..11d02f9aaf 100644 --- a/src/types.jl +++ b/src/types.jl @@ -93,7 +93,8 @@ struct RainVariable{T} end end -Base.@kwdef mutable struct RainVariables +Base.@kwdef mutable struct RainVariables{PS} + param_set::PS rain_model::String = "default_rain_model" max_supersaturation::Float64 C_drag::Float64 @@ -115,7 +116,7 @@ Base.@kwdef mutable struct RainVariables Env_QR::RainVariable Env_RainArea::RainVariable end -function RainVariables(namelist, Gr::Grid) +function RainVariables(namelist, Gr::Grid, param_set::APS) QR = RainVariable(Gr, "qr_mean", "kg/kg") # temporary variables for diagnostics to know where the rain is coming from @@ -147,7 +148,8 @@ function RainVariables(namelist, Gr::Grid) a_vent = parse_namelist(namelist, "microphysics", "a_vent"; default = 1.5) b_vent = parse_namelist(namelist, "microphysics", "b_vent"; default = 0.53) - return RainVariables(; + return RainVariables{typeof(param_set)}(; + param_set, rain_model, max_supersaturation, C_drag, @@ -863,7 +865,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2} minimum_area = 1e-5 # Create the class for rain - Rain = RainVariables(namelist, Gr) + Rain = RainVariables(namelist, Gr, param_set) # Create the updraft variable class (major diagnostic and prognostic variables) UpdVar = UpdraftVariables(n_updrafts, namelist, Gr)