Skip to content

Commit

Permalink
Merge pull request #7 from rgknox/xuchongang/fates_hydro_improve_read…
Browse files Browse the repository at this point in the history 
…y_to_merge

conflict resolutions
  • Loading branch information
@xuchongang
xuchongang authored Jan 14, 2019
2 parents 8971025 + 4f0310d commit 911b539
Show file tree
Hide file tree
Showing 13 changed files with 552 additions and 201 deletions.
31 changes: 17 additions & 14 deletions biogeochem/EDCohortDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ module EDCohortDynamicsMod
use EDTypesMod , only : min_npm2, min_nppatch
use EDTypesMod , only : min_n_safemath
use EDTypesMod , only : nlevleaf
use EDTypesMod , only : ican_upper
use FatesInterfaceMod , only : hlm_use_planthydro
use FatesInterfaceMod , only : hlm_parteh_mode
use FatesPlantHydraulicsMod, only : FuseCohortHydraulics
Expand Down Expand Up @@ -633,21 +634,25 @@ subroutine terminate_cohorts( currentSite, currentPatch, level )
endif ! if (.not.currentCohort%isnew .and. level == 2) then

if (terminate == 1) then

! preserve a record of the to-be-terminated cohort for mortality accounting
if (currentCohort%canopy_layer .eq. 1) then
levcan = 1
else
levcan = 2
endif
levcan = currentCohort%canopy_layer

currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) = &
currentSite%terminated_nindivs(currentCohort%size_class,currentCohort%pft,levcan) + currentCohort%n
!
currentSite%termination_carbonflux(levcan) = currentSite%termination_carbonflux(levcan) + &
currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)

if( hlm_use_planthydro == itrue ) then
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

currentSite%term_carbonflux_canopy = currentSite%term_carbonflux_canopy + &
currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
else
currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) = &
currentSite%term_nindivs_ustory(currentCohort%size_class,currentCohort%pft) + currentCohort%n

currentSite%term_carbonflux_ustory = currentSite%term_carbonflux_ustory + &
currentCohort%n * (struct_c+sapw_c+leaf_c+fnrt_c+store_c+repro_c)
end if

!put the litter from the terminated cohorts straight into the fragmenting pools
Expand Down Expand Up @@ -965,7 +970,6 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in)
currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn
currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn
currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn
currentCohort%fmort = (currentCohort%n*currentCohort%fmort + nextc%n*nextc%fmort)/newn
currentCohort%frmort = (currentCohort%n*currentCohort%frmort + nextc%n*nextc%frmort)/newn

! logging mortality, Yi Xu
Expand Down Expand Up @@ -1352,7 +1356,6 @@ subroutine copy_cohort( currentCohort,copyc )
! Mortality diagnostics
n%cmort = o%cmort
n%bmort = o%bmort
n%fmort = o%fmort
n%hmort = o%hmort
n%frmort = o%frmort

Expand Down
81 changes: 58 additions & 23 deletions biogeochem/EDPatchDynamicsMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@ module EDPatchDynamicsMod
use EDTypesMod , only : dtype_ifall
use EDTypesMod , only : dtype_ilog
use EDTypesMod , only : dtype_ifire
use EDTypesMod , only : ican_upper
use FatesInterfaceMod , only : hlm_use_planthydro
use FatesInterfaceMod , only : hlm_numSWb
use FatesInterfaceMod , only : bc_in_type
Expand Down Expand Up @@ -139,7 +140,6 @@ subroutine disturbance_rates( site_in, bc_in)
currentCohort%bmort = bmort
currentCohort%hmort = hmort
currentCohort%frmort = frmort
currentCohort%fmort = 0.0_r8 ! Fire mortality is initialized as zero, but may be changed

