Merge pull request #1128 from mpaiao/mpaiao-pr-allom-modes

New allometric modes

biogeochem/EDCanopyStructureMod.F90

@@ -22,6 +22,7 @@ module EDCanopyStructureMod
   use FatesAllometryMod     , only : tree_sai
   use EDTypesMod            , only : ed_site_type
   use FatesAllometryMod     , only : VegAreaLayer
+  use FatesAllometryMod     , only : CrownDepth
   use FatesPatchMod,          only : fates_patch_type
   use FatesCohortMod,         only : fates_cohort_type
   use EDParamsMod            , only : nclmax
@@ -1532,7 +1533,7 @@ subroutine leaf_area_profile( currentSite )
     real(r8) :: elai_layer,tlai_layer    ! leaf area per canopy area
     real(r8) :: esai_layer,tsai_layer    ! stem area per canopy area
     real(r8) :: vai_top,vai_bot          ! integrated top down veg area index at boundary of layer
-
+    real(r8) :: crown_depth              ! Current cohort's crown depth
     real(r8) :: layer_bottom_height,layer_top_height,lai,sai  ! Can be removed later
     !----------------------------------------------------------------------
 
@@ -1611,7 +1612,14 @@ subroutine leaf_area_profile( currentSite )
                    write(fates_log(), *) ' currentPatch%nrad(cl,ft): ', currentPatch%nrad(cl,ft)
                    call endrun(msg=errMsg(sourcefile, __LINE__))
                 end if
+
+                !---~---
+                !   Find current crown depth using the allometric function.
+                !---~---
+                call CrownDepth(currentCohort%height,currentCohort%pft,crown_depth)
+                !---~---
 
+
                 ! --------------------------------------------------------------------------
                 ! Whole layers.  Make a weighted average of the leaf area in each layer
                 ! before dividing it by the total area. Fill up layer for whole layers.
@@ -1625,12 +1633,10 @@ subroutine leaf_area_profile( currentSite )
                    ! is obscured by snow.
 
                    layer_top_height = currentCohort%height - &
-                        ( real(iv-1,r8)/currentCohort%NV * currentCohort%height *  &
-                        prt_params%crown_depth_frac(currentCohort%pft) )
-
+                        ( real(iv-1,r8)/currentCohort%NV * crown_depth )
+
                    layer_bottom_height = currentCohort%height - &
-                        ( real(iv,r8)/currentCohort%NV * currentCohort%height * &
-                        prt_params%crown_depth_frac(currentCohort%pft) )
+                        ( real(iv,r8)/currentCohort%NV * crown_depth )
 
                    fraction_exposed = 1.0_r8
                    if(currentSite%snow_depth  > layer_top_height)then

main/EDParamsMod.F90

@@ -77,6 +77,7 @@ module EDParamsMod
    real(r8),protected, public :: ED_val_patch_fusion_tol              ! minimum fraction in difference in profiles between patches
    real(r8),protected, public :: ED_val_canopy_closure_thresh         ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry
    integer,protected, public  :: stomatal_model                       ! switch for choosing between stomatal conductance models, 1 for Ball-Berry, 2 for Medlyn
+   integer,protected, public  :: dayl_switch                          ! switch for turning on or off day length factor scaling for photosynthetic parameters
    integer,protected, public  :: regeneration_model                   ! Switch for choosing between regeneration models:
                                                                       ! (1) for Fates default
                                                                       ! (2) for the Tree Recruitment Scheme (Hanbury-Brown et al., 2022)
@@ -166,6 +167,7 @@ module EDParamsMod
    character(len=param_string_length),parameter,public :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
    character(len=param_string_length),parameter,public :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"      
    character(len=param_string_length),parameter,public :: ED_name_stomatal_model= "fates_leaf_stomatal_model"
+   character(len=param_string_length),parameter,public :: ED_name_dayl_switch= "fates_daylength_factor_switch"
    character(len=param_string_length),parameter,public :: ED_name_regeneration_model= "fates_regeneration_model"
 
    character(len=param_string_length),parameter,public :: name_theta_cj_c3 = "fates_leaf_theta_cj_c3"
