Luna Bug #3, hard coded initial conditions

[wwieder]

Along with #958 and #953, this seems to be a bug in the LUNA code where initial conditions at the start of the growing season for Vcmax and Jmax are hard coded to the same value, globally

CTSM/src/biogeophys/LunaMod.F90

Line 475 in c5a9030

jmx25_z(p, z) = 85._r8

Extensive discussion can be found in #947, where @lmbirch89 has suggested using values from the last day of the previous year so that there's some climatological memory of appropriate LUNA parameters at the start of the growing season.

It likely makes sense to bring this into #961, where previous bugs were addressed to evaluate the impact on global simulations

[ekluzek]

I've lifted these changes from @lmbirch89 PR. And I've changed it so that it initializes the values to the previous hard-coded values for every year on startup. After running it saves the last values to the restart file and uses it.

Here are the changes to do that:

diff --git a/src/biogeophys/LunaMod.F90 b/src/biogeophys/LunaMod.F90
index 9087d1e1..6bf6010d 100644
--- a/src/biogeophys/LunaMod.F90
+++ b/src/biogeophys/LunaMod.F90
@@ -306,7 +306,9 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
     vcmx25_z      => photosyns_inst%vcmx25_z_patch                    , & ! Output: [real(r8) (:,:) ] patch leaf Vc,max25 (umol/m2 leaf/s) for canopy layer 
     jmx25_z       => photosyns_inst%jmx25_z_patch                     , & ! Output: [real(r8) (:,:) ] patch leaf Jmax25 (umol electron/m**2/s) for canopy layer
     pnlc_z        => photosyns_inst%pnlc_z_patch                      , & ! Output: [real(r8) (:,:) ] patch proportion of leaf nitrogen allocated for light capture for canopy layer 
-    enzs_z        => photosyns_inst%enzs_z_patch                        & ! Output: [real(r8) (:,:) ] enzyme decay status 1.0-fully active; 0-all decayed during stress
+    enzs_z        => photosyns_inst%enzs_z_patch                      , & ! Output: [real(r8) (:,:) ] enzyme decay status 1.0-fully active; 0-all decayed during stress
+    vcmx_prevyr   => photosyns_inst%vcmx_prevyr                       , & ! Output: [real(r8) (:,:) ] patch leaf Vc,max25 from previous year avg
+    jmx_prevyr    => photosyns_inst%jmx_prevyr                          & ! Output: [real(r8) (:,:) ] patch leaf Jmax25 from previous year avg
     )  
     !----------------------------------------------------------------------------------------------------------------------------------------------------------
     !set timestep
@@ -410,10 +412,12 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
                          chg = vcmx25_opt-vcmx25_z(p, z)
                          chg_constrn = min(abs(chg),vcmx25_z(p, z)*max_daily_pchg)
                          vcmx25_z(p, z)  = vcmx25_z(p, z)+sign(1.0_r8,chg)*chg_constrn
+                         vcmx_prevyr(p,z) = vcmx25_z(p,z)
                           
                          chg = jmx25_opt-jmx25_z(p, z)
                          chg_constrn = min(abs(chg),jmx25_z(p, z)*max_daily_pchg)
                          jmx25_z(p, z)  = jmx25_z(p, z)+sign(1.0_r8,chg)*chg_constrn 
+                         jmx_prevyr(p,z) = jmx25_z(p,z)
 
                          PNlc_z(p, z)= PNlcopt
 
@@ -472,8 +476,8 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
                 endif !if not C3 plants                   
          else
             do z = 1 , nrad(p)
-               jmx25_z(p, z) = 85._r8
-               vcmx25_z(p, z) = 50._r8
+               jmx25_z(p, z) = jmx_prevyr(p,z)
+               vcmx25_z(p, z) = vcmx_prevyr(p,z)
             end do
          endif !checking for LAI and LNC
      endif !the first daycheck 
diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90
index a111cab1..a0710358 100644
--- a/src/biogeophys/PhotosynthesisMod.F90
+++ b/src/biogeophys/PhotosynthesisMod.F90
@@ -183,6 +183,8 @@ module  PhotosynthesisMod
      ! LUNA specific variables
      real(r8), pointer, public  :: vcmx25_z_patch    (:,:) ! patch  leaf Vc,max25 (umol CO2/m**2/s) for canopy layer 
      real(r8), pointer, public  :: jmx25_z_patch     (:,:) ! patch  leaf Jmax25 (umol electron/m**2/s) for canopy layer 
+     real(r8), pointer, public  :: vcmx_prevyr       (:,:) ! patch  leaf Vc,max25 previous year running avg
+     real(r8), pointer, public  :: jmx_prevyr       (:,:) ! patch  leaf Jmax25 previous year running avg
      real(r8), pointer, public  :: pnlc_z_patch      (:,:) ! patch proportion of leaf nitrogen allocated for light capture for canopy layer
      real(r8), pointer, public  :: enzs_z_patch      (:,:) ! enzyme decay status 1.0-fully active; 0-all decayed during stress
      real(r8), pointer, public  :: fpsn24_patch      (:)   ! 24 hour mean patch photosynthesis (umol CO2/m**2 ground/day)
@@ -328,6 +330,8 @@ subroutine InitAllocate(this, bounds)
       ! statements.
       allocate(this%vcmx25_z_patch  (begp:endp,1:nlevcan)) ; this%vcmx25_z_patch    (:,:) = 30._r8
       allocate(this%jmx25_z_patch   (begp:endp,1:nlevcan)) ; this%jmx25_z_patch     (:,:) = 60._r8 
+      allocate(this%vcmx_prevyr     (begp:endp,1:nlevcan)) ; this%vcmx_prevyr       (:,:) = 85._r8
+      allocate(this%jmx_prevyr      (begp:endp,1:nlevcan)) ; this%jmx_prevyr        (:,:) = 50._r8
       allocate(this%pnlc_z_patch    (begp:endp,1:nlevcan)) ; this%pnlc_z_patch      (:,:) = 0.01_r8
       allocate(this%fpsn24_patch    (begp:endp))           ; this%fpsn24_patch      (:)   = nan
       allocate(this%enzs_z_patch    (begp:endp,1:nlevcan)) ; this%enzs_z_patch      (:,:) = 1._r8
@@ -833,6 +837,14 @@ subroutine Restart(this, bounds, ncid, flag)
          dim1name='pft', dim2name='levcan', switchdim=.true., &
          long_name='Maximum carboxylation rate at 25 celcius for canopy layers', units='umol CO2/m**2/s', &
          interpinic_flag='interp', readvar=readvar, data=this%jmx25_z_patch)
+      call restartvar(ncid=ncid, flag=flag, varname='vcmx_prevyr', xtype=ncd_double,  &
+         dim1name='pft', dim2name='levcan', switchdim=.true., &
+         long_name='avg carboxylation rate at 25 celsius for canopy layers', units='umol CO2/m**2/s', &
+         interpinic_flag='interp', readvar=readvar, data=this%vcmx_prevyr)
+      call restartvar(ncid=ncid, flag=flag, varname='jmx_prevyr', xtype=ncd_double,  &
+         dim1name='pft', dim2name='levcan', switchdim=.true., &
+         long_name='avg carboxylation rate at 25 celsius for canopy layers', units='umol CO2/m**2/s', &
+         interpinic_flag='interp', readvar=readvar, data=this%jmx_prevyr)
       call restartvar(ncid=ncid, flag=flag, varname='pnlc_z', xtype=ncd_double,  &
          dim1name='pft', dim2name='levcan', switchdim=.true., &
          long_name='proportion of leaf nitrogen allocated for light capture', units='unitless', &