LUNA day length factor

[wwieder]

Brief summary of bug

@lmbirch89 thinks this is a bug, based on the orig. LUNA paper.
@dlawrenncar agrees, but there are big answer changes that should be considered before bringing them into CLM5 physics

General bug information

CTSM version you are using: [output of git describe]
CLM5 release

Does this bug cause significantly incorrect results in the model's science? [Yes / No]
Yes

Configurations affected: [Fill this in if known.]
default CLM5

Details of bug

From Leah's paper:
"The function, f(daylegth), is a scaling factor that is based on the formulation in Bauerle et al. (2012). However, the denominator is set to 12, when it should be the maximum day length possible at a particular latitude (Bauerle et al., 2012). This may not be appropriate where some regions experience up to 24 hour day light in summer, which allows Equation 1 to be > 1, particularly around the summer solstice in June."

Important details of your setup / configuration so we can reproduce the bug

CTSM/src/biogeophys/LunaMod.F90

Line 900 in 384c726

JmaxCoef = Jmaxb1 * ((hourpd / 12.0_r8)**2.0_r8) * (1.0_r8 - exp(-relhExp * max(relh10 - minrelh, 0.0_r8) / &

Important output or errors that show the problem

[ekluzek]

The changes from @lmbirch89 are:

diff --git a/src/biogeophys/LunaMod.F90 b/src/biogeophys/LunaMod.F90
index 4947d1c1..4b19e64c 100644
--- a/src/biogeophys/LunaMod.F90
+++ b/src/biogeophys/LunaMod.F90
@@ -402,7 +402,7 @@ module LunaMod
                          PNcbold   = 0.0_r8                                     
                          call NitrogenAllocation(FNCa,forc_pbot10(p), relh10, CO2a10, O2a10, PARi10, PARimx10, rb10v, hourpd, &
                               tair10, tleafd10, tleafn10, &
-                              Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp, PNlcold, PNetold, PNrespold, &
+                              Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp, PNlcold, PNetold, PNrespold, dayl_factor(p), &
                               PNcbold, PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
                          vcmx25_opt= PNcbopt * FNCa * Fc25
                          jmx25_opt= PNetopt * FNCa * Fj25
@@ -791,7 +791,7 @@ end subroutine Clear24_Climate_LUNA
 !Use the LUNA model to calculate the Nitrogen partioning 
 subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PARimx10,rb10, hourpd, tair10, tleafd10, tleafn10, &
      Jmaxb0, Jmaxb1, Wc2Wjb0, relhExp,&
-     PNlcold, PNetold, PNrespold, PNcbold, &
+     PNlcold, PNetold, PNrespold, PNcbold, dayl_factor &
      PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
   implicit none
   real(r8), intent (in) :: FNCa                       !Area based functional nitrogen content (g N/m2 leaf)
@@ -814,6 +814,7 @@ subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PAR
   real(r8), intent (in) :: PNetold                    !old value of the proportion of nitrogen allocated to electron transport (unitless)
   real(r8), intent (in) :: PNrespold                  !old value of the proportion of nitrogen allocated to respiration (unitless)
   real(r8), intent (in) :: PNcbold                    !old value of the proportion of nitrogen allocated to carboxylation (unitless)  
+  real(r8), intent (in) :: dayl_factor                !daylight scale factor
   real(r8), intent (out):: PNstoreopt                 !optimal proportion of nitrogen for storage 
   real(r8), intent (out):: PNlcopt                    !optimal proportion of nitrogen for light capture 
   real(r8), intent (out):: PNetopt                    !optimal proportion of nitrogen for electron transport 
@@ -897,7 +898,7 @@ subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PAR
   tleafd10c = min(max(tleafd10, Trange1), Trange2)    !constrain the physiological range
   tleafn10c = min(max(tleafn10, Trange1), Trange2)    !constrain the physiological range
   ci = 0.7_r8 * CO2a10 
-  JmaxCoef = Jmaxb1 * ((hourpd / 12.0_r8)**2.0_r8) * (1.0_r8 - exp(-relhExp * max(relh10 - minrelh, 0.0_r8) / &
+  JmaxCoef = Jmaxb1 * dayl_factor * (1.0_r8 - exp(-relhExp * max(relh10 - minrelh, 0.0_r8) / &
        (1.0_r8 - minrelh)))
   do while (PNlcoldi .NE. PNlc .and. jj < 100)      
      Fc = VcmxTKattge(tair10, tleafd10c) * Fc25

Note, that dayl_ is set in CanopyFluxes as follows...

         ! calculate dayl_factor as the ratio of (current:max dayl)^2
         ! set a minimum of 0.01 (1%) for the dayl_factor
         dayl_factor(p)=min(1._r8,max(0.01_r8,(dayl(g)*dayl(g))/(max_dayl(g)*max_dayl(g))))

hourpd in Luna is dayl / 3600.0. So algebraically the two terms are equivalent when max_dayl is assumed to be 12 hours. The only difference then are the constraints of making sure it's between 0.01 and 1.

[ekluzek]

I checked the code history and this problem goes back to the very first version of Luna we had from 2015.