Merge pull request #752 from ckoven/variable_vai_bins

Variable VAI bin widths

biogeochem/EDCanopyStructureMod.F90

@@ -1476,8 +1476,8 @@ subroutine leaf_area_profile( currentSite )
 
     ! !USES:
 
-    use EDtypesMod           , only : area, dinc_ed, hitemax, n_hite_bins
-
+    use EDtypesMod           , only : area, dinc_vai, dlower_vai, hitemax, n_hite_bins
+  
     !
     ! !ARGUMENTS
     type(ed_site_type)     , intent(inout) :: currentSite
@@ -1508,8 +1508,6 @@ subroutine leaf_area_profile( currentSite )
 
     !----------------------------------------------------------------------
 
-
-
     smooth_leaf_distribution = 0
 
     ! Here we are trying to generate a profile of leaf area, indexed by 'z' and by pft
@@ -1572,7 +1570,7 @@ subroutine leaf_area_profile( currentSite )
              currentCohort%sai =  currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area
 
              ! Number of actual vegetation layers in this cohort's crown
-             currentCohort%nv =  ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)
+             currentCohort%nv =  count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
 
              currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV)
 
@@ -1736,15 +1734,16 @@ subroutine leaf_area_profile( currentSite )
 
                    if(iv==currentCohort%NV) then
                       remainder = (currentCohort%treelai + currentCohort%treesai) - &
-                           (dinc_ed*real(currentCohort%nv-1,r8))
-                      if(remainder > dinc_ed )then
+                           (dlower_vai(iv) - dinc_vai(iv))
+                      if(remainder > dinc_vai(iv) )then
                          write(fates_log(), *)'ED: issue with remainder', &
-                              currentCohort%treelai,currentCohort%treesai,dinc_ed, &
+                              currentCohort%treelai,currentCohort%treesai,dinc_vai(iv), & 
                               currentCohort%NV,remainder
+
                          call endrun(msg=errMsg(sourcefile, __LINE__))
                       endif
                    else
-                      remainder = dinc_ed
+                      remainder = dinc_vai(iv)
                    end if
 
                    currentPatch%tlai_profile(cl,ft,iv) = currentPatch%tlai_profile(cl,ft,iv) + &

biogeochem/EDPhysiologyMod.F90

@@ -37,8 +37,8 @@ module EDPhysiologyMod
   use EDTypesMod          , only : site_massbal_type
   use EDTypesMod          , only : numlevsoil_max
   use EDTypesMod          , only : numWaterMem
-  use EDTypesMod          , only : dl_sf, dinc_ed, area_inv
-  use EDTypesMod                , only : AREA
+  use EDTypesMod          , only : dl_sf, dinc_vai, dlower_vai, area_inv
+  use EDTypesMod          , only : AREA
   use FatesLitterMod      , only : ncwd
   use FatesLitterMod      , only : ndcmpy
   use FatesLitterMod      , only : ilabile
@@ -388,7 +388,7 @@ subroutine trim_canopy( currentSite )
     real(r8) :: sapw_c                ! sapwood carbon [kg]
     real(r8) :: store_c               ! storage carbon [kg]
     real(r8) :: struct_c              ! structure carbon [kg]
-    real(r8) :: leaf_inc              ! LAI-only portion of the vegetation increment of dinc_ed
+    real(r8) :: leaf_inc              ! LAI-only portion of the vegetation increment of dinc_vai
     real(r8) :: lai_canopy_above      ! the LAI in the canopy layers above the layer of interest
     real(r8) :: lai_layers_above      ! the LAI in the leaf layers, within the current canopy,
     ! above the leaf layer of interest
@@ -472,7 +472,7 @@ subroutine trim_canopy( currentSite )
                currentPatch%canopy_layer_tlai, currentCohort%treelai, &
                currentCohort%vcmax25top,0 )
 
-          currentCohort%nv      = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)
+          currentCohort%nv      = count((currentCohort%treelai+currentCohort%treesai) .gt. dlower_vai(:)) + 1
 
           if (currentCohort%nv > nlevleaf)then
              write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv, &
@@ -503,12 +503,12 @@ subroutine trim_canopy( currentSite )
           do z = 1, currentCohort%nv
 
              ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves)
-
-             leaf_inc    = dinc_ed * &
+             leaf_inc    = dinc_vai(z) * &
                   currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
-
+             
              ! Now calculate the cumulative top-down lai of the current layer's midpoint within the current cohort
-             lai_layers_above      = leaf_inc * (z-1)
+             lai_layers_above      = (dlower_vai(z) - dinc_vai(z)) * &
+                  currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
              lai_current           = min(leaf_inc, currentCohort%treelai - lai_layers_above)
              cumulative_lai_cohort = lai_layers_above + 0.5*lai_current
 

