From 7223c9b5fcb2f2c9bc8053a10a00500407cd61b7 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Thu, 9 Jan 2020 10:47:39 -0800 Subject: [PATCH 01/27] Squashed commit of the following: commit 0bbb2c78b57af37940b1e7740cce74074cbe22da Author: Jessica Needham Date: Thu Dec 19 15:53:46 2019 -0800 Squashed commit of the following: commit 194f0739b94f67a26a0280a73d1b52a7502d9f42 Author: Jessica Needham Date: Tue Dec 17 14:47:53 2019 -0800 [ Fix bug in error message] [Error message if user tries to have age dependent mortality and inventory init on. ] Fixes: [NGT-ED Github issue #] User interface changes?: [ No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing commit 34019b7b1613f73bf8abc2d2006069668ccfbc69 Author: Jessica Needham Date: Tue Dec 17 12:09:06 2019 -0800 [ Update FatesPFTIndexswapper to be compatible with python3] [Change iteritems to items ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:Doesn't change fates. No testing. commit ed84697d9c1a136cb7dcaa5f3d6e2badce2a89b9 Author: Jessica Needham Date: Tue Dec 17 11:44:47 2019 -0800 [Size- and age-dependent mortality functionality ] [These changes introduce two new mortality terms: M9 - size dependent and M10 age dependent mortality. They are logistic functions that gradually increase the rate of mortality with size or age. Each is controlled by two parameters a rate and an inflection point. To turn off size dependent mortality the inflection point can be set to _. To turn off age mortality set the fates_cohort_age_fusion_tol parameter to _. Cohort age will simply stay at 0. Otherwise, cohort age is tracked. Cohorts cannot be fused if they are too different in age. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes - five new parameters. User has to specify age and size mortality as either on or off by setting fates_mort_ip_size_senescence and fates_cohort_age_fusion_tol. If they decide to have size or age on, then sensible parameter values need to be selected for inflection point and rate.] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: No testing yet. --- biogeochem/EDCanopyStructureMod.F90 | 17 +- biogeochem/EDCohortDynamicsMod.F90 | 105 +++++++-- biogeochem/EDMortalityFunctionsMod.F90 | 89 +++++++- biogeochem/EDPatchDynamicsMod.F90 | 28 ++- biogeochem/EDPhysiologyMod.F90 | 4 +- main/EDInitMod.F90 | 8 +- main/EDMainMod.F90 | 34 ++- main/EDParamsMod.F90 | 44 +++- main/EDPftvarcon.F90 | 44 +++- main/EDTypesMod.F90 | 29 ++- main/FatesConstantsMod.F90 | 1 + main/FatesHistoryInterfaceMod.F90 | 271 ++++++++++++++++++++--- main/FatesHistoryVariableType.F90 | 15 +- main/FatesIODimensionsMod.F90 | 16 +- main/FatesIOVariableKindMod.F90 | 2 + main/FatesInterfaceMod.F90 | 49 +++- main/FatesInventoryInitMod.F90 | 4 +- main/FatesParametersInterface.F90 | 3 +- main/FatesRestartInterfaceMod.F90 | 70 +++++- main/FatesSizeAgeTypeIndicesMod.F90 | 35 +++ parameter_files/fates_params_default.cdl | 28 ++- tools/FatesPFTIndexSwapper.py | 7 +- tools/modify_fates_paramfile.py | 4 +- 23 files changed, 780 insertions(+), 127 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index bfa50250a9..8f473d440a 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -29,7 +29,8 @@ module EDCanopyStructureMod use FatesInterfaceMod , only : numpft use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use EDTypesMod , only : maxCohortsPerPatch - + use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : all_carbon_elements use PRTGenericMod, only : leaf_organ @@ -1258,9 +1259,11 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesInterfaceMod , only : bc_in_type use EDPatchDynamicsMod , only : set_patchno use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index + use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst - + use EDParamsMod , only : ED_val_cohort_age_fusion_tol + ! !ARGUMENTS integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) @@ -1279,8 +1282,9 @@ subroutine canopy_summarization( nsites, sites, bc_in ) real(r8) :: sapw_c ! sapwood carbon [kg] real(r8) :: store_c ! storage carbon [kg] real(r8) :: struct_c ! structure carbon [kg] - !---------------------------------------------------------------------- + !---------------------------------------------------------------------- + if ( debug ) then write(fates_log(),*) 'in canopy_summarization' endif @@ -1319,8 +1323,13 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ! Update the cohort's index within the size bin classes ! Update the cohort's index within the SCPF classification system call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & - currentCohort%size_class,currentCohort%size_by_pft_class) + currentCohort%size_class,currentCohort%size_by_pft_class) + if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + call coagetype_class_index(currentCohort%coage,currentCohort%pft, & + currentCohort%coage_class,currentCohort%coage_by_pft_class) + end if + call carea_allom(currentCohort%dbh,currentCohort%n,sites(s)%spread,& currentCohort%pft,currentCohort%c_area) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 0339aa4ea8..7129a4d034 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -36,6 +36,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : ican_upper use EDTypesMod , only : site_fluxdiags_type use EDTypesMod , only : num_elements + use EDParamsMod , only : ED_val_cohort_age_fusion_tol use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -50,6 +51,7 @@ module EDCohortDynamicsMod use FatesPlantHydraulicsMod, only : SavePreviousCompartmentVolumes use FatesPlantHydraulicsMod, only : ConstrainRecruitNumber use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index + use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : bfineroot use FatesAllometryMod , only : bsap_allom @@ -128,10 +130,10 @@ module EDCohortDynamicsMod !-------------------------------------------------------------------------------------! - subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & + + subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & prt, laimemory, status, recruitstatus,ctrim, & clayer, spread, bc_in) - ! ! !DESCRIPTION: ! create new cohort @@ -149,6 +151,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type), intent(inout), pointer :: patchptr + integer, intent(in) :: pft ! Cohort Plant Functional Type integer, intent(in) :: clayer ! canopy status of cohort ! (1 = canopy, 2 = understorey, etc.) @@ -159,6 +162,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & real(r8), intent(in) :: nn ! number of individuals in cohort ! per 'area' (10000m2 default) real(r8), intent(in) :: hite ! height: meters + real(r8), intent(in) :: coage ! cohort age in days real(r8), intent(in) :: dbh ! dbh: cm class(prt_vartypes),target :: prt ! The allocated PARTEH ! object @@ -169,6 +173,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & real(r8), intent(in) :: spread ! The community assembly effects how ! spread crowns are in horizontal space type(bc_in_type), intent(in) :: bc_in ! External boundary conditions + ! !LOCAL VARIABLES: type(ed_cohort_type), pointer :: new_cohort ! Pointer to New Cohort structure. @@ -180,7 +185,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & integer :: nlevsoi_hyd ! number of hydraulically active soil layers !---------------------------------------------------------------------- - + allocate(new_cohort) call nan_cohort(new_cohort) ! Make everything in the cohort not-a-number @@ -208,6 +213,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & new_cohort%n = nn new_cohort%hite = hite new_cohort%dbh = dbh + new_cohort%coage = coage new_cohort%canopy_trim = ctrim new_cohort%canopy_layer = clayer new_cohort%canopy_layer_yesterday = real(clayer, r8) @@ -218,10 +224,16 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, dbh, & ! leaf age fractions (which are defined by PARTEH) call UpdateCohortBioPhysRates(new_cohort) - call sizetype_class_index(new_cohort%dbh,new_cohort%pft, & + call sizetype_class_index(new_cohort%dbh, new_cohort%pft, & new_cohort%size_class,new_cohort%size_by_pft_class) - + ! If cohort age trackign is off we call this here once + ! just so everythin is in the first bin - + ! this makes it easier to copy and terminate cohorts later + ! we don't need to update this ever if cohort age tracking is off + call coagetype_class_index(new_cohort%coage, new_cohort%pft, & + new_cohort%coage_class,new_cohort%coage_by_pft_class) + ! This routine may be called during restarts, and at this point in the call sequence ! the actual cohort data is unknown, as this is really only used for allocation ! In these cases, testing if things like biomass are reasonable is pre-mature @@ -472,8 +484,12 @@ subroutine nan_cohort(cc_p) currentCohort%size_class = fates_unset_int ! size class index currentCohort%size_class_lasttimestep = fates_unset_int ! size class index currentCohort%size_by_pft_class = fates_unset_int ! size by pft classification index + currentCohort%coage_class = fates_unset_int ! cohort age class index + currentCohort%coage_by_pft_class = fates_unset_int ! cohort age by pft class index + currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) - currentCohort%dbh = nan ! 'diameter at breast height' in cm + currentCohort%dbh = nan ! 'diameter at breast height' in cm + currentCohort%coage = nan ! age of the cohort in years currentCohort%hite = nan ! height: meters currentCohort%laimemory = nan ! target leaf biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 @@ -583,6 +599,7 @@ subroutine zero_cohort(cc_p) currentcohort%ts_net_uptake(:) = 0._r8 currentcohort%fraction_crown_burned = 0._r8 currentCohort%size_class = 1 + currentCohort%coage_class = 1 currentCohort%seed_prod = 0._r8 currentCohort%size_class_lasttimestep = 0 currentcohort%npp_acc_hold = 0._r8 @@ -921,7 +938,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! Join similar cohorts to reduce total number ! ! !USES: - use EDParamsMod , only : ED_val_cohort_fusion_tol + use EDParamsMod , only : ED_val_cohort_size_fusion_tol + use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use FatesConstantsMod, only : days_per_year + ! ! !ARGUMENTS type (ed_site_type), intent(inout), target :: currentSite @@ -943,15 +963,19 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) integer :: nocohorts real(r8) :: newn real(r8) :: diff + real(r8) :: coage_diff real(r8) :: leaf_c_next ! Leaf carbon * plant density of current (for weighting) real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting) real(r8) :: leaf_c_target - real(r8) :: dynamic_fusion_tolerance + real(r8) :: dynamic_size_fusion_tolerance + real(r8) :: dynamic_age_fusion_tolerance real(r8) :: dbh real(r8) :: leaf_c ! leaf carbon [kg] integer :: largersc, smallersc, sc_i ! indices for tracking the growth flux caused by fusion real(r8) :: larger_n, smaller_n + integer :: oldercacls, youngercacls, cacls_i ! indices for tracking the age flux caused by fusion + real(r8) :: older_n, younger_n logical, parameter :: fuse_debug = .false. ! This debug is over-verbose ! and gets its own flag @@ -959,7 +983,9 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) !---------------------------------------------------------------------- !set initial fusion tolerance - dynamic_fusion_tolerance = ED_val_cohort_fusion_tol + dynamic_size_fusion_tolerance = ED_val_cohort_size_fusion_tol + ! set the cohort age fusion tolerance (this is in days - remains constant) + dynamic_age_fusion_tolerance = ED_val_cohort_age_fusion_tol !This needs to be a function of the canopy layer, because otherwise, at canopy closure !the number of cohorts doubles and very dissimilar cohorts are fused together @@ -994,7 +1020,20 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) !Criteria used to divide up the height continuum into different cohorts. - if (diff < dynamic_fusion_tolerance) then + if (diff < dynamic_size_fusion_tolerance) then + + ! Only fuse if the cohorts are within x years of each other + ! if they are the same age we make diff 0- to avoid errors divding by zero + !NB if cohort age tracking is off then the age of both should be 0 + ! and hence the age fusion criterion is met + if (currentCohort%coage .eq. nextc%coage) then + coage_diff = 0.0_r8 + else + coage_diff = abs((currentCohort%coage - nextc%coage)/ & + (0.5*(currentCohort%coage + nextc%coage))) + end if + + if (coage_diff <= dynamic_age_fusion_tolerance ) then ! Don't fuse a cohort with itself! if (.not.associated(currentCohort,nextc) ) then @@ -1025,6 +1064,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite + write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh write(fates_log(),*) 'pft:',currentCohort%pft,nextc%pft write(fates_log(),*) 'canopy_trim:',currentCohort%canopy_trim,nextc%canopy_trim @@ -1036,6 +1076,17 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) end do end if + ! new cohort age is weighted mean of two cohorts + currentCohort%coage = & + (currentCohort%coage * (currentCohort%n/(currentCohort%n + nextc%n))) + & + (nextc%coage * (nextc%n/(currentCohort%n + nextc%n))) + + ! update the cohort age again + if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + call coagetype_class_index(currentCohort%coage, currentCohort%pft, & + currentCohort%coage_class, currentCohort%coage_by_pft_class) + end if + ! Fuse all mass pools call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, & currentCohort%n/newn ) @@ -1167,7 +1218,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) - + + if(hlm_use_planthydro.eq.itrue) then call FuseCohortHydraulics(currentSite,currentCohort,nextc,bc_in,newn) endif @@ -1212,9 +1264,10 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! now that we've tracked the change flux. reset the memory of the prior timestep currentCohort%size_class_lasttimestep = currentCohort%size_class endif + ! Flux and biophysics variables have not been calculated for recruits we just default to - ! their initization values, which should be the same for eahc + ! their initization values, which should be the same for each if ( .not.currentCohort%isnew) then currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + & @@ -1245,6 +1298,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn + currentCohort%smort = (currentCohort%n*currentCohort%smort + nextc%n*nextc%smort)/newn + currentCohort%asmort = (currentCohort%n*currentCohort%asmort + nextc%n*nextc%asmort)/newn currentCohort%frmort = (currentCohort%n*currentCohort%frmort + nextc%n*nextc%frmort)/newn ! logging mortality, Yi Xu @@ -1316,7 +1371,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) endif !canopy layer endif !pft endif !index no. - endif !diff + endif ! cohort age diff + endif !diff nextc => nextnextc @@ -1350,8 +1406,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! !---------------------------------------------------------------------! - dynamic_fusion_tolerance = dynamic_fusion_tolerance * 1.1_r8 - + dynamic_size_fusion_tolerance = dynamic_size_fusion_tolerance * 1.1_r8 + dynamic_age_fusion_tolerance = dynamic_age_fusion_tolerance * 1.1_r8 !write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance else @@ -1359,13 +1415,14 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) iterate = 0 endif - if ( dynamic_fusion_tolerance .gt. 100._r8) then + if ( dynamic_size_fusion_tolerance .gt. 100._r8) then ! something has gone terribly wrong and we need to report what write(fates_log(),*) 'exceeded reasonable expectation of cohort fusion.' currentCohort => currentPatch%tallest nocohorts = 0 do while(associated(currentCohort)) - write(fates_log(),*) 'cohort ', nocohorts, currentCohort%dbh, currentCohort%canopy_layer, currentCohort%n + write(fates_log(),*) 'cohort ', nocohorts, currentCohort%dbh,& + currentCohort%coage, currentCohort%canopy_layer, currentCohort%n nocohorts = nocohorts + 1 currentCohort => currentCohort%shorter enddo @@ -1569,7 +1626,8 @@ subroutine copy_cohort( currentCohort,copyc ) ! VEGETATION STRUCTURE n%pft = o%pft n%n = o%n - n%dbh = o%dbh + n%dbh = o%dbh + n%coage = o%coage n%hite = o%hite n%laimemory = o%laimemory n%lai = o%lai @@ -1584,9 +1642,10 @@ subroutine copy_cohort( currentCohort,copyc ) n%excl_weight = o%excl_weight n%prom_weight = o%prom_weight n%size_class = o%size_class - n%size_class_lasttimestep = o%size_class_lasttimestep + n%size_class_lasttimestep = o%size_class_lasttimestep n%size_by_pft_class = o%size_by_pft_class - + n%coage_class = o%coage_class + n%coage_by_pft_class = o%coage_by_pft_class ! This transfers the PRT objects over. call n%prt%CopyPRTVartypes(o%prt) @@ -1638,6 +1697,8 @@ subroutine copy_cohort( currentCohort,copyc ) n%cmort = o%cmort n%bmort = o%bmort n%hmort = o%hmort + n%smort = o%smort + n%asmort = o%asmort n%frmort = o%frmort ! logging mortalities, Yi Xu @@ -1670,7 +1731,9 @@ subroutine copy_cohort( currentCohort,copyc ) n%size_class = o%size_class n%size_class_lasttimestep = o%size_class_lasttimestep n%size_by_pft_class = o%size_by_pft_class - + n%coage_class = o%coage_class + n%coage_by_pft_class = o%coage_by_pft_class + !Pointers n%taller => NULL() ! pointer to next tallest cohort n%shorter => NULL() ! pointer to next shorter cohort diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index eed631046a..0ea929c85a 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -19,6 +19,7 @@ module EDMortalityFunctionsMod use FatesInterfaceMod , only : hlm_use_planthydro use EDLoggingMortalityMod , only : LoggingMortality_frac use EDParamsMod , only : fates_mortality_disturbance_fraction + use EDParamsMod , only : ED_val_cohort_age_fusion_tol use FatesInterfaceMod , only : bc_in_type use PRTGenericMod, only : all_carbon_elements @@ -42,11 +43,11 @@ module EDMortalityFunctionsMod - subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort ) + subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmort ) ! ============================================================================ ! Calculate mortality rates from carbon storage, hydraulic cavitation, - ! background and freezing + ! background and freezing and size dependent senescence ! ============================================================================ use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm @@ -58,12 +59,18 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort ) real(r8),intent(out) :: cmort ! carbon starvation mortality real(r8),intent(out) :: hmort ! hydraulic failure mortality real(r8),intent(out) :: frmort ! freezing stress mortality + real(r8),intent(out) :: smort ! size dependent senescence term + real(r8),intent(out) :: asmort ! age dependent senescence term real(r8) :: frac ! relativised stored carbohydrate real(r8) :: leaf_c_target ! target leaf biomass kgC real(r8) :: store_c real(r8) :: hf_sm_threshold ! hydraulic failure soil moisture threshold real(r8) :: hf_flc_threshold ! hydraulic failure fractional loss of conductivity threshold + real(r8) :: mort_ip_size_senescence ! inflection point for increase in mortality with dbh + real(r8) :: mort_r_size_senescence ! rate of mortality increase with dbh in senesence term + real(r8) :: mort_ip_age_senescence ! inflection point for increase in mortality with age + real(r8) :: mort_r_age_senescence ! rate of mortality increase with age in senescence term real(r8) :: temp_dep_fraction ! Temp. function (freezing mortality) real(r8) :: temp_in_C ! Daily averaged temperature in Celcius real(r8) :: min_fmc_ag ! minimum fraction of maximum conductivity for aboveground @@ -75,12 +82,42 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort ) logical, parameter :: test_zero_mortality = .false. ! Developer test which ! may help to debug carbon imbalances ! and the like + + ! Size Dependent Senescence + ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh + mort_r_size_senescence = EDPftvarcon_inst%mort_r_size_senescence(cohort_in%pft) + mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) + + ! if the user has set the inflection point to be a huge number then + ! we zero smort (even though it would essentially be zero anyway) + if (mort_ip_size_senescence < 10000 ) then + smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & + (cohort_in%dbh - mort_ip_size_senescence) ) ) + else + smort = 0.0_r8 + end if + + ! if user has set cohort age fusion param > 0 we calculate age + ! dependent mortality + if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + ! Age Dependent Senescence + ! rate and inflection point define the change in mortality with age + mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) + mort_ip_age_senescence = EDPftvarcon_inst%mort_ip_age_senescence(cohort_in%pft) + asmort = 1.0_r8 / (1.0_r8 + exp(-1.0_r8 * mort_r_age_senescence * & + (cohort_in%coage - mort_ip_age_senescence ) ) ) + else + asmort = 0.0_r8 + end if + - if (hlm_use_ed_prescribed_phys .eq. ifalse) then + +if (hlm_use_ed_prescribed_phys .eq. ifalse) then ! 'Background' mortality (can vary as a function of ! density as in ED1.0 and ED2.0, but doesn't here for tractability) - bmort = EDPftvarcon_inst%bmort(cohort_in%pft) + + bmort = EDPftvarcon_inst%bmort(cohort_in%pft) ! Proxy for hydraulic failure induced mortality. hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft) @@ -141,9 +178,10 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort ) !mortality_rates = bmort + hmort + cmort - else ! i.e. hlm_use_ed_prescribed_phys is true + else ! i.e. hlm_use_ed_prescribed_phys is true + if ( cohort_in%canopy_layer .eq. 1) then - bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft) + bmort = EDPftvarcon_inst%prescribed_mortality_canopy(cohort_in%pft) else bmort = EDPftvarcon_inst%prescribed_mortality_understory(cohort_in%pft) endif @@ -157,6 +195,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort ) hmort = 0.0_r8 frmort = 0.0_r8 bmort = 0.0_r8 + smort = 0.0_r8 + asmort = 0.0_r8 end if return @@ -174,6 +214,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) ! ! !USES: + use FatesInterfaceMod, only : hlm_freq_day ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite @@ -185,6 +226,8 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) real(r8) :: bmort ! background mortality rate (fraction per year) real(r8) :: hmort ! hydraulic failure mortality rate (fraction per year) real(r8) :: frmort ! freezing mortality rate (fraction per year) + real(r8) :: smort ! size dependent senescence mortality rate (fraction per year) + real(r8) :: asmort ! age dependent senescence mortality rate (fraction per year) real(r8) :: dndt_logging ! Mortality rate (per day) associated with the a logging event integer :: ipft ! local copy of the pft index !---------------------------------------------------------------------- @@ -193,7 +236,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) ! Mortality for trees in the understorey. !if trees are in the canopy, then their death is 'disturbance'. This probably needs a different terminology - call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort) + call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort, asmort) call LoggingMortality_frac(ipft, currentCohort%dbh, currentCohort%canopy_layer, & currentCohort%lmort_direct, & currentCohort%lmort_collateral, & @@ -206,14 +249,36 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) if (currentCohort%canopy_layer > 1)then ! Include understory logging mortality rates not associated with disturbance dndt_logging = (currentCohort%lmort_direct + & - currentCohort%lmort_collateral + & - currentCohort%lmort_infra)/hlm_freq_day - currentCohort%dndt = -1.0_r8 * (cmort+hmort+bmort+frmort+dndt_logging) * currentCohort%n + currentCohort%lmort_collateral + & + currentCohort%lmort_infra)/hlm_freq_day + + ! this caps bmort so that daily mortality cannot exceed 0.995 + if ((cmort+hmort+bmort+frmort+smort+asmort + dndt_logging)*hlm_freq_day & + > 0.995_r8)then + bmort = (0.995_r8 - ((cmort+hmort+asmort+frmort+smort+dndt_logging)*hlm_freq_day)) & + /hlm_freq_day + endif + + currentCohort%dndt = -1.0_r8 * & + (cmort+hmort+bmort+frmort+smort+asmort + dndt_logging) & + * currentCohort%n else + + ! cap on bmort for canopy layers - bmort is adjusted so that total mortality + ! including disturbance fraction does not exceed 0.995 + if(((cmort+hmort+bmort+frmort+smort+asmort)*hlm_freq_day) > & + 0.995_r8 - fates_mortality_disturbance_fraction)then + bmort = (0.995_r8 - fates_mortality_disturbance_fraction - & + ((cmort+hmort+asmort+frmort+smort)*hlm_freq_day))/hlm_freq_day + endif + + ! Mortality from logging in the canopy is ONLY disturbance generating, don't ! update number densities via non-disturbance inducing death - currentCohort%dndt = -(1.0_r8 - fates_mortality_disturbance_fraction) & - * (cmort+hmort+bmort+frmort) * currentCohort%n + currentCohort%dndt= -(1.0_r8-fates_mortality_disturbance_fraction) & + * (cmort+hmort+bmort+frmort+smort+asmort+dndt_logging) * & + currentCohort%n + endif return diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 1f6501aa44..de4c64fc98 100644 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -158,6 +158,8 @@ subroutine disturbance_rates( site_in, bc_in) real(r8) :: bmort real(r8) :: hmort real(r8) :: frmort + real(r8) :: smort + real(r8) :: asmort real(r8) :: lmort_direct real(r8) :: lmort_collateral @@ -181,8 +183,8 @@ subroutine disturbance_rates( site_in, bc_in) ! Mortality for trees in the understorey. currentCohort%patchptr => currentPatch - call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort) - currentCohort%dmort = cmort+hmort+bmort+frmort + call mortality_rates(currentCohort,bc_in,cmort,hmort,bmort,frmort,smort,asmort) + currentCohort%dmort = cmort+hmort+bmort+frmort+smort+asmort call carea_allom(currentCohort%dbh,currentCohort%n,site_in%spread,currentCohort%pft, & currentCohort%c_area) @@ -191,6 +193,8 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%bmort = bmort currentCohort%hmort = hmort currentCohort%frmort = frmort + currentCohort%smort = smort + currentCohort%asmort = asmort call LoggingMortality_frac(currentCohort%pft, currentCohort%dbh, currentCohort%canopy_layer, & lmort_direct,lmort_collateral,lmort_infra,l_degrad ) @@ -287,6 +291,8 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction) end if currentCohort => currentCohort%taller enddo !currentCohort @@ -307,6 +313,8 @@ subroutine disturbance_rates( site_in, bc_in) currentCohort%bmort = currentCohort%bmort*(1.0_r8 - fates_mortality_disturbance_fraction) currentCohort%dmort = currentCohort%dmort*(1.0_r8 - fates_mortality_disturbance_fraction) currentCohort%frmort = currentCohort%frmort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%smort = currentCohort%smort*(1.0_r8 - fates_mortality_disturbance_fraction) + currentCohort%asmort = currentCohort%asmort*(1.0_r8 - fates_mortality_disturbance_fraction) currentCohort%lmort_direct = 0.0_r8 currentCohort%lmort_collateral = 0.0_r8 currentCohort%lmort_infra = 0.0_r8 @@ -621,6 +629,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = nan nc%bmort = nan nc%frmort = nan + nc%smort = nan + nc%asmort = nan nc%lmort_direct = nan nc%lmort_collateral = nan nc%lmort_infra = nan @@ -673,6 +683,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral @@ -697,6 +709,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral @@ -752,6 +766,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral @@ -825,6 +841,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort ! since these are the ones that weren't logged, @@ -884,6 +902,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral @@ -904,6 +924,8 @@ subroutine spawn_patches( currentSite, bc_in) nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort nc%frmort = currentCohort%frmort + nc%smort = currentCohort%smort + nc%asmort = currentCohort%asmort nc%dmort = currentCohort%dmort nc%lmort_direct = currentCohort%lmort_direct nc%lmort_collateral = currentCohort%lmort_collateral @@ -1958,7 +1980,7 @@ subroutine zero_patch(cp_p) ! FIRE - currentPatch%litter_moisture(:) = 0.0_r8 ! litter moisture + currentPatch%litter_moisture(:) = 0.0_r8 ! litter moisture currentPatch%fuel_eff_moist = 0.0_r8 ! average fuel moisture content of the ground fuel ! (incl. live grasses. omits 1000hr fuels) currentPatch%livegrass = 0.0_r8 ! total ag grass biomass in patch. 1=c3 grass, 2=c4 grass. gc/m2 diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 971eb48e9e..815c13875c 100644 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -1361,6 +1361,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) temp_cohort%canopy_trim = 0.8_r8 !starting with the canopy not fully expanded temp_cohort%pft = ft temp_cohort%hite = EDPftvarcon_inst%hgt_min(ft) + temp_cohort%coage = 0.0_r8 call h2d_allom(temp_cohort%hite,ft,temp_cohort%dbh) ! Initialize live pools @@ -1460,7 +1461,6 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! Initialize the PARTEH object, and determine the initial masses of all ! organs and elements. ! ----------------------------------------------------------------------------- - prt => null() call InitPRTObject(prt) @@ -1560,7 +1560,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in ) ! This initializes the cohort call create_cohort(currentSite,currentPatch, temp_cohort%pft, temp_cohort%n, & - temp_cohort%hite, temp_cohort%dbh, prt, & + temp_cohort%hite, temp_cohort%coage, temp_cohort%dbh, prt, & temp_cohort%laimemory, cohortstatus, recruitstatus, & temp_cohort%canopy_trim, currentPatch%NCL_p, currentSite%spread, bc_in) diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 9613d79957..285c5cae85 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -41,6 +41,7 @@ module EDInitMod use FatesInterfaceMod , only : numpft use FatesInterfaceMod , only : nleafage use FatesInterfaceMod , only : nlevsclass + use FatesInterfaceMod , only : nlevcoage use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : bagw_allom use FatesAllometryMod , only : bbgw_allom @@ -201,6 +202,7 @@ subroutine zero_site( site_in ) ! fusoin-induced growth flux of individuals site_in%growthflux_fusion(:,:) = 0._r8 +! site_in%ageflux_fusion(:,:) = 0._r8 ! demotion/promotion info site_in%demotion_rate(:) = 0._r8 @@ -468,6 +470,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) temp_cohort%pft = pft temp_cohort%n = EDPftvarcon_inst%initd(pft) * patch_in%area temp_cohort%hite = EDPftvarcon_inst%hgt_min(pft) + ! Calculate the plant diameter from height call h2d_allom(temp_cohort%hite,pft,temp_cohort%dbh) @@ -514,6 +517,9 @@ subroutine init_cohorts( site_in, patch_in, bc_in) if ( debug ) write(fates_log(),*) 'EDInitMod.F90 call create_cohort ' + temp_cohort%coage = 0.0_r8 + + ! -------------------------------------------------------------------------------- ! Initialize the mass of every element in every organ of the organ ! -------------------------------------------------------------------------------- @@ -581,7 +587,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in) call prt_obj%CheckInitialConditions() call create_cohort(site_in, patch_in, pft, temp_cohort%n, temp_cohort%hite, & - temp_cohort%dbh, prt_obj, temp_cohort%laimemory, cstatus, rstatus, & + temp_cohort%coage, temp_cohort%dbh, prt_obj, temp_cohort%laimemory, cstatus, rstatus, & temp_cohort%canopy_trim, 1, site_in%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index b11de373ed..a70932029d 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -1,3 +1,4 @@ + module EDMainMod ! =========================================================================== @@ -40,8 +41,10 @@ module EDMainMod use EDCohortDynamicsMod , only : UpdateCohortBioPhysRates use SFMainMod , only : fire_model use FatesSizeAgeTypeIndicesMod, only : get_age_class_index + use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use FatesLitterMod , only : litter_type use FatesLitterMod , only : ncwd + use EDtypesMod , only : ed_site_type use EDtypesMod , only : ed_patch_type use EDtypesMod , only : ed_cohort_type @@ -275,12 +278,14 @@ end subroutine ed_ecosystem_dynamics !-------------------------------------------------------------------------------! subroutine ed_integrate_state_variables(currentSite, bc_in ) ! + ! !DESCRIPTION: ! FIX(SPM,032414) refactor so everything goes through interface ! ! !USES: - ! + use EDParamsMod, only : ED_val_cohort_age_fusion_tol ! !ARGUMENTS: + type(ed_site_type) , intent(inout) :: currentSite type(bc_in_type) , intent(in) :: bc_in @@ -300,6 +305,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) real(r8) :: delta_dbh ! correction for dbh real(r8) :: delta_hite ! correction for hite + real(r8) :: current_npp ! place holder for calculating npp each year in prescribed physiology mode !----------------------------------------------------------------------- ! Set a pointer to this sites carbon12 mass balance @@ -355,13 +361,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) if (currentCohort%canopy_layer .eq. 1) then currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_canopy(ft) & * currentCohort%c_area / currentCohort%n - currentCohort%npp_acc = currentCohort%npp_acc_hold / hlm_days_per_year - ! for mass balancing currentCohort%gpp_acc = currentCohort%npp_acc currentCohort%resp_acc = 0._r8 - else currentCohort%npp_acc_hold = EDPftvarcon_inst%prescribed_npp_understory(ft) & * currentCohort%c_area / currentCohort%n @@ -377,7 +380,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%gpp_acc_hold = currentCohort%gpp_acc * real(hlm_days_per_year,r8) currentCohort%resp_acc_hold = currentCohort%resp_acc * real(hlm_days_per_year,r8) endif - + ! Conduct Maintenance Turnover (parteh) if(debug) call currentCohort%prt%CheckMassConservation(ft,3) if(any(currentSite%dstatus == [phen_dstat_moiston,phen_dstat_timeon])) then @@ -437,7 +440,21 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in) call updateSizeDepTreeHydStates(currentSite,currentCohort) end if + + ! if we are in age-dependent mortality mode + if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + ! update cohort age + currentCohort%coage = currentCohort%coage + hlm_freq_day + if(currentCohort%coage < 0.0_r8)then + write(fates_log(),*) 'negative cohort age?',currentCohort%coage + end if + + ! update cohort age class and age x pft class + call coagetype_class_index(currentCohort%coage, currentCohort%pft, & + currentCohort%coage_class,currentCohort%coage_by_pft_class) + end if + currentCohort => currentCohort%taller end do @@ -755,6 +772,8 @@ subroutine bypass_dynamics(currentSite) currentCohort%hmort = 0.0_r8 currentCohort%cmort = 0.0_r8 currentCohort%frmort = 0.0_r8 + currentCohort%smort = 0.0_r8 + currentCohort%asmort = 0.0_r8 currentCohort%dndt = 0.0_r8 currentCohort%dhdt = 0.0_r8 @@ -767,5 +786,8 @@ subroutine bypass_dynamics(currentSite) end subroutine bypass_dynamics - end module EDMainMod + + + + diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 986151a8fe..7c186dde77 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -40,7 +40,8 @@ module EDParamsMod real(r8),protected, public :: ED_val_phen_mindayson real(r8),protected, public :: ED_val_phen_ncolddayslim real(r8),protected, public :: ED_val_phen_coldtemp - real(r8),protected, public :: ED_val_cohort_fusion_tol + real(r8),protected, public :: ED_val_cohort_size_fusion_tol + real(r8),protected, public :: ED_val_cohort_age_fusion_tol real(r8),protected, public :: ED_val_patch_fusion_tol real(r8),protected, public :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry @@ -58,7 +59,7 @@ module EDParamsMod real(r8),protected,allocatable,public :: ED_val_history_sizeclass_bin_edges(:) real(r8),protected,allocatable,public :: ED_val_history_ageclass_bin_edges(:) real(r8),protected,allocatable,public :: ED_val_history_height_bin_edges(:) - + real(r8),protected,allocatable,public :: ED_val_history_coageclass_bin_edges(:) character(len=param_string_length),parameter,public :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac" character(len=param_string_length),parameter,public :: ED_name_comp_excln = "fates_comp_excln" @@ -79,7 +80,8 @@ module EDParamsMod character(len=param_string_length),parameter,public :: ED_name_phen_mindayson= "fates_phen_mindayson" character(len=param_string_length),parameter,public :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim" character(len=param_string_length),parameter,public :: ED_name_phen_coldtemp= "fates_phen_coldtemp" - character(len=param_string_length),parameter,public :: ED_name_cohort_fusion_tol= "fates_cohort_size_fusion_tol" + character(len=param_string_length),parameter,public :: ED_name_cohort_size_fusion_tol= "fates_cohort_size_fusion_tol" + character(len=param_string_length),parameter,public :: ED_name_cohort_age_fusion_tol = "fates_cohort_age_fusion_tol" character(len=param_string_length),parameter,public :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" character(len=param_string_length),parameter,public :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh" @@ -94,6 +96,7 @@ module EDParamsMod character(len=param_string_length),parameter,public :: ED_name_history_sizeclass_bin_edges= "fates_history_sizeclass_bin_edges" character(len=param_string_length),parameter,public :: ED_name_history_ageclass_bin_edges= "fates_history_ageclass_bin_edges" character(len=param_string_length),parameter,public :: ED_name_history_height_bin_edges= "fates_history_height_bin_edges" + character(len=param_string_length),parameter,public :: ED_name_history_coageclass_bin_edges = "fates_history_coageclass_bin_edges" ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE) ! ---------------------------------------------------------------------------------------------- @@ -180,7 +183,8 @@ subroutine FatesParamsInit() ED_val_phen_mindayson = nan ED_val_phen_ncolddayslim = nan ED_val_phen_coldtemp = nan - ED_val_cohort_fusion_tol = nan + ED_val_cohort_size_fusion_tol = nan + ED_val_cohort_age_fusion_tol = nan ED_val_patch_fusion_tol = nan ED_val_canopy_closure_thresh = nan hydr_kmax_rsurf1 = nan @@ -211,7 +215,10 @@ subroutine FatesRegisterParams(fates_params) use FatesParametersInterface, only : fates_parameters_type, dimension_name_scalar, dimension_shape_1d use FatesParametersInterface, only : dimension_name_history_size_bins, dimension_name_history_age_bins - use FatesParametersInterface, only : dimension_name_history_height_bins, dimension_shape_scalar + use FatesParametersInterface, only : dimension_name_history_height_bins + use FatesParametersInterface, only : dimension_name_history_coage_bins + use FatesParametersInterface, only : dimension_shape_scalar + implicit none @@ -221,6 +228,7 @@ subroutine FatesRegisterParams(fates_params) character(len=param_string_length), parameter :: dim_names_sizeclass(1) = (/dimension_name_history_size_bins/) character(len=param_string_length), parameter :: dim_names_ageclass(1) = (/dimension_name_history_age_bins/) character(len=param_string_length), parameter :: dim_names_height(1) = (/dimension_name_history_height_bins/) + character(len=param_string_length), parameter :: dim_names_coageclass(1) = (/dimension_name_history_coage_bins/) call FatesParamsInit() @@ -281,7 +289,10 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_phen_coldtemp, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - call fates_params%RegisterParameter(name=ED_name_cohort_fusion_tol, dimension_shape=dimension_shape_scalar, & + call fates_params%RegisterParameter(name=ED_name_cohort_size_fusion_tol, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + + call fates_params%RegisterParameter(name=ED_name_cohort_age_fusion_tol, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) call fates_params%RegisterParameter(name=ED_name_patch_fusion_tol, dimension_shape=dimension_shape_scalar, & @@ -295,7 +306,7 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=hydr_name_kmax_rsurf2, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - + call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -350,6 +361,10 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + + call fates_params%RegisterParameter(name=ED_name_history_coageclass_bin_edges, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names_coageclass) + end subroutine FatesRegisterParams @@ -423,8 +438,11 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_phen_coldtemp, & data=ED_val_phen_coldtemp) - call fates_params%RetreiveParameter(name=ED_name_cohort_fusion_tol, & - data=ED_val_cohort_fusion_tol) + call fates_params%RetreiveParameter(name=ED_name_cohort_size_fusion_tol, & + data=ED_val_cohort_size_fusion_tol) + + call fates_params%RetreiveParameter(name=ED_name_cohort_age_fusion_tol, & + data=ED_val_cohort_age_fusion_tol) call fates_params%RetreiveParameter(name=ED_name_patch_fusion_tol, & data=ED_val_patch_fusion_tol) @@ -493,7 +511,10 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameterAllocate(name=ED_name_history_height_bin_edges, & data=ED_val_history_height_bin_edges) - + + call fates_params%RetreiveParameterAllocate(name=ED_name_history_coageclass_bin_edges, & + data=ED_val_history_coageclass_bin_edges) + end subroutine FatesReceiveParams @@ -528,7 +549,8 @@ subroutine FatesReportParams(is_master) write(fates_log(),fmt0) 'ED_val_phen_mindayson = ',ED_val_phen_mindayson write(fates_log(),fmt0) 'ED_val_phen_ncolddayslim = ',ED_val_phen_ncolddayslim write(fates_log(),fmt0) 'ED_val_phen_coldtemp = ',ED_val_phen_coldtemp - write(fates_log(),fmt0) 'ED_val_cohort_fusion_tol = ',ED_val_cohort_fusion_tol + write(fates_log(),fmt0) 'ED_val_cohort_size_fusion_tol = ',ED_val_cohort_size_fusion_tol + write(fates_log(),fmt0) 'ED_val_cohort_age_fusion_tol = ',ED_val_cohort_age_fusion_tol write(fates_log(),fmt0) 'ED_val_patch_fusion_tol = ',ED_val_patch_fusion_tol write(fates_log(),fmt0) 'ED_val_canopy_closure_thresh = ',ED_val_canopy_closure_thresh write(fates_log(),fmt0) 'hydr_kmax_rsurf1 = ',hydr_kmax_rsurf1 diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 7a2fe0ef3a..8c1e67cc1c 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -94,6 +94,10 @@ module EDPftvarcon real(r8), allocatable :: maintresp_reduction_intercept(:) ! intercept of MR reduction as f(carbon storage), ! 0=no throttling, 1=max throttling real(r8), allocatable :: bmort(:) + real(r8), allocatable :: mort_ip_size_senescence(:) ! inflection point of dbh dependent senescence + real(r8), allocatable :: mort_r_size_senescence(:) ! rate of change in mortality with dbh + real(r8), allocatable :: mort_ip_age_senescence(:) ! inflection point of age dependent senescence + real(r8), allocatable :: mort_r_age_senescence(:) ! rate of change in mortality with age real(r8), allocatable :: mort_scalar_coldstress(:) real(r8), allocatable :: mort_scalar_cstarvation(:) real(r8), allocatable :: mort_scalar_hydrfailure(:) @@ -167,8 +171,8 @@ module EDPftvarcon ! Prescribed Physiology Mode Parameters real(r8), allocatable :: prescribed_npp_canopy(:) ! this is only for the special ! prescribed_physiology_mode - real(r8), allocatable :: prescribed_npp_understory(:) ! this is only for the special - ! prescribed_physiology_mode + real(r8), allocatable :: prescribed_npp_understory(:) ! this is only for the special + ! prescribed physiology mode real(r8), allocatable :: prescribed_mortality_canopy(:) ! this is only for the special ! prescribed_physiology_mode real(r8), allocatable :: prescribed_mortality_understory(:) ! this is only for the special @@ -709,6 +713,22 @@ subroutine Register_PFT(this, fates_params) call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_mort_r_size_senescence' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_mort_ip_size_senescence' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_mort_r_age_senescence' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_mort_ip_age_senescence' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_mort_scalar_coldstress' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -1209,6 +1229,22 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%bmort) + name = 'fates_mort_ip_size_senescence' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%mort_ip_size_senescence) + + name = 'fates_mort_r_size_senescence' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%mort_r_size_senescence) + + name = 'fates_mort_ip_age_senescence' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%mort_ip_age_senescence) + + name = 'fates_mort_r_age_senescence' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%mort_r_age_senescence) + name = 'fates_mort_scalar_coldstress' call fates_params%RetreiveParameterAllocate(name=name, & data=this%mort_scalar_coldstress) @@ -1938,6 +1974,10 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'grperc = ',EDPftvarcon_inst%grperc write(fates_log(),fmt0) 'c2b = ',EDPftvarcon_inst%c2b write(fates_log(),fmt0) 'bmort = ',EDPftvarcon_inst%bmort + write(fates_log(),fmt0) 'mort_ip_size_senescence = ', EDPftvarcon_inst%mort_ip_size_senescence + write(fates_log(),fmt0) 'mort_r_size_senescence = ', EDPftvarcon_inst%mort_r_size_senescence + write(fates_log(),fmt0) 'mort_ip_age_senescence = ', EDPftvarcon_inst%mort_ip_age_senescence + write(fates_log(),fmt0) 'mort_r_age_senescence = ', EDPftvarcon_inst%mort_r_age_senescence write(fates_log(),fmt0) 'mort_scalar_coldstress = ',EDPftvarcon_inst%mort_scalar_coldstress write(fates_log(),fmt0) 'mort_scalar_cstarvation = ',EDPftvarcon_inst%mort_scalar_cstarvation write(fates_log(),fmt0) 'mort_scalar_hydrfailure = ',EDPftvarcon_inst%mort_scalar_hydrfailure diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 6ed7598d63..f7508b3984 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -14,7 +14,8 @@ module EDTypesMod use FatesLitterMod, only : litter_type use FatesLitterMod, only : ncwd use FatesConstantsMod, only : n_anthro_disturbance_categories - + use FatesConstantsMod, only : days_per_year + implicit none private ! By default everything is private save @@ -22,7 +23,7 @@ module EDTypesMod integer, parameter, public :: maxPatchesPerSite = 14 ! maximum number of patches to live on a site integer, parameter, public :: maxPatchesPerSite_by_disttype(n_anthro_disturbance_categories) = & (/ 10, 4 /) !!! MUST SUM TO maxPatchesPerSite !!! - integer, parameter, public :: maxCohortsPerPatch = 100 ! maximum number of cohorts per patch + integer, parameter, public :: maxCohortsPerPatch = 300 ! maximum number of cohorts per patch integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers integer, parameter, public :: ican_upper = 1 ! Nominal index for the upper canopy @@ -39,6 +40,7 @@ module EDTypesMod integer, parameter, public :: max_nleafage = 4 ! This is the maximum number of leaf age pools, ! used for allocating scratch space + ! ------------------------------------------------------------------------------------- ! Radiation parameters ! These should be part of the radiation module, but since we only have one option @@ -213,6 +215,7 @@ module EDTypesMod integer :: pft ! pft number real(r8) :: n ! number of individuals in cohort per 'area' (10000m2 default) real(r8) :: dbh ! dbh: cm + real(r8) :: coage ! cohort age in years real(r8) :: hite ! height: meters integer :: indexnumber ! unique number for each cohort. (within clump?) real(r8) :: laimemory ! target leaf biomass- set from previous year: kGC per indiv @@ -239,12 +242,14 @@ module EDTypesMod ! this is used for history output. We maintain this in the main cohort memory ! because we don't want to continually re-calculate the cohort's position when ! performing size diagnostics at high-frequency calls + integer :: coage_class ! An index that indicates which age bin the cohort currently resides in + ! used for history output. integer :: size_by_pft_class ! An index that indicates the cohorts position of the joint size-class x functional ! type classification. We also maintain this in the main cohort memory ! because we don't want to continually re-calculate the cohort's position when ! performing size diagnostics at high-frequency calls - integer :: size_class_lasttimestep ! size class of the cohort at the end of the previous timestep (used for calculating growth flux) - + integer :: coage_by_pft_class ! An index that indicates the cohorts position of the join cohort age class x PFT + integer :: size_class_lasttimestep ! size class of the cohort at the last time step ! CARBON FLUXES @@ -316,7 +321,9 @@ module EDTypesMod real(r8) :: cmort ! carbon starvation mortality rate n/year real(r8) :: hmort ! hydraulic failure mortality rate n/year real(r8) :: frmort ! freezing mortality n/year - + real(r8) :: smort ! senesence mortality n/year + real(r8) :: asmort ! age senescence mortality n/year + ! Logging Mortality Rate ! Yi Xu & M. Huang real(r8) :: lmort_direct ! directly logging rate fraction /per logging activity @@ -747,9 +754,9 @@ module EDTypesMod real(r8), allocatable :: fmort_rate_crown(:,:) ! rate of individuals killed due to fire mortality ! from crown damage per year. on size x pft array - real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin - ! due to fusion in a given day. on size x pft array - + real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin + ! due to fusion in a given day. on size x pft array +! real(r8), allocatable :: ageflux_fusion(:,:) ! rate of individuals movign into a given age class bin due to fusion- age x pft @@ -958,6 +965,7 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%n = ', ccohort%n write(fates_log(),*) 'co%dbh = ', ccohort%dbh write(fates_log(),*) 'co%hite = ', ccohort%hite + write(fates_log(),*) 'co%coage = ', ccohort%coage write(fates_log(),*) 'co%laimemory = ', ccohort%laimemory write(fates_log(),*) 'leaf carbon = ', ccohort%prt%GetState(leaf_organ,all_carbon_elements) @@ -980,6 +988,8 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%prom_weight = ', ccohort%prom_weight write(fates_log(),*) 'co%size_class = ', ccohort%size_class write(fates_log(),*) 'co%size_by_pft_class = ', ccohort%size_by_pft_class + write(fates_log(),*) 'co%coage_class = ', ccohort%coage_class + write(fates_log(),*) 'co%coage_by_pft_class = ', ccohort%coage_by_pft_class write(fates_log(),*) 'co%gpp_acc_hold = ', ccohort%gpp_acc_hold write(fates_log(),*) 'co%gpp_acc = ', ccohort%gpp_acc write(fates_log(),*) 'co%gpp_tstep = ', ccohort%gpp_tstep @@ -1001,8 +1011,11 @@ subroutine dump_cohort(ccohort) write(fates_log(),*) 'co%c_area = ', ccohort%c_area write(fates_log(),*) 'co%cmort = ', ccohort%cmort write(fates_log(),*) 'co%bmort = ', ccohort%bmort + write(fates_log(),*) 'co%smort = ', ccohort%smort + write(fates_log(),*) 'co%asmort = ', ccohort%asmort write(fates_log(),*) 'co%hmort = ', ccohort%hmort write(fates_log(),*) 'co%frmort = ', ccohort%frmort + write(fates_log(),*) 'co%asmort = ', ccohort%asmort write(fates_log(),*) 'co%isnew = ', ccohort%isnew write(fates_log(),*) 'co%dndt = ', ccohort%dndt write(fates_log(),*) 'co%dhdt = ', ccohort%dhdt diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 3f8aaafc57..46d6aa8ea6 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -1,4 +1,5 @@ module FatesConstantsMod + ! This module is used to define global _immutable_ data. Everything in ! this module must have the parameter attribute. diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index b1870beee0..9fcb35023e 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -37,6 +37,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceMod , only : nlevheight use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_model_day + use FatesInterfaceMod , only : nlevcoage ! FIXME(bja, 2016-10) need to remove CLM dependancy use EDPftvarcon , only : EDPftvarcon_inst @@ -98,6 +99,7 @@ module FatesHistoryInterfaceMod ! ! For reference, some standardized abbreviations of the FATES dimensions are listed here: ! scls = size-class dimension + ! cacls = cohort age-class dimension ! pft = the pft dimension ! age = the age bin dimension ! height = the height bin dimension @@ -107,12 +109,14 @@ module FatesHistoryInterfaceMod ! we have create some "multiplexed" dimensions that combine two or more dimensions into a ! single dimension. Examples of these are the following: ! scpf = size class x PFT + ! cacpf = cohort age class x PFT ! cnlf = canopy layer x leaf layer ! cnlfpft = canopy layer x leaf layer x PFT ! scag = size class bin x age bin ! scagpft = size class bin x age bin x PFT ! agepft = age bin x PFT - ! + + ! A recipe for adding a new history variable to this module: ! (1) decide what time frequency it makes sense to update the variable at, and what dimension(s) ! you want to output the variable on @@ -292,6 +296,7 @@ module FatesHistoryInterfaceMod integer :: ih_growthflux_si_scpf integer :: ih_growthflux_fusion_si_scpf integer :: ih_ba_si_scpf + integer :: ih_agb_si_scpf integer :: ih_m1_si_scpf integer :: ih_m2_si_scpf integer :: ih_m3_si_scpf @@ -300,9 +305,13 @@ module FatesHistoryInterfaceMod integer :: ih_m6_si_scpf integer :: ih_m7_si_scpf integer :: ih_m8_si_scpf + integer :: ih_m9_si_scpf + integer :: ih_m10_si_scpf integer :: ih_crownfiremort_si_scpf integer :: ih_cambialfiremort_si_scpf + integer :: ih_m10_si_capf + integer :: ih_nplant_si_capf integer :: ih_ar_si_scpf integer :: ih_ar_grow_si_scpf @@ -346,6 +355,11 @@ module FatesHistoryInterfaceMod integer :: ih_m6_si_scls integer :: ih_m7_si_scls integer :: ih_m8_si_scls + integer :: ih_m9_si_scls + integer :: ih_m10_si_scls + + integer :: ih_m10_si_cacls + integer :: ih_nplant_si_cacls ! lots of non-default diagnostics for understanding canopy versus understory carbon balances integer :: ih_rdark_canopy_si_scls @@ -536,6 +550,7 @@ module FatesHistoryInterfaceMod type(iovar_map_type), pointer :: iovar_map(:) + !! THESE WERE EXPLICITLY PRIVATE WHEN TYPE WAS PUBLIC integer, private :: patch_index_, column_index_, levgrnd_index_, levscpf_index_ integer, private :: levscls_index_, levpft_index_, levage_index_ @@ -545,7 +560,9 @@ module FatesHistoryInterfaceMod integer, private :: levheight_index_ integer, private :: levelem_index_, levelpft_index_ integer, private :: levelcwd_index_, levelage_index_ + integer, private :: levcacls_index_, levcapf_index_ + contains procedure :: Init @@ -566,6 +583,8 @@ module FatesHistoryInterfaceMod procedure :: levgrnd_index procedure :: levscpf_index procedure :: levscls_index + procedure :: levcapf_index + procedure :: levcacls_index procedure :: levpft_index procedure :: levage_index procedure :: levfuel_index @@ -593,6 +612,8 @@ module FatesHistoryInterfaceMod procedure, private :: set_column_index procedure, private :: set_levgrnd_index procedure, private :: set_levscpf_index + procedure, private :: set_levcacls_index + procedure, private :: set_levcapf_index procedure, private :: set_levscls_index procedure, private :: set_levpft_index procedure, private :: set_levage_index @@ -625,6 +646,7 @@ subroutine Init(this, num_threads, fates_bounds) use FatesIODimensionsMod, only : patch, column, levgrnd, levscpf use FatesIODimensionsMod, only : levscls, levpft, levage + use FatesIODimensionsMod, only : levcacls, levcapf use FatesIODimensionsMod, only : levfuel, levcwdsc, levscag use FatesIODimensionsMod, only : levscagpft, levagepft use FatesIODimensionsMod, only : levcan, levcnlf, levcnlfpft @@ -666,6 +688,16 @@ subroutine Init(this, num_threads, fates_bounds) call this%dim_bounds(dim_count)%Init(levscls, num_threads, & fates_bounds%size_class_begin, fates_bounds%size_class_end) + dim_count = dim_count + 1 + call this%set_levcacls_index(dim_count) + call this%dim_bounds(dim_count)%Init(levcacls, num_threads, & + fates_bounds%coage_class_begin, fates_bounds%coage_class_end) + + dim_count = dim_count + 1 + call this%set_levcapf_index(dim_count) + call this%dim_bounds(dim_count)%Init(levcapf, num_threads, & + fates_bounds%coagepf_class_begin, fates_bounds%coagepf_class_end) + dim_count = dim_count + 1 call this%set_levpft_index(dim_count) call this%dim_bounds(dim_count)%Init(levpft, num_threads, & @@ -783,6 +815,14 @@ subroutine SetThreadBoundsEach(this, thread_index, thread_bounds) call this%dim_bounds(index)%SetThreadBounds(thread_index, & thread_bounds%size_class_begin, thread_bounds%size_class_end) + index = this%levcacls_index() + call this%dim_bounds(index)%SetThreadBounds(thread_index, & + thread_bounds%coage_class_begin, thread_bounds%coage_class_end) + + index = this%levcapf_index() + call this%dim_bounds(index)%SetThreadBounds(thread_index, & + thread_bounds%coagepf_class_begin, thread_bounds%coagepf_class_end) + index = this%levpft_index() call this%dim_bounds(index)%SetThreadBounds(thread_index, & thread_bounds%pft_class_begin, thread_bounds%pft_class_end) @@ -856,6 +896,7 @@ subroutine assemble_history_output_types(this) use FatesIOVariableKindMod, only : patch_r8, patch_ground_r8, patch_size_pft_r8 use FatesIOVariableKindMod, only : site_r8, site_ground_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 + use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8 @@ -888,6 +929,12 @@ subroutine assemble_history_output_types(this) call this%set_dim_indices(site_size_r8, 1, this%column_index()) call this%set_dim_indices(site_size_r8, 2, this%levscls_index()) + call this%set_dim_indices(site_coage_r8, 1, this%column_index()) + call this%set_dim_indices(site_coage_r8, 2, this%levcacls_index()) + + call this%set_dim_indices(site_coage_pft_r8, 1, this%column_index()) + call this%set_dim_indices(site_coage_pft_r8, 2, this%levcapf_index()) + call this%set_dim_indices(site_pft_r8, 1, this%column_index()) call this%set_dim_indices(site_pft_r8, 2, this%levpft_index()) @@ -1047,6 +1094,34 @@ integer function levscls_index(this) levscls_index = this%levscls_index_ end function levscls_index +!========================================================================= + subroutine set_levcacls_index(this, index) + implicit none + class(fates_history_interface_type), intent(inout) :: this + integer, intent(in) :: index + this%levcacls_index_ = index +end subroutine set_levcacls_index + +integer function levcacls_index(this) + implicit none + class(fates_history_interface_type), intent(in) :: this + levcacls_index = this%levcacls_index_ +end function levcacls_index + +!========================================================================= + subroutine set_levcapf_index(this, index) + implicit none + class(fates_history_interface_type), intent(inout) :: this + integer, intent(in) :: index + this%levcapf_index_ = index + end subroutine set_levcapf_index + +integer function levcapf_index(this) + implicit none + class(fates_history_interface_type), intent(in) :: this + levcapf_index = this%levcapf_index_ +end function levcapf_index + ! ======================================================================= subroutine set_levpft_index(this, index) implicit none @@ -1352,6 +1427,7 @@ subroutine init_dim_kinds_maps(this) use FatesIOVariableKindMod, only : patch_r8, patch_ground_r8, patch_size_pft_r8 use FatesIOVariableKindMod, only : site_r8, site_ground_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 + use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8 @@ -1395,6 +1471,14 @@ subroutine init_dim_kinds_maps(this) index = index + 1 call this%dim_kinds(index)%Init(site_size_r8, 2) + ! site x cohort age-class/pft + index = index + 1 + call this%dim_kinds(index)%Init(site_coage_pft_r8, 2) + + ! site x cohort age-class + index = index + 1 + call this%dim_kinds(index)%Init(site_coage_r8, 2) + ! site x pft index = index + 1 call this%dim_kinds(index)%Init(site_pft_r8, 2) @@ -1533,6 +1617,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) use FatesSizeAgeTypeIndicesMod, only : get_age_class_index use FatesSizeAgeTypeIndicesMod, only : get_height_index use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index + use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDTypesMod , only : nlevleaf use EDParamsMod, only : ED_val_history_height_bin_edges @@ -1561,6 +1646,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) integer :: cwd integer :: elcwd, elpft ! combined index of element and pft or cwd integer :: i_scpf,i_pft,i_scls ! iterators for scpf, pft, and scls dims + integer :: i_cacls, i_capf ! iterators for cohort age and cohort age x pft integer :: i_cwd,i_fuel ! iterators for cwd and fuel dims integer :: iscag ! size-class x age index integer :: iscagpft ! size-class x age x pft index @@ -1578,6 +1664,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8) :: n_perm2 ! individuals per m2 for the whole column real(r8) :: patch_scaling_scalar ! ratio of canopy to patch area for counteracting patch scaling real(r8) :: dbh ! diameter ("at breast height") + real(r8) :: coage ! cohort age real(r8) :: npp_partition_error ! a check that the NPP partitions sum to carbon allocation real(r8) :: frac_canopy_in_bin ! fraction of a leaf's canopy that is within a given height bin real(r8) :: binbottom,bintop ! edges of height bins @@ -1612,6 +1699,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8), parameter :: tiny = 1.e-5_r8 ! some small number real(r8), parameter :: reallytalltrees = 1000. ! some large number (m) + integer :: tmp + associate( hio_npatches_si => this%hvars(ih_npatches_si)%r81d, & hio_ncohorts_si => this%hvars(ih_ncohorts_si)%r81d, & hio_trimming_pa => this%hvars(ih_trimming_pa)%r81d, & @@ -1695,8 +1784,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_growthflux_si_scpf => this%hvars(ih_growthflux_si_scpf)%r82d, & hio_growthflux_fusion_si_scpf => this%hvars(ih_growthflux_fusion_si_scpf)%r82d, & hio_ba_si_scpf => this%hvars(ih_ba_si_scpf)%r82d, & + hio_agb_si_scpf => this%hvars(ih_agb_si_scpf)%r82d, & hio_nplant_si_scpf => this%hvars(ih_nplant_si_scpf)%r82d, & - + hio_nplant_si_capf => this%hvars(ih_nplant_si_capf)%r82d, & + hio_m1_si_scpf => this%hvars(ih_m1_si_scpf)%r82d, & hio_m2_si_scpf => this%hvars(ih_m2_si_scpf)%r82d, & hio_m3_si_scpf => this%hvars(ih_m3_si_scpf)%r82d, & @@ -1705,6 +1796,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m6_si_scpf => this%hvars(ih_m6_si_scpf)%r82d, & hio_m7_si_scpf => this%hvars(ih_m7_si_scpf)%r82d, & hio_m8_si_scpf => this%hvars(ih_m8_si_scpf)%r82d, & + hio_m9_si_scpf => this%hvars(ih_m9_si_scpf)%r82d, & + hio_m10_si_scpf => this%hvars(ih_m10_si_scpf)%r82d, & + hio_m10_si_capf => this%hvars(ih_m10_si_capf)%r82d, & + hio_crownfiremort_si_scpf => this%hvars(ih_crownfiremort_si_scpf)%r82d, & hio_cambialfiremort_si_scpf => this%hvars(ih_cambialfiremort_si_scpf)%r82d, & @@ -1717,8 +1812,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m5_si_scls => this%hvars(ih_m5_si_scls)%r82d, & hio_m6_si_scls => this%hvars(ih_m6_si_scls)%r82d, & hio_m7_si_scls => this%hvars(ih_m7_si_scls)%r82d, & - hio_m8_si_scls => this%hvars(ih_m8_si_scls)%r82d, & - hio_c13disc_si_scpf => this%hvars(ih_c13disc_si_scpf)%r82d, & + hio_m8_si_scls => this%hvars(ih_m8_si_scls)%r82d, & + hio_m9_si_scls => this%hvars(ih_m9_si_scls)%r82d, & + hio_m10_si_scls => this%hvars(ih_m10_si_scls)%r82d, & + hio_m10_si_cacls => this%hvars(ih_m10_si_cacls)%r82d, & + + hio_c13disc_si_scpf => this%hvars(ih_c13disc_si_scpf)%r82d, & + hio_cwd_elcwd => this%hvars(ih_cwd_elcwd)%r82d, & hio_cwd_ag_elem => this%hvars(ih_cwd_ag_elem)%r82d, & hio_cwd_bg_elem => this%hvars(ih_cwd_bg_elem)%r82d, & @@ -1728,6 +1828,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_agb_si_scls => this%hvars(ih_agb_si_scls)%r82d, & hio_biomass_si_scls => this%hvars(ih_biomass_si_scls)%r82d, & hio_nplant_si_scls => this%hvars(ih_nplant_si_scls)%r82d, & + hio_nplant_si_cacls => this%hvars(ih_nplant_si_cacls)%r82d, & hio_nplant_canopy_si_scls => this%hvars(ih_nplant_canopy_si_scls)%r82d, & hio_nplant_understory_si_scls => this%hvars(ih_nplant_understory_si_scls)%r82d, & hio_lai_canopy_si_scls => this%hvars(ih_lai_canopy_si_scls)%r82d, & @@ -1842,7 +1943,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) io_si = this%iovar_map(nc)%site_index(s) io_pa1 = this%iovar_map(nc)%patch1_index(s) io_soipa = io_pa1-1 - + ! Set trimming on the soil patch to 1.0 hio_trimming_pa(io_soipa) = 1.0_r8 @@ -1947,7 +2048,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) ft = ccohort%pft call sizetype_class_index(ccohort%dbh, ccohort%pft, ccohort%size_class, ccohort%size_by_pft_class) - + call coagetype_class_index(ccohort%coage, ccohort%pft, & + ccohort%coage_class, ccohort%coage_by_pft_class) + ! Increment the number of cohorts per site hio_ncohorts_si(io_si) = hio_ncohorts_si(io_si) + 1._r8 @@ -2103,14 +2206,20 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_npp_stor_si(io_si) = hio_npp_stor_si(io_si) + store_c_net_alloc * n_perm2 associate( scpf => ccohort%size_by_pft_class, & - scls => ccohort%size_class ) - - gpp_cached = hio_gpp_si_scpf(io_si,scpf) + scls => ccohort%size_class, & + cacls => ccohort%coage_class, & + capf => ccohort%coage_by_pft_class) + + + gpp_cached = hio_gpp_si_scpf(io_si,scpf) + hio_gpp_si_scpf(io_si,scpf) = hio_gpp_si_scpf(io_si,scpf) + & n_perm2*ccohort%gpp_acc_hold ! [kgC/m2/yr] hio_npp_totl_si_scpf(io_si,scpf) = hio_npp_totl_si_scpf(io_si,scpf) + & ccohort%npp_acc_hold *n_perm2 + + hio_npp_leaf_si_scpf(io_si,scpf) = hio_npp_leaf_si_scpf(io_si,scpf) + & leaf_c_net_alloc*n_perm2 hio_npp_fnrt_si_scpf(io_si,scpf) = hio_npp_fnrt_si_scpf(io_si,scpf) + & @@ -2155,7 +2264,11 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m7_si_scpf(io_si,scpf) = hio_m7_si_scpf(io_si,scpf) + & (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n hio_m8_si_scpf(io_si,scpf) = hio_m8_si_scpf(io_si,scpf) + ccohort%frmort*ccohort%n - + hio_m9_si_scpf(io_si,scpf) = hio_m9_si_scpf(io_si,scpf) + ccohort%smort*ccohort%n + hio_m10_si_scpf(io_si,scpf) = hio_m10_si_scpf(io_si,scpf) + ccohort%asmort*ccohort%n + + hio_m10_si_capf(io_si,capf) = hio_m10_si_capf(io_si,capf) + ccohort%asmort*ccohort%n + hio_m1_si_scls(io_si,scls) = hio_m1_si_scls(io_si,scls) + ccohort%bmort*ccohort%n hio_m2_si_scls(io_si,scls) = hio_m2_si_scls(io_si,scls) + ccohort%hmort*ccohort%n hio_m3_si_scls(io_si,scls) = hio_m3_si_scls(io_si,scls) + ccohort%cmort*ccohort%n @@ -2163,7 +2276,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n hio_m8_si_scls(io_si,scls) = hio_m8_si_scls(io_si,scls) + & ccohort%frmort*ccohort%n - + hio_m9_si_scls(io_si,scls) = hio_m9_si_scls(io_si,scls) + ccohort%smort*ccohort%n + + hio_m10_si_scls(io_si,scls) = hio_m10_si_scls(io_si,scls) + ccohort%asmort*ccohort%n + + hio_m10_si_cacls(io_si,cacls) = hio_m10_si_cacls(io_si,cacls)+ & + ccohort%asmort*ccohort%n + !C13 discrimination if(gpp_cached + ccohort%gpp_acc_hold > 0.0_r8)then hio_c13disc_si_scpf(io_si,scpf) = ((hio_c13disc_si_scpf(io_si,scpf) * gpp_cached) + & @@ -2174,10 +2293,15 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! number density [/ha] hio_nplant_si_scpf(io_si,scpf) = hio_nplant_si_scpf(io_si,scpf) + ccohort%n - - + ! number density along the cohort age dimension + hio_nplant_si_capf(io_si,capf) = hio_nplant_si_capf(io_si,capf) + ccohort%n + + ! number density by size and biomass hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & - total_c * ccohort%n * EDPftvarcon_inst%allom_agb_frac(ccohort%pft) * AREA_INV + total_c * ccohort%n * EDPftvarcon_inst%allom_agb_frac(ccohort%pft) * AREA_INV + + hio_agb_si_scpf(io_si,scpf) = hio_agb_si_scpf(io_si,scpf) + & + total_c * ccohort%n * EDPftvarcon_inst%allom_agb_frac(ccohort%pft) * AREA_INV hio_biomass_si_scls(io_si,scls) = hio_biomass_si_scls(io_si,scls) + & total_c * ccohort%n * AREA_INV @@ -2189,7 +2313,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nplant_si_scag(io_si,iscag) = hio_nplant_si_scag(io_si,iscag) + ccohort%n hio_nplant_si_scls(io_si,scls) = hio_nplant_si_scls(io_si,scls) + ccohort%n - + hio_nplant_si_cacls(io_si,cacls) = hio_nplant_si_cacls(io_si,cacls)+ccohort%n + ! update size, age, and PFT - indexed quantities iscagpft = get_sizeagepft_class_index(ccohort%dbh,cpatch%age,ccohort%pft) @@ -2207,12 +2332,12 @@ subroutine update_history_dyn(this,nc,nsites,sites) total_c * ccohort%n * AREA_INV - ! update SCPF/SCLS- and canopy/subcanopy- partitioned quantities if (ccohort%canopy_layer .eq. 1) then hio_nplant_canopy_si_scag(io_si,iscag) = hio_nplant_canopy_si_scag(io_si,iscag) + ccohort%n hio_mortality_canopy_si_scag(io_si,iscag) = hio_mortality_canopy_si_scag(io_si,iscag) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n hio_ddbh_canopy_si_scag(io_si,iscag) = hio_ddbh_canopy_si_scag(io_si,iscag) + & ccohort%ddbhdt*ccohort%n hio_bstor_canopy_si_scpf(io_si,scpf) = hio_bstor_canopy_si_scpf(io_si,scpf) + & @@ -2223,11 +2348,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_canopy_biomass_pa(io_pa) = hio_canopy_biomass_pa(io_pa) + n_density * total_c * g_per_kg !hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & - ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ! ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n hio_mortality_canopy_si_scpf(io_si,scpf) = hio_mortality_canopy_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + & - (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & + + (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort + & + ccohort%smort + ccohort%asmort) * ccohort%n + & + (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_nplant_canopy_si_scpf(io_si,scpf) = hio_nplant_canopy_si_scpf(io_si,scpf) + ccohort%n @@ -2252,12 +2380,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! sum of all mortality hio_mortality_canopy_si_scls(io_si,scls) = hio_mortality_canopy_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_canopy_mortality_carbonflux_si(io_si) = hio_canopy_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * & total_c * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_c * & ccohort%n * g_per_kg * ha_per_m2 @@ -2299,7 +2429,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) else hio_nplant_understory_si_scag(io_si,iscag) = hio_nplant_understory_si_scag(io_si,iscag) + ccohort%n hio_mortality_understory_si_scag(io_si,iscag) = hio_mortality_understory_si_scag(io_si,iscag) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n hio_ddbh_understory_si_scag(io_si,iscag) = hio_ddbh_understory_si_scag(io_si,iscag) + & ccohort%ddbhdt*ccohort%n hio_bstor_understory_si_scpf(io_si,scpf) = hio_bstor_understory_si_scpf(io_si,scpf) + & @@ -2309,11 +2440,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_understory_biomass_pa(io_pa) = hio_understory_biomass_pa(io_pa) + & n_density * total_c * g_per_kg !hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & - ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * ccohort%n + ! (ccohort%bmort + ccohort%hmort + ccohort%cmort + + ! ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n hio_mortality_understory_si_scpf(io_si,scpf) = hio_mortality_understory_si_scpf(io_si,scpf)+ & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & - (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n + & + (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_nplant_understory_si_scpf(io_si,scpf) = hio_nplant_understory_si_scpf(io_si,scpf) + ccohort%n @@ -2339,12 +2472,14 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! sum of all mortality hio_mortality_understory_si_scls(io_si,scls) = hio_mortality_understory_si_scls(io_si,scls) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort ) * ccohort%n + & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * ccohort%n + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * & ccohort%n * sec_per_day * days_per_year hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & - (ccohort%bmort + ccohort%hmort + ccohort%cmort + ccohort%frmort) * & + (ccohort%bmort + ccohort%hmort + ccohort%cmort + & + ccohort%frmort + ccohort%smort + ccohort%asmort) * & total_c * ccohort%n * g_per_kg * days_per_sec * years_per_day * ha_per_m2 + & (ccohort%lmort_direct + ccohort%lmort_collateral + ccohort%lmort_infra) * total_c * & ccohort%n * g_per_kg * ha_per_m2 @@ -2396,6 +2531,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) end do end if ccohort%size_class_lasttimestep = scls + + ! end associate else ! i.e. cohort%isnew @@ -2404,7 +2541,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) i_scpf = (ccohort%pft-1)*nlevsclass+1 hio_growthflux_si_scpf(io_si,i_scpf) = hio_growthflux_si_scpf(io_si,i_scpf) + ccohort%n * days_per_year ccohort%size_class_lasttimestep = 1 - ! + end if ! resolve some canopy area profiles, both total and of occupied leaves @@ -2573,10 +2710,22 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! while in this loop, pass the fusion-induced growth rate flux to history hio_growthflux_fusion_si_scpf(io_si,i_scpf) = hio_growthflux_fusion_si_scpf(io_si,i_scpf) + & sites(s)%growthflux_fusion(i_scls, i_pft) * days_per_year - end do end do ! + + ! do i_pft = 1, numpft + ! do i_cacls = 1, nlevcoage +! i_capf = (i_pft -1)*nlevcoage + i_cacls +! fix this! + ! i_capf = 1 + ! hio_ageflux_fusion_si_capf(io_si,i_capf) = hio_ageflux_fusion_si_capf(io_si,i_capf) + & + ! sites(s)%ageflux_fusion(i_capf, i_pft) * days_per_year + ! end do + ! end do + + + ! treat carbon flux from imort the same way hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & sites(s)%imort_carbonflux @@ -2611,8 +2760,10 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_m4_si_scpf(io_si,i_scpf) + & hio_m5_si_scpf(io_si,i_scpf) + & hio_m6_si_scpf(io_si,i_scpf) + & - hio_m7_si_scpf(io_si,i_scpf) + & - hio_m8_si_scpf(io_si,i_scpf) + hio_m7_si_scpf(io_si,i_scpf) + & + hio_m8_si_scpf(io_si,i_scpf) + & + hio_m9_si_scpf(io_si,i_scpf) + & + hio_m10_si_scpf(io_si,i_scpf) end do end do @@ -2748,6 +2899,9 @@ subroutine update_history_dyn(this,nc,nsites,sites) end do + + + ! pass demotion rates and associated carbon fluxes to history do i_scls = 1,nlevsclass hio_demotion_rate_si_scls(io_si,i_scls) = sites(s)%demotion_rate(i_scls) * days_per_year @@ -2964,7 +3118,7 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) ! Calculate index for the scpf class associate( scpf => ccohort%size_by_pft_class, & scls => ccohort%size_class ) - + ! scale up cohort fluxes to their patches hio_npp_pa(io_pa) = hio_npp_pa(io_pa) + & ccohort%npp_tstep * g_per_kg * n_density * per_dt_tstep @@ -3310,6 +3464,7 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) hio_btran_scpf => this%hvars(ih_btran_scpf)%r82d, & hio_h2oveg_si => this%hvars(ih_h2oveg_si)%r81d, & hio_nplant_si_scpf => this%hvars(ih_nplant_si_scpf)%r82d, & + hio_nplant_si_capf => this%hvars(ih_nplant_si_capf)%r82d, & hio_h2oveg_hydro_err_si => this%hvars(ih_h2oveg_hydro_err_si)%r81d ) ! Flush the relevant history variables @@ -3641,6 +3796,7 @@ subroutine define_history_vars(this, initialize_variables) use FatesIOVariableKindMod, only : patch_r8, patch_ground_r8, patch_size_pft_r8 use FatesIOVariableKindMod, only : site_r8, site_ground_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 + use FatesIOVariableKindMod, only : site_coage_pft_r8, site_coage_r8 use FatesIOVariableKindMod, only : site_height_r8 use FatesInterfaceMod , only : hlm_use_planthydro @@ -4509,11 +4665,21 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_ba_si_scpf ) + call this%set_history_var(vname='AGB_SCPF', units = 'kgC/m2', & + long='Aboveground biomass by pft/size', use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_agb_si_scpf ) + call this%set_history_var(vname='NPLANT_SCPF', units = 'N/ha', & long='stem number density by pft/size', use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_scpf ) + call this%set_history_var(vname='NPLANT_CAPF', units = 'N/ha', & + long='stem number density by pft/coage', use_default='inactive', & + avgflag='A', vtype=site_coage_pft_r8, hlms='CLM:ALM',flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_capf ) + call this%set_history_var(vname='M1_SCPF', units = 'N/ha/yr', & long='background mortality by pft/size', use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -4564,6 +4730,21 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m8_si_scpf ) + call this%set_history_var(vname='M9_SCPF', units = 'N/ha/yr', & + long='senescence mortality by pft/size',use_default='inactive', & + avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m9_si_scpf ) + + call this%set_history_var(vname='M10_SCPF', units = 'N/ha/yr', & + long='age senescence mortality by pft/size',use_default='inactive', & + avgflag='A', vtype =site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m10_si_scpf ) + + call this%set_history_var(vname='M10_CAPF',units='N/ha/yr', & + long='age senescence mortality by pft/cohort age',use_default='inactive', & + avgflag='A', vtype =site_coage_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index =ih_m10_si_capf ) + call this%set_history_var(vname='MORTALITY_CANOPY_SCPF', units = 'N/ha/yr', & long='total mortality of canopy plants by pft/size', use_default='inactive', & avgflag='A', vtype=site_size_pft_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -4727,7 +4908,7 @@ subroutine define_history_vars(this, initialize_variables) long='number of canopy plants by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_canopy_si_scls ) - + call this%set_history_var(vname='LAI_CANOPY_SCLS', units = 'm2/m2', & long='Leaf are index (LAI) by size class', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -4763,6 +4944,11 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_scls ) + call this%set_history_var(vname='NPLANT_CACLS', units = 'indiv/ha', & + long='number of plants by coage class', use_default='active', & + avgflag='A', vtype=site_coage_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_nplant_si_cacls ) + call this%set_history_var(vname='M1_SCLS', units = 'N/ha/yr', & long='background mortality by size', use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & @@ -4802,7 +4988,22 @@ subroutine define_history_vars(this, initialize_variables) long='freezing mortality by size',use_default='active', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m8_si_scls ) - + + call this%set_history_var(vname='M9_SCLS', units = 'N/ha/yr', & + long='senescence mortality by size',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m9_si_scls ) + + call this%set_history_var(vname='M10_SCLS', units = 'N/ha/yr', & + long='age senescence mortality by size',use_default='active', & + avgflag='A', vtype=site_coage_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m10_si_scls ) + + call this%set_history_var(vname='M10_CACLS', units = 'N/ha/yr', & + long='age senescence mortality by cohort age',use_default='active', & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m10_si_cacls ) + call this%set_history_var(vname='CARBON_BALANCE_CANOPY_SCLS', units = 'kg C / ha / yr', & long='CARBON_BALANCE for canopy plants by size class', use_default='inactive', & avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & diff --git a/main/FatesHistoryVariableType.F90 b/main/FatesHistoryVariableType.F90 index 41451f0615..83bb9fdf1f 100644 --- a/main/FatesHistoryVariableType.F90 +++ b/main/FatesHistoryVariableType.F90 @@ -7,6 +7,7 @@ module FatesHistoryVariableType use FatesIOVariableKindMod, only : patch_r8, patch_ground_r8, patch_size_pft_r8 use FatesIOVariableKindMod, only : site_r8, site_ground_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 + use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 use FatesIOVariableKindMod, only : site_height_r8, patch_int use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 @@ -57,7 +58,6 @@ subroutine Init(this, vname, units, long, use_default, & vtype, avgflag, flushval, upfreq, num_dim_kinds, dim_kinds, dim_bounds) - implicit none class(fates_history_variable_type), intent(inout) :: this @@ -133,6 +133,14 @@ subroutine Init(this, vname, units, long, use_default, & allocate(this%r82d(lb1:ub1, lb2:ub2)) this%r82d(:,:) = flushval + case(site_coage_r8) + allocate(this%r82d(lb1:ub1, lb2:ub2)) + this%r82d(:,:) = flushval + + case(site_coage_pft_r8) + allocate(this%r82d(lb1:ub1, lb2:ub2)) + this%r82d(:,:) = flushval + case(site_pft_r8) allocate(this%r82d(lb1:ub1, lb2:ub2)) this%r82d(:,:) = flushval @@ -257,7 +265,6 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) end subroutine GetBounds subroutine Flush(this, thread, dim_bounds, dim_kinds) - implicit none @@ -285,6 +292,10 @@ subroutine Flush(this, thread, dim_bounds, dim_kinds) this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_size_r8) this%r82d(lb1:ub1, lb2:ub2) = this%flushval + case(site_coage_pft_r8) + this%r82d(lb1:ub1, lb2:ub2) = this%flushval + case(site_coage_r8) + this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_pft_r8) this%r82d(lb1:ub1, lb2:ub2) = this%flushval case(site_age_r8) diff --git a/main/FatesIODimensionsMod.F90 b/main/FatesIODimensionsMod.F90 index 152a951e63..0b3c72967d 100644 --- a/main/FatesIODimensionsMod.F90 +++ b/main/FatesIODimensionsMod.F90 @@ -7,7 +7,9 @@ module FatesIODimensionsMod ! The following dimension names must be replicated in ! CLM/ALMs histFileMod.F90 and - + character(*), parameter, public :: levcapf = 'fates_levcapf' ! matches histFileMod + character(*), parameter, public :: levcacls = 'fates_levcacls' ! matches histFileMod + character(*), parameter, public :: cohort = 'cohort' ! matches clm_varcon character(*), parameter, public :: patch = 'patch' ! matches clm_varcon character(*), parameter, public :: column = 'column' ! matches clm_varcon @@ -30,7 +32,7 @@ module FatesIODimensionsMod character(*), parameter, public :: levelpft = 'fates_levelpft' character(*), parameter, public :: levelcwd = 'fates_levelcwd' character(*), parameter, public :: levelage = 'fates_levelage' - + ! patch = This is a structure that records where FATES patch boundaries ! on each thread point to in the host IO array, this structure ! is allocated by number of threads @@ -45,9 +47,15 @@ module FatesIODimensionsMod ! levscpf = This is a structure that records the boundaries for the ! number of size-class x pft dimension + ! levcapf = This is a structure that records the boundaries for the + ! number of cohort-age-class x pft dimension + ! levscls = This is a structure that records the boundaries for the ! number of size-class dimension + ! levcacls = This is a structure that records the boundaries for the + ! number of cohort age class dimension + ! levpft = This is a structure that records the boundaries for the ! number of pft dimension @@ -105,8 +113,12 @@ module FatesIODimensionsMod integer :: agepft_class_end integer :: sizepft_class_begin integer :: sizepft_class_end + integer :: coagepf_class_begin + integer :: coagepf_class_end integer :: size_class_begin integer :: size_class_end + integer :: coage_class_begin + integer :: coage_class_end integer :: pft_class_begin integer :: pft_class_end integer :: age_class_begin diff --git a/main/FatesIOVariableKindMod.F90 b/main/FatesIOVariableKindMod.F90 index 7d6ae688c3..a46b020370 100644 --- a/main/FatesIOVariableKindMod.F90 +++ b/main/FatesIOVariableKindMod.F90 @@ -18,6 +18,8 @@ module FatesIOVariableKindMod character(*), parameter, public :: site_ground_r8 = 'SI_GRND_R8' character(*), parameter, public :: site_size_pft_r8 = 'SI_SCPF_R8' character(*), parameter, public :: site_size_r8 = 'SI_SCLS_R8' + character(*), parameter, public :: site_coage_pft_r8 = 'SI_CAPF_R8' + character(*), parameter, public :: site_coage_r8 = 'SI_CACLS_R8' character(*), parameter, public :: patch_int = 'PA_INT' character(*), parameter, public :: cohort_r8 = 'CO_R8' character(*), parameter, public :: cohort_int = 'CO_INT' diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 0254c89123..351bd92f7b 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -201,6 +201,11 @@ module FatesInterfaceMod ! well. ! ------------------------------------------------------------------------------------- + real(r8), public, allocatable :: fates_hdim_levcoage(:) ! cohort age class lower bound dimension + integer , public, allocatable :: fates_hdim_pfmap_levcapf(:) ! map of pfts into cohort age class x pft dimension + integer , public, allocatable :: fates_hdim_camap_levcapf(:) ! map of cohort age class into cohort age x pft dimension + + real(r8), public, allocatable :: fates_hdim_levsclass(:) ! plant size class lower bound dimension integer , public, allocatable :: fates_hdim_pfmap_levscpf(:) ! map of pfts into size-class x pft dimension integer , public, allocatable :: fates_hdim_scmap_levscpf(:) ! map of size-class into size-class x pft dimension @@ -231,7 +236,6 @@ module FatesInterfaceMod integer , public, allocatable :: fates_hdim_cwdmap_levelcwd(:) ! map of cwds in the element x cwd dimension integer , public, allocatable :: fates_hdim_agemap_levelage(:) ! map of ages in the element x age dimension - ! ------------------------------------------------------------------------------------ ! DYNAMIC BOUNDARY CONDITIONS ! ------------------------------------------------------------------------------------ @@ -266,6 +270,7 @@ module FatesInterfaceMod integer, public, protected :: nlevsclass ! The total number of cohort size class bins output to history integer, public, protected :: nlevage ! The total number of patch age bins output to history integer, public, protected :: nlevheight ! The total number of height bins output to history + integer, public, protected :: nlevcoage ! The total number of cohort age bins output to history integer, public, protected :: nleafage ! The total number of leaf age classes ! ------------------------------------------------------------------------------------- @@ -1037,6 +1042,7 @@ subroutine set_fates_global_elements(use_fates) use EDParamsMod, only : ED_val_history_sizeclass_bin_edges, ED_val_history_ageclass_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges + use EDParamsMod, only : ED_val_history_coageclass_bin_edges use CLMFatesParamInterfaceMod , only : FatesReadParameters implicit none @@ -1099,6 +1105,7 @@ subroutine set_fates_global_elements(use_fates) nlevsclass = size(ED_val_history_sizeclass_bin_edges,dim=1) nlevage = size(ED_val_history_ageclass_bin_edges,dim=1) nlevheight = size(ED_val_history_height_bin_edges,dim=1) + nlevcoage = size(ED_val_history_coageclass_bin_edges,dim=1) ! do some checks on the size, age, and height bin arrays to make sure they make sense: ! make sure that all start at zero, and that both are monotonically increasing @@ -1114,6 +1121,10 @@ subroutine set_fates_global_elements(use_fates) write(fates_log(), *) 'height class bins specified in parameter file must start at zero' call endrun(msg=errMsg(sourcefile, __LINE__)) endif + if ( ED_val_history_coageclass_bin_edges(1) .ne. 0._r8 ) then + write(fates_log(), *) 'cohort age class bines specified in parameter file must start at zero' + call endrun(msg=errMsg(sourcefile, __LINE__)) + endif do i = 2,nlevsclass if ( (ED_val_history_sizeclass_bin_edges(i) - ED_val_history_sizeclass_bin_edges(i-1)) .le. 0._r8) then write(fates_log(), *) 'size class bins specified in parameter file must be monotonically increasing' @@ -1132,6 +1143,12 @@ subroutine set_fates_global_elements(use_fates) call endrun(msg=errMsg(sourcefile, __LINE__)) end if end do + do i = 2,nlevcoage + if ( (ED_val_history_coageclass_bin_edges(i) - ED_val_history_coageclass_bin_edges(i-1)) .le. 0._r8) then + write(fates_log(), *) 'cohort age class bins specified in parameter file must be monotonically increasing' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + end do ! Set Various Mapping Arrays used in history output as well ! These will not be used if use_ed or use_fates is false @@ -1163,10 +1180,11 @@ subroutine fates_history_maps use EDParamsMod, only : ED_val_history_sizeclass_bin_edges use EDParamsMod, only : ED_val_history_ageclass_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges + use EDParamsMod, only : ED_val_history_coageclass_bin_edges ! ------------------------------------------------------------------------------------------ ! This subroutine allocates and populates the variables - ! that define the mapping of variables in history files in multiplexed dimensions liked + ! that define the mapping of variables in history files in multiplexed dimensions like ! the "scpf" format ! back to ! their respective single component dimensions, like size-class "sc" and pft "pf" @@ -1181,6 +1199,7 @@ subroutine fates_history_maps integer :: ileaf integer :: iage integer :: iheight + integer :: icoage integer :: iel allocate( fates_hdim_levsclass(1:nlevsclass )) @@ -1191,6 +1210,9 @@ subroutine fates_history_maps allocate( fates_hdim_levcwdsc(1:NCWD )) allocate( fates_hdim_levage(1:nlevage )) allocate( fates_hdim_levheight(1:nlevheight )) + allocate( fates_hdim_levcoage(1:nlevcoage )) + allocate( fates_hdim_pfmap_levcapf(1:nlevcoage*numpft)) + allocate( fates_hdim_camap_levcapf(1:nlevcoage*numpft)) allocate( fates_hdim_levcan(nclmax)) allocate( fates_hdim_levelem(num_elements)) @@ -1219,6 +1241,7 @@ subroutine fates_history_maps fates_hdim_levsclass(:) = ED_val_history_sizeclass_bin_edges(:) fates_hdim_levage(:) = ED_val_history_ageclass_bin_edges(:) fates_hdim_levheight(:) = ED_val_history_height_bin_edges(:) + fates_hdim_levcoage(:) = ED_val_history_coageclass_bin_edges(:) ! make pft array do ipft=1,numpft @@ -1283,6 +1306,15 @@ subroutine fates_history_maps end do end do + i=0 + do ipft=1,numpft + do icoage=1,nlevcoage + i=i+1 + fates_hdim_pfmap_levcapf(i) = ipft + fates_hdim_camap_levcapf(i) = icoage + end do + end do + i=0 do ican=1,nclmax do ileaf=1,nlevleaf @@ -1396,7 +1428,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) ! ! RGK-2016 ! --------------------------------------------------------------------------------- - + use EDParamsMod , only : ED_val_cohort_age_fusion_tol + ! Arguments integer, optional, intent(in) :: ival real(r8), optional, intent(in) :: rval @@ -1507,6 +1540,16 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + if ( hlm_use_inventory_init.eq.1 .and. ED_val_cohort_age_fusion_tol > 0.0_r8) then + if (fates_global_verbose()) then + write(fates_log(), *) 'Fates inventory init cannot be used with age dependent mortality' + write(fates_log(), *) 'Set fates_cohort_age_fusion_tol to 0 or turn off inventory init' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if ( .not.((hlm_use_inventory_init.eq.1).or.(hlm_use_inventory_init.eq.0)) ) then if (fates_global_verbose()) then diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 8fc3524789..8457d45963 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -1112,8 +1112,10 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & call prt_obj%CheckInitialConditions() + ! Since spread is a canopy level calculation, we need to provide an initial guess here. - call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & +call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & + temp_cohort%coage, temp_cohort%dbh, & prt_obj, temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & 1, csite%spread, bc_in) diff --git a/main/FatesParametersInterface.F90 b/main/FatesParametersInterface.F90 index fa3843cb31..ebaad3fa7c 100644 --- a/main/FatesParametersInterface.F90 +++ b/main/FatesParametersInterface.F90 @@ -34,7 +34,8 @@ module FatesParametersInterface character(len=*), parameter, public :: dimension_name_history_size_bins = 'fates_history_size_bins' character(len=*), parameter, public :: dimension_name_history_age_bins = 'fates_history_age_bins' character(len=*), parameter, public :: dimension_name_history_height_bins = 'fates_history_height_bins' - + character(len=*), parameter, public :: dimension_name_history_coage_bins = 'fates_history_coage_bins' + ! Dimensions in the host namespace: character(len=*), parameter, public :: dimension_name_host_allpfts = 'allpfts' diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index f0513d7e5c..f86bdc429a 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -14,6 +14,8 @@ module FatesRestartInterfaceMod use FatesIODimensionsMod, only : fates_io_dimension_type use FatesIOVariableKindMod, only : fates_io_variable_kind_type use FatesRestartVariableMod, only : fates_restart_variable_type + use FatesInterfaceMod, only : nlevcoage + use FatesInterfaceMod, only : bc_in_type use FatesInterfaceMod, only : bc_out_type use FatesInterfaceMod, only : hlm_use_planthydro @@ -90,6 +92,7 @@ module FatesRestartInterfaceMod integer :: ir_canopy_trim_co integer :: ir_size_class_lasttimestep_co integer :: ir_dbh_co + integer :: ir_coage_co integer :: ir_g_sb_laweight_co integer :: ir_height_co integer :: ir_laimemory_co @@ -104,6 +107,8 @@ module FatesRestartInterfaceMod integer :: ir_hmort_co integer :: ir_cmort_co integer :: ir_frmort_co + integer :: ir_smort_co + integer :: ir_asmort_co !Logging integer :: ir_lmort_direct_co @@ -139,7 +144,6 @@ module FatesRestartInterfaceMod ! Site level - integer :: ir_watermem_siwm integer :: ir_vegtempmem_sitm integer :: ir_seed_bank_sift @@ -649,6 +653,10 @@ subroutine define_restart_vars(this, initialize_variables) long_name='ed cohort - diameter at breast height', units='cm', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_dbh_co ) + call this%set_restart_var(vname='fates_coage', vtype=cohort_r8, & + long_name='ed cohort - age in days', units='days', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_coage_co ) + call this%set_restart_var(vname='fates_height', vtype=cohort_r8, & long_name='ed cohort - plant height', units='m', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_height_co ) @@ -713,6 +721,16 @@ subroutine define_restart_vars(this, initialize_variables) units='/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_frmort_co ) + call this%set_restart_var(vname='fates_smort', vtype=cohort_r8, & + long_name='ed cohort - senescence mortality rate', & + units='/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_smort_co ) + + call this%set_restart_var(vname='fates_asmort', vtype=cohort_r8, & + long_name='ed cohort - age senescence mortality rate', & + units = '/year', flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_asmort_co ) + call this%set_restart_var(vname='fates_lmort_direct', vtype=cohort_r8, & long_name='ed cohort - directly logging mortality rate', & units='%/event', flushval = flushzero, & @@ -1413,10 +1431,13 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_si_lyr_shell ! site - layer x shell index integer :: io_idx_si_scpf ! each size-class x pft index within site integer :: io_idx_si_sc ! each size-class index within site + integer :: io_idx_si_capf ! each cohort age-class x pft index within site + integer :: io_idx_si_cacls ! each cohort age class index within site integer :: io_idx_si_cwd ! each site-cwd index integer :: io_idx_si_pft ! each site-pft index integer :: io_idx_si_vtmem ! indices for veg-temp memory at site + ! Some counters (for checking mostly) integer :: totalcohorts ! total cohort count on this thread (diagnostic) integer :: patchespersite ! number of patches per site @@ -1431,6 +1452,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: i_var ! loop counter for PRT variables integer :: i_pos ! loop counter for discrete PRT positions integer :: i_scls ! loop counter for size-class + integer :: i_cacls ! loop counter for cohort age class integer :: i_cwd ! loop counter for cwd integer :: i_pft ! loop counter for pft @@ -1460,6 +1482,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_seed_prod_co => this%rvars(ir_seed_prod_co)%r81d, & rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & rio_dbh_co => this%rvars(ir_dbh_co)%r81d, & + rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & @@ -1473,6 +1496,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & + rio_smort_co => this%rvars(ir_smort_co)%r81d, & + rio_asmort_co => this%rvars(ir_asmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & @@ -1540,7 +1565,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_lyr_shell = io_idx_co_1st io_idx_si_scpf = io_idx_co_1st io_idx_si_sc = io_idx_co_1st - + io_idx_si_capf = io_idx_co_1st + io_idx_si_cacls= io_idx_co_1st ! recruitment rate do i_pft = 1,numpft @@ -1661,6 +1687,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_seed_prod_co(io_idx_co) = ccohort%seed_prod rio_size_class_lasttimestep(io_idx_co) = ccohort%size_class_lasttimestep rio_dbh_co(io_idx_co) = ccohort%dbh + rio_coage_co(io_idx_co) = ccohort%coage rio_height_co(io_idx_co) = ccohort%hite rio_laimemory_co(io_idx_co) = ccohort%laimemory rio_g_sb_laweight_co(io_idx_co)= ccohort%g_sb_laweight @@ -1676,6 +1703,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_bmort_co(io_idx_co) = ccohort%bmort rio_hmort_co(io_idx_co) = ccohort%hmort rio_cmort_co(io_idx_co) = ccohort%cmort + rio_smort_co(io_idx_co) = ccohort%smort + rio_asmort_co(io_idx_co) = ccohort%asmort rio_frmort_co(io_idx_co) = ccohort%frmort !Logging @@ -1813,7 +1842,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_termnindiv_cano_siscpf(io_idx_si_scpf) = sites(s)%term_nindivs_canopy(i_scls,i_pft) rio_termnindiv_usto_siscpf(io_idx_si_scpf) = sites(s)%term_nindivs_ustory(i_scls,i_pft) rio_growflx_fusion_siscpf(io_idx_si_scpf) = sites(s)%growthflux_fusion(i_scls, i_pft) - + io_idx_si_scpf = io_idx_si_scpf + 1 end do @@ -1822,8 +1851,16 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_sc = io_idx_si_sc + 1 end do - + +! do i_cacls = 1, nlevcoage + ! do i_pft = 1, numpft + ! rio_ageflx_fusion_sicapf(io_idx_si_capf) = sites(s)%ageflux_fusion(i_cacls, i_pft) + ! io_idx_si_capf = io_idx_si_capf + 1 + ! end do + ! io_idx_si_cacls = io_idx_si_cacls + 1 + ! end do + rio_termcflux_cano_si(io_idx_si) = sites(s)%term_carbonflux_canopy rio_termcflux_usto_si(io_idx_si) = sites(s)%term_carbonflux_ustory rio_democflux_si(io_idx_si) = sites(s)%demotion_carbonflux @@ -2147,6 +2184,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_si_lyr_shell ! site - layer x shell index integer :: io_idx_si_scpf ! each size-class x pft index within site integer :: io_idx_si_sc ! each size-class index within site + integer :: io_idx_si_cacls ! each coage class index within site + integer :: io_idx_si_capf ! each cohort age class x pft index within site integer :: io_idx_si_cwd integer :: io_idx_si_pft @@ -2162,7 +2201,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: i_var ! loop counter for PRT variables integer :: i_pos ! loop counter for discrete PRT positions integer :: i_pft ! loop counter for pft - integer :: i_scls ! loop counter for size-class + integer :: i_scls ! loop counter for size-clas + integer :: i_cacls ! loop counter for cohort age class associate( rio_npatch_si => this%rvars(ir_npatch_si)%int1d, & rio_cd_status_si => this%rvars(ir_cd_status_si)%int1d, & @@ -2185,6 +2225,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_seed_prod_co => this%rvars(ir_seed_prod_co)%r81d, & rio_size_class_lasttimestep => this%rvars(ir_size_class_lasttimestep_co)%int1d, & rio_dbh_co => this%rvars(ir_dbh_co)%r81d, & + rio_coage_co => this%rvars(ir_coage_co)%r81d, & rio_g_sb_laweight_co => this%rvars(ir_g_sb_laweight_co)%r81d, & rio_height_co => this%rvars(ir_height_co)%r81d, & rio_laimemory_co => this%rvars(ir_laimemory_co)%r81d, & @@ -2198,6 +2239,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_bmort_co => this%rvars(ir_bmort_co)%r81d, & rio_hmort_co => this%rvars(ir_hmort_co)%r81d, & rio_cmort_co => this%rvars(ir_cmort_co)%r81d, & + rio_smort_co => this%rvars(ir_smort_co)%r81d, & + rio_asmort_co => this%rvars(ir_asmort_co)%r81d, & rio_frmort_co => this%rvars(ir_frmort_co)%r81d, & rio_lmort_direct_co => this%rvars(ir_lmort_direct_co)%r81d, & rio_lmort_collateral_co => this%rvars(ir_lmort_collateral_co)%r81d, & @@ -2254,7 +2297,9 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_scpf = io_idx_co_1st io_idx_si_sc = io_idx_co_1st - + io_idx_si_capf = io_idx_co_1st + io_idx_si_cacls= io_idx_co_1st + ! read seed_bank info(site-level, but PFT-resolved) do i_pft = 1,numpft sites(s)%recruitment_rate(i_pft) = rio_recrate_sift(io_idx_co_1st+i_pft-1) @@ -2348,6 +2393,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%seed_prod = rio_seed_prod_co(io_idx_co) ccohort%size_class_lasttimestep = rio_size_class_lasttimestep(io_idx_co) ccohort%dbh = rio_dbh_co(io_idx_co) + ccohort%coage = rio_coage_co(io_idx_co) ccohort%g_sb_laweight= rio_g_sb_laweight_co(io_idx_co) ccohort%hite = rio_height_co(io_idx_co) ccohort%laimemory = rio_laimemory_co(io_idx_co) @@ -2362,6 +2408,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ccohort%bmort = rio_bmort_co(io_idx_co) ccohort%hmort = rio_hmort_co(io_idx_co) ccohort%cmort = rio_cmort_co(io_idx_co) + ccohort%smort = rio_smort_co(io_idx_co) + ccohort%asmort = rio_asmort_co(io_idx_co) ccohort%frmort = rio_frmort_co(io_idx_co) !Logging @@ -2574,8 +2622,16 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_sc = io_idx_si_sc + 1 end do - + ! do i_cacls = i,nlevcoage + ! do i_pft = 1, numpft + ! sites(s)%ageflux_fusion(i_cacls, i_pft) = rio_ageflx_fusion_sicapf(io_idx_si_capf) + ! io_idx_si_capf = io_idx_si_capf + 1 + ! end do + ! io_idx_si_cacls = io_idx_si_cacls + 1 + ! end do + + sites(s)%term_carbonflux_canopy = rio_termcflux_cano_si(io_idx_si) sites(s)%term_carbonflux_ustory = rio_termcflux_usto_si(io_idx_si) sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) diff --git a/main/FatesSizeAgeTypeIndicesMod.F90 b/main/FatesSizeAgeTypeIndicesMod.F90 index 383621a55e..5683cc5302 100644 --- a/main/FatesSizeAgeTypeIndicesMod.F90 +++ b/main/FatesSizeAgeTypeIndicesMod.F90 @@ -4,9 +4,11 @@ module FatesSizeAgeTypeIndicesMod use FatesInterfaceMod, only : nlevsclass use FatesInterfaceMod, only : nlevage use FatesInterfaceMod, only : nlevheight + use FatesInterfaceMod, only : nlevcoage use EDParamsMod, only : ED_val_history_sizeclass_bin_edges use EDParamsMod, only : ED_val_history_ageclass_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges + use EDParamsMod, only : ED_val_history_coageclass_bin_edges implicit none private ! Modules are private by default @@ -19,6 +21,8 @@ module FatesSizeAgeTypeIndicesMod public :: get_height_index public :: get_sizeagepft_class_index public :: get_agepft_class_index + public :: coagetype_class_index + public :: get_coage_class_index contains @@ -83,6 +87,37 @@ function get_size_class_index(dbh) result(cohort_size_class) end function get_size_class_index + ! ===================================================================================== + + subroutine coagetype_class_index(coage,pft,coage_class,coage_by_pft_class) + + ! Arguments + real(r8),intent(in) :: coage + integer,intent(in) :: pft + integer,intent(out) :: coage_class + integer,intent(out) :: coage_by_pft_class + + coage_class = get_coage_class_index(coage) + + coage_by_pft_class = (pft-1)*nlevcoage+coage_class + + return + end subroutine coagetype_class_index + + ! ======================================================================================== + + function get_coage_class_index(coage) result(cohort_coage_class) + + real(r8), intent(in) :: coage + + integer :: cohort_coage_class + + cohort_coage_class = count(coage-ED_val_history_coageclass_bin_edges.ge.0.0_r8) + + end function get_coage_class_index + + + ! ===================================================================================== function get_height_index(height) result(cohort_height_index) diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 3573dd4feb..e6207ce634 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -25,6 +25,9 @@ variables: double fates_history_sizeclass_bin_edges(fates_history_size_bins) ; fates_history_sizeclass_bin_edges:units = "cm" ; fates_history_sizeclass_bin_edges:long_name = "Lower edges for DBH size class bins used in size-resolved cohort history output" ; + double fates_history_coageclass_bin_edges(fates_history_coage_bins) ; + fates_history_coageclass_bin_edges:units = "years" ; + fates_history_coageclass_bin_edges:long_name = "Lower edges for cohort age class bins used in cohort age resolved history output" ; char fates_pftname(fates_pft, fates_string_length) ; fates_pftname:units = "unitless - string" ; fates_pftname:long_name = "Description of plant type" ; @@ -301,6 +304,18 @@ variables: double fates_mort_bmort(fates_pft) ; fates_mort_bmort:units = "1/yr" ; fates_mort_bmort:long_name = "background mortality rate" ; + double fates_mort_ip_size_senescence(fates_pft); + fates_mort_ip_size_senescence:units = "dbh cm"; + fates_mort_ip_size_senescence:long_name ="Mortality dbh senescence inflection point"; + double fates_mort_r_size_senescence(fates_pft); + fates_mort_r_size_senescence:units = "mortality rate dbh^-1"; + fates_mort_r_size_senescence:long_name = "Mortality dbh senescence rate of change"; + double fates_mort_ip_age_senescence(fates_pft); + fates_mort_ip_age_senescence:units = "years"; + fates_mort_ip_age_senescence:long_name ="Mortality cohort age senescence inflection point"; + double fates_mort_r_age_senescence(fates_pft); + fates_mort_r_age_senescence:units = "mortality rate year^-1"; + fates_mort_r_age_senescence:long_name = "Mortality age senescence rate of change"; double fates_mort_freezetol(fates_pft) ; fates_mort_freezetol:units = "NA" ; fates_mort_freezetol:long_name = "minimum temperature tolerance (NOT USED)" ; @@ -531,7 +546,7 @@ variables: fates_canopy_closure_thresh:long_name = "tree canopy coverage at which crown area allometry changes from savanna to forest value" ; double fates_cohort_age_fusion_tol ; fates_cohort_age_fusion_tol:units = "unitless" ; - fates_cohort_age_fusion_tol:long_name = "minimum fraction in differece in cohort age between cohorts. 0 or _ implies functionality is turned off." ; + fates_cohort_age_fusion_tol:long_name = "minimum fraction in differece in cohort age between cohorts. 0 implies functionality is turned off." ; double fates_cohort_size_fusion_tol ; fates_cohort_size_fusion_tol:units = "unitless" ; fates_cohort_size_fusion_tol:long_name = "minimum fraction in difference in dbh between cohorts" ; @@ -684,6 +699,8 @@ data: fates_history_sizeclass_bin_edges = 0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100 ; + fates_history_coageclass_bin_edges = 0, 5 ; + fates_pftname = "broadleaf_evergreen_tropical_tree ", "needleleaf_evergreen_extratrop_tree ", @@ -984,6 +1001,14 @@ data: fates_mort_bmort = 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014 ; + fates_mort_ip_size_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; + + fates_mort_r_size_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; + + fates_mort_ip_age_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; + + fates_mort_r_age_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; + fates_mort_freezetol = 2.5, -55, -80, -30, 2.5, -30, -60, -10, -80, -80, -20, 2.5 ; @@ -1136,6 +1161,7 @@ data: fates_seed_suppl = 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ; + fates_senleaf_long_fdrought = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ; fates_smpsc = -255000, -255000, -255000, -255000, -255000, -255000, -255000, diff --git a/tools/FatesPFTIndexSwapper.py b/tools/FatesPFTIndexSwapper.py index 091ffc30ce..be467a6a83 100755 --- a/tools/FatesPFTIndexSwapper.py +++ b/tools/FatesPFTIndexSwapper.py @@ -10,12 +10,13 @@ # ======================================================================================= import numpy as np +from numpy import * import sys import getopt import code # For development: code.interact(local=locals()) from datetime import datetime from scipy.io import netcdf -import matplotlib.pyplot as plt +#import matplotlib.pyplot as plt # ======================================================================================= # Parameters @@ -144,7 +145,7 @@ def main(argv): fp_in = netcdf.netcdf_file(input_fname, 'r') - for key, value in sorted(fp_in.dimensions.iteritems()): + for key, value in sorted(fp_in.dimensions.items()): if(key==pft_dim_name): fp_out.createDimension(key,int(num_pft_out)) print('Creating Dimension: {}={}'.format(key,num_pft_out)) @@ -152,7 +153,7 @@ def main(argv): fp_out.createDimension(key,int(value)) print('Creating Dimension: {}={}'.format(key,value)) - for key, value in sorted(fp_in.variables.iteritems()): + for key, value in sorted(fp_in.variables.items()): print('Creating Variable: ',key) # code.interact(local=locals()) diff --git a/tools/modify_fates_paramfile.py b/tools/modify_fates_paramfile.py index c3cfa649fb..f2e1729de2 100755 --- a/tools/modify_fates_paramfile.py +++ b/tools/modify_fates_paramfile.py @@ -86,7 +86,7 @@ def main(): npft_file = var.shape[i] pftdim = i otherdimpresent = False - elif var.dimensions[i] in ['fates_history_age_bins','fates_history_size_bins','fates_history_height_bins','fates_NCWD','fates_litterclass','fates_leafage_class','fates_prt_organs','fates_hydr_organs','fates_variants']: + elif var.dimensions[i] in ['fates_history_age_bins','fates_history_size_bins','fates_history_coage_bins','fates_history_height_bins','fates_NCWD','fates_litterclass','fates_leafage_class','fates_prt_organs','fates_hydr_organs','fates_variants']: otherdimpresent = True otherdimname = var.dimensions[i] otherdimlength = var.shape[i] @@ -101,7 +101,7 @@ def main(): ### ok, we find ourselves in the situation where we need to rewrite the netcdf from scratch with its revised shape. # # first lets chech to make sure the dimension we are changing can be changed without breaking things. - plastic_dimensions_list = ['fates_history_age_bins','fates_history_size_bins','fates_history_height_bins','fates_leafage_class'] + plastic_dimensions_list = ['fates_history_age_bins','fates_history_size_bins','fates_history_coage_bins','fates_history_height_bins','fates_leafage_class'] if otherdimname not in plastic_dimensions_list: raise ValueError('asking to change the shape of a dimension, '+otherdimname+', that will probably break things') else: From 18bcd440e5d760b07a2c54db17e6ffa6d0d8f7db Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Thu, 9 Jan 2020 13:07:06 -0800 Subject: [PATCH 02/27] [ Fix bug in parameter file - cohort age classes must be defined ] [fates_history_coage_bins must be defined in the parameter file, even if age dependent mortality is off. This is because cohorts are still initialised with age 0 and there is a call to get cohort age class. The age is not incremented if age mortality is turned off but we still need this initial call. Cohort age bins must therefore be at least one bin wide. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No, just changing defaults] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Testing on my branch --- parameter_files/fates_params_default.cdl | 30 +++--------------------- 1 file changed, 3 insertions(+), 27 deletions(-) diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index e6207ce634..9219938b89 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -4,7 +4,7 @@ dimensions: fates_history_age_bins = 7 ; fates_history_height_bins = 6 ; fates_history_size_bins = 13 ; - fates_history_coage_bins = 1 ; + fates_history_coage_bins = 2 ; fates_hydr_organs = 4 ; fates_leafage_class = 1 ; fates_litterclass = 6 ; @@ -25,9 +25,6 @@ variables: double fates_history_sizeclass_bin_edges(fates_history_size_bins) ; fates_history_sizeclass_bin_edges:units = "cm" ; fates_history_sizeclass_bin_edges:long_name = "Lower edges for DBH size class bins used in size-resolved cohort history output" ; - double fates_history_coageclass_bin_edges(fates_history_coage_bins) ; - fates_history_coageclass_bin_edges:units = "years" ; - fates_history_coageclass_bin_edges:long_name = "Lower edges for cohort age class bins used in cohort age resolved history output" ; char fates_pftname(fates_pft, fates_string_length) ; fates_pftname:units = "unitless - string" ; fates_pftname:long_name = "Description of plant type" ; @@ -304,18 +301,6 @@ variables: double fates_mort_bmort(fates_pft) ; fates_mort_bmort:units = "1/yr" ; fates_mort_bmort:long_name = "background mortality rate" ; - double fates_mort_ip_size_senescence(fates_pft); - fates_mort_ip_size_senescence:units = "dbh cm"; - fates_mort_ip_size_senescence:long_name ="Mortality dbh senescence inflection point"; - double fates_mort_r_size_senescence(fates_pft); - fates_mort_r_size_senescence:units = "mortality rate dbh^-1"; - fates_mort_r_size_senescence:long_name = "Mortality dbh senescence rate of change"; - double fates_mort_ip_age_senescence(fates_pft); - fates_mort_ip_age_senescence:units = "years"; - fates_mort_ip_age_senescence:long_name ="Mortality cohort age senescence inflection point"; - double fates_mort_r_age_senescence(fates_pft); - fates_mort_r_age_senescence:units = "mortality rate year^-1"; - fates_mort_r_age_senescence:long_name = "Mortality age senescence rate of change"; double fates_mort_freezetol(fates_pft) ; fates_mort_freezetol:units = "NA" ; fates_mort_freezetol:long_name = "minimum temperature tolerance (NOT USED)" ; @@ -546,7 +531,7 @@ variables: fates_canopy_closure_thresh:long_name = "tree canopy coverage at which crown area allometry changes from savanna to forest value" ; double fates_cohort_age_fusion_tol ; fates_cohort_age_fusion_tol:units = "unitless" ; - fates_cohort_age_fusion_tol:long_name = "minimum fraction in differece in cohort age between cohorts. 0 implies functionality is turned off." ; + fates_cohort_age_fusion_tol:long_name = "minimum fraction in differece in cohort age between cohorts. _ implies functionality is turned off." ; double fates_cohort_size_fusion_tol ; fates_cohort_size_fusion_tol:units = "unitless" ; fates_cohort_size_fusion_tol:long_name = "minimum fraction in difference in dbh between cohorts" ; @@ -692,14 +677,13 @@ data: fates_history_ageclass_bin_edges = 0, 1, 2, 5, 10, 20, 50 ; - fates_history_coageclass_bin_edges = _ ; + fates_history_coageclass_bin_edges = 0, 5 ; fates_history_height_bin_edges = 0, 0.1, 0.3, 1, 3, 10 ; fates_history_sizeclass_bin_edges = 0, 5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100 ; - fates_history_coageclass_bin_edges = 0, 5 ; fates_pftname = "broadleaf_evergreen_tropical_tree ", @@ -1001,14 +985,6 @@ data: fates_mort_bmort = 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014, 0.014 ; - fates_mort_ip_size_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; - - fates_mort_r_size_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; - - fates_mort_ip_age_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; - - fates_mort_r_age_senescence = _, _, _, _, _, _, _, _, _, _, _, _ ; - fates_mort_freezetol = 2.5, -55, -80, -30, 2.5, -30, -60, -10, -80, -80, -20, 2.5 ; From 7b3490362a18f4dfa80355c3ff79cb12c1c1dfb6 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Thu, 9 Jan 2020 13:42:02 -0800 Subject: [PATCH 03/27] [ Update FatesPFTindexSwapper to be compatible with python 3] [ Iteritems to items ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: --- tools/FatesPFTIndexSwapper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/FatesPFTIndexSwapper.py b/tools/FatesPFTIndexSwapper.py index 091ffc30ce..d14b3f0b8c 100755 --- a/tools/FatesPFTIndexSwapper.py +++ b/tools/FatesPFTIndexSwapper.py @@ -144,7 +144,7 @@ def main(argv): fp_in = netcdf.netcdf_file(input_fname, 'r') - for key, value in sorted(fp_in.dimensions.iteritems()): + for key, value in sorted(fp_in.dimensions.items()): if(key==pft_dim_name): fp_out.createDimension(key,int(num_pft_out)) print('Creating Dimension: {}={}'.format(key,num_pft_out)) @@ -152,7 +152,7 @@ def main(argv): fp_out.createDimension(key,int(value)) print('Creating Dimension: {}={}'.format(key,value)) - for key, value in sorted(fp_in.variables.iteritems()): + for key, value in sorted(fp_in.variables.items()): print('Creating Variable: ',key) # code.interact(local=locals()) From cfb431c3fb1b3b13510002d555d66b3d3e7c7a92 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Thu, 9 Jan 2020 14:44:30 -0800 Subject: [PATCH 04/27] [ Change the logic for how the size and age mortality functions are switched on] [Default is for all params related to size and age mortality to be _. In FATES these get interpreted as a very large number. Therefore, if a user specifies param values for size or age mortality, we can use a condition that these values are < 1000 which will allow size and age mortality functionality to be turned on. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Tested on my branch. --- biogeochem/EDCanopyStructureMod.F90 | 2 +- biogeochem/EDCohortDynamicsMod.F90 | 2 +- biogeochem/EDMortalityFunctionsMod.F90 | 9 ++++----- main/EDMainMod.F90 | 2 +- main/FatesInterfaceMod.F90 | 4 ++-- 5 files changed, 9 insertions(+), 10 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 8f473d440a..fa7282aaac 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -1325,7 +1325,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) - if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then call coagetype_class_index(currentCohort%coage,currentCohort%pft, & currentCohort%coage_class,currentCohort%coage_by_pft_class) end if diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 7129a4d034..721e5c63a6 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1082,7 +1082,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) (nextc%coage * (nextc%n/(currentCohort%n + nextc%n))) ! update the cohort age again - if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then call coagetype_class_index(currentCohort%coage, currentCohort%pft, & currentCohort%coage_class, currentCohort%coage_by_pft_class) end if diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 0ea929c85a..1478e8312b 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -51,7 +51,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! ============================================================================ use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - + use FatesInterfaceMod , only : hlm_hio_ignore_val type (ed_cohort_type), intent(in) :: cohort_in type (bc_in_type), intent(in) :: bc_in @@ -88,8 +88,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor mort_r_size_senescence = EDPftvarcon_inst%mort_r_size_senescence(cohort_in%pft) mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) - ! if the user has set the inflection point to be a huge number then - ! we zero smort (even though it would essentially be zero anyway) + ! if param values have been set then calculate smort if (mort_ip_size_senescence < 10000 ) then smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & (cohort_in%dbh - mort_ip_size_senescence) ) ) @@ -97,9 +96,9 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor smort = 0.0_r8 end if - ! if user has set cohort age fusion param > 0 we calculate age + ! if user has set cohort age fusion param to not be _ we calculate age ! dependent mortality - if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then ! Age Dependent Senescence ! rate and inflection point define the change in mortality with age mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index a70932029d..b32e6a4236 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -442,7 +442,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) end if ! if we are in age-dependent mortality mode - if (ED_val_cohort_age_fusion_tol > 0.0_r8) then + if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then ! update cohort age currentCohort%coage = currentCohort%coage + hlm_freq_day if(currentCohort%coage < 0.0_r8)then diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 351bd92f7b..b971d11742 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1541,10 +1541,10 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( hlm_use_inventory_init.eq.1 .and. ED_val_cohort_age_fusion_tol > 0.0_r8) then + if ( hlm_use_inventory_init.eq.1 .and. ED_val_cohort_age_fusion_tol < 10000.0_r8) then if (fates_global_verbose()) then write(fates_log(), *) 'Fates inventory init cannot be used with age dependent mortality' - write(fates_log(), *) 'Set fates_cohort_age_fusion_tol to 0 or turn off inventory init' + write(fates_log(), *) 'Set fates_cohort_age_fusion_tol to _ or turn off inventory init' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if From 17a60ef90f10994ae76467d05ed3402ffa887afd Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Fri, 10 Jan 2020 14:20:49 -0800 Subject: [PATCH 05/27] [Fix bug in FatesHistoryInterfaceMod ] [M9_SCLS and M10_SCLS should both be written out on site_size dimension. M10_CACLS should be written out on site_coage dimension. These were mixed up.] Fixes: [NGT-ED Github issue #] User interface changes?: [ No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Running tests on my branch. --- main/FatesHistoryInterfaceMod.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 9fcb35023e..f8dd0f814a 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -4996,12 +4996,12 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='M10_SCLS', units = 'N/ha/yr', & long='age senescence mortality by size',use_default='active', & - avgflag='A', vtype=site_coage_r8, hlms='CLM:ALM', flushval=0.0_r8, & + avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m10_si_scls ) call this%set_history_var(vname='M10_CACLS', units = 'N/ha/yr', & long='age senescence mortality by cohort age',use_default='active', & - avgflag='A', vtype=site_size_r8, hlms='CLM:ALM', flushval=0.0_r8, & + avgflag='A', vtype=site_coage_r8, hlms='CLM:ALM', flushval=0.0_r8, & upfreq=1, ivar=ivar, initialize=initialize_variables, index = ih_m10_si_cacls ) call this%set_history_var(vname='CARBON_BALANCE_CANOPY_SCLS', units = 'kg C / ha / yr', & From 07898674a79ac3a056db932bc9efaaf7f4044d4d Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Fri, 17 Jan 2020 14:04:17 -0800 Subject: [PATCH 06/27] [Fix the small simple things like indentations and forgotten _r8s ] [ Fixing some of indentation issues and forgotten _r8s. Also delete commented out code that was not needed. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: No testing. --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- biogeochem/EDMortalityFunctionsMod.F90 | 8 +- main/EDInitMod.F90 | 1 - main/EDTypesMod.F90 | 1 - main/FatesHistoryInterfaceMod.F90 | 12 --- main/FatesInventoryInitMod.F90 | 116 ++++++++++++------------- main/FatesRestartInterfaceMod.F90 | 19 +--- 7 files changed, 64 insertions(+), 95 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 721e5c63a6..22a6444289 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1028,7 +1028,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! and hence the age fusion criterion is met if (currentCohort%coage .eq. nextc%coage) then coage_diff = 0.0_r8 - else + else coage_diff = abs((currentCohort%coage - nextc%coage)/ & (0.5*(currentCohort%coage + nextc%coage))) end if diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 1478e8312b..9a1c466209 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -84,12 +84,12 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! and the like ! Size Dependent Senescence - ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh + ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh mort_r_size_senescence = EDPftvarcon_inst%mort_r_size_senescence(cohort_in%pft) mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) - ! if param values have been set then calculate smort - if (mort_ip_size_senescence < 10000 ) then + ! if param values have been set then calculate smort + if (mort_ip_size_senescence < 10000.0_r8 ) then smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & (cohort_in%dbh - mort_ip_size_senescence) ) ) else @@ -108,7 +108,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor else asmort = 0.0_r8 end if - + if (hlm_use_ed_prescribed_phys .eq. ifalse) then diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 285c5cae85..086bc54525 100644 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -202,7 +202,6 @@ subroutine zero_site( site_in ) ! fusoin-induced growth flux of individuals site_in%growthflux_fusion(:,:) = 0._r8 -! site_in%ageflux_fusion(:,:) = 0._r8 ! demotion/promotion info site_in%demotion_rate(:) = 0._r8 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index f7508b3984..43a8599f42 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -756,7 +756,6 @@ module EDTypesMod real(r8), allocatable :: growthflux_fusion(:,:) ! rate of individuals moving into a given size class bin ! due to fusion in a given day. on size x pft array -! real(r8), allocatable :: ageflux_fusion(:,:) ! rate of individuals movign into a given age class bin due to fusion- age x pft diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index f8dd0f814a..be3cd50e6e 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -2713,18 +2713,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) end do end do ! - - ! do i_pft = 1, numpft - ! do i_cacls = 1, nlevcoage -! i_capf = (i_pft -1)*nlevcoage + i_cacls -! fix this! - ! i_capf = 1 - ! hio_ageflux_fusion_si_capf(io_si,i_capf) = hio_ageflux_fusion_si_capf(io_si,i_capf) + & - ! sites(s)%ageflux_fusion(i_capf, i_pft) * days_per_year - ! end do - ! end do - - ! treat carbon flux from imort the same way hio_understory_mortality_carbonflux_si(io_si) = hio_understory_mortality_carbonflux_si(io_si) + & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 8457d45963..67da02d0a8 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -1053,71 +1053,71 @@ subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & prt_obj => null() call InitPRTObject(prt_obj) - + do el = 1,num_elements - - element_id = element_list(el) - - ! If this is carbon12, then the initialization is straight forward - ! otherwise, we use stoichiometric ratios - select case(element_id) - case(carbon12_element) - - m_struct = b_struct - m_leaf = b_leaf - m_fnrt = b_fnrt - m_sapw = b_sapw - m_store = b_store - m_repro = 0._r8 - - case(nitrogen_element) - - m_struct = b_struct*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,struct_organ) - m_leaf = b_leaf*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,leaf_organ) - m_fnrt = b_fnrt*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,fnrt_organ) - m_sapw = b_sapw*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,sapw_organ) - m_store = b_store*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,store_organ) - m_repro = 0._r8 - - case(phosphorus_element) - - m_struct = b_struct*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,struct_organ) - m_leaf = b_leaf*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,leaf_organ) - m_fnrt = b_fnrt*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,fnrt_organ) - m_sapw = b_sapw*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,sapw_organ) - m_store = b_store*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,store_organ) - m_repro = 0._r8 - end select - - select case(hlm_parteh_mode) - case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) - - ! Equally distribute leaf mass into available age-bins - do iage = 1,nleafage - call SetState(prt_obj,leaf_organ, element_id,m_leaf/real(nleafage,r8),iage) - end do - - call SetState(prt_obj,fnrt_organ, element_id, m_fnrt) - call SetState(prt_obj,sapw_organ, element_id, m_sapw) - call SetState(prt_obj,store_organ, element_id, m_store) - call SetState(prt_obj,struct_organ, element_id, m_struct) - call SetState(prt_obj,repro_organ, element_id, m_repro) - - case default - write(fates_log(),*) 'Unspecified PARTEH module during inventory intitialization' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end select - + + element_id = element_list(el) + + ! If this is carbon12, then the initialization is straight forward + ! otherwise, we use stoichiometric ratios + select case(element_id) + case(carbon12_element) + + m_struct = b_struct + m_leaf = b_leaf + m_fnrt = b_fnrt + m_sapw = b_sapw + m_store = b_store + m_repro = 0._r8 + + case(nitrogen_element) + + m_struct = b_struct*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,struct_organ) + m_leaf = b_leaf*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,leaf_organ) + m_fnrt = b_fnrt*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,fnrt_organ) + m_sapw = b_sapw*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,sapw_organ) + m_store = b_store*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,store_organ) + m_repro = 0._r8 + + case(phosphorus_element) + + m_struct = b_struct*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,struct_organ) + m_leaf = b_leaf*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,leaf_organ) + m_fnrt = b_fnrt*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,fnrt_organ) + m_sapw = b_sapw*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,sapw_organ) + m_store = b_store*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,store_organ) + m_repro = 0._r8 + end select + + select case(hlm_parteh_mode) + case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) + + ! Equally distribute leaf mass into available age-bins + do iage = 1,nleafage + call SetState(prt_obj,leaf_organ, element_id,m_leaf/real(nleafage,r8),iage) + end do + + call SetState(prt_obj,fnrt_organ, element_id, m_fnrt) + call SetState(prt_obj,sapw_organ, element_id, m_sapw) + call SetState(prt_obj,store_organ, element_id, m_store) + call SetState(prt_obj,struct_organ, element_id, m_struct) + call SetState(prt_obj,repro_organ, element_id, m_repro) + + case default + write(fates_log(),*) 'Unspecified PARTEH module during inventory intitialization' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end select + end do call prt_obj%CheckInitialConditions() ! Since spread is a canopy level calculation, we need to provide an initial guess here. -call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & - temp_cohort%coage, temp_cohort%dbh, & - prt_obj, temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & - 1, csite%spread, bc_in) + call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, & + temp_cohort%coage, temp_cohort%dbh, & + prt_obj, temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & + 1, csite%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index f86bdc429a..7ebc67ca8c 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -1851,15 +1851,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_sc = io_idx_si_sc + 1 end do - - -! do i_cacls = 1, nlevcoage - ! do i_pft = 1, numpft - ! rio_ageflx_fusion_sicapf(io_idx_si_capf) = sites(s)%ageflux_fusion(i_cacls, i_pft) - ! io_idx_si_capf = io_idx_si_capf + 1 - ! end do - ! io_idx_si_cacls = io_idx_si_cacls + 1 - ! end do rio_termcflux_cano_si(io_idx_si) = sites(s)%term_carbonflux_canopy rio_termcflux_usto_si(io_idx_si) = sites(s)%term_carbonflux_ustory @@ -2623,15 +2614,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_sc = io_idx_si_sc + 1 end do - ! do i_cacls = i,nlevcoage - ! do i_pft = 1, numpft - ! sites(s)%ageflux_fusion(i_cacls, i_pft) = rio_ageflx_fusion_sicapf(io_idx_si_capf) - ! io_idx_si_capf = io_idx_si_capf + 1 - ! end do - ! io_idx_si_cacls = io_idx_si_cacls + 1 - ! end do - - + sites(s)%term_carbonflux_canopy = rio_termcflux_cano_si(io_idx_si) sites(s)%term_carbonflux_ustory = rio_termcflux_usto_si(io_idx_si) sites(s)%demotion_carbonflux = rio_democflux_si(io_idx_si) From d8e7935542222ca2333d7007c5b394a461538331 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Sun, 19 Jan 2020 19:23:51 -0800 Subject: [PATCH 07/27] [Change the switch for age dependent mortality ] [Make the switch for age-dependent mortality ip_age_senescence to match size-dependent mortality which is turned on by setting ip_size_mortality. Cohort age tracking is turned on independently by the flag 'cohort_age_tracking'. This is a new parameter in the parameter file. Users specify 0 for no age tracking and 1 for cohort age tracking. This makes it easier to have cohort age tracking independent of age-dependent mortality - might be useful in the future for lifetime accumulation of pathogen damage or something. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes - new parameter in the param file. Default 0 means age tracking of cohorts is off. If the user makes it 1 cohort age tracking is turned on. Age-dependent mortality is now turned on by setting the ip_age_senescence parameter.] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Runs as expected on lawrencium. --- biogeochem/EDCanopyStructureMod.F90 | 7 ++++--- biogeochem/EDCohortDynamicsMod.F90 | 17 +++++++++------- biogeochem/EDMortalityFunctionsMod.F90 | 11 +++++----- main/EDMainMod.F90 | 26 ++++++++++++++---------- main/EDParamsMod.F90 | 17 +++++++++++++--- parameter_files/fates_params_default.cdl | 7 ++++++- 6 files changed, 54 insertions(+), 31 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index fa7282aaac..5b538d29b5 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -29,7 +29,7 @@ module EDCanopyStructureMod use FatesInterfaceMod , only : numpft use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use EDTypesMod , only : maxCohortsPerPatch - use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use EDParamsMod , only : cohort_age_tracking use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : all_carbon_elements @@ -1262,7 +1262,8 @@ subroutine canopy_summarization( nsites, sites, bc_in ) use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst - use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use EDParamsMod , only : cohort_age_tracking + use FatesConstantsMod , only : itrue ! !ARGUMENTS integer , intent(in) :: nsites @@ -1325,7 +1326,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) - if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then + if (cohort_age_tracking) then call coagetype_class_index(currentCohort%coage,currentCohort%pft, & currentCohort%coage_class,currentCohort%coage_by_pft_class) end if diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 22a6444289..c69eacf85e 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -37,6 +37,7 @@ module EDCohortDynamicsMod use EDTypesMod , only : site_fluxdiags_type use EDTypesMod , only : num_elements use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use EDParamsMod , only : cohort_age_tracking use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -940,6 +941,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! !USES: use EDParamsMod , only : ED_val_cohort_size_fusion_tol use EDParamsMod , only : ED_val_cohort_age_fusion_tol + use EDParamsMod , only : cohort_age_tracking + use FatesConstantsMod , only : itrue use FatesConstantsMod, only : days_per_year ! @@ -1026,12 +1029,12 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! if they are the same age we make diff 0- to avoid errors divding by zero !NB if cohort age tracking is off then the age of both should be 0 ! and hence the age fusion criterion is met - if (currentCohort%coage .eq. nextc%coage) then - coage_diff = 0.0_r8 - else - coage_diff = abs((currentCohort%coage - nextc%coage)/ & - (0.5*(currentCohort%coage + nextc%coage))) - end if + if (abs(currentCohort%coage - nextc%coage) 0.995_r8)then - bmort = (0.995_r8 - ((cmort+hmort+asmort+frmort+smort+dndt_logging)*hlm_freq_day)) & + > 1.0_r8)then + bmort = (1.0_r8 - ((cmort+hmort+asmort+frmort+smort+dndt_logging)*hlm_freq_day)) & /hlm_freq_day endif @@ -266,8 +265,8 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) ! cap on bmort for canopy layers - bmort is adjusted so that total mortality ! including disturbance fraction does not exceed 0.995 if(((cmort+hmort+bmort+frmort+smort+asmort)*hlm_freq_day) > & - 0.995_r8 - fates_mortality_disturbance_fraction)then - bmort = (0.995_r8 - fates_mortality_disturbance_fraction - & + 1.0_r8 - fates_mortality_disturbance_fraction)then + bmort = (1.0_r8 - fates_mortality_disturbance_fraction - & ((cmort+hmort+asmort+frmort+smort)*hlm_freq_day))/hlm_freq_day endif diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index b32e6a4236..804c40986a 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -55,6 +55,7 @@ module EDMainMod use EDTypesMod , only : element_pos use EDTypesMod , only : phen_dstat_moiston use EDTypesMod , only : phen_dstat_timeon + use EDParamsMod , only : cohort_age_tracking use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : primaryforest, secondaryforest use FatesConstantsMod , only : nearzero @@ -283,7 +284,8 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! FIX(SPM,032414) refactor so everything goes through interface ! ! !USES: - use EDParamsMod, only : ED_val_cohort_age_fusion_tol + use EDParamsMod, only : cohort_age_tracking + use FatesConstantsMod, only : itrue ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: currentSite @@ -442,18 +444,20 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) end if ! if we are in age-dependent mortality mode - if (ED_val_cohort_age_fusion_tol < 10000.0_r8) then - ! update cohort age - currentCohort%coage = currentCohort%coage + hlm_freq_day - if(currentCohort%coage < 0.0_r8)then - write(fates_log(),*) 'negative cohort age?',currentCohort%coage - end if - - ! update cohort age class and age x pft class + write(fates_log(),*) 'cohort_age_tracking',cohort_age_tracking + if (cohort_age_tracking) then + ! update cohort age + currentCohort%coage = currentCohort%coage + hlm_freq_day + write(fates_log(),*) 'cohort age: ',currentCohort%coage + if(currentCohort%coage < 0.0_r8)then + write(fates_log(),*) 'negative cohort age?',currentCohort%coage + end if + + ! update cohort age class and age x pft class call coagetype_class_index(currentCohort%coage, currentCohort%pft, & - currentCohort%coage_class,currentCohort%coage_by_pft_class) + currentCohort%coage_class,currentCohort%coage_by_pft_class) end if - + currentCohort => currentCohort%taller end do diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 7c186dde77..5dfd785af8 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -48,7 +48,10 @@ module EDParamsMod logical,protected, public :: active_crown_fire ! flag, 1=active crown fire 0=no active crown fire character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire" - + + logical,protected, public :: cohort_age_tracking ! flag, 1=track cohort age 0=no cohort age tracking + character(len=param_string_length),parameter :: fates_name_cohort_age_tracking = "fates_cohort_age_tracking" + real(r8), protected, public :: cg_strikes ! fraction of cloud to ground lightning strikes (0-1) character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes" @@ -359,6 +362,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=fates_name_active_crown_fire, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=fates_name_cohort_age_tracking, dimension_shape=dimension_shape_scalar, & + dimension_names=dim_names_scalar) + call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -379,7 +385,7 @@ subroutine FatesReceiveParams(fates_params) class(fates_parameters_type), intent(inout) :: fates_params real(r8) :: active_crown_fire_real !Local temp to transfer real data in file - + real(r8) :: cohort_age_tracking_real ! Local temp to transfer real data in file call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, & data=fates_mortality_disturbance_fraction) @@ -495,10 +501,14 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=fates_name_q10_froz, & data=q10_froz) - call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, & + call fates_params%RetreiveParameter(name=fates_name_active_crown_fire, & data=active_crown_fire_real) active_crown_fire = (abs(active_crown_fire_real-1.0_r8) Date: Mon, 20 Jan 2020 15:28:07 -0800 Subject: [PATCH 08/27] [Update param file with metadata and add checks on param values ] [Add metadata to the param file so that it is clear which mortality params act as flags to turn on age and size dependent mortality. Suggest sensible values for these parameters. Add checks so that age mortality cannot be on if cohort age tracking is off. Add checks so that negative inflection points are not allowed and that rate params must be specified. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing yet. --- main/EDPftvarcon.F90 | 66 +++++++++++++++++++++++- parameter_files/fates_params_default.cdl | 8 +-- 2 files changed, 69 insertions(+), 5 deletions(-) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 8c1e67cc1c..9903b25e11 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2079,7 +2079,7 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! cannot have a structural biomass allometry intercept of 0, and a non-woody ! plant (grass) can't have a non-zero intercept... ! ----------------------------------------------------------------------------------- - + use EDParamsMod , only : cohort_age_tracking ! Argument logical, intent(in) :: is_master ! Only log if this is the master proc @@ -2154,7 +2154,71 @@ subroutine FatesCheckParams(is_master, parteh_mode) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + ! Check to see if cohort tracking is on if age-senescence is on + !---------------------------------------------------------------------------------- + + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 100000.0_r8 .and. & + cohort_age_tracking .neqv. .TRUE. ) ) then + + write(fates_log(),*) 'Age-dependent mortality cannot be on' + write(fates_log(),*) 'if cohort age tracking is off' + write(fates_log(),*) 'Set cohort_age_tracking to 1' + write(fates_log(),*) 'To turn on cohort age tracking or ' + write(fates_log(),*) 'Set mort_ip_age_senescence to _ ' + write(fates_log(),*) 'to turn off age-dependent mortality ' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + ! Check that parameter ranges for age-dependent mortality make sense + !----------------------------------------------------------------------------------- + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 100000.0_r8 .and. & + EDPftvarcon_inst%mort_r_age_senescence(ipft) > 100000.0_r8 ) ) then + + write(fates_log(),*) 'Age-dependent mortality is on' + write(fates_log(),*) 'Please also set mort_r_age_senescence' + write(fates_log(),*) 'Sensible values are between 0.03-0.06' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Check that parameter ranges for age-dependent mortality make sense + !----------------------------------------------------------------------------------- + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 0.0_r8 ) ) then + + write(fates_log(),*) 'Negative mort_ip_age_senescence means' + write(fates_log(),*) 'trees will die at a negative age.' + write(fates_log(),*) 'Sensible values are between 10 and 2000?' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Check that parameter ranges for size-dependent mortality make sense + !----------------------------------------------------------------------------------- + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 100000.0_r8 .and. & + EDPftvarcon_inst%mort_r_size_senescence(ipft) > 100000.0_r8 ) ) then + + write(fates_log(),*) 'Size-dependent mortality is on' + write(fates_log(),*) 'Please also set mort_r_size_senescence' + write(fates_log(),*) 'Sensible values are between 0.03-0.06' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Check that parameter ranges for size-dependent mortality make sense + !----------------------------------------------------------------------------------- + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 0.0_r8 ) ) then + write(fates_log(),*) 'Negative mort_ip_size_senescence means' + write(fates_log(),*) 'trees will die at a negative dbh.' + write(fates_log(),*) 'Sensible values are between 1 and 300?' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + ! Check to see if mature and base seed allocation is greater than 1 ! ---------------------------------------------------------------------------------- diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index 309a7e8e78..cfe0b20a7f 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -312,16 +312,16 @@ variables: fates_mort_hf_sm_threshold:long_name = "soil moisture (btran units) at which drought mortality begins for non-hydraulic model" ; double fates_mort_ip_age_senescence(fates_pft) ; fates_mort_ip_age_senescence:units = "years" ; - fates_mort_ip_age_senescence:long_name = "Mortality cohort age senescence inflection point" ; + fates_mort_ip_age_senescence:long_name = "Mortality cohort age senescence inflection point. If _ this mortality term is off. Setting this value turns on age dependent mortality. " ; double fates_mort_ip_size_senescence(fates_pft) ; fates_mort_ip_size_senescence:units = "dbh cm" ; - fates_mort_ip_size_senescence:long_name = "Mortality dbh senescence inflection point" ; + fates_mort_ip_size_senescence:long_name = "Mortality dbh senescence inflection point. If _ this mortality term is off. Setting this value turns on size dependent mortality" ; double fates_mort_r_age_senescence(fates_pft) ; fates_mort_r_age_senescence:units = "mortality rate year^-1" ; - fates_mort_r_age_senescence:long_name = "Mortality age senescence rate of change" ; + fates_mort_r_age_senescence:long_name = "Mortality age senescence rate of change. Sensible range is around 0.03-0.06. Larger values givesteeper mortality curves." ; double fates_mort_r_size_senescence(fates_pft) ; fates_mort_r_size_senescence:units = "mortality rate dbh^-1" ; - fates_mort_r_size_senescence:long_name = "Mortality dbh senescence rate of change" ; + fates_mort_r_size_senescence:long_name = "Mortality dbh senescence rate of change. Sensible range is around 0.03-0.06. Larger values give steeper mortality curves." ; double fates_mort_scalar_coldstress(fates_pft) ; fates_mort_scalar_coldstress:units = "1/yr" ; fates_mort_scalar_coldstress:long_name = "maximum mortality rate from cold stress" ; From 9e5ab3259e97097108694a5ac7cf757c3efed87d Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Mon, 20 Jan 2020 21:47:20 -0800 Subject: [PATCH 09/27] [ Fix bug in flags for mortality params ] [Fix some bugs so that error messages work if the user attempts to run size or age dependent mortality without the correct parameter combinations. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Brief testing. --- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- main/EDMainMod.F90 | 2 -- main/EDPftvarcon.F90 | 23 +++++++++++++---------- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index a0249d59c5..05b90663ea 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -86,7 +86,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! rate (r) and inflection point (ip) define the increase in mortality rate with dbh mort_r_size_senescence = EDPftvarcon_inst%mort_r_size_senescence(cohort_in%pft) mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) - + ! if param values have been set then calculate smort if (mort_ip_size_senescence < 10000.0_r8 ) then smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 804c40986a..74b0b8f8ef 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -444,11 +444,9 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) end if ! if we are in age-dependent mortality mode - write(fates_log(),*) 'cohort_age_tracking',cohort_age_tracking if (cohort_age_tracking) then ! update cohort age currentCohort%coage = currentCohort%coage + hlm_freq_day - write(fates_log(),*) 'cohort age: ',currentCohort%coage if(currentCohort%coage < 0.0_r8)then write(fates_log(),*) 'negative cohort age?',currentCohort%coage end if diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 9903b25e11..4757d6328e 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2158,8 +2158,11 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check to see if cohort tracking is on if age-senescence is on !---------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 100000.0_r8 .and. & - cohort_age_tracking .neqv. .TRUE. ) ) then + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 10000.0_r8 ) .and. & + ( cohort_age_tracking .neqv. .TRUE. ) ) then + + write(fates_log(),*) 'ip_age = ', EDPftvarcon_inst%mort_ip_age_senescence(ipft) + write(fates_log(),*) 'cohort_age_tracking = ', cohort_age_tracking write(fates_log(),*) 'Age-dependent mortality cannot be on' write(fates_log(),*) 'if cohort age tracking is off' @@ -2174,8 +2177,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for age-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 100000.0_r8 .and. & - EDPftvarcon_inst%mort_r_age_senescence(ipft) > 100000.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 10000.0_r8 ) .and. & + ( EDPftvarcon_inst%mort_r_age_senescence(ipft) > 10000.0_r8 ) ) then write(fates_log(),*) 'Age-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_age_senescence' @@ -2197,8 +2200,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for size-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 100000.0_r8 .and. & - EDPftvarcon_inst%mort_r_size_senescence(ipft) > 100000.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 10000.0_r8 ) .and. & + ( EDPftvarcon_inst%mort_r_size_senescence(ipft) > 10000.0_r8 ) ) then write(fates_log(),*) 'Size-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_size_senescence' @@ -2217,13 +2220,13 @@ subroutine FatesCheckParams(is_master, parteh_mode) write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if - - + + ! Check to see if mature and base seed allocation is greater than 1 ! ---------------------------------------------------------------------------------- if ( ( EDPftvarcon_inst%seed_alloc(ipft) + & - EDPftvarcon_inst%seed_alloc_mature(ipft)) > 1.0_r8 ) then + EDPftvarcon_inst%seed_alloc_mature(ipft)) > 1.0_r8 ) then write(fates_log(),*) 'The sum of seed allocation from base and mature trees may' write(fates_log(),*) ' not exceed 1.' @@ -2232,7 +2235,7 @@ subroutine FatesCheckParams(is_master, parteh_mode) write(fates_log(),*) ' seed_alloc_mature: ',EDPftvarcon_inst%seed_alloc_mature(ipft) write(fates_log(),*) ' Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) - + end if ! Check if woody plants have a structural biomass (agb) intercept From b01c1363872b01023d5e6c9781752ee355bbb370 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Tue, 21 Jan 2020 11:13:59 -0800 Subject: [PATCH 10/27] [In EDCohortDynamicsMod increase maxCohortsPerPatch to 300 if cohort age tracking is on ] [If cohort age tracking is on then cohorts can only fuse if size and age are similar enough. This means we end up with many more cohorts and maxCohortsPerPatch needs to be increased. I added a local variable maxCohortsPerPatch_age_tracking to EDCohortDynamicsMod.F90 that is used in place of maxCohortsPerPatch if age tracking is on. Not sure this is the best way to go about it but it seems to be working. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Quick test on lobata. It compiles and runs. --- biogeochem/EDCohortDynamicsMod.F90 | 52 +++++++++++++++++++++++------- 1 file changed, 40 insertions(+), 12 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index c69eacf85e..419dc0e919 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -944,6 +944,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) use EDParamsMod , only : cohort_age_tracking use FatesConstantsMod , only : itrue use FatesConstantsMod, only : days_per_year + use EDTypesMod , only : maxCohortsPerPatch ! ! !ARGUMENTS @@ -971,7 +972,8 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) real(r8) :: leaf_c_curr ! Leaf carbon * plant density of next (for weighting) real(r8) :: leaf_c_target real(r8) :: dynamic_size_fusion_tolerance - real(r8) :: dynamic_age_fusion_tolerance + real(r8) :: dynamic_age_fusion_tolerance + integer :: maxCohortsPerPatch_age_tracking real(r8) :: dbh real(r8) :: leaf_c ! leaf carbon [kg] @@ -990,6 +992,12 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! set the cohort age fusion tolerance (this is in days - remains constant) dynamic_age_fusion_tolerance = ED_val_cohort_age_fusion_tol + if ( cohort_age_tracking ) then + maxCohortsPerPatch_age_tracking = 300 + end if + + + !This needs to be a function of the canopy layer, because otherwise, at canopy closure !the number of cohorts doubles and very dissimilar cohorts are fused together !because c_area and biomass are non-linear with dbh, this causes several mass inconsistancies @@ -999,6 +1007,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) iterate = 1 fusion_took_place = 0 + !---------------------------------------------------------------------! ! Keep doing this until nocohorts <= maxcohorts ! !---------------------------------------------------------------------! @@ -1404,20 +1413,39 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) currentCohort => currentCohort%shorter enddo - if (nocohorts > maxCohortsPerPatch) then - iterate = 1 - !---------------------------------------------------------------------! - ! Making profile tolerance larger means that more fusion will happen ! - !---------------------------------------------------------------------! - dynamic_size_fusion_tolerance = dynamic_size_fusion_tolerance * 1.1_r8 - dynamic_age_fusion_tolerance = dynamic_age_fusion_tolerance * 1.1_r8 - !write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance - else + if ( cohort_age_tracking ) then + if ( nocohorts > maxCohortsPerPatch_age_tracking ) then + iterate = 1 + !---------------------------------------------------------------------! + ! Making profile tolerance larger means that more fusion will happen ! + !---------------------------------------------------------------------! + dynamic_size_fusion_tolerance = dynamic_size_fusion_tolerance * 1.1_r8 + dynamic_age_fusion_tolerance = dynamic_age_fusion_tolerance * 1.1_r8 + !write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance - iterate = 0 - endif + else + + iterate = 0 + endif + else + + if (nocohorts > maxCohortsPerPatch) then + iterate = 1 + !---------------------------------------------------------------------! + ! Making profile tolerance larger means that more fusion will happen ! + !---------------------------------------------------------------------! + dynamic_size_fusion_tolerance = dynamic_size_fusion_tolerance * 1.1_r8 + !write(fates_log(),*) 'maxcohorts exceeded',dynamic_fusion_tolerance + + else + + iterate = 0 + endif + end if + + if ( dynamic_size_fusion_tolerance .gt. 100._r8) then ! something has gone terribly wrong and we need to report what write(fates_log(),*) 'exceeded reasonable expectation of cohort fusion.' From e45153c7271edcef6c4edce82976c7a80e5eb040 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Tue, 21 Jan 2020 13:48:45 -0800 Subject: [PATCH 11/27] [ More detailed flags for mortality parameter checks in EDPftvarcon ] [Check all size and age mortality parameters are positive and turned on in the correct combinations. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Appears to work as expected on lobata. --- main/EDPftvarcon.F90 | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 4757d6328e..9834f38751 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2189,11 +2189,13 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for age-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 0.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 0.0_r8 ) .or. & + (EDPftvarcon_inst%mort_r_age_senescence(ipft) < 0.0_r8 ) ) then - write(fates_log(),*) 'Negative mort_ip_age_senescence means' - write(fates_log(),*) 'trees will die at a negative age.' - write(fates_log(),*) 'Sensible values are between 10 and 2000?' + write(fates_log(),*) 'Either mort_ip_age_senescence or mort_r_age_senescece' + write(fates_log(),*) 'is negative which makes no biological sense.' + write(fates_log(),*) 'Sensible values for ip are between 10 and 2000?' + write(fates_log(),*) 'Sensible values for r are between 0.03-0.06' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -2212,11 +2214,13 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for size-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 0.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 0.0_r8 ) .or. & + ( EDPftvarcon_inst%mort_r_size_senescence(ipft) < 0.0_r8 ) ) then - write(fates_log(),*) 'Negative mort_ip_size_senescence means' - write(fates_log(),*) 'trees will die at a negative dbh.' - write(fates_log(),*) 'Sensible values are between 1 and 300?' + write(fates_log(),*) 'Either mort_ip_size_senescence or mort_r_size_senescence' + write(fates_log(),*) 'is negative which makes no biological sense.' + write(fates_log(),*) 'Sensible values for ip are between 1 and 300?' + write(fates_log(),*) 'Sensible values for r are between 0.03-0.06' write(fates_log(),*) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if From fb52c507da612490a3a5ea53dd820f1d44599e40 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Wed, 22 Jan 2020 15:27:42 -0800 Subject: [PATCH 12/27] [ Change the cap on mortality ] [ Change the cap on mortality so that if annual mortality of both background and size or age dependent mortality gets above 1 background mortality is reduced. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes (describe what changes), No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: 200 year run on Lawrencium worked with both size and age dependent mortality. --- biogeochem/EDMortalityFunctionsMod.F90 | 16 ++++++---------- 1 file changed, 6 insertions(+), 10 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 05b90663ea..fc547683c4 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -250,11 +250,9 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) currentCohort%lmort_collateral + & currentCohort%lmort_infra)/hlm_freq_day - ! this caps bmort so that daily mortality cannot exceed 0.995 - if ((cmort+hmort+bmort+frmort+smort+asmort + dndt_logging)*hlm_freq_day & - > 1.0_r8)then - bmort = (1.0_r8 - ((cmort+hmort+asmort+frmort+smort+dndt_logging)*hlm_freq_day)) & - /hlm_freq_day + ! this caps bmort so that daily mortality cannot exceed 1 + if (( bmort + smort + asmort ) > 1.0_r8 ) then + bmort = (1.0_r8 - (asmort + smort )) endif currentCohort%dndt = -1.0_r8 * & @@ -263,11 +261,9 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) else ! cap on bmort for canopy layers - bmort is adjusted so that total mortality - ! including disturbance fraction does not exceed 0.995 - if(((cmort+hmort+bmort+frmort+smort+asmort)*hlm_freq_day) > & - 1.0_r8 - fates_mortality_disturbance_fraction)then - bmort = (1.0_r8 - fates_mortality_disturbance_fraction - & - ((cmort+hmort+asmort+frmort+smort)*hlm_freq_day))/hlm_freq_day + ! including disturbance fraction does not exceed 1 + if (( bmort + smort + asmort) > 1.0_r8 - fates_mortality_disturbance_fraction ) then + bmort = (1.0_r8 - fates_mortality_disturbance_fraction - ( asmort + smort) ) endif From ab6fc2008175f6da0354375ed47a0eab77712006 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Thu, 23 Jan 2020 12:21:37 -0800 Subject: [PATCH 13/27] [ Remove cap on mortality ] [ This cap was not actually coming into effect. ] Fixes: [NGT-ED Github issue #] User interface changes?: No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: No testing --- biogeochem/EDMortalityFunctionsMod.F90 | 9 --------- 1 file changed, 9 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index fc547683c4..6e8acdbee0 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -250,21 +250,12 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) currentCohort%lmort_collateral + & currentCohort%lmort_infra)/hlm_freq_day - ! this caps bmort so that daily mortality cannot exceed 1 - if (( bmort + smort + asmort ) > 1.0_r8 ) then - bmort = (1.0_r8 - (asmort + smort )) - endif currentCohort%dndt = -1.0_r8 * & (cmort+hmort+bmort+frmort+smort+asmort + dndt_logging) & * currentCohort%n else - ! cap on bmort for canopy layers - bmort is adjusted so that total mortality - ! including disturbance fraction does not exceed 1 - if (( bmort + smort + asmort) > 1.0_r8 - fates_mortality_disturbance_fraction ) then - bmort = (1.0_r8 - fates_mortality_disturbance_fraction - ( asmort + smort) ) - endif ! Mortality from logging in the canopy is ONLY disturbance generating, don't From dd28a9b10476e56126976bc6a8074ec88bd7aa06 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Wed, 29 Jan 2020 17:20:29 -0800 Subject: [PATCH 14/27] [Change check on initialised params to use the check_initialized that is defined in PRTGenericMod ] [PRTGenericMod contains a parameter that is used to check if parameters have been initialised i.e. are not _ in the param file. Use this variable for places where we check if mortality parameters are set. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:no testing --- biogeochem/EDMortalityFunctionsMod.F90 | 7 ++++--- main/EDPftvarcon.F90 | 13 +++++++------ main/FatesInterfaceMod.F90 | 4 ++-- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 6e8acdbee0..69b1af37b9 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -51,7 +51,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesInterfaceMod , only : hlm_hio_ignore_val - + use PRTGenericMod , only : check_initialized + type (ed_cohort_type), intent(in) :: cohort_in type (bc_in_type), intent(in) :: bc_in real(r8),intent(out) :: bmort ! background mortality : Fraction per year @@ -88,7 +89,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) ! if param values have been set then calculate smort - if (mort_ip_size_senescence < 10000.0_r8 ) then + if ( mort_ip_size_senescence < check_initialized ) then smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & (cohort_in%dbh - mort_ip_size_senescence) ) ) else @@ -97,7 +98,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! if user has set cohort age fusion param to not be _ we calculate age ! dependent mortality - if (mort_ip_age_senescence < 10000.0_r8) then + if ( mort_ip_age_senescence < check_initialized ) then ! Age Dependent Senescence ! rate and inflection point define the change in mortality with age mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index 9834f38751..a496091766 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2079,7 +2079,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! cannot have a structural biomass allometry intercept of 0, and a non-woody ! plant (grass) can't have a non-zero intercept... ! ----------------------------------------------------------------------------------- - use EDParamsMod , only : cohort_age_tracking + use EDParamsMod , only : cohort_age_tracking + use PRTGenericMod , only : check_initialized ! Argument logical, intent(in) :: is_master ! Only log if this is the master proc @@ -2158,7 +2159,7 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check to see if cohort tracking is on if age-senescence is on !---------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 10000.0_r8 ) .and. & + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < check_initialized ) .and. & ( cohort_age_tracking .neqv. .TRUE. ) ) then write(fates_log(),*) 'ip_age = ', EDPftvarcon_inst%mort_ip_age_senescence(ipft) @@ -2177,8 +2178,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for age-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 10000.0_r8 ) .and. & - ( EDPftvarcon_inst%mort_r_age_senescence(ipft) > 10000.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < check_initialized ) .and. & + ( EDPftvarcon_inst%mort_r_age_senescence(ipft) > check_initialized ) ) then write(fates_log(),*) 'Age-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_age_senescence' @@ -2202,8 +2203,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for size-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 10000.0_r8 ) .and. & - ( EDPftvarcon_inst%mort_r_size_senescence(ipft) > 10000.0_r8 ) ) then + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < check_initialized ) .and. & + ( EDPftvarcon_inst%mort_r_size_senescence(ipft) > check_initialized ) ) then write(fates_log(),*) 'Size-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_size_senescence' diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b971d11742..b92f0cb26a 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1541,10 +1541,10 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( hlm_use_inventory_init.eq.1 .and. ED_val_cohort_age_fusion_tol < 10000.0_r8) then + if ( hlm_use_inventory_init.eq.1 .and. cohort_age_tracking .eqv. 1) then if (fates_global_verbose()) then write(fates_log(), *) 'Fates inventory init cannot be used with age dependent mortality' - write(fates_log(), *) 'Set fates_cohort_age_fusion_tol to _ or turn off inventory init' + write(fates_log(), *) 'Set cohort_age_tracking to 0 or turn off inventory init' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if From 63c2958706d15e9b263fc0fe8edd78c46390aee7 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Wed, 11 Mar 2020 10:58:44 -0700 Subject: [PATCH 15/27] [ Fix some minor bugs/make changes suggested by Rosie Fisher ] [Change units on cohort age to days in EDCohortDynamics, add units in comments in other places in the code, add a missing _r8. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: no testing. --- biogeochem/EDCohortDynamicsMod.F90 | 10 +++++----- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 419dc0e919..b8c5fc016a 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -490,7 +490,7 @@ subroutine nan_cohort(cc_p) currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm - currentCohort%coage = nan ! age of the cohort in years + currentCohort%coage = nan ! age of the cohort in days currentCohort%hite = nan ! height: meters currentCohort%laimemory = nan ! target leaf biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 @@ -987,9 +987,9 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) !---------------------------------------------------------------------- - !set initial fusion tolerance + !set initial fusion tolerance (in cm) dynamic_size_fusion_tolerance = ED_val_cohort_size_fusion_tol - ! set the cohort age fusion tolerance (this is in days - remains constant) + ! set the cohort age fusion tolerance (in days) dynamic_age_fusion_tolerance = ED_val_cohort_age_fusion_tol if ( cohort_age_tracking ) then @@ -1028,7 +1028,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) do while (associated(nextc)) nextnextc => nextc%shorter - diff = abs((currentCohort%dbh - nextc%dbh)/(0.5*(currentCohort%dbh + nextc%dbh))) + diff = abs((currentCohort%dbh - nextc%dbh)/(0.5_r8*(currentCohort%dbh + nextc%dbh))) !Criteria used to divide up the height continuum into different cohorts. @@ -1042,7 +1042,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) coage_diff = 0.0_r8 else coage_diff = abs((currentCohort%coage - nextc%coage)/ & - (0.5*(currentCohort%coage + nextc%coage))) + (0.5_r8*(currentCohort%coage + nextc%coage))) end if if (coage_diff <= dynamic_age_fusion_tolerance ) then diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 69b1af37b9..f25d37b11a 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -46,7 +46,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor ! ============================================================================ ! Calculate mortality rates from carbon storage, hydraulic cavitation, - ! background and freezing and size dependent senescence + ! background and freezing and size and age dependent senescence ! ============================================================================ use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm From 774b105f5aa94cd01d2213abb2e866b0a3ee936a Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Wed, 11 Mar 2020 12:44:23 -0700 Subject: [PATCH 16/27] [ Cohort age units to years ] [Units for cohort age back to years ] Fixes: [NGT-ED Github issue #] User interface changes?: [ No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing. --- biogeochem/EDCohortDynamicsMod.F90 | 6 +++--- biogeochem/EDMortalityFunctionsMod.F90 | 3 +-- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index b8c5fc016a..f07071d827 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -163,7 +163,7 @@ subroutine create_cohort(currentSite, patchptr, pft, nn, hite, coage, dbh, & real(r8), intent(in) :: nn ! number of individuals in cohort ! per 'area' (10000m2 default) real(r8), intent(in) :: hite ! height: meters - real(r8), intent(in) :: coage ! cohort age in days + real(r8), intent(in) :: coage ! cohort age in years real(r8), intent(in) :: dbh ! dbh: cm class(prt_vartypes),target :: prt ! The allocated PARTEH ! object @@ -490,7 +490,7 @@ subroutine nan_cohort(cc_p) currentCohort%n = nan ! number of individuals in cohort per 'area' (10000m2 default) currentCohort%dbh = nan ! 'diameter at breast height' in cm - currentCohort%coage = nan ! age of the cohort in days + currentCohort%coage = nan ! age of the cohort in years currentCohort%hite = nan ! height: meters currentCohort%laimemory = nan ! target leaf biomass- set from previous year: kGC per indiv currentCohort%lai = nan ! leaf area index of cohort m2/m2 @@ -989,7 +989,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) !set initial fusion tolerance (in cm) dynamic_size_fusion_tolerance = ED_val_cohort_size_fusion_tol - ! set the cohort age fusion tolerance (in days) + ! set the cohort age fusion tolerance (in fraction of years) dynamic_age_fusion_tolerance = ED_val_cohort_age_fusion_tol if ( cohort_age_tracking ) then diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index f25d37b11a..70fa8de984 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -96,8 +96,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor smort = 0.0_r8 end if - ! if user has set cohort age fusion param to not be _ we calculate age - ! dependent mortality + ! if param values have been set then calculate asmort if ( mort_ip_age_senescence < check_initialized ) then ! Age Dependent Senescence ! rate and inflection point define the change in mortality with age From c83725ab1466f4d62f643ba3f5bc37b628c40014 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Wed, 11 Mar 2020 13:59:05 -0700 Subject: [PATCH 17/27] [ Bug fix - .eqv. to .TRUE. for cohort_age_tracking check ] [In checks that age mortality and cohort age tracking are compatible the check on cohort age tracking needs to be a logical .TRUE. or .FALSE. Fixed bug where it was .eqv. 0 ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing. --- main/FatesInterfaceMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b92f0cb26a..ff73b1d5de 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1428,8 +1428,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) ! ! RGK-2016 ! --------------------------------------------------------------------------------- - use EDParamsMod , only : ED_val_cohort_age_fusion_tol - + use EDParamsMod , only : cohort_age_tracking + ! Arguments integer, optional, intent(in) :: ival real(r8), optional, intent(in) :: rval @@ -1541,7 +1541,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( hlm_use_inventory_init.eq.1 .and. cohort_age_tracking .eqv. 1) then + if ( hlm_use_inventory_init.eq.1 .and. cohort_age_tracking .eqv. .TRUE.) then if (fates_global_verbose()) then write(fates_log(), *) 'Fates inventory init cannot be used with age dependent mortality' write(fates_log(), *) 'Set cohort_age_tracking to 0 or turn off inventory init' From aef198f7d7f81a53398e4f8eb16c3f338f75ba08 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Wed, 11 Mar 2020 19:59:34 -0700 Subject: [PATCH 18/27] [Remove the cohort age tracking parameter and make cohort age tracking a 'hlm_use_'flag. ] [Make cohort age tracking an option that is turned on or off with a flag 'hlm_use_cohort_age_tracking'. Remove the parameter 'cohort_age_tracking' that used to serve this purpose. Also created a variable fates_check_param_set which we can compare mortality parameters against to see if the user has left them '_' meaning size and age mortality is off, or set them to real numbers meaning size or age mortality is on. ] Fixes: [NGT-ED Github issue #] User interface changes?: [Yes - one less parameter, but a new namelist option.] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: Brief testing on lobata. Seems to work as expected. --- biogeochem/EDCanopyStructureMod.F90 | 6 ++-- biogeochem/EDCohortDynamicsMod.F90 | 10 +++--- biogeochem/EDMortalityFunctionsMod.F90 | 6 ++-- main/EDMainMod.F90 | 6 ++-- main/EDParamsMod.F90 | 12 ------- main/EDPftvarcon.F90 | 44 +++++++----------------- main/FatesConstantsMod.F90 | 3 ++ main/FatesInterfaceMod.F90 | 9 +++-- parameter_files/fates_params_default.cdl | 5 --- 9 files changed, 36 insertions(+), 65 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 5b538d29b5..acf7a9edd0 100644 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -26,10 +26,10 @@ module EDCanopyStructureMod use FatesGlobals , only : endrun => fates_endrun use FatesInterfaceMod , only : hlm_days_per_year use FatesInterfaceMod , only : hlm_use_planthydro + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use FatesInterfaceMod , only : numpft use FatesPlantHydraulicsMod, only : UpdateH2OVeg,InitHydrCohort, RecruitWaterStorage use EDTypesMod , only : maxCohortsPerPatch - use EDParamsMod , only : cohort_age_tracking use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : all_carbon_elements @@ -1257,12 +1257,12 @@ subroutine canopy_summarization( nsites, sites, bc_in ) ! --------------------------------------------------------------------------------- use FatesInterfaceMod , only : bc_in_type + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use EDPatchDynamicsMod , only : set_patchno use FatesSizeAgeTypeIndicesMod, only : sizetype_class_index use FatesSizeAgeTypeIndicesMod, only : coagetype_class_index use EDtypesMod , only : area use EDPftvarcon , only : EDPftvarcon_inst - use EDParamsMod , only : cohort_age_tracking use FatesConstantsMod , only : itrue ! !ARGUMENTS @@ -1326,7 +1326,7 @@ subroutine canopy_summarization( nsites, sites, bc_in ) call sizetype_class_index(currentCohort%dbh,currentCohort%pft, & currentCohort%size_class,currentCohort%size_by_pft_class) - if (cohort_age_tracking) then + if (hlm_use_cohort_age_tracking .eq. itrue) then call coagetype_class_index(currentCohort%coage,currentCohort%pft, & currentCohort%coage_class,currentCohort%coage_by_pft_class) end if diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index f07071d827..d75194ce99 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -9,6 +9,7 @@ module EDCohortDynamicsMod use FatesInterfaceMod , only : hlm_freq_day use FatesInterfaceMod , only : bc_in_type use FatesInterfaceMod , only : hlm_use_planthydro + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : fates_unset_int use FatesConstantsMod , only : itrue,ifalse @@ -37,7 +38,6 @@ module EDCohortDynamicsMod use EDTypesMod , only : site_fluxdiags_type use EDTypesMod , only : num_elements use EDParamsMod , only : ED_val_cohort_age_fusion_tol - use EDParamsMod , only : cohort_age_tracking use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_parteh_mode use FatesPlantHydraulicsMod, only : FuseCohortHydraulics @@ -941,7 +941,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! !USES: use EDParamsMod , only : ED_val_cohort_size_fusion_tol use EDParamsMod , only : ED_val_cohort_age_fusion_tol - use EDParamsMod , only : cohort_age_tracking + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use FatesConstantsMod , only : itrue use FatesConstantsMod, only : days_per_year use EDTypesMod , only : maxCohortsPerPatch @@ -992,7 +992,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) ! set the cohort age fusion tolerance (in fraction of years) dynamic_age_fusion_tolerance = ED_val_cohort_age_fusion_tol - if ( cohort_age_tracking ) then + if ( hlm_use_cohort_age_tracking .eq. itrue) then maxCohortsPerPatch_age_tracking = 300 end if @@ -1094,7 +1094,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) (nextc%coage * (nextc%n/(currentCohort%n + nextc%n))) ! update the cohort age again - if (cohort_age_tracking) then + if (hlm_use_cohort_age_tracking .eq.itrue) then call coagetype_class_index(currentCohort%coage, currentCohort%pft, & currentCohort%coage_class, currentCohort%coage_by_pft_class) end if @@ -1414,7 +1414,7 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) enddo - if ( cohort_age_tracking ) then + if ( hlm_use_cohort_age_tracking .eq.itrue) then if ( nocohorts > maxCohortsPerPatch_age_tracking ) then iterate = 1 !---------------------------------------------------------------------! diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 70fa8de984..71b2b2b29c 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -51,7 +51,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesInterfaceMod , only : hlm_hio_ignore_val - use PRTGenericMod , only : check_initialized + use FatesConstantsMod, only : fates_check_param_set type (ed_cohort_type), intent(in) :: cohort_in type (bc_in_type), intent(in) :: bc_in @@ -89,7 +89,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor mort_ip_size_senescence = EDPftvarcon_inst%mort_ip_size_senescence(cohort_in%pft) ! if param values have been set then calculate smort - if ( mort_ip_size_senescence < check_initialized ) then + if ( mort_ip_size_senescence < fates_check_param_set ) then smort = 1.0_r8 / ( 1.0_r8 + exp( -1.0_r8 * mort_r_size_senescence * & (cohort_in%dbh - mort_ip_size_senescence) ) ) else @@ -97,7 +97,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor end if ! if param values have been set then calculate asmort - if ( mort_ip_age_senescence < check_initialized ) then + if ( mort_ip_age_senescence < fates_check_param_set ) then ! Age Dependent Senescence ! rate and inflection point define the change in mortality with age mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 74b0b8f8ef..4641aec3e8 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -15,6 +15,7 @@ module EDMainMod use FatesInterfaceMod , only : hlm_current_month use FatesInterfaceMod , only : hlm_current_day use FatesInterfaceMod , only : hlm_use_planthydro + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use FatesInterfaceMod , only : hlm_reference_date use FatesInterfaceMod , only : hlm_use_ed_prescribed_phys use FatesInterfaceMod , only : hlm_use_ed_st3 @@ -55,7 +56,6 @@ module EDMainMod use EDTypesMod , only : element_pos use EDTypesMod , only : phen_dstat_moiston use EDTypesMod , only : phen_dstat_timeon - use EDParamsMod , only : cohort_age_tracking use FatesConstantsMod , only : itrue,ifalse use FatesConstantsMod , only : primaryforest, secondaryforest use FatesConstantsMod , only : nearzero @@ -284,7 +284,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! FIX(SPM,032414) refactor so everything goes through interface ! ! !USES: - use EDParamsMod, only : cohort_age_tracking + use FatesInterfaceMod, only : hlm_use_cohort_age_tracking use FatesConstantsMod, only : itrue ! !ARGUMENTS: @@ -444,7 +444,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) end if ! if we are in age-dependent mortality mode - if (cohort_age_tracking) then + if (hlm_use_cohort_age_tracking .eq. itrue) then ! update cohort age currentCohort%coage = currentCohort%coage + hlm_freq_day if(currentCohort%coage < 0.0_r8)then diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 5dfd785af8..d9c7a4e050 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -49,9 +49,6 @@ module EDParamsMod logical,protected, public :: active_crown_fire ! flag, 1=active crown fire 0=no active crown fire character(len=param_string_length),parameter :: fates_name_active_crown_fire = "fates_fire_active_crown_fire" - logical,protected, public :: cohort_age_tracking ! flag, 1=track cohort age 0=no cohort age tracking - character(len=param_string_length),parameter :: fates_name_cohort_age_tracking = "fates_cohort_age_tracking" - real(r8), protected, public :: cg_strikes ! fraction of cloud to ground lightning strikes (0-1) character(len=param_string_length),parameter :: fates_name_cg_strikes="fates_fire_cg_strikes" @@ -361,9 +358,6 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=fates_name_active_crown_fire, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) - - call fates_params%RegisterParameter(name=fates_name_cohort_age_tracking, dimension_shape=dimension_shape_scalar, & - dimension_names=dim_names_scalar) call fates_params%RegisterParameter(name=fates_name_cg_strikes, dimension_shape=dimension_shape_scalar, & dimension_names=dim_names_scalar) @@ -385,7 +379,6 @@ subroutine FatesReceiveParams(fates_params) class(fates_parameters_type), intent(inout) :: fates_params real(r8) :: active_crown_fire_real !Local temp to transfer real data in file - real(r8) :: cohort_age_tracking_real ! Local temp to transfer real data in file call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, & data=fates_mortality_disturbance_fraction) @@ -505,10 +498,6 @@ subroutine FatesReceiveParams(fates_params) data=active_crown_fire_real) active_crown_fire = (abs(active_crown_fire_real-1.0_r8) check_initialized ) ) then + if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < fates_check_param_set ) .and. & + ( EDPftvarcon_inst%mort_r_age_senescence(ipft) > fates_check_param_set ) ) then write(fates_log(),*) 'Age-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_age_senescence' @@ -2203,8 +2185,8 @@ subroutine FatesCheckParams(is_master, parteh_mode) ! Check that parameter ranges for size-dependent mortality make sense !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < check_initialized ) .and. & - ( EDPftvarcon_inst%mort_r_size_senescence(ipft) > check_initialized ) ) then + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < fates_check_param_set ) .and. & + ( EDPftvarcon_inst%mort_r_size_senescence(ipft) > fates_check_param_set ) ) then write(fates_log(),*) 'Size-dependent mortality is on' write(fates_log(),*) 'Please also set mort_r_size_senescence' diff --git a/main/FatesConstantsMod.F90 b/main/FatesConstantsMod.F90 index 46d6aa8ea6..a53a1d3e33 100644 --- a/main/FatesConstantsMod.F90 +++ b/main/FatesConstantsMod.F90 @@ -21,6 +21,9 @@ module FatesConstantsMod ! Used to initialize and test unset r8s real(fates_r8), parameter, public :: fates_unset_r8 = -1.e36_fates_r8 + ! Used to check if a parameter was specified in the parameter file (or left as _) + real(fates_r8), parameter, public :: fates_check_param_set = 9.9e32_fates_r8 + ! Integer equivalent of true (in case some compilers dont auto convert) integer, parameter, public :: itrue = 1 diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index ff73b1d5de..2605a4d2c3 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -138,6 +138,9 @@ module FatesInterfaceMod ! 1 = TRUE, 0 = FALSE ! THIS IS CURRENTLY NOT SUPPORTED + integer, public, protected :: hlm_use_cohort_age_tracking ! This flag signals whether or not to use + ! cohort age tracking. 1 = TRUE, 0 = FALSE + integer, public, protected :: hlm_use_ed_st3 ! This flag signals whether or not to use ! (ST)atic (ST)and (ST)ructure mode (ST3) ! Essentially, this gives us the ability @@ -1428,7 +1431,6 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) ! ! RGK-2016 ! --------------------------------------------------------------------------------- - use EDParamsMod , only : cohort_age_tracking ! Arguments integer, optional, intent(in) :: ival @@ -1462,6 +1464,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_parteh_mode = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int + hlm_use_cohort_age_tracking = unset_int hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int @@ -1541,10 +1544,10 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - if ( hlm_use_inventory_init.eq.1 .and. cohort_age_tracking .eqv. .TRUE.) then + if ( hlm_use_inventory_init.eq.1 .and. hlm_use_cohort_age_tracking .eq.1) then if (fates_global_verbose()) then write(fates_log(), *) 'Fates inventory init cannot be used with age dependent mortality' - write(fates_log(), *) 'Set cohort_age_tracking to 0 or turn off inventory init' + write(fates_log(), *) 'Set hlm_use_cohort_age_tracking to 0 or turn off inventory init' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index cfe0b20a7f..45f8721b69 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -529,9 +529,6 @@ variables: double fates_canopy_closure_thresh ; fates_canopy_closure_thresh:units = "unitless" ; fates_canopy_closure_thresh:long_name = "tree canopy coverage at which crown area allometry changes from savanna to forest value" ; - double fates_cohort_age_tracking ; - fates_cohort_age_tracking:units = "0 or 1" ; - fates_cohort_age_tracking:long_name = "flag, 1=track cohort age, 0=no cohort age tracking" ; double fates_cohort_age_fusion_tol ; fates_cohort_age_fusion_tol:units = "unitless" ; fates_cohort_age_fusion_tol:long_name = "minimum fraction in differece in cohort age between cohorts." ; @@ -1229,8 +1226,6 @@ data: fates_canopy_closure_thresh = 0.8 ; - fates_cohort_age_tracking = 0 ; - fates_cohort_age_fusion_tol = _ ; fates_cohort_size_fusion_tol = 0.08 ; From 9345e67b890bc35a95f27ebd2e7a3206ca05cf40 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Thu, 12 Mar 2020 13:27:33 -0700 Subject: [PATCH 19/27] [Fix bug with hlm_use_cohort_age_tracking flag ] [After being unset make sure flag hlm_use_cohort_age_tracking gets given a 0 or 1. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing. --- main/FatesInterfaceMod.F90 | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 2605a4d2c3..f6cca89ada 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1522,6 +1522,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + if ( .not.((hlm_use_ed_st3.eq.1).or.(hlm_use_ed_st3.eq.0)) ) then if (fates_global_verbose()) then write(fates_log(), *) 'The FATES namelist stand structure flag must be 0 or 1, exiting' @@ -1655,6 +1657,14 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + if(hlm_use_cohort_age_tracking .eq. unset_int) then + if (fates_global_verbose()) then + write(fates_log(), *) 'switch for cohort_age_tracking unset: hlm_use_cohort_age_tracking, exiting' + end if + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + if (fates_global_verbose()) then write(fates_log(), *) 'Checked. All control parameters sent to FATES.' end if @@ -1737,6 +1747,13 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' end if + case('use_cohort_age_tracking') + hlm_use_cohort_age_tracking = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering hlm_use_cohort_age_tracking= ',ival,' to FATES' + end if + + case('use_logging') hlm_use_logging = ival if (fates_global_verbose()) then From 54e8dfb7c27133fdbc261f91f89de08cb5ba71d3 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Tue, 17 Mar 2020 14:28:11 -0700 Subject: [PATCH 20/27] [Add one additional test to ensure age-dependent mortality is not on without cohort age tracking ] [Add a test in FatesInterfaceMod to ensuer that if age-dependent mortality parameters are set then hlm_use_cohort_age_tracking is also set to true. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:Complies and then crashes with the expected error message. --- main/EDPftvarcon.F90 | 26 +++++++++++++------------- main/FatesInterfaceMod.F90 | 17 +++++++++++++++++ 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index c00383c07c..b4654f2c13 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -2178,19 +2178,6 @@ subroutine FatesCheckParams(is_master, parteh_mode) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - ! Check that parameter ranges for age-dependent mortality make sense - !----------------------------------------------------------------------------------- - if ( ( EDPftvarcon_inst%mort_ip_age_senescence(ipft) < 0.0_r8 ) .or. & - (EDPftvarcon_inst%mort_r_age_senescence(ipft) < 0.0_r8 ) ) then - - write(fates_log(),*) 'Either mort_ip_age_senescence or mort_r_age_senescece' - write(fates_log(),*) 'is negative which makes no biological sense.' - write(fates_log(),*) 'Sensible values for ip are between 10 and 2000?' - write(fates_log(),*) 'Sensible values for r are between 0.03-0.06' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - ! Check that parameter ranges for size-dependent mortality make sense !----------------------------------------------------------------------------------- if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < fates_check_param_set ) .and. & @@ -2217,6 +2204,19 @@ subroutine FatesCheckParams(is_master, parteh_mode) end if + ! Check that parameter ranges for size-dependent mortality make sense + !----------------------------------------------------------------------------------- + if ( ( EDPftvarcon_inst%mort_ip_size_senescence(ipft) < 0.0_r8 ) .or. & + ( EDPftvarcon_inst%mort_r_size_senescence(ipft) < 0.0_r8 ) ) then + + write(fates_log(),*) 'Either mort_ip_size_senescence or mort_r_size_senescence' + write(fates_log(),*) 'is negative which makes no biological sense.' + write(fates_log(),*) 'Sensible values for ip are between 1 and 300?' + write(fates_log(),*) 'Sensible values for r are between 0.03-0.06' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + ! Check to see if mature and base seed allocation is greater than 1 ! ---------------------------------------------------------------------------------- diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index f6cca89ada..2a6aa169ed 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1046,6 +1046,7 @@ subroutine set_fates_global_elements(use_fates) use EDParamsMod, only : ED_val_history_sizeclass_bin_edges, ED_val_history_ageclass_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges use EDParamsMod, only : ED_val_history_coageclass_bin_edges + use FatesConstantsMod, only : fates_check_param_set use CLMFatesParamInterfaceMod , only : FatesReadParameters implicit none @@ -1089,6 +1090,22 @@ subroutine set_fates_global_elements(use_fates) nleafage = size(EDPftvarcon_inst%leaf_long,dim=2) end if + + ! Check that if age-dependent mortality is on, cohort age tracking is also on + !----------------------------------------------------------------------------------- + if ( ( ANY(EDPftvarcon_inst%mort_ip_age_senescence < fates_check_param_set )) .and. & + (hlm_use_cohort_age_tracking .eq. ifalse ) ) then + + write(fates_log(),*) 'Age dependent mortality cannot be on if' + write(fates_log(),*) 'cohort age tracking is off.' + write(fates_log(),*) 'Set hlm_use_cohort_age_tracking = .true.' + write(fates_log(),*) 'in FATES namelist options' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if + + + ! These values are used to define the restart file allocations and general structure ! of memory for the cohort arrays From fc733f9da5c381a5ce9eab71113a0b4c99cbb9ca Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Wed, 18 Mar 2020 10:12:58 -0700 Subject: [PATCH 21/27] [ Fix bug in error message for cohort age tracking flag ] [Make sure the test for cohort age tracking being properly set is only called AFTER hlm_use_cohort_age_tracking has been initialised. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary: No testing. --- main/FatesInterfaceMod.F90 | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 2a6aa169ed..b3088d7a1c 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1046,7 +1046,6 @@ subroutine set_fates_global_elements(use_fates) use EDParamsMod, only : ED_val_history_sizeclass_bin_edges, ED_val_history_ageclass_bin_edges use EDParamsMod, only : ED_val_history_height_bin_edges use EDParamsMod, only : ED_val_history_coageclass_bin_edges - use FatesConstantsMod, only : fates_check_param_set use CLMFatesParamInterfaceMod , only : FatesReadParameters implicit none @@ -1059,6 +1058,8 @@ subroutine set_fates_global_elements(use_fates) ! first read the non-PFT parameters call FatesReadParameters() + write(fates_log(),*)'hlm_use_cohort_age_tracking = ', hlm_use_cohort_age_tracking + ! Identify the number of PFTs by evaluating a pft array ! Using wood density as that is not expected to be deprecated any time soon @@ -1090,22 +1091,6 @@ subroutine set_fates_global_elements(use_fates) nleafage = size(EDPftvarcon_inst%leaf_long,dim=2) end if - - ! Check that if age-dependent mortality is on, cohort age tracking is also on - !----------------------------------------------------------------------------------- - if ( ( ANY(EDPftvarcon_inst%mort_ip_age_senescence < fates_check_param_set )) .and. & - (hlm_use_cohort_age_tracking .eq. ifalse ) ) then - - write(fates_log(),*) 'Age dependent mortality cannot be on if' - write(fates_log(),*) 'cohort age tracking is off.' - write(fates_log(),*) 'Set hlm_use_cohort_age_tracking = .true.' - write(fates_log(),*) 'in FATES namelist options' - write(fates_log(),*) 'Aborting' - call endrun(msg=errMsg(sourcefile, __LINE__)) - end if - - - ! These values are used to define the restart file allocations and general structure ! of memory for the cohort arrays @@ -1448,7 +1433,9 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) ! ! RGK-2016 ! --------------------------------------------------------------------------------- - + use FatesConstantsMod, only : fates_check_param_set + + ! Arguments integer, optional, intent(in) :: ival real(r8), optional, intent(in) :: rval @@ -1539,6 +1526,17 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) call endrun(msg=errMsg(sourcefile, __LINE__)) end if + + if ( ( ANY(EDPftvarcon_inst%mort_ip_age_senescence < fates_check_param_set )) .and. & + (hlm_use_cohort_age_tracking .eq. ifalse ) ) then + + write(fates_log(),*) 'Age dependent mortality cannot be on if' + write(fates_log(),*) 'cohort age tracking is off.' + write(fates_log(),*) 'Set hlm_use_cohort_age_tracking = .true.' + write(fates_log(),*) 'in FATES namelist options' + write(fates_log(),*) 'Aborting' + call endrun(msg=errMsg(sourcefile, __LINE__)) + end if if ( .not.((hlm_use_ed_st3.eq.1).or.(hlm_use_ed_st3.eq.0)) ) then From eed4545771bf1b490ba15ebb123d2d7160d12451 Mon Sep 17 00:00:00 2001 From: Jessica Needham Date: Wed, 18 Mar 2020 15:47:16 -0700 Subject: [PATCH 22/27] [ifalse to 0 in a test of hlm_use_cohort_age_flag ] [In the test to ensure cohort age tracking is on if age mortality is on, the test needs to be that cohort_age_tracking is 0, not ifalse. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:no testing --- main/FatesInterfaceMod.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index b3088d7a1c..9291bdbd03 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1058,8 +1058,6 @@ subroutine set_fates_global_elements(use_fates) ! first read the non-PFT parameters call FatesReadParameters() - write(fates_log(),*)'hlm_use_cohort_age_tracking = ', hlm_use_cohort_age_tracking - ! Identify the number of PFTs by evaluating a pft array ! Using wood density as that is not expected to be deprecated any time soon @@ -1528,7 +1526,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if ( ( ANY(EDPftvarcon_inst%mort_ip_age_senescence < fates_check_param_set )) .and. & - (hlm_use_cohort_age_tracking .eq. ifalse ) ) then + (hlm_use_cohort_age_tracking .eq.0 ) ) then write(fates_log(),*) 'Age dependent mortality cannot be on if' write(fates_log(),*) 'cohort age tracking is off.' From 22eb568540c7b2f0c5c79dfdf8d8ce75bb6a9b7b Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Thu, 19 Mar 2020 16:10:01 -0700 Subject: [PATCH 23/27] [ Update parameter file] [Cohort age fusion can be set at 0.08 as default. It will not be used if age tracking is off. It will get adjusted if age tracking is on. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing. --- parameter_files/fates_params_default.cdl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parameter_files/fates_params_default.cdl b/parameter_files/fates_params_default.cdl index ffa9e6633c..067ca4155f 100644 --- a/parameter_files/fates_params_default.cdl +++ b/parameter_files/fates_params_default.cdl @@ -1234,7 +1234,7 @@ data: fates_canopy_closure_thresh = 0.8 ; - fates_cohort_age_fusion_tol = _ ; + fates_cohort_age_fusion_tol = 0.08 ; fates_cohort_size_fusion_tol = 0.08 ; From 337ec1307fa73407222daaffc2f1fffd7f587163 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Fri, 20 Mar 2020 12:52:51 -0700 Subject: [PATCH 24/27] [ Remove checks on cohort age bin edges ] [If cohort age tracking is off these bin edges can be unset so we don't need a check on this. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:no testing --- main/FatesInterfaceMod.F90 | 4 ---- 1 file changed, 4 deletions(-) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 9291bdbd03..8bc84cd8fe 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1124,10 +1124,6 @@ subroutine set_fates_global_elements(use_fates) write(fates_log(), *) 'height class bins specified in parameter file must start at zero' call endrun(msg=errMsg(sourcefile, __LINE__)) endif - if ( ED_val_history_coageclass_bin_edges(1) .ne. 0._r8 ) then - write(fates_log(), *) 'cohort age class bines specified in parameter file must start at zero' - call endrun(msg=errMsg(sourcefile, __LINE__)) - endif do i = 2,nlevsclass if ( (ED_val_history_sizeclass_bin_edges(i) - ED_val_history_sizeclass_bin_edges(i-1)) .le. 0._r8) then write(fates_log(), *) 'size class bins specified in parameter file must be monotonically increasing' From f4afdf23540728305dbecd63e03571e889dc2f0b Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 24 Mar 2020 12:45:36 -0700 Subject: [PATCH 25/27] Fixed order of operations on mortality condition for age mortality, conditioned some history variables when age is off --- biogeochem/EDMortalityFunctionsMod.F90 | 9 ++++++-- main/FatesHistoryInterfaceMod.F90 | 29 ++++++++++++++++---------- 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index c27113f364..08c3120ae6 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -98,11 +98,16 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort,smort,asmor end if ! if param values have been set then calculate asmort + + + + mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) + mort_ip_age_senescence = EDPftvarcon_inst%mort_ip_age_senescence(cohort_in%pft) + if ( mort_ip_age_senescence < fates_check_param_set ) then ! Age Dependent Senescence ! rate and inflection point define the change in mortality with age - mort_r_age_senescence = EDPftvarcon_inst%mort_r_age_senescence(cohort_in%pft) - mort_ip_age_senescence = EDPftvarcon_inst%mort_ip_age_senescence(cohort_in%pft) + asmort = 1.0_r8 / (1.0_r8 + exp(-1.0_r8 * mort_r_age_senescence * & (cohort_in%coage - mort_ip_age_senescence ) ) ) else diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 31f7f0b15c..4a1692dc5d 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -29,6 +29,7 @@ module FatesHistoryInterfaceMod use FatesInterfaceMod , only : hlm_hio_ignore_val use FatesInterfaceMod , only : hlm_use_planthydro use FatesInterfaceMod , only : hlm_use_ed_st3 + use FatesInterfaceMod , only : hlm_use_cohort_age_tracking use FatesInterfaceMod , only : numpft use FatesInterfaceMod , only : hlm_freq_day use EDParamsMod , only : ED_val_comp_excln @@ -2274,10 +2275,15 @@ subroutine update_history_dyn(this,nc,nsites,sites) (ccohort%lmort_direct+ccohort%lmort_collateral+ccohort%lmort_infra) * ccohort%n hio_m8_si_scpf(io_si,scpf) = hio_m8_si_scpf(io_si,scpf) + ccohort%frmort*ccohort%n hio_m9_si_scpf(io_si,scpf) = hio_m9_si_scpf(io_si,scpf) + ccohort%smort*ccohort%n - hio_m10_si_scpf(io_si,scpf) = hio_m10_si_scpf(io_si,scpf) + ccohort%asmort*ccohort%n - - hio_m10_si_capf(io_si,capf) = hio_m10_si_capf(io_si,capf) + ccohort%asmort*ccohort%n - + + if (hlm_use_cohort_age_tracking .eq.itrue) then + hio_m10_si_scpf(io_si,scpf) = hio_m10_si_scpf(io_si,scpf) + ccohort%asmort*ccohort%n + hio_m10_si_capf(io_si,capf) = hio_m10_si_capf(io_si,capf) + ccohort%asmort*ccohort%n + hio_m10_si_scls(io_si,scls) = hio_m10_si_scls(io_si,scls) + ccohort%asmort*ccohort%n + hio_m10_si_cacls(io_si,cacls) = hio_m10_si_cacls(io_si,cacls)+ & + ccohort%asmort*ccohort%n + end if + hio_m1_si_scls(io_si,scls) = hio_m1_si_scls(io_si,scls) + ccohort%bmort*ccohort%n hio_m2_si_scls(io_si,scls) = hio_m2_si_scls(io_si,scls) + ccohort%hmort*ccohort%n hio_m3_si_scls(io_si,scls) = hio_m3_si_scls(io_si,scls) + ccohort%cmort*ccohort%n @@ -2287,10 +2293,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ccohort%frmort*ccohort%n hio_m9_si_scls(io_si,scls) = hio_m9_si_scls(io_si,scls) + ccohort%smort*ccohort%n - hio_m10_si_scls(io_si,scls) = hio_m10_si_scls(io_si,scls) + ccohort%asmort*ccohort%n - - hio_m10_si_cacls(io_si,cacls) = hio_m10_si_cacls(io_si,cacls)+ & - ccohort%asmort*ccohort%n + !C13 discrimination if(gpp_cached + ccohort%gpp_acc_hold > 0.0_r8)then @@ -2302,9 +2305,13 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! number density [/ha] hio_nplant_si_scpf(io_si,scpf) = hio_nplant_si_scpf(io_si,scpf) + ccohort%n + ! number density along the cohort age dimension - hio_nplant_si_capf(io_si,capf) = hio_nplant_si_capf(io_si,capf) + ccohort%n - + if (hlm_use_cohort_age_tracking .eq.itrue) then + hio_nplant_si_capf(io_si,capf) = hio_nplant_si_capf(io_si,capf) + ccohort%n + hio_nplant_si_cacls(io_si,cacls) = hio_nplant_si_cacls(io_si,cacls)+ccohort%n + end if + ! number density by size and biomass hio_agb_si_scls(io_si,scls) = hio_agb_si_scls(io_si,scls) + & total_c * ccohort%n * EDPftvarcon_inst%allom_agb_frac(ccohort%pft) * AREA_INV @@ -2322,7 +2329,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_nplant_si_scag(io_si,iscag) = hio_nplant_si_scag(io_si,iscag) + ccohort%n hio_nplant_si_scls(io_si,scls) = hio_nplant_si_scls(io_si,scls) + ccohort%n - hio_nplant_si_cacls(io_si,cacls) = hio_nplant_si_cacls(io_si,cacls)+ccohort%n + ! update size, age, and PFT - indexed quantities From ee90e0f83d4c8cca81acb865aff5b784fcf1a8e0 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Tue, 24 Mar 2020 15:32:52 -0700 Subject: [PATCH 26/27] [ Fix bug in mortality change in number density - logging should only be understory ] [Remove where I accidently added a dndt_logging. ] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: [Names] Spotted by rgknox Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing. --- biogeochem/EDMortalityFunctionsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDMortalityFunctionsMod.F90 b/biogeochem/EDMortalityFunctionsMod.F90 index 08c3120ae6..2a99cfb485 100644 --- a/biogeochem/EDMortalityFunctionsMod.F90 +++ b/biogeochem/EDMortalityFunctionsMod.F90 @@ -266,7 +266,7 @@ subroutine Mortality_Derivative( currentSite, currentCohort, bc_in) ! Mortality from logging in the canopy is ONLY disturbance generating, don't ! update number densities via non-disturbance inducing death currentCohort%dndt= -(1.0_r8-fates_mortality_disturbance_fraction) & - * (cmort+hmort+bmort+frmort+smort+asmort+dndt_logging) * & + * (cmort+hmort+bmort+frmort+smort+asmort) * & currentCohort%n endif From 10809b93fc63716bfe048d46e665f76cbf4531a5 Mon Sep 17 00:00:00 2001 From: JessicaNeedham Date: Thu, 26 Mar 2020 13:17:25 -0700 Subject: [PATCH 27/27] [Change maxCohortsPerPatch back to 100] [This was meant to go back to 100 now that we have a local variable for when cohort age tracking is on. Also fix alignment in EDCohortDynamics Fuse_cohorts.] Fixes: [NGT-ED Github issue #] User interface changes?: [No] Code review: Found by rgknox Test suite: [suite name, machine, compilers] Test baseline: Test namelist changes: Test answer changes: [bit for bit, roundoff, climate changing] Test summary:No testing --- biogeochem/EDCohortDynamicsMod.F90 | 668 ++++++++++++++--------------- main/EDTypesMod.F90 | 2 +- 2 files changed, 335 insertions(+), 335 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 699adfc710..09af1236a8 100644 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -1055,353 +1055,353 @@ subroutine fuse_cohorts(currentSite, currentPatch, bc_in) if (coage_diff <= dynamic_age_fusion_tolerance ) then - ! Don't fuse a cohort with itself! - if (.not.associated(currentCohort,nextc) ) then - - if (currentCohort%pft == nextc%pft) then - - ! check cohorts in same c. layer. before fusing - - if (currentCohort%canopy_layer == nextc%canopy_layer) then - - ! Note: because newly recruited cohorts that have not experienced - ! a day yet will have un-known flux quantities or change rates - ! we don't want them fusing with non-new cohorts. We allow them - ! to fuse with other new cohorts to keep the total number of cohorts - ! down. - - if( currentCohort%isnew.eqv.nextc%isnew ) then - - newn = currentCohort%n + nextc%n - - fusion_took_place = 1 - - if ( fuse_debug .and. currentCohort%isnew ) then - write(fates_log(),*) 'Fusing Two Cohorts' - write(fates_log(),*) 'newn: ',newn - write(fates_log(),*) 'Cohort I, Cohort II' - write(fates_log(),*) 'n:',currentCohort%n,nextc%n - write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew - write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory - write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite - write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage - write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh - write(fates_log(),*) 'pft:',currentCohort%pft,nextc%pft - write(fates_log(),*) 'canopy_trim:',currentCohort%canopy_trim,nextc%canopy_trim - write(fates_log(),*) 'canopy_layer_yesterday:', & - currentCohort%canopy_layer_yesterday,nextc%canopy_layer_yesterday - do i=1, nlevleaf - write(fates_log(),*) 'leaf level: ',i,'year_net_uptake', & - currentCohort%year_net_uptake(i),nextc%year_net_uptake(i) - end do - end if - - ! new cohort age is weighted mean of two cohorts - currentCohort%coage = & - (currentCohort%coage * (currentCohort%n/(currentCohort%n + nextc%n))) + & - (nextc%coage * (nextc%n/(currentCohort%n + nextc%n))) - - ! update the cohort age again - if (hlm_use_cohort_age_tracking .eq.itrue) then - call coagetype_class_index(currentCohort%coage, currentCohort%pft, & - currentCohort%coage_class, currentCohort%coage_by_pft_class) - end if - - ! Fuse all mass pools - call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, & - currentCohort%n/newn ) - - ! Leaf biophysical rates (use leaf mass weighting) - ! ----------------------------------------------------------------- - call UpdateCohortBioPhysRates(currentCohort) - - currentCohort%laimemory = (currentCohort%n*currentCohort%laimemory & - + nextc%n*nextc%laimemory)/newn - - currentCohort%sapwmemory = (currentCohort%n*currentCohort%sapwmemory & - + nextc%n*nextc%sapwmemory)/newn - - currentCohort%structmemory = (currentCohort%n*currentCohort%structmemory & - + nextc%n*nextc%structmemory)/newn - - currentCohort%canopy_trim = (currentCohort%n*currentCohort%canopy_trim & - + nextc%n*nextc%canopy_trim)/newn - - ! c13disc_acc calculation; weighted mean by GPP - if ((currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) .eq. 0.0_r8) then - currentCohort%c13disc_acc = 0.0_r8 - else - currentCohort%c13disc_acc = (currentCohort%n * currentCohort%gpp_acc * currentCohort%c13disc_acc + & - nextc%n * nextc%gpp_acc * nextc%c13disc_acc)/ & - (currentCohort%n * currentCohort%gpp_acc + nextc%n * nextc%gpp_acc) - endif - - select case(cohort_fusion_conservation_method) - ! - ! ----------------------------------------------------------------- - ! Because cohort fusion is an unavoidable but non-physical process, - ! and because of the various nonlinear allometric relationships, - ! it isn't possible to simultaneously conserve all of the allometric - ! relationships during cohort fusion. We will always conserve carbon, - ! but there are choices to made about what else to conserve or not. - ! In particular, there is a choice to be made of conservation amongst - ! the number density, stem diameter, and crown area. Below, - ! some different conservation relationships can be chosen during fusion. - ! ----------------------------------------------------------------- - ! - case(conserve_crownarea_and_number_not_dbh) - ! - ! ----------------------------------------------------------------- - ! conserve total crown area during the fusion step, and then calculate - ! dbh of the fused cohort as that which conserves both crown area and - ! the dbh to crown area allometry. dbh will be updated in the next - ! growth step in the (likely) event that dbh to structural iomass - ! allometry is exceeded. if using a capped crown area allometry and - ! above the cap, then calculate as the weighted average of fusing - ! cohorts' dbh + ! Don't fuse a cohort with itself! + if (.not.associated(currentCohort,nextc) ) then + + if (currentCohort%pft == nextc%pft) then + + ! check cohorts in same c. layer. before fusing + + if (currentCohort%canopy_layer == nextc%canopy_layer) then + + ! Note: because newly recruited cohorts that have not experienced + ! a day yet will have un-known flux quantities or change rates + ! we don't want them fusing with non-new cohorts. We allow them + ! to fuse with other new cohorts to keep the total number of cohorts + ! down. + + if( currentCohort%isnew.eqv.nextc%isnew ) then + + newn = currentCohort%n + nextc%n + + fusion_took_place = 1 + + if ( fuse_debug .and. currentCohort%isnew ) then + write(fates_log(),*) 'Fusing Two Cohorts' + write(fates_log(),*) 'newn: ',newn + write(fates_log(),*) 'Cohort I, Cohort II' + write(fates_log(),*) 'n:',currentCohort%n,nextc%n + write(fates_log(),*) 'isnew:',currentCohort%isnew,nextc%isnew + write(fates_log(),*) 'laimemory:',currentCohort%laimemory,nextc%laimemory + write(fates_log(),*) 'hite:',currentCohort%hite,nextc%hite + write(fates_log(),*) 'coage:',currentCohort%coage,nextc%coage + write(fates_log(),*) 'dbh:',currentCohort%dbh,nextc%dbh + write(fates_log(),*) 'pft:',currentCohort%pft,nextc%pft + write(fates_log(),*) 'canopy_trim:',currentCohort%canopy_trim,nextc%canopy_trim + write(fates_log(),*) 'canopy_layer_yesterday:', & + currentCohort%canopy_layer_yesterday,nextc%canopy_layer_yesterday + do i=1, nlevleaf + write(fates_log(),*) 'leaf level: ',i,'year_net_uptake', & + currentCohort%year_net_uptake(i),nextc%year_net_uptake(i) + end do + end if + + ! new cohort age is weighted mean of two cohorts + currentCohort%coage = & + (currentCohort%coage * (currentCohort%n/(currentCohort%n + nextc%n))) + & + (nextc%coage * (nextc%n/(currentCohort%n + nextc%n))) + + ! update the cohort age again + if (hlm_use_cohort_age_tracking .eq.itrue) then + call coagetype_class_index(currentCohort%coage, currentCohort%pft, & + currentCohort%coage_class, currentCohort%coage_by_pft_class) + end if + + ! Fuse all mass pools + call currentCohort%prt%WeightedFusePRTVartypes(nextc%prt, & + currentCohort%n/newn ) + + ! Leaf biophysical rates (use leaf mass weighting) ! ----------------------------------------------------------------- - ! - - call carea_allom(currentCohort%dbh,currentCohort%n, & - currentSite%spread,currentCohort%pft,& - currentCohort%c_area,inverse=.false.) - - call carea_allom(nextc%dbh,nextc%n, & - currentSite%spread,nextc%pft,& - nextc%c_area,inverse=.false.) - - currentCohort%c_area = currentCohort%c_area + nextc%c_area - - ! - call carea_allom(dbh,newn,currentSite%spread,currentCohort%pft,& - currentCohort%c_area,inverse=.true.) - ! - if (abs(dbh-fates_unset_r8) nextc%shorter - tallerCohort => nextc%taller - - if (.not. associated(tallerCohort)) then - currentPatch%tallest => shorterCohort - if(associated(shorterCohort)) shorterCohort%taller => null() - else - tallerCohort%shorter => shorterCohort - endif - - if (.not. associated(shorterCohort)) then - currentPatch%shortest => tallerCohort - if(associated(tallerCohort)) tallerCohort%shorter => null() - else - shorterCohort%taller => tallerCohort - endif - - ! At this point, nothing should be pointing to current Cohort - ! update hydraulics quantities that are functions of hite & biomasses - ! deallocate the hydro structure of nextc - if (hlm_use_planthydro.eq.itrue) then - call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & - currentCohort%pft,currentCohort%c_area) - leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) - currentCohort%treelai = tree_lai(leaf_c, & - currentCohort%pft, currentCohort%c_area, currentCohort%n, & - currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & - currentCohort%vcmax25top ) - call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in) - endif - - call DeallocateCohort(nextc) - deallocate(nextc) - nullify(nextc) - - endif ! if( currentCohort%isnew.eqv.nextc%isnew ) then - endif !canopy layer - endif !pft - endif !index no. - endif ! cohort age diff - endif !diff - + ! + ! it is possible that fusion has caused cohorts separated by at least two size bin deltas to join. + ! so slightly complicated to keep track of because the resulting cohort could be in one of the old bins or in between + ! structure as a loop to handle the general case + ! + ! first the positive growth case + do sc_i = smallersc + 1, currentCohort%size_class + currentSite%growthflux_fusion(sc_i, currentCohort%pft) = & + currentSite%growthflux_fusion(sc_i, currentCohort%pft) + smaller_n + end do + ! + ! next the negative growth case + do sc_i = currentCohort%size_class + 1, largersc + currentSite%growthflux_fusion(sc_i, currentCohort%pft) = & + currentSite%growthflux_fusion(sc_i, currentCohort%pft) - larger_n + end do + ! now that we've tracked the change flux. reset the memory of the prior timestep + currentCohort%size_class_lasttimestep = currentCohort%size_class + endif + + + ! Flux and biophysics variables have not been calculated for recruits we just default to + ! their initization values, which should be the same for each + + if ( .not.currentCohort%isnew) then + currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + & + nextc%n*nextc%seed_prod)/newn + currentCohort%gpp_acc = (currentCohort%n*currentCohort%gpp_acc + & + nextc%n*nextc%gpp_acc)/newn + currentCohort%npp_acc = (currentCohort%n*currentCohort%npp_acc + & + nextc%n*nextc%npp_acc)/newn + currentCohort%resp_acc = (currentCohort%n*currentCohort%resp_acc + & + nextc%n*nextc%resp_acc)/newn + currentCohort%resp_acc_hold = & + (currentCohort%n*currentCohort%resp_acc_hold + & + nextc%n*nextc%resp_acc_hold)/newn + currentCohort%npp_acc_hold = & + (currentCohort%n*currentCohort%npp_acc_hold + & + nextc%n*nextc%npp_acc_hold)/newn + currentCohort%gpp_acc_hold = & + (currentCohort%n*currentCohort%gpp_acc_hold + & + nextc%n*nextc%gpp_acc_hold)/newn + + currentCohort%dmort = (currentCohort%n*currentCohort%dmort + & + nextc%n*nextc%dmort)/newn + + currentCohort%fire_mort = (currentCohort%n*currentCohort%fire_mort + & + nextc%n*nextc%fire_mort)/newn + + ! mortality diagnostics + currentCohort%cmort = (currentCohort%n*currentCohort%cmort + nextc%n*nextc%cmort)/newn + currentCohort%hmort = (currentCohort%n*currentCohort%hmort + nextc%n*nextc%hmort)/newn + currentCohort%bmort = (currentCohort%n*currentCohort%bmort + nextc%n*nextc%bmort)/newn + currentCohort%smort = (currentCohort%n*currentCohort%smort + nextc%n*nextc%smort)/newn + currentCohort%asmort = (currentCohort%n*currentCohort%asmort + nextc%n*nextc%asmort)/newn + currentCohort%frmort = (currentCohort%n*currentCohort%frmort + nextc%n*nextc%frmort)/newn + + ! logging mortality, Yi Xu + currentCohort%lmort_direct = (currentCohort%n*currentCohort%lmort_direct + & + nextc%n*nextc%lmort_direct)/newn + currentCohort%lmort_collateral = (currentCohort%n*currentCohort%lmort_collateral + & + nextc%n*nextc%lmort_collateral)/newn + currentCohort%lmort_infra = (currentCohort%n*currentCohort%lmort_infra + & + nextc%n*nextc%lmort_infra)/newn + currentCohort%l_degrad = (currentCohort%n*currentCohort%l_degrad + & + nextc%n*nextc%l_degrad)/newn + + ! biomass and dbh tendencies + currentCohort%ddbhdt = (currentCohort%n*currentCohort%ddbhdt + & + nextc%n*nextc%ddbhdt)/newn + + do i=1, nlevleaf + if (currentCohort%year_net_uptake(i) == 999._r8 .or. nextc%year_net_uptake(i) == 999._r8) then + currentCohort%year_net_uptake(i) = & + min(nextc%year_net_uptake(i),currentCohort%year_net_uptake(i)) + else + currentCohort%year_net_uptake(i) = (currentCohort%n*currentCohort%year_net_uptake(i) + & + nextc%n*nextc%year_net_uptake(i))/newn + endif + enddo + + end if !(currentCohort%isnew) + + currentCohort%n = newn + + ! Set pointers and remove the current cohort from the list + + shorterCohort => nextc%shorter + tallerCohort => nextc%taller + + if (.not. associated(tallerCohort)) then + currentPatch%tallest => shorterCohort + if(associated(shorterCohort)) shorterCohort%taller => null() + else + tallerCohort%shorter => shorterCohort + endif + + if (.not. associated(shorterCohort)) then + currentPatch%shortest => tallerCohort + if(associated(tallerCohort)) tallerCohort%shorter => null() + else + shorterCohort%taller => tallerCohort + endif + + ! At this point, nothing should be pointing to current Cohort + ! update hydraulics quantities that are functions of hite & biomasses + ! deallocate the hydro structure of nextc + if (hlm_use_planthydro.eq.itrue) then + call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread, & + currentCohort%pft,currentCohort%c_area) + leaf_c = currentCohort%prt%GetState(leaf_organ, carbon12_element) + currentCohort%treelai = tree_lai(leaf_c, & + currentCohort%pft, currentCohort%c_area, currentCohort%n, & + currentCohort%canopy_layer, currentPatch%canopy_layer_tlai, & + currentCohort%vcmax25top ) + call updateSizeDepTreeHydProps(currentSite,currentCohort, bc_in) + endif + + call DeallocateCohort(nextc) + deallocate(nextc) + nullify(nextc) + + endif ! if( currentCohort%isnew.eqv.nextc%isnew ) then + endif !canopy layer + endif !pft + endif !index no. + endif ! cohort age diff + endif !diff + nextc => nextnextc - + enddo !end checking nextc cohort loop ! Ususally we always point to the next cohort. But remember ... diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index d95f2f6baa..913c1a8ef5 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -23,7 +23,7 @@ module EDTypesMod integer, parameter, public :: maxPatchesPerSite = 14 ! maximum number of patches to live on a site integer, parameter, public :: maxPatchesPerSite_by_disttype(n_anthro_disturbance_categories) = & (/ 10, 4 /) !!! MUST SUM TO maxPatchesPerSite !!! - integer, parameter, public :: maxCohortsPerPatch = 300 ! maximum number of cohorts per patch + integer, parameter, public :: maxCohortsPerPatch = 100 ! maximum number of cohorts per patch integer, parameter, public :: nclmax = 2 ! Maximum number of canopy layers integer, parameter, public :: ican_upper = 1 ! Nominal index for the upper canopy