Merge pull request #1 from rgknox/fates_master_medlyn

merge medlyn changes up to master
ckoven · Jun 8, 2020 · 7f84a4b · 7f84a4b
2 parents b636984 + 5e9e1ee
commit 7f84a4b
Show file tree

Hide file tree

Showing 15 changed files with 229 additions and 219 deletions.
diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90
@@ -64,7 +64,6 @@ module EDCohortDynamicsMod
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : ForceDBH
   use FatesAllometryMod  , only : tree_lai, tree_sai
-  use FatesAllometryMod  , only : i_biomass_rootprof_context 
   use FatesAllometryMod    , only : set_root_fraction
   use PRTGenericMod,          only : prt_carbon_allom_hyp   
   use PRTGenericMod,          only : prt_cnp_flex_allom_hyp
@@ -847,8 +846,7 @@ subroutine SendCohortToLitter(csite,cpatch,ccohort,nplant)
 
     plant_dens = nplant/cpatch%area
 
-    call set_root_fraction(csite%rootfrac_scr, pft, csite%zi_soil, &
-         icontext = i_biomass_rootprof_context)
+    call set_root_fraction(csite%rootfrac_scr, pft, csite%zi_soil)
 
     do el=1,num_elements
 

diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90
@@ -55,7 +55,6 @@ module EDLoggingMortalityMod
    use PRTGenericMod     , only : sapw_organ, struct_organ, leaf_organ
    use PRTGenericMod     , only : fnrt_organ, store_organ, repro_organ
    use FatesAllometryMod , only : set_root_fraction
-   use FatesAllometryMod , only : i_biomass_rootprof_context
 
    implicit none
    private
@@ -409,9 +408,7 @@ subroutine logging_litter_fluxes(currentSite, currentPatch, newPatch, patch_site
             ! derived from the current patch, so we need to multiply by patch_areadis/np%area
             ! ----------------------------------------------------------------------------------------
 
-            call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, &
-                  icontext = i_biomass_rootprof_context)
-
+            call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil)
 
             ag_wood = (direct_dead+indirect_dead) * (struct_m + sapw_m ) * &
                   EDPftvarcon_inst%allom_agb_frac(currentCohort%pft)

diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -53,7 +53,6 @@ module EDPatchDynamicsMod
   use EDParamsMod          , only : fates_mortality_disturbance_fraction
   use FatesAllometryMod    , only : carea_allom
   use FatesAllometryMod    , only : set_root_fraction
-  use FatesAllometryMod    , only : i_biomass_rootprof_context
   use FatesConstantsMod    , only : g_per_kg
   use FatesConstantsMod    , only : ha_per_m2
   use FatesConstantsMod    , only : days_per_sec
@@ -1531,8 +1530,7 @@ subroutine fire_litter_fluxes(currentSite, currentPatch, newPatch, patch_site_ar
 
              site_mass%burn_flux_to_atm = site_mass%burn_flux_to_atm + burned_mass
 
-             call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, &
-                   icontext = i_biomass_rootprof_context)
+             call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil)
 
              ! Contribution of dead trees to root litter (no root burn flux to atm)
              do dcmpy=1,ndcmpy
@@ -1741,8 +1739,8 @@ subroutine mortality_litter_fluxes(currentSite, currentPatch, newPatch, patch_si
           ag_wood = num_dead * (struct_m + sapw_m) * EDPftvarcon_inst%allom_agb_frac(pft)
           bg_wood = num_dead * (struct_m + sapw_m) * (1.0_r8-EDPftvarcon_inst%allom_agb_frac(pft))
 
-          call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, &
-                icontext = i_biomass_rootprof_context)
+          call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil)
+
 
           do c=1,ncwd
 

diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90
@@ -81,7 +81,6 @@ module EDPhysiologyMod
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : CheckIntegratedAllometries
   use FatesAllometryMod, only : set_root_fraction
-  use FatesAllometryMod, only : i_biomass_rootprof_context 
 
   use PRTGenericMod, only : prt_carbon_allom_hyp
   use PRTGenericMod, only : prt_cnp_flex_allom_hyp
@@ -1797,8 +1796,7 @@ subroutine CWDInput( currentSite, currentPatch, litt)
     do while(associated(currentCohort))
       pft = currentCohort%pft        
 
