Photosynthesis Refactors

Merge branch 'rgknox-more-photo-refactors'

Plant hydraulics (when it is one day merged in, not today) will enable
photosynthesis to experience hydraulic limitations at the cohort
scale. Currently photosynthesis is calculated at the leaf-layer and
pft scale, but we need to have code that can also handle unique
calculations on every cohort's leaf layer. These refactors of
photosynthesis enable this optioning. In order to do this, a good deal
of code modularization was implemented. In particular, the various
processes that are calculated at the leaf layer, such as the
calculation of specific biophysical rates (such as vcmax, jmax, tpu,
etc), photosynthesis and leaf maintenance respiration where all sent
to their own subroutines. Some efforts were also made to improve
variable naming so they are more descriptive, and to add some
explanatory comments.

The changes do impact things like order of operations. The changes
should not create climate changing results. They should not be b4b
either.

Single site 50 year regression tests at the 1x1 brazil site indicate
that there are no discernible differences in gpp, npp, ar and canopy
structure. See plots below.  The standard ed, clm_short_45 and
clm_short_50 were also compared against baselines, with the
expectation that differences would be present, but small. Results
indicate that is true.

User interface changes?: No

Code review:

Testing:
  rgknox:
    Test suite: edTest,clm_short_45,clm_short_50: lawrencium-lr3, intel
    Test baseline: 5c5928f (b4b on ed tests not expected pass)
    Test namelist changes: none
    Test answer changes: b4b on clm_short_45 and clm_short_50
    Test summary: all Pass

  andre:
    Test suite: ed - yellowstone intel, pgi, gnu
                     hobart nag
    Test baseline: d116008
    Test namelist changes: none
    Test answer changes: yes for ED compsets, bit for bit for other clm
    Test summary: all functionality pass, bit for bit as above

    Test suite: clm_short - yellowstone intel, pgi, gnu
    Test baseline: clm4_5_12_r195
    Test namelist changes: none
    Test answer changes: yes for ED compsets, bit for bit for other clm
    Test summary: all tests pass

components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90

@@ -246,7 +246,6 @@ subroutine allocate_live_biomass(cc_p,mode)
        currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
             currentcohort%laimemory)*leaf_frac
 
-
     else ! Leaves are on (leaves_off_switch==1)
 
     !the purpose of this section is to figure out the root and stem biomass when the leaves are off
@@ -1038,13 +1037,15 @@ subroutine copy_cohort( currentCohort,copyc )
     n%gpp_acc_hold    = o%gpp_acc_hold
     n%gpp_acc         = o%gpp_acc
     n%gpp_tstep       = o%gpp_tstep
+
     n%npp_acc_hold    = o%npp_acc_hold
     n%npp_tstep       = o%npp_tstep
 
     if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ia ',o%npp_acc
     if ( DEBUG ) write(fates_log(),*) 'EDcohortDyn Ib ',o%resp_acc
 
-    n%npp_acc_hold    = o%npp_acc_hold
+    n%npp_acc         = o%npp_acc
+
     n%resp_tstep      = o%resp_tstep
     n%resp_acc        = o%resp_acc
     n%resp_acc_hold   = o%resp_acc_hold

components/clm/src/ED/main/EDTypesMod.F90

@@ -149,6 +149,16 @@ module EDTypesMod
   ! HLM will interpret that the value should not be included in the average.
   real(r8) :: cp_hio_ignore_val
 
+
+
+  ! Module switches (this will be read in one day)
+  ! This variable only exists now to serve as a place holder
+  !!!!!!!!!! THIS SHOULD NOT BE SET TO TRUE !!!!!!!!!!!!!!!!!
+  logical,parameter :: use_fates_plant_hydro = .false.
+
+
+
+
   !************************************
   !** COHORT type structure          **
   !************************************
@@ -660,4 +670,10 @@ subroutine set_root_fraction( this , depth_gl)
 
   end subroutine set_root_fraction
 