biogeochem/FatesAllometryMod.F90

@@ -95,7 +95,7 @@ module FatesAllometryMod
   use shr_log_mod      , only : errMsg => shr_log_errMsg
   use FatesGlobals     , only : fates_log
   use FatesGlobals     , only : endrun => fates_endrun
-  use EDTypesMod       , only : nlevleaf, dinc_ed
+  use EDTypesMod       , only : nlevleaf, dinc_vai
   use EDTypesMod       , only : nclmax
 
 
@@ -731,7 +731,7 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, &
 
     tree_sai   =  prt_params%allom_sai_scaler(pft) * target_lai
 
-    if( (treelai + tree_sai) > (nlevleaf*dinc_ed) )then
+    if( (treelai + tree_sai) > (sum(dinc_vai)) )then
 
        call h_allom(dbh,pft,h)
 
@@ -740,7 +740,8 @@ real(r8) function tree_sai( pft, dbh, canopy_trim, c_area, nplant, cl, &
        write(fates_log(),*) 'sai: ',tree_sai
        write(fates_log(),*) 'target_lai: ',target_lai
        write(fates_log(),*) 'lai+sai: ',treelai+tree_sai
-       write(fates_log(),*) 'nlevleaf,dinc_ed,nlevleaf*dinc_ed :',nlevleaf,dinc_ed,nlevleaf*dinc_ed
+       write(fates_log(),*) 'dinc_vai:',dinc_vai
+       write(fates_log(),*) 'nlevleaf,sum(dinc_vai):',nlevleaf,sum(dinc_vai)
        write(fates_log(),*) 'pft: ',pft
        write(fates_log(),*) 'call id: ',call_id
        write(fates_log(),*) 'n: ',nplant

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -117,7 +117,8 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     use EDTypesMod        , only : ed_cohort_type
     use EDTypesMod        , only : ed_site_type
     use EDTypesMod        , only : maxpft
-    use EDTypesMod        , only : dinc_ed
+    use EDTypesMod        , only : dinc_vai
+    use EDTypesMod        , only : dlower_vai
     use FatesInterfaceTypesMod , only : bc_in_type
     use FatesInterfaceTypesMod , only : bc_out_type
     use EDCanopyStructureMod, only : calc_areaindex
@@ -227,9 +228,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     ! over each cohort x layer.
     real(r8) :: cohort_eleaf_area  ! This is the effective leaf area [m2] reported by each cohort
     real(r8) :: lnc_top            ! Leaf nitrogen content per unit area at canopy top [gN/m2]
-    real(r8) :: lmr25top           ! canopy top leaf maint resp rate at 25C
-    ! for this plant or pft (umol CO2/m**2/s)
-    real(r8) :: leaf_inc           ! LAI-only portion of the vegetation increment of dinc_ed
+    real(r8) :: lmr25top           ! canopy top leaf maint resp rate at 25C 
+                                   ! for this plant or pft (umol CO2/m**2/s)
+    real(r8) :: leaf_inc           ! LAI-only portion of the vegetation increment of dinc_vai
     real(r8) :: lai_canopy_above   ! the LAI in the canopy layers above the layer of interest
     real(r8) :: lai_layers_above   ! the LAI in the leaf layers, within the current canopy,
     ! above the leaf layer of interest
@@ -419,31 +420,33 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                  (hlm_parteh_mode .ne. prt_carbon_allom_hyp )   ) then
 
                                if (hlm_use_planthydro.eq.itrue ) then
-
+                                  
                                   stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentCohort%co_hydr%btran )
-                                  btran_eff = currentCohort%co_hydr%btran
+                                  btran_eff = currentCohort%co_hydr%btran 
 
-                                  ! dinc_ed is the total vegetation area index of each "leaf" layer
+                                  ! dinc_vai(:) is the total vegetation area index of each "leaf" layer
                                   ! we convert to the leaf only portion of the increment
                                   ! ------------------------------------------------------
-                                  leaf_inc    = dinc_ed * &
+                                  leaf_inc    = dinc_vai(iv) * &
                                        currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
 
                                   ! Now calculate the cumulative top-down lai of the current layer's midpoint
-                                  lai_canopy_above  = sum(currentPatch%canopy_layer_tlai(1:cl-1))
-                                  lai_layers_above  = leaf_inc * (iv-1)
+                                  lai_canopy_above  = sum(currentPatch%canopy_layer_tlai(1:cl-1)) 
+
+                                  lai_layers_above  = (dlower_vai(iv) - dinc_vai(iv)) * &
+                                       currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
                                   lai_current       = min(leaf_inc, currentCohort%treelai - lai_layers_above)
