added rcemipii insolation and changed the 'idealized_insolation' argu…

…ment to 'insolation'

config/default_configs/default_config.yml

@@ -102,9 +102,9 @@ dt_rad:
 idealized_clouds:
   help: "Use idealized clouds in radiation model [`false` (default), `true`]"
   value: false
-idealized_insolation:
-  help: "Use idealized insolation in radiation model [`false`, `true` (default)]"
-  value: true
+insolation:
+  help: "Insolation used in radiation model [`idealized` (default), `timevarying`, `rcemipii`]"
+  value: "idealized"
 add_isothermal_boundary_layer:
   help: "Add an isothermal boundary layer above the domain top for radiation [`false`, `true` (default)]"
   value: true

config/gpu_configs/gpu_aquaplanet_dyamond_ss_1process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: equil
 precip_model: 1M
 rad: allskywithclear
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: 1hours
 vert_diff: FriersonDiffusion
 implicit_diffusion: true

config/gpu_configs/gpu_aquaplanet_dyamond_ss_2process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/gpu_configs/gpu_aquaplanet_dyamond_ss_4process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/gpu_configs/gpu_aquaplanet_dyamond_ws_1process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/gpu_configs/gpu_aquaplanet_dyamond_ws_2process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/gpu_configs/gpu_aquaplanet_dyamond_ws_4process.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/longrun_configs/amip_target_diagedmf.yml

@@ -12,7 +12,7 @@ override_τ_precip: false
 rad: "allskywithclear"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
-idealized_insolation: false
+insolation: "timevarying"
 surface_setup: "DefaultMoninObukhov"
 turbconv: "diagnostic_edmfx"
 ode_algo: ARS343

config/longrun_configs/longrun_aquaplanet_dyamond.yml

@@ -7,7 +7,7 @@ moist: "equil"
 precip_model: "0M" 
 cloud_model: "grid_scale"
 rad: "allskywithclear" 
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 vert_diff: "true" 
 surface_setup: "DefaultMoninObukhov" 

config/longrun_configs/longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_tvinsol_0M_slabocean.yml

@@ -14,7 +14,7 @@ vert_diff: "FriersonDiffusion"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 rad: "clearsky"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 albedo_model: "RegressionFunctionAlbedo"
 prognostic_surface: "PrognosticSurfaceTemperature"

config/model_configs/aquaplanet_diagedmf.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: equil 
 surface_setup: DefaultMoninObukhov 
 rad: allskywithclear
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours
 turbconv: diagnostic_edmfx 

config/model_configs/aquaplanet_progedmf.yml

@@ -9,7 +9,7 @@ dz_top: 3000.0
 moist: equil
 surface_setup: DefaultMoninObukhov
 rad: allskywithclear
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours
 turbconv: prognostic_edmfx

config/model_configs/aquaplanet_rhoe_equil_clearsky_tvinsol_0M_slabocean.yml

@@ -8,7 +8,7 @@ vert_diff: "FriersonDiffusion"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 rad: "clearsky"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 prognostic_surface: "PrognosticSurfaceTemperature"
 check_conservation: true

config/model_configs/gpu_aquaplanet_dyamond.yml

@@ -8,7 +8,7 @@ dz_top: 3000.0
 moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 dt_rad: "1hours"
 dt_cloud_fraction: "1hours"
 vert_diff: "FriersonDiffusion"

config/model_configs/rcemipii_sphere_diagnostic_edmfx.yml

@@ -4,6 +4,7 @@ rad: allskywithclear
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
+insolation: rcemipii
 prognostic_tke: true
 edmfx_upwinding: first_order
 edmfx_entr_model: "Generalized"

config/model_configs/single_column_radiative_equilibrium_allsky_idealized_clouds.yml

@@ -9,7 +9,7 @@ dz_top: 10000.0
 idealized_clouds: true
 dz_bottom: 100.0
 config: "column"
-idealized_insolation: false
+insolation: "timevarying"
 z_max: 70000.0
 dt_save_to_sol: "30hours"
 rad: "allskywithclear"

config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_raw_zonallyasymmetric.yml

@@ -9,7 +9,7 @@ t_end: "15hours"
 non_orographic_gravity_wave: true
 dz_bottom: 300.0
 vert_diff: "true"
-idealized_insolation: false
+insolation: "timevarying"
 z_max: 55000.0
 precip_model: "0M"
 cloud_model: "grid_scale"

config/model_configs/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.yml

@@ -11,7 +11,7 @@ moist: "equil"
 precip_model: "1M"
 rad: "allskywithclear"
 aerosol_radiation: true
-idealized_insolation: false
+insolation: "timevarying"
 rayleigh_sponge: true
 non_orographic_gravity_wave: true
 orographic_gravity_wave: "gfdl_restart"

