Merge pull request #440 from xuchongang/xuchongang/fates_hydro_improve_ready_to_merge

Improve the plant hydro codes to accomodate dynamic forest

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -25,7 +25,7 @@ module EDCanopyStructureMod
   use FatesInterfaceMod     , only : hlm_days_per_year
   use FatesInterfaceMod     , only : hlm_use_planthydro
   use FatesInterfaceMod     , only : numpft
-  use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort
+  use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage
 
   use PRTGenericMod,          only : leaf_organ
   use PRTGenericMod,          only : all_carbon_elements
@@ -503,6 +503,7 @@ subroutine DemoteFromLayer(currentSite,currentPatch,i_lyr)
                   ! demoted to the understory
 
                   allocate(copyc)
+
                   call InitPRTCohort(copyc)
                   if( hlm_use_planthydro.eq.itrue ) then
                      call InitHydrCohort(currentSite,copyc)
@@ -850,6 +851,7 @@ subroutine PromoteIntoLayer(currentSite,currentPatch,i_lyr)
                   elseif ( cc_gain > nearzero .and. cc_gain < currentCohort%c_area) then
 
                      allocate(copyc)
+
                      call InitPRTCohort(copyc)
                      if( hlm_use_planthydro.eq.itrue ) then
                         call InitHydrCohort(CurrentSite,copyc)
@@ -1759,6 +1761,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
      ! If hydraulics is turned on, update the amount of water bound in vegetation
      if (hlm_use_planthydro.eq.itrue) then
+        call RecruitWaterStorage(nsites,sites,bc_out)
         call UpdateH2OVeg(nsites,sites,bc_out)
      end if
 

biogeochem/EDCohortDynamicsMod.F90

@@ -8,6 +8,7 @@ module EDCohortDynamicsMod
   use FatesGlobals          , only : fates_log
   use FatesInterfaceMod     , only : hlm_freq_day
   use FatesInterfaceMod     , only : bc_in_type
+  use FatesInterfaceMod     , only : hlm_use_planthydro
   use FatesConstantsMod     , only : r8 => fates_r8
   use FatesConstantsMod     , only : fates_unset_int
   use FatesConstantsMod     , only : itrue,ifalse
@@ -31,6 +32,9 @@ module EDCohortDynamicsMod
   use FatesPlantHydraulicsMod, only : InitHydrCohort
   use FatesPlantHydraulicsMod, only : DeallocateHydrCohort
   use FatesPlantHydraulicsMod, only : AccumulateMortalityWaterStorage
+  use FatesPlantHydraulicsMod, only : UpdateTreeHydrNodes
+  use FatesPlantHydraulicsMod, only : UpdateTreeHydrLenVolCond
+  use FatesPlantHydraulicsMod, only : SavePreviousCompartmentVolumes
   use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index
   use FatesAllometryMod  , only : bleaf
   use FatesAllometryMod  , only : bfineroot
@@ -98,6 +102,13 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine
     !
     ! !DESCRIPTION:
     ! create new cohort
+    ! There are 4 places this is called
+    ! 1) Initializing new cohorts at the beginning of a cold-start simulation
+    ! 2) Initializing new recruits during dynamics
+    ! 3) Initializing new cohorts at the beginning of a inventory read
+    ! 4) Initializing new cohorts during restart
+    !
+    ! It is assumed that in the first 3, this is called with a reasonable amount of starter information.
     !
     ! !USES:
     !
@@ -125,7 +136,8 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine
     ! !LOCAL VARIABLES:
     type(ed_cohort_type), pointer :: new_cohort         ! Pointer to New Cohort structure.
     type(ed_cohort_type), pointer :: storesmallcohort 
-    type(ed_cohort_type), pointer :: storebigcohort   
+    type(ed_cohort_type), pointer :: storebigcohort  
+    integer :: nlevsoi_hyd                      ! number of hydraulically active soil layers 
     integer :: tnull,snull                      ! are the tallest and shortest cohorts allocate
     !----------------------------------------------------------------------
 
@@ -242,12 +254,31 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, bleaf, bfine
     new_cohort%isnew = .true.
 
     if( hlm_use_planthydro.eq.itrue ) then
