implements TOA emdf model

cleans up files and adjusts config forcing pathing and height

Debugging

working scm adding insolation and cos_zenith

updates diagnostics and removes global variable call

updates config link

config/model_configs/prognostic_edmfx_gcmdriven_fullcolumn.yml

@@ -0,0 +1,48 @@
+initial_condition: "GCM"
+external_forcing: "GCM"
+external_forcing_file: "/central/groups/esm/zhaoyi/GCMForcedLES/forcing/corrected/HadGEM2-A_amip.2004-2008.07.nc" 
+surface_setup: "GCM"
+turbconv: "prognostic_edmfx"
+implicit_diffusion: true
+implicit_sgs_advection: false
+approximate_linear_solve_iters: 2
+max_newton_iters_ode: 1
+edmfx_upwinding: first_order
+edmfx_entr_model: "Generalized"
+edmfx_detr_model: "Generalized"
+edmfx_sgs_mass_flux: true
+edmfx_sgs_diffusive_flux: true
+edmfx_nh_pressure: true
+edmfx_filter: true
+prognostic_tke: true
+moist: "equil"
+config: "column"
+z_max: 40e3
+z_elem: 60
+z_stretch: true
+dz_bottom: 30
+dz_top: 3000
+perturb_initstate: false
+dt: "10secs"
+t_end: "48hours"
+dt_save_state_to_disk: "6hours"
+call_cloud_diagnostics_per_stage : true
+output_dir: output/gcm_driven_scm
+toml: [toml/prognostic_edmfx_bomex.toml]
+netcdf_output_at_levels: true
+netcdf_interpolation_num_points: [2, 2, 60]
+output_default_diagnostics: false
+rad: allskywithclear
+insolation: "gcmdriven"
+diagnostics:
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
+    period: 10mins
+  - short_name: [arup, waup, taup, thetaaup, haup, husup, hurup, clwup, cliup, waen, taen, thetaaen, haen, husen, huren, clwen, clien, tke]
+    period: 10mins
+  - short_name: [entr, detr, lmix, bgrad, strain, edt, evu]
+    period: 10mins
+  - short_name: [rlut, rlutcs, rsut, rsutcs]
+    period: 10mins
+  - reduction_time: max
+    short_name: tke
+    period: 10mins

src/cache/cache.jl

@@ -214,7 +214,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
     precipitation = precipitation_cache(Y, atmos)
     subsidence = subsidence_cache(Y, atmos)
     large_scale_advection = large_scale_advection_cache(Y, atmos)
-    external_forcing = external_forcing_cache(Y, atmos)
+    external_forcing = external_forcing_cache(Y, atmos, params)
     edmf_coriolis = edmf_coriolis_cache(Y, atmos)
     forcing = forcing_cache(Y, atmos)
     non_orographic_gravity_wave = non_orographic_gravity_wave_cache(Y, atmos)

src/callbacks/callbacks.jl

@@ -190,6 +190,13 @@ function set_insolation_variables!(Y, p, t, ::RCEMIPIIInsolation)
     rrtmgp_model.weighted_irradiance .= FT(551.58)
 end
 
+function set_insolation_variables!(Y, p, t, ::GCMDrivenInsolation)
+    FT = Spaces.undertype(axes(Y.c))
+    (; rrtmgp_model) = p.radiation
+    rrtmgp_model.cos_zenith .= p.external_forcing.cos_zenith
+    rrtmgp_model.weighted_irradiance .= p.external_forcing.insolation
+end
+
 function set_insolation_variables!(Y, p, t, ::IdealizedInsolation)
     FT = Spaces.undertype(axes(Y.c))
     bottom_coords = Fields.coordinate_field(Spaces.level(Y.c, 1))

src/initial_conditions/InitialConditions.jl

@@ -15,7 +15,8 @@ import ..PrognosticEDMFX
 import ..DiagnosticEDMFX
 import ..n_mass_flux_subdomains
 import ..gcm_driven_profile
-import ..gcm_driven_reference
+import ..gcm_height
+import ..gcm_driven_profile_tmean
 
 import Thermodynamics.TemperatureProfiles:
     DecayingTemperatureProfile, DryAdiabaticProfile
