Merge pull request #402 from kovenock/kovenock_slaprofile

SLA Profile with Canopy Depth

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -1016,8 +1016,10 @@ subroutine canopy_summarization( nsites, sites, bc_in )
 
              call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,&
                   currentCohort%pft,currentCohort%c_area)
-             currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, &
-                  currentCohort%pft, currentCohort%c_area, currentCohort%n )
+
+             currentCohort%treelai = tree_lai(currentCohort%bl,             &
+                  currentCohort%pft, currentCohort%c_area, currentCohort%n, &
+                  currentCohort%canopy_layer, currentPatch%canopy_layer_tlai )
 
              canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
 
@@ -1087,7 +1089,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
     !
     ! The following patch level diagnostics are updated here:
     ! 
-    ! currentPatch%canopy_layer_tai(cl)    ! TAI of each canopy layer
+    ! currentPatch%canopy_layer_tlai(cl)   ! total leaf area index of canopy layer
     ! currentPatch%ncan(cl,ft)             ! number of vegetation layers needed
     !                                      ! in this patch's pft/canopy-layer 
     ! currentPatch%nrad(cl,ft)             ! same as ncan, but does not include
@@ -1157,7 +1159,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
        ! calculate tree lai and sai.
        ! --------------------------------------------------------------------------------
 
-       currentPatch%canopy_layer_tai(:)         = 0._r8
+       currentPatch%canopy_layer_tlai(:)        = 0._r8
        currentPatch%ncan(:,:)                   = 0 
        currentPatch%nrad(:,:)                   = 0 
        patch_lai                                = 0._r8
@@ -1176,38 +1178,38 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
 
        if (currentPatch%total_canopy_area > nearzero ) then
 
-       currentCohort => currentPatch%shortest
+
+       currentCohort => currentPatch%tallest
        do while(associated(currentCohort)) 
 
           ft = currentCohort%pft
           cl = currentCohort%canopy_layer
 
-          currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, currentCohort%pft, &
-               currentCohort%c_area, currentCohort%n )
-          currentCohort%treesai = tree_sai(currentCohort%dbh, currentCohort%pft, currentCohort%canopy_trim, &
-               currentCohort%c_area, currentCohort%n)
+          ! Calculate LAI of layers above
+          ! Note that the canopy_layer_lai is also calculated in this loop
+          ! but since we go top down in terms of plant size, we should be okay
+
+          currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%pft, currentCohort%c_area, &
+                                           currentCohort%n, currentCohort%canopy_layer,               &
+                                           currentPatch%canopy_layer_tlai )    
+
+          currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, &
+                                           currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
+                                           currentPatch%canopy_layer_tlai, currentCohort%treelai )  
 
           currentCohort%lai =  currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area 
           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 =  ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)  
 
           currentPatch%ncan(cl,ft) = max(currentPatch%ncan(cl,ft),currentCohort%NV)
 
           patch_lai = patch_lai + currentCohort%lai
 
-!          currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + &
-!                currentCohort%lai + currentCohort%sai
+          currentPatch%canopy_layer_tlai(cl) = currentPatch%canopy_layer_tlai(cl) + currentCohort%lai
 
-          do cl = 1,nclmax-1
-             if(currentCohort%canopy_layer == cl)then
-                currentPatch%canopy_layer_tai(cl) = currentPatch%canopy_layer_tai(cl) + &
-                     currentCohort%lai + currentCohort%sai
-             endif
-          enddo
-
-          currentCohort => currentCohort%taller 
+          currentCohort => currentCohort%shorter 
 
        enddo !currentCohort
 
@@ -1384,7 +1386,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
 
                    if(iv==currentCohort%NV) then
                       remainder = (currentCohort%treelai + currentCohort%treesai) - &
-                            (dinc_ed*dble(currentCohort%NV-1.0_r8))
+                            (dinc_ed*dble(currentCohort%nv-1.0_r8))
                       if(remainder > dinc_ed )then
                          write(fates_log(), *)'ED: issue with remainder', &
                                currentCohort%treelai,currentCohort%treesai,dinc_ed, & 