+
+  ! =====================================================================================
+
+
+
+
 end module EDTypesMod

components/clm/src/ED/main/FatesInterfaceMod.F90

@@ -191,14 +191,13 @@ module FatesInterfaceMod
       ! Shaded canopy resistance [s/m]
       real(r8), allocatable :: rssha_pa(:)
 
-      ! Canopy conductance [mmol m-2 s-1]
-      real(r8), allocatable :: gccanopy_pa(:)
+      ! leaf photosynthesis (umol CO2 /m**2/ s)
+      ! (NOT CURRENTLY USED, PLACE-HOLDER)
+      !real(r8), allocatable :: psncanopy_pa(:)
 
-      ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s)
-      real(r8), allocatable :: psncanopy_pa(:)
-
-      ! patch sunlit leaf maintenance respiration rate (umol CO2/m**2/s) 
-      real(r8), allocatable :: lmrcanopy_pa(:)
+      ! leaf maintenance respiration rate (umol CO2/m**2/s) 
+      ! (NOT CURRENTLY USED, PLACE-HOLDER)
+      !real(r8), allocatable :: lmrcanopy_pa(:)
 
       ! Canopy Radiation Boundaries
       ! ---------------------------------------------------------------------------------
@@ -404,11 +403,9 @@ subroutine allocate_bcout(bc_out)
       allocate(bc_out%btran_pa(maxPatchesPerCol))
 
       ! Photosynthesis
+
       allocate(bc_out%rssun_pa(maxPatchesPerCol))
       allocate(bc_out%rssha_pa(maxPatchesPerCol))
-      allocate(bc_out%gccanopy_pa(maxPatchesPerCol))
-      allocate(bc_out%lmrcanopy_pa(maxPatchesPerCol))
-      allocate(bc_out%psncanopy_pa(maxPatchesPerCol))
 
       ! Canopy Radiation
       allocate(bc_out%albd_parb(maxPatchesPerCol,cp_numSWb))
@@ -481,9 +478,6 @@ subroutine zero_bcs(this,s)
 
       this%bc_out(s)%rssun_pa(:)     = 0.0_r8
       this%bc_out(s)%rssha_pa(:)     = 0.0_r8
-      this%bc_out(s)%gccanopy_pa(:)  = 0.0_r8
-      this%bc_out(s)%psncanopy_pa(:) = 0.0_r8
-      this%bc_out(s)%lmrcanopy_pa(:) = 0.0_r8
 
       this%bc_out(s)%albd_parb(:,:) = 0.0_r8
       this%bc_out(s)%albi_parb(:,:) = 0.0_r8

components/clm/src/ED/main/FatesParameterDerivedMod.F90

@@ -0,0 +1,116 @@
+module FatesParameterDerivedMod
+
+  ! -------------------------------------------------------------------------------------
+  ! This module contains all procedures types and settings for any quantities that are
+  ! statically derived from static model parameters.  These are unchanging quantities
+  ! and are based off of simple relationships from parameters that the user can
+  ! vary.  This should be called once, and early in the model initialization call
+  ! sequence immediately after FATES parameters are read in.
+  ! -------------------------------------------------------------------------------------
+
+  use FatesConstantsMod, only : r8 => fates_r8
+  use FatesConstantsMod, only : umolC_to_kgC
+  use FatesConstantsMod, only : g_per_kg
+
+  type param_derived_type
+
+     real(r8), allocatable :: vcmax25top(:) ! canopy top: maximum rate of carboxylation
+                                            ! at 25C (umol CO2/m**2/s)
+     real(r8), allocatable :: jmax25top(:)  ! canopy top: maximum electron transport 
+                                            ! rate at 25C (umol electrons/m**2/s)
+     real(r8), allocatable :: tpu25top(:)   ! canopy top: triose phosphate utilization
+                                            ! rate at 25C (umol CO2/m**2/s)
+     real(r8), allocatable :: kp25top(:)    ! canopy top: initial slope of CO2 response
+                                            ! curve (C4 plants) at 25C
+     real(r8), allocatable :: lmr25top(:)   ! canopy top: leaf maintenance respiration
+                                            ! rate at 25C (umol CO2/m**2/s)
+   contains
+
+     procedure :: Init
+     procedure :: InitAllocate
+
+  end type param_derived_type
+
+  type(param_derived_type) :: param_derived
+
+contains
+
+  subroutine InitAllocate(this,maxpft)
+
+    class(param_derived_type), intent(inout) :: this
+    integer, intent(in)                                 :: maxpft
+
+    allocate(this%vcmax25top(maxpft))
+    allocate(this%jmax25top(maxpft))
+    allocate(this%tpu25top(maxpft))
+    allocate(this%kp25top(maxpft))
+    allocate(this%lmr25top(maxpft))
+
+    return
+  end subroutine InitAllocate
+
+  ! =====================================================================================
+
+  subroutine Init(this,maxpft)
+
+    use pftconMod        , only: pftcon
+
+    class(param_derived_type), intent(inout) :: this
+    integer, intent(in)                      :: maxpft
+
+    ! local variables
+    integer  :: ft                 ! pft index
+    real(r8) :: lnc                ! leaf N concentration (gN leaf/m^2)
+
+    associate( &
+         slatop    => pftcon%slatop , & ! specific leaf area at top of canopy, 
+                                        ! projected area basis [m^2/gC]
+         fnitr     => pftcon%fnitr  , & ! foliage nitrogen limitation factor (-)
+         leafcn    => pftcon%leafcn )   ! leaf C:N (gC/gN)
+
+      call this%InitAllocate(maxpft)
+
+      do ft = 1,maxpft
+
+         ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf)
+         lnc  = 1._r8 / (slatop(ft) * leafcn(ft))
+
+         ! at the moment in ED we assume that there is no active N cycle. 
+         ! This should change, of course. FIX(RF,032414) Sep2011. 
+         ! fudge - shortcut using fnitr as a proxy for vcmax... 
+         this%vcmax25top(ft) = fnitr(ft)
+
+         ! Parameters derived from vcmax25top. 
+         ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593
+         ! used jmax25 = 1.97 vcmax25, from Wullschleger (1993) Journal of 
+         ! Experimental Botany 44:907-920.  Here use a factor "1.67", from 
+         ! Medlyn et al (2002) Plant, Cell and Environment 25:1167-1179
+
+         ! RF - copied this from the CLM trunk code, but where did it come from, 
+         ! and how can we make these consistant? 
+         ! jmax25top(ft) =  &
+         ! (2.59_r8 - 0.035_r8*min(max((t10(p)-tfrzc),11._r8),35._r8)) * vcmax25top(ft)
+
+         this%jmax25top(ft) = 1.67_r8   * this%vcmax25top(ft)
+         this%tpu25top(ft)  = 0.167_r8  * this%vcmax25top(ft)
+         this%kp25top(ft)   = 20000._r8 * this%vcmax25top(ft)
+
+         ! Leaf maintenance respiration to match the base rate used in CN
+         ! but with the new temperature functions for C3 and C4 plants.
+         !
+         !
+         ! CN respiration has units:  g C / g N [leaf] / s. This needs to be
+         ! converted from g C / g N [leaf] / s to umol CO2 / m**2 [leaf] / s
+         !
+         ! Then scale this value at the top of the canopy for canopy depth
+
+         this%lmr25top(ft) = 2.525e-6_r8 * (1.5_r8 ** ((25._r8 - 20._r8)/10._r8))
+         this%lmr25top(ft) = this%lmr25top(ft) * lnc / (umolC_to_kgC * g_per_kg)
+
+      end do !ft 
+    end associate
+    return
+  end subroutine Init
+
+
+end module FatesParameterDerivedMod