-      call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil, &
-            icontext = i_biomass_rootprof_context)
+      call set_root_fraction(currentSite%rootfrac_scr, pft, currentSite%zi_soil)
 
       leaf_m_turnover   = currentCohort%prt%GetTurnover(leaf_organ,element_id)
       store_m_turnover  = currentCohort%prt%GetTurnover(store_organ,element_id)

diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90
@@ -128,22 +128,6 @@ module FatesAllometryMod
   character(len=*), parameter :: sourcefile = __FILE__
 
 
-  ! The code will call the wrapper routine "set_root_fraction"
-  ! in at least two different context.  In one context it will query
-  ! set_root_fraction to describe the depth profile of hydraulicly
-  ! active roots.  In the other context, it will ask the wrapper
-  ! to define the profile of roots as litter.  We allow these
-  ! two contexts to differ.  While not fully implemented, the use
-  ! will have control parameters to choose from different relationships
-  ! in these two contexts.  The calling function, therefore
-  ! has to tell the wrapper function which context (water or biomass)
-  ! is being querried.  So that we don't have to do messy string
-  ! parsing, we have two pre-defined flags.
-
-  integer, parameter, public :: i_hydro_rootprof_context  = 1
-  integer, parameter, public :: i_biomass_rootprof_context = 2
-
-
   ! If testing b4b with older versions, do not remove sapwood
   ! Our old methods with saldarriaga did not remove sapwood from the
   ! bdead pool.  But newer allometries are providing total agb
@@ -1977,7 +1961,8 @@ end subroutine carea_2pwr
 
   ! =========================================================================
 
-  subroutine set_root_fraction(root_fraction, ft, zi, icontext )
+  subroutine set_root_fraction(root_fraction, ft, zi)
+
     !
     ! !DESCRIPTION:
     !  Calculates the fractions of the root biomass in each layer for each pft. 
@@ -1992,7 +1977,6 @@ subroutine set_root_fraction(root_fraction, ft, zi, icontext )
     real(r8),intent(inout) :: root_fraction(:) ! Normalized profile
     integer, intent(in)    :: ft               ! functional typpe
     real(r8),intent(in)    :: zi(0:)            ! Center of depth [m]
-    integer,intent(in)     :: icontext
 
     ! Parameters
     !
@@ -2009,8 +1993,8 @@ subroutine set_root_fraction(root_fraction, ft, zi, icontext )
     ! exponential.
     ! All methods return a normalized profile.
 
-    integer, parameter :: exponential_1p_profile_type = 1
-    integer, parameter :: jackson_beta_profile_type   = 2
+    integer, parameter :: jackson_beta_profile_type   = 1
+    integer, parameter :: exponential_1p_profile_type = 2
     integer, parameter :: exponential_2p_profile_type = 3
 
     integer :: root_profile_type
@@ -2027,40 +2011,21 @@ subroutine set_root_fraction(root_fraction, ft, zi, icontext )
        call endrun(msg=errMsg(sourcefile, __LINE__))
     end if
 
-    if(icontext == i_hydro_rootprof_context) then
-
-       root_profile_type = exponential_2p_profile_type
-
-    else if(icontext == i_biomass_rootprof_context) then
-
-       root_profile_type = jackson_beta_profile_type
-
-    else
-       write(fates_log(),*) 'An undefined context for calculating root profiles was provided'
-       write(fates_log(),*) 'There are only two contexts, hydraulic and biomass, pick one.'
-       write(fates_log(),*) 'Aborting'
-       call endrun(msg=errMsg(sourcefile, __LINE__))
-    end if
-
 
-    select case(root_profile_type)
+    select case(nint(EDPftvarcon_inst%fnrt_prof_mode(ft)))
     case ( exponential_1p_profile_type ) 
-       call exponential_1p_root_profile(root_fraction, ft, zi) 
+       call exponential_1p_root_profile(root_fraction, ft, zi, EDPftvarcon_inst%fnrt_prof_a(ft)) 
     case ( jackson_beta_profile_type )
-       call jackson_beta_root_profile(root_fraction, ft, zi)
+       call jackson_beta_root_profile(root_fraction, ft, zi, EDPftvarcon_inst%fnrt_prof_a(ft))
     case ( exponential_2p_profile_type ) 
