wip - running shortwave and longwave solver with developing interface

Author: bischtob

src/Radiation/gray_breezy.jl

@@ -1,8 +1,7 @@
 using RRTMGP.Optics: GrayOpticalThicknessSchneider2004
-using RRTMGP.Parameters: RRTMGPParameters
 using RRTMGP.AtmosphericStates: GrayAtmosphericState
-using RRTMGP.RTE: TwoStreamLWRTE
-using RRTMGP.RTESolver: solve_lw!
+using RRTMGP.RTE: TwoStreamLWRTE, TwoStreamSWRTE
+using RRTMGP.RTESolver: solve_lw!, solve_sw!
 
 using Oceananigans
 using Oceananigans: field
@@ -17,14 +16,26 @@ DA = array_type(architecture)
 Nx, Ny, Nz = 1, 1, 64
 Lx, Ly, Lz = 1, 20_000_000, 1_000 # domain size in meters
 Nbnd, Ngpt = 1, 1 # gray model
+deg2rad = FT(π / 180)
+
 lw_inc_flux = nothing # no incoming longwave flux
-sfc_emissivity = 1
+sw_inc_flux = FT(1407.679) # no incoming shortwave flux
+sw_inc_flux_diffuse = nothing
+zenith_angle = FT(52.95) # zenith angle in degrees
+sfc_emissivity = FT(1) # surface emissivity
+albedo_direct = FT(0.1) # surface albedo (direct)
+albedo_diffuse = FT(0.1) # surface albedo (diffuse)
+
 planet_radius = 6_371_000
 p_surface = 100_000 # surface pressure (Pa)
 p_top = 9_000 # top of atmosphere pressure / emission level (Pa)
 lat_center = 0
 optical_properties = GrayOpticalThicknessSchneider2004(FT)
-
+#     cos_zenith .= cos(deg2rad * 52.95) # corresponding to ~52.95 deg zenith angle    
+#     toa_flux .= FT(1407.679)
+#     sfc_alb_direct .= FT(0.1)
+#     sfc_alb_diffuse .= FT(0.1)
+#     inc_flux_diffuse = nothing
 # Grid setup, we also need to say where on the planet our box is located
 grid = RectilinearGrid(
     architecture,
@@ -55,10 +66,28 @@ atmospheric_state = GrayAtmosphericState(
 
 # Set up radiation model
 sfc_emission = DA{FT}(undef, Nbnd, Nx, Ny)
+sfc_alb_direct = DA{FT}(undef, Nbnd, Nx, Ny)
+sfc_alb_diffuse = DA{FT}(undef, Nbnd, Nx, Ny)
+cos_zenith = DA{FT}(undef, Nx, Ny)
+toa_flux = DA{FT}(undef, Nx, Ny)
+lw_toa_inc_flux = nothing
+inc_flux_diffuse = nothing
 fill!(sfc_emission, FT(sfc_emissivity))
+fill!(sfc_alb_direct, FT(albedo_direct))
+fill!(sfc_alb_diffuse, FT(albedo_diffuse))
+fill!(cos_zenith, FT(cos(deg2rad * zenith_angle)))
+fill!(toa_flux, FT(sw_inc_flux))
 SLVLW = TwoStreamLWRTE
-slv_lw = SLVLW(grid; sfc_emission, lw_inc_flux)
+SLVSW = TwoStreamSWRTE
+lw_params = (; sfc_emission, lw_inc_flux)
+sw_params = (; cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse)
+slv_lw = SLVLW(grid; lw_params...)
+slv_sw = SLVSW(grid; sw_params...)
 
-# solve the radiation model
-solve_lw!(slv_lw, atmospheric_state)
+# Solve the LW radiation model
+function update_radative_fluxes!(slv_lw, slv_sw, atmospheric_state)
+    solve_lw!(slv_lw, atmospheric_state)
+    solve_sw!(slv_sw, atmospheric_state)
+end
 
+update_radative_fluxes!(slv_lw, slv_sw, atmospheric_state)

src/Radiation/rrtmgp_interface.jl

@@ -95,15 +95,16 @@ function RRTMGP.RTE.TwoStreamLWRTE(
     FT = eltype(grid)
     Nx, Ny = grid.Nx, grid.Ny
     Ncols = Nx * Ny
+    Nbnd = size(sfc_emission, 1)
     grid_params = RRTMGPGridParams(grid)
     params = RRTMGPParameters(eltype(grid))
     (; context) = grid_params
     op = TwoStream(grid_params)
     src = SourceLW2Str(grid_params; params)
-    if lw_inc_flux !== nothing
-        lw_inc_flux = reshape(lw_inc_flux, 1, Ncols)
-    end
-    bcs = LwBCs(reshape(sfc_emission, 1, Ncols), lw_inc_flux)
+    bcs = LwBCs(
+        reshape(sfc_emission, Nbnd, Ncols), 
+        lw_inc_flux isa Nothing ? nothing : reshape(lw_inc_flux, 1, Ncols)
+    )
     fluxb = FluxLW(grid_params)
     flux = FluxLW(grid_params)
     return TwoStreamLWRTE(context, op, src, bcs, fluxb, flux)
@@ -134,11 +135,20 @@ function RRTMGP.RTE.TwoStreamSWRTE(
     inc_flux_diffuse,
     sfc_alb_diffuse,
 )
+    Nx, Ny = grid.Nx, grid.Ny
+    Ncols = Nx * Ny
+    Nbnd = size(sfc_alb_direct, 1)
     grid_params = RRTMGPGridParams(grid)
     (; context) = grid_params
     op = TwoStream(grid_params)
     src = SourceSW2Str(grid_params)
-    bcs = SwBCs(cos_zenith, toa_flux, sfc_alb_direct, inc_flux_diffuse, sfc_alb_diffuse)
+    bcs = SwBCs(
+        reshape(cos_zenith, Ncols),
+        reshape(toa_flux, Ncols),
+        reshape(sfc_alb_direct, Nbnd, Ncols), 
+        inc_flux_diffuse isa Nothing ? nothing : reshape(inc_flux_diffuse, Nbnd, Ncols), 
+        reshape(sfc_alb_diffuse, Nbnd, Ncols),
+    )
     fluxb = FluxSW(grid_params)
     flux = FluxSW(grid_params)
     return TwoStreamSWRTE(context, op, src, bcs, fluxb, flux)