NGEET · rgknox · Jun 21, 2022 · Apr 27, 2022 · Apr 27, 2022 · Apr 27, 2022
diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90
diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90
@@ -2008,7 +2008,7 @@ subroutine CrownDepth(height,ft,crown_depth)
     ! Original FATES crown depth heigh used for hydraulics
     ! crown_depth               = min(height,0.1_r8)
 
-    crown_depth = prt_params%crown(ft) * height
+    crown_depth = prt_params%crown_depth_frac(ft) * height
 
 
     return

diff --git a/biogeophys/FatesPlantHydraulicsMod.F90 b/biogeophys/FatesPlantHydraulicsMod.F90
@@ -48,7 +48,7 @@ module FatesPlantHydraulicsMod
   use EDParamsMod       , only : hydr_psi0
   use EDParamsMod       , only : hydr_psicap
   use EDParamsMod       , only : hydr_htftype_node
-  use EDParamsMod       , only : hydr_solver_type
+  use EDParamsMod       , only : hydr_solver
 
   use EDTypesMod        , only : ed_site_type
   use EDTypesMod        , only : ed_patch_type
@@ -1468,7 +1468,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_none)
 
         csite_hydr%nlevrhiz = bc_in(s)%nlevsoil
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver_type,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
 
         do j=1,csite_hydr%nlevrhiz
            csite_hydr%map_r2s(j,1) = j
@@ -1480,7 +1480,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_combine12)
 
         csite_hydr%nlevrhiz     = max(1,bc_in(s)%nlevsoil-1)
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver_type,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
 
         csite_hydr%map_r2s(1,1) = 1
         j_bc                 = min(2,bc_in(s)%nlevsoil) ! this protects 1 soil layer
@@ -1498,7 +1498,7 @@ subroutine InitHydrSites(sites,bc_in)
      case(rhizlayer_aggmeth_balN)
 
         csite_hydr%nlevrhiz = min(aggN,bc_in(s)%nlevsoil)
-        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver_type,bc_in(s)%nlevsoil)
+        call sites(s)%si_hydr%InitHydrSite(numpft,nlevsclass,hydr_solver,bc_in(s)%nlevsoil)
 
         ntoagg = int(ceiling(real(bc_in(s)%nlevsoil)/real(csite_hydr%nlevrhiz)-nearzero))
 
@@ -2627,21 +2627,21 @@ subroutine hydraulics_bc ( nsites, sites, bc_in, bc_out, dtime)
               ! from leaf to the current soil layer.  This does NOT
               ! update cohort%th_*
 
-              if(hydr_solver_type == hydr_solver_2DNewton) then
+              if(hydr_solver == hydr_solver_2DNewton) then
 
                  call MatSolve2D(csite_hydr,ccohort,ccohort_hydr, &
                       dtime,qflx_tran_veg_indiv, &
                       sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, &
                       dth_layershell_col)
 
-              elseif(hydr_solver_type == hydr_solver_2DPicard) then
+              elseif(hydr_solver == hydr_solver_2DPicard) then
 
                  call PicardSolve2D(csite_hydr,ccohort,ccohort_hydr, &
                       dtime,qflx_tran_veg_indiv, &
                       sapflow,rootuptake(1:nlevrhiz),wb_err_plant,dwat_plant, &
                       dth_layershell_col,csite_hydr%num_nodes)
 
-              elseif(hydr_solver_type == hydr_solver_1DTaylor ) then
+              elseif(hydr_solver == hydr_solver_1DTaylor ) then
 
                  ! ---------------------------------------------------------------------------------
                  ! Approach: do nlevsoi_hyd sequential solutions to Richards' equation,
