Sign consistency with R_n

CliMA · Oct 7, 2024 · 0b68e7c · 0b68e7c
1 parent 8a4cfba
commit 0b68e7c
Show file tree

Hide file tree

Showing 17 changed files with 71 additions and 44 deletions.
diff --git a/docs/make.jl b/docs/make.jl
@@ -52,6 +52,7 @@ include("list_diagnostics.jl")
 pages = Any[
     "Home" => "index.md",
     "Getting Started" => "getting_started.md",
+    "Variables used in ClimaLand" => "glossary.md",
     "Tutorials" => tutorials,
     "Standalone models" => standalone_models,
     "Diagnostics" => diagnostics,

diff --git a/docs/src/glossary.md b/docs/src/glossary.md
@@ -0,0 +1,36 @@
+# Definition of variables and terms
+
+## Sign and Units conventions
+
+### Basic units
+All variables are in standard units: `kg` for mass, `m` for length, and
+seconds `s` for time.  Zenith angle is reported in radians. We sometimes
+use the units of moles and number of photons in the Canopy model where noted.
+
+Currently, dates are used only for interfacing with a calendar
+(for forcing data, for output, for the insolation), and the simulation time
+is in units of seconds from a reference date.
+
+### Units for fluxes
+- Energy fluxes have units of W/m^2.
+- Water fluxes have units of m^3/m^2/s = m/s. We are planning to change this to a mass flux.
+- Carbon fluxes have units of moles CO2/m^2/s or kg C/m^2/s, depending on the context. We are planning to change this to be consistent throughout.
+
+### Sign conventions for fluxes
+The majority of fluxes follow the convention that upwards is positive and
+downwards is negative. So, we have that:
+- Sensible, latent, and vapor fluxes are positive if towards the atmosphere.
+- Snowmelt is negative.
+- Root extraction is positive if the canopy is taking water from the soil.
+- The ground heat flux between soil and snow is positive if the snow is warming.
+
+Some fluxes, however, are represented only with scalars, with a direction
+built in, such as:
+- Precipitation is positive by definition, but always is downwards.
+- Radiative fluxes marked with `_d` are downwelling and positive by
+definition.
+- Radiative fluxes marked with `_u` are upwelling and positive by definition.
+- Net radiation is defined to be `R_n = SW_d + LW_d - SW_u - LW_u`, where
+`SW` and `LW indicate the total short and longwave fluxes. Therefore, the
+ClimaLand convention is that `R_n` is positive if the land is gaining energy.
+- Runoff is defined to be positive.
diff --git a/docs/tutorials/standalone/Bucket/bucket_tutorial.jl b/docs/tutorials/standalone/Bucket/bucket_tutorial.jl
@@ -47,11 +47,11 @@
 # ``
 
 # ``
-# (1-σ) (R_n+ SHF + LHF)_{soil} + σG_{undersnow} = -κ_{soil} \frac{\partial T}{\partial z}|_{z = z_{sfc}}
+# (1-σ) (SHF + LHF - R_n)_{soil} + σG_{undersnow} = -κ_{soil} \frac{\partial T}{\partial z}|_{z = z_{sfc}}
 # ``
 
 # ``
-# G_{undersnow} = (R_n+ SHF + LHF)_{snow} - F_{intosnow}
+# G_{undersnow} = (SHF + LHF - R_n)_{snow} - F_{intosnow}
 # ``
 
 # ``
@@ -60,7 +60,7 @@
 
 
 # ``
-# R_n = -(1-α)*SW↓ -LW↓ + σ_{SB} T_{sfc}^4
+# R_n = (1-α)*SW↓ + LW↓ - σ_{SB} T_{sfc}^4
 # ``
 
 # where the water fluxes are : `I` the infiltration as defined in [1], `P_liq` (m/s) the
@@ -419,7 +419,7 @@ plot(
     legend = :bottomright,
 )
 plot!(sol.t ./ 86400, turbulent_energy_flux, label = "Turbulent fluxes")
-plot!(sol.t ./ 86400, R_n .+ turbulent_energy_flux, label = "Net flux")
+plot!(sol.t ./ 86400, turbulent_energy_flux .- R_n, label = "Net flux")
 savefig("energy_f.png")
 # ![](energy_f.png)
 

diff --git a/experiments/integrated/performance/conservation/ozark_conservation.jl b/experiments/integrated/performance/conservation/ozark_conservation.jl
@@ -224,8 +224,7 @@ for float_type in (Float32, Float64)
         ]
         # Radiation is computed in LW and SW components
         # with positive numbers indicating the soil gaining energy.
-        soil_Rn =
-            -1 .* [parent(sv.saveval[k].soil.R_n)[1] for k in 2:length(sol.t)]
+        soil_Rn = [parent(sv.saveval[k].soil.R_n)[1] for k in 2:length(sol.t)]
         # Root sink term: a positive root extraction is a sink term for soil; add minus sign
         root_sink_energy = [
             sum(-1 .* sv.saveval[k].root_energy_extraction) for
@@ -239,7 +238,7 @@ for float_type in (Float32, Float64)
         # d[∫Idz] = [-(F_sfc - F_bot) + ∫Sdz]dt = -ΔF dt + ∫Sdz dt
         # N.B. in ClimaCore, sum(field) -> integral
         rhs_soil_energy =
-            -(LHF .+ SHF .+ soil_Rn .- soil_bottom_flux) .+ root_sink_energy
+            -(LHF .+ SHF .- soil_Rn .- soil_bottom_flux) .+ root_sink_energy
 
         net_soil_energy_storage =
             [sum(sol.u[k].soil.ρe_int)[1] for k in 1:length(sol.t)]

diff --git a/experiments/standalone/Bucket/bucket_era5.jl b/experiments/standalone/Bucket/bucket_era5.jl
@@ -333,7 +333,7 @@ F_sfc = [
     Array(
         Remapping.interpolate(
             remapper,
-            saved_values.saveval[k].bucket.R_n .+
+            -1 .* saved_values.saveval[k].bucket.R_n .+
             saved_values.saveval[k].bucket.turbulent_fluxes.lhf .+
             saved_values.saveval[k].bucket.turbulent_fluxes.shf,
         ),

diff --git a/experiments/standalone/Bucket/global_bucket_temporalmap.jl b/experiments/standalone/Bucket/global_bucket_temporalmap.jl
@@ -271,7 +271,7 @@ F_sfc = [
     Array(
         Remapping.interpolate(
             remapper,
-            saved_values.saveval[k].bucket.R_n .+
+            -1 .* saved_values.saveval[k].bucket.R_n .+
             saved_values.saveval[k].bucket.turbulent_fluxes.lhf .+
             saved_values.saveval[k].bucket.turbulent_fluxes.shf,
         ),

diff --git a/src/integrated/soil_canopy_model.jl b/src/integrated/soil_canopy_model.jl
@@ -379,7 +379,6 @@ function lsm_radiant_energy_fluxes!(
     ϵ_canopy = p.canopy.radiative_transfer.ϵ # this takes into account LAI/SAI
     @. LW_d_canopy = ((1 - ϵ_canopy) * LW_d + ϵ_canopy * _σ * T_canopy^4) # double checked
     @. LW_u_soil = ϵ_soil * _σ * p.T_ground^4 + (1 - ϵ_soil) * LW_d_canopy # double checked
-    # This is a sign inconsistency. Here Rn is positive if towards soil. X_X
     @. R_net_soil += ϵ_soil * LW_d_canopy - ϵ_soil * _σ * p.T_ground^4 # double checked
     @. LW_net_canopy =
         ϵ_canopy * LW_d - 2 * ϵ_canopy * _σ * T_canopy^4 + ϵ_canopy * LW_u_soil

diff --git a/src/integrated/soil_snow_model.jl b/src/integrated/soil_snow_model.jl
@@ -203,13 +203,11 @@ function make_update_boundary_fluxes(
         ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
         @. p.atmos_energy_flux =
             (1 - p.snow.snow_cover_fraction) * (
-                p.soil.turbulent_fluxes.lhf +
-                p.soil.turbulent_fluxes.shf +
+                p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf -
                 p.soil.R_n
             ) +
             p.snow.snow_cover_fraction * (
-                p.snow.turbulent_fluxes.lhf +
-                p.snow.turbulent_fluxes.shf +
+                p.snow.turbulent_fluxes.lhf + p.snow.turbulent_fluxes.shf -
                 p.snow.R_n
             ) +
             p.drivers.P_snow * ρe_falling_snow
@@ -359,8 +357,7 @@ function snow_boundary_fluxes!(
     @. p.snow.total_energy_flux =
         P_snow * ρe_falling_snow +
         (
-            p.snow.turbulent_fluxes.lhf +
-            p.snow.turbulent_fluxes.shf +
+            p.snow.turbulent_fluxes.lhf + p.snow.turbulent_fluxes.shf -
             p.snow.R_n - p.snow.energy_runoff - p.ground_heat_flux
         ) * p.snow.snow_cover_fraction
 end
@@ -407,7 +404,7 @@ function soil_boundary_fluxes!(
 
     @. p.soil.top_bc.heat =
         (1 - p.snow.snow_cover_fraction) * (
-            p.soil.R_n +
+            -p.soil.R_n +
             p.soil.turbulent_fluxes.lhf +
             p.soil.turbulent_fluxes.shf
         ) +

diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl
@@ -461,10 +461,7 @@ function net_radiation(
     T_sfc = surface_temperature(model, Y, p, t)
     α_sfc = surface_albedo(model, Y, p)
     ϵ_sfc = surface_emissivity(model, Y, p)
-    # Recall that the user passed the LW and SW downwelling radiation,
-    # where positive values indicate toward surface, so we need a negative sign out front
-    # in order to inidicate positive R_n  = towards atmos.
-    R_n = @.(-(1 - α_sfc) * SW_d - ϵ_sfc * (LW_d - _σ * T_sfc^4))
+    R_n = @.((1 - α_sfc) * SW_d + ϵ_sfc * (LW_d - _σ * T_sfc^4)) # positive if the land absorbs energy
     return R_n
 end
 

diff --git a/src/standalone/Bucket/Bucket.jl b/src/standalone/Bucket/Bucket.jl
@@ -449,7 +449,7 @@ function make_update_aux(model::BucketModel{FT}) where {FT}
         # The below must be adjusted to compute F_sfc over snow and over soil
         # if we want the snow cover fraction to be intermediate between 0 and 1.
         @. p.bucket.F_sfc = (
-            p.bucket.turbulent_fluxes.shf .+ p.bucket.turbulent_fluxes.lhf +
+            p.bucket.turbulent_fluxes.shf .+ p.bucket.turbulent_fluxes.lhf -
             p.bucket.R_n
         ) # Eqn (21)
         _T_freeze = LP.T_freeze(model.parameters.earth_param_set)

diff --git a/src/standalone/Snow/boundary_fluxes.jl b/src/standalone/Snow/boundary_fluxes.jl
@@ -122,8 +122,7 @@ function snow_boundary_fluxes!(
     @. p.snow.total_energy_flux =
         P_snow * ρe_falling_snow +
         (
-            p.snow.turbulent_fluxes.lhf +
-            p.snow.turbulent_fluxes.shf +
+            p.snow.turbulent_fluxes.lhf + p.snow.turbulent_fluxes.shf -
             p.snow.R_n - p.snow.energy_runoff
         ) * p.snow.snow_cover_fraction
 end
diff --git a/src/standalone/Soil/boundary_conditions.jl b/src/standalone/Soil/boundary_conditions.jl
@@ -800,7 +800,7 @@ function soil_boundary_fluxes!(
     @. p.soil.top_bc.water =
         p.soil.infiltration + p.soil.turbulent_fluxes.vapor_flux_liq
     @. p.soil.top_bc.heat =
-        p.soil.R_n + p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf
+        p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf - p.soil.R_n
 end
 
 """

diff --git a/test/integrated/soil_snow.jl b/test/integrated/soil_snow.jl
@@ -192,13 +192,11 @@ _LH_f0 = FT(LP.LH_f0(earth_param_set))
 _ρ_liq = FT(LP.ρ_cloud_liq(earth_param_set))
 ρe_falling_snow = -_LH_f0 * _ρ_liq # per unit vol of liquid water
 @test p.atmos_energy_flux == @. (1 - p.snow.snow_cover_fraction) * (
-             p.soil.turbulent_fluxes.lhf +
-             p.soil.turbulent_fluxes.shf +
+             p.soil.turbulent_fluxes.lhf + p.soil.turbulent_fluxes.shf -
              p.soil.R_n
          ) +
          p.snow.snow_cover_fraction * (
-             p.snow.turbulent_fluxes.lhf +
-             p.snow.turbulent_fluxes.shf +
+             p.snow.turbulent_fluxes.lhf + p.snow.turbulent_fluxes.shf -
              p.snow.R_n
          ) +
          p.drivers.P_snow * ρe_falling_snow

diff --git a/test/standalone/Bucket/snow_bucket_tests.jl b/test/standalone/Bucket/snow_bucket_tests.jl
@@ -150,9 +150,9 @@ for FT in (Float32, Float64)
                 ),
             ) .== FT(0.0),
         )
-        F_sfc = @. p.bucket.turbulent_fluxes.lhf +
-           p.bucket.turbulent_fluxes.shf +
-           p.bucket.R_n
+        F_sfc =
+            @. p.bucket.turbulent_fluxes.lhf + p.bucket.turbulent_fluxes.shf -
+               p.bucket.R_n
         partitioned_fluxes = @. ClimaLand.Bucket.partition_snow_surface_fluxes(
             Y.bucket.σS,
             p.bucket.T_sfc,
@@ -252,7 +252,7 @@ for FT in (Float32, Float64)
                     snow_cover_fraction.(Y.bucket.σS),
                     p.bucket.turbulent_fluxes.vapor_flux,
                     p.bucket.turbulent_fluxes.lhf .+
-                    p.bucket.turbulent_fluxes.shf .+ p.bucket.R_n,
+                    p.bucket.turbulent_fluxes.shf .- p.bucket.R_n,
                     _ρLH_f0,
                     _T_freeze,
                 )
@@ -262,7 +262,7 @@ for FT in (Float32, Float64)
             G = partitioned_fluxes.G_under_snow
             F_sfc =
                 p.bucket.turbulent_fluxes.lhf .+
-                p.bucket.turbulent_fluxes.shf .+ p.bucket.R_n .-
+                p.bucket.turbulent_fluxes.shf .- p.bucket.R_n .-
                 _ρLH_f0 .* FT(snow_precip(t0))
             F_water_sfc =
                 FT(liquid_precip(t0)) + FT(snow_precip(t0)) .+
@@ -366,7 +366,7 @@ for FT in (Float32, Float64)
                     snow_cover_fraction.(Y.bucket.σS),
                     p.bucket.turbulent_fluxes.vapor_flux,
                     p.bucket.turbulent_fluxes.lhf .+
-                    p.bucket.turbulent_fluxes.shf .+ p.bucket.R_n,
+                    p.bucket.turbulent_fluxes.shf .- p.bucket.R_n,
                     _ρLH_f0,
                     _T_freeze,
                 )
@@ -376,7 +376,7 @@ for FT in (Float32, Float64)
             G = partitioned_fluxes.G_under_snow
             F_sfc =
                 p.bucket.turbulent_fluxes.lhf .+
-                p.bucket.turbulent_fluxes.shf .+ p.bucket.R_n .-
+                p.bucket.turbulent_fluxes.shf .- p.bucket.R_n .-
                 _ρLH_f0 .* FT(snow_precip(t0))
             F_water_sfc =
                 FT(liquid_precip(t0)) + FT(snow_precip(t0)) .+
@@ -475,9 +475,9 @@ for FT in (Float32, Float64)
             return σS > eps(typeof(σS)) ? typeof(σS)(1.0) : typeof(σS)(0.0)
         end
         σ_snow = snow_cover_fraction.(Y.bucket.σS)
-        F_sfc = @. p.bucket.turbulent_fluxes.lhf +
-           p.bucket.turbulent_fluxes.shf +
-           p.bucket.R_n
+        F_sfc =
+            @. p.bucket.turbulent_fluxes.lhf + p.bucket.turbulent_fluxes.shf -
+               p.bucket.R_n
         partitioned_fluxes =
             partition_snow_surface_fluxes.(
                 Y.bucket.σS,

diff --git a/test/standalone/Bucket/soil_bucket_tests.jl b/test/standalone/Bucket/soil_bucket_tests.jl
@@ -215,7 +215,7 @@ for FT in (Float32, Float64)
             F_water_sfc = FT(precip(t0)) .+ p.bucket.turbulent_fluxes.vapor_flux
             F_sfc =
                 p.bucket.turbulent_fluxes.lhf .+
-                p.bucket.turbulent_fluxes.shf .+ p.bucket.R_n
+                p.bucket.turbulent_fluxes.shf .- p.bucket.R_n
             surface_space = model.domain.space.surface
             A_sfc = sum(ones(surface_space))
 

diff --git a/test/standalone/Snow/snow.jl b/test/standalone/Snow/snow.jl
@@ -84,7 +84,7 @@ import ClimaLand.Parameters as LP
 
     # Check if aux update occurred correctly
     @test p.snow.R_n ==
-          @. (-(1 - α_snow) * 20.0f0 - ϵ_snow * (20.0f0 - _σ * p.snow.T_sfc^4))
+          @. ((1 - α_snow) * 20.0f0 + ϵ_snow * (20.0f0 - _σ * p.snow.T_sfc^4))
     @test p.snow.R_n == ClimaLand.net_radiation(
         model.boundary_conditions.radiation,
         model,
@@ -143,8 +143,9 @@ import ClimaLand.Parameters as LP
     )
     @test dY.snow.S == net_water_fluxes
     @test dY.snow.U == @.(
-        -p.snow.turbulent_fluxes.shf - p.snow.turbulent_fluxes.lhf -
-        p.snow.R_n + p.snow.energy_runoff
+        -p.snow.turbulent_fluxes.shf - p.snow.turbulent_fluxes.lhf +
+        p.snow.R_n +
+        p.snow.energy_runoff
     )
 
 

diff --git a/test/standalone/Soil/climate_drivers.jl b/test/standalone/Soil/climate_drivers.jl
@@ -229,7 +229,7 @@ for FT in (Float32, Float64)
 
             expected_water_flux = @. FT(precip(t)) .+ conditions.vapor_flux_liq
             @test computed_water_flux == expected_water_flux
-            expected_energy_flux = @. R_n + conditions.lhf + conditions.shf
+            expected_energy_flux = @. conditions.lhf + conditions.shf - R_n
             @test computed_energy_flux == expected_energy_flux
 
             # Test soil resistances for liquid water