use net radiative fluxes [skip ci]

juliasloan25 · juliasloan25 · commit e56e5ca45107 · 2025-03-14T14:06:42.000-07:00
diff --git a/bucket_restart/output_active b/bucket_restart/output_active
@@ -0,0 +1 @@
+output_0000
diff --git a/experiments/long_runs/data_paper_plots.jl b/experiments/long_runs/data_paper_plots.jl
@@ -74,42 +74,46 @@ function get_sim_var_dict(simdir)
             return sim_var_swd
         end
 
-    # # Compute LWN = LWD - LWU
-    # sim_var_dict["lwn"] =
-    #     () -> begin
-    #         lwd = sim_var_dict["lwd"]()
-    #         lwu = sim_var_dict["lwu"]()
-    #         lwn = lwd - lwu
-    #         lwn.attributes["long_name"] = "Net longwave radiation"
-    #         lwn.attributes["units"] = "W m⁻²"
-    #         return lwn
-    #     end
-
-    # # Compute SWN = SWD - SWU
-    # sim_var_dict["swn"] =
-    #     () -> begin
-    #         swd = sim_var_dict["swd"]()
-    #         swu = sim_var_dict["swu"]()
-    #         swn = swd - swu
-    #         swn.attributes["long_name"] = "Net shortwave radiation"
-    #         swn.attributes["units"] = "W m⁻²"
-    #         return swn
-    #     end
+    # Compute LWN = LWD - LWU
+    sim_var_dict["lwn"] =
+        () -> begin
+            lwd = sim_var_dict["lwd"]()
+            lwu = sim_var_dict["lwu"]()
+            lwn = lwd - lwu
+            lwn.attributes["long_name"] = "Net longwave radiation"
+            lwn.attributes["units"] = "W m⁻²"
+            lwn.attributes["start_date"] = lwu.attributes["start_date"]
+            return lwn
+        end
+
+    # Compute SWN = SWD - SWU
+    sim_var_dict["swn"] =
+        () -> begin
+            swd = sim_var_dict["swd"]()
+            swu = sim_var_dict["swu"]()
+            swn = swd - swu
+            swn.attributes["long_name"] = "Net shortwave radiation"
+            swn.attributes["units"] = "W m⁻²"
+            swn.attributes["start_date"] = swu.attributes["start_date"]
+            return swn
+        end
 
     return sim_var_dict
 end
 
 function get_obs_var_dict()
+    # contains monthly mslhf, msshf, msuwlwrf, msuwswrf
     era5_data_path = joinpath(
         ClimaLand.Artifacts.era5_monthly_averages_2008_folder_path(),
         "era5_monthly_surface_fluxes_200801-200812.nc",
     )
-    # era5_land_forcing_data_path = joinpath(
-    #     ClimaLand.Artifacts.era5_land_forcing_data2008_folder_path(
-    #         lowres = false,
-    #     ),
-    #     "era5_2008_0.9x1.25.nc",
-    # )
+    # contains hourly msdwlwrf, msdwswrf
+    era5_land_forcing_data_path = joinpath(
+        ClimaLand.Artifacts.era5_land_forcing_data2008_folder_path(
+            lowres = false,
+        ),
+        "era5_2008_1.0x1.0.nc",
+    )
 
     # Dict for loading in observational data
     obs_var_dict = Dict{String, Any}()
@@ -132,11 +136,11 @@ function get_obs_var_dict()
             )
 
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data =
                 -1 * cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
             # # relabel longitude values from [0, 359] to [-180, 179]