call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, &
lmort_direct,lmort_collateral,lmort_infra )
Expand Down Expand Up @@ -216,7 +216,6 @@ subroutine disturbance_rates( site_in, bc_in)
currentCohort => currentPatch%shortest
do while(associated(currentCohort))
if(currentCohort%canopy_layer == 1)then
currentCohort%fmort = 0.0_r8
currentCohort%cmort = currentCohort%cmort*(1.0_r8 - fates_mortality_disturbance_fraction)
currentCohort%hmort = currentCohort%hmort*(1.0_r8 - fates_mortality_disturbance_fraction)
currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction)
Expand Down Expand Up @@ -267,7 +266,6 @@ subroutine disturbance_rates( site_in, bc_in)
currentCohort%lmort_direct = 0.0_r8
currentCohort%lmort_collateral = 0.0_r8
currentCohort%lmort_infra = 0.0_r8
currentCohort%fmort = 0.0_r8
end if
currentCohort => currentCohort%taller
enddo !currentCohort
Expand Down Expand Up @@ -323,12 +321,14 @@ 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
integer :: levcan ! canopy level
real(r8) :: leaf_c ! leaf carbon [kg]
real(r8) :: fnrt_c ! fineroot carbon [kg]
real(r8) :: sapw_c ! sapwood carbon [kg]
real(r8) :: store_c ! storage carbon [kg]
real(r8) :: struct_c ! structure carbon [kg]
real(r8) :: total_c ! total carbon of plant [kg]
real(r8) :: fnrt_c ! fineroot carbon [kg]
real(r8) :: sapw_c ! sapwood carbon [kg]
real(r8) :: store_c ! storage carbon [kg]
real(r8) :: struct_c ! structure carbon [kg]
real(r8) :: total_c ! total carbon of plant [kg]

!---------------------------------------------------------------------

storesmallcohort => null() ! storage of the smallest cohort for insertion routine
Expand All @@ -351,12 +351,17 @@ subroutine spawn_patches( currentSite, bc_in)
call endrun(msg=errMsg(sourcefile, __LINE__))
end if

site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate
! Only create new patches that have non-negligible amount of land
if((currentPatch%area*currentPatch%disturbance_rate) > nearzero ) then
site_areadis = site_areadis + currentPatch%area * currentPatch%disturbance_rate
end if

currentPatch => currentPatch%older

enddo ! end loop over patches. sum area disturbed for all patches.

if (site_areadis > 0.0_r8) then
if (site_areadis > nearzero) then

cwd_ag_local = 0.0_r8
cwd_bg_local = 0.0_r8
leaf_litter_local = 0.0_r8
Expand All @@ -378,6 +383,8 @@ subroutine spawn_patches( currentSite, bc_in)
! This is the amount of patch area that is disturbed, and donated by the donor
patch_site_areadis = currentPatch%area * currentPatch%disturbance_rate

if (patch_site_areadis > nearzero) then

call average_patch_properties(currentPatch, new_patch, patch_site_areadis)

if (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. &
Expand Down Expand Up @@ -442,7 +449,6 @@ subroutine spawn_patches( currentSite, bc_in)
nc%cmort = nan ! The mortality diagnostics are set to nan because the cohort should dissappear
nc%hmort = nan
nc%bmort = nan
nc%fmort = nan
nc%frmort = nan
nc%lmort_direct = nan
nc%lmort_collateral = nan
Expand Down Expand Up @@ -485,7 +491,6 @@ subroutine spawn_patches( currentSite, bc_in)
! number density in EDCLMLink, and the number density of this new patch is donated
! so with the number density must come the effective mortality rates.

nc%fmort = 0.0_r8 ! Should had also been zero in the donor
nc%cmort = currentCohort%cmort
nc%hmort = currentCohort%hmort
nc%bmort = currentCohort%bmort
Expand All @@ -510,7 +515,6 @@ subroutine spawn_patches( currentSite, bc_in)
! Those remaining in the existing
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)

nc%fmort = 0.0_r8
nc%cmort = currentCohort%cmort
nc%hmort = currentCohort%hmort
nc%bmort = currentCohort%bmort
Expand All @@ -533,11 +537,40 @@ subroutine spawn_patches( currentSite, bc_in)
! loss of individuals from source patch due to area shrinking
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)

levcan = currentCohort%canopy_layer

if(levcan==ican_upper) then

! before changing number densities, track total rate of trees that died
! due to fire, as well as from each fire mortality term
currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) = &
currentSite%fmort_rate_canopy(currentCohort%size_class, currentCohort%pft) + &
nc%n * currentCohort%fire_mort / hlm_freq_day

currentSite%fmort_carbonflux_canopy = currentSite%fmort_carbonflux_canopy + &
(nc%n * currentCohort%fire_mort) * &
total_c * g_per_kg * days_per_sec * ha_per_m2