-       call InitHydrCohort(CurrentSite,new_cohort)
-       call updateSizeDepTreeHydProps(CurrentSite,new_cohort, bc_in) 
-       call initTreeHydStates(CurrentSite,new_cohort, bc_in)
+
+       nlevsoi_hyd = currentSite%si_hydr%nlevsoi_hyd
+
+       ! This allocates array spaces
+       call InitHydrCohort(currentSite,new_cohort)
+
+       ! This calculates node heights
+       call UpdateTreeHydrNodes(new_cohort%co_hydr,new_cohort%pft, &
+                                new_cohort%hite,nlevsoi_hyd,bc_in)
+
+       ! This calculates volumes, lengths and max conductances
+       call UpdateTreeHydrLenVolCond(new_cohort,nlevsoi_hyd,bc_in)
+       
+       ! Since this is a newly initialized plant, we set the previous compartment-size
+       ! equal to the ones we just calculated.
+       call SavePreviousCompartmentVolumes(new_cohort%co_hydr)
+       
+       ! This comes up with starter suctions and then water contents
+       ! based on the soil values
+       call initTreeHydStates(currentSite,new_cohort, bc_in)
+
        if(recruitstatus==1)then
-          new_cohort%co_hydr%is_newly_recuited = .true.
+          new_cohort%co_hydr%is_newly_recruited = .true.
        endif
+
     endif
 
     call insert_cohort(new_cohort, patchptr%tallest, patchptr%shortest, tnull, snull, &
@@ -607,6 +638,9 @@ subroutine terminate_cohorts( currentSite, currentPatch, level )
           ! preserve a record of the to-be-terminated cohort for mortality accounting
           levcan = currentCohort%canopy_layer
 
+          if( hlm_use_planthydro == itrue ) &
+             call AccumulateMortalityWaterStorage(currentSite,currentCohort,currentCohort%n)
+
           if(levcan==ican_upper) then
              currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) = &
                    currentSite%term_nindivs_canopy(currentCohort%size_class,currentCohort%pft) + currentCohort%n
@@ -734,6 +768,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
      real(r8) :: newn
      real(r8) :: diff
      real(r8) :: dynamic_fusion_tolerance
+     real(r8) :: leaf_c             ! leaf carbon [kg]
 
      integer  :: largersc, smallersc, sc_i        ! indices for tracking the growth flux caused by fusion
      real(r8) :: larger_n, smaller_n
@@ -850,8 +885,11 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
 
                                 call sizetype_class_index(currentCohort%dbh,currentCohort%pft, &
                                       currentCohort%size_class,currentCohort%size_by_pft_class)
+				      
 
-                                if(hlm_use_planthydro.eq.itrue) call FuseCohortHydraulics(currentSite,currentCohort,nextc,bc_in,newn)
+                                if(hlm_use_planthydro.eq.itrue) then			  					  				  
+				    call FuseCohortHydraulics(currentSite,currentCohort,nextc,bc_in,newn)				    
+				 endif
 
                                 ! recent canopy history
                                 currentCohort%canopy_layer_yesterday  = (currentCohort%n*currentCohort%canopy_layer_yesterday  + &
@@ -980,7 +1018,18 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
                                 endif
 
                                 ! At this point, nothing should be pointing to current Cohort
-                                if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(nextc)
+				! update hydraulics quantities that are functions of hite & biomasses
+				! deallocate the hydro structure of nextc
+                                if (hlm_use_planthydro.eq.itrue) then				    
+				    call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, &
+				          currentCohort%pft,currentCohort%c_area)
+                                    leaf_c   = currentCohort%prt%GetState(leaf_organ, all_carbon_elements)
+                                    currentCohort%treelai = tree_lai(leaf_c,             &
+                                       currentCohort%pft, currentCohort%c_area, currentCohort%n, &
+                                       currentCohort%canopy_layer, currentPatch%canopy_layer_tlai )			    
+				   call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in)  				   
+				   call DeallocateHydrCohort(nextc)
+				endif
 
                                 ! Deallocate the cohort's PRT structure
                                 call nextc%prt%DeallocatePRTVartypes()

biogeochem/EDMortalityFunctionsMod.F90

