NGEET · rgknox · Jan 16, 2024 · May 15, 2023 · May 15, 2023 · May 16, 2023
diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90
diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90
@@ -2773,8 +2773,9 @@ subroutine fuse_2_patches(csite, dp, rp)
     rp%zstar                = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area
     rp%c_stomata            = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
     rp%c_lblayer            = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area
-    rp%radiation_error      = (dp%radiation_error*dp%area + rp%radiation_error*rp%area) * inv_sum_area
-
+    rp%rad_error(1)         = (dp%rad_error(1)*dp%area + rp%rad_error(1)*rp%area) * inv_sum_area
+    rp%rad_error(2)         = (dp%rad_error(2)*dp%area + rp%rad_error(2)*rp%area) * inv_sum_area
+
     rp%area = rp%area + dp%area !THIS MUST COME AT THE END!
 
     !insert donor cohorts into recipient patch

diff --git a/biogeochem/FatesAllometryMod.F90 b/biogeochem/FatesAllometryMod.F90
@@ -92,12 +92,13 @@ module FatesAllometryMod
   use FatesConstantsMod, only : calloc_abs_error
   use FatesConstantsMod, only : fates_unset_r8
   use FatesConstantsMod, only : itrue
+  use FatesConstantsMod, only : nearzero
   use shr_log_mod      , only : errMsg => shr_log_errMsg
   use FatesGlobals     , only : fates_log
   use FatesGlobals     , only : endrun => fates_endrun
   use FatesGlobals     , only : FatesWarn,N2S,A2S,I2S
-  use EDParamsMod      , only : nlevleaf, dinc_vai
-  use EDParamsMod      , only : nclmax
+  use EDParamsMod      , only : nlevleaf,dinc_vai,dlower_vai
+  use EDParamsMod       , only : nclmax
   use DamageMainMod    , only : GetCrownReduction
 
   implicit none
@@ -125,7 +126,8 @@ module FatesAllometryMod
   public :: set_root_fraction  ! Generic wrapper to calculate normalized
                                ! root profiles
   public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode
-
+  public :: VegAreaLayer
+
   logical         , parameter :: verbose_logging = .false.
   character(len=*), parameter :: sourcefile = __FILE__
 
@@ -2436,7 +2438,8 @@ real(r8) function decay_coeff_kn(pft,vcmax25top)
   end function decay_coeff_kn
 
   ! =====================================================================================
-subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )
+
+  subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )
 
      ! =========================================================================
      ! This subroutine estimates the diameter based on either the structural biomass
@@ -2586,6 +2589,182 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d
      return
   end subroutine ForceDBH
 
