Skip to content

Commit

Permalink
Merge #80
Browse files Browse the repository at this point in the history
80: merge turbulence and TurbProgTKE r=charleskawczynski a=yairchn

this PR deletes Turbulence.jl and moves the relevant functions to TurbulencePrognosticTKE.jl 

Co-authored-by: yairchn <yairchn@caltech.edu>
  • Loading branch information
bors[bot] and yairchn committed Aug 4, 2021
2 parents 4b0ef91 + 4ec5d64 commit 46e8471
Show file tree
Hide file tree
Showing 5 changed files with 88 additions and 117 deletions.
2 changes: 1 addition & 1 deletion integration_tests/utils/main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ function Simulation1d(namelist)
Ref = TurbulenceConvection.ReferenceState(Gr, param_set)
GMV = TurbulenceConvection.GridMeanVariables(namelist, Gr, Ref, param_set)
Case = Cases.CasesFactory(namelist, Gr, Ref)
Turb = TurbulenceConvection.ParameterizationFactory(namelist, Gr, Ref, param_set)
Turb = TurbulenceConvection.EDMF_PrognosticTKE(namelist, Gr, Ref, param_set)
TS = TurbulenceConvection.TimeStepping(namelist)
Stats = TurbulenceConvection.NetCDFIO_Stats(namelist, Gr)
return Simulation1d(Gr, Ref, GMV, Case, Turb, TS, Stats)
Expand Down
62 changes: 0 additions & 62 deletions src/Turbulence.jl

This file was deleted.

