Merge pull request #961 from ekluzek/lunadaylbugs

Fixes for the LUNA dayl bugs

Author: ekluzek

bld/namelist_files/namelist_defaults_ctsm.xml

@@ -393,8 +393,8 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <!-- The default filenames are given relative to the root directory
      for the CLM2 data in the CESM distribution -->
 <!-- Plant function types (relative to {csmdata}) -->
-<paramfile phys="clm5_0">lnd/clm2/paramdata/clm5_params.c200614.nc</paramfile>
-<paramfile phys="clm4_5">lnd/clm2/paramdata/clm45_params.c200614.nc</paramfile>
+<paramfile phys="clm5_0">lnd/clm2/paramdata/clm5_params.c200624.nc</paramfile>
+<paramfile phys="clm4_5">lnd/clm2/paramdata/clm45_params.c200624.nc</paramfile>
 
 <!-- ================================================================== -->
 <!-- FATES default parameter file                                       -->
@@ -442,7 +442,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <use_fertilizer use_crop=".true." phys="clm4_5">.true.</use_fertilizer>
 <use_fertilizer use_crop=".true." phys="clm5_0">.true.</use_fertilizer>
 
-<jmaxb1 use_luna=".true." phys="clm5_0">0.093563</jmaxb1>
+<jmaxb1 use_luna=".true." phys="clm5_0">0.17</jmaxb1>
 
 <use_grainproduct>.false.</use_grainproduct>
 <use_grainproduct use_crop=".true." phys="clm4_5">.false.</use_grainproduct> <!-- 1-year grain product pool default to off for clm45 if crop is turned on -->

doc/ChangeLog

@@ -1,4 +1,136 @@
 ===============================================================
