From 271c185e39c612decf8afb1c066b581b591c8a32 Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Sat, 10 Nov 2018 18:52:37 -0600 Subject: [PATCH 01/11] Add notes from Marcos Nearly all of what's been added isn't actually implemented, but it helps when tracking comparisons --- ED/src/dynamics/disturbance.f90 | 9 +- ED/src/dynamics/forestry.f90 | 21 +++- ED/src/init/landuse_init.f90 | 187 +++++++++++++++++++++++++++----- ED/src/memory/disturb_coms.f90 | 18 +++ ED/src/memory/ed_max_dims.F90 | 5 +- ED/src/memory/ed_state_vars.F90 | 24 ++++ 6 files changed, 230 insertions(+), 34 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index bc3671c28..18969a5e6 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -46,6 +46,7 @@ subroutine apply_disturbances(cgrid) use ed_misc_coms , only : current_time & ! intent(in) , ibigleaf ! ! intent(in) use disturb_coms , only : min_patch_area & ! intent(in) + , treefall_hite_threshold & ! intent(in) ! Added from Marcos , mature_harvest_age & ! intent(in) , plantation_year & ! intent(in) , plantation_rotation & ! intent(in) @@ -138,9 +139,9 @@ subroutine apply_disturbances(cgrid) !------------------------------------------------------------------------------------! ! nnsp_ble stands for the new number of site patches for big leaf. Each site ! ! may have additional mypfts, except if land use is 1 (cropland/pasture) or 2 ! - ! (forest plantation), in which case only one PFT is allowed. ! + ! (forest plantation), or 8 (cropland), in which case only one PFT is allowed. ! !------------------------------------------------------------------------------------! - nnsp_ble = 2 + (n_dist_types - 2) * mypfts + nnsp_ble = 3 + (n_dist_types - 3) * mypfts !------------------------------------------------------------------------------------! @@ -1530,7 +1531,7 @@ subroutine site_disturbance_rates(year, cgrid) cpoly%secondary_harvest_target(isi) = 0.0 !------------------------------------------------------------------------! - else + else ! if landuse(12) > 0 !------------------------------------------------------------------------! ! Logging based on target biomass, leave the disturbance rates as ! @@ -1578,7 +1579,7 @@ subroutine site_disturbance_rates(year, cgrid) !------------------------------------------------------------------------! end if !------------------------------------------------------------------------------! - end if + end if ! End ianth+disturb=1 !------------------------------------------------------------------------------! end do siteloop diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 4bee2eff8..3747b3f8a 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -26,7 +26,10 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) use disturb_coms , only : ianth_disturb & ! intent(in) , lutime & ! intent(in) , min_patch_area & ! intent(in) - , min_harvest_biomass ! ! intent(in) + , min_harvest_biomass & ! intent(in) + , plantation_rotation & ! intent(in) + , mature_harvest_age ! ! intent(in) +! , min_oldgrowth ! ! intent(in) ! New; need to add use fuse_fiss_utils , only : terminate_patches ! ! subroutine use ed_max_dims , only : n_pft & ! intent(in) , n_dbh ! ! intent(in) @@ -49,6 +52,9 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) real :: secondary_harvest_target real :: site_harvest_target real :: site_harvestable_agb +! real :: site_hvmax_btimber +! real :: site_hvpot_btimber + real :: area_mature_primary real :: hvagb_mature_primary real :: area_mature_secondary @@ -59,6 +65,9 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) real :: lambda_mature_secondary real :: lambda_mature_plantation real :: harvest_deficit +! real, dimension(onsp) :: pat_hvmax_btimber +! real, dimension(onsp) :: pat_hvpot_btimber + !----- Local constants. ----------------------------------------------------------------! logical , parameter :: print_detail = .false. !---------------------------------------------------------------------------------------! @@ -73,6 +82,11 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) csite => cpoly%site(isi) !---------------------------------------------------------------------------------------! + !----- Initialise book keeping variables. -------------------------------------------! + pat_hvmax_btimber(:) = 0.0 + pat_hvpot_btimber(:) = 0.0 + !------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! @@ -94,10 +108,15 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) ! demand is for trees above a minimum size (which is very common practice in the ! ! tropics). ! !---------------------------------------------------------------------------------------! +! site_hvmax_btimber = 0. +! site_hvpot_btimber = 0. + site_harvestable_agb = 0. hpat_loop: do ipa=1,onsp cpatch => csite%patch(ipa) + ilu = csite%dist_type(ipa) + is_oldgrowth = csite%age(ipa) >= min_oldgrowth(ilu) !----- Select the minimum DBH depending on the forest category. ---------------------! select case(ilu) diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index b819065de..6c2f91772 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -19,6 +19,12 @@ subroutine landuse_init , num_lu_trans & ! intent(in) , ianth_disturb & ! intent(in) , lu_database ! ! intent(in) +! , sl_pft & ! intent(in) +! , sl_scale & ! intent(in) +! , sl_nyrs & ! intent(in) +! , sl_yr_first & ! intent(in) +! , sl_mindbh_harvest & ! intent(in) +! , sl_prob_harvest ! ! intent(in) use ed_misc_coms , only : iyeara & ! intent(in) , iyearz ! ! intent(in) use grid_coms , only : ngrids ! ! intent(in) @@ -26,6 +32,8 @@ subroutine landuse_init , huge_lu & ! intent(in) , n_pft & ! intent(in) , maxlist ! ! intent(in) + ! use detailed_coms , only : idetailed ! ! intent(in) ! Maybe not needed + implicit none !----- Local variables -----------------------------------------------------------------! @@ -42,6 +50,8 @@ subroutine landuse_init character(len=6) :: hform character(len=str_len) :: lu_name character(len=str_len) :: cdum +! character(len=13) :: hifmt ! if writing LU table +! character(len=15) :: hffmt ! if writing LU table integer :: nharvest integer , dimension(n_pft) :: harvest_pft integer :: nf @@ -63,6 +73,7 @@ subroutine landuse_init real , dimension(n_pft) :: harvprob_1ary real , dimension(n_pft) :: mindbh_2ary real , dimension(n_pft) :: harvprob_2ary + real , dimension(num_lu_trans) :: landuse_now real :: lu_area real :: lu_area_i real :: wlon @@ -76,14 +87,19 @@ subroutine landuse_init real , parameter :: huge_dbh = huge(1.) !----- External function. --------------------------------------------------------------! real , external :: dist_gc +! real , external :: solid_area ! Need only for writing LU table !---------------------------------------------------------------------------------------! !----- Finding number of simulation years ----------------------------------------------! sim_years = iyearz-iyeara+1 + !------ Decide whether to write lu settings. ----------------------------------------! + ! write_lu_settings = btest(idetailed,6) .and. ianth_disturb /= 0 + !------------------------------------------------------------------------------------! + !----- Crashing the run if the user set up a very long run... --------------------------! - if (ianth_disturb == 1 .and. sim_years > max_lu_years) then + if (ianth_disturb /= 0 .and. sim_years > max_lu_years) then write (unit=*,fmt='(a,1x,i5)') 'IYEARA (From namelist) :',iyeara write (unit=*,fmt='(a,1x,i5)') 'IYEARZ (From namelist) : ',iyearz write (unit=*,fmt='(a,1x,i5)') 'MAX_LU_YEARS (From disturb_coms.f90): ',max_lu_years @@ -112,6 +128,31 @@ subroutine landuse_init cpoly => cgrid%polygon(ipy) select case (ianth_disturb) + case (0) + !------------------------------------------------------------------------------! + ! Anthropogenic disturbance is not used this time, allocate only a single ! + ! landuse year. ! + !------------------------------------------------------------------------------! + allocate(cpoly%clutimes(1,cpoly%nsites)) + !------------------------------------------------------------------------------! + + !----- Set the parameters in a way that no logging/ploughing will happen. -----! + do isi = 1,cpoly%nsites + cpoly%num_landuse_years(isi) = 1 + cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh + cpoly%probharv_primary (1:n_pft,isi) = 0. + cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh + cpoly%probharv_secondary(1:n_pft,isi) = 0. + cpoly%clutimes(1,isi)%landuse_year = iyeara + cpoly%clutimes(1,isi)%landuse(1:num_lu_trans) = 0.0 + end do + !------------------------------------------------------------------------------! + + + !----- No plantations. --------------------------------------------------------! + cpoly%plantation(:) = 0 + !------------------------------------------------------------------------------! + case (1) !------------------------------------------------------------------------------! @@ -145,6 +186,11 @@ subroutine landuse_init harvprob_1ary(1:n_pft) = 0. mindbh_2ary(1:n_pft) = huge_dbh harvprob_2ary(1:n_pft) = 0. + + !----- Initialise plantation patches if plantation information is available. --! + cpoly%plantation(:) = 0 + call read_plantation_fractions(cpoly,cgrid%lon(ipy),cgrid%lat(ipy),igr) + !----- Define the format for the header. --------------------------------------! write(hform,fmt='(a,i3.3,a)') '(a',str_len,')' @@ -334,6 +380,11 @@ subroutine landuse_init !----- Determine the number of disturbance years. -----------------------! cpoly%num_landuse_years(isi) = cpoly%num_landuse_years(1) + !----- PFT-dependent harvest characteristics. ------------------------! + cpoly%mindbh_primary (:,isi) = cpoly%mindbh_primary (:,1) + cpoly%probharv_primary (:,isi) = cpoly%probharv_primary (:,1) + cpoly%mindbh_secondary (:,isi) = cpoly%mindbh_secondary (:,1) + cpoly%probharv_secondary(:,isi) = cpoly%probharv_secondary(:,1) !----- Disturbances. ----------------------------------------------------! do iyear = 1,cpoly%num_landuse_years(isi) @@ -385,33 +436,113 @@ subroutine landuse_init close(unit=12,status='keep') !------------------------------------------------------------------------------! - !----- Initialise plantation patches if plantation information is available. --! - cpoly%plantation(:) = 0 - call read_plantation_fractions(cpoly,cgrid%lon(ipy),cgrid%lat(ipy),igr) - case (0) - !------------------------------------------------------------------------------! - ! Anthropogenic disturbance is not used this time, allocate only a single ! - ! landuse year. ! - !------------------------------------------------------------------------------! - allocate(cpoly%clutimes(1,cpoly%nsites)) - !------------------------------------------------------------------------------! - - !----- Set the parameters in a way that no logging/ploughing will happen. -----! - do isi = 1,cpoly%nsites - cpoly%num_landuse_years(isi) = 1 - cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh - cpoly%probharv_primary (1:n_pft,isi) = 0. - cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh - cpoly%probharv_secondary(1:n_pft,isi) = 0. - cpoly%clutimes(1,isi)%landuse_year = iyeara - cpoly%clutimes(1,isi)%landuse(1:num_lu_trans) = 0.0 - end do - !------------------------------------------------------------------------------! - - - !----- No plantations. --------------------------------------------------------! - cpoly%plantation(:) = 0 - !------------------------------------------------------------------------------! +! case (2) +! !---------------------------------------------------------------------------! +! ! Make the land use data based on ED2IN. ! +! ! Work with the first site, then copy the data to the others. ! +! !---------------------------------------------------------------------------! +! isi = 1 +! csite => cpoly%site(isi) +! !---------------------------------------------------------------------------! +! +! +! !----- No plantations. -----------------------------------------------------! +! cpoly%plantation(:) = 0 +! !---------------------------------------------------------------------------! +! +! +! +! !----- Determine the number of disturbance years. --------------------------! +! cpoly%num_landuse_years(isi) = sim_years +! !---------------------------------------------------------------------------! +! +! +! +! !----- File exists, allocate the maximum number of years. ------------------! +! allocate(cpoly%clutimes(sim_years,cpoly%nsites)) +! !---------------------------------------------------------------------------! +! +! +! !----- Initialise the PFT-dependent arrays. --------------------------------! +! cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh +! cpoly%probharv_primary (1:n_pft,isi) = 0. +! cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh +! cpoly%probharv_secondary(1:n_pft,isi) = 0. +! !---------------------------------------------------------------------------! +! +! +! +! !------ Find the number of PFT that can be harvested. ----------------------! +! nharvest = count(sl_pft >= 1 .and. sl_pft <= n_pft) +! !---------------------------------------------------------------------------! +! +! !----- Fill the arrays with the appropriate PFT. ---------------------------! +! harvloop_two: do h=1,nharvest +! ipft = sl_pft(h) +! if (ipft >= 1 .and. ipft <= n_pft) then +! cpoly%mindbh_harvest(ipft,isi) = sl_mindbh_harvest(h) +! cpoly%prob_harvest (ipft,isi) = sl_prob_harvest (h) +! end if +! end do harvloop_two +! !---------------------------------------------------------------------------! +! +! +! +! !---------------------------------------------------------------------------! +! ! Fill in the disturbance matrices and biomass target. ! +! !---------------------------------------------------------------------------! +! iyear = 0 +! do yd_this = iyeara,iyearz +! iyear = iyear + 1 +! clutime => cpoly%clutimes(iyear,isi) +! +! clutime%landuse_year = yd_this +! clutime%landuse(1:num_lu_trans) = 0. +! +! !------------------------------------------------------------------------! +! ! Decide whether to include logging disturbance in this year. ! +! !------------------------------------------------------------------------! +! if (yd_this >= sl_yr_first) then +! if ( (sl_scale == 1) .or. (mod(yd_this-sl_yr_first,sl_nyrs) == 0) ) & +! then +! clutime%landuse(12) = -1.0 +! clutime%landuse(14) = -1.0 +! end if +! end if +! !------------------------------------------------------------------------! +! end do +! !---------------------------------------------------------------------------! +! +! +! !---------------------------------------------------------------------------! +! ! Copy the information from the first site to the others. ! +! !---------------------------------------------------------------------------! +! siteloop_two: do isi = 2,cpoly%nsites +! csite => cpoly%site(isi) +! +! !----- Determine the number of disturbance years. -----------------------! +! cpoly%num_landuse_years(isi) = cpoly%num_landuse_years(1) +! !------------------------------------------------------------------------! +! +! +! !----- PFT-dependent harvest characteristics. ---------------------------! +! cpoly%mindbh_harvest(:,isi) = cpoly%mindbh_harvest(:,1) +! cpoly%prob_harvest (:,isi) = cpoly%prob_harvest (:,1) +! !------------------------------------------------------------------------! +! +! +! +! !----- Disturbances. ----------------------------------------------------! +! do iyear = 1,cpoly%num_landuse_years(isi) +! clutime => cpoly%clutimes(iyear,isi) +! onelutime => cpoly%clutimes(iyear,1) +! clutime%landuse_year = onelutime%landuse_year +! clutime%landuse(1:num_lu_trans) = onelutime%landuse(1:num_lu_trans) +! end do +! !------------------------------------------------------------------------! +! +! end do siteloop_two +! end select end do polyloop diff --git a/ED/src/memory/disturb_coms.f90 b/ED/src/memory/disturb_coms.f90 index 7b6b7a3e0..4d25c1a65 100644 --- a/ED/src/memory/disturb_coms.f90 +++ b/ED/src/memory/disturb_coms.f90 @@ -59,6 +59,16 @@ module disturb_coms !----- Dimensionless parameter controlling speed of fire spread. -----------------------! real :: fire_parameter +! +! ! FROM MARCOS +! !----- Fractions of fast and structural carbon and nitrogen lost through combustion. ---! +! real :: f_combusted_fast_c +! real :: f_combusted_struct_c +! real :: f_combusted_fast_n +! real :: f_combusted_struct_n +! !---- Maximum height for non-grass cohort to be considered part of fuel. ---------------! +! real :: fuel_height_max + !---------------------------------------------------------------------------------------! ! Anthropogenic disturbance. 1 means that anthropogenic disturbances will be ! ! included, whereas 0 means that it won't. ! @@ -113,6 +123,14 @@ module disturb_coms real :: treefall_hite_threshold !----- Cut-off for different fire survivorship. ----------------------------------------! real :: fire_hite_threshold +! +! ! MARCOS -- Not implemented +! !---------------------------------------------------------------------------------------! +! ! Minimum age above which we disregard the disturbance type (land use) and assume ! +! ! old growth, thus allowing patch fusion to occur. ! +! !---------------------------------------------------------------------------------------! +! real, dimension(n_dist_types) :: min_oldgrowth +! !=======================================================================================! !=======================================================================================! diff --git a/ED/src/memory/ed_max_dims.F90 b/ED/src/memory/ed_max_dims.F90 index a63b7f739..bf82ee03f 100644 --- a/ED/src/memory/ed_max_dims.F90 +++ b/ED/src/memory/ed_max_dims.F90 @@ -179,8 +179,11 @@ module ed_max_dims ! 4 -- Fire. ! ! 5 -- Forest regrowth. ! ! 6 -- Logged forest. ! + ! MARCOS Added ! + ! 7 -- Understory Thin -- Logged forest (skid trail + road). ! + ! 8 -- Cropland. ! !---------------------------------------------------------------------------------------! - integer, parameter :: n_dist_types = 6 + integer, parameter :: n_dist_types = 8 ! Changed from 6 to include understory thin !---------------------------------------------------------------------------------------! diff --git a/ED/src/memory/ed_state_vars.F90 b/ED/src/memory/ed_state_vars.F90 index 7d1d6739a..05f71cf41 100644 --- a/ED/src/memory/ed_state_vars.F90 +++ b/ED/src/memory/ed_state_vars.F90 @@ -1030,6 +1030,10 @@ module ed_state_vars !!4. Fire.\n !!5. Forest regrowth.\n !!6. Logged forest.\n + ! ! MARCOS -- Not Implemented + ! !!7. Skid trails and roads.\n + ! !!8. Croplands.\n + real , pointer,dimension(:) :: fast_soil_C ! Date: Sat, 10 Nov 2018 20:28:41 -0600 Subject: [PATCH 02/11] Initial commits for harvest option based on area and mindbh --- ED/src/dynamics/disturbance.f90 | 63 +-- ED/src/dynamics/forestry.f90 | 652 +++++++++++++++++++++----------- ED/src/dynamics/mortality.f90 | 18 +- ED/src/init/landuse_init.f90 | 1 + ED/src/io/ed_opspec.F90 | 6 +- ED/src/memory/disturb_coms.f90 | 7 +- 6 files changed, 489 insertions(+), 258 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index bc3671c28..6a9f0d079 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -108,6 +108,7 @@ subroutine apply_disturbances(cgrid) real , dimension(n_pft,n_dbh) :: initial_agb real , dimension(n_pft,n_dbh) :: initial_basal_area real , dimension(n_pft) :: mindbh_harvest + real , dimension(n_pft) :: harvprob real :: pot_area_remain real :: area_fac real :: orig_area @@ -657,6 +658,7 @@ subroutine apply_disturbances(cgrid) disturbed = act_area_loss(ipa,new_lu) > tiny_num dist_path = 0 mindbh_harvest(1:n_pft) = huge(1.) + harvprob (1:n_pft) = 0. !---------------------------------------------------------------------! @@ -697,8 +699,10 @@ subroutine apply_disturbances(cgrid) !----- Logging. ---------------------------------------------------! if (mature_primary) then mindbh_harvest(1:n_pft) = cpoly%mindbh_primary(1:n_pft,isi) + harvprob (1:n_pft) = cpoly%harvprob_primary(1:n_pft,isi) else if (mature_secondary) then mindbh_harvest(1:n_pft) = cpoly%mindbh_secondary(1:n_pft,isi) + harvprob (1:n_pft) = cpoly%harvprob_secondary(1:n_pft,isi) end if !------------------------------------------------------------------! end select @@ -725,7 +729,7 @@ subroutine apply_disturbances(cgrid) end if !----- Find the mortality associated with disturbance. ------------! call disturbance_mortality(csite,ipa,dist_rate,new_lu & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) !------------------------------------------------------------------! @@ -743,9 +747,9 @@ subroutine apply_disturbances(cgrid) area_fac = act_area_loss(ipa,new_lu) / csite%area(onsp+new_lu) call increment_patch_vars(csite,onsp+new_lu,ipa,area_fac) call insert_survivors(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) call accum_dist_litt(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) !---------------------------------------------------------------! case (1) !---------------------------------------------------------------! @@ -760,9 +764,9 @@ subroutine apply_disturbances(cgrid) area_fac = act_area_loss(ipa,new_lu)/csite%area(onsp+new_lu) call increment_patch_vars(csite,onsp+new_lu,ipa,area_fac) call insert_survivors(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) call accum_dist_litt(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) !------------------------------------------------------------! case default !------------------------------------------------------------! @@ -811,9 +815,9 @@ subroutine apply_disturbances(cgrid) if (same_pft) then call increment_patch_vars(csite,npa,ipa,area_fac) call insert_survivors(csite,npa,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) call accum_dist_litt(csite,npa,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob) end if !---------------------------------------------------------! end do @@ -1361,13 +1365,13 @@ subroutine site_disturbance_rates(year, cgrid) !------------------------------------------------------------------------------! - - !------------------------------------------------------------------------------! - ! Disturbance that creates new "burnt patches". Only non-cultivated ! - ! lands may suffer this disturbance. ! - !------------------------------------------------------------------------------! - cpoly%disturbance_rates(4,3:6,isi) = fire_disturbance_rate - !------------------------------------------------------------------------------! +! REDUNDANT WITH ABOVE +! !------------------------------------------------------------------------------! +! ! Disturbance that creates new "burnt patches". Only non-cultivated ! +! ! lands may suffer this disturbance. ! +! !------------------------------------------------------------------------------! +! cpoly%disturbance_rates(4,3:6,isi) = fire_disturbance_rate +! !------------------------------------------------------------------------------! @@ -1443,12 +1447,13 @@ subroutine site_disturbance_rates(year, cgrid) !------------------------------------------------------------------------------! ! Harvesting (either plantation -> plantation or logging) when a biomass ! - ! target does not exist (e.g. SimAmazonia)). Convert the ! + ! target does not exist (e.g. SimAmazonia)). Convert the ! ! harvest probability of being cut given that the DBH exceeds the minimum DBH. ! ! This is done only when anthropogenic disturbance is on and we are not seek- ! ! ing the biomass target, otherwise we set it to zero. ! !------------------------------------------------------------------------------! - if (ianth_disturb == 1) then + select case (ianth_disturb) + case (1) if (clutime%landuse(12) <= 0) then !----- Loop over all patches, and find the harvest probability. ---------! @@ -1578,7 +1583,19 @@ subroutine site_disturbance_rates(year, cgrid) !------------------------------------------------------------------------! end if !------------------------------------------------------------------------------! - end if + case (2) + !------------------------------------------------------------------------! + ! Logging based on grid fraction and tree size. Harvest creates new ! + ! patches and target biomass harvest set to zero. ! + !------------------------------------------------------------------------! + cpoly%disturbance_rates(6,6,isi) = clutime%landuse(13) + cpoly%disturbance_rates(2,2,isi) = clutime%landuse(15) + + cpoly%primary_harvest_target (isi) = 0.0 + cpoly%secondary_harvest_target (isi) = 0.0 + !------------------------------------------------------------------------! + + end select ! End ianth=1 !------------------------------------------------------------------------------! end do siteloop @@ -2700,7 +2717,7 @@ end subroutine increment_patch_vars ! This subroutine will populate the disturbed patch with the cohorts that were ! ! disturbed but did not go extinct. ! !---------------------------------------------------------------------------------------! - subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest) + subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob) use update_derived_props_module use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure @@ -2715,6 +2732,7 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest integer , intent(in) :: np integer , intent(in) :: cp real , dimension(n_pft), intent(in) :: mindbh_harvest + real , dimension(n_pft), intent(in) :: harvprob real , intent(in) :: area_fac !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch @@ -2744,7 +2762,7 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest mask(:) = .false. survivalloop: do ico = 1,cpatch%ncohorts - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) & + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2780,7 +2798,7 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest cohortloop: do ico = 1,cpatch%ncohorts - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) & + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2826,7 +2844,7 @@ end subroutine insert_survivors !=======================================================================================! ! This subroutine updates the litter pools after a disturbance takes place. ! !---------------------------------------------------------------------------------------! - subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest) + subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob) use ed_state_vars, only : sitetype & ! structure , patchtype & ! structure , polygontype ! ! structure @@ -2845,6 +2863,7 @@ subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest) integer , intent(in) :: np integer , intent(in) :: cp real , dimension(n_pft), intent(in) :: mindbh_harvest + real , dimension(n_pft), intent(in) :: harvprob integer , intent(in) :: new_lu real , intent(in) :: area_fac integer , intent(in) :: dist_path @@ -2892,7 +2911,7 @@ subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest) !---------------------------------------------------------------------------------! !----- Find survivorship. --------------------------------------------------------! - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) !---------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 4bee2eff8..4178e664e 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -63,232 +63,432 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) logical , parameter :: print_detail = .false. !---------------------------------------------------------------------------------------! - - !----- Nothing to do here in case anthropogenic disturbance is turned off. -------------! - if (ianth_disturb == 0) return - !---------------------------------------------------------------------------------------! - - - !----- Link to the current site. -------------------------------------------------------! - csite => cpoly%site(isi) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Set biomass targets based on current rates and unapplied harvest from previous ! - ! years (memory). These are in kgC/m2. In case DBH-based logging is applied, then ! - ! harvest target is set to zero. ! - !---------------------------------------------------------------------------------------! - primary_harvest_target = cpoly%primary_harvest_target (isi) & - + cpoly%primary_harvest_memory (isi) - secondary_harvest_target = cpoly%secondary_harvest_target(isi) & - + cpoly%primary_harvest_memory (isi) - site_harvest_target = primary_harvest_target + secondary_harvest_target - !---------------------------------------------------------------------------------------! - - - !---------------------------------------------------------------------------------------! - ! Find total harvestable biomass density in kgC/m2. The harvestable biomass is ! - ! not necessarily the same as the total site biomass, because it is possible that the ! - ! demand is for trees above a minimum size (which is very common practice in the ! - ! tropics). ! - !---------------------------------------------------------------------------------------! - site_harvestable_agb = 0. - hpat_loop: do ipa=1,onsp - cpatch => csite%patch(ipa) - ilu = csite%dist_type(ipa) - - !----- Select the minimum DBH depending on the forest category. ---------------------! - select case(ilu) - case (1) - !----- Agriculture. No timber harvesting here. ----------------------------------! - mindbh_harvest(:) = huge(1.) - !---------------------------------------------------------------------------------! - case (2) - !----- Forest plantation. Usually all biomass is cleared. -----------------------! - mindbh_harvest(:) = 0. - !---------------------------------------------------------------------------------! - case (4,5,6) - !----- Secondary vegetation. -----------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_secondary(:,isi) - !---------------------------------------------------------------------------------! - case (3) - !----- Primary vegetation. -------------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_primary (:,isi) - !---------------------------------------------------------------------------------! - end select - !------------------------------------------------------------------------------------! - - !------------------------------------------------------------------------------------! - ! Go over each cohort, seek harvestable biomass. ! - !------------------------------------------------------------------------------------! - hcoh_loop: do ico=1,cpatch%ncohorts - ipft = cpatch%pft(ico) - if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then - !----- Cohort is harvestable. -------------------------------------------------! - harvestable_agb(ipa) = harvestable_agb(ipa) & - + cpatch%nplant(ico) * cpatch%agb(ico) - !------------------------------------------------------------------------------! - end if - !---------------------------------------------------------------------------------! - end do hcoh_loop - !------------------------------------------------------------------------------------! - - !----- Update site biomass only when the patch has sufficient biomass. --------------! - if (harvestable_agb(ipa) >= min_harvest_biomass) then - site_harvestable_agb = site_harvestable_agb & - + harvestable_agb(ipa) * csite%area(ipa) - end if - !------------------------------------------------------------------------------------! - end do hpat_loop - !---------------------------------------------------------------------------------------! - - - - !=======================================================================================! - !=======================================================================================! - ! Find out whether any harvest can occur at this site. For harvest to occur, the ! - ! site must meet two criteria: ! - ! a. Harvestable biomass must be greater than min_harvest_biomass; ! - ! b. Target biomass must be greater than min_harvest_biomass. ! - ! ! - ! In case one or both criteria are not met, we add the harvest target to memory. ! - !---------------------------------------------------------------------------------------! - if ( site_harvestable_agb < min_harvest_biomass .or. & - site_harvest_target <= site_harvestable_agb * min_patch_area ) then - cpoly%primary_harvest_target (isi) = 0.0 - cpoly%secondary_harvest_target(isi) = 0.0 - cpoly%primary_harvest_memory (isi) = primary_harvest_target - cpoly%secondary_harvest_memory(isi) = secondary_harvest_target - !------------------------------------------------------------------------------------! - - - !------------------------------------------------------------------------------------! - ! Print the inventory and the target. ! - !------------------------------------------------------------------------------------! - if (print_detail) then - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a)' ) '---------------------------------------------' - write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST WON''T OCCUR THIS YEAR' - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi - write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & - , primary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & - , secondary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & - , site_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & - , site_harvestable_agb - write (unit=*,fmt='(a)' ) '---------------------------------------------' - write (unit=*,fmt='(a)' ) ' ' - end if - !------------------------------------------------------------------------------------! - return - end if - !---------------------------------------------------------------------------------------! - - - - !------ Compute current stocks of agb in mature forests. -------------------------------! - call inventory_mat_forests(cpoly,isi,onsp,harvestable_agb & - ,area_mature_primary ,hvagb_mature_primary & - ,area_mature_secondary ,hvagb_mature_secondary & - ,area_mature_plantation,hvagb_mature_plantation) - !---------------------------------------------------------------------------------------! - - - - !------ Compute the mature-forest harvest rates. ---------------------------------------! - call mat_forest_harv_rates(hvagb_mature_primary,hvagb_mature_secondary & - ,hvagb_mature_plantation,primary_harvest_target & - ,secondary_harvest_target,lambda_mature_primary & - ,lambda_mature_secondary,lambda_mature_plantation & - ,harvest_deficit) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Compute the area lost by mature patches through harvest. ! - !---------------------------------------------------------------------------------------! - call area_harvest_mature(cpoly,isi,onsp,harvestable_agb,pot_area_harv & - ,lambda_mature_primary,lambda_mature_secondary & - ,lambda_mature_plantation) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Compute the area lost by immature patches to meet the biomass demand. ! - !---------------------------------------------------------------------------------------! - call area_harvest_immature(cpoly,isi,onsp,harvestable_agb,pot_area_harv,harvest_deficit) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Print the inventory and the target. ! - !---------------------------------------------------------------------------------------! - if (print_detail) then - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST RATES' - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi - write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & - , primary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & - , secondary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & - , site_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & - , site_harvestable_agb - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PRIMARY) = ' & - , hvagb_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (SECONDARY) = ' & - , hvagb_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PLANTATION) = ' & - , hvagb_mature_plantation - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PRIMARY) = ' & - , area_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (SECONDARY) = ' & - , area_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PLANTATION) = ' & - , area_mature_plantation - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PRIMARY) = ' & - , lambda_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (SECONDARY) = ' & - , lambda_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PLANTATION) = ' & - , lambda_mature_plantation - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a,1x,es12.5)') ' HARVEST DEFICIT = ', harvest_deficit - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(5(a,1x))' ) ' IPA',' LU',' AGE',' AREA' & - ,' HV_AREA' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - do ipa=1,onsp - write (unit=*,fmt='(2(i5,1x),3(f12.7,1x))') ipa,csite%dist_type(ipa) & - ,csite%age(ipa),csite%area(ipa) & - ,pot_area_harv(ipa) - end do - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(a)' ) ' ' - end if - !---------------------------------------------------------------------------------------! - - - !---------------------------------------------------------------------------------------! - ! Reset the primary forest memory, and assign any remaining deficit to the ! - ! secondary forest memory. ! - !---------------------------------------------------------------------------------------! - cpoly%primary_harvest_memory (isi) = 0.0 - cpoly%secondary_harvest_memory(isi) = harvest_deficit - !---------------------------------------------------------------------------------------! - + select case (ianth_disturb) + !----- Nothing to do because not basing harvest off of biomass -------------------------! + case (0) ! Nothing to do because anthropogenic disturbance is turned off + return + case (2) ! Harvest is based on size + area + return + !---------------------------------------------------------------------------------------! + + case (1) + !----- Link to the current site. -------------------------------------------------------! + csite => cpoly%site(isi) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Set biomass targets based on current rates and unapplied harvest from previous ! + ! years (memory). These are in kgC/m2. In case DBH-based logging is applied, then ! + ! harvest target is set to zero. ! + !---------------------------------------------------------------------------------------! + primary_harvest_target = cpoly%primary_harvest_target (isi) & + + cpoly%primary_harvest_memory (isi) + secondary_harvest_target = cpoly%secondary_harvest_target(isi) & + + cpoly%primary_harvest_memory (isi) + site_harvest_target = primary_harvest_target + secondary_harvest_target + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Find total harvestable biomass density in kgC/m2. The harvestable biomass is ! + ! not necessarily the same as the total site biomass, because it is possible that the ! + ! demand is for trees above a minimum size (which is very common practice in the ! + ! tropics). ! + !---------------------------------------------------------------------------------------! + site_harvestable_agb = 0. + hpat_loop: do ipa=1,onsp + cpatch => csite%patch(ipa) + ilu = csite%dist_type(ipa) + + !----- Select the minimum DBH depending on the forest category. ---------------------! + select case(ilu) + case (1) + !----- Agriculture. No timber harvesting here. ----------------------------------! + mindbh_harvest(:) = huge(1.) + !---------------------------------------------------------------------------------! + case (2) + !----- Forest plantation. Usually all biomass is cleared. -----------------------! + mindbh_harvest(:) = 0. + !---------------------------------------------------------------------------------! + case (4,5,6) + !----- Secondary vegetation. -----------------------------------------------------! + mindbh_harvest(:) = cpoly%mindbh_secondary(:,isi) + !---------------------------------------------------------------------------------! + case (3) + !----- Primary vegetation. -------------------------------------------------------! + mindbh_harvest(:) = cpoly%mindbh_primary (:,isi) + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Go over each cohort, seek harvestable biomass. ! + !------------------------------------------------------------------------------------! + hcoh_loop: do ico=1,cpatch%ncohorts + ipft = cpatch%pft(ico) + if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then + !----- Cohort is harvestable. -------------------------------------------------! + harvestable_agb(ipa) = harvestable_agb(ipa) & + + cpatch%nplant(ico) * cpatch%agb(ico) + !------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------! + end do hcoh_loop + !------------------------------------------------------------------------------------! + + !----- Update site biomass only when the patch has sufficient biomass. --------------! + if (harvestable_agb(ipa) >= min_harvest_biomass) then + site_harvestable_agb = site_harvestable_agb & + + harvestable_agb(ipa) * csite%area(ipa) + end if + !------------------------------------------------------------------------------------! + end do hpat_loop + !---------------------------------------------------------------------------------------! + + + + !=======================================================================================! + !=======================================================================================! + ! Find out whether any harvest can occur at this site. For harvest to occur, the ! + ! site must meet two criteria: ! + ! a. Harvestable biomass must be greater than min_harvest_biomass; ! + ! b. Target biomass must be greater than min_harvest_biomass. ! + ! ! + ! In case one or both criteria are not met, we add the harvest target to memory. ! + !---------------------------------------------------------------------------------------! + if ( site_harvestable_agb < min_harvest_biomass .or. & + site_harvest_target <= site_harvestable_agb * min_patch_area ) then + cpoly%primary_harvest_target (isi) = 0.0 + cpoly%secondary_harvest_target(isi) = 0.0 + cpoly%primary_harvest_memory (isi) = primary_harvest_target + cpoly%secondary_harvest_memory(isi) = secondary_harvest_target + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Print the inventory and the target. ! + !------------------------------------------------------------------------------------! + if (print_detail) then + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '---------------------------------------------' + write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST WON''T OCCUR THIS YEAR' + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi + write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & + , primary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & + , secondary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & + , site_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & + , site_harvestable_agb + write (unit=*,fmt='(a)' ) '---------------------------------------------' + write (unit=*,fmt='(a)' ) ' ' + end if + !------------------------------------------------------------------------------------! + return + end if + !---------------------------------------------------------------------------------------! + + + + !------ Compute current stocks of agb in mature forests. -------------------------------! + call inventory_mat_forests(cpoly,isi,onsp,harvestable_agb & + ,area_mature_primary ,hvagb_mature_primary & + ,area_mature_secondary ,hvagb_mature_secondary & + ,area_mature_plantation,hvagb_mature_plantation) + !---------------------------------------------------------------------------------------! + + + + !------ Compute the mature-forest harvest rates. ---------------------------------------! + call mat_forest_harv_rates(hvagb_mature_primary,hvagb_mature_secondary & + ,hvagb_mature_plantation,primary_harvest_target & + ,secondary_harvest_target,lambda_mature_primary & + ,lambda_mature_secondary,lambda_mature_plantation & + ,harvest_deficit) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Compute the area lost by mature patches through harvest. ! + !---------------------------------------------------------------------------------------! + call area_harvest_mature(cpoly,isi,onsp,harvestable_agb,pot_area_harv & + ,lambda_mature_primary,lambda_mature_secondary & + ,lambda_mature_plantation) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Compute the area lost by immature patches to meet the biomass demand. ! + !---------------------------------------------------------------------------------------! + call area_harvest_immature(cpoly,isi,onsp,harvestable_agb,pot_area_harv,harvest_deficit) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Print the inventory and the target. ! + !---------------------------------------------------------------------------------------! + if (print_detail) then + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST RATES' + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi + write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & + , primary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & + , secondary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & + , site_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & + , site_harvestable_agb + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PRIMARY) = ' & + , hvagb_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (SECONDARY) = ' & + , hvagb_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PLANTATION) = ' & + , hvagb_mature_plantation + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PRIMARY) = ' & + , area_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (SECONDARY) = ' & + , area_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PLANTATION) = ' & + , area_mature_plantation + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PRIMARY) = ' & + , lambda_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (SECONDARY) = ' & + , lambda_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PLANTATION) = ' & + , lambda_mature_plantation + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,es12.5)') ' HARVEST DEFICIT = ', harvest_deficit + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(5(a,1x))' ) ' IPA',' LU',' AGE',' AREA' & + ,' HV_AREA' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + do ipa=1,onsp + write (unit=*,fmt='(2(i5,1x),3(f12.7,1x))') ipa,csite%dist_type(ipa) & + ,csite%age(ipa),csite%area(ipa) & + ,pot_area_harv(ipa) + end do + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' ' + end if + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Reset the primary forest memory, and assign any remaining deficit to the ! + ! secondary forest memory. ! + !---------------------------------------------------------------------------------------! + cpoly%primary_harvest_memory (isi) = 0.0 + cpoly%secondary_harvest_memory(isi) = harvest_deficit + !---------------------------------------------------------------------------------------! + case (2) + !----- Link to the current site. -------------------------------------------------------! + csite => cpoly%site(isi) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Set biomass targets based on current rates and unapplied harvest from previous ! + ! years (memory). These are in kgC/m2. In case DBH-based logging is applied, then ! + ! harvest target is set to zero. ! + !---------------------------------------------------------------------------------------! + primary_harvest_target = 0 + secondary_harvest_target = 0 + site_harvest_target = primary_harvest_target + secondary_harvest_target + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Find total harvestable biomass density in kgC/m2. The harvestable biomass is ! + ! not necessarily the same as the total site biomass, because it is possible that the ! + ! demand is for trees above a minimum size (which is very common practice in the ! + ! tropics). ! + !---------------------------------------------------------------------------------------! + site_harvestable_agb = 0. + hpat_loop: do ipa=1,onsp + cpatch => csite%patch(ipa) + ilu = csite%dist_type(ipa) + + !----- Select the minimum DBH depending on the forest category. ---------------------! + select case(ilu) + case (1) + !----- Agriculture. No timber harvesting here. ----------------------------------! + mindbh_harvest(:) = huge(1.) + !---------------------------------------------------------------------------------! + case (2) + !----- Forest plantation. Usually all biomass is cleared. -----------------------! + mindbh_harvest(:) = 0. + !---------------------------------------------------------------------------------! + case (4,5,6) + !----- Secondary vegetation. -----------------------------------------------------! + mindbh_harvest(:) = cpoly%mindbh_secondary(:,isi) + !---------------------------------------------------------------------------------! + case (3) + !----- Primary vegetation. -------------------------------------------------------! + mindbh_harvest(:) = cpoly%mindbh_primary (:,isi) + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! + + !------------------------------------------------------------------------------------! + ! Go over each cohort, seek harvestable biomass. ! + !------------------------------------------------------------------------------------! + hcoh_loop: do ico=1,cpatch%ncohorts + ipft = cpatch%pft(ico) + if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then + !----- Cohort is harvestable. -------------------------------------------------! + harvestable_agb(ipa) = harvestable_agb(ipa) & + + cpatch%nplant(ico) * cpatch%agb(ico) + !------------------------------------------------------------------------------! + end if + !---------------------------------------------------------------------------------! + end do hcoh_loop + !------------------------------------------------------------------------------------! + + !----- Update site biomass only when the patch has sufficient biomass. --------------! + if (harvestable_agb(ipa) >= min_harvest_biomass) then + site_harvestable_agb = site_harvestable_agb & + + harvestable_agb(ipa) * csite%area(ipa) + end if + !------------------------------------------------------------------------------------! + end do hpat_loop + !---------------------------------------------------------------------------------------! + + + + !=======================================================================================! + !=======================================================================================! + ! Find out whether any harvest can occur at this site. For harvest to occur, the ! + ! site must meet two criteria: ! + ! a. Harvestable biomass must be greater than min_harvest_biomass; ! + ! b. Target biomass must be greater than min_harvest_biomass. ! + ! ! + ! In case one or both criteria are not met, we add the harvest target to memory. ! + !---------------------------------------------------------------------------------------! + if ( site_harvestable_agb < min_harvest_biomass .or. & + site_harvest_target <= site_harvestable_agb * min_patch_area ) then + cpoly%primary_harvest_target (isi) = 0.0 + cpoly%secondary_harvest_target(isi) = 0.0 + cpoly%primary_harvest_memory (isi) = primary_harvest_target + cpoly%secondary_harvest_memory(isi) = secondary_harvest_target + !------------------------------------------------------------------------------------! + + return + end if + !---------------------------------------------------------------------------------------! + + + + !------ Compute current stocks of agb in mature forests. -------------------------------! + call inventory_mat_forests(cpoly,isi,onsp,harvestable_agb & + ,area_mature_primary ,hvagb_mature_primary & + ,area_mature_secondary ,hvagb_mature_secondary & + ,area_mature_plantation,hvagb_mature_plantation) + !---------------------------------------------------------------------------------------! + + + + !------ Compute the mature-forest harvest rates. ---------------------------------------! + call mat_forest_harv_rates(hvagb_mature_primary,hvagb_mature_secondary & + ,hvagb_mature_plantation,primary_harvest_target & + ,secondary_harvest_target,lambda_mature_primary & + ,lambda_mature_secondary,lambda_mature_plantation & + ,harvest_deficit) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Compute the area lost by mature patches through harvest. ! + !---------------------------------------------------------------------------------------! + call area_harvest_mature(cpoly,isi,onsp,harvestable_agb,pot_area_harv & + ,lambda_mature_primary,lambda_mature_secondary & + ,lambda_mature_plantation) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Compute the area lost by immature patches to meet the biomass demand. ! + !---------------------------------------------------------------------------------------! + call area_harvest_immature(cpoly,isi,onsp,harvestable_agb,pot_area_harv,harvest_deficit) + !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! Print the inventory and the target. ! + !---------------------------------------------------------------------------------------! + if (print_detail) then + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST RATES' + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi + write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & + , primary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & + , secondary_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & + , site_harvest_target + write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & + , site_harvestable_agb + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PRIMARY) = ' & + , hvagb_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (SECONDARY) = ' & + , hvagb_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PLANTATION) = ' & + , hvagb_mature_plantation + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PRIMARY) = ' & + , area_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (SECONDARY) = ' & + , area_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PLANTATION) = ' & + , area_mature_plantation + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PRIMARY) = ' & + , lambda_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (SECONDARY) = ' & + , lambda_mature_secondary + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PLANTATION) = ' & + , lambda_mature_plantation + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,es12.5)') ' HARVEST DEFICIT = ', harvest_deficit + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(5(a,1x))' ) ' IPA',' LU',' AGE',' AREA' & + ,' HV_AREA' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + do ipa=1,onsp + write (unit=*,fmt='(2(i5,1x),3(f12.7,1x))') ipa,csite%dist_type(ipa) & + ,csite%age(ipa),csite%area(ipa) & + ,pot_area_harv(ipa) + end do + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' ' + end if + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! + ! Reset the primary forest memory, and assign any remaining deficit to the ! + ! secondary forest memory. ! + !---------------------------------------------------------------------------------------! + cpoly%primary_harvest_memory (isi) = 0.0 + cpoly%secondary_harvest_memory(isi) = 0.0 + !---------------------------------------------------------------------------------------! + + end select return end subroutine find_harvest_area !==========================================================================================! diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index c75294002..b4e1e08d7 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -109,7 +109,7 @@ end subroutine mortality_rates ! disturbance. ! !---------------------------------------------------------------------------------------! subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path & - ,mindbh_harvest) + ,mindbh_harvest,harvprob) use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure use ed_max_dims , only : n_pft ! ! intent(in) @@ -121,6 +121,7 @@ subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path integer , intent(in) :: new_lu integer , intent(in) :: dist_path real , dimension(n_pft), intent(in) :: mindbh_harvest + real , dimension(n_pft), intent(in) :: harvprob !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch integer :: ico @@ -129,7 +130,7 @@ subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path cpatch => csite%patch(ipa) do ico=1,cpatch%ncohorts - f_survival = survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) + f_survival = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) cpatch%mort_rate(5,ico) = cpatch%mort_rate(5,ico) & - log( f_survival & + (1.0 - f_survival) * exp(- disturbance_rate) ) @@ -159,13 +160,15 @@ end subroutine disturbance_mortality ! -- mindbh_harvest: minimum DBH for harvesting (selective logging and forest ! ! plantantions). If the tree DBH is greater than mindbh_harvest, the tree may ! ! be harvested, otherwise it may be damaged by logging but not harvested. ! + ! -- harvprob: the probability of harvest for a trees meeting mindbh (ianth=2) ! ! -- cpatch: current patch. ! ! -- ico: index for current cohort. ! !---------------------------------------------------------------------------------------! - real function survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) + real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) use ed_state_vars, only : patchtype ! ! structure use disturb_coms , only : treefall_hite_threshold & ! intent(in) - , fire_hite_threshold ! ! intent(in) + , fire_hite_threshold & ! intent(in) + , ianth_disturb ! ! intent(in) use pft_coms , only : treefall_s_ltht & ! intent(in) , treefall_s_gtht & ! intent(in) , fire_s_ltht & ! intent(in) @@ -176,6 +179,7 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) !----- Arguments. -------------------------------------------------------------------! type(patchtype) , target :: cpatch real , dimension(n_pft), intent(in) :: mindbh_harvest + real , dimension(n_pft), intent(in) :: harvprob integer , intent(in) :: ico integer , intent(in) :: new_lu integer , intent(in) :: dist_path @@ -280,7 +284,11 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,cpatch,ico) ! applied. ! !---------------------------------------------------------------------------------! if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then - survivorship = 0.0 + if (ianth_disturb == 2) then + survivorship = 1 - harvprob(ipft) + else + survivorship = 0.0 + end if else survivorship = treefall_s_ltht(ipft) end if diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index b819065de..b6b075d19 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -311,6 +311,7 @@ subroutine landuse_init if ( clutime%landuse(14) > 0. ) then clutime%landuse(14) = lu_area_i * clutime%landuse(14) end if + ! Note: These two columns aren't really used right now clutime%landuse(16) = lu_area_i * clutime%landuse(16) clutime%landuse(18) = lu_area_i * clutime%landuse(18) end do diff --git a/ED/src/io/ed_opspec.F90 b/ED/src/io/ed_opspec.F90 index 4275a9b84..4f87639e9 100644 --- a/ED/src/io/ed_opspec.F90 +++ b/ED/src/io/ed_opspec.F90 @@ -2067,9 +2067,9 @@ subroutine ed_opspec_misc end if end select - if (ianth_disturb < 0 .or. ianth_disturb > 1) then + if (ianth_disturb < 0 .or. ianth_disturb > 2) then write (reason,fmt='(a,1x,i4,a)') & - 'Invalid IANTH_DISTURB, it must be between 0 and 1. Yours is set to',ianth_disturb,'...' + 'Invalid IANTH_DISTURB, it must be between 0 and 2. Yours is set to',ianth_disturb,'...' call opspec_fatal(reason,'opspec_misc') ifaterr = ifaterr +1 end if @@ -2128,7 +2128,7 @@ subroutine ed_opspec_misc end do pftloop !----- Checking whether the user choice for agriculture and plantation make sense. -----! - if (ianth_disturb == 1) then + if (ianth_disturb > 0) then !------ Checking the plantation PFT. It must be a tree PFT. ------------------------! select case (plantation_stock) case (2,3,4,6,7,8,9,10,11,17) diff --git a/ED/src/memory/disturb_coms.f90 b/ED/src/memory/disturb_coms.f90 index 7b6b7a3e0..f015e930e 100644 --- a/ED/src/memory/disturb_coms.f90 +++ b/ED/src/memory/disturb_coms.f90 @@ -60,8 +60,11 @@ module disturb_coms real :: fire_parameter !---------------------------------------------------------------------------------------! - ! Anthropogenic disturbance. 1 means that anthropogenic disturbances will be ! - ! included, whereas 0 means that it won't. ! + ! Anthropogenic disturbance. ! + ! 0. No anthropogenic disturbance ! + ! 1. Anthropogenic disturbances; harvest based on biomass ! + ! 2. Anthropogenic disturbances; forest harvest base don land area ! + ! 3. Marcos' single-site scheme (not implemented) ! !---------------------------------------------------------------------------------------! integer :: ianth_disturb !---------------------------------------------------------------------------------------! From 6b68f445772096a7529f0991ac29cbeb6147d4df Mon Sep 17 00:00:00 2001 From: crollinson Date: Sat, 10 Nov 2018 20:43:17 -0600 Subject: [PATCH 03/11] remove finding target biomass for ianth_disturb=2 --- ED/src/dynamics/forestry.f90 | 197 ----------------------------------- 1 file changed, 197 deletions(-) diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 2371d39e2..193f0141e 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -300,203 +300,6 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) cpoly%primary_harvest_memory (isi) = 0.0 cpoly%secondary_harvest_memory(isi) = harvest_deficit !---------------------------------------------------------------------------------------! - case (2) - !----- Link to the current site. -------------------------------------------------------! - csite => cpoly%site(isi) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Set biomass targets based on current rates and unapplied harvest from previous ! - ! years (memory). These are in kgC/m2. In case DBH-based logging is applied, then ! - ! harvest target is set to zero. ! - !---------------------------------------------------------------------------------------! - primary_harvest_target = 0 - secondary_harvest_target = 0 - site_harvest_target = primary_harvest_target + secondary_harvest_target - !---------------------------------------------------------------------------------------! - - - !---------------------------------------------------------------------------------------! - ! Find total harvestable biomass density in kgC/m2. The harvestable biomass is ! - ! not necessarily the same as the total site biomass, because it is possible that the ! - ! demand is for trees above a minimum size (which is very common practice in the ! - ! tropics). ! - !---------------------------------------------------------------------------------------! - site_harvestable_agb = 0. - hpat_loop: do ipa=1,onsp - cpatch => csite%patch(ipa) - ilu = csite%dist_type(ipa) - - !----- Select the minimum DBH depending on the forest category. ---------------------! - select case(ilu) - case (1) - !----- Agriculture. No timber harvesting here. ----------------------------------! - mindbh_harvest(:) = huge(1.) - !---------------------------------------------------------------------------------! - case (2) - !----- Forest plantation. Usually all biomass is cleared. -----------------------! - mindbh_harvest(:) = 0. - !---------------------------------------------------------------------------------! - case (4,5,6) - !----- Secondary vegetation. -----------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_secondary(:,isi) - !---------------------------------------------------------------------------------! - case (3) - !----- Primary vegetation. -------------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_primary (:,isi) - !---------------------------------------------------------------------------------! - end select - !------------------------------------------------------------------------------------! - - !------------------------------------------------------------------------------------! - ! Go over each cohort, seek harvestable biomass. ! - !------------------------------------------------------------------------------------! - hcoh_loop: do ico=1,cpatch%ncohorts - ipft = cpatch%pft(ico) - if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then - !----- Cohort is harvestable. -------------------------------------------------! - harvestable_agb(ipa) = harvestable_agb(ipa) & - + cpatch%nplant(ico) * cpatch%agb(ico) - !------------------------------------------------------------------------------! - end if - !---------------------------------------------------------------------------------! - end do hcoh_loop - !------------------------------------------------------------------------------------! - - !----- Update site biomass only when the patch has sufficient biomass. --------------! - if (harvestable_agb(ipa) >= min_harvest_biomass) then - site_harvestable_agb = site_harvestable_agb & - + harvestable_agb(ipa) * csite%area(ipa) - end if - !------------------------------------------------------------------------------------! - end do hpat_loop - !---------------------------------------------------------------------------------------! - - - - !=======================================================================================! - !=======================================================================================! - ! Find out whether any harvest can occur at this site. For harvest to occur, the ! - ! site must meet two criteria: ! - ! a. Harvestable biomass must be greater than min_harvest_biomass; ! - ! b. Target biomass must be greater than min_harvest_biomass. ! - ! ! - ! In case one or both criteria are not met, we add the harvest target to memory. ! - !---------------------------------------------------------------------------------------! - if ( site_harvestable_agb < min_harvest_biomass .or. & - site_harvest_target <= site_harvestable_agb * min_patch_area ) then - cpoly%primary_harvest_target (isi) = 0.0 - cpoly%secondary_harvest_target(isi) = 0.0 - cpoly%primary_harvest_memory (isi) = primary_harvest_target - cpoly%secondary_harvest_memory(isi) = secondary_harvest_target - !------------------------------------------------------------------------------------! - - return - end if - !---------------------------------------------------------------------------------------! - - - - !------ Compute current stocks of agb in mature forests. -------------------------------! - call inventory_mat_forests(cpoly,isi,onsp,harvestable_agb & - ,area_mature_primary ,hvagb_mature_primary & - ,area_mature_secondary ,hvagb_mature_secondary & - ,area_mature_plantation,hvagb_mature_plantation) - !---------------------------------------------------------------------------------------! - - - - !------ Compute the mature-forest harvest rates. ---------------------------------------! - call mat_forest_harv_rates(hvagb_mature_primary,hvagb_mature_secondary & - ,hvagb_mature_plantation,primary_harvest_target & - ,secondary_harvest_target,lambda_mature_primary & - ,lambda_mature_secondary,lambda_mature_plantation & - ,harvest_deficit) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Compute the area lost by mature patches through harvest. ! - !---------------------------------------------------------------------------------------! - call area_harvest_mature(cpoly,isi,onsp,harvestable_agb,pot_area_harv & - ,lambda_mature_primary,lambda_mature_secondary & - ,lambda_mature_plantation) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Compute the area lost by immature patches to meet the biomass demand. ! - !---------------------------------------------------------------------------------------! - call area_harvest_immature(cpoly,isi,onsp,harvestable_agb,pot_area_harv,harvest_deficit) - !---------------------------------------------------------------------------------------! - - - - !---------------------------------------------------------------------------------------! - ! Print the inventory and the target. ! - !---------------------------------------------------------------------------------------! - if (print_detail) then - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST RATES' - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi - write (unit=*,fmt='(a,1x,es12.5)') ' PRIMARY TARGET AGB = ' & - , primary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' SECONDARY TARGET AGB = ' & - , secondary_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL TARGET AGB = ' & - , site_harvest_target - write (unit=*,fmt='(a,1x,es12.5)') ' TOTAL HARVESTABLE AGB = ' & - , site_harvestable_agb - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PRIMARY) = ' & - , hvagb_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (SECONDARY) = ' & - , hvagb_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AGB (PLANTATION) = ' & - , hvagb_mature_plantation - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PRIMARY) = ' & - , area_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (SECONDARY) = ' & - , area_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV AREA (PLANTATION) = ' & - , area_mature_plantation - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PRIMARY) = ' & - , lambda_mature_primary - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (SECONDARY) = ' & - , lambda_mature_secondary - write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PLANTATION) = ' & - , lambda_mature_plantation - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a,1x,es12.5)') ' HARVEST DEFICIT = ', harvest_deficit - write (unit=*,fmt='(a)' ) ' ' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(5(a,1x))' ) ' IPA',' LU',' AGE',' AREA' & - ,' HV_AREA' - write (unit=*,fmt='(a)' ) '------------------------------------------------' - do ipa=1,onsp - write (unit=*,fmt='(2(i5,1x),3(f12.7,1x))') ipa,csite%dist_type(ipa) & - ,csite%age(ipa),csite%area(ipa) & - ,pot_area_harv(ipa) - end do - write (unit=*,fmt='(a)' ) '------------------------------------------------' - write (unit=*,fmt='(a)' ) ' ' - end if - !---------------------------------------------------------------------------------------! - - - !---------------------------------------------------------------------------------------! - ! Reset the primary forest memory, and assign any remaining deficit to the ! - ! secondary forest memory. ! - !---------------------------------------------------------------------------------------! - cpoly%primary_harvest_memory (isi) = 0.0 - cpoly%secondary_harvest_memory(isi) = 0.0 - !---------------------------------------------------------------------------------------! - end select return end subroutine find_harvest_area From 1f0c0c8305db8a397a883bf22b04a5288b2ec153 Mon Sep 17 00:00:00 2001 From: crollinson Date: Sat, 10 Nov 2018 20:50:09 -0600 Subject: [PATCH 04/11] typo correct in harvest probability assignment --- ED/src/dynamics/disturbance.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index ec5fbd689..77e6ea94c 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -700,10 +700,10 @@ subroutine apply_disturbances(cgrid) !----- Logging. ---------------------------------------------------! if (mature_primary) then mindbh_harvest(1:n_pft) = cpoly%mindbh_primary(1:n_pft,isi) - harvprob (1:n_pft) = cpoly%harvprob_primary(1:n_pft,isi) + harvprob (1:n_pft) = cpoly%probharv_primary(1:n_pft,isi) else if (mature_secondary) then mindbh_harvest(1:n_pft) = cpoly%mindbh_secondary(1:n_pft,isi) - harvprob (1:n_pft) = cpoly%harvprob_secondary(1:n_pft,isi) + harvprob (1:n_pft) = cpoly%probharv_secondary(1:n_pft,isi) end if !------------------------------------------------------------------! end select From 326e6316af9ee25b024ce0c596cc7ecd806e456b Mon Sep 17 00:00:00 2001 From: crollinson Date: Sun, 11 Nov 2018 13:26:02 -0600 Subject: [PATCH 05/11] some more places to add ianth = 2 --- ED/src/dynamics/disturbance.f90 | 2 +- ED/src/init/ed_params.f90 | 4 ++-- ED/src/init/landuse_init.f90 | 4 ++-- ED/src/io/ed_read_ed10_20_history.f90 | 8 ++++++++ ED/src/io/ed_read_ed21_history.F90 | 14 +++++++------- 5 files changed, 20 insertions(+), 12 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index 77e6ea94c..016ee0130 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -222,7 +222,7 @@ subroutine apply_disturbances(cgrid) !------------------------------------------------------------------------------! ! Find the area to be harvested when biomass targets have been ! - ! established. ! + ! established. Returns zero if ianth != 1 ! !------------------------------------------------------------------------------! call find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) !------------------------------------------------------------------------------! diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index f3ebd2749..82926e84b 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -148,12 +148,12 @@ subroutine load_ed_ecosystem_params() ! Warn the user in case the PFT choice for agriculture or forest plantation was ! ! inconsistent. ! !---------------------------------------------------------------------------------------! - if (count(include_pft_ag) == 0 .and. ianth_disturb == 1) then + if (count(include_pft_ag) == 0 .and. ianth_disturb >= 1) then call warning ('PFT defined in agri_stock is not included in include_these_pft,'// & ' your croplands will be barren and not very profitable...' & ,'load_ecosystem_params','ed_params.f90') end if - if (count(include_pft_fp) == 0 .and. ianth_disturb == 1) then + if (count(include_pft_fp) == 0 .and. ianth_disturb >= 1) then call warning ('PFT defined in plantation_stock is not listed in include_these_pft,'//& ' your forest plantation will be barren and not very profitable ...' & ,'load_ecosystem_params','ed_params.f90') diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index 80d81947a..5d6114bef 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -116,7 +116,7 @@ subroutine landuse_init !------------------------------------------------------------------------------------! ! Find the list of disturbance rate files. ! !------------------------------------------------------------------------------------! - if (ianth_disturb == 1) then + if (ianth_disturb >= 1) then call ed_filelist(full_list,lu_database(igr),nflist) call ed1_fileinfo('.lu',nflist,full_list,nfllu,lu_list,llon_list,llat_list) end if @@ -153,7 +153,7 @@ subroutine landuse_init cpoly%plantation(:) = 0 !------------------------------------------------------------------------------! - case (1) + case (1,2) !------------------------------------------------------------------------------! ! Comput the distance between the current polygon and all the files. ! diff --git a/ED/src/io/ed_read_ed10_20_history.f90 b/ED/src/io/ed_read_ed10_20_history.f90 index 1d2c52b13..d10a5b905 100644 --- a/ED/src/io/ed_read_ed10_20_history.f90 +++ b/ED/src/io/ed_read_ed10_20_history.f90 @@ -444,6 +444,10 @@ subroutine read_ed10_ed20_history_file !----- Anthropogenic, assume logging. --------------------------! csite%dist_type (ip2) = 6 !---------------------------------------------------------------! + case (2) + !----- Anthropogenic, assume logging. --------------------------! + csite%dist_type (ip2) = 6 + !---------------------------------------------------------------! end select !------------------------------------------------------------------! case (2) @@ -521,6 +525,10 @@ subroutine read_ed10_ed20_history_file !----- Anthropogenic, assume logging. -----------------------------! csite%dist_type (ip) = 6 !------------------------------------------------------------------! + case (2) + !----- Anthropogenic, assume logging. -----------------------------! + csite%dist_type (ip) = 6 + !------------------------------------------------------------------! end select !---------------------------------------------------------------------! case (2) diff --git a/ED/src/io/ed_read_ed21_history.F90 b/ED/src/io/ed_read_ed21_history.F90 index 6b041d6f8..6d26bc298 100644 --- a/ED/src/io/ed_read_ed21_history.F90 +++ b/ED/src/io/ed_read_ed21_history.F90 @@ -513,7 +513,7 @@ subroutine read_ed21_history_file !---------------------------------------------------------------! if (plantation(ipa) == 0 .and. ianth_disturb == 0) then csite%dist_type(ipa) = 5 - else if (plantation(ipa) == 0 .and. ianth_disturb == 1) then + else if (plantation(ipa) == 0 .and. ianth_disturb >= 1) then csite%dist_type(ipa) = 6 else csite%dist_type(ipa) = 2 @@ -1127,7 +1127,7 @@ subroutine read_ed21_history_unstruct slat_rscl(:) = 100. nlat_rscl(:) = -100. inquire(file=trim(lu_rescale_file(igr)),exist=exists) - rescale_glob = ianth_disturb == 1 .and. exists + rescale_glob = ianth_disturb >= 1 .and. exists nrescale = 0 if (rescale_glob) then open (unit=13,file=trim(lu_rescale_file(igr)),status='old',action='read') @@ -1147,7 +1147,7 @@ subroutine read_ed21_history_unstruct end do readrescale rescale_glob = nrescale > 0 close (unit=13,status='keep') - elseif (ianth_disturb == 1 .and. len_trim(lu_rescale_file(igr)) > 0) then + elseif (ianth_disturb >= 1 .and. len_trim(lu_rescale_file(igr)) > 0) then write (unit=*,fmt='(a)') '-------------------------------------------------------' write (unit=*,fmt='(a)') ' WARNING! WARNING! WARNING! WARNING! WARNING! WARNING! ' write (unit=*,fmt='(a)') ' WARNING! WARNING! WARNING! WARNING! WARNING! WARNING! ' @@ -1675,7 +1675,7 @@ subroutine read_ed21_history_unstruct !---------------------------------------------------------------! if (plantation(ipa) == 0 .and. ianth_disturb == 0) then csite%dist_type(ipa) = 5 - else if (plantation(ipa) == 0 .and. ianth_disturb == 1) then + else if (plantation(ipa) == 0 .and. ianth_disturb >= 1) then csite%dist_type(ipa) = 6 else csite%dist_type(ipa) = 2 @@ -2324,7 +2324,7 @@ subroutine read_ed21_polyclone slat_rscl(:) = 100. nlat_rscl(:) = -100. inquire(file=trim(lu_rescale_file(igr)),exist=exists) - rescale_glob = ianth_disturb == 1 .and. exists + rescale_glob = ianth_disturb >= 1 .and. exists nrescale = 0 if (rescale_glob) then open (unit=13,file=trim(lu_rescale_file(igr)),status='old',action='read') @@ -2344,7 +2344,7 @@ subroutine read_ed21_polyclone end do readrescale rescale_glob = nrescale > 0 close (unit=13,status='keep') - elseif (ianth_disturb == 1 .and. len_trim(lu_rescale_file(igr)) > 0) then + elseif (ianth_disturb >= 1 .and. len_trim(lu_rescale_file(igr)) > 0) then write (unit=*,fmt='(a)') '-------------------------------------------------------' write (unit=*,fmt='(a)') ' WARNING! WARNING! WARNING! WARNING! WARNING! WARNING! ' write (unit=*,fmt='(a)') ' WARNING! WARNING! WARNING! WARNING! WARNING! WARNING! ' @@ -2811,7 +2811,7 @@ subroutine read_ed21_polyclone !---------------------------------------------------------------! if (plantation(ipa) == 0 .and. ianth_disturb == 0) then csite%dist_type(ipa) = 5 - else if (plantation(ipa) == 0 .and. ianth_disturb == 1) then + else if (plantation(ipa) == 0 .and. ianth_disturb >= 1) then csite%dist_type(ipa) = 6 else csite%dist_type(ipa) = 2 From 9a391050649e8b50534bba32f87b2f5fe4e0fc4e Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Wed, 11 Sep 2019 12:10:18 -0500 Subject: [PATCH 06/11] Update Forestry: names changes; add above/below harvest Gradually bringing in some ideas from Marcos to allow for more robust management schemes. Based on Marcos's work; now do not distinguish betwen primary & secondary forest during management & shift some of the terminology to better specify plantation parameters. Also, add new harvest probability to allow for understory thinning. This should have similar function as Marcos's skid trails and collatoral damage idea, but also in a more targeted manner that mirrors some potential management strategies like a targeted understory thin to remove late-successional (maple) stems & promote oak regeneration. (I've probably missed stuff and this won't compile, but it's a start) --- ED/src/dynamics/disturbance.f90 | 78 +++++++++++++-------- ED/src/dynamics/forestry.f90 | 14 +--- ED/src/dynamics/mortality.f90 | 22 ++++-- ED/src/init/landuse_init.f90 | 120 ++++++++++++++++++++------------ ED/src/memory/ed_state_vars.F90 | 22 +++--- 5 files changed, 154 insertions(+), 102 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index 016ee0130..6c8856a56 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -109,7 +109,8 @@ subroutine apply_disturbances(cgrid) real , dimension(n_pft,n_dbh) :: initial_agb real , dimension(n_pft,n_dbh) :: initial_basal_area real , dimension(n_pft) :: mindbh_harvest - real , dimension(n_pft) :: harvprob + real , dimension(n_pft) :: harvprob_g + real , dimension(n_pft) :: harvprob_l real :: pot_area_remain real :: area_fac real :: orig_area @@ -659,7 +660,8 @@ subroutine apply_disturbances(cgrid) disturbed = act_area_loss(ipa,new_lu) > tiny_num dist_path = 0 mindbh_harvest(1:n_pft) = huge(1.) - harvprob (1:n_pft) = 0. + harvprob_g (1:n_pft) = 0. + harvprob_l (1:n_pft) = 0. !---------------------------------------------------------------------! @@ -699,11 +701,13 @@ subroutine apply_disturbances(cgrid) case (6) !----- Logging. ---------------------------------------------------! if (mature_primary) then - mindbh_harvest(1:n_pft) = cpoly%mindbh_primary(1:n_pft,isi) - harvprob (1:n_pft) = cpoly%probharv_primary(1:n_pft,isi) + mindbh_harvest(1:n_pft) = cpoly%mindbh_harvest(1:n_pft,isi) + harvprob_g (1:n_pft) = cpoly%prob_harvest_g(1:n_pft,isi) + harvprob_l (1:n_pft) = cpoly%prob_harvest_l(1:n_pft,isi) else if (mature_secondary) then - mindbh_harvest(1:n_pft) = cpoly%mindbh_secondary(1:n_pft,isi) - harvprob (1:n_pft) = cpoly%probharv_secondary(1:n_pft,isi) + mindbh_harvest(1:n_pft) = cpoly%mindbh_harvest(1:n_pft,isi) + harvprob_g (1:n_pft) = cpoly%prob_harvest_g(1:n_pft,isi) + harvprob_l (1:n_pft) = cpoly%prob_harvest_l(1:n_pft,isi) end if !------------------------------------------------------------------! end select @@ -730,7 +734,7 @@ subroutine apply_disturbances(cgrid) end if !----- Find the mortality associated with disturbance. ------------! call disturbance_mortality(csite,ipa,dist_rate,new_lu & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !------------------------------------------------------------------! @@ -748,9 +752,9 @@ subroutine apply_disturbances(cgrid) area_fac = act_area_loss(ipa,new_lu) / csite%area(onsp+new_lu) call increment_patch_vars(csite,onsp+new_lu,ipa,area_fac) call insert_survivors(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) call accum_dist_litt(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !---------------------------------------------------------------! case (1) !---------------------------------------------------------------! @@ -765,9 +769,9 @@ subroutine apply_disturbances(cgrid) area_fac = act_area_loss(ipa,new_lu)/csite%area(onsp+new_lu) call increment_patch_vars(csite,onsp+new_lu,ipa,area_fac) call insert_survivors(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) call accum_dist_litt(csite,onsp+new_lu,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !------------------------------------------------------------! case default !------------------------------------------------------------! @@ -816,9 +820,9 @@ subroutine apply_disturbances(cgrid) if (same_pft) then call increment_patch_vars(csite,npa,ipa,area_fac) call insert_survivors(csite,npa,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) call accum_dist_litt(csite,npa,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest,harvprob) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) end if !---------------------------------------------------------! end do @@ -1476,12 +1480,21 @@ subroutine site_disturbance_rates(year, cgrid) !------------------------------------------------------------------! if (csite%age(ipa) > plantation_rotation) then cohortloop_02: do ico=1,cpatch%ncohorts - ipft = cpatch%pft(ico) - weight = cpatch%nplant(ico) * cpatch%basarea(ico) & - * csite%area(ipa) - pharvest (ilu) = pharvest(ilu) & - + cpoly%probharv_secondary(ipft,isi) * weight - sumweight(ilu) = sumweight(ilu) + weight + ipft = cpatch%pft(ico) + if (cpatch%dbh(ico) >= cpoly%mindbh_harvest(ipft,isi)) then + weight = cpatch%nplant(ico) * cpatch%basarea(ico) & + * csite%area(ipa) + pharvest (ilu) = pharvest(ilu) & + + cpoly%prob_harvest_g(ipft,isi) & + * weight + sumweight(ilu) = sumweight(ilu) + weight + else + weight = cpatch%nplant(ico) * cpatch%basarea(ico) & + * csite%area(ipa) + pharvest (ilu) = pharvest(ilu) & + + cpoly%prob_harvest_l(ipft,isi) & + * weight + sumweight(ilu) = sumweight(ilu) + weight end do cohortloop_02 !---------------------------------------------------------------! end if @@ -1494,13 +1507,20 @@ subroutine site_disturbance_rates(year, cgrid) if (csite%age(ipa) > mature_harvest_age) then cohortloop_06: do ico=1,cpatch%ncohorts ipft = cpatch%pft(ico) - if (cpatch%dbh(ico) >= cpoly%mindbh_secondary(ipft,isi)) then + if (cpatch%dbh(ico) >= cpoly%mindbh_harvest(ipft,isi)) then weight = cpatch%nplant(ico) * cpatch%basarea(ico) & * csite%area(ipa) pharvest (ilu) = pharvest(ilu) & - + cpoly%probharv_secondary(ipft,isi) & + + cpoly%prob_harvest_g(ipft,isi) & * weight sumweight(ilu) = sumweight(ilu) + weight + else + weight = cpatch%nplant(ico) * cpatch%basarea(ico) & + * csite%area(ipa) + pharvest (ilu) = pharvest(ilu) & + + cpoly%prob_harvest_l(ipft,isi) & + * weight + sumweight(ilu) = sumweight(ilu) + weight end if end do cohortloop_06 !---------------------------------------------------------------! @@ -2718,7 +2738,7 @@ end subroutine increment_patch_vars ! This subroutine will populate the disturbed patch with the cohorts that were ! ! disturbed but did not go extinct. ! !---------------------------------------------------------------------------------------! - subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob) + subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob_g, harvprob_l) use update_derived_props_module use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure @@ -2733,7 +2753,8 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest integer , intent(in) :: np integer , intent(in) :: cp real , dimension(n_pft), intent(in) :: mindbh_harvest - real , dimension(n_pft), intent(in) :: harvprob + real , dimension(n_pft), intent(in) :: harvprob_g + real , dimension(n_pft), intent(in) :: harvprob_l real , intent(in) :: area_fac !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch @@ -2763,7 +2784,7 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest mask(:) = .false. survivalloop: do ico = 1,cpatch%ncohorts - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) & + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g, harvprob_l,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2799,7 +2820,7 @@ subroutine insert_survivors(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest cohortloop: do ico = 1,cpatch%ncohorts - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) & + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g, harvprob_l,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2845,7 +2866,7 @@ end subroutine insert_survivors !=======================================================================================! ! This subroutine updates the litter pools after a disturbance takes place. ! !---------------------------------------------------------------------------------------! - subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob) + subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest,harvprob_g, harvprob_l) use ed_state_vars, only : sitetype & ! structure , patchtype & ! structure , polygontype ! ! structure @@ -2864,7 +2885,8 @@ subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest, integer , intent(in) :: np integer , intent(in) :: cp real , dimension(n_pft), intent(in) :: mindbh_harvest - real , dimension(n_pft), intent(in) :: harvprob + real , dimension(n_pft), intent(in) :: harvprob_g + real , dimension(n_pft), intent(in) :: harvprob_l integer , intent(in) :: new_lu real , intent(in) :: area_fac integer , intent(in) :: dist_path @@ -2912,7 +2934,7 @@ subroutine accum_dist_litt(csite,np,cp,new_lu,area_fac,dist_path,mindbh_harvest, !---------------------------------------------------------------------------------! !----- Find survivorship. --------------------------------------------------------! - survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) + survival_fac = survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g,harvprob_l,cpatch,ico) !---------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 193f0141e..ddbd7a9f3 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -117,17 +117,9 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) !----- Agriculture. No timber harvesting here. ----------------------------------! mindbh_harvest(:) = huge(1.) !---------------------------------------------------------------------------------! - case (2) - !----- Forest plantation. Usually all biomass is cleared. -----------------------! - mindbh_harvest(:) = 0. - !---------------------------------------------------------------------------------! - case (4,5,6) - !----- Secondary vegetation. -----------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_secondary(:,isi) - !---------------------------------------------------------------------------------! - case (3) - !----- Primary vegetation. -------------------------------------------------------! - mindbh_harvest(:) = cpoly%mindbh_primary (:,isi) + case (2:6) + !----- Harvesting of some sort. ALl use same variable now -----------------------! + mindbh_harvest(:) = cpoly%mindbh_harvest(:,isi) !---------------------------------------------------------------------------------! end select !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index b4e1e08d7..110a058dc 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -109,7 +109,7 @@ end subroutine mortality_rates ! disturbance. ! !---------------------------------------------------------------------------------------! subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path & - ,mindbh_harvest,harvprob) + ,mindbh_harvest,harvprob_g,harvprob_l) use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure use ed_max_dims , only : n_pft ! ! intent(in) @@ -121,7 +121,8 @@ subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path integer , intent(in) :: new_lu integer , intent(in) :: dist_path real , dimension(n_pft), intent(in) :: mindbh_harvest - real , dimension(n_pft), intent(in) :: harvprob + real , dimension(n_pft), intent(in) :: harvprob_g + real , dimension(n_pft), intent(in) :: harvprob_l !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch integer :: ico @@ -164,7 +165,7 @@ end subroutine disturbance_mortality ! -- cpatch: current patch. ! ! -- ico: index for current cohort. ! !---------------------------------------------------------------------------------------! - real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) + real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g,harvprob_l,cpatch,ico) use ed_state_vars, only : patchtype ! ! structure use disturb_coms , only : treefall_hite_threshold & ! intent(in) , fire_hite_threshold & ! intent(in) @@ -179,7 +180,8 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) !----- Arguments. -------------------------------------------------------------------! type(patchtype) , target :: cpatch real , dimension(n_pft), intent(in) :: mindbh_harvest - real , dimension(n_pft), intent(in) :: harvprob + real , dimension(n_pft), intent(in) :: harvprob_g + real , dimension(n_pft), intent(in) :: harvprob_l integer , intent(in) :: ico integer , intent(in) :: new_lu integer , intent(in) :: dist_path @@ -277,7 +279,12 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) case (6) !---------------------------------------------------------------------------------! - ! Logging. At this point a single pathway exists: cohorts above threshold ! + ! Logging. ! + ! NEW: Survivorship = inverse fraction of amount removed above and below the hite ! + ! threshold. This allows for understory thinning as a management strategy ! + ! ! + !---------------------------------------------------------------------------------! + ! OLD: At this point a single pathway exists: cohorts above threshold ! ! are completely removed, and small cohorts have the same survivorship as small ! ! cohorts at a treefall site. Both could be re-visited in the future, and ! ! different pathways for conventional and reduced-impact logging could be ! @@ -285,12 +292,13 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) !---------------------------------------------------------------------------------! if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then if (ianth_disturb == 2) then - survivorship = 1 - harvprob(ipft) + survivorship = 1 - harvprob_g(ipft) else survivorship = 0.0 end if else - survivorship = treefall_s_ltht(ipft) + survivorship = harvprob_l(ipft) + ! survivorship = treefall_s_ltht(ipft) ! OLD version end if !---------------------------------------------------------------------------------! end select diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index 5d6114bef..845519b16 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -69,10 +69,12 @@ subroutine landuse_init integer :: yd_this integer :: yd_last logical :: inside - real , dimension(n_pft) :: mindbh_1ary - real , dimension(n_pft) :: harvprob_1ary - real , dimension(n_pft) :: mindbh_2ary - real , dimension(n_pft) :: harvprob_2ary + real , dimension(n_pft) :: mindbh_slog + real , dimension(n_pft) :: harvprob_slog_g + real , dimension(n_pft) :: harvprob_slog_l + real , dimension(n_pft) :: mindbh_fplt + real , dimension(n_pft) :: harvprob_fplt_g + real , dimension(n_pft) :: harvprob_fplt_l real , dimension(num_lu_trans) :: landuse_now real :: lu_area real :: lu_area_i @@ -139,10 +141,9 @@ subroutine landuse_init !----- Set the parameters in a way that no logging/ploughing will happen. -----! do isi = 1,cpoly%nsites cpoly%num_landuse_years(isi) = 1 - cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh - cpoly%probharv_primary (1:n_pft,isi) = 0. - cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh - cpoly%probharv_secondary(1:n_pft,isi) = 0. + cpoly%mindbh_harvest (1:n_pft,isi) = huge_dbh + cpoly%prob_havest_g (1:n_pft,isi) = 0. + cpoly%prob_havest_l (1:n_pft,isi) = 0. cpoly%clutimes(1,isi)%landuse_year = iyeara cpoly%clutimes(1,isi)%landuse(1:num_lu_trans) = 0.0 end do @@ -182,10 +183,12 @@ subroutine landuse_init ! harvesting. The actual variables will be read in the following block. ! !------------------------------------------------------------------------------! harvest_pft(1:n_pft) = -1 - mindbh_1ary(1:n_pft) = huge_dbh - harvprob_1ary(1:n_pft) = 0. - mindbh_2ary(1:n_pft) = huge_dbh - harvprob_2ary(1:n_pft) = 0. + mindbh_slog(1:n_pft) = huge_dbh + harvprob_slog_g(1:n_pft) = 0. + harvprob_slog_l(1:n_pft) = 0. + mindbh_fplt(1:n_pft) = huge_dbh + harvprob_fplt_g(1:n_pft) = 0. + harvprob_fplt_l(1:n_pft) = 0. !----- Initialise plantation patches if plantation information is available. --! cpoly%plantation(:) = 0 @@ -239,19 +242,27 @@ subroutine landuse_init read (unit=12,fmt=hform) cdum cdum = cdum(hoff:) - read (cdum, fmt=*) (mindbh_1ary(h) ,h=1,nharvest) + read (cdum, fmt=*) (mindbh_slog(h) ,h=1,nharvest) read (unit=12,fmt=hform) cdum cdum = cdum(hoff:) - read (cdum, fmt=*) (harvprob_1ary(h) ,h=1,nharvest) + read (cdum, fmt=*) (harvprob_slog_g(h) ,h=1,nharvest) read (unit=12,fmt=hform) cdum cdum = cdum(hoff:) - read (cdum, fmt=*) (mindbh_2ary(h) ,h=1,nharvest) + read (cdum, fmt=*) (harvprob_slog_l(h) ,h=1,nharvest) read (unit=12,fmt=hform) cdum cdum = cdum(hoff:) - read (cdum, fmt=*) (harvprob_2ary(h) ,h=1,nharvest) + read (cdum, fmt=*) (mindbh_fplt(h) ,h=1,nharvest) + + read (unit=12,fmt=hform) cdum + cdum = cdum(hoff:) + read (cdum, fmt=*) (harvprob_fplt_g(h) ,h=1,nharvest) + + read (unit=12,fmt=hform) cdum + cdum = cdum(hoff:) + read (cdum, fmt=*) (harvprob_fplt_l(h) ,h=1,nharvest) else !---------------------------------------------------------------------------! ! No specific PFT information was given, this is likely to be a case in ! @@ -262,8 +273,8 @@ subroutine landuse_init nopftloop: do ipft=1,n_pft h=h+1 harvest_pft(h) = ipft - mindbh_1ary(1:n_pft) = 0. - mindbh_2ary(1:n_pft) = 0. + mindbh_slog(1:n_pft) = 0. + mindbh_fplt(1:n_pft) = 0. end do nopftloop end if read (unit=12,fmt=*) @@ -311,23 +322,48 @@ subroutine landuse_init !----- Initialise the PFT-dependent arrays. --------------------------------! - cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh - cpoly%probharv_primary (1:n_pft,isi) = 0. - cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh - cpoly%probharv_secondary(1:n_pft,isi) = 0. - - !----- Fill the arrays with the appropriate PFT. ---------------------------! - do ipft=1,n_pft - harvloop: do h=1,nharvest - if (harvest_pft(h) == ipft) then - cpoly%mindbh_primary (ipft,isi) = mindbh_1ary (h) - cpoly%probharv_primary (ipft,isi) = harvprob_1ary(h) - cpoly%mindbh_secondary (ipft,isi) = mindbh_2ary (h) - cpoly%probharv_secondary(ipft,isi) = harvprob_2ary(h) - exit harvloop - end if - end do harvloop - end do + cpoly%mindbh_harvest(1:n_pft,isi) = huge_dbh + cpoly%prob_harvest_g(1:n_pft,isi) = 0. + cpoly%prob_harvest_l(1:n_pft,isi) = 0. + + ! From Marcos + !----- Fill the arrays with the appropriate PFT. ------------------------! + select case(cpoly%plantation(isi)) + case (0) + harvloop_slog: do h=1,nharvest + ipft = harvest_pft(h) + if (ipft >= 1 .and. ipft <= n_pft) then + cpoly%mindbh_harvest(ipft,isi) = mindbh_slog (h) + cpoly%prob_harvest_g (ipft,isi) = harvprob_slog_g(h) + cpoly%prob_harvest_l (ipft,isi) = harvprob_slog_l(h) + end if + end do harvloop_slog + case (1) + harvloop_fplt: do h=1,nharvest + ipft = harvest_pft(h) + if (ipft >= 1 .and. ipft <= n_pft) then + cpoly%mindbh_harvest(ipft,isi) = mindbh_fplt (h) + cpoly%prob_harvest_g (ipft,isi) = harvprob_fplt_g(h) + cpoly%prob_harvest_l (ipft,isi) = harvprob_fplt_l(h) + end if + end do harvloop_fplt + end select + !------------------------------------------------------------------------! + +! +! ! Original +! !----- Fill the arrays with the appropriate PFT. ---------------------------! +! do ipft=1,n_pft +! harvloop: do h=1,nharvest +! if (harvest_pft(h) == ipft) then +! cpoly%mindbh_primary (ipft,isi) = mindbh_slog (h) +! cpoly%probharv_primary (ipft,isi) = harvprob_slog(h) +! cpoly%mindbh_secondary (ipft,isi) = mindbh_fplt (h) +! cpoly%probharv_secondary(ipft,isi) = harvprob_fplt(h) +! exit harvloop +! end if +! end do harvloop +! end do !----- Padding disturbances with zero before first available lu year. ------! @@ -382,10 +418,9 @@ subroutine landuse_init cpoly%num_landuse_years(isi) = cpoly%num_landuse_years(1) !----- PFT-dependent harvest characteristics. ------------------------! - cpoly%mindbh_primary (:,isi) = cpoly%mindbh_primary (:,1) - cpoly%probharv_primary (:,isi) = cpoly%probharv_primary (:,1) - cpoly%mindbh_secondary (:,isi) = cpoly%mindbh_secondary (:,1) - cpoly%probharv_secondary(:,isi) = cpoly%probharv_secondary(:,1) + cpoly%mindbh_harvest(:,isi) = cpoly%mindbh_harvest(:,1) + cpoly%prob_harvest_g(:,isi) = cpoly%prob_harvest_g(:,1) + cpoly%prob_harvest_l(:,isi) = cpoly%prob_harvest_l(:,1) !----- Disturbances. ----------------------------------------------------! do iyear = 1,cpoly%num_landuse_years(isi) @@ -422,10 +457,9 @@ subroutine landuse_init !----- Set the parameters in a way that no logging/ploughing will happen. --! do isi = 1,cpoly%nsites cpoly%num_landuse_years(isi) = 1 - cpoly%mindbh_primary (1:n_pft,isi) = huge_dbh - cpoly%probharv_primary (1:n_pft,isi) = 0. - cpoly%mindbh_secondary (1:n_pft,isi) = huge_dbh - cpoly%probharv_secondary(1:n_pft,isi) = 0. + cpoly%mindbh_harvest(1:n_pft,isi) = huge_dbh + cpoly%prob_harvest_g(1:n_pft,isi) = 0. + cpoly%prob_harvest_l(1:n_pft,isi) = 0. cpoly%clutimes(1,isi)%landuse_year = iyeara cpoly%clutimes(1,isi)%landuse(1:num_lu_trans) = 0.0 end do diff --git a/ED/src/memory/ed_state_vars.F90 b/ED/src/memory/ed_state_vars.F90 index 05f71cf41..c0e118441 100644 --- a/ED/src/memory/ed_state_vars.F90 +++ b/ED/src/memory/ed_state_vars.F90 @@ -1964,12 +1964,10 @@ module ed_state_vars type(lutime), pointer,dimension(:,:) :: clutimes !(luyears,nsites) ! Date: Wed, 11 Sep 2019 12:30:48 -0500 Subject: [PATCH 07/11] Fix typos that prevent compilation --- ED/src/dynamics/disturbance.f90 | 1 + ED/src/dynamics/mortality.f90 | 2 +- ED/src/init/landuse_init.f90 | 6 +++--- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index 6c8856a56..f5f8fb59c 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -1495,6 +1495,7 @@ subroutine site_disturbance_rates(year, cgrid) + cpoly%prob_harvest_l(ipft,isi) & * weight sumweight(ilu) = sumweight(ilu) + weight + end if end do cohortloop_02 !---------------------------------------------------------------! end if diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index 110a058dc..20dc8707f 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -131,7 +131,7 @@ subroutine disturbance_mortality(csite,ipa,disturbance_rate,new_lu,dist_path cpatch => csite%patch(ipa) do ico=1,cpatch%ncohorts - f_survival = survivorship(new_lu,dist_path,mindbh_harvest,harvprob,cpatch,ico) + f_survival = survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g,harvprob_l,cpatch,ico) cpatch%mort_rate(5,ico) = cpatch%mort_rate(5,ico) & - log( f_survival & + (1.0 - f_survival) * exp(- disturbance_rate) ) diff --git a/ED/src/init/landuse_init.f90 b/ED/src/init/landuse_init.f90 index 845519b16..2347dbda3 100644 --- a/ED/src/init/landuse_init.f90 +++ b/ED/src/init/landuse_init.f90 @@ -142,8 +142,8 @@ subroutine landuse_init do isi = 1,cpoly%nsites cpoly%num_landuse_years(isi) = 1 cpoly%mindbh_harvest (1:n_pft,isi) = huge_dbh - cpoly%prob_havest_g (1:n_pft,isi) = 0. - cpoly%prob_havest_l (1:n_pft,isi) = 0. + cpoly%prob_harvest_g (1:n_pft,isi) = 0. + cpoly%prob_harvest_l (1:n_pft,isi) = 0. cpoly%clutimes(1,isi)%landuse_year = iyeara cpoly%clutimes(1,isi)%landuse(1:num_lu_trans) = 0.0 end do @@ -707,4 +707,4 @@ end subroutine read_plantation_fractions !==========================================================================================! !==========================================================================================! -end module landuse_init_module \ No newline at end of file +end module landuse_init_module From 064f39200bb06886406cf905bb3e4629bc656c7c Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Thu, 12 Sep 2019 12:38:50 -0500 Subject: [PATCH 08/11] fix above/below bug --- ED/src/dynamics/mortality.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index 20dc8707f..eb0bf36a6 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -293,11 +293,11 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g,harvprob_l if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then if (ianth_disturb == 2) then survivorship = 1 - harvprob_g(ipft) - else - survivorship = 0.0 - end if + ! else + ! survivorship = 0.0 + ! end if else - survivorship = harvprob_l(ipft) + survivorship = 1 - harvprob_l(ipft) ! survivorship = treefall_s_ltht(ipft) ! OLD version end if !---------------------------------------------------------------------------------! From 0ed954081ac859a8b4275754d5729ebbd3da1e74 Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Thu, 12 Sep 2019 13:49:36 -0500 Subject: [PATCH 09/11] Bug fix to compile --- ED/src/dynamics/mortality.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index eb0bf36a6..c40cea663 100644 --- a/ED/src/dynamics/mortality.f90 +++ b/ED/src/dynamics/mortality.f90 @@ -291,7 +291,7 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,harvprob_g,harvprob_l ! applied. ! !---------------------------------------------------------------------------------! if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then - if (ianth_disturb == 2) then + !if (ianth_disturb == 2) then survivorship = 1 - harvprob_g(ipft) ! else ! survivorship = 0.0 From 4ca39eb5064e9a53046dd016943f09f9c4ff3fd5 Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Mon, 16 Sep 2019 11:56:04 -0500 Subject: [PATCH 10/11] Quick forestry fixes Reset number of LU dimensions to 6 since we're not adding new types; fix references in harvest folder so that we continue to repurpose secondary forest for plantation and primary for 'natural' forest. --- ED/src/dynamics/disturbance.f90 | 5 +++-- ED/src/memory/ed_max_dims.F90 | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index f5f8fb59c..2a5c8e037 100644 --- a/ED/src/dynamics/disturbance.f90 +++ b/ED/src/dynamics/disturbance.f90 @@ -1610,11 +1610,12 @@ subroutine site_disturbance_rates(year, cgrid) ! Logging based on grid fraction and tree size. Harvest creates new ! ! patches and target biomass harvest set to zero. ! !------------------------------------------------------------------------! - cpoly%disturbance_rates(6,6,isi) = clutime%landuse(13) - cpoly%disturbance_rates(2,2,isi) = clutime%landuse(15) + cpoly%disturbance_rates(6,6,isi) = clutime%landuse(15) + cpoly%disturbance_rates(2,2,isi) = clutime%landuse(13) cpoly%primary_harvest_target (isi) = 0.0 cpoly%secondary_harvest_target (isi) = 0.0 + !------------------------------------------------------------------------! end select ! End ianth select diff --git a/ED/src/memory/ed_max_dims.F90 b/ED/src/memory/ed_max_dims.F90 index bf82ee03f..c46366cf0 100644 --- a/ED/src/memory/ed_max_dims.F90 +++ b/ED/src/memory/ed_max_dims.F90 @@ -183,7 +183,7 @@ module ed_max_dims ! 7 -- Understory Thin -- Logged forest (skid trail + road). ! ! 8 -- Cropland. ! !---------------------------------------------------------------------------------------! - integer, parameter :: n_dist_types = 8 ! Changed from 6 to include understory thin + integer, parameter :: n_dist_types = 6 ! Changed from 6 to include understory thin !---------------------------------------------------------------------------------------! From 9dcb2934e77fa1ad325b802b4fecc64c58a32841 Mon Sep 17 00:00:00 2001 From: Christy Rollinson Date: Mon, 16 Sep 2019 11:57:10 -0500 Subject: [PATCH 11/11] Forestry fix: need to allocate potential harvest area if we don't define the potential harvest area in this routine, harvest ends up geting rescaled. This is clunky code, but was trying to make it consistent with the older, biomass-based harvest schemes. --- ED/src/dynamics/forestry.f90 | 66 ++++++++++++++++++++++++++++++++++-- 1 file changed, 63 insertions(+), 3 deletions(-) diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index ddbd7a9f3..dd38c13ac 100644 --- a/ED/src/dynamics/forestry.f90 +++ b/ED/src/dynamics/forestry.f90 @@ -53,8 +53,7 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) real :: site_harvest_target real :: site_harvestable_agb ! real :: site_hvmax_btimber -! real :: site_hvpot_btimber - +! real :: site_hvpot_btimber real :: area_mature_primary real :: hvagb_mature_primary real :: area_mature_secondary @@ -77,7 +76,68 @@ subroutine find_harvest_area(cpoly,isi,onsp,harvestable_agb,pot_area_harv) case (0) ! Nothing to do because anthropogenic disturbance is turned off return case (2) ! Harvest is based on size + area - return + + !----- Link to the current site. -------------------------------------------------------! + csite => cpoly%site(isi) + !---------------------------------------------------------------------------------------! + + cpoly%primary_harvest_memory (isi) = 0.0 + cpoly%secondary_harvest_memory(isi) = 0.0 + + lambda_mature_plantation = cpoly%disturbance_rates(2,2,isi) + lambda_mature_primary = cpoly%disturbance_rates(6,6,isi) + + !---------------------------------------------------------------------------------------! + ! Loop over patches. ! + !---------------------------------------------------------------------------------------! + patch_loop: do ipa=1,onsp + + !------------------------------------------------------------------------------------! + ! Find out whether to harvest this patch. ! + !------------------------------------------------------------------------------------! + select case (csite%dist_type(ipa)) + case (2) + !----- Forest plantation. --------------------------------------------------------! + if ( csite%age(ipa) > plantation_rotation ) then + pot_area_harv(ipa) = csite%area(ipa) * lambda_mature_plantation + end if + !---------------------------------------------------------------------------------! + case (6) + !----- Primary/Secondary forest. -------------------------------------------------! + if ( csite%age(ipa) > mature_harvest_age ) then + pot_area_harv(ipa) = csite%area(ipa) * cpoly%disturbance_rates(6,6,isi) + end if + case default + !----- Agriculture. Do not log. -------------------------------------------------! + continue + !---------------------------------------------------------------------------------! + end select + !------------------------------------------------------------------------------------! + end do patch_loop + !---------------------------------------------------------------------------------------! + + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' FORESTRY. HARVEST RATES' + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a,1x,i5)') ' ISI = ',isi + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PRIMARY) = ' & + , lambda_mature_primary + write (unit=*,fmt='(a,1x,es12.5)') ' HV LAMBDA (PLANTATION) = ' & + , lambda_mature_plantation + write (unit=*,fmt='(a)' ) ' ' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(5(a,1x))' ) ' IPA',' LU',' AGE',' AREA' & + ,' HV_AREA' + write (unit=*,fmt='(a)' ) '------------------------------------------------' + do ipa=1,onsp + write (unit=*,fmt='(2(i5,1x),3(f12.7,1x))') ipa,csite%dist_type(ipa) & + ,csite%age(ipa),csite%area(ipa) & + ,pot_area_harv(ipa) + end do + write (unit=*,fmt='(a)' ) '------------------------------------------------' + write (unit=*,fmt='(a)' ) ' ' + !---------------------------------------------------------------------------------------! case (1)