-                                  cumulative_lai    = lai_canopy_above + lai_layers_above + 0.5*lai_current
+                                  cumulative_lai    = lai_canopy_above + lai_layers_above + 0.5*lai_current 
 
                                   leaf_psi = currentCohort%co_hydr%psi_ag(1)
 
                                else
 
-                                  stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) )
+                                  stomatal_intercept_btran = max( cf/rsmax0,stomatal_intercept(ft)*currentPatch%btran_ft(ft) ) 
 
                                   btran_eff = currentPatch%btran_ft(ft)
                                   ! For consistency sake, we use total LAI here, and not exposed
-                                  ! if the plant is under-snow, it will be effectively dormant for
+                                  ! if the plant is under-snow, it will be effectively dormant for 
                                   ! the purposes of nscaler
 
                                   cumulative_lai = sum(currentPatch%canopy_layer_tlai(1:cl-1))  + &
@@ -460,7 +463,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
 
 
                                ! Bonan et al (2011) JGR, 116, doi:10.1029/2010JG001593 used
-                               ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al
+                               ! kn = 0.11. Here, derive kn from vcmax25 as in Lloyd et al 
                                ! (2010) Biogeosciences, 7, 1833-1859
 
                                kn = decay_coeff_kn(ft,currentCohort%vcmax25top)
@@ -479,7 +482,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                select case(hlm_parteh_mode)
                                case (prt_carbon_allom_hyp)
 
-                                 lnc_top  = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft)
+                                  lnc_top  = prt_params%nitr_stoich_p1(ft,prt_params%organ_param_id(leaf_organ))/slatop(ft)
 
                                case (prt_cnp_flex_allom_hyp)
 

main/EDParamsMod.F90

@@ -39,6 +39,8 @@ module EDParamsMod
 
    real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
    real(r8),protected, public :: ED_val_comp_excln
+   real(r8),protected, public :: ED_val_vai_top_bin_width
+   real(r8),protected, public :: ED_val_vai_width_increase_factor
    real(r8),protected, public :: ED_val_init_litter
    real(r8),protected, public :: ED_val_nignitions
    real(r8),protected, public :: ED_val_understorey_death
@@ -91,14 +93,14 @@ module EDParamsMod
    ! 1  = Christofferson et al. 2016 (TFS),   2 = Van Genuchten 1980
    integer, protected,allocatable,public :: hydr_htftype_node(:) 
 
-   character(len=param_string_length),parameter,public :: ED_name_vai_top_bin_width = "fates_vai_top_bin_width"
-   character(len=param_string_length),parameter,public :: ED_name_vai_width_increase_factor = "fates_vai_width_increase_factor"
    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"
    character(len=param_string_length),parameter,public :: ED_name_hydr_htftype_node = "fates_hydr_htftype_node"
    character(len=param_string_length),parameter,public :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac"
    character(len=param_string_length),parameter,public :: ED_name_comp_excln = "fates_comp_excln"
+   character(len=param_string_length),parameter,public :: ED_name_vai_top_bin_width = "fates_vai_top_bin_width"
+   character(len=param_string_length),parameter,public :: ED_name_vai_width_increase_factor = "fates_vai_width_increase_factor"
    character(len=param_string_length),parameter,public :: ED_name_init_litter = "fates_init_litter"
    character(len=param_string_length),parameter,public :: ED_name_nignitions = "fates_fire_nignitions"
    character(len=param_string_length),parameter,public :: ED_name_understorey_death = "fates_mort_understorey_death"
@@ -216,6 +218,8 @@ subroutine FatesParamsInit()
     maintresp_model                       = -9
     fates_mortality_disturbance_fraction  = nan
     ED_val_comp_excln                     = nan
+    ED_val_vai_top_bin_width              = nan
+    ED_val_vai_width_increase_factor      = nan
     ED_val_init_litter                    = nan
     ED_val_nignitions                     = nan
     ED_val_understorey_death              = nan
@@ -283,12 +287,6 @@ subroutine FatesRegisterParams(fates_params)
 
     call FatesParamsInit()
 
