diff --git a/integration_tests/utils/Cases.jl b/integration_tests/utils/Cases.jl
index 7561034e1..20b89a5a2 100644
--- a/integration_tests/utils/Cases.jl
+++ b/integration_tests/utils/Cases.jl
@@ -133,7 +133,7 @@ get_radiation_type(::LES_driven_SCM) = TC.RadiationLES
 ##### Default CasesBase behavior:
 #####
 
-initialize_radiation(self::CasesBase, grid, state, gm, param_set) = initialize(self.Rad, grid)
+initialize_radiation(self::CasesBase, grid, state, gm, param_set) = initialize(self.Rad, grid, state)
 
 function TC.initialize_io(self::CasesBase, Stats::NetCDFIO_Stats, ::BaseCase)
     add_ts(Stats, "Tsurface")
@@ -155,7 +155,7 @@ TC.update_forcing(self::CasesBase, grid, state, gm, TS::TimeStepping, param_set)
 TC.update_radiation(self::CasesBase, grid, state, gm, TS::TimeStepping, param_set) =
     update(self.Rad, grid, state, gm, param_set)
 
-initialize_forcing(self::CasesBase, grid::Grid, state, gm, param_set) = initialize(self.Fo, grid)
+initialize_forcing(self::CasesBase, grid::Grid, state, gm, param_set) = initialize(self.Fo, grid, state)
 
 #####
 ##### Soares
@@ -368,13 +368,14 @@ function initialize_surface(self::CasesBase{Bomex}, grid::Grid, state, param_set
 end
 
 function initialize_forcing(self::CasesBase{Bomex}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    initialize(self.Fo, grid, state)
     p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         # Geostrophic velocity profiles. vg = 0
-        self.Fo.ug[k] = -10.0 + (1.8e-3) * z
+        aux_gm.ug[k] = -10.0 + (1.8e-3) * z
         ts = TC.thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
         Π = TD.exner(ts)
         # Set large-scale cooling
@@ -383,7 +384,7 @@ function initialize_forcing(self::CasesBase{Bomex}, grid::Grid, state, gm, param
         else
             (-2.0 / (3600 * 24.0) + (z - 1500.0) * (0.0 - -2.0 / (3600 * 24.0)) / (3000.0 - 1500.0)) * Π
         end
-        self.Fo.dTdt[k] = dTdt
+        aux_gm.dTdt[k] = dTdt
 
         # Set large-scale drying
         dqtdt = if z <= 300.0
@@ -393,7 +394,7 @@ function initialize_forcing(self::CasesBase{Bomex}, grid::Grid, state, gm, param
         else
             0
         end
-        self.Fo.dqtdt[k] = dqtdt
+        aux_gm.dqtdt[k] = dqtdt
 
         #Set large scale subsidence
         subsidence = if z <= 1500.0
@@ -403,7 +404,7 @@ function initialize_forcing(self::CasesBase{Bomex}, grid::Grid, state, gm, param
         else
             0
         end
-        self.Fo.subsidence[k] = subsidence
+        aux_gm.subsidence[k] = subsidence
     end
     return nothing
 end
@@ -499,24 +500,25 @@ function initialize_surface(self::CasesBase{life_cycle_Tan2018}, grid::Grid, sta
     self.Sur.bflux = life_cycle_buoyancy_flux(param_set)
 end
 function initialize_forcing(self::CasesBase{life_cycle_Tan2018}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    initialize(self.Fo, grid, state)
     p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         # Geostrophic velocity profiles. vg = 0
-        self.Fo.ug[k] = -10.0 + (1.8e-3) * z
+        aux_gm.ug[k] = -10.0 + (1.8e-3) * z
         ts = TC.thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
         Π = TD.exner(ts)
         # Set large-scale cooling
-        self.Fo.dTdt[k] = if z <= 1500.0
+        aux_gm.dTdt[k] = if z <= 1500.0
             (-2.0 / (3600 * 24.0)) * Π
         else
             (-2.0 / (3600 * 24.0) + (z - 1500.0) * (0.0 - -2.0 / (3600 * 24.0)) / (3000.0 - 1500.0)) * Π
         end
 
         # Set large-scale drying
-        self.Fo.dqtdt[k] = if z <= 300.0
+        aux_gm.dqtdt[k] = if z <= 300.0
             -1.2e-8   #kg/(kg * s)
         elseif z > 300.0 && z <= 500.0
             -1.2e-8 + (z - 300.0) * (0.0 - -1.2e-8) / (500.0 - 300.0) #kg/(kg * s)
@@ -525,7 +527,7 @@ function initialize_forcing(self::CasesBase{life_cycle_Tan2018}, grid::Grid, sta
         end
 
         #Set large scale subsidence
-        self.Fo.subsidence[k] = if z <= 1500.0
+        aux_gm.subsidence[k] = if z <= 1500.0
             0.0 + z * (-0.65 / 100.0 - 0.0) / (1500.0 - 0.0)
         elseif z > 1500.0 && z <= 2100.0
             -0.65 / 100 + (z - 1500.0) * (0.0 - -0.65 / 100.0) / (2100.0 - 1500.0)
@@ -632,28 +634,29 @@ function initialize_surface(self::CasesBase{Rico}, grid::Grid, state, param_set)
 end
 
 function initialize_forcing(self::CasesBase{Rico}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    initialize(self.Fo, grid, state)
     p0_c = TC.center_ref_state(state).p0
     prog_gm = TC.center_prog_grid_mean(state)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
         z = grid.zc[k]
         ts = TC.thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
         Π = TD.exner(ts)
         # Geostrophic velocity profiles
-        self.Fo.ug[k] = -9.9 + 2.0e-3 * z
-        self.Fo.vg[k] = -3.8
+        aux_gm.ug[k] = -9.9 + 2.0e-3 * z
+        aux_gm.vg[k] = -3.8
         # Set large-scale cooling
-        self.Fo.dTdt[k] = (-2.5 / (3600.0 * 24.0)) * Π
+        aux_gm.dTdt[k] = (-2.5 / (3600.0 * 24.0)) * Π
 
         # Set large-scale moistening
-        self.Fo.dqtdt[k] = if z <= 2980.0
+        aux_gm.dqtdt[k] = if z <= 2980.0
             (-1.0 + 1.3456 / 2980.0 * z) / 86400.0 / 1000.0   #kg/(kg * s)
         else
             0.3456 / 86400.0 / 1000.0
         end
 
         #Set large scale subsidence
-        self.Fo.subsidence[k] = if z <= 2260.0
+        aux_gm.subsidence[k] = if z <= 2260.0
             -(0.005 / 2260.0) * z
         else
             -0.005
@@ -739,7 +742,6 @@ function initialize_profiles(self::CasesBase{TRMM_LBA}, grid::Grid, gm, state)
 
     zc_in = grid.zc
     molmass_ratio = CPP.molmass_ratio(param_set)
-    aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
 
     parent(prog_gm.u) .= pyinterp(zc_in, z_in, u_in)
@@ -779,8 +781,7 @@ function initialize_surface(self::CasesBase{TRMM_LBA}, grid::Grid, state, param_
 end
 
 function initialize_forcing(self::CasesBase{TRMM_LBA}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
-    self.Fo.dTdt = TC.center_field(grid)
+
     self.rad_time = range(10, 360; length = 36) .* 60
     #! format: off
     z_in         = arr_type([42.5, 200.92, 456.28, 743, 1061.08, 1410.52, 1791.32, 2203.48, 2647,3121.88, 3628.12,
@@ -911,8 +912,9 @@ function initialize_forcing(self::CasesBase{TRMM_LBA}, grid::Grid, state, gm, pa
     self.rad .= A # store matrix in self
     ind1 = Int(trunc(10.0 / 600.0)) + 1
     ind2 = Int(ceil(10.0 / 600.0))
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
-        self.Fo.dTdt[k] = if 10 % 600.0 == 0
+        aux_gm.dTdt[k] = if 10 % 600.0 == 0
             self.rad[ind1, k]
         else
             (self.rad[ind2, k] - self.rad[ind1, k]) / (self.rad_time[ind2] - self.rad_time[ind1]) * (10.0) +
@@ -933,12 +935,13 @@ end
 
 function TC.update_forcing(self::CasesBase{TRMM_LBA}, grid, state, gm, TS::TimeStepping, param_set)
 
+    aux_gm = TC.center_aux_grid_mean(state)
     ind2 = Int(ceil(TS.t / 600.0)) + 1
     ind1 = Int(trunc(TS.t / 600.0)) + 1
     rad_time = self.rad_time
     rad = self.rad
     @inbounds for k in real_center_indices(grid)
-        self.Fo.dTdt[k] = if grid.zc[k] >= 22699.48
+        aux_gm.dTdt[k] = if grid.zc[k] >= 22699.48
             0.0
         else
             if TS.t < 600.0 # first 10 min use the radiative forcing of t=10min (as in the paper)
@@ -980,12 +983,10 @@ end
 function initialize_profiles(self::CasesBase{ARM_SGP}, grid::Grid, gm, state)
     # ARM_SGP inputs
     p0 = TC.center_ref_state(state).p0
-    aux_gm = TC.center_aux_grid_mean(state)
-    prog_gm = TC.center_prog_grid_mean(state)
-    #! format: off
     param_set = TC.parameter_set(gm)
-    aux_gm = TC.center_aux_grid_mean(state)
     prog_gm = TC.center_prog_grid_mean(state)
+    aux_gm = TC.center_aux_grid_mean(state)
+    #! format: off
     z_in = arr_type([0.0, 50.0, 350.0, 650.0, 700.0, 1300.0, 2500.0, 5500.0 ]) #LES z is in meters
     Theta_in = arr_type([299.0, 301.5, 302.5, 303.53, 303.7, 307.13, 314.0, 343.2]) # K
     r_in = arr_type([15.2,15.17,14.98,14.8,14.7,13.5,3.0,3.0])/1000 # qt should be in kg/kg
@@ -1030,10 +1031,10 @@ function initialize_surface(self::CasesBase{ARM_SGP}, grid::Grid, state, param_s
 end
 
 function initialize_forcing(self::CasesBase{ARM_SGP}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
-        self.Fo.ug[k] = 10.0
-        self.Fo.vg[k] = 0.0
+        aux_gm.ug[k] = 10.0
+        aux_gm.vg[k] = 0.0
     end
     return nothing
 end
@@ -1058,6 +1059,7 @@ function TC.update_surface(self::CasesBase{ARM_SGP}, grid, state, gm, TS::TimeSt
 end
 
 function TC.update_forcing(self::CasesBase{ARM_SGP}, grid, state, gm, TS::TimeStepping, param_set)
+    aux_gm = TC.center_aux_grid_mean(state)
     p0_c = TC.center_ref_state(state).p0
     t_in = arr_type([0.0, 3.0, 6.0, 9.0, 12.0, 14.5]) .* 3600.0 #LES time is in sec
     AT_in = arr_type([0.0, 0.0, 0.0, -0.08, -0.016, -0.016]) ./ 3600.0 # Advective forcing for theta [K/h] converted to [K/sec]
@@ -1070,14 +1072,14 @@ function TC.update_forcing(self::CasesBase{ARM_SGP}, grid, state, gm, TS::TimeSt
         ts = TC.thermo_state_pθq(param_set, p0_c[k], prog_gm.θ_liq_ice[k], prog_gm.q_tot[k])
         Π = TD.exner(ts)
         z = grid.zc[k]
-        self.Fo.dTdt[k] = if z <= 1000.0
+        aux_gm.dTdt[k] = if z <= 1000.0
             dTdt
         elseif z > 1000.0 && z <= 2000.0
             dTdt * (1 - (z - 1000.0) / 1000.0)
         else
             0
         end
-        self.Fo.dqtdt[k] = if z <= 1000.0
+        aux_gm.dqtdt[k] = if z <= 1000.0
             dqtdt * Π
         elseif z > 1000.0 && z <= 2000.0
             dqtdt * Π * (1 - (z - 1000.0) / 1000.0)
@@ -1169,7 +1171,8 @@ function initialize_surface(self::CasesBase{GATE_III}, grid::Grid, state, param_
 end
 
 function initialize_forcing(self::CasesBase{GATE_III}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+
+    aux_gm = TC.center_aux_grid_mean(state)
     #LES z is in meters
     #! format: off
     z_in     = arr_type([ 0.0,   0.5,  1.0,  1.5,   2.0,   2.5,    3.0,   3.5,   4.0,   4.5,   5.0,   5.5,   6.0,
@@ -1194,8 +1197,12 @@ function initialize_forcing(self::CasesBase{GATE_III}, grid::Grid, state, gm, pa
 
     Qtend_in = r_tend_in ./ (1 .+ r_tend_in) # convert mixing ratio to specific humidity
 
-    self.Fo.dqtdt = pyinterp(grid.zc, z_in, Qtend_in)
-    self.Fo.dTdt = pyinterp(grid.zc, z_in, Ttend_in) + pyinterp(grid.zc, z_in, RAD_in)
+    dqtdt = pyinterp(grid.zc, z_in, Qtend_in)
+    dTdt = pyinterp(grid.zc, z_in, Ttend_in) + pyinterp(grid.zc, z_in, RAD_in)
+    for k in TC.real_center_indices(grid)
+        aux_gm.dqtdt[k] = dqtdt[k]
+        aux_gm.dTdt[k] = dTdt[k]
+    end
 end
 
 #####
@@ -1286,29 +1293,30 @@ function initialize_surface(self::CasesBase{DYCOMS_RF01}, grid::Grid, state, par
         )
 end
 function initialize_forcing(self::CasesBase{DYCOMS_RF01}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    aux_gm = TC.center_aux_grid_mean(state)
 
     # geostrophic velocity profiles
-    self.Fo.ug .= 7.0
-    self.Fo.vg .= -5.5
+    parent(aux_gm.ug) .= 7.0
+    parent(aux_gm.vg) .= -5.5
 
     # large scale subsidence
     divergence = 3.75e-6    # divergence is defined twice: here and in __init__ of ForcingDYCOMS_RF01 class
     # To be able to have self.Fo.divergence available here,
     # we would have to change the signature of ForcingBase class
     @inbounds for k in real_center_indices(grid)
-        self.Fo.subsidence[k] = -grid.zc[k] * divergence
+        aux_gm.subsidence[k] = -grid.zc[k] * divergence
     end
 
     # no large-scale drying
-    self.Fo.dqtdt .= 0.0 #kg/(kg * s)
+    parent(aux_gm.dqtdt) .= 0 #kg/(kg * s)
 end
 
 function initialize_radiation(self::CasesBase{DYCOMS_RF01}, grid::Grid, state, gm, param_set)
-    initialize(self.Rad, grid)
+    initialize(self.Rad, grid, state)
+    aux_gm = TC.center_aux_grid_mean(state)
 
     # no large-scale drying
-    self.Rad.dqtdt .= 0.0 #kg/(kg * s)
+    parent(aux_gm.dqtdt_rad) .= 0 #kg/(kg * s)
 
     # Radiation based on eq. 3 in Stevens et. al., (2005)
     # cloud-top cooling + cloud-base warming + cooling in free troposphere
@@ -1376,11 +1384,12 @@ function initialize_surface(self::CasesBase{GABLS}, grid::Grid, state, param_set
 end
 
 function initialize_forcing(self::CasesBase{GABLS}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    initialize(self.Fo, grid, state)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
         # Geostrophic velocity profiles.
-        self.Fo.ug[k] = 8.0
-        self.Fo.vg[k] = 0.0
+        aux_gm.ug[k] = 8.0
+        aux_gm.vg[k] = 0.0
     end
     return nothing
 end
@@ -1455,11 +1464,11 @@ function initialize_surface(self::CasesBase{SP}, grid::Grid, state, param_set)
 end
 
 function initialize_forcing(self::CasesBase{SP}, grid::Grid, state, gm, param_set)
-    initialize(self.Fo, grid)
+    aux_gm = TC.center_aux_grid_mean(state)
     @inbounds for k in real_center_indices(grid)
         # Geostrophic velocity profiles. vg = 0
-        self.Fo.ug[k] = 1.0
-        self.Fo.vg[k] = 0.0
+        aux_gm.ug[k] = 1.0
+        aux_gm.vg[k] = 0.0
     end
 end
 
@@ -1653,9 +1662,9 @@ function initialize_surface(self::CasesBase{LES_driven_SCM}, grid::Grid, state,
 end
 
 initialize_forcing(self::CasesBase{LES_driven_SCM}, grid::Grid, state, gm, param_set) =
-    initialize(self.Fo, grid, self.LESDat)
+    initialize(self.Fo, grid, state, self.LESDat)
 
 initialize_radiation(self::CasesBase{LES_driven_SCM}, grid::Grid, state, gm, param_set) =
-    initialize(self.Rad, grid, self.LESDat)
+    initialize(self.Rad, grid, state, self.LESDat)
 
 end
diff --git a/integration_tests/utils/main.jl b/integration_tests/utils/main.jl
index e1033c23f..75810d4f0 100644
--- a/integration_tests/utils/main.jl
+++ b/integration_tests/utils/main.jl
@@ -61,6 +61,23 @@ cent_aux_vars_gm(FT) = (;
     H_third_m = FT(0),
     W_third_m = FT(0),
     QT_third_m = FT(0),
+    # From RadiationBase
+    dTdt_rad = FT(0), # horizontal advection temperature tendency
+    dqtdt_rad = FT(0), # horizontal advection moisture tendency
+    # From ForcingBase
+    subsidence = FT(0), #Large-scale subsidence
+    dTdt = FT(0), #Large-scale temperature tendency
+    dqtdt = FT(0), #Large-scale moisture tendency
+    dTdt_hadv = FT(0), #Horizontal advection of temperature
+    dTdt_nudge = FT(0), #Temperature tendency due to relaxation to large-scale
+    dTdt_fluc = FT(0), #Vertical turbulent advection of temperature
+    dqtdt_hadv = FT(0), #Horizontal advection of moisture
+    dqtdt_nudge = FT(0), #Moisture tendency due to relaxation to large-scale
+    dqtdt_fluc = FT(0), #Vertical turbulent advection of moisture
+    u_nudge = FT(0), #Reference u profile for relaxation tendency
+    v_nudge = FT(0), #Reference v profile for relaxation tendency
+    ug = FT(0), #Geostrophic u velocity
+    vg = FT(0), #Geostrophic v velocity
 )
 cent_aux_vars_en_2m(FT) = (;
     dissipation = FT(0),
@@ -169,7 +186,7 @@ cent_aux_vars_edmf(FT, n_up) = (;
 cent_aux_vars(FT, n_up) = (; aux_vars_ref_state(FT)..., cent_aux_vars_gm(FT)..., cent_aux_vars_edmf(FT, n_up)...)
 
 # Face only
-face_aux_vars_gm(FT) = (; massflux_s = FT(0), diffusive_flux_s = FT(0), total_flux_s = FT(0))
+face_aux_vars_gm(FT) = (; massflux_s = FT(0), diffusive_flux_s = FT(0), total_flux_s = FT(0), f_rad = FT(0))
 face_aux_vars_up(FT) = (;
     w = FT(0),
     nh_pressure = FT(0),
diff --git a/src/Fields.jl b/src/Fields.jl
index ac6a55253..269289f00 100644
--- a/src/Fields.jl
+++ b/src/Fields.jl
@@ -7,32 +7,6 @@ function center_field(grid::Grid, nu::Int)
     return zeros(nu, grid.nz)
 end
 
-function face_field(grid::Grid)
-    return zeros(grid.nz + 1)
-end
-
-function face_field(grid::Grid, nu::Int)
-    return zeros(nu, grid.nz + 1)
-end
-
-function field(grid::Grid, loc::String)
-    @assert loc == "full" || loc == "half"
-    if loc == "full"
-        return face_field(grid)
-    else
-        return center_field(grid)
-    end
-end
-
-function field(grid::Grid, loc::String, nu::Int)
-    @assert loc == "full" || loc == "half"
-    if loc == "full"
-        return face_field(grid, nu)
-    else
-        return center_field(grid, nu)
-    end
-end
-
 # A complementary struct to ClimaCore's `PlusHalf` type.
 struct Cent{I <: Integer}
     i::I
diff --git a/src/Forcing.jl b/src/Forcing.jl
index a6c828645..76bac5f2e 100644
--- a/src/Forcing.jl
+++ b/src/Forcing.jl
@@ -1,32 +1,33 @@
-function initialize(self::ForcingBase, grid, ::ForcingBaseType)
-    self.subsidence = center_field(grid)
-    self.dTdt = center_field(grid)
-    self.dqtdt = center_field(grid)
-    self.ug = center_field(grid)
-    self.vg = center_field(grid)
-    return
-end
-
-initialize(self::ForcingBase{ForcingNone}, grid) = initialize(self, grid, ForcingBaseType())
-initialize(self::ForcingBase{ForcingStandard}, grid) = initialize(self, grid, ForcingBaseType())
-initialize(self::ForcingBase{ForcingDYCOMS_RF01}, grid) = initialize(self, grid, ForcingBaseType())
+initialize(self::ForcingBase, grid, state) = nothing
 
-function initialize(self::ForcingBase{ForcingLES}, grid, LESDat::LESData)
-    initialize(self, grid, ForcingBaseType())
-    NC.Dataset(LESDat.les_filename, "r") do data
+function initialize(self::ForcingBase{ForcingLES}, grid, state, LESDat::LESData)
+    aux_gm = center_aux_grid_mean(state)
+    nt = NC.Dataset(LESDat.les_filename, "r") do data
         imin = LESDat.imin
         imax = LESDat.imax
 
         zc_les = Array(get_nc_data(data, "zc"))
 
-        self.dTdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax))
-        self.dTdt_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_nudge", imin, imax))
-        self.dTdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax))
-        self.dqtdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax))
-        self.dqtdt_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_nudge", imin, imax))
-        self.dqtdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_fluc", imin, imax))
-        self.subsidence = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "ls_subsidence", imin, imax))
-        self.u_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "u_mean", imin, imax))
-        self.v_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "v_mean", imin, imax))
+        dTdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_hadv", imin, imax))
+        dTdt_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_nudge", imin, imax))
+        dTdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dtdt_fluc", imin, imax))
+        dqtdt_hadv = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_hadv", imin, imax))
+        dqtdt_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_nudge", imin, imax))
+        dqtdt_fluc = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "dqtdt_fluc", imin, imax))
+        subsidence = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "ls_subsidence", imin, imax))
+        u_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "u_mean", imin, imax))
+        v_nudge = pyinterp(grid.zc, zc_les, mean_nc_data(data, "profiles", "v_mean", imin, imax))
+        (; dTdt_hadv, dTdt_nudge, dTdt_fluc, dqtdt_hadv, dqtdt_nudge, dqtdt_fluc, subsidence, u_nudge, v_nudge)
+    end
+    for k in real_center_indices(grid)
+        aux_gm.dTdt_hadv[k] = nt.dTdt_hadv[k]
+        aux_gm.dTdt_nudge[k] = nt.dTdt_nudge[k]
+        aux_gm.dTdt_fluc[k] = nt.dTdt_fluc[k]
+        aux_gm.dqtdt_hadv[k] = nt.dqtdt_hadv[k]
+        aux_gm.dqtdt_nudge[k] = nt.dqtdt_nudge[k]
+        aux_gm.dqtdt_fluc[k] = nt.dqtdt_fluc[k]
+        aux_gm.subsidence[k] = nt.subsidence[k]
+        aux_gm.u_nudge[k] = nt.u_nudge[k]
+        aux_gm.v_nudge[k] = nt.v_nudge[k]
     end
 end
diff --git a/src/NetCDFIO.jl b/src/NetCDFIO.jl
index ecd00a6f5..c0c502bd1 100644
--- a/src/NetCDFIO.jl
+++ b/src/NetCDFIO.jl
@@ -1,10 +1,4 @@
 
-face_fields_list() = ()
-
-function is_face_field(var_name)
-    return var_name in face_fields_list()
-end
-
 # TODO: remove `vars` hack that avoids https://github.com/Alexander-Barth/NCDatasets.jl/issues/135
 
 mutable struct NetCDFIO_Stats
@@ -127,14 +121,6 @@ function close_files(self::NetCDFIO_Stats)
     close(self.root_grp)
 end
 
-#####
-##### Profile-specific
-#####
-
-# TODO: depricate
-add_profile(self::NetCDFIO_Stats, var_name::String) =
-    add_field(self, var_name; dims = is_face_field(var_name) ? ("zf", "t") : ("zc", "t"), group = "profiles")
-
 #####
 ##### Generic field
 #####
@@ -164,10 +150,6 @@ end
 # Field wrapper
 write_field(self::NetCDFIO_Stats, var_name::String, data; group) = write_field(self, var_name, vec(data); group = group)
 
-# TODO: depricate
-write_profile(self::NetCDFIO_Stats, var_name::String, data::T) where {T <: AbstractArray{Float64, 1}} =
-    write_field(self, var_name, data; group = "profiles")
-
 function write_field(self::NetCDFIO_Stats, var_name::String, data::T; group) where {T <: AbstractArray{Float64, 1}}
     if group == "profiles"
         # Hack to avoid https://github.com/Alexander-Barth/NCDatasets.jl/issues/135
diff --git a/src/Radiation.jl b/src/Radiation.jl
index 10cfd1717..487d345f6 100644
--- a/src/Radiation.jl
+++ b/src/Radiation.jl
@@ -1,18 +1,7 @@
-
-function initialize(self::RadiationBase, grid, ::RadiationBaseType)
-    self.dTdt = center_field(grid)
-    self.dqtdt = center_field(grid)
-    return
-end
-
 update(self::RadiationBase, grid, state, gm::GridMeanVariables, param_set) = nothing
+initialize(self::RadiationBase{RadiationNone}, grid, state) = nothing
 
-initialize(self::RadiationBase{RadiationNone}, grid) = initialize(self, grid, RadiationBaseType())
-
-function initialize(self::RadiationBase{RadiationDYCOMS_RF01}, grid)
-    initialize(self, grid, RadiationBaseType())
-
-
+function initialize(self::RadiationBase{RadiationDYCOMS_RF01}, grid, state)
     self.divergence = 3.75e-6  # divergence is defined twice: here and in initialize_forcing method of DYCOMS_RF01 case class
     # where it is used to initialize large scale subsidence
     self.alpha_z = 1.0
@@ -20,7 +9,6 @@ function initialize(self::RadiationBase{RadiationDYCOMS_RF01}, grid)
     self.F0 = 70.0
     self.F1 = 22.0
     self.divergence = 3.75e-6
-    self.f_rad = face_field(grid)
     return
 end
 
@@ -32,6 +20,7 @@ function calculate_radiation(self::RadiationBase{RadiationDYCOMS_RF01}, grid, st
     ρ0_f = face_ref_state(state).ρ0
     ρ0_c = center_ref_state(state).ρ0
     aux_gm = center_aux_grid_mean(state)
+    aux_gm_f = face_aux_grid_mean(state)
     prog_gm = center_prog_grid_mean(state)
     # find zi (level of 8.0 g/kg isoline of qt)
     # TODO: report bug: zi and ρ_i are not initialized
@@ -67,21 +56,21 @@ function calculate_radiation(self::RadiationBase{RadiationDYCOMS_RF01}, grid, st
     prob = ODE.ODEProblem(integrand, 0.0, z_span, params; dt = grid.Δz)
     sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12)
     q_1 = sol.(vec(grid.zf))
-    self.f_rad .= self.F0 .* exp.(-q_0)
-    self.f_rad .+= self.F1 .* exp.(-q_1)
+    parent(aux_gm_f.f_rad) .= self.F0 .* exp.(-q_0)
+    parent(aux_gm_f.f_rad) .+= self.F1 .* exp.(-q_1)
 
     # cooling in free troposphere
     @inbounds for k in real_face_indices(grid)
         if grid.zf[k] > zi
             cbrt_z = cbrt(grid.zf[k] - zi)
-            self.f_rad[k] += ρ_i * cp_d * self.divergence * self.alpha_z * (cbrt_z^4 / 4 + zi * cbrt_z)
+            aux_gm_f.f_rad[k] += ρ_i * cp_d * self.divergence * self.alpha_z * (cbrt_z^4 / 4 + zi * cbrt_z)
         end
     end
 
     @inbounds for k in real_center_indices(grid)
-        f_rad_dual = dual_faces(self.f_rad, grid, k)
+        f_rad_dual = dual_faces(aux_gm_f.f_rad, grid, k)
         ∇f_rad = ∇_staggered(f_rad_dual, grid)
-        self.dTdt[k] = -∇f_rad / ρ0_c[k] / cp_d
+        aux_gm.dTdt_rad[k] = -∇f_rad / ρ0_c[k] / cp_d
     end
 
     return
@@ -90,17 +79,20 @@ end
 update(self::RadiationBase{RadiationDYCOMS_RF01}, grid, state, gm::GridMeanVariables, param_set) =
     calculate_radiation(self, grid, state, gm, param_set)
 
-function initialize(self::RadiationBase{RadiationLES}, grid, LESDat::LESData)
-    initialize(self, grid, RadiationBaseType())
+function initialize(self::RadiationBase{RadiationLES}, grid, state, LESDat::LESData)
     # load from LES
-    NC.Dataset(LESDat.les_filename, "r") do data
+    aux_gm = center_aux_grid_mean(state)
+    dTdt = NC.Dataset(LESDat.les_filename, "r") do data
         imin = LESDat.imin
         imax = LESDat.imax
 
         # interpolate here
         zc_les = Array(get_nc_data(data, "zc"))
         meandata = mean_nc_data(data, "profiles", "dtdt_rad", imin, imax)
-        self.dTdt = pyinterp(grid.zc, zc_les, meandata)
+        pyinterp(grid.zc, zc_les, meandata)
+    end
+    @inbounds for k in real_center_indices(grid)
+        aux_gm.dTdt_rad[k] = dTdt[k]
     end
     return
 end
diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl
index 3d2abb42f..cb76d6e83 100755
--- a/src/Turbulence_PrognosticTKE.jl
+++ b/src/Turbulence_PrognosticTKE.jl
@@ -21,6 +21,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
     tendencies_gm = center_tendencies_grid_mean(state)
     param_set = parameter_set(gm)
     prog_gm = center_prog_grid_mean(state)
+    aux_gm = center_aux_grid_mean(state)
     aux_en = center_aux_environment(state)
     aux_en_f = face_aux_environment(state)
     aux_up = center_aux_updrafts(state)
@@ -42,11 +43,11 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         Π = TD.exner(ts)
 
         if Case.Fo.apply_coriolis
-            tendencies_gm.u[k] -= Case.Fo.coriolis_param * (Case.Fo.vg[k] - prog_gm.v[k])
-            tendencies_gm.v[k] += Case.Fo.coriolis_param * (Case.Fo.ug[k] - prog_gm.u[k])
+            tendencies_gm.u[k] -= Case.Fo.coriolis_param * (aux_gm.vg[k] - prog_gm.v[k])
+            tendencies_gm.v[k] += Case.Fo.coriolis_param * (aux_gm.ug[k] - prog_gm.u[k])
         end
         if rad_type(Case.Rad) <: Union{RadiationDYCOMS_RF01, RadiationLES}
-            tendencies_gm.θ_liq_ice[k] += Case.Rad.dTdt[k] / Π
+            tendencies_gm.θ_liq_ice[k] += aux_gm.dTdt_rad[k] / Π
         end
         H_cut = ccut_downwind(prog_gm.θ_liq_ice, grid, k)
         q_tot_cut = ccut_downwind(prog_gm.q_tot, grid, k)
@@ -54,32 +55,32 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
         ∇q_tot = c∇_downwind(q_tot_cut, grid, k; bottom = FreeBoundary(), top = SetGradient(0))
 
         if force_type(Case.Fo) <: ForcingDYCOMS_RF01
-            tendencies_gm.q_tot[k] += Case.Fo.dqtdt[k]
+            tendencies_gm.q_tot[k] += aux_gm.dqtdt[k]
             # Apply large-scale subsidence tendencies
-            tendencies_gm.θ_liq_ice[k] -= ∇H * Case.Fo.subsidence[k]
-            tendencies_gm.q_tot[k] -= ∇q_tot * Case.Fo.subsidence[k]
+            tendencies_gm.θ_liq_ice[k] -= ∇H * aux_gm.subsidence[k]
+            tendencies_gm.q_tot[k] -= ∇q_tot * aux_gm.subsidence[k]
         end
 
         if force_type(Case.Fo) <: ForcingStandard
             if Case.Fo.apply_subsidence
-                tendencies_gm.θ_liq_ice[k] -= ∇H * Case.Fo.subsidence[k]
-                tendencies_gm.q_tot[k] -= ∇q_tot * Case.Fo.subsidence[k]
+                tendencies_gm.θ_liq_ice[k] -= ∇H * aux_gm.subsidence[k]
+                tendencies_gm.q_tot[k] -= ∇q_tot * aux_gm.subsidence[k]
             end
-            tendencies_gm.θ_liq_ice[k] += Case.Fo.dTdt[k] / Π
-            tendencies_gm.q_tot[k] += Case.Fo.dqtdt[k]
+            tendencies_gm.θ_liq_ice[k] += aux_gm.dTdt[k] / Π
+            tendencies_gm.q_tot[k] += aux_gm.dqtdt[k]
         end
 
         if force_type(Case.Fo) <: ForcingLES
-            H_horz_adv = Case.Fo.dTdt_hadv[k] / Π
-            H_nudge = Case.Fo.dTdt_nudge[k] / Π
-            H_fluc = Case.Fo.dTdt_fluc[k] / Π
+            H_horz_adv = aux_gm.dTdt_hadv[k] / Π
+            H_nudge = aux_gm.dTdt_nudge[k] / Π
+            H_fluc = aux_gm.dTdt_fluc[k] / Π
 
-            gm_U_nudge_k = (Case.Fo.u_nudge[k] - prog_gm.u[k]) / Case.Fo.nudge_tau
-            gm_V_nudge_k = (Case.Fo.v_nudge[k] - prog_gm.v[k]) / Case.Fo.nudge_tau
+            gm_U_nudge_k = (aux_gm.u_nudge[k] - prog_gm.u[k]) / Case.Fo.nudge_tau
+            gm_V_nudge_k = (aux_gm.v_nudge[k] - prog_gm.v[k]) / Case.Fo.nudge_tau
             if Case.Fo.apply_subsidence
                 # Apply large-scale subsidence tendencies
-                gm_H_subsidence_k = -∇H * Case.Fo.subsidence[k]
-                gm_QT_subsidence_k = -∇q_tot * Case.Fo.subsidence[k]
+                gm_H_subsidence_k = -∇H * aux_gm.subsidence[k]
+                gm_QT_subsidence_k = -∇q_tot * aux_gm.subsidence[k]
             else
                 gm_H_subsidence_k = 0.0
                 gm_QT_subsidence_k = 0.0
@@ -87,7 +88,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, state, Case, gm,
 
             tendencies_gm.θ_liq_ice[k] += H_horz_adv + H_nudge + H_fluc + gm_H_subsidence_k
             tendencies_gm.q_tot[k] +=
-                Case.Fo.dqtdt_hadv[k] + Case.Fo.dqtdt_nudge[k] + gm_QT_subsidence_k + Case.Fo.dqtdt_fluc[k]
+                aux_gm.dqtdt_hadv[k] + aux_gm.dqtdt_nudge[k] + gm_QT_subsidence_k + aux_gm.dqtdt_fluc[k]
 
             tendencies_gm.u[k] += gm_U_nudge_k
             tendencies_gm.v[k] += gm_V_nudge_k
diff --git a/src/closures/entr_detr.jl b/src/closures/entr_detr.jl
index 15a875a40..3c26b46be 100644
--- a/src/closures/entr_detr.jl
+++ b/src/closures/entr_detr.jl
@@ -13,7 +13,6 @@ Cohen et al. (JAMES, 2020), given:
 """
 function entr_detr(param_set, εδ_model::MoistureDeficitEntr)
 
-    l = zeros(2)
     γ_lim = ICP.area_limiter_scale(param_set)
     β_lim = ICP.area_limiter_power(param_set)
     c_ε = CPEDMF.c_ε(param_set)
diff --git a/src/diagnostics.jl b/src/diagnostics.jl
index aa7c3fca7..6cd68c3d7 100644
--- a/src/diagnostics.jl
+++ b/src/diagnostics.jl
@@ -142,6 +142,10 @@ function io_dictionary_aux(state)
         "ed_length_scheme" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).mls),
         "mixing_length_ratio" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).ml_ratio),
         "entdet_balance_length" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_turbconv(state).l_entdet),
+
+        "rad_dTdt" => (; dims = ("zc", "t"), group = "profiles", field = center_aux_grid_mean(state).dTdt_rad),
+        "rad_flux" => (; dims = ("zf", "t"), group = "profiles", field = face_aux_grid_mean(state).f_rad),
+
     )
     return io_dict
 end
diff --git a/src/io.jl b/src/io.jl
index 85340f1a8..237a39e90 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -5,30 +5,6 @@ initialize_io(self::ForcingBase, Stats) = nothing
 initialize_io(self::RadiationBase, Stats::NetCDFIO_Stats) = nothing
 io(self::RadiationBase, grid, state, Stats::NetCDFIO_Stats) = nothing
 
-function initialize_io(self::RadiationBase{RadiationLES}, Stats::NetCDFIO_Stats)
-    add_profile(Stats, "rad_dTdt")
-    add_profile(Stats, "rad_flux")
-    return
-end
-
-function io(self::RadiationBase{RadiationLES}, Stats::NetCDFIO_Stats)
-    write_profile(Stats, "rad_dTdt", self.dTdt)
-    write_profile(Stats, "rad_flux", self.f_rad)
-    return
-end
-
-function initialize_io(self::RadiationBase{RadiationDYCOMS_RF01}, Stats::NetCDFIO_Stats)
-    add_profile(Stats, "rad_dTdt")
-    add_profile(Stats, "rad_flux")
-    return
-end
-
-function io(self::RadiationBase{RadiationDYCOMS_RF01}, Stats::NetCDFIO_Stats)
-    write_profile(Stats, "rad_dTdt", self.dTdt)
-    write_profile(Stats, "rad_flux", self.f_rad)
-    return
-end
-
 function initialize_io(en::EnvironmentVariables, Stats::NetCDFIO_Stats)
     add_ts(Stats, "env_cloud_base")
     add_ts(Stats, "env_cloud_top")
diff --git a/src/types.jl b/src/types.jl
index 0c22bccd1..c2270c888 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -351,38 +351,12 @@ LES-driven forcing
 $(DocStringExtensions.FIELDS)
 """
 Base.@kwdef mutable struct ForcingBase{T}
-    "Large-scale subsidence"
-    subsidence::AbstractArray{Float64, 1} = zeros(1)
-    "Large-scale temperature tendency"
-    dTdt::AbstractArray{Float64, 1} = zeros(1)
-    "Large-scale moisture tendency"
-    dqtdt::AbstractArray{Float64, 1} = zeros(1)
-    "Horizontal advection of temperature"
-    dTdt_hadv::AbstractArray{Float64, 1} = zeros(1)
-    "Temperature tendency due to relaxation to large-scale"
-    dTdt_nudge::AbstractArray{Float64, 1} = zeros(1)
-    "Vertical turbulent advection of temperature"
-    dTdt_fluc::AbstractArray{Float64, 1} = zeros(1)
-    "Horizontal advection of moisture"
-    dqtdt_hadv::AbstractArray{Float64, 1} = zeros(1)
-    "Moisture tendency due to relaxation to large-scale"
-    dqtdt_nudge::AbstractArray{Float64, 1} = zeros(1)
-    "Vertical turbulent advection of moisture"
-    dqtdt_fluc::AbstractArray{Float64, 1} = zeros(1)
-    "Reference u profile for relaxation tendency"
-    u_nudge::AbstractArray{Float64, 1} = zeros(1)
-    "Reference v profile for relaxation tendency"
-    v_nudge::AbstractArray{Float64, 1} = zeros(1)
     "Boolean specifying whether Coriolis forcing is applied"
     apply_coriolis::Bool = false
     "Boolean specifying whether subsidence forcing is applied"
     apply_subsidence::Bool = false
     "Coriolis parameter"
     coriolis_param::Float64 = 0
-    "Geostrophic u velocity"
-    ug::AbstractArray{Float64, 1} = zeros(1)
-    "Geostrophic v velocity"
-    vg::AbstractArray{Float64, 1} = zeros(1)
     "Momentum relaxation timescale"
     nudge_tau::Float64 = 0.0
     "Conversion function from forcing to prognostic"
@@ -392,14 +366,11 @@ end
 force_type(::ForcingBase{T}) where {T} = T
 
 Base.@kwdef mutable struct RadiationBase{T}
-    dTdt::AbstractArray{Float64, 1} = zeros(1) # horizontal advection temperature tendency
-    dqtdt::AbstractArray{Float64, 1} = zeros(1) # horizontal advection moisture tendency
     divergence::Float64 = 0
     alpha_z::Float64 = 0
     kappa::Float64 = 0
     F0::Float64 = 0
     F1::Float64 = 0
-    f_rad::AbstractArray{Float64, 1} = zeros(1)
 end
 
 rad_type(::RadiationBase{T}) where {T} = T