@@ -338,6 +340,7 @@ subroutine FatesParamsInit()
     ED_val_patch_fusion_tol               = nan
     ED_val_canopy_closure_thresh          = nan
     stomatal_model                        = -9
+    dayl_switch                           = -9
     regeneration_model                    = -9
     stomatal_assim_model                  = -9
     max_cohort_per_patch                  = -9
@@ -430,9 +433,8 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
-    ! Temporary until we add parameter to file
-    !call fates_params%RegisterParameter(name=ED_name_mort_cstarvation_model, dimension_shape=dimension_shape_scalar, &
-    !     dimension_names=dim_names_scalar)
+    call fates_params%RegisterParameter(name=ED_name_mort_cstarvation_model, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
 
     call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
@@ -493,7 +495,10 @@ subroutine FatesRegisterParams(fates_params)
 
     call fates_params%RegisterParameter(name=ED_name_stomatal_model, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
-
+
+    call fates_params%RegisterParameter(name=ED_name_dayl_switch, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=ED_name_regeneration_model, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -648,12 +653,10 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetrieveParameter(name=ED_name_mort_disturb_frac, &
           data=fates_mortality_disturbance_fraction)
 
-    ! Temporary until we add parameter to file 
-    !call fates_params%RetrieveParameter(name=ED_name_mort_cstarvation_model, &
-    !     data=tmpreal)
-    !    mort_cstarvation_model = nint(tmpreal)
-    mort_cstarvation_model = cstarvation_model_lin
-
+    call fates_params%RetrieveParameter(name=ED_name_mort_cstarvation_model, &
+         data=tmpreal)
+    mort_cstarvation_model = nint(tmpreal)
+
     call fates_params%RetrieveParameter(name=ED_name_comp_excln, &
          data=ED_val_comp_excln)
 
@@ -715,6 +718,10 @@ subroutine FatesReceiveParams(fates_params)
          data=tmpreal)
     stomatal_model = nint(tmpreal)
 
+    call fates_params%RetrieveParameter(name=ED_name_dayl_switch, &
+         data=tmpreal)
+    dayl_switch = nint(tmpreal)
+
     call fates_params%RetrieveParameter(name=ED_name_regeneration_model, &
          data=tmpreal)
     regeneration_model = nint(tmpreal)
@@ -884,6 +891,7 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh
         write(fates_log(),fmt0) 'regeneration_model = ',regeneration_model      
         write(fates_log(),fmt0) 'stomatal_model = ',stomatal_model
+        write(fates_log(),fmt0) 'dayl_switch = ',dayl_switch
         write(fates_log(),fmt0) 'stomatal_assim_model = ',stomatal_assim_model            
         write(fates_log(),fmt0) 'hydro_kmax_rsurf1 = ',hydr_kmax_rsurf1
         write(fates_log(),fmt0) 'hydro_kmax_rsurf2 = ',hydr_kmax_rsurf2  

main/EDPftvarcon.F90

@@ -92,6 +92,13 @@ module EDPftvarcon
 
      real(r8), allocatable :: maintresp_leaf_ryan1991_baserate(:)  ! leaf maintenance respiration per Ryan et al 1991
 
+     real(r8), allocatable :: leafn_vert_scaler_coeff1(:)  ! Coefficient one for decrease of leaf N through the canopy
+     real(r8), allocatable :: leafn_vert_scaler_coeff2(:)  ! Coefficient two for decrease of leaf N through the canopy 
+
+     real(r8), allocatable :: maintresp_leaf_vert_scaler_coeff1(:) ! leaf maintenance respiration decrease through the canopy param 1
+                                                                   ! only with Atkin et al. 2017 respiration model
+     real(r8), allocatable :: maintresp_leaf_vert_scaler_coeff2(:) ! leaf maintenance respiration decrease through the canopy param 2
+                                                                   ! only with Atkin et al. 2017 respiraiton model 
      real(r8), allocatable :: bmort(:)
      real(r8), allocatable :: mort_ip_size_senescence(:) ! inflection point of dbh dependent senescence
      real(r8), allocatable :: mort_r_size_senescence(:)  ! rate of change in mortality with dbh
@@ -471,6 +478,22 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
         dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    name = 'fates_leafn_vert_scaler_coeff1'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+        dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_leafn_vert_scaler_coeff2'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_maintresp_leaf_vert_scaler_coeff1'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+        dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
+    name = 'fates_maintresp_leaf_vert_scaler_coeff2'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_prescribed_npp_canopy'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -912,6 +935,22 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%maintresp_leaf_ryan1991_baserate)
 
+    name = 'fates_leafn_vert_scaler_coeff1'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%leafn_vert_scaler_coeff1)
+
+    name = 'fates_leafn_vert_scaler_coeff2'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%leafn_vert_scaler_coeff2)
+
+    name = 'fates_maintresp_leaf_vert_scaler_coeff1'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%maintresp_leaf_vert_scaler_coeff1)
+
+    name = 'fates_maintresp_leaf_vert_scaler_coeff2'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%maintresp_leaf_vert_scaler_coeff2)
+
     name = 'fates_prescribed_npp_canopy'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%prescribed_npp_canopy)