config/model_configs/sphere_aquaplanet_rhoe_nonequilmoist_allsky.yml

@@ -10,7 +10,7 @@ approximate_linear_solve_iters: 2
 moist: "nonequil"
 precip_model: "nothing"
 rad: "allskywithclear"
-idealized_insolation: false
+insolation: "timevarying"
 rayleigh_sponge: true
 surface_setup: "DefaultMoninObukhov"
 toml: [toml/sphere_aquaplanet_rhoe_equilmoist_allsky_gw_res.toml]

config/perf_configs/flame_perf_gw.yml

@@ -18,6 +18,6 @@ rayleigh_sponge: true
 toml:
   - "toml/flame_perf_gw.toml"
 orographic_gravity_wave: "raw_topo"
-idealized_insolation: false
+insolation: "timevarying"
 non_orographic_gravity_wave: true
 z_max: 55000.0

src/callbacks/callbacks.jl

@@ -59,7 +59,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
 
     (; ᶜts, cloud_diagnostics_tuple, sfc_conditions) = p.precomputed
     (; params) = p
-    (; idealized_insolation, idealized_h2o, idealized_clouds) = p.radiation
+    (; idealized_h2o, idealized_clouds) = p.radiation
     (; ᶠradiation_flux, radiation_model) = p.radiation
     (; radiation_mode) = p.atmos
 
@@ -128,9 +128,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
         end
     end
 
