From a17f8b71e8f6503187e75d009f673dbfa4a5c98a Mon Sep 17 00:00:00 2001 From: Julian Schmitt Date: Tue, 23 Jul 2024 14:47:40 -0700 Subject: [PATCH] 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 --- .../prognostic_edmfx_gcmdriven_fullcolumn.yml | 48 ++++++++++++++++ src/cache/cache.jl | 2 +- src/callbacks/callbacks.jl | 7 +++ src/initial_conditions/InitialConditions.jl | 4 +- src/initial_conditions/initial_conditions.jl | 29 ++++------ .../forcing/external_forcing.jl | 56 ++++++++++-------- src/solver/model_getters.jl | 4 +- src/solver/types.jl | 1 + src/surface_conditions/SurfaceConditions.jl | 1 - src/surface_conditions/surface_setups.jl | 8 +-- src/utils/read_gcm_driven_scm_data.jl | 57 +++++-------------- 11 files changed, 126 insertions(+), 91 deletions(-) create mode 100644 config/model_configs/prognostic_edmfx_gcmdriven_fullcolumn.yml diff --git a/config/model_configs/prognostic_edmfx_gcmdriven_fullcolumn.yml b/config/model_configs/prognostic_edmfx_gcmdriven_fullcolumn.yml new file mode 100644 index 00000000000..f6b977e2b97 --- /dev/null +++ b/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 diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 56b7a2bb498..6c1feb4e016 100644 --- a/src/cache/cache.jl +++ b/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) diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 622ee4e1d2c..a381e45b061 100644 --- a/src/callbacks/callbacks.jl +++ b/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)) diff --git a/src/initial_conditions/InitialConditions.jl b/src/initial_conditions/InitialConditions.jl index f88fcfa5a31..73e0e581518 100644 --- a/src/initial_conditions/InitialConditions.jl +++ b/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") diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index 97915ca3047..4c26362e838 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/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 diff --git a/src/prognostic_equations/forcing/external_forcing.jl b/src/prognostic_equations/forcing/external_forcing.jl index 822d41b4ae1..c4fc2560337 100644 --- a/src/prognostic_equations/forcing/external_forcing.jl +++ b/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) diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 8749eb86bbe..83352050b0e 100644 --- a/src/solver/model_getters.jl +++ b/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 diff --git a/src/solver/types.jl b/src/solver/types.jl index e4f36e58611..f26a0a2c2e6 100644 --- a/src/solver/types.jl +++ b/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 diff --git a/src/surface_conditions/SurfaceConditions.jl b/src/surface_conditions/SurfaceConditions.jl index b8b6ba75a5e..c14565f8560 100644 --- a/src/surface_conditions/SurfaceConditions.jl +++ b/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 diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl index 79c0c2e0008..68b3921f81b 100644 --- a/src/surface_conditions/surface_setups.jl +++ b/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 diff --git a/src/utils/read_gcm_driven_scm_data.jl b/src/utils/read_gcm_driven_scm_data.jl index e65f5ee5887..a9e4d89add9 100644 --- a/src/utils/read_gcm_driven_scm_data.jl +++ b/src/utils/read_gcm_driven_scm_data.jl @@ -3,9 +3,9 @@ import NCDatasets as NC import StatsBase: mean """ - gcm_driven_profile_tmean(ds, varname; imin = 793) + gcm_driven_profile_tmean(ds, varname) -Extract time-averaged data (`imin:end`) for `varname` from the "profile" group in the GCM-driven dataset `ds` +Extract time-averaged data for `varname` from the "profile" group in the GCM-driven dataset `ds` Returns a 1D ("z",) `Vector{FT}` of the time-averaged data. @@ -14,8 +14,8 @@ Returns a 1D ("z",) `Vector{FT}` of the time-averaged data. This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -function gcm_driven_profile_tmean(ds, varname; imin = 793) - vec(mean(gcm_driven_profile(ds, varname)[:, imin:end]; dims = 2)) +function gcm_driven_profile_tmean(ds, varname) + vec(mean(gcm_driven_profile(ds, varname), dims = 2)) end """ @@ -23,7 +23,7 @@ end Extract `varname` from the "profile" group in the GCM-driven dataset `ds` -Returns a 2D ("z", "t") `NCDatasets.Variable` object. +Returns a 2D ("z", "t") `Matrix` object. !!! note @@ -31,63 +31,34 @@ Returns a 2D ("z", "t") `NCDatasets.Variable` object. If this is not the case, "garbage" data may be returned with no warning. """ function gcm_driven_profile(ds, varname) - varname ∈ ("z", "z_half", "t") && - throw(ArgumentError("This method does not support access to $varname")) - _gcm_driven_variable(ds.group["profiles"], varname, ("z", "t")) + ds[varname][:,:] end """ gcm_driven_reference(ds, varname) -Extract `varname` from the "reference" group in the GCM-driven dataset `ds` +Extract `height` from the GCM-driven dataset `ds` -Returns a 1D ("z",) `NCDatasets.Variable` object. +Returns a 1D ("z",) `vec` object. !!! note This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -function gcm_driven_reference(ds, varname) - dimnames = endswith(varname, "_full") ? ("z_full",) : ("z",) - _gcm_driven_variable(ds.group["reference"], varname, dimnames) +function gcm_height(ds) + vec(mean(ds["zg"][:,:], dims=2)) end - """ - gcm_driven_timeseries(ds, varname) - -Extract `varname` from the "timeseries" group in the GCM-driven dataset `ds` +# gcm_driven_timeseries(ds, varname) -Returns a 1D ("t",) `NCDatasets.Variable` object. +# Get `varname` from the dataset `ds` and return values !!! note This method currently assumes the underlying data is `Float64`. If this is not the case, "garbage" data may be returned with no warning. """ -gcm_driven_timeseries(ds, varname) = - _gcm_driven_variable(ds.group["timeseries"], varname, ("t",)) - -""" - _gcm_driven_variable(ds, varname, dimnames, [FT=Float64]) - -Fetch a variable (`varname`) from a GCM-driven SCM NetCDF dataset (or subgroup) `ds`. - -Returns an `NCDatasets.Variable` object. - -!!! note - - It is critical to provide _correct_ data type `FT` (e.g. `Float64`) and dimension names `dimnames` - associated with the variable in the dataset, - otherwise "garbage" data is returned with no warning (wrong `FT`), or an error is thrown (wrong `dimnames`). -""" -function _gcm_driven_variable( - ds, - varname, - dimnames::NTuple{N, String}, - FT = Float64, -) where {N} - dimids = NC.nc_inq_dimid.(ds.ncid, dimnames) - varid = NC.nc_inq_varid(ds.ncid, varname) - NC.Variable{FT, N, typeof(ds)}(ds, varid, dimids) +function gcm_driven_timeseries(ds, varname) + ds[varname][:] end