Merge pull request #349 from ckoven/parprof

Some radiation profile updates: clumping, SAI changes, and diagnostics

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -10,8 +10,8 @@ module EDCanopyStructureMod
   use EDPftvarcon           , only : EDPftvarcon_inst
   use FatesAllometryMod     , only : carea_allom
   use EDCohortDynamicsMod   , only : copy_cohort, terminate_cohorts, fuse_cohorts
-  use EDCohortDynamicsMod   , only : tree_lai
-  use EDCohortDynamicsMod   , only : tree_sai
+  use FatesAllometryMod     , only : tree_lai
+  use FatesAllometryMod     , only : tree_sai
   use EDtypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type, ncwd
   use EDTypesMod            , only : nclmax
   use EDTypesMod            , only : nlevleaf
@@ -930,9 +930,10 @@ subroutine canopy_summarization( nsites, sites, bc_in )
              call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
                                        currentCohort%size_class,currentCohort%size_by_pft_class)
 
-             call carea_allom(currentCohort%dbh,currentCohort%n, &
-                              sites(s)%spread,currentCohort%pft,currentCohort%c_area)
-             currentCohort%treelai = tree_lai(currentCohort)
+             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 )
 
              canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
 
@@ -1097,8 +1098,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
           ft = currentCohort%pft
           cl = currentCohort%canopy_layer
 
-          currentCohort%treelai = tree_lai(currentCohort)    
-          currentCohort%treesai = tree_sai(currentCohort)
+          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%lai =  currentCohort%treelai *currentCohort%c_area/currentPatch%total_canopy_area 
           currentCohort%sai =  currentCohort%treesai *currentCohort%c_area/currentPatch%total_canopy_area  
 

biogeochem/EDCohortDynamicsMod.F90

@@ -36,6 +36,7 @@ module EDCohortDynamicsMod
   use FatesAllometryMod  , only : h_allom
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : StructureResetOfDH
+  use FatesAllometryMod  , only : tree_lai, tree_sai
   ! CIME globals
   use shr_log_mod           , only : errMsg => shr_log_errMsg
   !
@@ -51,8 +52,6 @@ module EDCohortDynamicsMod
   public :: sort_cohorts
   public :: copy_cohort
   public :: count_cohorts
-  public :: tree_lai
-  public :: tree_sai
 
   logical, parameter :: DEBUG  = .false. ! local debug flag
 
@@ -153,9 +152,11 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     ! Assign canopy extent and depth
     call carea_allom(new_cohort%dbh,new_cohort%n,new_cohort%siteptr%spread,new_cohort%pft,new_cohort%c_area)
 
-    new_cohort%treelai = tree_lai(new_cohort)
+    new_cohort%treelai = tree_lai(new_cohort%bl, new_cohort%status_coh, new_cohort%pft, &
+         new_cohort%c_area, new_cohort%n)
     new_cohort%lai     = new_cohort%treelai * new_cohort%c_area/patchptr%area
-    new_cohort%treesai = 0.0_r8 !FIX(RF,032414)   
+    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
@@ -1232,82 +1233,6 @@ function count_cohorts( currentPatch ) result ( backcount )
     endif
 
   end function count_cohorts
-
-  ! =====================================================================================
-
-  real(r8) function tree_lai( cohort_in )
-
-    ! ============================================================================
-    !  LAI of individual trees is a function of the total leaf area and the total canopy area.   
-    ! ============================================================================
-
-    type(ed_cohort_type), intent(inout) :: cohort_in       
-
-    real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground.
-    real(r8) :: slat               ! the sla of the top leaf layer. m2/kgC
-
-    if( cohort_in%bl  <  0._r8 .or. cohort_in%pft  ==  0 ) then
-       write(fates_log(),*) 'problem in treelai',cohort_in%bl,cohort_in%pft
-    endif
-
-    if( cohort_in%status_coh  ==  2 ) then ! are the leaves on? 
-       slat = 1000.0_r8 * EDPftvarcon_inst%slatop(cohort_in%pft) ! m2/g to m2/kg
-       leafc_per_unitarea = cohort_in%bl/(cohort_in%c_area/cohort_in%n) !KgC/m2
-       if(leafc_per_unitarea > 0.0_r8)then
-          tree_lai = leafc_per_unitarea * slat  !kg/m2 * m2/kg = unitless LAI 
-       else
-          tree_lai = 0.0_r8
-       endif
-    else
-       tree_lai = 0.0_r8
-    endif !status
-    cohort_in%treelai = tree_lai
-
-    ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it
-    ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a 
-    ! huge error 
-    if(cohort_in%treelai > nlevleaf*dinc_ed)then
-       write(fates_log(),*) 'too much lai' , cohort_in%treelai , cohort_in%pft , nlevleaf * dinc_ed
-    endif
-
-    return
-
-  end function tree_lai
-  
-  ! ============================================================================
-
-  real(r8) function tree_sai( cohort_in )
-
-    ! ============================================================================
-    !  SAI of individual trees is a function of the total dead biomass per unit canopy area.   
-    ! ============================================================================
-
-    type(ed_cohort_type), intent(inout) :: cohort_in       
-
-    real(r8) :: bdead_per_unitarea ! KgC of leaf per m2 area of ground.
-    real(r8) :: sai_scaler     
-
-    sai_scaler = EDPftvarcon_inst%allom_sai_scaler(cohort_in%pft) 
-
-    if( cohort_in%bdead  <  0._r8 .or. cohort_in%pft  ==  0 ) then
-       write(fates_log(),*) 'problem in treesai',cohort_in%bdead,cohort_in%pft
-    endif
-
-    bdead_per_unitarea = cohort_in%bdead/(cohort_in%c_area/cohort_in%n) !KgC/m2
-    tree_sai = bdead_per_unitarea * sai_scaler !kg/m2 * m2/kg = unitless LAI 
-   
-    cohort_in%treesai = tree_sai
-
-    ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it
-    ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a 
-    ! huge error 
-    if(cohort_in%treesai > nlevleaf*dinc_ed)then
-       write(fates_log(),*) 'too much sai' , cohort_in%treesai , cohort_in%pft , nlevleaf * dinc_ed
-    endif
-
-    return
-
-  end function tree_sai
 
   ! ============================================================================
 

biogeochem/EDPhysiologyMod.F90

@@ -17,8 +17,8 @@ module EDPhysiologyMod
   use FatesInterfaceMod, only    : bc_in_type
   use EDCohortDynamicsMod , only : zero_cohort
   use EDCohortDynamicsMod , only : create_cohort, sort_cohorts
-  use EDCohortDynamicsMod , only : tree_lai
-  use EDCohortDynamicsMod , only : tree_sai
+  use FatesAllometryMod   , only : tree_lai
+  use FatesAllometryMod   , only : tree_sai
 
   use EDTypesMod          , only : numWaterMem
   use EDTypesMod          , only : dl_sf, dinc_ed
@@ -185,7 +185,10 @@ subroutine trim_canopy( currentSite )
           trimmed = 0    
           ipft = currentCohort%pft
           call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area)
-          currentCohort%treelai = tree_lai(currentCohort)    
+          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)
           if (currentCohort%nv > nlevleaf)then
              write(fates_log(),*) 'nv > nlevleaf',currentCohort%nv,currentCohort%treelai,currentCohort%treesai, &

biogeochem/FatesAllometryMod.F90

@@ -86,9 +86,13 @@ module FatesAllometryMod
   use EDPFTvarcon      , only : EDPftvarcon_inst
   use FatesConstantsMod, only : r8 => fates_r8
   use FatesConstantsMod, only : i4 => fates_int
+  use FatesConstantsMod, only : g_per_kg 
+  use FatesConstantsMod, only : cm2_per_m2
+  use FatesConstantsMod, only : kg_per_Megag
   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
 
   implicit none
 
@@ -98,6 +102,8 @@ module FatesAllometryMod
   public :: bagw_allom    ! Generic AGWB (above grnd. woody bio) wrapper
   public :: blmax_allom   ! Generic maximum leaf biomass wrapper
   public :: bleaf         ! Generic actual leaf biomass wrapper
+  public :: tree_lai      ! Calculate tree-level LAI from actual leaf biomass
+  public :: tree_sai      ! Calculate tree-level SAI from target leaf biomass
   public :: bsap_allom    ! Generic sapwood wrapper
   public :: bbgw_allom    ! Generic coarse root wrapper
   public :: bfineroot     ! Generic actual fine root biomass wrapper
@@ -494,6 +500,87 @@ subroutine bleaf(d,ipft,canopy_trim,bl,dbldd)
     return
   end subroutine bleaf
 