-       call exponential_2p_root_profile(root_fraction, ft, zi)
+       call exponential_2p_root_profile(root_fraction, ft, zi, & 
+             EDPftvarcon_inst%fnrt_prof_a(ft),EDPftvarcon_inst%fnrt_prof_b(ft))
     case default
        write(fates_log(),*) 'An undefined root profile type was specified'
        write(fates_log(),*) 'Aborting'
        call endrun(msg=errMsg(sourcefile, __LINE__))
     end select
 
-!    if( abs(sum(root_fraction)-1.0_r8) > 1.e-9_r8 ) then
-!        write(fates_log(),*) 'Root fractions should add up to 1'
-!        write(fates_log(),*) root_fraction
-!        call endrun(msg=errMsg(sourcefile, __LINE__))
-!    end if
 
     correction = 1._r8 - sum(root_fraction)
     corr_id = maxloc(root_fraction)
@@ -2073,27 +2038,49 @@ end subroutine set_root_fraction
 
   ! =====================================================================================
 
-  subroutine exponential_2p_root_profile(root_fraction, ft, zi )
+  subroutine exponential_2p_root_profile(root_fraction, ft, zi, a, b)
     !
     ! !ARGUMENTS
     real(r8),intent(out) :: root_fraction(:)
     integer,intent(in)   :: ft
     real(r8),intent(in)  :: zi(0:)
+    real(r8),intent(in)  :: a     ! Exponential shape parameter a
+    real(r8),intent(in)  :: b     ! Exponential shape parameter b
 
     ! Locals
     integer  :: nlevsoil    ! Number of soil layers
     integer  :: lev         ! soil layer index
     real(r8) :: sum_rootfr  ! sum of root fraction for normalization
-
+
+
+    ! Original default parameters: 
+    !
+    ! broadleaf_evergreen_tropical_tree
+    ! needleleaf_evergreen_extratrop_tree
+    ! needleleaf_colddecid_extratrop_tree
+    ! broadleaf_evergreen_extratrop_tree
+    ! broadleaf_hydrodecid_tropical_tree
+    ! broadleaf_colddecid_extratrop_tree
+    ! broadleaf_evergreen_extratrop_shrub
+    ! broadleaf_hydrodecid_extratrop_shrub
+    ! broadleaf_colddecid_extratrop_shrub
+    ! arctic_c3_grass
+    ! cool_c3_grass
+    ! c4_grass
+    !
+    ! a = 7, 7, 7, 7, 6, 6, 7, 7, 7, 11, 11, 11 ;
+    ! b = 1, 2, 2, 1, 2, 2, 1.5, 1.5, 1.5, 2, 2, 2 ;
+
+
     nlevsoil = ubound(zi,1)
 
     sum_rootfr = 0.0_r8
     do lev = 1, nlevsoil
        root_fraction(lev) = .5_r8*( &
-             exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev-1))  &
-             + exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev-1))  &
-             - exp(-EDPftvarcon_inst%roota_par(ft) * zi(lev))    &
-             - exp(-EDPftvarcon_inst%rootb_par(ft) * zi(lev)))
+             exp(-a * zi(lev-1))  &
+             + exp(-b * zi(lev-1))  &
+             - exp(-a * zi(lev))    &
+             - exp(-b * zi(lev)))
 
        sum_rootfr = sum_rootfr + root_fraction(lev)
     end do
@@ -2106,30 +2093,32 @@ end subroutine exponential_2p_root_profile
 
   ! =====================================================================================
 
-  subroutine exponential_1p_root_profile(root_fraction, ft, zi)
+  subroutine exponential_1p_root_profile(root_fraction, ft, zi, a)
 
     !
     ! !ARGUMENTS
     real(r8),intent(out) :: root_fraction(:)
     integer,intent(in)   :: ft
     real(r8),intent(in)  :: zi(0:)
-
+    real(r8),intent(in)  :: a     ! Exponential shape parameter a
+
     !
     ! LOCAL VARIABLES:
     integer :: lev         ! soil depth layer index
     integer :: nlevsoil    ! number of soil layers
     real(r8) :: depth      ! Depth to middle of layer [m]
     real(r8) :: sum_rootfr ! sum of rooting profile for normalization
 