+Tag name:  ctsm1.0.dev102
+Originator(s):  erik (Erik Kluzek,UCAR/TSS,303-497-1326)
+Date: Fri Jun 26 01:32:04 MDT 2020
+One-line Summary: Some important fixes for LUNA in clm5_0, and small urban issue in clm5_0
+
+Purpose of changes
+------------------
+
+Fix some important issues that were found in LUNA for the arctic. These changes were large enough
+that we had to re-tune the clm5_0 parameters to give a similar simulation. Leah Birch noticed 
+that the formulations in the code disagreed with the paper. Correcting this caused arctic plants
+to be less productive. In order to compensate parameters for leafcn, and slatop were adjusted for 
+broadleaf_deciduous_boreal_shrub and c3_arctic_grass for clm5_0.
+
+The total vegetation carbon threshold on exit_spinup (from accelerated decomposition (AD) mode) was 
+increased from 0.1 to 1.0. Without this change some PFT's can die out in exit spinup even when
+they have significant carbon stores from the AD spinup mode.
+
+There is an implicit assumption in the urban building energy model that building width equals 
+street width. However, this assumption can/should be relaxed and building width can be derived 
+from the morphology dataset. 
+
+
+Bugs fixed or introduced
+------------------------
+
+Issues fixed (include CTSM Issue #):  #738 #803 #953 #958 #1056
+ Fixes #1056 -- Remove allpfts dimension
+ Fixes #953 -- Incorrect formula in LUNA uses day time rather than average of day/night
+ Fixes #958 -- LUNA day length factor missing
+ Fixes #738 -- Totvegc threshold increased from 0.1 to 1.
+ Fixes #803 -- No longer assume building width
+
+Known bugs introduced in this tag (include github issue ID): #1060
+   #1060 -- Some code cleanup of luna bug fixes
+
+Significant changes to scientifically-supported configurations
+--------------------------------------------------------------
+
+Does this tag change answers significantly for any of the following physics configurations?
+(Details of any changes will be given in the "Answer changes" section below.)
+
+    [Put an [X] in the box for any configuration with significant answer changes.]
+
+[x] clm5_0
+
+[x] ctsm5_0-nwp
+
+[ ] clm4_5
+
+Notes of particular relevance for users
+---------------------------------------
+
+Caveats for users (e.g., need to interpolate initial conditions): clm5_0 climate adjusted for arctic plants
+
+Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none
+
+Changes made to namelist defaults (e.g., changed parameter values):
+
+Changes to the datasets (e.g., parameter, surface or initial files): New parameter files
+
+Substantial timing or memory changes: None
+
+Notes of particular relevance for developers: (including Code reviews and testing)
+---------------------------------------------
+NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide
+
+Caveats for developers (e.g., code that is duplicated that requires double maintenance): We 
+  chose some bad names for variables, so we'll clean this up in a seperate step (see #1060)
+
+Changes to tests or testing: None
+
+Code reviewed by: self, olyson, bill-sacks, wweider
+
+
+CTSM testing: regular
+
+ [PASS means all tests PASS and OK means tests PASS other than expected fails.]
+
+  build-namelist tests:
+
+    cheyenne - PASS (266 compare differently because of update in paramsfile and jmaxb1)
+
+  python testing (see instructions in python/README.md; document testing done):
+
+    cheyenne - PASS
+
+  regular tests (aux_clm):
+
+    cheyenne ---- PASS
+    izumi ------- PASS
+
+If the tag used for baseline comparisons was NOT the previous tag, note that here:
+
+
+Answer changes
+--------------
+
+Changes answers relative to baseline:  Yes for clm5_0
+
+  If a tag changes answers relative to baseline comparison the
+  following should be filled in (otherwise remove this section):
+
+  Summarize any changes to answers, i.e.,
+    - what code configurations: clm5_0 when Luna on as well as smaller changes for urban
+    - what platforms/compilers: all
+    - nature of change: climate changing, but tuned to respond similar to clm5.0
+
+   If this tag changes climate describe the run(s) done to evaluate the new
+   climate (put details of the simulations in the experiment database)
+     luna changes: oleson case clm50_ctsm10d089_2deg_GSWP3V1_luna3_jmaxb1-0.17_slatopA_leafcnA_hist
+
+   URL for LMWG diagnostics output used to validate new climate:
+
+    urban changes:
+
+
+http://webext.cgd.ucar.edu/I2000/clm50_ctsm10d098_1deg_GSWP3V1_CON_FIXBUILDENERGY_2000/lnd/clm50_ctsm10d098_1deg_GSWP3V1_CON_FIXBUILDENERGY_2000.1991_2010-clm50_ctsm10d098_1deg_GSWP3V1_CON_2000.1991_2010/setsIndex.html
+	
+
+Detailed list of changes
+------------------------
+
+List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): None
+
+Pull Requests that document the changes (include PR ids): #1034 and #962
+(https://github.com/ESCOMP/ctsm/pull)
+   #1034 -- Building energy fix
+   #962 --- Fixes for the LUNA dayl bugs 
+
+===============================================================
+===============================================================
 Tag name:  ctsm1.0.dev101
 Originator(s):  erik (Erik Kluzek,UCAR/TSS,303-497-1326), oleson (Keith Oleson)
 Date: Wed Jun 17 23:51:22 MDT 2020

doc/ChangeSum

@@ -1,5 +1,6 @@
 Tag                   Who      Date  Summary
 ============================================================================================================================
+  ctsm1.0.dev102 erik/ole 06/26/2020 Some important fixes for LUNA in clm5_0, and small urban issue in clm5_0
   ctsm1.0.dev101 ole/erik 06/17/2020 Changes from Keith to bring a list of variables to the parameter file
   ctsm1.0.dev100    sacks 06/09/2020 Deallocate memory after running init_interp
   ctsm1.0.dev099    sacks 06/08/2020 Update cime, needed for izumi machine updates

src/biogeophys/LunaMod.F90

@@ -319,7 +319,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
@@ -345,7 +347,7 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
          hourpd = dayl(g) / 3600._r8             
          tleafd10 = t_veg10_day(p) - tfrz
          tleafn10 = t_veg10_night(p) - tfrz
-         tleaf10  = (dayl(g)*tleafd10 +(86400._r8-dayl(g)) * tleafd10)/86400._r8 	     
+         tleaf10  = (dayl(g)*tleafd10 +(86400._r8-dayl(g)) * tleafn10)/86400._r8
          tair10 = t10(p)- tfrz
          relh10 = min(1.0_r8, rh10_p(p))  
 	 rb10v = rb10_p(p)	     
@@ -415,18 +417,20 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, &
                          PNcbold   = 0.0_r8                                     
                          call NitrogenAllocation(FNCa,forc_pbot10(p), relh10, CO2a10, O2a10, PARi10, PARimx10, rb10v, hourpd, &
                               tair10, tleafd10, tleafn10, &
-                              Jmaxb1, PNlcold, PNetold, PNrespold, &
-                              PNcbold, PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
+                              Jmaxb1, PNlcold, PNetold, PNrespold, PNcbold, dayl_factor(p), &
+                              PNstoreopt, PNlcopt, PNetopt, PNrespopt, PNcbopt)
                          vcmx25_opt= PNcbopt * FNCa * Fc25
                          jmx25_opt= PNetopt * FNCa * Fj25
 
                          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
 
@@ -485,8 +489,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 
@@ -803,7 +807,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, &
-     Jmaxb1, PNlcold, PNetold, PNrespold, PNcbold, &
+     Jmaxb1, 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)
@@ -823,6 +827,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 
@@ -904,7 +909,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(-params_inst%relhExp * max(relh10 - &
+  JmaxCoef = Jmaxb1 * dayl_factor * (1.0_r8 - exp(-params_inst%relhExp * max(relh10 - &
       params_inst%minrelh, 0.0_r8) / (1.0_r8 - params_inst%minrelh)))
   do while (PNlcoldi .NE. PNlc .and. jj < 100)      
      Fc = VcmxTKattge(tair10, tleafd10c) * Fc25

src/biogeophys/PhotosynthesisMod.F90

@@ -200,6 +200,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)
@@ -344,6 +346,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
@@ -888,6 +892,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', &

src/soilbiogeochem/SoilBiogeochemPrecisionControlMod.F90

@@ -43,7 +43,7 @@ subroutine SoilBiogeochemPrecisionControlInit( soilbiogeochem_carbonstate_inst,
     type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst
     !
     ! !LOCAL VARIABLES:
-    real(r8), parameter :: totvegcthresh = 0.1_r8   ! Total vegetation carbon threshold to zero out decomposition pools
+    real(r8), parameter :: totvegcthresh = 1.0_r8   ! Total vegetation carbon threshold to zero out decomposition pools
     !-----------------------------------------------------------------------
     ccrit    =  1.e-8_r8              ! critical carbon state value for truncation (gC/m2)
     ncrit    =  1.e-8_r8              ! critical nitrogen state value for truncation (gN/m2)