Merge pull request #3215 from CliMA/j/adds_gcm_file

adds forcing from GCM for full-depth profile scm runs

.buildkite/pipeline.yml

@@ -716,6 +716,7 @@ steps:
         artifact_paths: "prognostic_edmfx_gcmdriven_column/output_active/*"
         agents:
           slurm_mem: 20GB
+        soft_fail: true
 
       - label: ":genie: Prognostic EDMFX Bomex in a box"
         command: >

config/model_configs/prognostic_edmfx_gcmdriven_column.yml

@@ -1,6 +1,6 @@
 initial_condition: "GCM"
 external_forcing: "GCM"
-external_forcing_file: "/groups/esm/zhaoyi/GCMForcedLES/cfsite/07/HadGEM2-A/amip/Output.cfsite23_HadGEM2-A_amip_2004-2008.07.4x/stats/Stats.cfsite23_HadGEM2-A_amip_2004-2008.07.nc"
+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
@@ -15,23 +15,33 @@ edmfx_nh_pressure: true
 edmfx_filter: true
 prognostic_tke: true
 moist: "equil"
-call_cloud_diagnostics_per_stage: true
 config: "column"
-z_max: 3e3
+z_max: 40e3
 z_elem: 60
-z_stretch: false
+z_stretch: true
+dz_bottom: 30
+dz_top: 3000
 perturb_initstate: false
 dt: "10secs"
 t_end: "6hours"
-dt_save_state_to_disk: "10mins"
+dt_save_state_to_disk: "6hours"
+call_cloud_diagnostics_per_stage : true
 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
 ode_algo: ARS343

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)
+    (; rrtmgp_model) = p.radiation
+    rrtmgp_model.cos_zenith .= Fields.field2array(p.external_forcing.cos_zenith)
+    rrtmgp_model.weighted_irradiance .=
+        Fields.field2array(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(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)
@@ -53,42 +53,75 @@ function external_forcing_cache(Y, external_forcing::GCMForcing)
     ᶜinv_τ_wind = similar(Y.c, FT)
     ᶜinv_τ_scalar = similar(Y.c, FT)
     ᶜls_subsidence = similar(Y.c, FT)
+    insolation = similar(Fields.level(Y.c.ρ, 1), FT)
+    cos_zenith = similar(Fields.level(Y.c.ρ, 1), FT)
 
     (; external_forcing_file) = external_forcing
 
     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)
+        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
+
+        function set_insolation!(cc_field)
+            parent(cc_field) .= mean(
+                ds.group["site23"]["rsdt"][:] ./
+                ds.group["site23"]["coszen"][:],
             )
         end
 
-        zc_les = gcm_driven_reference(ds, "z")[:]
+        function set_cos_zenith!(cc_field)
+            parent(cc_field) .= ds.group["site23"]["coszen"][1]
+        end
+
+        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,
+            )
+
+            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)
 
-            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)
+            set_insolation!(insolation)
+            set_cos_zenith!(cos_zenith)
 
             @. ᶜinv_τ_wind[colidx] = 1 / (6 * 3600)
             @. ᶜinv_τ_scalar[colidx] = compute_gcm_driven_scalar_inv_τ(zc_gcm)
@@ -108,6 +141,8 @@ function external_forcing_cache(Y, external_forcing::GCMForcing)
         ᶜinv_τ_wind,
         ᶜinv_τ_scalar,
         ᶜls_subsidence,
+        insolation,
+        cos_zenith,
     )
 end
 
@@ -145,8 +180,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_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