@@ -28,6 +29,7 @@ import AtmosphericProfilesLibrary as APL
 import SciMLBase
 import Interpolations as Intp
 import NCDatasets as NC
+import Statistics: mean
 
 include("local_state.jl")
 include("atmos_state.jl")

src/initial_conditions/initial_conditions.jl

@@ -1154,9 +1154,11 @@ function (initial_condition::GCMDriven)(params)
     thermo_params = CAP.thermodynamics_params(params)
 
     # Read forcing file
-    z_gcm = gcm_z(external_forcing_file)
+    z_gcm = NC.NCDataset(external_forcing_file) do ds
+        vec(gcm_height(ds.group["site23"]))
+    end
     vars = gcm_initial_conditions(external_forcing_file)
-    θ, u, v, q_tot, ρ₀ = map(vars) do value
+    T, u, v, q_tot, ρ₀ = map(vars) do value
         Intp.extrapolate(
             Intp.interpolate((z_gcm,), value, Intp.Gridded(Intp.Linear())),
             Intp.Flat(),
@@ -1169,10 +1171,10 @@ function (initial_condition::GCMDriven)(params)
         return LocalState(;
             params,
             geometry = local_geometry,
-            thermo_state = ts = TD.PhaseEquil_ρθq(
+            thermo_state = ts = TD.PhaseEquil_ρTq( 
                 thermo_params,
-                FT(ρ₀(z)),
-                FT(θ(z)),
+                FT(ρ₀(z)), 
+                FT(T(z)),
                 FT(q_tot(z)),
             ),
             velocity = Geometry.UVVector(FT(u(z)), FT(v(z))),
@@ -1182,22 +1184,15 @@ function (initial_condition::GCMDriven)(params)
     return local_state
 end
 
-# function gcm_z(external_forcing_file, FT::DataType)
-function gcm_z(external_forcing_file)
-    NC.NCDataset(external_forcing_file) do ds
-        gcm_driven_reference(ds, "z")[:]
-    end
-end
-
 # function gcm_initial_conditions(external_forcing_file, FT)
 function gcm_initial_conditions(external_forcing_file)
     NC.NCDataset(external_forcing_file) do ds
         (  # TODO: Cast to CuVector for GPU compatibility
-            gcm_driven_profile(ds, "thetali_mean")[:, 1], # 1 is initial time index
-            gcm_driven_profile(ds, "u_mean")[:, 1],
-            gcm_driven_profile(ds, "v_mean")[:, 1],
-            gcm_driven_profile(ds, "qt_mean")[:, 1],
-            gcm_driven_reference(ds, "rho0")[:],
+            gcm_driven_profile_tmean(ds.group["site23"], "ta"),
+            gcm_driven_profile_tmean(ds.group["site23"], "ua"),
+            gcm_driven_profile_tmean(ds.group["site23"], "va"),
+            gcm_driven_profile_tmean(ds.group["site23"], "hus"),
+            vec(mean( 1 ./ ds.group["site23"]["alpha"][:, :], dims = 2)), # convert alpha to rho using rho=1/alpha, take average profile
         )
     end
 end

src/prognostic_equations/forcing/external_forcing.jl

@@ -35,11 +35,11 @@ function compute_gcm_driven_scalar_inv_τ(z::FT) where {FT}
     end
 end
 
-external_forcing_cache(Y, atmos::AtmosModel) =
-    external_forcing_cache(Y, atmos.external_forcing)
+external_forcing_cache(Y, atmos::AtmosModel, params) =
+    external_forcing_cache(Y, atmos.external_forcing, params)
 
-external_forcing_cache(Y, external_forcing::Nothing) = (;)
-function external_forcing_cache(Y, external_forcing::GCMForcing)
+external_forcing_cache(Y, external_forcing::Nothing, params) = (;)
+function external_forcing_cache(Y, external_forcing::GCMForcing, params)
     FT = Spaces.undertype(axes(Y.c))
     ᶜdTdt_fluc = similar(Y.c, FT)
     ᶜdqtdt_fluc = similar(Y.c, FT)
@@ -56,39 +56,47 @@ function external_forcing_cache(Y, external_forcing::GCMForcing)
 
     (; external_forcing_file) = external_forcing
 
+    insolation, cos_zenith = [0f1], [0f1]
     NC.Dataset(external_forcing_file, "r") do ds
-        function setvar!(cc_field, varname, colidx, zc_gcm, zc_les)
+
+        function setvar!(cc_field, varname, colidx, zc_gcm, zc_forcing)
             parent(cc_field[colidx]) .= interp_vertical_prof(
                 zc_gcm,
-                zc_les,
-                gcm_driven_profile_tmean(ds, varname),  # TODO: time-varying tendencies
+                zc_forcing,
+                gcm_driven_profile_tmean(ds.group["site23"], varname), 
             )
         end
 
-        function setnudgevar!(cc_field, varname, colidx, zc_gcm, zc_les)
+        insolation[1] = mean(ds.group["site23"]["rsdt"][:] ./ ds.group["site23"]["coszen"][:]) #similar(Fields.level(Y.c.ρ, 1), FT) #similar(Spaces.level(Y.c.ρ, 1), FT)
+        cos_zenith[1] = ds.group["site23"]["coszen"][1]
+
+        function setvar_subsidence!(cc_field, varname, colidx, zc_gcm, zc_forcing, params)
             parent(cc_field[colidx]) .= interp_vertical_prof(
                 zc_gcm,
-                zc_les,
-                gcm_driven_profile(ds, varname)[:, 1],
+                zc_forcing,
+                gcm_driven_profile_tmean(ds.group["site23"] ,varname) .* 
+                            .-(gcm_driven_profile_tmean(ds.group["site23"] ,"alpha")) ./ 
+                            CAP.grav(params), 
             )
         end
 
-        zc_les = gcm_driven_reference(ds, "z")[:]
+
+        zc_forcing = gcm_height(ds.group["site23"])
         Fields.bycolumn(axes(Y.c)) do colidx
 
             zc_gcm = Fields.coordinate_field(Y.c).z[colidx]
 
-            setvar!(ᶜdTdt_fluc, "dtdt_fluc", colidx, zc_gcm, zc_les)
-            setvar!(ᶜdqtdt_fluc, "dqtdt_fluc", colidx, zc_gcm, zc_les)
-            setvar!(ᶜdTdt_hadv, "dtdt_hadv", colidx, zc_gcm, zc_les)
-            setvar!(ᶜdqtdt_hadv, "dqtdt_hadv", colidx, zc_gcm, zc_les)
-            setvar!(ᶜdTdt_rad, "dtdt_rad", colidx, zc_gcm, zc_les)
-            setvar!(ᶜls_subsidence, "ls_subsidence", colidx, zc_gcm, zc_les)
+            # setvar!(ᶜdTdt_fluc, "dtdt_fluc", colidx, zc_gcm, zc_forcing) #TODO: add these forcings back in
+            # setvar!(ᶜdqtdt_fluc, "dqtdt_fluc", colidx, zc_gcm, zc_forcing) #TODO: add these forcings back in
+            setvar!(ᶜdTdt_hadv, "tntha", colidx, zc_gcm, zc_forcing)
+            setvar!(ᶜdqtdt_hadv, "tnhusha", colidx, zc_gcm, zc_forcing)
+            setvar!(ᶜdTdt_rad, "tntr", colidx, zc_gcm, zc_forcing)
+            setvar_subsidence!(ᶜls_subsidence, "wap", colidx, zc_gcm, zc_forcing, params)
 
-            setnudgevar!(ᶜT_nudge, "temperature_mean", colidx, zc_gcm, zc_les)
-            setnudgevar!(ᶜqt_nudge, "qt_mean", colidx, zc_gcm, zc_les)
-            setnudgevar!(ᶜu_nudge, "u_mean", colidx, zc_gcm, zc_les)
-            setnudgevar!(ᶜv_nudge, "v_mean", colidx, zc_gcm, zc_les)
+            setvar!(ᶜT_nudge, "ta", colidx, zc_gcm, zc_forcing)
+            setvar!(ᶜqt_nudge, "hus", colidx, zc_gcm, zc_forcing)
+            setvar!(ᶜu_nudge, "ua", colidx, zc_gcm, zc_forcing)
+            setvar!(ᶜv_nudge, "va", colidx, zc_gcm, zc_forcing)
 
             @. ᶜinv_τ_wind[colidx] = 1 / (6 * 3600)
             @. ᶜinv_τ_scalar[colidx] = compute_gcm_driven_scalar_inv_τ(zc_gcm)
@@ -108,6 +116,8 @@ function external_forcing_cache(Y, external_forcing::GCMForcing)
         ᶜinv_τ_wind,
         ᶜinv_τ_scalar,
         ᶜls_subsidence,
+        insolation, 
+        cos_zenith,
     )
 end
 
@@ -145,8 +155,8 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::GCMForcing)
 
     ᶜdTdt_sum = p.scratch.ᶜtemp_scalar
     ᶜdqtdt_sum = p.scratch.ᶜtemp_scalar_2
-    @. ᶜdTdt_sum = ᶜdTdt_hadv + ᶜdTdt_fluc + ᶜdTdt_rad + ᶜdTdt_nudging
-    @. ᶜdqtdt_sum = ᶜdqtdt_hadv + ᶜdqtdt_fluc + ᶜdqtdt_nudging
+    @. ᶜdTdt_sum = ᶜdTdt_hadv + ᶜdTdt_nudging #  + ᶜdTdt_rad + ᶜdTdt_fluc remove nudging for now - TODO add back later
+    @. ᶜdqtdt_sum = ᶜdqtdt_hadv + ᶜdqtdt_nudging # + ᶜdqtdt_fluc remove nudging for now - TODO add back later
 
     T_0 = TD.Parameters.T_0(thermo_params)
     Lv_0 = TD.Parameters.LH_v0(thermo_params)

src/solver/model_getters.jl

@@ -49,13 +49,15 @@ end
 
 function get_insolation_form(parsed_args)
     insolation = parsed_args["insolation"]
-    @assert insolation in ("idealized", "timevarying", "rcemipii")
+    @assert insolation in ("idealized", "timevarying", "rcemipii", "gcmdriven")
     return if insolation == "idealized"
         IdealizedInsolation()
     elseif insolation == "timevarying"
         TimeVaryingInsolation()
     elseif insolation == "rcemipii"
         RCEMIPIIInsolation()
+    elseif insolation == "gcmdriven"
+        GCMDrivenInsolation()
     end
 end
 

src/solver/types.jl

@@ -32,6 +32,7 @@ abstract type AbstractInsolation end
 struct IdealizedInsolation <: AbstractInsolation end
 struct TimeVaryingInsolation <: AbstractInsolation end
 struct RCEMIPIIInsolation <: AbstractInsolation end
+struct GCMDrivenInsolation <: AbstractInsolation end
 
 abstract type AbstractSurfaceTemperature end
 struct PrescribedSurfaceTemperature <: AbstractSurfaceTemperature end

src/surface_conditions/SurfaceConditions.jl

@@ -9,7 +9,6 @@ import ..RCEMIPIISST
 import ..PrognosticSurfaceTemperature
 import ..PrescribedSurfaceTemperature
 import ..gcm_driven_timeseries
-import ..gcm_driven_reference
 
 import ..CT1, ..CT2, ..C12, ..CT12, ..C3
 import ..unit_basis_vector_data, ..projected_vector_data

src/surface_conditions/surface_setups.jl

@@ -263,12 +263,12 @@ function (surface_setup::GCMDriven)(params)
     return SurfaceState(; parameterization, T)
 end
 
-function gcm_surface_conditions(external_forcing_file; imin = 793)
+function gcm_surface_conditions(external_forcing_file)
     NC.NCDataset(external_forcing_file) do ds
         (
-            mean(gcm_driven_timeseries(ds, "surface_temperature")[imin:end]),
-            mean(gcm_driven_timeseries(ds, "lhf_surface_mean")[imin:end]),
-            mean(gcm_driven_timeseries(ds, "shf_surface_mean")[imin:end]),
+            mean(gcm_driven_timeseries(ds.group["site23"], "ts")),
+            mean(gcm_driven_timeseries(ds.group["site23"], "hfls")),
+            mean(gcm_driven_timeseries(ds.group["site23"], "hfss")),
         )
     end
 end