Adds Morel like MultiBandPhotosyntheticallyActiveRadiation

docs/src/model_components/light.md

@@ -2,21 +2,33 @@
 
 Nearly all BGC models require some model of the attenuation of PAR through the water. Usually this depends on the concentration of chlorophyll in the water (in phytoplankton), and may depend on the concentration of coloured dissolved organic matter or particulates.
 
-We currently have two models of light attenuation, a two band model by [Karleskind2011](@citet) and the more widely used three band model by [Morel1988](@citet). As the light level is diagnostic of the phytoplankton concentration these models are implemented with the light level as various auxiliary fields which are updated with callbacks within the biogeochemical model.
+We have two models implemented, a two band model by [Karleskind2011](@citet), and a more generic "multi band" model which can have the PAR split into arbitary many wavelength bands, but default to the widely used three band model by [Morel1988](@citet). As the light level is diagnostic of the phytoplankton concentration these models are implemented with the light level as various auxiliary fields which are updated within the biogeochemical model.
 
 Models requiring light attenuation models will set these up automatically, for example [LOBSTER](@ref LOBSTER) sets `light_attenuation_model = TwoBandPhotosyntheticallyActiveRadiation()`. You may choose others. Additionally, you can pass the surface PAR as a function of horizontal position and time. The default for LOBSTER is `(x, y, t) -> 100*max(0.0, cos(t*π/(12hours)))`.
 
-## Model equations (for the two band model)
+## The multi band model
+The surface intensity is split into multiple bands (usually with equal weight, but users may specify custom weights), and the attenuation of each band (i) is computed from the radiative transfer equation:
+```math
+\frac{\partial PAR^i}{\partial z} = PAR\^i (k^w(i) + \chi(i)Chl^{e(i)}),
+```
+where ``Chl`` is the concentration of chlorophyll, ``k^w(i)`` is the band specific water attenuation coefficient, ``\chi(i)`` the chlorophyll attenuation coefficient, and ``e(i)`` the chlorophyll exponent.
 
-Light attenuation is calculated by integrating attenuation (from the surface). The $PAR$ is considered as two components attenuated at different rates. At depth $z$ the total $PAR$ is given by:
+The water concentration of chlorophyll is returned by a function `chlorophyll` with arguments `biogeochemistry` and `model`. For the `LOBSTER` model this returns a constant ratio of the phytoplankton concentration, but may be different for other models.
 
-$PAR = \frac{PAR_0}{2} \left[\exp\left(k_rz + \chi_r\int_{z=0}^z Chl_r dz\right) + \exp\left(k_bz + \chi_b\int_{z=0}^z Chl_b dz\right)\right],$
+## The two band model
+Light attenuation is calculated by integrating attenuation (from the surface). The ``PAR`` is considered as two components attenuated at different rates. At depth $z$ the total $PAR$ is given by:
 
-where $PAR_0$ is the surface value, $k_r$ and $k_b$ are the red and blue attenuation coefficients of water, $\chi_r$ and $\chi_b$ are the red and blue chlorophyll attenuation coefficients, and $Chl_r$ and $Chl_b$ are the red and blue chlorophyll pigment concentrations. The chlorophyll pigment concentration is derived from the phytoplankton concentration where it is assumed that the pigment concentration is given by:
+```math
+PAR = \frac{PAR_0}{2} \left[\exp\left(k_rz + \chi_r\int_{z=0}^z Chl_r dz\right) + \exp\left(k_bz + \chi_b\int_{z=0}^z Chl_b dz\right)\right],
+```
 
-$Chl = PR_{Chl:P},$
+where ``PAR_0`` is the surface value, ``k_r`` and ``k_b`` are the red and blue attenuation coefficients of water, ``\chi_r`` and ``\chi_b`` are the red and blue chlorophyll attenuation coefficients, and ``Chl_r`` and ``Chl_b`` are the red and blue chlorophyll pigment concentrations. The chlorophyll pigment concentration is derived from the phytoplankton concentration where it is assumed that the pigment concentration is given by:
 