-    call fates_params%RegisterParameter(name=ED_name_vai_top_bin_width, dimension_shape=dimension_shape_scalar, &
-         dimension_names=dim_names_scalar)
-
-    call fates_params%RegisterParameter(name=ED_name_vai_width_increase_factor, dimension_shape=dimension_shape_scalar, &
-         dimension_names=dim_names_scalar)
-
     call fates_params%RegisterParameter(name=ED_name_photo_temp_acclim_timescale, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -310,6 +308,12 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
+    call fates_params%RegisterParameter(name=ED_name_vai_top_bin_width, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
+    call fates_params%RegisterParameter(name=ED_name_vai_width_increase_factor, dimension_shape=dimension_shape_scalar, &
+         dimension_names=dim_names_scalar)
+
     call fates_params%RegisterParameter(name=ED_name_init_litter, dimension_shape=dimension_shape_scalar, &
          dimension_names=dim_names_scalar)
 
@@ -463,12 +467,6 @@ subroutine FatesReceiveParams(fates_params)
     real(r8) :: tmpreal ! local real variable for changing type on read
     real(r8), allocatable :: hydr_htftype_real(:)
 
-    call fates_params%RetreiveParameter(name=ED_name_vai_top_bin_width, &
-         data=vai_top_bin_width)
-
-    call fates_params%RetreiveParameter(name=ED_name_vai_width_increase_factor, &
-         data=vai_width_increase_factor)
-
     call fates_params%RetreiveParameter(name=ED_name_photo_temp_acclim_timescale, &
          data=photo_temp_acclim_timescale)
 
@@ -486,6 +484,12 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=ED_name_comp_excln, &
          data=ED_val_comp_excln)
 
+    call fates_params%RetreiveParameter(name=ED_name_vai_top_bin_width, &
+         data=ED_val_vai_top_bin_width)
+
+    call fates_params%RetreiveParameter(name=ED_name_vai_width_increase_factor, &
+         data=ED_val_vai_width_increase_factor)
+
     call fates_params%RetreiveParameter(name=ED_name_init_litter, &
          data=ED_val_init_litter)
 
@@ -654,6 +658,8 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmti) 'hydr_htftype_node = ',hydr_htftype_node
         write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
         write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln
+        write(fates_log(),fmt0) 'ED_val_vai_top_bin_width = ',ED_val_vai_top_bin_width
+        write(fates_log(),fmt0) 'ED_val_vai_width_increase_factor = ',ED_val_vai_width_increase_factor
         write(fates_log(),fmt0) 'ED_val_init_litter = ',ED_val_init_litter
         write(fates_log(),fmt0) 'ED_val_nignitions = ',ED_val_nignitions
         write(fates_log(),fmt0) 'ED_val_understorey_death = ',ED_val_understorey_death

main/EDTypesMod.F90

@@ -18,6 +18,7 @@ module EDTypesMod
   use FatesLitterMod,        only : ncwd
   use FatesConstantsMod,     only : n_anthro_disturbance_categories
   use FatesConstantsMod,     only : days_per_year
+  use FatesConstantsMod,     only : fates_unset_r8
   use FatesRunningMeanMod,   only : rmean_type
   use FatesInterfaceTypesMod,only : bc_in_type
   use FatesInterfaceTypesMod,only : bc_out_type
@@ -38,7 +39,6 @@ module EDTypesMod
                                                           ! to understory layers (all layers that
                                                           ! are not the top canopy layer)
 
-  integer, parameter, public :: nlevleaf = 30             ! number of leaf layers in canopy layer
   integer, parameter, public :: maxpft = 16               ! maximum number of PFTs allowed
                                                           ! the parameter file may determine that fewer
                                                           ! are used, but this helps allocate scratch
@@ -60,6 +60,10 @@ module EDTypesMod
   integer, parameter, public :: idirect   = 1             ! This is the array index for direct radiation
   integer, parameter, public :: idiffuse  = 2             ! This is the array index for diffuse radiation
 
+  ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code
+  integer, parameter, public :: nlevleaf = 30             ! number of leaf+stem layers in canopy layer
+  real(r8), public :: dinc_vai(nlevleaf)   = fates_unset_r8 ! VAI bin widths array
+  real(r8), public :: dlower_vai(nlevleaf) = fates_unset_r8 ! lower edges of VAI bins
 
   ! TODO: we use this cp_maxSWb only because we have a static array q(size=2) of
   ! land-ice abledo for vis and nir.  This should be a parameter, which would
@@ -125,8 +129,6 @@ module EDTypesMod
 
   ! BIOLOGY/BIOGEOCHEMISTRY        
   integer , parameter, public :: num_vegtemp_mem      = 10         ! Window of time over which we track temp for cold sensecence (days)
-  real(r8), parameter, public :: dinc_ed              = 1.0_r8     ! size of VAI bins (LAI+SAI)  [CHANGE THIS NAME WITH NEXT INTERFACE
-                                                           ! UPDATE]
   integer , parameter, public :: N_DIST_TYPES         = 3          ! Disturbance Modes 1) tree-fall, 2) fire, 3) logging
   integer , parameter, public :: dtype_ifall          = 1          ! index for naturally occuring tree-fall generated event
   integer , parameter, public :: dtype_ifire          = 2          ! index for fire generated disturbance event