else
currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) = &
currentSite%fmort_rate_ustory(currentCohort%size_class, currentCohort%pft) + &
nc%n * currentCohort%fire_mort / hlm_freq_day

currentSite%fmort_carbonflux_ustory = currentSite%fmort_carbonflux_ustory + &
(nc%n * currentCohort%fire_mort) * &
total_c * g_per_kg * days_per_sec * ha_per_m2
end if

currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) = &
currentSite%fmort_rate_cambial(currentCohort%size_class, currentCohort%pft) + &
nc%n * currentCohort%cambial_mort / hlm_freq_day
currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) = &
currentSite%fmort_rate_crown(currentCohort%size_class, currentCohort%pft) + &
nc%n * currentCohort%crownfire_mort / hlm_freq_day

! loss of individual from fire in new patch.
nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort)

nc%fmort = currentCohort%fire_mort/hlm_freq_day

nc%cmort = currentCohort%cmort
nc%hmort = currentCohort%hmort
nc%bmort = currentCohort%bmort
Expand All @@ -546,7 +579,8 @@ subroutine spawn_patches( currentSite, bc_in)
nc%lmort_direct = currentCohort%lmort_direct
nc%lmort_collateral = currentCohort%lmort_collateral
nc%lmort_infra = currentCohort%lmort_infra



! Logging is the dominant disturbance
elseif (currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifall) .and. &
currentPatch%disturbance_rates(dtype_ilog) > currentPatch%disturbance_rates(dtype_ifire)) then ! Logging
Expand All @@ -567,7 +601,6 @@ subroutine spawn_patches( currentSite, bc_in)
nc%hmort = nan
nc%bmort = nan
nc%frmort = nan
nc%fmort = nan
nc%lmort_direct = nan
nc%lmort_collateral = nan
nc%lmort_infra = nan
Expand Down Expand Up @@ -609,7 +642,6 @@ subroutine spawn_patches( currentSite, bc_in)
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)


nc%fmort = 0.0_r8
nc%cmort = currentCohort%cmort
nc%hmort = currentCohort%hmort
nc%bmort = currentCohort%bmort
Expand All @@ -629,7 +661,6 @@ subroutine spawn_patches( currentSite, bc_in)
currentCohort%n = currentCohort%n * (1._r8 - patch_site_areadis/currentPatch%area)

! No grass impact mortality imposed on the newly created patch
nc%fmort = 0.0_r8
nc%cmort = currentCohort%cmort
nc%hmort = currentCohort%hmort
nc%bmort = currentCohort%bmort
Expand Down Expand Up @@ -682,10 +713,6 @@ subroutine spawn_patches( currentSite, bc_in)
enddo ! currentCohort
call sort_cohorts(currentPatch)

!zero disturbance accumulators
currentPatch%disturbance_rate = 0._r8
currentPatch%disturbance_rates = 0._r8

!update area of donor patch
currentPatch%area = currentPatch%area - patch_site_areadis

Expand All @@ -698,6 +725,12 @@ subroutine spawn_patches( currentSite, bc_in)
call terminate_cohorts(currentSite, currentPatch, 2)
call sort_cohorts(currentPatch)

end if ! if (patch_site_areadis > nearzero) then

!zero disturbance rate trackers
currentPatch%disturbance_rate = 0._r8
currentPatch%disturbance_rates = 0._r8

currentPatch => currentPatch%younger

enddo ! currentPatch patch loop.
Expand All @@ -723,9 +756,11 @@ subroutine spawn_patches( currentSite, bc_in)

endif !end new_patch area


call check_patch_area(currentSite)
call set_patchno(currentSite)

return
end subroutine spawn_patches

! ============================================================================
Expand Down
2 changes: 1 addition & 1 deletion biogeophys/FatesPlantRespPhotosynthMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -938,7 +938,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl, & ! in

if ( parsun_lsl <= 0._r8 ) then ! night time

anet_av_out = 0._r8
anet_av_out = -lmr
psn_out = 0._r8
rstoma_out = min(rsmax0, 1._r8/bbb * cf)

Expand Down
34 changes: 12 additions & 22 deletions fire/SFMainMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -158,14 +158,15 @@ end subroutine fire_danger_index
subroutine charecteristics_of_fuel ( currentSite )
!*****************************************************************