-            # obs_var.dims["longitude"] = obs_var.dims["longitude"] .- 180
+            # obs_var.dims[ClimaAnalysis.longitude_name(obs_var)] = obs_var.dims[ClimaAnalysis.longitude_name(obs_var)] .- 180
             obs_var_pos = ClimaAnalysis.OutputVar(
                 attribs,
                 obs_var.dims,
@@ -165,7 +169,7 @@ function get_obs_var_dict()
                 "start_date" => start_date,
             )
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data =
                 -1 * cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
             obs_var_pos = ClimaAnalysis.OutputVar(
@@ -190,7 +194,7 @@ function get_obs_var_dict()
             obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
 
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data = cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
 
             # Manually set attributes
@@ -222,7 +226,7 @@ function get_obs_var_dict()
             obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
 
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data = cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
             # Manually set attributes
             attribs = Dict(
@@ -245,15 +249,15 @@ function get_obs_var_dict()
     obs_var_dict["lwd"] =
         (start_date) -> begin
             obs_var = ClimaAnalysis.OutputVar(
-                era5_data_path,
+                era5_land_forcing_data_path,
                 "msdwlwrf",
                 new_start_date = start_date,
                 shift_by = Dates.firstdayofmonth,
             ) # units (W/m²)
             obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
 
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data = cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
             # Manually set attributes
             attribs = Dict(
@@ -276,15 +280,15 @@ function get_obs_var_dict()
     obs_var_dict["swd"] =
         (start_date) -> begin
             obs_var = ClimaAnalysis.OutputVar(
-                era5_data_path,
+                era5_land_forcing_data_path,
                 "msdwswrf",
                 new_start_date = start_date,
                 shift_by = Dates.firstdayofmonth,
             ) # units (W/m²)
             obs_var = ClimaAnalysis.replace(obs_var, missing => NaN)
 
             # Copy data at lon = 0 to lon = 360 to avoid white lines
-            push!(obs_var.dims["longitude"], 360)
+            push!(obs_var.dims[ClimaAnalysis.longitude_name(obs_var)], 360)
             new_data = cat(obs_var.data, obs_var.data[[1], :, :], dims = 1)
             # Manually set attributes
             attribs = Dict(
@@ -304,27 +308,81 @@ function get_obs_var_dict()
         end
 
 
-    # # Compute LWN = LWD - LWU
-    # obs_var_dict["lwn"] =
-    #     (start_date) -> begin
-    #         lwd = obs_var_dict["lwd"](start_date)
-    #         lwu = obs_var_dict["lwu"](start_date)
-    #         lwn = lwd - lwu
-    #         lwn.attributes["long_name"] = "Net longwave radiation"
-    #         lwn.attributes["units"] = "W m⁻²"
-    #         return lwu
-    #     end
-
-    # # Compute SWN = SWD - SWU
-    # obs_var_dict["swn"] =
-    #     (start_date) -> begin
-    #         swd = obs_var_dict["swd"](start_date)
-    #         swu = obs_var_dict["swu"](start_date)
-    #         swn = swd - swu
-    #         swn.attributes["long_name"] = "Net shortwave radiation"
-    #         swn.attributes["units"] = "W m⁻²"
-    #         return swn
-    #     end
+    # Compute LWN = LWD - LWU
+    obs_var_dict["lwn"] =
+        (start_date) -> begin
+            # LWD is very large because it's hourly, so it takes a while to load
+            lwd = obs_var_dict["lwd"](start_date)
+            lwu = obs_var_dict["lwu"](start_date)
+
+            monthly_times = lwu.dims["time"]
+            lwd_monthly_data = similar(lwu.data)
+            # average hourly LWD data to monthly so we can compare with monthly LWU
+            for (idx, time) in enumerate(monthly_times)
+                # collect all hourly time indices in this month
+                time_idxs = findall(x -> (x == time), lwd.dims["time"])
+
+                # average hourly data to monthly
+                lwd_monthly_data[:, :, idx] =
+                    mean(lwd.data[:, :, time_idxs], dims = 3)
+            end
+
+            # Compute LWN = LWD - LWU
+            lwn_data = lwd_monthly_data - lwu.data
+
+            attribs = Dict(
+                "short_name" => "lwn",
+                "long_name" => "Net longwave radiation",
+                "units" => "W m⁻²",
+                "start_date" => start_date,
+            )
+
+            lwn = ClimaAnalysis.OutputVar(
+                attribs,
+                lwu.dims,
+                lwu.dim_attributes,
+                lwn_data,
+            )
+            return lwn
+        end
+
+    # Compute SWN = SWD - SWU
+    obs_var_dict["swn"] =
+        (start_date) -> begin
+            # SWD is very large because it's hourly, so it takes a while to load
+            swd = obs_var_dict["swd"](start_date)
+            swu = obs_var_dict["swu"](start_date)
+
+            monthly_times = swu.dims["time"]
+            swd_monthly_data = similar(swu.data)
+            # average hourly SWD data to monthly so we can compare with monthly SWU
+            for (idx, time) in enumerate(monthly_times)
+                # collect all hourly time indices in this month
+                time_idxs = findall(x -> (x == time), swd.dims["time"])
+
+                # average hourly data to monthly
+                swd_monthly_data[:, :, idx] =
+                    mean(swd.data[:, :, time_idxs], dims = 3)
+            end
+
+            # Compute SWN = SWD - SWU
+            swn_data = swd_monthly_data - swu.data
+
+            attribs = Dict(
+                "short_name" => "swn",
+                "long_name" => "Net shortwave radiation",
+                "units" => "W m⁻²",
+                "start_date" => start_date,
+            )
+
+            swn = ClimaAnalysis.OutputVar(
+                attribs,
+                swu.dims,
+                swu.dim_attributes,
+                swn_data,
+            )
+            return swn
+        end
 
     return obs_var_dict
 end
diff --git a/experiments/long_runs/paper_plots.jl b/experiments/long_runs/paper_plots.jl
@@ -27,12 +27,14 @@ root_path = joinpath(pwd(), "snowy_land_longrun_gpu")
 outdir = "snowy_land_longrun_gpu-3785"
 root_path = outdir
 
-short_names = ["lhf", "shf", "lwu", "swu"]#, "lwn", "swn"]
+short_names = ["lhf", "shf", "lwn", "swn"]#, "lwu", "swu"]
 title_stubs = Dict(
     "lhf" => "LE",#"Latent heat flux",
     "shf" => "H",#"Sensible heat flux",
     "lwu" => "LWᵤ",#"Upward longwave radiation",
     "swu" => "SWᵤ",#"Upward shortwave radiation",
+    "lwn" => "LWₙ",#"Net longwave radiation",
+    "swn" => "SWₙ",#"Net shortwave radiation",
 )
 # Define levels for contour colorbars
 levels_dict = Dict(
@@ -92,8 +94,6 @@ function make_paper_figures(
 
         obs_var = obs_var_dict[short_name](sim_var.attributes["start_date"])
 
-
-
         sim_var_global_average = zeros(12)
         obs_var_global_average = zeros(12)
         for i in 1:12
@@ -510,7 +510,8 @@ function make_paper_figures(
     end
     save_name = joinpath(root_path, "combined_figures.png")
     save_name =
-        plot_bias ? joinpath(root_path, "combined_figures_bias.png") : save_name
+        plot_bias ? joinpath(root_path, "combined_figures_net-rad_bias.png") :
+        save_name
     save_name =
         plot_seasonal ? joinpath(root_path, "combined_figures_seasonal.png") :
         save_name