@@ -4942,7 +4942,7 @@ subroutine MatSolve2D(csite_hydr,cohort,cohort_hydr, &
 
 
    ! This NaN's the scratch arrays
-   call csite_hydr%FlushSiteScratch(hydr_solver_type)
+   call csite_hydr%FlushSiteScratch(hydr_solver)
 
    ! This is the maximum number of iterations needed for this cohort
    ! (each soil layer has a different number, this saves the max)
@@ -5725,7 +5725,7 @@ subroutine PicardSolve2D(csite_hydr,cohort,cohort_hydr, &
 
 
     ! This NaN's the scratch arrays
-    call csite_hydr%FlushSiteScratch(hydr_solver_type)
+    call csite_hydr%FlushSiteScratch(hydr_solver)
 
     ! This is the maximum number of iterations needed for this cohort
     ! (each soil layer has a different number, this saves the max)

diff --git a/biogeophys/FatesPlantRespPhotosynthMod.F90 b/biogeophys/FatesPlantRespPhotosynthMod.F90
@@ -50,7 +50,7 @@ module FATESPlantRespPhotosynthMod
   use PRTGenericMod,     only : store_organ
   use PRTGenericMod,     only : repro_organ
   use PRTGenericMod,     only : struct_organ
-  use EDParamsMod,       only : ED_val_base_mr_20, stomatal_model
+  use EDParamsMod,       only : ED_val_base_mr_20, stomatal_model, stomatal_assim_model
   use PRTParametersMod,  only : prt_params
   use EDPftvarcon         , only : EDPftvarcon_inst
 
@@ -93,12 +93,9 @@ module FATESPlantRespPhotosynthMod
 
   ! Alternatively, Gross Assimilation can be used to estimate
   ! leaf co2 partial pressure and therefore conductance. The default
-  !is to use anet
-  logical, parameter :: use_agross = .false.
-
-
-
-
+  ! is to use anet
+  integer, parameter :: net_assim_model = 1
+  integer, parameter :: gross_assim_model = 2
 
 
 contains
@@ -1152,7 +1149,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
               ! using anet in calculating gs this is version B  
               anet = agross  - lmr
 
-              if (use_agross) then
+              if ( stomatal_assim_model == gross_assim_model ) then
                  if ( stomatal_model == medlyn_model ) then
                     write (fates_log(),*) 'Gross Assimilation conductance is incompatible with the Medlyn model'
                     call endrun(msg=errMsg(sourcefile, __LINE__))

diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90
@@ -93,16 +93,6 @@ module EDParamsMod
    ! 1  = Christofferson et al. 2016 (TFS),   2 = Van Genuchten 1980
    integer, protected,allocatable,public :: hydr_htftype_node(:)
 
-   ! Switch that defines which hydraulic solver to use
-   ! 1 = Taylor solution that solves plant fluxes with 1 layer
-   !     sequentially placing solution on top of previous layer solves
-   ! 2 = Newton-Raphson solution that solves all fluxes in a plant and
-   !     the soil simultaneously, 2D: soil x (root + shell)
-   ! 3 = Picard solution that solves all fluxes in a plant and
-   !     the soil simultaneously, 2D: soil x (root + shell)
-
-   integer, parameter, public :: hydr_solver_type = 1    ! 1 = hydr_solver_1DTaylor
-
    character(len=param_string_length),parameter,public :: ED_name_photo_temp_acclim_timescale = "fates_photo_temp_acclim_timescale"
    character(len=param_string_length),parameter,public :: name_photo_tempsens_model = "fates_photo_tempsens_model"
    character(len=param_string_length),parameter,public :: name_maintresp_model = "fates_maintresp_model"
@@ -160,10 +150,27 @@ module EDParamsMod
    real(r8),protected,public :: hydr_psicap        !  sapwood water potential at which capillary reserves exhausted (MPa)
    character(len=param_string_length),parameter,public :: hydr_name_psicap = "fates_hydr_psicap"
 
+
+   ! Switch that defines which hydraulic solver to use
+   ! 1 = Taylor solution that solves plant fluxes with 1 layer
+   !     sequentially placing solution on top of previous layer solves
+   ! 2 = Picard solution that solves all fluxes in a plant and
+   !     the soil simultaneously, 2D: soil x (root + shell)
+   ! 3 = Newton-Raphson (Deprecated) solution that solves all fluxes in a plant and
+   !     the soil simultaneously, 2D: soil x (root + shell)
+
+   integer,protected,public :: hydr_solver        !  switch designating hydraulics numerical solver
+   character(len=param_string_length),parameter,public :: hydr_name_solver = "fates_hydr_solver"
+
    !Soil BGC parameters, mostly used for testing FATES when not coupled to the dynamics bgc hlm
    ! ----------------------------------------------------------------------------------------------
    real(r8),protected,public :: bgc_soil_salinity ! site-level soil salinity for FATES when not coupled to dynamic soil BGC of salinity
    character(len=param_string_length),parameter,public :: bgc_name_soil_salinity= "fates_soil_salinity"      
+
+   ! Switch designating whether to use net or gross assimilation in the stomata model
+   integer, protected, public :: stomatal_assim_model
+   character(len=param_string_length), parameter, public :: stomatal_assim_name = "fates_stomatal_assim"
+
 
    ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE)
    ! ----------------------------------------------------------------------------------------------
@@ -250,10 +257,12 @@ subroutine FatesParamsInit()
     ED_val_patch_fusion_tol               = nan
     ED_val_canopy_closure_thresh          = nan
     stomatal_model                        = -9
+    stomatal_assim_model                  = -9
     hydr_kmax_rsurf1                      = nan
     hydr_kmax_rsurf2                      = nan
     hydr_psi0                             = nan
     hydr_psicap                           = nan
+    hydr_solver                           = -9
     bgc_soil_salinity                     = nan
     logging_dbhmin                        = nan
     logging_dbhmax                        = nan
@@ -383,7 +392,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=stomatal_assim_name, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf1, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -561,6 +573,10 @@ subroutine FatesReceiveParams(fates_params)
          data=tmpreal)
     stomatal_model = nint(tmpreal)
 