@@ -16,6 +16,7 @@ module EDMortalityFunctionsMod
    use FatesInterfaceMod     , only : bc_in_type
    use FatesInterfaceMod     , only : hlm_use_ed_prescribed_phys
    use FatesInterfaceMod     , only : hlm_freq_day
+   use FatesInterfaceMod     , only : hlm_use_planthydro
    use EDLoggingMortalityMod , only : LoggingMortality_frac
    use EDParamsMod           , only : fates_mortality_disturbance_fraction
    use FatesInterfaceMod     , only : bc_in_type
@@ -62,8 +63,14 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     real(r8) :: leaf_c_target      ! target leaf biomass kgC
     real(r8) :: store_c
     real(r8) :: hf_sm_threshold    ! hydraulic failure soil moisture threshold 
+    real(r8) :: hf_flc_threshold   ! hydraulic failure fractional loss of conductivity threshold
     real(r8) :: temp_dep_fraction  ! Temp. function (freezing mortality)
     real(r8) :: temp_in_C          ! Daily averaged temperature in Celcius
+    real(r8) :: min_fmc_ag         ! minimum fraction of maximum conductivity for aboveground
+    real(r8) :: min_fmc_tr         ! minimum fraction of maximum conductivity for transporting root
+    real(r8) :: min_fmc_ar         ! minimum fraction of maximum conductivity for absorbing root
+    real(r8) :: min_fmc            ! minimum fraction of maximum conductivity for whole plant
+    real(r8) :: flc                ! fractional loss of conductivity 
     real(r8), parameter :: frost_mort_buffer = 5.0_r8  ! 5deg buffer for freezing mortality
     logical, parameter :: test_zero_mortality = .false. ! Developer test which
                                                         ! may help to debug carbon imbalances
@@ -77,12 +84,28 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
 
     ! Proxy for hydraulic failure induced mortality. 
     hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft)
-
-    if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= hf_sm_threshold)then 
+    hf_flc_threshold = EDPftvarcon_inst%hf_flc_threshold(cohort_in%pft)
+    if(hlm_use_planthydro.eq.itrue)then
+     !note the flc is set as the fraction of max conductivity in hydro
+     min_fmc_ag = minval(cohort_in%co_hydr%flc_ag(:))
+     min_fmc_tr = minval(cohort_in%co_hydr%flc_troot(:))
+     min_fmc_ar = minval(cohort_in%co_hydr%flc_aroot(:))
+     min_fmc = min(min_fmc_ag, min_fmc_tr)
+     min_fmc = min(min_fmc, min_fmc_ar)
+     flc = 1.0_r8-min_fmc
+     if(flc >= hf_flc_threshold .and. hf_flc_threshold < 1.0_r8 )then 
+       hmort = (flc-hf_flc_threshold)/(1.0_r8-hf_flc_threshold) * &
+           EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft)
+     else
+       hmort = 0.0_r8
+     endif      
+    else
+     if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= hf_sm_threshold)then 
        hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft)
      else
        hmort = 0.0_r8
-     endif 
+     endif
+    endif 
 
     ! Carbon Starvation induced mortality.
     if ( cohort_in%dbh  >  0._r8 ) then

biogeochem/EDPhysiologyMod.F90

@@ -13,12 +13,13 @@ module EDPhysiologyMod
   use FatesInterfaceMod, only    : hlm_day_of_year
   use FatesInterfaceMod, only    : numpft
   use FatesInterfaceMod, only    : hlm_use_planthydro
+  use FatesInterfaceMod, only    : hlm_parteh_mode
   use FatesConstantsMod, only    : r8 => fates_r8
   use FatesConstantsMod, only    : nearzero
   use EDPftvarcon      , only    : EDPftvarcon_inst
   use FatesInterfaceMod, only    : bc_in_type
   use EDCohortDynamicsMod , only : zero_cohort
-  use EDCohortDynamicsMod , only : create_cohort, sort_cohorts
+  use EDCohortDynamicsMod , only : create_cohort, sort_cohorts,InitPRTCohort
   use FatesAllometryMod   , only : tree_lai
   use FatesAllometryMod   , only : tree_sai
   use FatesAllometryMod   , only : decay_coeff_kn