@@ -1737,6 +1776,10 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'hydro_vg_alpha_node  = ',EDPftvarcon_inst%hydr_vg_alpha_node
         write(fates_log(),fmt0) 'hydro_vg_m_node  = ',EDPftvarcon_inst%hydr_vg_m_node
         write(fates_log(),fmt0) 'hydro_vg_n_node  = ',EDPftvarcon_inst%hydr_vg_n_node
+        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff1 = ',EDPftvarcon_inst%leafn_vert_scaler_coeff1
+        write(fates_log(),fmt0) 'leafn_vert_scaler_coeff2 = ',EDPftvarcon_inst%leafn_vert_scaler_coeff2
+        write(fates_log(),fmt0) 'maintresp_leaf_vert_scaler_coeff1 = ',EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff1
+        write(fates_log(),fmt0) 'maintresp_leaf_vert_scaler_coeff2 = ',EDPftvarcon_inst%maintresp_leaf_vert_scaler_coeff2
         write(fates_log(),*) '-------------------------------------------------'
 
      end if
@@ -1845,7 +1888,6 @@ subroutine FatesCheckParams(is_master)
                     call endrun(msg=errMsg(sourcefile, __LINE__))
                  end if
               end if
-
            end if
         end if
 
@@ -1871,7 +1913,25 @@ subroutine FatesCheckParams(is_master)
               call endrun(msg=errMsg(sourcefile, __LINE__))
            end if
         end if
-
+
+        ! We are using a simple phosphatase model right now. There is
+        ! no critical value (lambda) , and there is no preferential uptake (alpha).
+        ! Make sure these parameters are both set to 0.
+
+        if ((hlm_phosphorus_spec>0) .and. (trim(hlm_nu_com).eq.'ECA')) then
+           if (any(abs(EDPftvarcon_inst%eca_lambda_ptase(:)) > nearzero ) ) then
+              write(fates_log(),*) 'Critical Values for phosphatase in ECA are not'
+              write(fates_log(),*) 'enabled right now. Please set fates_eca_lambda_ptase = 0'
+              write(fates_log(),*) ' Aborting'
+              call endrun(msg=errMsg(sourcefile, __LINE__))
+           end if
+           if (any(abs(EDPftvarcon_inst%eca_alpha_ptase(:)) > nearzero ) ) then
+              write(fates_log(),*) 'There is no preferential plant uptake of P from phosphatase'
+              write(fates_log(),*) 'enabled right now. Please set fates_eca_alpha_ptase = 0'
+              write(fates_log(),*) ' Aborting'
+              call endrun(msg=errMsg(sourcefile, __LINE__))
+           end if
+        end if
 
      case (prt_carbon_allom_hyp)
         ! No additional checks needed for now.
@@ -1915,6 +1975,9 @@ subroutine FatesCheckParams(is_master)
         end if
      end if
 
+
+
+
      do ipft = 1,npft
 
         ! xl must be between -0.4 and 0.6 according to Bonan (2019) doi:10.1017/9781107339217 pg. 238