+    call fates_params%RetreiveParameter(name=stomatal_assim_name, &
+         data=tmpreal)
+    stomatal_assim_model = nint(tmpreal)
+
     call fates_params%RetreiveParameter(name=hydr_name_kmax_rsurf1, &
           data=hydr_kmax_rsurf1)
 
@@ -572,7 +588,11 @@ subroutine FatesReceiveParams(fates_params)
 
     call fates_params%RetreiveParameter(name=hydr_name_psicap, &
           data=hydr_psicap)
-
+
+    call fates_params%RetreiveParameter(name=hydr_name_solver, &
+         data=tmpreal)
+    hydr_solver = nint(tmpreal)
+
     call fates_params%RetreiveParameter(name=bgc_name_soil_salinity, &
           data=bgc_soil_salinity)	  
 
@@ -689,11 +709,13 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'ED_val_cohort_age_fusion_tol = ',ED_val_cohort_age_fusion_tol
         write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol
         write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh
-        write(fates_log(),fmt0) 'stomatal_model = ',stomatal_model      
+        write(fates_log(),fmt0) 'stomatal_model = ',stomatal_model
+        write(fates_log(),fmt0) 'stomatal_assim_model = ',stomatal_assim_model            
         write(fates_log(),fmt0) 'hydr_kmax_rsurf1 = ',hydr_kmax_rsurf1
         write(fates_log(),fmt0) 'hydr_kmax_rsurf2 = ',hydr_kmax_rsurf2  
         write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0
         write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap
+        write(fates_log(),fmt0) 'hydr_solver = ',hydr_solver
         write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity
         write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin
         write(fates_log(),fmt0) 'logging_dbhmax = ',logging_dbhmax

diff --git a/main/FatesHydraulicsMemMod.F90 b/main/FatesHydraulicsMemMod.F90
@@ -14,8 +14,8 @@ module FatesHydraulicsMemMod
    ! Define the various different solver options for hydraulics
 
    integer, parameter, public :: hydr_solver_1DTaylor = 1
-   integer, parameter, public :: hydr_solver_2DNewton = 2
-   integer, parameter, public :: hydr_solver_2DPicard = 3
+   integer, parameter, public :: hydr_solver_2DNewton = 3
+   integer, parameter, public :: hydr_solver_2DPicard = 2
 
    ! Number of soil layers for indexing cohort fine root quanitities
    ! NOTE: The hydraulics code does have some capacity to run a single soil

diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl
@@ -201,9 +201,9 @@ variables:
 	double fates_fire_bark_scaler(fates_pft) ;
 		fates_fire_bark_scaler:units = "fraction" ;
 		fates_fire_bark_scaler:long_name = "the thickness of a cohorts bark as a fraction of its dbh" ;
-	double fates_fire_crown_depth_frac(fates_pft) ;
-		fates_fire_crown_depth_frac:units = "fraction" ;
-		fates_fire_crown_depth_frac:long_name = "the depth of a cohorts crown as a fraction of its height" ;
+	double fates_crown_depth_frac(fates_pft) ;
+		fates_crown_depth_frac:units = "fraction" ;
+		fates_crown_depth_frac:long_name = "the depth of a cohorts crown as a fraction of its height" ;
 	double fates_fire_crown_kill(fates_pft) ;
 		fates_fire_crown_kill:units = "NA" ;
 		fates_fire_crown_kill:long_name = "fire parameter, see equation 22 in Thonicke et al 2010" ;
@@ -656,6 +656,9 @@ variables:
 	double fates_hydr_psicap ;
 		fates_hydr_psicap:units = "MPa" ;
 		fates_hydr_psicap:long_name = "sapwood water potential at which capillary reserves exhausted" ;