use SFParamsMod, only : SF_val_alpha_FMC, SF_val_SAV, SF_val_FBD
use SFParamsMod, only : SF_val_drying_ratio, SF_val_SAV, SF_val_FBD

type(ed_site_type), intent(in), target :: currentSite

type(ed_patch_type), pointer :: currentPatch
type(ed_cohort_type), pointer :: currentCohort

real(r8) timeav_swc
real(r8) timeav_swc
real(r8) alpha_FMC(nfsc) ! Relative fuel moisture adjusted per drying ratio
real(r8) fuel_moisture(nfsc) ! Scaled moisture content of small litter fuels.
real(r8) MEF(nfsc) ! Moisture extinction factor of fuels integer n

Expand Down Expand Up @@ -235,13 +236,15 @@ subroutine charecteristics_of_fuel ( currentSite )
! Equation 6 in Thonicke et al. 2010. across leaves,twig, small branch, and large branch
! dead leaves and twigs included in 1hr pool per Thonicke (2010)
! Calculate fuel moisture for trunks to hold value for fuel consumption
fuel_moisture(dl_sf:tr_sf) = exp(-1.0_r8 * SF_val_alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI)
alpha_FMC(dl_sf:tr_sf) = SF_val_SAV(dl_sf:tr_sf)/SF_val_drying_ratio

fuel_moisture(dl_sf:tr_sf) = exp(-1.0_r8 * alpha_FMC(dl_sf:tr_sf) * currentSite%acc_NI)

if(write_SF == itrue)then
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'fm ',fuel_moisture
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'csa ',currentSite%acc_NI
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC
if ( hlm_masterproc == itrue ) write(fates_log(),*) 'sfv ',alpha_FMC
endif
! FIX(RF,032414): needs refactoring.
! average water content !is this the correct metric?
Expand Down Expand Up @@ -408,8 +411,8 @@ subroutine rate_of_spread ( currentSite )
use SFParamsMod, only : SF_val_miner_total, &
SF_val_part_dens, &
SF_val_miner_damp, &
SF_val_fuel_energy, &
SF_val_wind_max
SF_val_fuel_energy

use FatesInterfaceMod, only : hlm_current_day, hlm_current_month

type(ed_site_type), intent(in), target :: currentSite
Expand All @@ -426,7 +429,6 @@ subroutine rate_of_spread ( currentSite )
real(r8) beta_ratio ! ratio of beta/beta_op
real(r8) a_beta ! dummy variable for product of a* beta_ratio for react_v_opt equation
real(r8) a,b,c,e ! function of fuel sav
real(r8) wind_elev_fire !wind speed (m/min) at elevevation relevant for fire

logical,parameter :: debug_windspeed = .false. !for debugging

Expand Down Expand Up @@ -493,21 +495,9 @@ subroutine rate_of_spread ( currentSite )

! Equation A5 in Thonicke et al. 2010
! phi_wind (unitless)
! convert wind_elev_fire from m/min to ft/min for Rothermel ROS eqn
! wind max per Lasslop et al 2014 to linearly reduce ROS for high wind speeds
!OLD! phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e))
if (currentPatch%effect_wspeed .le. SF_val_wind_max) then
wind_elev_fire = currentPatch%effect_wspeed
phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e))
if (debug_windspeed) write(fates_log(),*) 'SF wind LESS max ', currentPatch%effect_wspeed
if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day
else
!max condition 225 ft/min (FIREMIP Rabin table A10 JSBACH-Spitfire) convert to 68.577 m/min
wind_elev_fire = max(0.0_r8,(68.577_r8-0.5_r8*currentPatch%effect_wspeed))
phi_wind = c * ((3.281_r8*wind_elev_fire)**b)*(beta_ratio**(-e))
if (debug_windspeed) write(fates_log(),*) 'SF wind GREATER max ', currentPatch%effect_wspeed
if (debug_windspeed) write(fates_log(),*) 'month and day', hlm_current_month, hlm_current_day
endif
! convert current_wspeed (wind at elev relevant to fire) from m/min to ft/min for Rothermel ROS eqn
phi_wind = c * ((3.281_r8*currentPatch%effect_wspeed)**b)*(beta_ratio**(-e))


! ---propagating flux----
! Equation A2 in Thonicke et al.2010 and Eq. 42 Rothermal 1972
Expand Down
Loading

0 comments on commit 911b539

Please sign in to comment.