@@ -39,6 +40,12 @@ module EDPhysiologyMod
   use EDParamsMod           , only : fates_mortality_disturbance_fraction
 
   use FatesPlantHydraulicsMod  , only : AccumulateMortalityWaterStorage
+  use FatesPlantHydraulicsMod  , only : updateSizeDepTreeHydProps
+  use FatesPlantHydraulicsMod  , only : initTreeHydStates
+  use FatesPlantHydraulicsMod  , only : InitHydrCohort
+  use FatesPlantHydraulicsMod  , only : ConstrainRecruitNumber
+  use FatesPlantHydraulicsMod  , only : DeallocateHydrCohort
+  
   use FatesConstantsMod     , only : itrue,ifalse
   use FatesConstantsMod     , only : calloc_abs_error
 
@@ -54,7 +61,8 @@ module EDPhysiologyMod
   use FatesAllometryMod  , only : carea_allom
   use FatesAllometryMod  , only : CheckIntegratedAllometries
   use FatesAllometryMod  , only : StructureResetOfDH
-
+  
+  use PRTGenericMod, only : prt_carbon_allom_hyp
   use PRTGenericMod, only : leaf_organ
   use PRTGenericMod, only : all_carbon_elements
   use PRTGenericMod, only : carbon12_element
@@ -971,6 +979,10 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
 
     allocate(temp_cohort) ! create temporary cohort
     call zero_cohort(temp_cohort)
+    if( hlm_use_planthydro.eq.itrue ) then
+	call InitHydrCohort(CurrentSite,temp_cohort)
+    endif
+    call InitPRTCohort(temp_cohort)
 
     do ft = 1,numpft
 
@@ -1014,28 +1026,64 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
        else
           ! prescribed recruitment rates. number per sq. meter per year
           temp_cohort%n        = currentPatch%area * EDPftvarcon_inst%prescribed_recruitment(ft) * hlm_freq_day
-          ! modify the carbon balance accumulators to take into account the different way of defining recruitment
-          ! add prescribed rates as an input C flux, and the recruitment that would have otherwise occured as an output flux
-          ! (since the carbon associated with them effectively vanishes)
-          currentSite%flux_in = currentSite%flux_in + temp_cohort%n * &
-                (b_store + b_leaf + b_fineroot + b_sapwood + b_dead)
-          currentSite%flux_out = currentSite%flux_out + currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day
        endif
 
        if (temp_cohort%n > 0.0_r8 )then
-          if ( debug ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
-          call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
+          if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort '
+	  !constrain the number of individual based on rhyzosphere water availability
+	  if( hlm_use_planthydro.eq.itrue ) then
+	      call carea_allom(temp_cohort%dbh,temp_cohort%n,currentSite%spread, &
+				          ft,temp_cohort%c_area)
+	      if(associated(currentPatch%shortest)) then
+	         temp_cohort%canopy_layer =  currentPatch%shortest%canopy_layer
+	      else
+	         temp_cohort%canopy_layer = 1
+	      endif		
+	      temp_cohort%pft = ft
+	      select case(hlm_parteh_mode)
+                case (prt_carbon_allom_hyp)
+
+                  call SetState(temp_cohort%prt,leaf_organ, carbon12_element, b_leaf)
+                  call SetState(temp_cohort%prt,fnrt_organ, carbon12_element, b_fineroot)
+                  call SetState(temp_cohort%prt,sapw_organ, carbon12_element, b_sapwood)
+                  call SetState(temp_cohort%prt,store_organ, carbon12_element, b_store)
+                  call SetState(temp_cohort%prt,struct_organ, carbon12_element, b_dead)
+                  call SetState(temp_cohort%prt,repro_organ , carbon12_element, 0.0_r8)
+
+              end select
+              temp_cohort%treelai = tree_lai(b_leaf, ft,&
+				             temp_cohort%c_area,temp_cohort%n, &
+                                             temp_cohort%canopy_layer,currentPatch%canopy_layer_tlai)	      
+              call updateSizeDepTreeHydProps(CurrentSite,temp_cohort, bc_in) 
+              call initTreeHydStates(CurrentSite,temp_cohort, bc_in)
+ 	      call ConstrainRecruitNumber(currentSite,temp_cohort, bc_in)
+	  endif
+	  if(temp_cohort%n > 0.0_r8) then
+            call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, &
                 b_leaf, b_fineroot, b_sapwood, b_dead, b_store, &  
                 temp_cohort%laimemory, cohortstatus,recruitstatus, temp_cohort%canopy_trim, currentPatch%NCL_p, &
                 currentSite%spread, bc_in)
