Merge pull request #279 from ckoven/sitelevelspread

Modified crown allometric spread logic

biogeochem/EDCanopyStructureMod.F90

@@ -802,59 +802,51 @@ subroutine canopy_spread( currentSite )
     !  Calculates the spatial spread of tree canopies based on canopy closure.                             
     !
     ! !USES:
-    use EDParamsMod , only : ED_val_maxspread, ED_val_minspread 
+    use EDTypesMod        , only : AREA
+    use EDParamsMod, only : ED_val_canopy_closure_thresh    
     !
     ! !ARGUMENTS    
     type (ed_site_type), intent(inout), target :: currentSite
     !
     ! !LOCAL VARIABLES:
     type (ed_cohort_type), pointer :: currentCohort
     type (ed_patch_type) , pointer :: currentPatch
-    real(r8) :: arealayer(nlevleaf) ! Amount of canopy in each layer. 
+    real(r8) :: sitelevel_canopyarea  ! Amount of canopy in top layer at the site level
     real(r8) :: inc                   ! Arbitrary daily incremental change in canopy area 
     integer  :: z
     !----------------------------------------------------------------------
 
-    inc = 0.005_r8
+    inc = 0.05_r8
 
     currentPatch => currentSite%oldest_patch
 
+    sitelevel_canopyarea = 0.0_r8   
     do while (associated(currentPatch))
 
-       !calculate canopy area in each canopy storey...
-       arealayer = 0.0_r8   
+       !calculate canopy area in each patch...
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort))
           currentCohort%c_area = c_area(currentCohort) 
-          if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
-             arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area
+          if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then
+             sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area
           endif
           currentCohort => currentCohort%shorter
        enddo
 
-       !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
-       do z = 1,nclmax  
-
-          if(arealayer(z)/currentPatch%area > 0.9_r8)then
-             currentPatch%spread(z) = currentPatch%spread(z) - inc
-          else 
-             currentPatch%spread(z) = currentPatch%spread(z) + inc 
-          endif
-          if(currentPatch%spread(z) >= ED_val_maxspread)then 
-             currentPatch%spread(z) = ED_val_maxspread
-          endif
-          if(currentPatch%spread(z) <=  ED_val_minspread)then
-             currentPatch%spread(z) = ED_val_minspread
-          endif
-        enddo !z
-        !write(fates_log(),*) 'spread',currentPatch%spread(1:2)
-        !currentPatch%spread(:) = ED_val_maxspread
-        !FIX(RF,033114) spread is off
-        !write(fates_log(),*) 'canopy_spread',currentPatch%area,currentPatch%spread(1:2)
-        currentPatch => currentPatch%younger
+       currentPatch => currentPatch%younger
 
     enddo !currentPatch
 
+    !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner
+    if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then
+       currentSite%spread = currentSite%spread - inc
+    else 
+       currentSite%spread = currentSite%spread + inc 
+    endif
+
+    ! put within bounds to make sure it stays between 0 and 1
+    currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8)
+
   end subroutine canopy_spread
 
 

biogeochem/EDGrowthFunctionsMod.F90

@@ -225,14 +225,13 @@ real(r8) function c_area( cohort_in )
     ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. 
     ! ============================================================================
 
-    use EDParamsMod              , only : ED_val_grass_spread
     use EDTypesMod               , only : nclmax
 
     type(ed_cohort_type), intent(in) :: cohort_in       
 
     real(r8) :: dbh ! Tree diameter at breat height. cm. 
     real(r8) :: crown_area_to_dbh_exponent
-    integer  :: can_layer_index
+    real(r8) :: spreadterm
 
     ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth,
     ! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero)
@@ -244,9 +243,8 @@ real(r8) function c_area( cohort_in )
        write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%max_dbh
        write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody
        write(fates_log(),*) 'z_area 4',cohort_in%n
-       write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread
+       write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread
        write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer
-       write(fates_log(),*) 'z_area 7',ED_val_grass_spread
     end if
 
     dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft))
@@ -260,15 +258,12 @@ real(r8) function c_area( cohort_in )
     ! So we allow layer index exceedence here and force it down to max.
     ! (rgk/cdk 05/2017)
     ! ----------------------------------------------------------------------------------
