From 206df5272e7392d8cf12a18c4c7d12c49afa4980 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Tue, 26 Sep 2017 17:42:50 -0700 Subject: [PATCH 01/10] change to canopy spread logic; move to site level and also make a pft trait --- biogeochem/EDCanopyStructureMod.F90 | 50 ++++++++++++++--------------- biogeochem/EDGrowthFunctionsMod.F90 | 13 ++++---- biogeochem/EDPatchDynamicsMod.F90 | 13 ++------ main/EDInitMod.F90 | 12 ++++--- main/EDParamsMod.F90 | 50 +++++++---------------------- main/EDPftvarcon.F90 | 20 ++++++++++++ main/EDTypesMod.F90 | 4 ++- main/FatesHistoryInterfaceMod.F90 | 12 +++---- main/FatesInventoryInitMod.F90 | 5 +-- main/FatesRestartInterfaceMod.F90 | 35 ++++++-------------- 10 files changed, 93 insertions(+), 121 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index ddab8e74af..3a58ce73ae 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -756,7 +756,8 @@ subroutine canopy_spread( currentSite ) ! Calculates the spatial spread of tree canopies based on canopy closure. ! ! !USES: - use EDParamsMod , only : ED_val_maxspread, ED_val_minspread + use EDTypesMod , only : AREA + use EDParamsMod, only : ED_val_canopy_closure_thresh ! ! !ARGUMENTS type (ed_site_type), intent(inout), target :: currentSite @@ -769,46 +770,43 @@ subroutine canopy_spread( currentSite ) integer :: z !---------------------------------------------------------------------- - inc = 0.005_r8 + inc = 0.05_r8 currentPatch => currentSite%oldest_patch + arealayer = 0.0_r8 do while (associated(currentPatch)) - !calculate canopy area in each canopy storey... - arealayer = 0.0_r8 + !calculate canopy area in each patch... currentCohort => currentPatch%tallest do while (associated(currentCohort)) currentCohort%c_area = c_area(currentCohort) - if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then - arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area + if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then + arealayer = arealayer + currentCohort%c_area endif currentCohort => currentCohort%shorter enddo - !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner - do z = 1,nclmax - - if(arealayer(z)/currentPatch%area > 0.9_r8)then - currentPatch%spread(z) = currentPatch%spread(z) - inc - else - currentPatch%spread(z) = currentPatch%spread(z) + inc - endif - if(currentPatch%spread(z) >= ED_val_maxspread)then - currentPatch%spread(z) = ED_val_maxspread - endif - if(currentPatch%spread(z) <= ED_val_minspread)then - currentPatch%spread(z) = ED_val_minspread - endif - enddo !z - !write(fates_log(),*) 'spread',currentPatch%spread(1:2) - !currentPatch%spread(:) = ED_val_maxspread - !FIX(RF,033114) spread is off - !write(fates_log(),*) 'canopy_spread',currentPatch%area,currentPatch%spread(1:2) - currentPatch => currentPatch%younger + currentPatch => currentPatch%younger enddo !currentPatch + !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner + if( arealayer/AREA .gt. ED_val_canopy_closure_thresh ) then + currentSite%spread = currentSite%spread - inc + else + currentSite%spread = currentSite%spread + inc + endif + if(currentSite%spread >= 1._r8)then + currentSite%spread = 1._r8 + endif + if(currentSite%spread <= 0._r8)then + currentSite%spread = 0._r8 + endif + + ! put within bounds to make sure it stays between 0 and 1 + currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8) + end subroutine canopy_spread diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 index 47bbb9591c..7abb19a19a 100755 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -225,13 +225,13 @@ real(r8) function c_area( cohort_in ) ! ============================================================================ use EDParamsMod , only : ED_val_grass_spread + use EDParamsMod , only : ED_val_maxspread, ED_val_minspread use EDTypesMod , only : nclmax type(ed_cohort_type), intent(in) :: cohort_in real(r8) :: dbh ! Tree diameter at breat height. cm. real(r8) :: crown_area_to_dbh_exponent - integer :: can_layer_index ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth, ! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero) @@ -259,12 +259,13 @@ real(r8) function c_area( cohort_in ) ! So we allow layer index exceedence here and force it down to max. ! (rgk/cdk 05/2017) ! ---------------------------------------------------------------------------------- - - can_layer_index = min(cohort_in%canopy_layer,nclmax) - - if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then + + if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then + ! apply site-level spread elasticity to the cohort crown allometry term + spreadterm = cohort_in%siteptr%spread * ED_val_maxspread + (1._r8 - cohort_in%siteptr%spread * ED_val_minspread) + ! c_area = 3.142_r8 * cohort_in%n * & - (cohort_in%patchptr%spread(can_layer_index)*dbh)**crown_area_to_dbh_exponent + (spreadterm * dbh)**crown_area_to_dbh_exponent else c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent end if diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 533073b0b8..ae6b13c716 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -187,7 +187,7 @@ subroutine spawn_patches( currentSite, bc_in) ! ! !USES: - use EDParamsMod , only : ED_val_maxspread, ED_val_understorey_death + use EDParamsMod , only : ED_val_understorey_death use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts ! @@ -211,7 +211,6 @@ subroutine spawn_patches( currentSite, bc_in) real(r8) :: leaf_litter_local(maxpft) ! initial value of leaf litter. KgC/m2 real(r8) :: cwd_ag_local(ncwd) ! initial value of above ground coarse woody debris. KgC/m2 real(r8) :: cwd_bg_local(ncwd) ! initial value of below ground coarse woody debris. KgC/m2 - real(r8) :: spread_local(nclmax) ! initial value of canopy spread parameter.no units !--------------------------------------------------------------------- storesmallcohort => null() ! storage of the smallest cohort for insertion routine @@ -239,7 +238,6 @@ subroutine spawn_patches( currentSite, bc_in) cwd_bg_local = 0.0_r8 leaf_litter_local = 0.0_r8 root_litter_local = 0.0_r8 - spread_local(1:nclmax) = ED_val_maxspread age = 0.0_r8 allocate(new_patch) @@ -250,7 +248,7 @@ subroutine spawn_patches( currentSite, bc_in) ! call zero_patch(new_patch) call create_patch(currentSite, new_patch, age, site_areadis, & - spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & + cwd_ag_local, cwd_bg_local, leaf_litter_local, & root_litter_local) new_patch%tallest => null() @@ -537,8 +535,6 @@ subroutine average_patch_properties( currentPatch, newPatch, patch_site_areadis enddo - newPatch%spread = newPatch%spread + currentPatch%spread * patch_site_areadis/newPatch%area - end subroutine average_patch_properties ! ============================================================================ @@ -868,7 +864,7 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat end subroutine mortality_litter_fluxes ! ============================================================================ - subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_local,cwd_bg_local, & + subroutine create_patch(currentSite, new_patch, age, areap,cwd_ag_local,cwd_bg_local, & leaf_litter_local,root_litter_local) ! ! !DESCRIPTION: @@ -885,7 +881,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ real(r8), intent(in) :: cwd_bg_local(:) ! initial value of below ground coarse woody debris. KgC/m2 real(r8), intent(in) :: root_litter_local(:)! initial value of root litter. KgC/m2 real(r8), intent(in) :: leaf_litter_local(:)! initial value of leaf litter. KgC/m2 - real(r8), intent(in) :: spread_local(:) ! initial value of canopy spread parameter.no units ! ! !LOCAL VARIABLES: !--------------------------------------------------------------------- @@ -915,7 +910,6 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_ new_patch%age = age new_patch%age_class = 1 new_patch%area = areap - new_patch%spread = spread_local new_patch%cwd_ag = cwd_ag_local new_patch%cwd_bg = cwd_bg_local new_patch%leaf_litter = leaf_litter_local @@ -1012,7 +1006,6 @@ subroutine zero_patch(cp_p) currentPatch%nrad(:,:) = 999 ! number of exposed leaf layers for each canopy layer and pft currentPatch%ncan(:,:) = 999 ! number of total leaf layers for each canopy layer and pft currentPatch%lai = nan ! leaf area index of patch - currentPatch%spread(:) = nan ! dynamic ratio of dbh to canopy area. currentPatch%pft_agb_profile(:,:) = nan ! DISTURBANCE diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 73bda67f5d..b2bde7517e 100755 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -108,6 +108,9 @@ subroutine zero_site( site_in ) site_in%leaf_litter_diagnostic_input_carbonflux(:) = 0._r8 site_in%root_litter_diagnostic_input_carbonflux(:) = 0._r8 + ! canopy spread + site_in%spread = 0._r8 + end subroutine zero_site ! ============================================================================ @@ -116,6 +119,7 @@ subroutine set_site_properties( nsites, sites) ! !DESCRIPTION: ! ! !USES: + use EDParamsMod , only : ED_val_maxspread ! ! !ARGUMENTS @@ -134,6 +138,7 @@ subroutine set_site_properties( nsites, sites) real(r8) :: watermem integer :: dleafoff integer :: dleafon + real(r8) :: spread_local !---------------------------------------------------------------------- if ( hlm_is_restart == ifalse ) then @@ -184,7 +189,7 @@ subroutine set_site_properties( nsites, sites) sites(s)%frac_burnt = 0.0_r8 sites(s)%old_stock = 0.0_r8 - + sites(s)%spread_local = 1.0_r8 end do return @@ -201,7 +206,6 @@ subroutine init_patches( nsites, sites, bc_in) ! - use EDParamsMod , only : ED_val_maxspread use FatesPlantHydraulicsMod, only : updateSizeDepRhizHydProps use FatesInventoryInitMod, only : initialize_sites_by_inventory @@ -215,7 +219,6 @@ subroutine init_patches( nsites, sites, bc_in) integer :: s real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) - real(r8) :: spread_local(nclmax) real(r8) :: leaf_litter_local(maxpft) real(r8) :: root_litter_local(maxpft) real(r8) :: age !notional age of this patch @@ -228,7 +231,6 @@ subroutine init_patches( nsites, sites, bc_in) cwd_bg_local(:) = 0.0_r8 !ED_val_init_litter leaf_litter_local(:) = 0.0_r8 root_litter_local(:) = 0.0_r8 - spread_local(:) = ED_val_maxspread age = 0.0_r8 ! --------------------------------------------------------------------------------------------- @@ -263,7 +265,7 @@ subroutine init_patches( nsites, sites, bc_in) ! make new patch... call create_patch(sites(s), newp, age, AREA, & - spread_local, cwd_ag_local, cwd_bg_local, leaf_litter_local, & + cwd_ag_local, cwd_bg_local, leaf_litter_local, & root_litter_local) call init_cohorts(newp, bc_in(s)) diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index 34e4cfa807..34a7f28b61 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -2,6 +2,7 @@ module EDParamsMod ! ! module that deals with reading the ED parameter file ! +canopyclosure_threshold use FatesParametersInterface, only : param_string_length use FatesGlobals , only : fates_log @@ -23,11 +24,8 @@ module EDParamsMod real(r8),protected :: ED_size_diagnostic_scale ! Flag to switch between a linear and exponential ! scale on the plant size axis in diagnostics (NOT USED YET) real(r8),protected :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance - real(r8),protected :: ED_val_grass_spread real(r8),protected :: ED_val_comp_excln real(r8),protected :: ED_val_stress_mort - real(r8),protected :: ED_val_maxspread ! maximum ratio of dbh to canopy area (cm/m2) - real(r8),protected :: ED_val_minspread ! minimum ratio of dbh to canopy area (cm/m2) real(r8),protected :: ED_val_init_litter real(r8),protected :: ED_val_nignitions real(r8),protected :: ED_val_understorey_death @@ -47,14 +45,12 @@ module EDParamsMod real(r8),protected :: ED_val_phen_coldtemp real(r8),protected :: ED_val_cohort_fusion_tol real(r8),protected :: ED_val_patch_fusion_tol + real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry character(len=param_string_length),parameter :: ED_name_size_diagnostic_scale = "fates_size_diagnostic_scale" character(len=param_string_length),parameter :: ED_name_mort_disturb_frac = "fates_mort_disturb_frac" - character(len=param_string_length),parameter :: ED_name_grass_spread = "fates_grass_spread" character(len=param_string_length),parameter :: ED_name_comp_excln = "fates_comp_excln" character(len=param_string_length),parameter :: ED_name_stress_mort = "fates_stress_mort" - character(len=param_string_length),parameter :: ED_name_maxspread = "fates_maxspread" - character(len=param_string_length),parameter :: ED_name_minspread = "fates_minspread" character(len=param_string_length),parameter :: ED_name_init_litter = "fates_init_litter" character(len=param_string_length),parameter :: ED_name_nignitions = "fates_nignitions" character(len=param_string_length),parameter :: ED_name_understorey_death = "fates_understorey_death" @@ -73,7 +69,8 @@ module EDParamsMod character(len=param_string_length),parameter :: ED_name_phen_ncolddayslim= "fates_phen_ncolddayslim" character(len=param_string_length),parameter :: ED_name_phen_coldtemp= "fates_phen_coldtemp" character(len=param_string_length),parameter :: ED_name_cohort_fusion_tol= "fates_cohort_fusion_tol" - character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" + character(len=param_string_length),parameter :: ED_name_patch_fusion_tol= "fates_patch_fusion_tol" + character(len=param_string_length),parameter :: ED_name_canopy_closure_thresh= "fates_canopy_closure_thresh" ! Hydraulics Control Parameters (ONLY RELEVANT WHEN USE_FATES_HYDR = TRUE) ! ---------------------------------------------------------------------------------------------- @@ -121,11 +118,8 @@ subroutine FatesParamsInit() ED_size_diagnostic_scale = nan fates_mortality_disturbance_fraction = nan - ED_val_grass_spread = nan ED_val_comp_excln = nan ED_val_stress_mort = nan - ED_val_maxspread = nan - ED_val_minspread = nan ED_val_init_litter = nan ED_val_nignitions = nan ED_val_understorey_death = nan @@ -145,6 +139,7 @@ subroutine FatesParamsInit() ED_val_phen_coldtemp = nan ED_val_cohort_fusion_tol = nan ED_val_patch_fusion_tol = nan + ED_val_canopy_closure_thresh = nan hydr_psi0 = nan hydr_psicap = nan @@ -179,27 +174,15 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_mort_disturb_frac, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_grass_spread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_grass_spread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_comp_excln, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) call fates_params%RegisterParameter(name=ED_name_stress_mort, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_maxspread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - - call fates_params%RegisterParameter(name=ED_name_minspread, dimension_shape=dimension_shape_1d, & - dimension_names=dim_names) - call fates_params%RegisterParameter(name=ED_name_init_litter, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -257,6 +240,9 @@ subroutine FatesRegisterParams(fates_params) call fates_params%RegisterParameter(name=ED_name_patch_fusion_tol, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) + call fates_params%RegisterParameter(name=ED_name_canopy_closure_thresh, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names) + call fates_params%RegisterParameter(name=hydr_name_psi0, dimension_shape=dimension_shape_1d, & dimension_names=dim_names) @@ -297,27 +283,15 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_mort_disturb_frac, & data=fates_mortality_disturbance_fraction) - call fates_params%RetreiveParameter(name=ED_name_grass_spread, & - data=ED_val_grass_spread) - call fates_params%RetreiveParameter(name=ED_name_comp_excln, & data=ED_val_comp_excln) - call fates_params%RetreiveParameter(name=ED_name_grass_spread, & - data=ED_val_grass_spread) - call fates_params%RetreiveParameter(name=ED_name_comp_excln, & data=ED_val_comp_excln) call fates_params%RetreiveParameter(name=ED_name_stress_mort, & data=ED_val_stress_mort) - call fates_params%RetreiveParameter(name=ED_name_maxspread, & - data=ED_val_maxspread) - - call fates_params%RetreiveParameter(name=ED_name_minspread, & - data=ED_val_minspread) - call fates_params%RetreiveParameter(name=ED_name_init_litter, & data=ED_val_init_litter) @@ -375,6 +349,9 @@ subroutine FatesReceiveParams(fates_params) call fates_params%RetreiveParameter(name=ED_name_patch_fusion_tol, & data=ED_val_patch_fusion_tol) + call fates_params%RetreiveParameter(name=ED_name_canopy_closure_thresh, & + data=ED_val_canopy_closure_thresh) + call fates_params%RetreiveParameter(name=hydr_name_psi0, & data=hydr_psi0) @@ -412,13 +389,9 @@ subroutine FatesReportParams(is_master) write(fates_log(),*) '----------- FATES Scalar Parameters -----------------' write(fates_log(),fmt0) 'ED_size_diagnostic_scale = ',ED_size_diagnostic_scale write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction - write(fates_log(),fmt0) 'ED_val_grass_spread = ',ED_val_grass_spread write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln - write(fates_log(),fmt0) 'ED_val_grass_spread = ',ED_val_grass_spread write(fates_log(),fmt0) 'ED_val_comp_excln = ', ED_val_comp_excln write(fates_log(),fmt0) 'ED_val_stress_mort = ',ED_val_stress_mort - write(fates_log(),fmt0) 'ED_val_maxspread = ',ED_val_maxspread - write(fates_log(),fmt0) 'ED_val_minspread = ',ED_val_minspread write(fates_log(),fmt0) 'ED_val_init_litter = ',ED_val_init_litter write(fates_log(),fmt0) 'ED_val_nignitions = ',ED_val_nignitions write(fates_log(),fmt0) 'ED_val_understorey_death = ',ED_val_understorey_death @@ -438,6 +411,7 @@ subroutine FatesReportParams(is_master) 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_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_psi0 = ',hydr_psi0 write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin diff --git a/main/EDPftvarcon.F90 b/main/EDPftvarcon.F90 index d4b9fdf18a..dfa3dcf0a0 100644 --- a/main/EDPftvarcon.F90 +++ b/main/EDPftvarcon.F90 @@ -120,6 +120,8 @@ module EDPftvarcon real(r8), allocatable :: allom_d2bl_slascaler(:) ! real(r8), allocatable :: allom_blca_expnt_diff(:) ! Any difference in the exponent between the leaf ! biomass and crown area scaling + real(r8), allocatable :: allom_d2ca_coefficient_max(:) ! upper (savanna) value for crown area to dbh coefficient + real(r8), allocatable :: allom_d2ca_coefficient_min(:) ! lower (closed-canopy forest) value for crown area to dbh coefficient real(r8), allocatable :: allom_agb1(:) ! Parameter 1 for agb allometry real(r8), allocatable :: allom_agb2(:) ! Parameter 2 for agb allometry real(r8), allocatable :: allom_agb3(:) ! Parameter 3 for agb allometry @@ -462,6 +464,14 @@ 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_allom_d2ca_coefficient_max' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + + name = 'fates_allom_d2ca_coefficient_min' + call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & + dimension_names=dim_names, lower_bounds=dim_lower_bound) + name = 'fates_allom_d2bl_slascaler' call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, & dimension_names=dim_names, lower_bounds=dim_lower_bound) @@ -830,6 +840,14 @@ subroutine Receive_PFT(this, fates_params) call fates_params%RetreiveParameterAllocate(name=name, & data=this%allom_blca_expnt_diff) + name = 'fates_allom_d2ca_coefficient_max' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%allom_d2ca_coefficient_max) + + name = 'fates_allom_d2ca_coefficient_min' + call fates_params%RetreiveParameterAllocate(name=name, & + data=this%allom_d2ca_coefficient_min) + name = 'fates_allom_d2bl_slascaler' call fates_params%RetreiveParameterAllocate(name=name, & data=this%allom_d2bl_slascaler) @@ -1394,6 +1412,8 @@ subroutine FatesReportPFTParams(is_master) write(fates_log(),fmt0) 'allom_sai_scaler = ',EDPftvarcon_inst%allom_sai_scaler write(fates_log(),fmt0) 'allom_d2bl_slascaler = ',EDPftvarcon_inst%allom_d2bl_slascaler write(fates_log(),fmt0) 'allom_blca_expnt_diff = ',EDPftvarcon_inst%allom_blca_expnt_diff + write(fates_log(),fmt0) 'allom_d2ca_coefficient_max = ',EDPftvarcon_inst%allom_d2ca_coefficient_max + write(fates_log(),fmt0) 'allom_d2ca_coefficient_min = ',EDPftvarcon_inst%allom_d2ca_coefficient_min write(fates_log(),fmt0) 'allom_agb1 = ',EDPftvarcon_inst%allom_agb1 write(fates_log(),fmt0) 'allom_agb2 = ',EDPftvarcon_inst%allom_agb2 write(fates_log(),fmt0) 'allom_agb3 = ',EDPftvarcon_inst%allom_agb3 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index e0b7f2b44c..a6631cea6c 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -298,7 +298,6 @@ module EDTypesMod integer :: ncl_p ! Number of occupied canopy layers ! LEAF ORGANIZATION - real(r8) :: spread(nclmax) ! dynamic ratio of dbh to canopy area: cm/m2 real(r8) :: pft_agb_profile(maxpft,n_dbh_bins) ! binned above ground biomass, for patch fusion: KgC/m2 real(r8) :: canopy_layer_lai(nclmax) ! lai that is shading this canopy layer: m2/m2 real(r8) :: total_canopy_area ! area that is covered by vegetation : m2 @@ -538,6 +537,9 @@ module EDTypesMod real(r8) :: leaf_litter_diagnostic_input_carbonflux(1:maxpft) ! diagnostic flux to AG litter [kg C / m2 / yr] real(r8) :: root_litter_diagnostic_input_carbonflux(1:maxpft) ! diagnostic flux to BG litter [kg C / m2 / yr] + ! Canopy Spread + real(r8) :: spread ! dynamic canopy allometric term [unitless] + end type ed_site_type contains diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 0f93936e2b..bc4f2d4ba0 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -41,7 +41,6 @@ module FatesHistoryInterfaceMod integer, private :: ih_trimming_pa integer, private :: ih_area_plant_pa integer, private :: ih_area_treespread_pa - integer, private :: ih_canopy_spread_pa integer, private :: ih_nesterov_fire_danger_pa integer, private :: ih_spitfire_ROS_pa integer, private :: ih_effect_wspeed_pa @@ -121,6 +120,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_promotion_carbonflux_si integer, private :: ih_canopy_mortality_carbonflux_si integer, private :: ih_understory_mortality_carbonflux_si + integer, private :: ih_canopy_spread_si ! Indices to (site x scpf) variables integer, private :: ih_nplant_si_scpf @@ -1305,6 +1305,8 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! The seed bank is a site level variable hio_seed_bank_si(io_si) = sum(sites(s)%seed_bank) * g_per_kg + hio_canopy_spread_si(io_si) = sites(s)%spread + ipa = 0 cpatch => sites(s)%oldest_patch do while(associated(cpatch)) @@ -1328,7 +1330,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_zstar_si_age(io_si,cpatch%age_class) = hio_zstar_si_age(io_si,cpatch%age_class) & + cpatch%zstar * cpatch%area * AREA_INV endif - + ccohort => cpatch%shortest do while(associated(ccohort)) @@ -1686,8 +1688,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) g_per_kg * patch_scaling_scalar * years_per_day * days_per_sec - hio_canopy_spread_pa(io_pa) = cpatch%spread(1) - do i_cwd = 1, ncwd hio_cwd_ag_si_cwdsc(io_si, i_cwd) = hio_cwd_ag_si_cwdsc(io_si, i_cwd) + & cpatch%CWD_AG(i_cwd)*cpatch%area * AREA_INV * g_per_kg @@ -2506,8 +2506,8 @@ subroutine define_history_vars(this, initialize_variables) call this%set_history_var(vname='CANOPY_SPREAD', units='0-1', & long='Scaling factor between tree basal area and canopy area', & use_default='active', & - avgflag='A', vtype=patch_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & - ivar=ivar, initialize=initialize_variables, index = ih_canopy_spread_pa) + avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_canopy_spread_si) call this%set_history_var(vname='PFTbiomass', units='gC/m2', & long='total PFT level biomass', use_default='active', & diff --git a/main/FatesInventoryInitMod.F90 b/main/FatesInventoryInitMod.F90 index 52a75e480c..1f063a2cb7 100644 --- a/main/FatesInventoryInitMod.F90 +++ b/main/FatesInventoryInitMod.F90 @@ -77,7 +77,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) use EDTypesMod, only : nclmax use EDTypesMod, only : maxpft use EDTypesMod, only : ncwd - use EDParamsMod, only : ED_val_maxspread use EDPatchDynamicsMod, only : create_patch use EDPatchDynamicsMod, only : fuse_patches use EDCohortDynamicsMod, only : fuse_cohorts @@ -104,7 +103,6 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) character(len=line_strlen) :: header_str ! large string for whole lines real(r8) :: age_init ! dummy value for creating a patch real(r8) :: area_init ! dummy value for creating a patch - real(r8) :: spread_init(nclmax) ! dummy value for creating a patch real(r8) :: cwd_ag_init(ncwd) ! dummy value for creating a patch real(r8) :: cwd_bg_init(ncwd) ! dummy value for creating a patch real(r8) :: leaf_litter_init(maxpft) ! dummy value for creating a patch @@ -247,13 +245,12 @@ subroutine initialize_sites_by_inventory(nsites,sites,bc_in) age_init = 0.0_r8 area_init = 0.0_r8 - spread_init(:) = ED_val_maxspread cwd_ag_init(:) = 0.0_r8 cwd_bg_init(:) = 0.0_r8 leaf_litter_init(:) = 0.0_r8 root_litter_init(:) = 0.0_r8 - call create_patch(sites(s), newpatch, age_init, area_init, spread_init, & + call create_patch(sites(s), newpatch, age_init, area_init, & cwd_ag_init, cwd_bg_init, & leaf_litter_init, root_litter_init ) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index 347bfc9052..78d165364a 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -114,7 +114,7 @@ module FatesRestartInterfaceMod integer, private :: ir_leaf_litter_in_paft integer, private :: ir_root_litter_in_paft integer, private :: ir_seed_bank_sift - integer, private :: ir_spread_pacl + integer, private :: ir_spread_si integer, private :: ir_livegrass_pa integer, private :: ir_age_pa integer, private :: ir_area_pa @@ -815,7 +815,7 @@ subroutine define_restart_vars(this, initialize_variables) call this%set_restart_var(vname='fates_spread', vtype=cohort_r8, & long_name='dynamic ratio of dbh to canopy area, by patch x canopy-layer', & units='cm/m2', flushval = flushzero, & - hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_spread_pacl ) + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_spread_si ) call this%set_restart_var(vname='fates_livegrass', vtype=cohort_r8, & long_name='total AGB from grass, by patch', & @@ -961,7 +961,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_cl ! each canopy layer class within each patch (pa_cl) integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation integer :: io_idx_si_wmem ! each water memory class within each site @@ -1046,7 +1045,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites) rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & rio_root_litter_in_paft => this%rvars(ir_root_litter_in_paft)%r81d, & rio_seed_bank_sift => this%rvars(ir_seed_bank_sift)%r81d, & - rio_spread_pacl => this%rvars(ir_spread_pacl)%r81d, & + rio_spread_si => this%rvars(ir_spread_si)%r81d, & rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & @@ -1078,7 +1077,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st @@ -1086,6 +1084,9 @@ subroutine set_restart_vectors(this,nc,nsites,sites) do i = 1,numpft rio_seed_bank_sift(io_idx_co_1st+i-1) = sites(s)%seed_bank(i) end do + + ! canopy spread term + rio_spread_si(io_idx_si) = sites(s)%spread cpatch => sites(s)%oldest_patch @@ -1200,11 +1201,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 end do - do i = 1,nclmax ! nclmax currently 2 - rio_spread_pacl(io_idx_pa_cl) = cpatch%spread(i) - io_idx_pa_cl = io_idx_pa_cl + 1 - end do - if ( DEBUG ) write(fates_log(),*) 'CLTV io_idx_pa_sunz 1 ',io_idx_pa_sunz if ( DEBUG ) write(fates_log(),*) 'CLTV 1186 ',nlevleaf,numpft,nclmax @@ -1232,7 +1228,6 @@ subroutine set_restart_vectors(this,nc,nsites,sites) ! reset counters so that they are all advanced evenly. io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_co = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st @@ -1312,7 +1307,6 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) use EDGrowthFunctionsMod, only : Dbh use EDCohortDynamicsMod, only : create_cohort use EDInitMod, only : zero_site - use EDParamsMod, only : ED_val_maxspread use EDPatchDynamicsMod, only : create_patch use EDPftvarcon, only : EDPftvarcon_inst @@ -1329,7 +1323,6 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) type(ed_cohort_type), allocatable :: temp_cohort real(r8) :: cwd_ag_local(ncwd) real(r8) :: cwd_bg_local(ncwd) - real(r8) :: spread_local(nclmax) real(r8) :: leaf_litter_local(maxpft) real(r8) :: root_litter_local(maxpft) real(r8) :: patch_age @@ -1348,7 +1341,6 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) cwd_bg_local(:) = 0.0_r8 leaf_litter_local(:) = 0.0_r8 root_litter_local(:) = 0.0_r8 - spread_local(:) = ED_val_maxspread patch_age = 0.0_r8 ! ---------------------------------------------------------------------------------- @@ -1396,7 +1388,7 @@ subroutine create_patchcohort_structure(this, nc, nsites, sites, bc_in) ! make new patch call create_patch(sites(s), newp, patch_age, area, & - spread_local, cwd_ag_local, cwd_bg_local, & + cwd_ag_local, cwd_bg_local, & leaf_litter_local, root_litter_local) newp%siteptr => sites(s) @@ -1537,7 +1529,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) integer :: io_idx_co ! cohort index integer :: io_idx_pa_pft ! each pft within each patch (pa_pft) integer :: io_idx_pa_cwd ! each cwd class within each patch (pa_cwd) - integer :: io_idx_pa_cl ! each canopy layer class within each patch (pa_cl) integer :: io_idx_pa_sunz ! index for the combined dimensions for radiation integer :: io_idx_si_wmem ! each water memory class within each site @@ -1616,7 +1607,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites) rio_leaf_litter_in_paft => this%rvars(ir_leaf_litter_in_paft)%r81d, & rio_root_litter_in_paft => this%rvars(ir_root_litter_in_paft)%r81d, & rio_seed_bank_sift => this%rvars(ir_seed_bank_sift)%r81d, & - rio_spread_pacl => this%rvars(ir_spread_pacl)%r81d, & + rio_spread_si => this%rvars(ir_spread_si)%r81d, & rio_livegrass_pa => this%rvars(ir_livegrass_pa)%r81d, & rio_age_pa => this%rvars(ir_age_pa)%r81d, & rio_area_pa => this%rvars(ir_area_pa)%r81d, & @@ -1637,7 +1628,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_co = io_idx_co_1st io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st io_idx_si_wmem = io_idx_co_1st @@ -1645,6 +1635,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites) do i = 1,numpft sites(s)%seed_bank(i) = rio_seed_bank_sift(io_idx_co_1st+i-1) enddo + + sites(s)%spread = rio_spread_si(io_idx_si) ! Perform a check on the number of patches per site patchespersite = 0 @@ -1725,7 +1717,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) cpatch%root_litter(:) = 0.0_r8 cpatch%leaf_litter_in(:) = 0.0_r8 cpatch%root_litter_in(:) = 0.0_r8 - cpatch%spread(:) = 0.0_r8 ! ! deal with patch level fields here @@ -1760,11 +1751,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_pa_cwd = io_idx_pa_cwd + 1 enddo - do i = 1,nclmax ! nclmax currently 2 - cpatch%spread(i) = rio_spread_pacl(io_idx_pa_cl) - io_idx_pa_cl = io_idx_pa_cl + 1 - enddo - if ( DEBUG ) write(fates_log(),*) 'CVTL io_idx_pa_sunz 1 ',io_idx_pa_sunz do k = 1,nlevleaf ! nlevleaf currently 40 @@ -1790,7 +1776,6 @@ subroutine get_restart_vectors(this, nc, nsites, sites) ! and max the number of allowed cohorts per patch io_idx_pa_pft = io_idx_co_1st io_idx_pa_cwd = io_idx_co_1st - io_idx_pa_cl = io_idx_co_1st io_idx_co = io_idx_co_1st io_idx_pa_sunz = io_idx_co_1st From 9b64aeb8af9918eff9e76d0d77d6e106a69fc6e7 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 26 Sep 2017 21:52:43 -0600 Subject: [PATCH 02/10] bugfixes and cleanup on prior --- biogeochem/EDCanopyStructureMod.F90 | 8 +------- biogeochem/EDGrowthFunctionsMod.F90 | 20 +++++++------------- main/EDInitMod.F90 | 4 +--- main/EDParamsMod.F90 | 1 - main/FatesHistoryInterfaceMod.F90 | 2 +- 5 files changed, 10 insertions(+), 25 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 95b0e1809d..7e4608b264 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -811,7 +811,7 @@ subroutine canopy_spread( currentSite ) ! !LOCAL VARIABLES: type (ed_cohort_type), pointer :: currentCohort type (ed_patch_type) , pointer :: currentPatch - real(r8) :: arealayer(nlevleaf) ! Amount of canopy in each layer. + real(r8) :: arealayer ! Amount of canopy in each layer. real(r8) :: inc ! Arbitrary daily incremental change in canopy area integer :: z !---------------------------------------------------------------------- @@ -843,12 +843,6 @@ subroutine canopy_spread( currentSite ) else currentSite%spread = currentSite%spread + inc endif - if(currentSite%spread >= 1._r8)then - currentSite%spread = 1._r8 - endif - if(currentSite%spread <= 0._r8)then - currentSite%spread = 0._r8 - endif ! put within bounds to make sure it stays between 0 and 1 currentSite%spread = max(min(currentSite%spread, 1._r8), 0._r8) diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 index 7abb19a19a..b2c357a2d0 100755 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -224,14 +224,13 @@ real(r8) function c_area( cohort_in ) ! Function of DBH (cm) canopy spread (m/cm) and number of individuals. ! ============================================================================ - use EDParamsMod , only : ED_val_grass_spread - use EDParamsMod , only : ED_val_maxspread, ED_val_minspread use EDTypesMod , only : nclmax type(ed_cohort_type), intent(in) :: cohort_in real(r8) :: dbh ! Tree diameter at breat height. cm. real(r8) :: crown_area_to_dbh_exponent + real(r8) :: spreadterm ! default is to use the same exponent as the dbh to bleaf exponent so that per-plant canopy depth remains invariant during growth, ! but allowed to vary via the allom_blca_expnt_diff term (which has default value of zero) @@ -243,9 +242,8 @@ real(r8) function c_area( cohort_in ) write(fates_log(),*) 'z_area 2',EDPftvarcon_inst%max_dbh write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody write(fates_log(),*) 'z_area 4',cohort_in%n - write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread + write(fates_log(),*) 'z_area 5',cohort_in%siteptr%spread write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer - write(fates_log(),*) 'z_area 7',ED_val_grass_spread end if dbh = min(cohort_in%dbh,EDPftvarcon_inst%max_dbh(cohort_in%pft)) @@ -260,15 +258,11 @@ real(r8) function c_area( cohort_in ) ! (rgk/cdk 05/2017) ! ---------------------------------------------------------------------------------- - if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then - ! apply site-level spread elasticity to the cohort crown allometry term - spreadterm = cohort_in%siteptr%spread * ED_val_maxspread + (1._r8 - cohort_in%siteptr%spread * ED_val_minspread) - ! - c_area = 3.142_r8 * cohort_in%n * & - (spreadterm * dbh)**crown_area_to_dbh_exponent - else - c_area = 3.142_r8 * cohort_in%n * (ED_val_grass_spread*dbh)**crown_area_to_dbh_exponent - end if + ! apply site-level spread elasticity to the cohort crown allometry term + spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + & + (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft) + ! + c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent end function c_area diff --git a/main/EDInitMod.F90 b/main/EDInitMod.F90 index 1b5c060d8d..1df4d1c772 100755 --- a/main/EDInitMod.F90 +++ b/main/EDInitMod.F90 @@ -122,7 +122,6 @@ subroutine set_site_properties( nsites, sites) ! !DESCRIPTION: ! ! !USES: - use EDParamsMod , only : ED_val_maxspread ! ! !ARGUMENTS @@ -141,7 +140,6 @@ subroutine set_site_properties( nsites, sites) real(r8) :: watermem integer :: dleafoff integer :: dleafon - real(r8) :: spread_local !---------------------------------------------------------------------- if ( hlm_is_restart == ifalse ) then @@ -192,7 +190,7 @@ subroutine set_site_properties( nsites, sites) sites(s)%frac_burnt = 0.0_r8 sites(s)%old_stock = 0.0_r8 - sites(s)%spread_local = 1.0_r8 + sites(s)%spread = 1.0_r8 end do return diff --git a/main/EDParamsMod.F90 b/main/EDParamsMod.F90 index e8bc151afb..d57b2b5879 100644 --- a/main/EDParamsMod.F90 +++ b/main/EDParamsMod.F90 @@ -2,7 +2,6 @@ module EDParamsMod ! ! module that deals with reading the ED parameter file ! -canopyclosure_threshold use FatesParametersInterface, only : param_string_length use FatesGlobals , only : fates_log diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 22e123df19..36d5d5db33 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -1146,7 +1146,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_trimming_pa => this%hvars(ih_trimming_pa)%r81d, & hio_area_plant_pa => this%hvars(ih_area_plant_pa)%r81d, & hio_area_treespread_pa => this%hvars(ih_area_treespread_pa)%r81d, & - hio_canopy_spread_pa => this%hvars(ih_canopy_spread_pa)%r81d, & + hio_canopy_spread_si => this%hvars(ih_canopy_spread_si)%r81d, & hio_biomass_si_pft => this%hvars(ih_biomass_si_pft)%r82d, & hio_leafbiomass_si_pft => this%hvars(ih_leafbiomass_si_pft)%r82d, & hio_storebiomass_si_pft => this%hvars(ih_storebiomass_si_pft)%r82d, & From cba6a07eb1f6eb32d6ec7c4df4a39e87f821b7b9 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 28 Sep 2017 11:43:19 -0600 Subject: [PATCH 03/10] fixed a bug in the restart code --- main/FatesRestartInterfaceMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index a93b8775a9..37f69eb030 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -843,7 +843,7 @@ subroutine define_restart_vars(this, initialize_variables) units='kgC/m2/year', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_seed_bank_sift ) - call this%set_restart_var(vname='fates_spread', vtype=cohort_r8, & + call this%set_restart_var(vname='fates_spread', vtype=site_r8, & long_name='dynamic ratio of dbh to canopy area, by patch x canopy-layer', & units='cm/m2', flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_spread_si ) From e372f3a8e45a8fd0c7175c99ed0004da8015a2ce Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 28 Sep 2017 14:23:43 -0600 Subject: [PATCH 04/10] added new variable to track total biomass within each patch age bin --- main/FatesHistoryInterfaceMod.F90 | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 5afa5cff17..682d8a400f 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -249,6 +249,7 @@ module FatesHistoryInterfaceMod integer, private :: ih_ncl_si_age integer, private :: ih_npatches_si_age integer, private :: ih_zstar_si_age + integer, private :: ih_biomass_si_age ! Indices to hydraulics variables @@ -1267,6 +1268,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_ncl_si_age => this%hvars(ih_ncl_si_age)%r82d, & hio_npatches_si_age => this%hvars(ih_npatches_si_age)%r82d, & hio_zstar_si_age => this%hvars(ih_zstar_si_age)%r82d, & + hio_biomass_si_age => this%hvars(ih_biomass_si_age)%r82d, & hio_litter_moisture_si_fuel => this%hvars(ih_litter_moisture_si_fuel)%r82d, & hio_cwd_ag_si_cwdsc => this%hvars(ih_cwd_ag_si_cwdsc)%r82d, & hio_cwd_bg_si_cwdsc => this%hvars(ih_cwd_bg_si_cwdsc)%r82d, & @@ -1397,6 +1399,11 @@ subroutine update_history_dyn(this,nc,nsites,sites) hio_biomass_si_pft(io_si, ft) = hio_biomass_si_pft(io_si, ft) + & (ccohort%n * AREA_INV) * ccohort%b * g_per_kg + ! update total biomass per age bin + hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) & + + ccohort%b * ccohort%n + + ! Site by Size-Class x PFT (SCPF) ! ------------------------------------------------------------------------ @@ -2624,6 +2631,12 @@ subroutine define_history_vars(this, initialize_variables) avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & ivar=ivar, initialize=initialize_variables, index = ih_zstar_si_age ) + call this%set_history_var(vname='BIOMASS_BY_AGE', units='m', & + long='Total Biomass within a given patch age bin (kg C)', & + use_default='inactive', & + avgflag='A', vtype=site_age_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=1, & + ivar=ivar, initialize=initialize_variables, index = ih_biomass_si_age ) + ! Fire Variables call this%set_history_var(vname='FIRE_NESTEROV_INDEX', units='none', & From 75016e101fc36f2b578aca8d9fce203f07841eb5 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 28 Sep 2017 15:12:54 -0600 Subject: [PATCH 05/10] unit fix on new biomass per patch age variable --- main/FatesHistoryInterfaceMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 682d8a400f..6a2b8ca08c 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -1401,7 +1401,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) ! update total biomass per age bin hio_biomass_si_age(io_si,cpatch%age_class) = hio_biomass_si_age(io_si,cpatch%age_class) & - + ccohort%b * ccohort%n + + ccohort%b * ccohort%n * AREA_INV ! Site by Size-Class x PFT (SCPF) From 7c06febbdd23e4254b711e668a0bf71c04fd6575 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 2 Oct 2017 15:35:45 -0600 Subject: [PATCH 06/10] removed pi-like term from the crown allometry calculation, and moving it into the allom_d2ca_coefficient parameters --- biogeochem/EDGrowthFunctionsMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 index 7845729ef5..30c30048b7 100755 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -263,7 +263,7 @@ real(r8) function c_area( cohort_in ) spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + & (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft) ! - c_area = 3.142_r8 * cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent + c_area = cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent end function c_area From 06e3c885ee88d3b298ef758fda21a998e5b3564c Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 4 Oct 2017 12:51:44 -0600 Subject: [PATCH 07/10] pulled crown area : dbh coefficient out of exponent parentheses, and changed the name of arealayer var in canopy_spread subroutine --- biogeochem/EDCanopyStructureMod.F90 | 8 ++++---- biogeochem/EDGrowthFunctionsMod.F90 | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/biogeochem/EDCanopyStructureMod.F90 b/biogeochem/EDCanopyStructureMod.F90 index 7e4608b264..626e0e9988 100755 --- a/biogeochem/EDCanopyStructureMod.F90 +++ b/biogeochem/EDCanopyStructureMod.F90 @@ -811,7 +811,7 @@ subroutine canopy_spread( currentSite ) ! !LOCAL VARIABLES: type (ed_cohort_type), pointer :: currentCohort type (ed_patch_type) , pointer :: currentPatch - real(r8) :: arealayer ! Amount of canopy in each layer. + real(r8) :: sitelevel_canopyarea ! Amount of canopy in top layer at the site level real(r8) :: inc ! Arbitrary daily incremental change in canopy area integer :: z !---------------------------------------------------------------------- @@ -820,7 +820,7 @@ subroutine canopy_spread( currentSite ) currentPatch => currentSite%oldest_patch - arealayer = 0.0_r8 + sitelevel_canopyarea = 0.0_r8 do while (associated(currentPatch)) !calculate canopy area in each patch... @@ -828,7 +828,7 @@ subroutine canopy_spread( currentSite ) do while (associated(currentCohort)) currentCohort%c_area = c_area(currentCohort) if(EDPftvarcon_inst%woody(currentCohort%pft) .eq. 1 .and. currentCohort%canopy_layer .eq. 1 ) then - arealayer = arealayer + currentCohort%c_area + sitelevel_canopyarea = sitelevel_canopyarea + currentCohort%c_area endif currentCohort => currentCohort%shorter enddo @@ -838,7 +838,7 @@ subroutine canopy_spread( currentSite ) enddo !currentPatch !If the canopy area is approaching closure, squash the tree canopies and make them taller and thinner - if( arealayer/AREA .gt. ED_val_canopy_closure_thresh ) then + if( sitelevel_canopyarea/AREA .gt. ED_val_canopy_closure_thresh ) then currentSite%spread = currentSite%spread - inc else currentSite%spread = currentSite%spread + inc diff --git a/biogeochem/EDGrowthFunctionsMod.F90 b/biogeochem/EDGrowthFunctionsMod.F90 index 30c30048b7..1de0c041e9 100755 --- a/biogeochem/EDGrowthFunctionsMod.F90 +++ b/biogeochem/EDGrowthFunctionsMod.F90 @@ -263,7 +263,7 @@ real(r8) function c_area( cohort_in ) spreadterm = cohort_in%siteptr%spread * EDPftvarcon_inst%allom_d2ca_coefficient_max(cohort_in%pft) + & (1._r8 - cohort_in%siteptr%spread) * EDPftvarcon_inst%allom_d2ca_coefficient_min(cohort_in%pft) ! - c_area = cohort_in%n * (spreadterm * dbh)**crown_area_to_dbh_exponent + c_area = cohort_in%n * spreadterm * dbh ** crown_area_to_dbh_exponent end function c_area From aea25a4118d86932c829a7061290b92942bb3ea8 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 4 Oct 2017 12:45:01 -0700 Subject: [PATCH 08/10] Turned on the NL switch controlling logging. --- biogeochem/EDLoggingMortalityMod.F90 | 3 +-- main/FatesInterfaceMod.F90 | 10 +++------- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index 7d2afb870a..f3636ad8c1 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -78,8 +78,7 @@ subroutine IsItLoggingTime(is_master,currentSite) logging_time = .false. icode = int(logging_event_code) -! if(hlm_use_logging.eq.ifalse) return ! Don't turn on until fates-clm adds - ! this to the interface (RGK 08-2017) + if(hlm_use_logging.eq.ifalse) return if(icode .eq. 1) then ! Logging is turned off diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index a7bbcb7147..0189876fb1 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -1124,11 +1124,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) hlm_use_vertsoilc = unset_int hlm_use_spitfire = unset_int hlm_use_planthydro = unset_int - hlm_use_logging = ifalse ! (RGK: to allow backwards compatibility - ! defaulting to FALSE, this will allow - ! the call from the HLM to not necessarily - ! be forced to exist. Will set this to unset - ! along with other non backwards compatible changes + hlm_use_logging = unset_int hlm_use_ed_st3 = unset_int hlm_use_ed_prescribed_phys = unset_int hlm_use_inventory_init = unset_int @@ -1172,7 +1168,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if ( .not.((hlm_use_logging .eq.1).or.(hlm_use_logging.eq.0)) ) then if (fates_global_verbose()) then - write(fates_log(), *) 'The FATES namelist planthydro flag must be 0 or 1, exiting' + write(fates_log(), *) 'The FATES namelist use_logging flag must be 0 or 1, exiting' end if call endrun(msg=errMsg(sourcefile, __LINE__)) end if @@ -1411,7 +1407,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) case('use_logging') hlm_use_logging = ival if (fates_global_verbose()) then - write(fates_log(),*) 'Transfering hlm_use_planthydro= ',ival,' to FATES' + write(fates_log(),*) 'Transfering hlm_use_logging= ',ival,' to FATES' end if case('use_ed_st3') From 2d63622e9b513c5731719d166119a1dc9becc264 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 4 Oct 2017 12:55:09 -0700 Subject: [PATCH 09/10] logging NL switch, syntax fixes. --- biogeochem/EDLoggingMortalityMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/biogeochem/EDLoggingMortalityMod.F90 b/biogeochem/EDLoggingMortalityMod.F90 index f3636ad8c1..4f87fcaaed 100644 --- a/biogeochem/EDLoggingMortalityMod.F90 +++ b/biogeochem/EDLoggingMortalityMod.F90 @@ -36,6 +36,7 @@ module EDLoggingMortalityMod use FatesInterfaceMod , only : hlm_model_day use FatesInterfaceMod , only : hlm_day_of_year use FatesInterfaceMod , only : hlm_days_per_year + use FatesInterfaceMod , only : hlm_use_logging use FatesConstantsMod , only : itrue,ifalse use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log From 4f6a3429b2d316b3a7b2ab854931af9f1055df6c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 5 Oct 2017 15:57:04 -0700 Subject: [PATCH 10/10] bug fixes: dimension structure incorrectly defined as a class instead of type, unused hvar pointers and weird association statement. --- main/FatesHistoryInterfaceMod.F90 | 19 ++++++------------- main/FatesHistoryVariableType.F90 | 9 ++++++--- main/FatesRestartVariableType.F90 | 2 +- 3 files changed, 13 insertions(+), 17 deletions(-) diff --git a/main/FatesHistoryInterfaceMod.F90 b/main/FatesHistoryInterfaceMod.F90 index 5afa5cff17..08828774fb 100644 --- a/main/FatesHistoryInterfaceMod.F90 +++ b/main/FatesHistoryInterfaceMod.F90 @@ -869,20 +869,17 @@ subroutine flush_hvars(this,nc,upfreq_in) class(fates_history_interface_type) :: this integer,intent(in) :: nc integer,intent(in) :: upfreq_in - integer :: ivar - type(fates_history_variable_type),pointer :: hvar integer :: lb1,ub1,lb2,ub2 do ivar=1,ubound(this%hvars,1) - associate( hvar => this%hvars(ivar) ) - if (hvar%upfreq == upfreq_in) then ! Only flush variables with update on dynamics step - call hvar%Flush(nc, this%dim_bounds, this%dim_kinds) - end if - end associate + if (this%hvars(ivar)%upfreq == upfreq_in) then ! Only flush variables with update on dynamics step + call this%hvars(ivar)%flush(nc, this%dim_bounds, this%dim_kinds) + + end if end do - end subroutine flush_hvars +end subroutine flush_hvars ! ===================================================================================== @@ -915,7 +912,6 @@ subroutine set_history_var(this, vname, units, long, use_default, avgflag, vtype ! not used ! locals - type(fates_history_variable_type), pointer :: hvar integer :: ub1, lb1, ub2, lb2 ! Bounds for allocating the var integer :: ityp @@ -1110,7 +1106,7 @@ subroutine update_history_dyn(this,nc,nsites,sites) use EDTypesMod , only : nlevleaf ! Arguments - class(fates_history_interface_type) :: this + class(fates_history_interface_type) :: this integer , intent(in) :: nc ! clump index integer , intent(in) :: nsites type(ed_site_type) , intent(inout), target :: sites(nsites) @@ -1135,7 +1131,6 @@ subroutine update_history_dyn(this,nc,nsites,sites) real(r8) :: patch_scaling_scalar ! ratio of canopy to patch area for counteracting patch scaling real(r8) :: dbh ! diameter ("at breast height") - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort @@ -1888,7 +1883,6 @@ subroutine update_history_prod(this,nc,nsites,sites,dt_tstep) real(r8), parameter :: tiny = 1.e-5_r8 ! some small number integer :: ipa2 ! patch incrementer integer :: cnlfpft_indx, cnlf_indx, ipft, ican, ileaf ! more iterators and indices - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort real(r8) :: per_dt_tstep ! Time step in frequency units (/s) @@ -2246,7 +2240,6 @@ subroutine update_history_hydraulics(this,nc,nsites,sites,dt_tstep) integer :: ipa2 ! patch incrementer integer :: iscpf ! index of the scpf group - type(fates_history_variable_type),pointer :: hvar type(ed_patch_type),pointer :: cpatch type(ed_cohort_type),pointer :: ccohort type(ed_cohort_hydr_type), pointer :: ccohort_hydr diff --git a/main/FatesHistoryVariableType.F90 b/main/FatesHistoryVariableType.F90 index e76531062c..86b9a4f875 100644 --- a/main/FatesHistoryVariableType.F90 +++ b/main/FatesHistoryVariableType.F90 @@ -89,7 +89,7 @@ subroutine Init(this, vname, units, long, use_default, & call dim_kinds(dk_index)%set_active() call this%GetBounds(0, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) - + ! NOTE(rgk, 2016-09) currently, all array spaces are flushed each ! time the update is called. The flush here on the initialization ! may be redundant, but will prevent issues in the future if we @@ -167,7 +167,7 @@ subroutine Init(this, vname, units, long, use_default, & end subroutine Init ! ===================================================================================== - + subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) use FatesIODimensionsMod, only : fates_io_dimension_type @@ -176,7 +176,7 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) class(fates_history_variable_type), intent(inout) :: this integer, intent(in) :: thread - class(fates_io_dimension_type), intent(in) :: dim_bounds(:) + type(fates_io_dimension_type), intent(in) :: dim_bounds(:) type(fates_io_variable_kind_type), intent(in) :: dim_kinds(:) integer, intent(out) :: lb1 integer, intent(out) :: ub1 @@ -205,14 +205,17 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) ub2 = dim_bounds(d_index)%upper_bound end if else + d_index = dim_kinds(this%dim_kinds_index)%dim1_index lb1 = dim_bounds(d_index)%clump_lower_bound(thread) ub1 = dim_bounds(d_index)%clump_upper_bound(thread) + if(ndims>1)then d_index = dim_kinds(this%dim_kinds_index)%dim2_index lb2 = dim_bounds(d_index)%clump_lower_bound(thread) ub2 = dim_bounds(d_index)%clump_upper_bound(thread) end if + end if end subroutine GetBounds diff --git a/main/FatesRestartVariableType.F90 b/main/FatesRestartVariableType.F90 index 40648fb4c6..48152adf7f 100644 --- a/main/FatesRestartVariableType.F90 +++ b/main/FatesRestartVariableType.F90 @@ -118,7 +118,7 @@ subroutine GetBounds(this, thread, dim_bounds, dim_kinds, lb1, ub1, lb2, ub2) class(fates_restart_variable_type), intent(inout) :: this integer, intent(in) :: thread - class(fates_io_dimension_type), intent(in) :: dim_bounds(:) + type(fates_io_dimension_type), intent(in) :: dim_bounds(:) type(fates_io_variable_kind_type), intent(in) :: dim_kinds(:) integer, intent(out) :: lb1 integer, intent(out) :: ub1