-    if !idealized_insolation && !(p.atmos.surface_albedo isa CouplerAlbedo)
-        set_insolation_variables!(Y, p, t)
-    end
+    set_insolation_variables!(Y, p, t, p.atmos.insolation)
 
     if !idealized_clouds && !(
         radiation_mode isa RRTMGPI.GrayRadiation ||
@@ -176,50 +174,81 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
     return nothing
 end
 
-function set_insolation_variables!(Y, p, t)
+#Uniform insolation, magnitudes from Wing et al. (2018)
+#Summary will show incorrect insolation magnitude, but the magnitude is correct above the isothermal boundary layer
+function set_insolation_variables!(Y, p, t, ::RCEMIPIIInsolation)
+    FT = Spaces.undertype(axes(Y.c))
+    (; radiation_model) = p.radiation
+    radiation_model.cos_zenith .= cosd(FT(42.05))
+    radiation_model.weighted_irradiance .= FT(551.58)
+end
 
+function set_insolation_variables!(Y, p, t, ::IdealizedInsolation)
     FT = Spaces.undertype(axes(Y.c))
-    params = p.params
-    insolation_params = CAP.insolation_params(params)
-    (; insolation_tuple, radiation_model) = p.radiation
-
-    current_datetime = p.start_date + Dates.Second(round(Int, t)) # current time
-    max_zenith_angle = FT(π) / 2 - eps(FT)
-    irradiance = FT(CAP.tot_solar_irrad(params))
-    au = FT(CAP.astro_unit(params))
-    date0 = DateTime("2000-01-01T11:58:56.816")
-    d, δ, η_UTC =
-        FT.(
-            Insolation.helper_instantaneous_zenith_angle(
-                current_datetime,
-                date0,
-                insolation_params,
-            )
-        )
     bottom_coords = Fields.coordinate_field(Spaces.level(Y.c, 1))
     if eltype(bottom_coords) <: Geometry.LatLongZPoint
-        cos_zenith =
-            Fields.array2field(radiation_model.cos_zenith, axes(bottom_coords))
-        weighted_irradiance = Fields.array2field(
-            radiation_model.weighted_irradiance,
-            axes(bottom_coords),
-        )
-        @. insolation_tuple = instantaneous_zenith_angle(
-            d,
-            δ,
-            η_UTC,
-            bottom_coords.long,
-            bottom_coords.lat,
-        ) # the tuple is (zenith angle, azimuthal angle, earth-sun distance)
-        @. cos_zenith = cos(min(first(insolation_tuple), max_zenith_angle))
-        @. weighted_irradiance = irradiance * (au / last(insolation_tuple))^2
+        latitude = Fields.field2array(bottom_coords.lat)
     else
-        # assume that the latitude and longitude are both 0 for flat space
-        insolation_tuple = instantaneous_zenith_angle(d, δ, η_UTC, FT(0), FT(0))
-        radiation_model.cos_zenith .=
-            cos(min(first(insolation_tuple), max_zenith_angle))
-        radiation_model.weighted_irradiance .=
-            irradiance * (au / last(insolation_tuple))^2
+        latitude = Fields.field2array(zero(bottom_coords.z)) # flat space is on Equator
+    end
+    (; radiation_model) = p.radiation
+    # perpetual equinox with no diurnal cycle
+    radiation_model.cos_zenith .= cos(FT(π) / 3)
+    weighted_irradiance =
+        @. 1360 * (1 + FT(1.2) / 4 * (1 - 3 * sind(latitude)^2)) /
+           (4 * cos(FT(π) / 3))
+    radiation_model.weighted_irradiance .= weighted_irradiance
+end
+
+function set_insolation_variables!(Y, p, t, ::TimeVaryingInsolation)
+    if !(p.atmos.surface_albedo isa CouplerAlbedo)
+        FT = Spaces.undertype(axes(Y.c))
+        params = p.params
+        insolation_params = CAP.insolation_params(params)
+        (; insolation_tuple, radiation_model) = p.radiation
+
+        current_datetime = p.start_date + Dates.Second(round(Int, t)) # current time
+        max_zenith_angle = FT(π) / 2 - eps(FT)
+        irradiance = FT(CAP.tot_solar_irrad(params))
+        au = FT(CAP.astro_unit(params))
+        date0 = DateTime("2000-01-01T11:58:56.816")
+        d, δ, η_UTC =
+            FT.(
+                Insolation.helper_instantaneous_zenith_angle(
+                    current_datetime,
+                    date0,
+                    insolation_params,
+                )
+            )
+        bottom_coords = Fields.coordinate_field(Spaces.level(Y.c, 1))
+        if eltype(bottom_coords) <: Geometry.LatLongZPoint
+            cos_zenith = Fields.array2field(
+                radiation_model.cos_zenith,
+                axes(bottom_coords),
+            )
+            weighted_irradiance = Fields.array2field(
+                radiation_model.weighted_irradiance,
+                axes(bottom_coords),
+            )
+            @. insolation_tuple = instantaneous_zenith_angle(
+                d,
+                δ,
+                η_UTC,
+                bottom_coords.long,
+                bottom_coords.lat,
+            ) # the tuple is (zenith angle, azimuthal angle, earth-sun distance)
+            @. cos_zenith = cos(min(first(insolation_tuple), max_zenith_angle))
+            @. weighted_irradiance =
+                irradiance * (au / last(insolation_tuple))^2
+        else
+            # assume that the latitude and longitude are both 0 for flat space
+            insolation_tuple =
+                instantaneous_zenith_angle(d, δ, η_UTC, FT(0), FT(0))
+            radiation_model.cos_zenith .=
+                cos(min(first(insolation_tuple), max_zenith_angle))
+            radiation_model.weighted_irradiance .=
+                irradiance * (au / last(insolation_tuple))^2
+        end
     end
 end
 

src/parameterized_tendencies/radiation/RRTMGPInterface.jl

@@ -14,27 +14,23 @@ using NVTX
 abstract type AbstractRRTMGPMode end
 struct GrayRadiation <: AbstractRRTMGPMode
     idealized_h2o::Bool
-    idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
 end
 struct ClearSkyRadiation <: AbstractRRTMGPMode
     idealized_h2o::Bool
-    idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool
 end
 struct AllSkyRadiation <: AbstractRRTMGPMode
     idealized_h2o::Bool
-    idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool
 end
 struct AllSkyRadiationWithClearSkyDiagnostics <: AbstractRRTMGPMode
     idealized_h2o::Bool
-    idealized_insolation::Bool
     idealized_clouds::Bool
     add_isothermal_boundary_layer::Bool
     aerosol_radiation::Bool

src/parameterized_tendencies/radiation/radiation.jl

@@ -40,7 +40,7 @@ function radiation_model_cache(
 )
     context = ClimaComms.context(axes(Y.c))
     device = context.device
-    (; idealized_h2o, idealized_insolation, idealized_clouds) = radiation_mode
+    (; idealized_h2o, idealized_clouds) = radiation_mode
     if !(radiation_mode isa RRTMGPI.GrayRadiation)
         (; aerosol_radiation) = radiation_mode
         if aerosol_radiation && !("SO4" in aerosol_names)
@@ -207,14 +207,7 @@ function radiation_model_cache(
             )
         end
 
-        if idealized_insolation # perpetual equinox with no diurnal cycle
-            cos_zenith = cos(FT(π) / 3)
-            weighted_irradiance =
-                @. 1360 * (1 + FT(1.2) / 4 * (1 - 3 * sind(latitude)^2)) /
-                   (4 * cos_zenith)
-        else
-            cos_zenith = weighted_irradiance = NaN # initialized in callback
-        end
+        cos_zenith = weighted_irradiance = NaN # initialized in callback
 
         radiation_model = RRTMGPI.RRTMGPModel(
             rrtmgp_params,
@@ -237,7 +230,6 @@ function radiation_model_cache(
         )
     end
     return (;
-        idealized_insolation,
         orbital_data,
         idealized_h2o,
         idealized_clouds,