1 change: 0 additions & 1 deletion src/TurbulenceConvection.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,6 @@ include("Variables.jl")
include("EDMF_Rain.jl")
include("EDMF_Environment.jl")
include("EDMF_Updrafts.jl")
include("Turbulence.jl")
include("Turbulence_PrognosticTKE.jl")
include("Forcing.jl")
include("Radiation.jl")
Expand Down
68 changes: 57 additions & 11 deletions src/Turbulence_PrognosticTKE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -251,7 +251,7 @@ function update(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBas
grid = get_grid(self)
update_inversion(self, GMV, Case.inversion_option)
compute_pressure_plume_spacing(self, GMV, Case)
self.wstar = get_wstar(Case.Sur.bflux, self.base.zi)
self.wstar = get_wstar(Case.Sur.bflux, self.zi)
if TS.nstep == 0
decompose_environment(self, GMV, "values")

Expand Down Expand Up @@ -317,7 +317,19 @@ function update(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBas
GMV.cloud_cover = min(self.EnvVar.cloud_cover + sum(self.UpdVar.cloud_cover), 1)
# Back out the tendencies of the grid mean variables for the whole timestep
# by differencing GMV.new and GMV.values
update(self.base, GMV, Case, TS)
update_turbulence(self, GMV, Case, TS)
return
end

function update_turbulence(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase, TS::TimeStepping)

@inbounds for k in real_center_indicies(self.Gr)
GMV.H.tendencies[k] += (GMV.H.new[k] - GMV.H.values[k]) * TS.dti
GMV.QT.tendencies[k] += (GMV.QT.new[k] - GMV.QT.values[k]) * TS.dti
GMV.U.tendencies[k] += (GMV.U.new[k] - GMV.U.values[k]) * TS.dti
GMV.V.tendencies[k] += (GMV.V.new[k] - GMV.V.values[k]) * TS.dti
end

return
end

Expand Down Expand Up @@ -368,17 +380,12 @@ function compute_prognostic_updrafts(
return
end

function update_inversion(self, GMV::GridMeanVariables, option)
update_inversion(self.base, GMV, option)
return
end

function compute_mixing_length(self, obukhov_length, ustar, GMV::GridMeanVariables)

grid = get_grid(self)
ref_state = reference_state(self)
gw = grid.gw
tau = get_mixing_tau(self.base.zi, self.wstar)
tau = get_mixing_tau(self.zi, self.wstar)
l = pyzeros(3)
m_eps = 1.0e-9 # Epsilon to avoid zero
@inbounds for k in real_center_indicies(grid)
Expand Down Expand Up @@ -578,7 +585,7 @@ end
function set_updraft_surface_bc(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase)

update_inversion(self, GMV, Case.inversion_option)
self.wstar = get_wstar(Case.Sur.bflux, self.base.zi)
self.wstar = get_wstar(Case.Sur.bflux, self.zi)

gw = get_grid(self).gw
dzi = get_grid(self).dzi
Expand Down Expand Up @@ -955,7 +962,7 @@ function compute_turbulent_entrainment(self::EDMF_PrognosticTKE, GMV::GridMeanVa

grid = get_grid(self)
ref_state = reference_state(self)
tau = get_mixing_tau(self.base.zi, self.wstar)
tau = get_mixing_tau(self.zi, self.wstar)

@inbounds for i in xrange(self.n_updrafts)
@inbounds for k in real_center_indicies(grid)
Expand Down Expand Up @@ -1914,7 +1921,7 @@ function initialize_covariance(self::EDMF_PrognosticTKE, GMV::GridMeanVariables,

ws = self.wstar
us = Case.Sur.ustar
zs = self.base.zi
zs = self.zi
grid = get_grid(self)

reset_surface_covariance(self, GMV, Case)
Expand Down Expand Up @@ -2498,3 +2505,42 @@ function GMV_third_m(
end
return
end


# Update the diagnosis of the inversion height, using the maximum temperature gradient method
function update_inversion(self::EDMF_PrognosticTKE, GMV::GridMeanVariables, option)
theta_rho = center_field(self.Gr)
∇θ_liq_max = 0.0
k_fi = first_center(self.Gr)

@inbounds for k in real_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


if option == "theta_rho"
@inbounds for k in real_center_indicies(self.Gr)
if theta_rho[k] > theta_rho[k_fi]
self.zi = self.Gr.z_half[k]
break
end
end
elseif option == "thetal_maxgrad"

@inbounds for k in real_center_indicies(self.Gr)
∇θ_liq = ∇_upwind(GMV.THL.values, self.Gr, k)
if ∇θ_liq > ∇θ_liq_max
∇θ_liq_max = ∇θ_liq
self.zi = self.Gr.z[k]
end
end
elseif option == "critical_Ri"
self.zi = get_inversion(theta_rho, GMV.U.values, GMV.V.values, self.Gr, Ri_bulk_crit(self))

else
error("INVERSION HEIGHT OPTION NOT RECOGNIZED")
end

return
end
72 changes: 30 additions & 42 deletions src/types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -643,28 +643,6 @@ struct EnvironmentThermodynamics{A1}
end
end

mutable struct ParameterizationBase{T}
turbulence_tendency::T
Gr::Grid
Ref::ReferenceState
KM::VariableDiagnostic
KH::VariableDiagnostic
prandtl_number::Float64
Ri_bulk_crit::Float64
zi::Float64
# A base class common to all turbulence parameterizations
function ParameterizationBase(namelist, Gr::Grid, Ref::ReferenceState)
turbulence_tendency = center_field(Gr)
KM = VariableDiagnostic(Gr, "half", "scalar", "sym", "diffusivity", "m^2/s") # eddy viscosity
KH = VariableDiagnostic(Gr, "half", "scalar", "sym", "viscosity", "m^2/s") # eddy diffusivity
# get values from namelist
prandtl_number = namelist["turbulence"]["prandtl_number_0"]
Ri_bulk_crit = namelist["turbulence"]["Ri_bulk_crit"]

return new{typeof(turbulence_tendency)}(turbulence_tendency, Gr, Ref, KM, KH, prandtl_number, Ri_bulk_crit, 0)
end
end

# SurfaceMoninObukhovDry:
# Needed for dry cases (qt=0). They have to be initialized with nonzero qtg for the
# reference profiles. This surface subroutine sets the latent heat flux to zero
Expand Down Expand Up @@ -773,15 +751,15 @@ Base.@kwdef mutable struct CasesBase{T}
shf0::Float64 = 0
end

struct SimilarityED{PS}
param_set::PS
base::ParameterizationBase
extrapolate_buoyancy::Bool
end

mutable struct EDMF_PrognosticTKE{PS, A1, A2}
param_set::PS
base::ParameterizationBase
turbulence_tendency::A1
Gr::Grid
Ref::ReferenceState
KM::VariableDiagnostic
KH::VariableDiagnostic
Ri_bulk_crit::Float64
zi::Float64
n_updrafts::Int
use_const_plume_spacing::Bool
entr_detr_fp::Function
Expand Down Expand Up @@ -870,8 +848,13 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
detr_surface_bc::Float64
dt_upd::Float64
function EDMF_PrognosticTKE(namelist, Gr::Grid, Ref::ReferenceState, param_set::PS) where {PS}
# Initialize the base parameterization class
base = ParameterizationBase(namelist, Gr, Ref)
turbulence_tendency = center_field(Gr)
KM = VariableDiagnostic(Gr, "half", "scalar", "sym", "diffusivity", "m^2/s") # eddy viscosity
KH = VariableDiagnostic(Gr, "half", "scalar", "sym", "viscosity", "m^2/s") # eddy diffusivity
# get values from namelist
prandtl_number = namelist["turbulence"]["prandtl_number_0"]
Ri_bulk_crit = namelist["turbulence"]["Ri_bulk_crit"]
zi = 0.0

# Set the number of updrafts (1)
n_updrafts = parse_namelist(namelist, "turbulence", "EDMF_PrognosticTKE", "updraft_number"; default = 1)
Expand Down Expand Up @@ -1073,7 +1056,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
# Added by Ignacio : Length scheme in use (mls), and smooth min effect (ml_ratio)
# Variable Prandtl number initialized as neutral value.
ones_vec = center_field(Gr)
prandtl_nvec = base.prandtl_number .* ones_vec
prandtl_nvec = prandtl_number .* ones_vec
mls = center_field(Gr)
ml_ratio = center_field(Gr)
l_entdet = center_field(Gr)
Expand All @@ -1086,7 +1069,13 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
A2 = typeof(horizontal_KM)
return new{PS, A1, A2}(
param_set,
base,
turbulence_tendency,
Gr,
Ref,
KM,
KH,
Ri_bulk_crit,
zi,
n_updrafts,
use_const_plume_spacing,
entr_detr_fp,
Expand Down Expand Up @@ -1165,7 +1154,7 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
diffusive_flux_v,
massflux_tke,
prandtl_nvec,
base.prandtl_number,
prandtl_number,
mls,
ml_ratio,
l_entdet,
Expand All @@ -1177,13 +1166,12 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2}
)
end
end
get_grid(edmf::EDMF_PrognosticTKE) = edmf.base.Gr
get_grid(edmf::EDMF_PrognosticTKE) = edmf.Gr
get_grid(obj) = obj.Gr
reference_state(edmf::EDMF_PrognosticTKE) = edmf.base.Ref
prandtl_number(edmf::EDMF_PrognosticTKE) = edmf.base.prandtl_number
turbulence_tendency(edmf::EDMF_PrognosticTKE) = edmf.base.turbulence_tendency
diffusivity_m(edmf::EDMF_PrognosticTKE) = edmf.base.KM
diffusivity_h(edmf::EDMF_PrognosticTKE) = edmf.base.KH
Ri_bulk_crit(edmf::EDMF_PrognosticTKE) = edmf.base.Ri_bulk_crit
Ri_bulk_crit(base::ParameterizationBase) = base.Ri_bulk_crit
reference_state(edmf::EDMF_PrognosticTKE) = edmf.Ref
prandtl_number(edmf::EDMF_PrognosticTKE) = edmf.prandtl_number
turbulence_tendency(edmf::EDMF_PrognosticTKE) = edmf.turbulence_tendency
diffusivity_m(edmf::EDMF_PrognosticTKE) = edmf.KM
diffusivity_h(edmf::EDMF_PrognosticTKE) = edmf.KH
Ri_bulk_crit(edmf::EDMF_PrognosticTKE) = edmf.Ri_bulk_crit
parameter_set(obj) = obj.param_set

0 comments on commit 46e8471

Please sign in to comment.