+   double fates_hydr_solver ;
+	   fates_hydr_solver:units = "unitless" ;
+		fates_hydr_solver:long_name = "switch designating which numerical solver for plant hydraulics, 1 = 1D taylor, 2 = 2D Picard" ;
 	double fates_init_litter ;
 		fates_init_litter:units = "NA" ;
 		fates_init_litter:long_name = "Initialization value for litter pool in cold-start (NOT USED)" ;
@@ -743,6 +746,9 @@ variables:
 	double fates_soil_salinity ;
 		fates_soil_salinity:units = "ppt" ;
 		fates_soil_salinity:long_name = "soil salinity used for model when not coupled to dynamic soil salinity" ;
+	double fates_stomatal_assim ;
+	   fates_stomatal_assim:units = "unitless" ;
+		fates_stomatal_assim:long_name = "a switch designating whether to use net (1) or gross (2) assimilation in the stomatal model" ;
 	double fates_theta_cj_c3 ;
 		fates_theta_cj_c3:units = "unitless" ;
 		fates_theta_cj_c3:long_name = "Empirical curvature parameter for ac, aj photosynthesis co-limitation in c3 plants" ;
@@ -939,7 +945,7 @@ data:
  fates_fire_bark_scaler = 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07,
     0.07, 0.07, 0.07, 0.07 ;
 
- fates_fire_crown_depth_frac = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.95, 0.95,
+ fates_crown_depth_frac = 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.95, 0.95,
     0.95, 1, 1, 1 ;
 
  fates_fire_crown_kill = 0.775, 0.775, 0.775, 0.775, 0.775, 0.775, 0.775,
@@ -1421,6 +1427,8 @@ data:
 
  fates_hydr_psicap = -0.6 ;
 
+ fates_hydr_solver = 1 ;
+
  fates_init_litter = 0.05 ;
 
  fates_leaf_stomatal_model = 1 ;
@@ -1479,6 +1487,8 @@ data:
 
  fates_soil_salinity = 0.4 ;
 
+ fates_stomatal_assim = 1 ;
+
  fates_theta_cj_c3 = 0.999 ;
 
  fates_theta_cj_c4 = 0.999 ;

diff --git a/parteh/PRTParametersMod.F90 b/parteh/PRTParametersMod.F90
@@ -101,7 +101,7 @@ module PRTParametersMod
      real(r8), allocatable :: c2b(:)                        ! Carbon to biomass multiplier [kg/kgC]
      real(r8), allocatable :: wood_density(:)               ! wood density  g cm^-3  ...
      real(r8), allocatable :: woody(:)                      ! Does the plant have wood?      (1=yes, 0=no)
-     real(r8), allocatable :: crown(:)                      ! fraction of the height of the plant
+     real(r8), allocatable :: crown_depth_frac(:)           ! fraction of the height of the plant
                                                             ! that is occupied by crown
      real(r8), allocatable :: slamax(:)                     ! Maximum specific leaf area of plant (at bottom) [m2/gC]
      real(r8), allocatable :: slatop(:)                     ! Specific leaf area at canopy top [m2/gC]

diff --git a/parteh/PRTParamsFATESMod.F90 b/parteh/PRTParamsFATESMod.F90
@@ -179,7 +179,7 @@ subroutine PRTRegisterPFT(fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
-    name = 'fates_fire_crown_depth_frac'
+    name = 'fates_crown_depth_frac'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
@@ -435,9 +435,9 @@ subroutine PRTReceivePFT(fates_params)
     call fates_params%RetreiveParameterAllocate(name=name, &
          data=prt_params%fnrt_prof_mode)
 
-    name = 'fates_fire_crown_depth_frac'
+    name = 'fates_crown_depth_frac'
     call fates_params%RetreiveParameterAllocate(name=name, &
-         data=prt_params%crown)
+         data=prt_params%crown_depth_frac)
 
     name = 'fates_woody'
     call fates_params%RetreiveParameterAllocate(name=name, &
@@ -911,7 +911,7 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'prt_phos_stoich_p2 = ',prt_params%phos_stoich_p2
         write(fates_log(),fmt0) 'prt_alloc_priority = ',prt_params%alloc_priority
         write(fates_log(),fmt0) 'woody = ',prt_params%woody
-        write(fates_log(),fmt0) 'crown = ',prt_params%crown
+        write(fates_log(),fmt0) 'crown_depth_frac = ',prt_params%crown_depth_frac
         write(fates_log(),fmt0) 'roota_par = ',prt_params%fnrt_prof_a
         write(fates_log(),fmt0) 'rootb_par = ',prt_params%fnrt_prof_b
         write(fates_log(),fmt0) 'fnrt_prof_mode = ',prt_params%fnrt_prof_mode