-
-    can_layer_index = min(cohort_in%canopy_layer,nclmax)
-
-    if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then 
-       c_area = 3.142_r8 * cohort_in%n * &
-            (cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent
-    else
-       c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent
-    end if
+
+    ! apply site-level spread elasticity to the cohort crown allometry term
+    spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + &
+         (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft)
+    !
+    c_area = cohort_in%n * spreadterm * dbh ** crown_area_to_dbh_exponent
 
   end function c_area
 

biogeochem/EDPatchDynamicsMod.F90

@@ -276,7 +276,7 @@ subroutine spawn_patches( currentSite, bc_in)
     !
     ! !USES:
 
-    use EDParamsMod         , only : ED_val_maxspread, ED_val_understorey_death
+    use EDParamsMod         , only : ED_val_understorey_death
     use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts 
 
     !
@@ -300,7 +300,6 @@ subroutine spawn_patches( currentSite, bc_in)
     real(r8) :: leaf_litter_local(maxpft)    ! initial value of leaf litter. KgC/m2
     real(r8) :: cwd_ag_local(ncwd)           ! initial value of above ground coarse woody debris. KgC/m2
     real(r8) :: cwd_bg_local(ncwd)           ! initial value of below ground coarse woody debris. KgC/m2
-    real(r8) :: spread_local(nclmax)         ! initial value of canopy spread parameter.no units 
     !---------------------------------------------------------------------
 
     storesmallcohort => null() ! storage of the smallest cohort for insertion routine
@@ -328,12 +327,11 @@ subroutine spawn_patches( currentSite, bc_in)
        cwd_bg_local = 0.0_r8
        leaf_litter_local = 0.0_r8
        root_litter_local = 0.0_r8
-       spread_local(1:nclmax) = ED_val_maxspread
        age = 0.0_r8
 
        allocate(new_patch)
        call create_patch(currentSite, new_patch, age, site_areadis, &
-            spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, &
+            cwd_ag_local, cwd_bg_local, leaf_litter_local, &
             root_litter_local)
 
        new_patch%tallest  => null()
@@ -748,8 +746,6 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis
 
     enddo
 
-    newPatch%spread = newPatch%spread + currentPatch%spread * patch_site_areadis/newPatch%area    
-
   end subroutine average_patch_properties
 
   ! ============================================================================
@@ -1079,7 +1075,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
   end subroutine mortality_litter_fluxes
 
   ! ============================================================================
-  subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_local,cwd_bg_local, &
+  subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_local, &
        leaf_litter_local,root_litter_local)
     !
     ! !DESCRIPTION:
@@ -1096,7 +1092,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
     real(r8), intent(in) :: cwd_bg_local(:)     ! initial value of below ground coarse woody debris. KgC/m2
     real(r8), intent(in) :: root_litter_local(:)! initial value of root litter. KgC/m2
     real(r8), intent(in) :: leaf_litter_local(:)! initial value of leaf litter. KgC/m2
-    real(r8), intent(in) :: spread_local(:)     ! initial value of canopy spread parameter.no units 
     !
     ! !LOCAL VARIABLES:
     !---------------------------------------------------------------------
@@ -1126,7 +1121,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
     new_patch%age                = age   
     new_patch%age_class          = 1
     new_patch%area               = areap 
-    new_patch%spread             = spread_local
     new_patch%cwd_ag             = cwd_ag_local
     new_patch%cwd_bg             = cwd_bg_local
     new_patch%leaf_litter        = leaf_litter_local
@@ -1221,7 +1215,6 @@ subroutine zero_patch(cp_p)
     currentPatch%nrad(:,:)                  = 999    ! number of exposed leaf layers for each canopy layer and pft
     currentPatch%ncan(:,:)                  = 999    ! number of total leaf layers for each canopy layer and pft
     currentPatch%lai                        = nan    ! leaf area index of patch
-    currentPatch%spread(:)                  = nan    ! dynamic ratio of dbh to canopy area.
     currentPatch%pft_agb_profile(:,:)       = nan    
 
     ! DISTURBANCE 

main/EDInitMod.F90

