diff --git a/main/ChecksBalancesMod.F90 b/main/ChecksBalancesMod.F90 index 325a1089e8..f041389fd3 100644 --- a/main/ChecksBalancesMod.F90 +++ b/main/ChecksBalancesMod.F90 @@ -302,7 +302,7 @@ subroutine CheckIntegratedMassPools(site) select case(element_list(el)) case(carbon12_element) - net_uptake = diag%npp + site_mass%net_root_uptake*area_inv + net_uptake = (site_mass%gpp_acc - site_mass%aresp_acc + site_mass%net_root_uptake)*area_inv case(nitrogen_element) net_uptake = site_mass%net_root_uptake*area_inv case(phosphorus_element) diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 7a5bd21840..8d75256e50 100644 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -198,7 +198,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) call ZeroBCOutCarbonFluxes(bc_out) ! Zero mass balance - call TotalBalanceCheck(currentSite, 0) + call TotalBalanceCheck(currentSite, 0, is_restarting=.false.) ! We do not allow phenology while in ST3 mode either, it is hypothetically ! possible to allow this, but we have not plugged in the litter fluxes @@ -263,7 +263,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) currentPatch => currentPatch%younger enddo - call TotalBalanceCheck(currentSite,1) + call TotalBalanceCheck(currentSite,1,is_restarting=.false.) currentPatch => currentSite%oldest_patch do while (associated(currentPatch)) @@ -286,7 +286,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) end if - call TotalBalanceCheck(currentSite,2) + call TotalBalanceCheck(currentSite,2,is_restarting=.false.) !********************************************************************************* ! Patch dynamics sub-routines: fusion, new patch creation (spwaning), termination. @@ -304,7 +304,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) call spawn_patches(currentSite, bc_in) - call TotalBalanceCheck(currentSite,3) + call TotalBalanceCheck(currentSite,3,is_restarting=.false.) ! fuse on the spawned patches. call fuse_patches(currentSite, bc_in ) @@ -319,14 +319,14 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, bc_out) end if ! SP has changes in leaf carbon but we don't expect them to be in balance. - call TotalBalanceCheck(currentSite,4) + call TotalBalanceCheck(currentSite,4,is_restarting=.false.) ! kill patches that are too small call terminate_patches(currentSite, bc_in) end if ! Final instantaneous mass balance check - call TotalBalanceCheck(currentSite,5) + call TotalBalanceCheck(currentSite,5,is_restarting=.false.) end subroutine ed_ecosystem_dynamics @@ -668,11 +668,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentSite%mass_balance(element_pos(carbon12_element))%net_root_uptake - & currentCohort%daily_c_efflux*currentCohort%n - ! Save NPP diagnostic for flux accounting [kg/m2/day] - - currentSite%flux_diags%npp = currentSite%flux_diags%npp + & - currentCohort%npp_acc_hold/real( hlm_days_per_year,r8) * currentCohort%n * area_inv - ! And simultaneously add the input fluxes to mass balance accounting site_cmass%gpp_acc = site_cmass%gpp_acc + & currentCohort%gpp_acc * currentCohort%n @@ -682,8 +677,6 @@ subroutine ed_integrate_state_variables(currentSite, bc_in, bc_out ) currentCohort%resp_excess_hold*currentCohort%n + & currentCohort%resp_g_acc_hold*currentCohort%n/real( hlm_days_per_year,r8) - - call currentCohort%prt%CheckMassConservation(ft,5) ! Update the leaf biophysical rates based on proportion of leaf @@ -848,10 +841,7 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) call set_patchno(currentSite,.true.,1) end if - ! Pass site-level mass fluxes to output boundary conditions - ! [kg/site/day] * [site/m2 day/sec] = [kgC/m2/s] - bc_out%gpp_site = site_cmass%gpp_acc * area_inv / sec_per_day - bc_out%ar_site = site_cmass%aresp_acc * area_inv / sec_per_day + if(hlm_use_sp.eq.ifalse .and. (.not.is_restarting))then call canopy_spread(currentSite) @@ -860,13 +850,13 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) site_cmass%aresp_acc = 0._r8 end if - call TotalBalanceCheck(currentSite,6) + call TotalBalanceCheck(currentSite,6,is_restarting=is_restarting) if(hlm_use_sp.eq.ifalse .and. (.not.is_restarting) )then call canopy_structure(currentSite, bc_in) endif - call TotalBalanceCheck(currentSite,final_check_id) + call TotalBalanceCheck(currentSite,final_check_id,is_restarting=is_restarting) ! Update recruit L2FRs based on new canopy position call SetRecruitL2FR(currentSite) @@ -927,26 +917,34 @@ subroutine ed_update_site( currentSite, bc_in, bc_out, is_restarting ) ! Set boundary condition to HLM for carbon loss to atm from fires and grazing ! [kgC/ha/day]*[ha/m2]*[day/s] = [kg/m2/s] - bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * ha_per_m2 * days_per_sec - bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * ha_per_m2 * days_per_sec - - + bc_out%fire_closs_to_atm_si = site_cmass%burn_flux_to_atm * area_inv * days_per_sec + bc_out%grazing_closs_to_atm_si = site_cmass%herbivory_flux_out * area_inv * days_per_sec + bc_out%gpp_site = site_cmass%gpp_acc * area_inv * days_per_sec + bc_out%ar_site = site_cmass%aresp_acc * area_inv * days_per_sec end subroutine ed_update_site !-------------------------------------------------------------------------------! - subroutine TotalBalanceCheck (currentSite, call_index ) + subroutine TotalBalanceCheck (currentSite, call_index, is_restarting ) ! ! !DESCRIPTION: ! This routine looks at the mass flux in and out of the FATES and compares it to ! the change in total stocks (states). ! Fluxes in are NPP. Fluxes out are decay of CWD and litter into SOM pools. + ! Note: If the model is restarting, it is assumed that the mass stocks + ! that were saved in the restart file are the "old" stocks, and they + ! should equal the stocks that are currently in the sites. However, + ! the fluxes on a restart are saved so that they can inform the HLM + ! on the next day, so we have to modify our mass balance check to ignore + ! fluxes on restarts. ! ! !ARGUMENTS: type(ed_site_type) , intent(inout) :: currentSite integer , intent(in) :: call_index + logical , intent(in) :: is_restarting + ! ! !LOCAL VARIABLES: type(site_massbal_type),pointer :: site_mass @@ -987,32 +985,35 @@ subroutine TotalBalanceCheck (currentSite, call_index ) change_in_stock = 0.0_r8 - ! Loop through the number of elements in the system - do el = 1, num_elements + do_elem_loop: do el = 1, num_elements site_mass => currentSite%mass_balance(el) call SiteMassStock(currentSite,el,total_stock,biomass_stock,litter_stock,seed_stock) change_in_stock = total_stock - site_mass%old_stock - - flux_in = site_mass%seed_in + & - site_mass%net_root_uptake + & - site_mass%gpp_acc + & - site_mass%flux_generic_in + & - site_mass%patch_resize_err - - flux_out = sum(site_mass%wood_product_harvest(:)) + & - sum(site_mass%wood_product_landusechange(:)) + & - site_mass%burn_flux_to_atm + & - site_mass%seed_out + & - site_mass%flux_generic_out + & - site_mass%frag_out + & - site_mass%aresp_acc + & - site_mass%herbivory_flux_out - + if(is_restarting) then + flux_in = 0._r8 + flux_out = 0._r8 + else + flux_in = site_mass%seed_in + & + site_mass%net_root_uptake + & + site_mass%gpp_acc + & + site_mass%flux_generic_in + & + site_mass%patch_resize_err + + flux_out = sum(site_mass%wood_product_harvest(:)) + & + sum(site_mass%wood_product_landusechange(:)) + & + site_mass%burn_flux_to_atm + & + site_mass%seed_out + & + site_mass%flux_generic_out + & + site_mass%frag_out + & + site_mass%aresp_acc + & + site_mass%herbivory_flux_out + end if + net_flux = flux_in - flux_out error = abs(net_flux - change_in_stock) @@ -1110,12 +1111,13 @@ subroutine TotalBalanceCheck (currentSite, call_index ) ! This is the last check of the sequence, where we update our total ! error check and the final fates stock - if(call_index == final_check_id) then - site_mass%old_stock = total_stock - site_mass%err_fates = net_flux - change_in_stock + if(call_index == final_check_id .and. .not.is_restarting) then + site_mass%old_stock = total_stock + site_mass%err_fates = net_flux - change_in_stock end if - end do + end do do_elem_loop + end if ! not SP mode end subroutine TotalBalanceCheck diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 90be4df5ec..46c65218b1 100644 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -217,18 +217,15 @@ module EDTypesMod type, public :: site_fluxdiags_type - ! This is for all diagnostics that are uniform over all elements (C,N,P) + ! These are site level flux diagnostics that are not used + ! in mass balance checks. We use these structures + ! to inform the history output. These values are not + ! zero'd when dynamics are completed. These values + ! are zero'd on cold-starts, and on restarts prior to the read + ! This is for all diagnostics that are uniform over all elements (C,N,P) type(elem_diag_type), pointer :: elem(:) - ! This variable is slated as to-do, but the fluxdiags type needs - ! to be refactored first. Currently this type is allocated - ! by chemical species (ie C, N or P). GPP is C, but not N or P (RGK 0524) - ! Previous day GPP [kgC/m2/year], partitioned by size x pft - !real(r8),allocatable :: gpp_prev_scpf(:) - - real(r8) :: npp ! kg m-2 day-1 - ! Nutrient Flux Diagnostics real(r8) :: resp_excess ! plant carbon respired due to carbon overflow @@ -678,7 +675,6 @@ subroutine ZeroFluxDiags(this) end do - this%npp = 0._r8 this%resp_excess = 0._r8 this%nh4_uptake = 0._r8 this%no3_uptake = 0._r8 @@ -694,11 +690,6 @@ subroutine ZeroFluxDiags(this) this%p_uptake_scpf(:) = 0._r8 this%p_efflux_scpf(:) = 0._r8 - ! We don't zero gpp_prev_scpf because this is not - ! incremented like others, it is assigned at the end - ! of the daily history write process - - return end subroutine ZeroFluxDiags diff --git a/main/FatesRestartInterfaceMod.F90 b/main/FatesRestartInterfaceMod.F90 index fb8deaced4..82835aa27d 100644 --- a/main/FatesRestartInterfaceMod.F90 +++ b/main/FatesRestartInterfaceMod.F90 @@ -114,6 +114,8 @@ module FatesRestartInterfaceMod integer :: ir_landuse_config_si integer :: ir_gpp_acc_si integer :: ir_aresp_acc_si + integer :: ir_herbivory_flux_out_si + integer :: ir_burn_flux_to_atm_si integer :: ir_ncohort_pa integer :: ir_canopy_layer_co @@ -1199,6 +1201,15 @@ subroutine define_restart_vars(this, initialize_variables) units='kg/ha', veclength=num_elements, flushval = flushzero, & hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_errfates_mbal) + call this%RegisterCohortVector(symbol_base='herbivory_flux_out', vtype=site_r8, & + long_name_base='Mass flux of herbivory losses at the site level', & + units='kg/ha/day', veclength=num_elements, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_herbivory_flux_out_si) + + call this%RegisterCohortVector(symbol_base='burn_flux_to_atm', vtype=site_r8, & + long_name_base='Mass flux of burn loss to the atmosphere at site level', & + units='kg/ha/day', veclength=num_elements, flushval = flushzero, & + hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_burn_flux_to_atm_si) ! Time integrated mass balance accounting [kg/m2] call this%RegisterCohortVector(symbol_base='fates_liveveg_intflux', vtype=site_r8, & @@ -2543,7 +2554,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites) do i_lu_donor = 1, n_landuse_cats do i_lu_receiver = 1, n_landuse_cats do i_dist = 1, n_dist_types - rio_disturbance_rates_siluludi(io_idx_si_luludi) = sites(s)%disturbance_rates(i_dist,i_lu_donor, i_lu_receiver) + rio_disturbance_rates_siluludi(io_idx_si_luludi) = & + sites(s)%disturbance_rates(i_dist,i_lu_donor, i_lu_receiver) io_idx_si_luludi = io_idx_si_luludi + 1 end do end do @@ -2557,19 +2569,30 @@ subroutine set_restart_vectors(this,nc,nsites,sites) io_idx_si_scpf = io_idx_co_1st do i_cwd=1,ncwd - this%rvars(ir_cwdagin_flxdg+el-1)%r81d(io_idx_si_cwd) = sites(s)%flux_diags%elem(el)%cwd_ag_input(i_cwd) - this%rvars(ir_cwdbgin_flxdg+el-1)%r81d(io_idx_si_cwd) = sites(s)%flux_diags%elem(el)%cwd_bg_input(i_cwd) + this%rvars(ir_cwdagin_flxdg+el-1)%r81d(io_idx_si_cwd) = & + sites(s)%flux_diags%elem(el)%cwd_ag_input(i_cwd) + this%rvars(ir_cwdbgin_flxdg+el-1)%r81d(io_idx_si_cwd) = & + sites(s)%flux_diags%elem(el)%cwd_bg_input(i_cwd) io_idx_si_cwd = io_idx_si_cwd + 1 end do do i_pft=1,numpft - this%rvars(ir_leaflittin_flxdg+el-1)%r81d(io_idx_si_pft) = sites(s)%flux_diags%elem(el)%surf_fine_litter_input(i_pft) - this%rvars(ir_rootlittin_flxdg+el-1)%r81d(io_idx_si_pft) = sites(s)%flux_diags%elem(el)%root_litter_input(i_pft) - this%rvars(ir_woodprod_harvest_mbal+el-1)%r81d(io_idx_si_pft) = sites(s)%mass_balance(el)%wood_product_harvest(i_pft) - this%rvars(ir_woodprod_landusechange_mbal+el-1)%r81d(io_idx_si_pft) = sites(s)%mass_balance(el)%wood_product_landusechange(i_pft) + this%rvars(ir_leaflittin_flxdg+el-1)%r81d(io_idx_si_pft) = & + sites(s)%flux_diags%elem(el)%surf_fine_litter_input(i_pft) + this%rvars(ir_rootlittin_flxdg+el-1)%r81d(io_idx_si_pft) = & + sites(s)%flux_diags%elem(el)%root_litter_input(i_pft) + this%rvars(ir_woodprod_harvest_mbal+el-1)%r81d(io_idx_si_pft) = & + sites(s)%mass_balance(el)%wood_product_harvest(i_pft) + this%rvars(ir_woodprod_landusechange_mbal+el-1)%r81d(io_idx_si_pft) = & + sites(s)%mass_balance(el)%wood_product_landusechange(i_pft) io_idx_si_pft = io_idx_si_pft + 1 end do + this%rvars(ir_herbivory_flux_out_si+el-1)%r81d(io_idx_si) = & + sites(s)%mass_balance(el)%herbivory_flux_out + this%rvars(ir_burn_flux_to_atm_si+el-1)%r81d(io_idx_si) = & + sites(s)%mass_balance(el)%burn_flux_to_atm + this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%old_stock this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si) = sites(s)%mass_balance(el)%err_fates @@ -3607,6 +3630,11 @@ subroutine get_restart_vectors(this, nc, nsites, sites) io_idx_si_pft = io_idx_si_pft + 1 end do + sites(s)%mass_balance(el)%herbivory_flux_out = & + this%rvars(ir_herbivory_flux_out_si+el-1)%r81d(io_idx_si) + sites(s)%mass_balance(el)%burn_flux_to_atm = & + this%rvars(ir_burn_flux_to_atm_si+el-1)%r81d(io_idx_si) + sites(s)%mass_balance(el)%old_stock = this%rvars(ir_oldstock_mbal+el-1)%r81d(io_idx_si) sites(s)%mass_balance(el)%err_fates = this%rvars(ir_errfates_mbal+el-1)%r81d(io_idx_si)