+  ! =========================================================================
+
+  subroutine VegAreaLayer(tree_lai,tree_sai,tree_height,iv,nv,pft,snow_depth, & 
+       vai_top,vai_bot, elai_layer,esai_layer,tlai_layer,tsai_layer)
+
+    ! -----------------------------------------------------------------------------------
+    ! This routine returns the exposed leaf and stem areas (m2 of leaf and stem) per m2 of
+    ! ground inside the crown, for the leaf-layer specified.
+    ! -----------------------------------------------------------------------------------
+
+    real(r8),intent(in) :: tree_lai         ! the in-crown leaf area index for the plant
+                                            ! [m2 leaf/m2 crown footprint]
+    real(r8),intent(in) :: tree_sai         ! the in-crown stem area index for the plant
+                                            ! [m2 stem/m2 crown footprint]
+    real(r8),intent(in) :: tree_height      ! the height of the plant [m]
+    integer,intent(in)  :: iv               ! vegetation layer index
+    integer,intent(in)  :: nv               ! this plants total number of veg layers
+    integer,intent(in)  :: pft              ! plant functional type index
+    real(r8),intent(in) :: snow_depth       ! the depth of snow on the ground [m]
+    real(r8),intent(out) :: vai_top
+    real(r8),intent(out) :: vai_bot          ! the VAI of the bin top and bottom
+    real(r8),intent(out) :: elai_layer       ! exposed leaf area index of the layer
+    real(r8),intent(out) :: esai_layer       ! exposed stem area index of the layer
+    real(r8),optional,intent(out) :: tlai_layer       ! total leaf area index of the layer
+    real(r8),optional,intent(out) :: tsai_layer       ! total stem area index of the layer
+
+                                 ! [m2 of leaf in bin / m2 crown footprint]
+    real(r8) :: tree_vai         ! the in-crown veg area index for the plant
+    real(r8) :: fraction_exposed ! fraction of the veg media that is above snow
+    real(r8) :: layer_top_height ! Physical height of the layer top relative to ground [m]
+    real(r8) :: layer_bot_height ! Physical height of the layer bottom relative to ground [m]
+    real(r8) :: tlai,tsai        ! temporary total area indices [m2/m2]
+    real(r8) :: fleaf            ! fraction of biomass in layer that is leaf
+    real(r8) :: remainder        ! old-method: remainder of biomass in last bin
+    integer, parameter :: layer_height_const_depth = 1 ! constant physical depth assumption
+    integer, parameter :: layer_height_const_lad   = 2 ! constant leaf area depth assumption
+    integer, parameter :: layer_height_method = layer_height_const_depth
+
+
+    logical, parameter :: preserve_b4b = .true.
+
+
+    if_preserve_b4b: if(preserve_b4b) then
+
+       if( (tree_lai+tree_sai) > nearzero)then
+
+          ! lai = tree_lai * currentCohort%c_area/currentPatch%total_canopy_area
+          ! sai = tree_sai * currentCohort%c_area/currentPatch%total_canopy_area
+
+          ! See issue: https://github.com/NGEET/fates/issues/899
+          ! fleaf = currentCohort%treelai / (currentCohort%treelai + currentCohort%treesai)
+          fleaf = tree_lai / (tree_lai+tree_sai)
+       else
+          fleaf = 0._r8
+       endif
+
+       if(iv==0)then  ! Whole plant (not just specific layer)
+
+          layer_top_height = tree_height
+          layer_bot_height = tree_height - ( tree_height * prt_params%crown_depth_frac(pft) )
+
+       else
+
+          layer_top_height = tree_height - ( real(iv-1,r8)/nv * tree_height *  &
+               prt_params%crown_depth_frac(pft) )
+
+          layer_bot_height = tree_height - ( real(iv,r8)/nv * tree_height * &
+               prt_params%crown_depth_frac(pft) )
+       end if
+
+       fraction_exposed = 1.0_r8
+       if(snow_depth  > layer_top_height)then
+          fraction_exposed = 0._r8
+       endif
+       if(snow_depth < layer_bot_height)then
+          fraction_exposed = 1._r8
+       endif
+       if(snow_depth >= layer_bot_height .and. &
+            snow_depth <= layer_top_height) then !only partly hidden...
+          fraction_exposed =  1._r8 - max(0._r8,(min(1.0_r8,(snow_depth -layer_bot_height)/ &
+               (layer_top_height-layer_bot_height ))))
+       endif
+
+       if(iv==0) then
+          remainder = tree_lai+tree_sai
+       elseif(iv==nv) then
+          remainder = (tree_lai + tree_sai) - &
+               (dlower_vai(iv) - dinc_vai(iv))
+          if(remainder > dinc_vai(iv) )then
+             write(fates_log(), *)'ED: issue with remainder', &
+                  tree_lai,tree_sai,dinc_vai(iv), & 
+                  nv,remainder
+             call endrun(msg=errMsg(sourcefile, __LINE__))
+          endif
+       else
+          remainder = dinc_vai(iv)
+       end if
+
+       if(present(tlai_layer)) tlai_layer = remainder * fleaf  !* currentCohort%c_area/currentPatch%total_canopy_area
+       if(present(tsai_layer)) tsai_layer = remainder * (1._r8-fleaf) !* currentCohort%c_area/currentPatch%total_canopy_area
+       elai_layer = fraction_exposed * remainder * fleaf
+       esai_layer = fraction_exposed * remainder * (1._r8-fleaf)
+
+
+       return
+
+    end if if_preserve_b4b
+
+
+    tree_vai = tree_lai + tree_sai
+
+    if_any_vai: if(tree_vai>0._r8)then
+
+       if(iv==0)then
+          vai_top = 0.0
+          vai_bot = tree_vai
+       else
+
+          if(iv>1)then
+             vai_top = dlower_vai(iv) - dinc_vai(iv)
+          else
+             vai_top = 0._r8
+          end if
+
+          if(iv<nv) then
+             vai_bot = dlower_vai(iv)
+          else
+             vai_bot = tree_vai
+          end if
+       end if
+
+       if(layer_height_method .eq. layer_height_const_depth)then
+          if(iv==0)then
+             layer_top_height = tree_height
+             layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft))
+          else
+             layer_top_height = tree_height*(1._r8 - real(iv-1,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
+             layer_bot_height = tree_height*(1._r8 - real(iv,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
+          end if
+       else
+          layer_top_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_top/tree_vai)
+          layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_bot/tree_vai)
+       end if
+
+       fraction_exposed =  1._r8 - max(0._r8,(min(1._r8, (snow_depth-layer_bot_height)/(layer_top_height-layer_bot_height))))
+
+       tlai = (vai_bot-vai_top) * tree_lai / tree_vai
+       tsai = (vai_bot-vai_top) * tree_sai / tree_vai
+
+       if(present(tlai_layer)) tlai_layer = tlai
+       if(present(tsai_layer)) tsai_layer = tsai
+
+       elai_layer = fraction_exposed * tlai
+       esai_layer = fraction_exposed * tsai
+
+       ! Update the vai at the bottom to be removed/decreased if there is no exposure
+       ! set the vai top and bottoms to the snow layer if below
+       !vai_top = min(vai_top,fraction_exposed*tree_vai)
+
+       vai_bot = vai_top + fraction_exposed*(vai_bot-vai_top)
+
+    else
+
+       if(present(tlai_layer)) tlai_layer = 0._r8
+       if(present(tsai_layer)) tsai_layer = 0._r8
+       elai_layer = 0._r8
+       esai_layer = 0._r8
+       vai_bot = 0._r8
+       vai_top = 0._r8
+
+    end if if_any_vai
+
+
+    return
+  end subroutine VegAreaLayer
+
   ! =========================================================================
 
   subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx)

diff --git a/biogeochem/FatesCohortMod.F90 b/biogeochem/FatesCohortMod.F90
@@ -200,6 +200,8 @@ module FatesCohortMod
 
     !---------------------------------------------------------------------------
 
+    integer :: twostr_col  ! The column index in the two-stream solution that this cohort is part of
+
     ! RESPIRATION COMPONENTS
     real(r8) :: rdark            ! dark respiration [kgC/indiv/s]
     real(r8) :: resp_g_tstep     ! growth respiration [kgC/indiv/timestep]