-where the ratio is constant and found in [Parameters](@ref parameters). The red and blue pigment concentrations are then found as $Chl_r = \left(\frac{Chl}{r_\text{pig}}\right)^{e_r}$ and $Chl_b = \left(\frac{Chl}{r_\text{pig}}\right)^{e_b}$. 
+```math
+Chl = PR_{Chl:P},
+```
+
+where the ratio is constant and found in [Parameters](@ref parameters). The red and blue pigment concentrations are then found as ``Chl_r = \left(\frac{Chl}{r_\text{pig}}\right)^{e_r}`` and ``Chl_b = \left(\frac{Chl}{r_\text{pig}}\right)^{e_b}``. 
 
 ### Parameter variable names
 

src/Light/2band.jl

@@ -83,7 +83,7 @@ end
                                                chlorophyll_blue_exponent::FT = 0.674,
                                                pigment_ratio::FT = 0.7,
                                                phytoplankton_chlorophyll_ratio::FT = 1.31,
-                                               surface_PAR::SPAR = (x, y, t) -> 100 * max(0.0, cos(t * π / 12hours)))
+                                               surface_PAR::SPAR = default_surface_PAR)
 
 Keyword Arguments
 ==================
@@ -103,7 +103,7 @@ function TwoBandPhotosyntheticallyActiveRadiation(; grid,
                                                     chlorophyll_blue_exponent::FT = 0.674,
                                                     pigment_ratio::FT = 0.7,
                                                     phytoplankton_chlorophyll_ratio::FT = 1.31,
-                                                    surface_PAR::SPAR = (x, y, t) -> 100 * max(0.0, cos(t * π / 12hours))) where {FT, SPAR} # mgChl/mol N
+                                                    surface_PAR::SPAR = default_surface_PAR) where {FT, SPAR} # mgChl/mol N
 
     field = CenterField(grid; boundary_conditions = 
                             regularize_field_boundary_conditions(

src/Light/Light.jl

@@ -1,12 +1,16 @@
 """
-Light attenuation by chlorophyll as described by [Karleskind2011](@citet) (implemented as twoBand) and [Morel1988](@citet).
+Light attenuation by chlorophyll as described by [Karleskind2011](@citet) (implemented as TwoBand) and [Morel1988](@citet) (as MultiBand).
 """
 module Light
 
-export TwoBandPhotosyntheticallyActiveRadiation, update_PAR!, PrescribedPhotosyntheticallyActiveRadiation
+export TwoBandPhotosyntheticallyActiveRadiation, 
+       PrescribedPhotosyntheticallyActiveRadiation, 
+       MultiBandPhotosyntheticallyActiveRadiation
+
+using Adapt
 
 using KernelAbstractions, Oceananigans.Units
-using Oceananigans.Architectures: device, architecture
+using Oceananigans.Architectures: device, architecture, on_architecture
 using Oceananigans.Utils: launch!
 using Oceananigans: Center, Face, fields
 using Oceananigans.Grids: node, znodes
@@ -20,6 +24,7 @@ using Oceananigans.BoundaryConditions: fill_halo_regions!,
                                        domain_boundary_indices,
                                        RightBoundary
 using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid
+using OceanBioME: chlorophyll
 
 import Adapt: adapt_structure, adapt
 import Base: show, summary
@@ -28,7 +33,7 @@ import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields, update_bio
 import Oceananigans.BoundaryConditions: _fill_top_halo!
 
 include("2band.jl")
-include("morel.jl")
+include("multi_band.jl")
 include("prescribed.jl")
 
 default_surface_PAR(x, y, t) = default_surface_PAR(t)

src/Light/morel_coefficients.jl

@@ -0,0 +1,31 @@
+const MOREL_λ = [350, 355, 360, 365, 370, 375, 380, 385, 390, 395, 
+                 400, 405, 410, 415, 420, 425, 430, 435, 440, 445, 
+                 450, 455, 460, 465, 470, 475, 480, 485, 490, 495, 
+                 500, 505, 510, 515, 520, 525, 530, 535, 540, 545, 
+                 550, 555, 560, 565, 570, 575, 580, 585, 590, 595, 
+                 600, 605, 610, 615, 620, 625, 630, 635, 640, 645, 
+                 650, 655, 660, 665, 670, 675, 680, 685, 690, 695, 700]
+
+const MOREL_kʷ = [0.0271 , 0.0238 , 0.0216 , 0.0188 , 0.0177 , 0.01595, 0.0151 , 0.01376, 0.01271, 0.01208,
+                  0.01042, 0.0089 , 0.00812, 0.00765, 0.00758, 0.00768, 0.0077 , 0.00792, 0.00885, 0.0099 ,
+                  0.01148, 0.01182, 0.01188, 0.01211, 0.01251, 0.0132 , 0.01444, 0.01526, 0.0166 , 0.01885,
+                  0.02188, 0.02701, 0.03385, 0.0409 , 0.04214, 0.04287, 0.04454, 0.0463 , 0.04846, 0.05212,
+                  0.05746, 0.06053, 0.0628 , 0.06507, 0.07034, 0.07801, 0.09038, 0.11076, 0.13584, 0.16792,
+                  0.2231 , 0.25838, 0.26506, 0.26843, 0.27612, 0.284  , 0.29218, 0.30176, 0.31134, 0.32553,
+                  0.34052, 0.3715 , 0.41048, 0.42947, 0.43946, 0.44844, 0.46543, 0.48642, 0.5164 , 0.55939, 0.62438]
+
+const MOREL_e = [0.778  , 0.767  , 0.756  , 0.737  , 0.72   , 0.7    , 0.685  , 0.673  , 0.67   , 0.66   ,
+                 0.64358, 0.64776, 0.65175, 0.65555, 0.65917, 0.66259, 0.66583, 0.66889, 0.67175, 0.67443,
+                 0.67692, 0.67923, 0.68134, 0.68327, 0.68501, 0.68657, 0.68794, 0.68903, 0.68955, 0.68947,
+                 0.6888 , 0.68753, 0.68567, 0.6832 , 0.68015, 0.67649, 0.67224, 0.66739, 0.66195, 0.65591,
+                 0.64927, 0.64204, 0.64   , 0.63   , 0.623  , 0.615  , 0.61   , 0.614  , 0.618  , 0.622  ,
+                 0.626  , 0.63   , 0.634  , 0.638  , 0.642  , 0.647  , 0.653  , 0.658  , 0.663  , 0.667  ,
+                 0.672  , 0.677  , 0.682  , 0.687  , 0.695  , 0.697  , 0.693  , 0.665  , 0.64   , 0.62   , 0.6    ]
+
+const MOREL_χ = [0.153  , 0.149  , 0.144  , 0.14   , 0.136  , 0.131  , 0.127  , 0.123  , 0.119  , 0.118  ,
+                 0.11748, 0.12066, 0.12259, 0.12326, 0.12269, 0.12086, 0.11779, 0.11372, 0.10963, 0.1056 ,
+                 0.10165, 0.09776, 0.09393, 0.09018, 0.08649, 0.08287, 0.07932, 0.07584, 0.07242, 0.06907,
+                 0.06579, 0.06257, 0.05943, 0.05635, 0.05341, 0.05072, 0.04829, 0.04611, 0.04419, 0.04253,
+                 0.04111, 0.03996, 0.039  , 0.0375 , 0.036  , 0.034  , 0.033  , 0.0328 , 0.0325 , 0.033  ,
+                 0.034  , 0.035  , 0.036  , 0.0375 , 0.0385 , 0.04   , 0.042  , 0.043  , 0.044  , 0.0445 ,
+                 0.045  , 0.046  , 0.0475 , 0.049  , 0.0515 , 0.052  , 0.0505 , 0.044  , 0.039  , 0.034  , 0.03   ]