@@ -1484,7 +1486,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
                        enddo !currentCohort
                        call endrun(msg=errMsg(sourcefile, __LINE__))
                     end if
-                end do
+                 end do
 
                 do ft = 1,numpft
                    do iv = 1,currentPatch%ncan(cl,ft)

biogeochem/EDCohortDynamicsMod.F90

@@ -21,7 +21,6 @@ module EDCohortDynamicsMod
   use EDTypesMod            , only : min_npm2, min_nppatch
   use EDTypesMod            , only : min_n_safemath
   use EDTypesMod            , only : nlevleaf
-  use EDTypesMod            , only : dinc_ed
   use FatesInterfaceMod      , only : hlm_use_planthydro
   use FatesPlantHydraulicsMod, only : FuseCohortHydraulics
   use FatesPlantHydraulicsMod, only : CopyCohortHydraulics
@@ -152,11 +151,16 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine
     ! Assign canopy extent and depth
     call carea_allom(new_cohort%dbh,new_cohort%n,spread,new_cohort%pft,new_cohort%c_area)
 
-    new_cohort%treelai = tree_lai(new_cohort%bl, new_cohort%status_coh, new_cohort%pft, &
-         new_cohort%c_area, new_cohort%n)
+    new_cohort%treelai = tree_lai(new_cohort%bl, new_cohort%pft, new_cohort%c_area,    &
+                                  new_cohort%n, new_cohort%canopy_layer,               &
+                                  patchptr%canopy_layer_tlai )    
+
+    new_cohort%treesai = tree_sai(new_cohort%pft, new_cohort%dbh, new_cohort%canopy_trim,   &
+                                  new_cohort%c_area, new_cohort%n, new_cohort%canopy_layer, &
+                                  patchptr%canopy_layer_tlai, new_cohort%treelai )  
+
     new_cohort%lai     = new_cohort%treelai * new_cohort%c_area/patchptr%area
-    new_cohort%treesai = tree_sai(new_cohort%dbh, new_cohort%pft, new_cohort%canopy_trim, &
-         new_cohort%c_area, new_cohort%n)
+
 
     ! Put cohort at the right place in the linked list
     storebigcohort   => patchptr%tallest

biogeochem/EDPatchDynamicsMod.F90

@@ -1267,7 +1267,7 @@ subroutine zero_patch(cp_p)
     currentPatch%age                        = nan                          
     currentPatch%age_class                  = 1
     currentPatch%area                       = nan                                           
-    currentPatch%canopy_layer_tai(:)        = nan               
+    currentPatch%canopy_layer_tlai(:)       = nan               
     currentPatch%total_canopy_area          = nan
 
     currentPatch%tlai_profile(:,:,:)        = nan 
@@ -1347,7 +1347,7 @@ subroutine zero_patch(cp_p)
     currentPatch%burnt_frac_litter(:)       = 0.0_r8 
     currentPatch%btran_ft(:)                = 0.0_r8
 
-    currentPatch%canopy_layer_tai(:)        = 0.0_r8
+    currentPatch%canopy_layer_tlai(:)       = 0.0_r8
 
     currentPatch%seeds_in(:)                = 0.0_r8
     currentPatch%seed_decay(:)              = 0.0_r8

biogeochem/EDPhysiologyMod.F90

@@ -20,6 +20,7 @@ module EDPhysiologyMod
   use EDCohortDynamicsMod , only : create_cohort, sort_cohorts
   use FatesAllometryMod   , only : tree_lai
   use FatesAllometryMod   , only : tree_sai
+  use FatesAllometryMod   , only : decay_coeff_kn
 
   use EDTypesMod          , only : numWaterMem
   use EDTypesMod          , only : dl_sf, dinc_ed
@@ -162,40 +163,59 @@ subroutine trim_canopy( currentSite )
     ! Canopy trimming / leaf optimisation. Removes leaves in negative annual carbon balance. 
     !
     ! !USES:
-    !
-    !
+
     ! !ARGUMENTS    
     type (ed_site_type),intent(inout), target :: currentSite
     !
     ! !LOCAL VARIABLES:
     type (ed_cohort_type) , pointer :: currentCohort
     type (ed_patch_type)  , pointer :: currentPatch
 
-    integer  :: z          ! leaf layer
-    integer  :: ipft       ! pft index
-    integer  :: trimmed    ! was this layer trimmed in this year? If not expand the canopy. 
-    real(r8) :: tar_bl     ! target leaf biomass       (leaves flushed, trimmed)
-    real(r8) :: tar_bfr    ! target fine-root biomass  (leaves flushed, trimmed)
-    real(r8) :: bfr_per_bleaf ! ratio of fine root per leaf biomass
+    integer  :: z               ! leaf layer
+    integer  :: ipft            ! pft index
+    logical  :: trimmed         ! was this layer trimmed in this year? If not expand the canopy. 
+    real(r8) :: tar_bl          ! target leaf biomass       (leaves flushed, trimmed)
+    real(r8) :: tar_bfr         ! target fine-root biomass  (leaves flushed, trimmed)
+    real(r8) :: bfr_per_bleaf   ! ratio of fine root per leaf biomass
+    real(r8) :: sla_levleaf     ! sla at leaf level z
+    real(r8) :: nscaler_levleaf ! nscaler value at leaf level z
+    integer  :: cl              ! canopy layer index
+    real(r8) :: kn              ! nitrogen decay coefficient
+    real(r8) :: sla_max         ! Observational constraint on how large sla (m2/gC) can become
+
+    real(r8) :: leaf_inc           ! LAI-only portion of the vegetation increment of dinc_ed
+    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
+    real(r8) :: lai_current        ! the LAI in the current leaf layer
+    real(r8) :: cumulative_lai     ! the cumulative LAI, top down, to the leaf layer of interest
 
     !----------------------------------------------------------------------
 
     currentPatch => currentSite%youngest_patch
-
     do while(associated(currentPatch))
+       
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort)) 
-          trimmed = 0    
+
+          trimmed = .false.
           ipft = currentCohort%pft
           call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area)
-          currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%status_coh, currentCohort%pft, &
-               currentCohort%c_area, currentCohort%n )    
-          currentCohort%treesai = tree_sai(currentCohort%dbh, currentCohort%pft, currentCohort%canopy_trim, &
-               currentCohort%c_area, currentCohort%n)    
-          currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)
+          currentCohort%treelai = tree_lai(currentCohort%bl, currentCohort%pft, currentCohort%c_area, &
+                                           currentCohort%n, currentCohort%canopy_layer,               &
+                                           currentPatch%canopy_layer_tlai )    
+
+          currentCohort%treesai = tree_sai(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_trim, &
+                                           currentCohort%c_area, currentCohort%n, currentCohort%canopy_layer, &
+                                           currentPatch%canopy_layer_tlai, currentCohort%treelai )  
+
+          currentCohort%nv      = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed)
+
           if (currentCohort%nv > nlevleaf)then
-             write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv,currentCohort%treelai,currentCohort%treesai, &
-                  currentCohort%c_area,currentCohort%n,currentCohort%bl
+             write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv, &
+                   currentCohort%treelai,currentCohort%treesai, &
+                   currentCohort%c_area,currentCohort%n,currentCohort%bl
+             call endrun(msg=errMsg(sourcefile, __LINE__))
           endif
 
           call bleaf(currentcohort%dbh,ipft,currentcohort%canopy_trim,tar_bl)
@@ -206,35 +226,68 @@ subroutine trim_canopy( currentSite )
              bfr_per_bleaf = tar_bfr/tar_bl
           endif
 
+          ! Identify current canopy layer (cl)
+          cl = currentCohort%canopy_layer
+          
+          ! PFT-level maximum SLA value, even if under a thick canopy (same units as slatop)
+          sla_max = EDPftvarcon_inst%slamax(ipft)
+
           !Leaf cost vs netuptake for each leaf layer. 
-          do z = 1,nlevleaf
-             if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer. 
+          do z = 1, currentCohort%nv
+
+             ! Calculate the cumulative total vegetation area index (no snow occlusion, stems and leaves)
+
+             leaf_inc    = dinc_ed * &
+                   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 * (z-1)
+             lai_current       = min(leaf_inc, currentCohort%treelai - lai_layers_above)
+             cumulative_lai    = lai_canopy_above + lai_layers_above + 0.5*lai_current
+             
+             if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer.
+             
+                   
+                ! Calculate sla_levleaf following the sla profile with overlying leaf area
+                ! Scale for leaf nitrogen profile
+                kn = decay_coeff_kn(ipft)
+                ! Nscaler value at leaf level z
+                nscaler_levleaf = exp(-kn * cumulative_lai)
+                ! Sla value at leaf level z after nitrogen profile scaling (m2/gC)
+                sla_levleaf = EDPftvarcon_inst%slatop(ipft)/nscaler_levleaf
+
+                if(sla_levleaf > sla_max)then
+                   sla_levleaf = sla_max
+                end if
+                   
                 !Leaf Cost kgC/m2/year-1
                 !decidous costs. 
                 if (EDPftvarcon_inst%season_decid(ipft) == 1.or. &
                      EDPftvarcon_inst%stress_decid(ipft) == 1)then 
 
-
-                   currentCohort%leaf_cost =  1._r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8)
+                   ! Leaf cost at leaf level z accounting for sla profile (kgC/m2)
+                   currentCohort%leaf_cost =  1._r8/(sla_levleaf*1000.0_r8)
 
                    if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
                       ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment
                       ! to the leaf increment; otherwise do not.
                       currentCohort%leaf_cost = currentCohort%leaf_cost + &
-                           1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
+                           1.0_r8/(sla_levleaf*1000.0_r8) * &
                            bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
                    endif
 
                    currentCohort%leaf_cost = currentCohort%leaf_cost * &
                          (EDPftvarcon_inst%grperc(ipft) + 1._r8)
                 else !evergreen costs
-                   currentCohort%leaf_cost = 1.0_r8/(EDPftvarcon_inst%slatop(ipft)* &
+                   ! Leaf cost at leaf level z accounting for sla profile
+                   currentCohort%leaf_cost = 1.0_r8/(sla_levleaf* &
                         EDPftvarcon_inst%leaf_long(ipft)*1000.0_r8) !convert from sla in m2g-1 to m2kg-1
                    if ( int(EDPftvarcon_inst%allom_fmode(ipft)) .eq. 1 ) then
                       ! if using trimmed leaf for fine root biomass allometry, add the cost of the root increment
                       ! to the leaf increment; otherwise do not.
                       currentCohort%leaf_cost = currentCohort%leaf_cost + &
-                           1.0_r8/(EDPftvarcon_inst%slatop(ipft)*1000.0_r8) * &
+                           1.0_r8/(sla_levleaf*1000.0_r8) * &
                            bfr_per_bleaf / EDPftvarcon_inst%root_long(ipft)
                    endif
                    currentCohort%leaf_cost = currentCohort%leaf_cost * &
@@ -256,15 +309,15 @@ subroutine trim_canopy( currentSite )
                             currentCohort%laimemory = currentCohort%laimemory * &
                                   (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) 
                          endif
-                         trimmed = 1
+                         trimmed = .true.
                       endif
                    endif
                 endif
              endif !leaf activity? 
           enddo !z
 
           currentCohort%year_net_uptake(:) = 999.0_r8
-          if (trimmed == 0.and.currentCohort%canopy_trim < 1.0_r8)then
+          if ( (.not.trimmed) .and.currentCohort%canopy_trim < 1.0_r8)then
              currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft)
           endif 
 
@@ -2416,4 +2469,9 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out)
 
     end subroutine flux_into_litter_pools
 
+    ! ===================================================================================
+
+
+
+
 end module EDPhysiologyMod