@@ -111,6 +111,9 @@ subroutine zero_site( site_in )
     ! Resources management (logging/harvesting, etc)
     site_in%resources_management%trunk_product_site  = 0.0_r8
 
+    ! canopy spread
+    site_in%spread = 0._r8
+
   end subroutine zero_site
 
   ! ============================================================================
@@ -187,7 +190,7 @@ subroutine set_site_properties( nsites, sites)
        sites(s)%frac_burnt = 0.0_r8
        sites(s)%old_stock  = 0.0_r8
 
-
+       sites(s)%spread     = 1.0_r8
     end do
 
     return
@@ -204,7 +207,6 @@ subroutine init_patches( nsites, sites, bc_in)
      !
 
 
-     use EDParamsMod            , only : ED_val_maxspread
      use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps 
      use FatesInventoryInitMod,   only : initialize_sites_by_inventory
 
@@ -218,7 +220,6 @@ subroutine init_patches( nsites, sites, bc_in)
      integer  :: s
      real(r8) :: cwd_ag_local(ncwd)
      real(r8) :: cwd_bg_local(ncwd)
-     real(r8) :: spread_local(nclmax)
      real(r8) :: leaf_litter_local(maxpft)
      real(r8) :: root_litter_local(maxpft)
      real(r8) :: age !notional age of this patch
@@ -237,7 +238,6 @@ subroutine init_patches( nsites, sites, bc_in)
      cwd_bg_local(:)      = 0.0_r8 !ED_val_init_litter
      leaf_litter_local(:) = 0.0_r8
      root_litter_local(:) = 0.0_r8
-     spread_local(:)      = ED_val_maxspread
      age                  = 0.0_r8
      ! ---------------------------------------------------------------------------------------------
 
@@ -276,7 +276,7 @@ subroutine init_patches( nsites, sites, bc_in)
 
            ! make new patch...
            call create_patch(sites(s), newp, age, AREA, &
-                 spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local,  &
+                 cwd_ag_local, cwd_bg_local, leaf_litter_local,  &
                  root_litter_local) 
 
            call init_cohorts(newp, bc_in(s))

main/EDParamsMod.F90

@@ -23,11 +23,8 @@ module EDParamsMod
    real(r8),protected :: ED_size_diagnostic_scale             ! Flag to switch between a linear and exponential
                                                               ! scale on the plant size axis in diagnostics (NOT USED YET)
    real(r8),protected :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
-   real(r8),protected :: ED_val_grass_spread
    real(r8),protected :: ED_val_comp_excln
    real(r8),protected :: ED_val_stress_mort
-   real(r8),protected :: ED_val_maxspread                  ! maximum ratio of dbh to canopy area (cm/m2)
-   real(r8),protected :: ED_val_minspread                  ! minimum ratio of dbh to canopy area (cm/m2)
    real(r8),protected :: ED_val_init_litter
    real(r8),protected :: ED_val_nignitions
    real(r8),protected :: ED_val_understorey_death
@@ -47,14 +44,12 @@ module EDParamsMod
    real(r8),protected :: ED_val_phen_coldtemp
    real(r8),protected :: ED_val_cohort_fusion_tol
    real(r8),protected :: ED_val_patch_fusion_tol
+   real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry
 
    character(len=param_string_length),parameter :: ED_name_size_diagnostic_scale = "fates_size_diagnostic_scale"
    character(len=param_string_length),parameter :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac"
-   character(len=param_string_length),parameter :: ED_name_grass_spread = "fates_grass_spread"
    character(len=param_string_length),parameter :: ED_name_comp_excln = "fates_comp_excln"
    character(len=param_string_length),parameter :: ED_name_stress_mort = "fates_stress_mort"
-   character(len=param_string_length),parameter :: ED_name_maxspread = "fates_maxspread"
-   character(len=param_string_length),parameter :: ED_name_minspread = "fates_minspread"
    character(len=param_string_length),parameter :: ED_name_init_litter = "fates_init_litter"
    character(len=param_string_length),parameter :: ED_name_nignitions = "fates_nignitions"
    character(len=param_string_length),parameter :: ED_name_understorey_death = "fates_understorey_death"