-    real(r8), parameter :: rootprof_exp  = 3.  ! how steep profile is
+    ! Typical default parameter is a = 3.  
+    ! how steep profile is
     ! for root C inputs (1/ e-folding depth) (1/m)
 
     nlevsoil = ubound(zi,1)
 
     ! define rooting profile from exponential parameters
     sum_rootfr = 0.0_r8
     do lev = 1,  nlevsoil
-       root_fraction(lev) = exp(-rootprof_exp * 0.5*(zi(lev)+zi(lev-1)) )
+       root_fraction(lev) = exp(-a * 0.5*(zi(lev)+zi(lev-1)) )
        sum_rootfr = sum_rootfr + root_fraction(lev)
     end do
 
@@ -2142,33 +2131,32 @@ end subroutine exponential_1p_root_profile
 
   ! =====================================================================================
 
-  subroutine jackson_beta_root_profile(root_fraction, ft, zi)
+  subroutine jackson_beta_root_profile(root_fraction, ft, zi, a)
+
+    ! -----------------------------------------------------------------------------------
+    ! use beta distribution parameter from Jackson et al., 1996
+    ! -----------------------------------------------------------------------------------
 
-
     ! !ARGUMENTS
     real(r8),intent(out) :: root_fraction(:) ! fraction of root mass in each soil layer
     integer,intent(in)   :: ft               ! functional type
     real(r8),intent(in)  :: zi(0:)           ! depth of layer interfaces 0-nlevsoil
-
+    real(r8),intent(in)  :: a                ! Exponential shape parameter a
+
     !
     ! LOCAL VARIABLES:
     integer :: lev         ! soil depth layer index
     integer :: nlevsoil    ! number of soil layers
     real(r8) :: sum_rootfr ! sum of rooting profile, for normalization 
 
-    ! Note cdk 2016/08 we actually want to use the carbon index here rather than the water index.  
-    ! Doing so will be answer changing though so perhaps easiest to do this in steps.
-    integer, parameter :: rooting_profile_varindex_water = 1
+    ! Original defaults in fates, a = 0.976 (all Pfts)
 
     nlevsoil = ubound(zi,1)
-    ! use beta distribution parameter from Jackson et al., 1996
+
     sum_rootfr = 0.0_r8
     do lev = 1, nlevsoil
        root_fraction(lev) = &
-             ( EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & 
-             ( zi(lev-1)*100._r8) - &
-             EDPftvarcon_inst%rootprof_beta(ft, rooting_profile_varindex_water) ** & 
-             ( zi(lev)*100._r8) )
+             ( a ** ( zi(lev-1)*100._r8) - a ** ( zi(lev)*100._r8) )
        sum_rootfr = sum_rootfr + root_fraction(lev)
     end do
 

diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90
@@ -19,7 +19,6 @@ module EDBtranMod
    use FatesInterfaceTypesMod , only : hlm_use_planthydro
    use FatesGlobals      , only : fates_log
    use FatesAllometryMod , only : set_root_fraction
-   use FatesAllometryMod , only : i_hydro_rootprof_context
 
    !
    implicit none
@@ -140,8 +139,7 @@ subroutine btran_ed( nsites, sites, bc_in, bc_out)
 
               do ft = 1,numpft
 
-                  call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, &
-                      icontext = i_hydro_rootprof_context)
+                  call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil ) 
 
                  cpatch%btran_ft(ft) = 0.0_r8
                  do j = 1,bc_in(s)%nlevsoil

diff --git a/biogeophys/FatesBstressMod.F90 b/biogeophys/FatesBstressMod.F90
@@ -19,7 +19,6 @@ module FatesBstressMod
    use FatesGlobals      , only : fates_log
    use EDBtranMod        , only : check_layer_water
    use FatesAllometryMod , only : set_root_fraction
-   use FatesAllometryMod , only : i_hydro_rootprof_context
 
    implicit none
    private
@@ -69,8 +68,7 @@ subroutine btran_sal_stress_fates( nsites, sites, bc_in)
               do ft = 1,numpft
                  cpatch%bstress_sal_ft(ft) = 0.0_r8
 
-                 call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil, &
-                       icontext = i_hydro_rootprof_context)
+                 call set_root_fraction(sites(s)%rootfrac_scr, ft, sites(s)%zi_soil )
 
                  do j = 1,bc_in(s)%nlevsoil