diff --git a/ED/src/dynamics/disturbance.f90 b/ED/src/dynamics/disturbance.f90 index bc3671c28..2a5c8e037 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) @@ -108,6 +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_g + real , dimension(n_pft) :: harvprob_l real :: pot_area_remain real :: area_fac real :: orig_area @@ -138,9 +141,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 !------------------------------------------------------------------------------------! @@ -220,7 +223,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) !------------------------------------------------------------------------------! @@ -657,6 +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_g (1:n_pft) = 0. + harvprob_l (1:n_pft) = 0. !---------------------------------------------------------------------! @@ -696,9 +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) + 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) + 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 @@ -725,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) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !------------------------------------------------------------------! @@ -743,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) + ,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) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !---------------------------------------------------------------! case (1) !---------------------------------------------------------------! @@ -760,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) + ,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) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) !------------------------------------------------------------! case default !------------------------------------------------------------! @@ -811,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) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) call accum_dist_litt(csite,npa,ipa,new_lu,area_fac & - ,dist_path,mindbh_harvest) + ,dist_path,mindbh_harvest,harvprob_g,harvprob_l) end if !---------------------------------------------------------! end do @@ -1361,13 +1370,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 +1452,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. ---------! @@ -1470,12 +1480,22 @@ 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 if end do cohortloop_02 !---------------------------------------------------------------! end if @@ -1488,13 +1508,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 !---------------------------------------------------------------! @@ -1530,7 +1557,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 +1605,20 @@ 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(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 !------------------------------------------------------------------------------! end do siteloop @@ -2700,7 +2740,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_g, harvprob_l) use update_derived_props_module use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure @@ -2715,6 +2755,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_g + real , dimension(n_pft), intent(in) :: harvprob_l real , intent(in) :: area_fac !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch @@ -2744,7 +2786,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_g, harvprob_l,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2780,7 +2822,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_g, harvprob_l,cpatch,ico) & * area_fac n_survivors = cpatch%nplant(ico) * survival_fac @@ -2826,7 +2868,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_g, harvprob_l) use ed_state_vars, only : sitetype & ! structure , patchtype & ! structure , polygontype ! ! structure @@ -2845,6 +2887,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_g + real , dimension(n_pft), intent(in) :: harvprob_l integer , intent(in) :: new_lu real , intent(in) :: area_fac integer , intent(in) :: dist_path @@ -2892,7 +2936,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_g,harvprob_l,cpatch,ico) !---------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/forestry.f90 b/ED/src/dynamics/forestry.f90 index 4bee2eff8..dd38c13ac 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,8 @@ 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,236 +64,295 @@ 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. !---------------------------------------------------------------------------------------! + 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 - !----- 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 - !------------------------------------------------------------------------------------! + !----- 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) - !------------------------------------------------------------------------------------! - ! 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 - !---------------------------------------------------------------------------------------! - + !---------------------------------------------------------------------------------------! + ! 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) + !----- 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:6) + !----- Harvesting of some sort. ALl use same variable now -----------------------! + mindbh_harvest(:) = cpoly%mindbh_harvest(:,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 + !---------------------------------------------------------------------------------------! + end select return end subroutine find_harvest_area !==========================================================================================! diff --git a/ED/src/dynamics/mortality.f90 b/ED/src/dynamics/mortality.f90 index c75294002..c40cea663 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_g,harvprob_l) use ed_state_vars, only : sitetype & ! structure , patchtype ! ! structure use ed_max_dims , only : n_pft ! ! intent(in) @@ -121,6 +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_g + real , dimension(n_pft), intent(in) :: harvprob_l !----- Local variables. -------------------------------------------------------------! type(patchtype) , pointer :: cpatch integer :: ico @@ -129,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,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) ) @@ -159,13 +161,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_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) + , 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 +180,8 @@ 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_g + real , dimension(n_pft), intent(in) :: harvprob_l integer , intent(in) :: ico integer , intent(in) :: new_lu integer , intent(in) :: dist_path @@ -273,16 +279,26 @@ real function survivorship(new_lu,dist_path,mindbh_harvest,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 ! ! applied. ! !---------------------------------------------------------------------------------! if (cpatch%dbh(ico) >= mindbh_harvest(ipft)) then - survivorship = 0.0 + !if (ianth_disturb == 2) then + survivorship = 1 - harvprob_g(ipft) + ! else + ! survivorship = 0.0 + ! end if else - survivorship = treefall_s_ltht(ipft) + survivorship = 1 - harvprob_l(ipft) + ! survivorship = treefall_s_ltht(ipft) ! OLD version end if !---------------------------------------------------------------------------------! end select diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index a776e232c..f6372cd04 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 b819065de..2347dbda3 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 @@ -59,10 +69,13 @@ 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 real :: wlon @@ -76,14 +89,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 @@ -100,7 +118,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 @@ -112,7 +130,31 @@ subroutine landuse_init cpoly => cgrid%polygon(ipy) select case (ianth_disturb) - case (1) + 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_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 + !------------------------------------------------------------------------------! + + + !----- No plantations. --------------------------------------------------------! + cpoly%plantation(:) = 0 + !------------------------------------------------------------------------------! + + case (1,2) !------------------------------------------------------------------------------! ! Comput the distance between the current polygon and all the files. ! @@ -141,10 +183,17 @@ 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 + 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,')' @@ -193,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_slog_g(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_l(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=*) (mindbh_fplt(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=*) (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 ! @@ -216,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=*) @@ -265,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. ------! @@ -311,6 +393,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 @@ -334,6 +417,10 @@ 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_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) @@ -370,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 @@ -385,33 +471,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 @@ -541,4 +707,4 @@ end subroutine read_plantation_fractions !==========================================================================================! !==========================================================================================! -end module landuse_init_module \ No newline at end of file +end module landuse_init_module 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/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 diff --git a/ED/src/memory/disturb_coms.f90 b/ED/src/memory/disturb_coms.f90 index 7b6b7a3e0..1d5b4ff3c 100644 --- a/ED/src/memory/disturb_coms.f90 +++ b/ED/src/memory/disturb_coms.f90 @@ -59,9 +59,22 @@ 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. ! + ! 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 !---------------------------------------------------------------------------------------! @@ -113,6 +126,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..c46366cf0 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 = 6 ! 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..c0e118441 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 !