components/clm/src/biogeophys/PhotosynthesisMod.F90

@@ -162,10 +162,6 @@ module  PhotosynthesisMod
      real(r8), pointer, public  :: rssun_patch       (:)   ! patch sunlit stomatal resistance (s/m)
      real(r8), pointer, public  :: rssha_patch       (:)   ! patch shaded stomatal resistance (s/m)
 
-     ! ED specific variables
-!     real(r8), pointer, public  :: psncanopy_patch   (:)   ! patch sunlit leaf photosynthesis (umol CO2 /m**2/ s) (ED specific)
-!     real(r8), pointer, public  :: lmrcanopy_patch   (:)   ! sunlit leaf maintenance respiration rate (umol CO2/m**2/s) (ED specific)
-
      ! 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 

components/clm/src/utils/clmfates_interfaceMod.F90

@@ -40,7 +40,6 @@ module CLMFatesInterfaceMod
    use TemperatureType   , only : temperature_type
    use EnergyFluxType    , only : energyflux_type
    use SoilStateType     , only : soilstate_type
-   use PhotosynthesisMod , only : photosyns_type
    use clm_varctl        , only : iulog, use_ed
    use clm_varcon        , only : tfrz
    use clm_varcon        , only : spval
@@ -50,6 +49,7 @@ module CLMFatesInterfaceMod
                                   nlevsoi,           &
                                   nlevdecomp,        &
                                   nlevdecomp_full
+   use PhotosynthesisMod , only : photosyns_type
    use atm2lndType       , only : atm2lnd_type
    use SurfaceAlbedoType , only : surfalb_type
    use SolarAbsorbedType , only : solarabs_type
@@ -100,7 +100,7 @@ module CLMFatesInterfaceMod
    use EDBtranMod            , only : btran_ed, &
                                       get_active_suction_layers
    use EDCanopyStructureMod  , only : canopy_summarization, update_hlm_dynamics
-   use EDPhotosynthesisMod   , only : Photosynthesis_ED
+   use FatesPlantRespPhotosynthMod, only : FatesPlantRespPhotosynthDrive
    use EDAccumulateFluxesMod , only : AccumulateFluxes_ED
    use EDPhysiologyMod       , only : flux_into_litter_pools
 
@@ -192,8 +192,9 @@ subroutine init(this, bounds_proc, use_ed)
       ! is not turned on
       ! ---------------------------------------------------------------------------------
 
-     use FatesInterfaceMod, only : FatesInterfaceInit
-
+      use FatesInterfaceMod, only : FatesInterfaceInit
+      use EDTypesMod       , only : numpft_ed
+      use FatesParameterDerivedMod, only : param_derived
 
       implicit none
 
@@ -215,10 +216,9 @@ subroutine init(this, bounds_proc, use_ed)
          ! This involves to stages
          ! 1) allocate the vectors
          ! 2) add the history variables defined in clm_inst to the history machinery
-!         call this%fates2hlm%Init(bounds_proc)
-
          call EDecophysconInit( EDpftvarcon_inst, numpft )
-
+         call param_derived%Init(numpft_ed)
+
       end if
 
       if(DEBUG)then
@@ -1341,11 +1341,11 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, &
 
       ! Call photosynthesis
 
-      call  Photosynthesis_ED (this%fates(nc)%nsites, &
-           this%fates(nc)%sites,  &
-           this%fates(nc)%bc_in,  &
-           this%fates(nc)%bc_out, &
-           dtime)
+      call FatesPlantRespPhotosynthDrive (this%fates(nc)%nsites, &
+                                this%fates(nc)%sites,  &
+                                this%fates(nc)%bc_in,  &
+                                this%fates(nc)%bc_out, &
+                                dtime)
 
       ! Perform a double check to see if all patches on naturally vegetated columns
       ! were activated for photosynthesis