-
-
-          ! keep track of how many individuals were recruited for passing to history
-          currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
+            ! keep track of how many individuals were recruited for passing to history
+            currentSite%recruitment_rate(ft) = currentSite%recruitment_rate(ft) + temp_cohort%n
+	    ! modify the carbon balance accumulators to take into account the different way of defining recruitment
+            ! add prescribed rates as an input C flux, and the recruitment that would have otherwise occured as an output flux
+            ! (since the carbon associated with them effectively vanishes)
+	    ! check the water for hydraulics
+	    if (hlm_use_ed_prescribed_phys .ne. ifalse .and. EDPftvarcon_inst%prescribed_recruitment(ft) .ge. 0. ) then
+              currentSite%flux_in = currentSite%flux_in + temp_cohort%n * &
+                (b_store + b_leaf + b_fineroot + b_sapwood + b_dead)
+              currentSite%flux_out = currentSite%flux_out + currentPatch%area * currentPatch%seed_germination(ft)*hlm_freq_day
+	    endif 
+	    
+	  endif 
 
        endif
     enddo  !pft loop
-
+    !deallocate the temporatory cohort
+    if (hlm_use_planthydro.eq.itrue) call DeallocateHydrCohort(temp_cohort)
+    !Deallocate the cohort's PRT structure
+    call temp_cohort%prt%DeallocatePRTVartypes()
+    deallocate(temp_cohort%prt)
     deallocate(temp_cohort) ! delete temporary cohort
 
   end subroutine recruitment
@@ -1173,7 +1221,9 @@ subroutine CWD_Input( currentSite, currentPatch)
                 hlm_freq_day * currentPatch%area
 
           if( hlm_use_planthydro == itrue ) then
-             call AccumulateMortalityWaterStorage(currentSite,currentCohort,dead_n)
+             !call AccumulateMortalityWaterStorage(currentSite,currentCohort,dead_n)
+             call AccumulateMortalityWaterStorage(currentSite,currentCohort,&
+                                                  -1.0_r8 * currentCohort%dndt * hlm_freq_day)
           end if
 
 

biogeochem/FatesAllometryMod.F90

@@ -117,6 +117,7 @@ module FatesAllometryMod
   public :: decay_coeff_kn
   public :: StructureResetOfDH ! Method to set DBH to sync with structure biomass
   public :: CheckIntegratedAllometries
+  public :: CrownDepth
   public :: set_root_fraction  ! Generic wrapper to calculate normalized
                                ! root profiles
 
@@ -1843,6 +1844,33 @@ subroutine h2d_martcano(h,p1,p2,p3,d,dddh)
     return
   end subroutine h2d_martcano
 
+  ! =====================================================================================
+
+
+  subroutine CrownDepth(height,crown_depth)
+
+    ! -----------------------------------------------------------------------------------
+    ! This routine returns the depth of a plant's crown.  Which is the length
+    ! from the bottom of the crown to the top in the vertical dimension.
+    ! 
+    ! This code may be used as a wrapper if different hypotheses are wished to be
+    ! optioned.
+    ! -----------------------------------------------------------------------------------
+    
+    real(r8),intent(in)  :: height   ! The height of the plant   [m]
+    real(r8),intent(out) :: crown_depth ! The depth of the crown [m]
+    
+    ! Alternative Hypothesis:
+    ! crown depth from Poorter, Bongers & Bongers
+    ! crown_depth = exp(-1.169_r8)*cCohort%hite**1.098_r8   
+    
+    crown_depth               = min(height,0.1_r8)
+    
+    return
+ end subroutine CrownDepth
+
+
+
   ! =============================================================================
   ! Specific diameter to crown area allometries
   ! =============================================================================