Remove const g parameter -> use CLIMAParameters

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

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,

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,

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

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")

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