+  ! =====================================================================================
+  
+  real(r8) function tree_lai( bl, status_coh, pft, c_area, n )
+
+    ! ============================================================================
+    !  LAI of individual trees is a function of the total leaf area and the total canopy area.   
+    ! ============================================================================
+
+    real(r8), intent(in) :: bl            ! plant leaf biomass [kg]     
+    integer, intent(in)  :: status_coh    ! growth status of plant  (2 = leaves on , 1 = leaves off)
+    integer, intent(in)  :: pft
+    real(r8), intent(in) :: c_area        ! areal extent of canopy (m2)
+    real(r8), intent(in) :: n             ! number of individuals in cohort per 'area' (10000m2 default)
+
+    real(r8) :: leafc_per_unitarea ! KgC of leaf per m2 area of ground.
+    real(r8) :: slat               ! the sla of the top leaf layer. m2/kgC
+
+    if( bl  <  0._r8 .or. pft  ==  0 ) then
+       write(fates_log(),*) 'problem in treelai',bl,pft
+    endif
+
+    slat = g_per_kg * EDPftvarcon_inst%slatop(pft) ! m2/g to m2/kg
+    leafc_per_unitarea = bl/(c_area/n) !KgC/m2
+    if(leafc_per_unitarea > 0.0_r8)then
+       tree_lai = leafc_per_unitarea * slat  !kg/m2 * m2/kg = unitless LAI 
+    else
+       tree_lai = 0.0_r8
+    endif
+
+
+    ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it
+    ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a 
+    ! huge error 
+    if(tree_lai > nlevleaf*dinc_ed)then
+       write(fates_log(),*) 'too much lai' , tree_lai , pft , nlevleaf * dinc_ed
+       write(fates_log(),*) 'Aborting'
+       call endrun(msg=errMsg(sourcefile, __LINE__))
+    endif
+
+    return
+
+  end function tree_lai
+
+  ! ============================================================================
+
+  real(r8) function tree_sai( dbh, pft, canopy_trim, c_area, n )
+
+    ! ============================================================================
+    !  SAI of individual trees is a function of the target leaf biomass
+    ! ============================================================================
+
+    real(r8),intent(in)  :: dbh
+    integer, intent(in)  :: pft
+    real(r8),intent(in)  :: canopy_trim
+    real(r8), intent(in) :: c_area        ! areal extent of canopy (m2)
+    real(r8), intent(in) :: n             ! number of individuals in cohort per 'area' (10000m2 default)
+
+    real(r8) :: leafc_per_unitarea ! KgC of target leaf per m2 area of ground.
+    real(r8) :: sai_scaler     
+    real(r8) :: b_leaf
+
+    sai_scaler = g_per_kg * EDPftvarcon_inst%allom_sai_scaler(pft)  ! m2/g to m2/kg
+
+    call bleaf(dbh,pft,canopy_trim,b_leaf)
+
+    leafc_per_unitarea = b_leaf/(c_area/n) !KgC/m2
+
+    tree_sai = leafc_per_unitarea * sai_scaler !kg/m2 * m2/kg = unitless SAI 
+
+    ! here, if the LAI exceeeds the maximum size of the possible array, then we have no way of accomodating it
+    ! at the moments nlevleaf default is 40, which is very large, so exceeding this would clearly illustrate a 
+    ! huge error 
+    if(tree_sai > nlevleaf*dinc_ed)then
+       write(fates_log(),*) 'too much sai' , tree_sai , pft , nlevleaf * dinc_ed
+       write(fates_log(),*) 'Aborting'
+       call endrun(msg=errMsg(sourcefile, __LINE__))
+    endif
+
+    return
+
+  end function tree_sai
 
   ! ============================================================================
   ! Generic sapwood biomass interface
@@ -804,9 +891,7 @@ end subroutine bbgw_const
 
  subroutine bsap_deprecated(d,h,dhdd,bleaf,dbleafdd,ipft,bsap,dbsapdd)
 
-    use FatesConstantsMod, only : g_per_kg
-    use FatesConstantsMod, only : cm2_per_m2
-    use FatesConstantsMod, only : kg_per_Megag
+
 
     ! -------------------------------------------------------------------------
     ! -------------------------------------------------------------------------
@@ -853,10 +938,6 @@ end subroutine bsap_deprecated
 
   subroutine bsap_dlinear(d,h,dhdd,bleaf,dbleafdd,ipft,bsap,dbsapdd)
 
-    use FatesConstantsMod, only : g_per_kg
-    use FatesConstantsMod, only : cm2_per_m2
-    use FatesConstantsMod, only : kg_per_Megag
-    
     ! -------------------------------------------------------------------------
     ! Calculate sapwood biomass based on leaf area to sapwood area
     ! proportionality.  In this function, the leaftosapwood area is a function
@@ -1685,7 +1766,6 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h )
      ! T
      ! ============================================================================
 
-
      use FatesConstantsMod     , only : calloc_abs_error
      ! Arguments