@@ -73,7 +68,8 @@ module EDParamsMod
    character(len=param_string_length),parameter :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim"   
    character(len=param_string_length),parameter :: ED_name_phen_coldtemp= "fates_phen_coldtemp"   
    character(len=param_string_length),parameter :: ED_name_cohort_fusion_tol= "fates_cohort_fusion_tol"   
-   character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"   
+   character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol"
+   character(len=param_string_length),parameter :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh"      
 
    ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE)
    ! ----------------------------------------------------------------------------------------------
@@ -121,11 +117,8 @@ subroutine FatesParamsInit()
 
     ED_size_diagnostic_scale              = nan
     fates_mortality_disturbance_fraction  = nan
-    ED_val_grass_spread                   = nan
     ED_val_comp_excln                     = nan
     ED_val_stress_mort                    = nan
-    ED_val_maxspread                      = nan
-    ED_val_minspread                      = nan
     ED_val_init_litter                    = nan
     ED_val_nignitions                     = nan
     ED_val_understorey_death              = nan
@@ -145,6 +138,7 @@ subroutine FatesParamsInit()
     ED_val_phen_coldtemp                  = nan
     ED_val_cohort_fusion_tol              = nan
     ED_val_patch_fusion_tol               = nan
+    ED_val_canopy_closure_thresh               = nan    
 
     hydr_psi0                             = nan
     hydr_psicap                           = nan
@@ -179,21 +173,12 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
-    call fates_params%RegisterParameter(name=ED_name_grass_spread, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names)
-
     call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
     call fates_params%RegisterParameter(name=ED_name_stress_mort, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
-    call fates_params%RegisterParameter(name=ED_name_maxspread, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names)
-
-    call fates_params%RegisterParameter(name=ED_name_minspread, dimension_shape=dimension_shape_1d, &
-         dimension_names=dim_names)
-
     call fates_params%RegisterParameter(name=ED_name_init_litter, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
@@ -251,6 +236,9 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=ED_name_patch_fusion_tol, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
+    call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names)
+
     call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
@@ -291,21 +279,12 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, &
           data=fates_mortality_disturbance_fraction)
 
-    call fates_params%RetreiveParameter(name=ED_name_grass_spread, &
-         data=ED_val_grass_spread)
-
     call fates_params%RetreiveParameter(name=ED_name_comp_excln, &
          data=ED_val_comp_excln)
 
     call fates_params%RetreiveParameter(name=ED_name_stress_mort, &
          data=ED_val_stress_mort)
 
-    call fates_params%RetreiveParameter(name=ED_name_maxspread, &
-         data=ED_val_maxspread)
-
-    call fates_params%RetreiveParameter(name=ED_name_minspread, &
-         data=ED_val_minspread)
-
     call fates_params%RetreiveParameter(name=ED_name_init_litter, &
          data=ED_val_init_litter)
 
@@ -363,6 +342,9 @@ subroutine FatesReceiveParams(fates_params)
     call fates_params%RetreiveParameter(name=ED_name_patch_fusion_tol, &
          data=ED_val_patch_fusion_tol)
 
+    call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, &
+         data=ED_val_canopy_closure_thresh)
+
     call fates_params%RetreiveParameter(name=hydr_name_psi0, &
           data=hydr_psi0)
 
@@ -400,11 +382,8 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),*) '-----------  FATES Scalar Parameters -----------------'
         write(fates_log(),fmt0) 'ED_size_diagnostic_scale = ',ED_size_diagnostic_scale
         write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
-        write(fates_log(),fmt0) 'ED_val_grass_spread = ',ED_val_grass_spread
-        write(fates_log(),fmt0) 'ED_val_comp_excln = ', ED_val_comp_excln
+        write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln
         write(fates_log(),fmt0) 'ED_val_stress_mort = ',ED_val_stress_mort
-        write(fates_log(),fmt0) 'ED_val_maxspread = ',ED_val_maxspread
-        write(fates_log(),fmt0) 'ED_val_minspread = ',ED_val_minspread
         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
@@ -424,6 +403,7 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'ED_val_phen_coldtemp = ',ED_val_phen_coldtemp
         write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_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) 'hydr_psi0 = ',hydr_psi0
         write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap
         write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin