diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index 288c8c9dd9..5b383735ea 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -733,7 +733,7 @@ subroutine fuse_cohorts(patchptr) nextc%n*nextc%year_net_uptake(i))/newn endif enddo - + currentCohort%n = newn !remove fused cohort from the list nextc%taller%shorter => nextnextc diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 5bfdf1c71a..58c4656c06 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -200,8 +200,11 @@ subroutine spawn_patches( currentSite ) ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. currentPatch => currentSite%youngest_patch + + ! zero site-level fire fluxes currentSite%cwd_ag_burned = 0.0_r8 currentSite%leaf_litter_burned = 0.0_r8 + currentSite%total_burn_flux_to_atm = 0.0_r8 site_areadis = 0.0_r8 do while(associated(currentPatch)) @@ -548,12 +551,14 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) burned_litter = new_patch%cwd_ag(c) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(c+1) !kG/m2/day new_patch%cwd_ag(c) = new_patch%cwd_ag(c) - burned_litter currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day enddo do p = 1,numpft_ed burned_litter = new_patch%leaf_litter(p) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(dg_sf) new_patch%leaf_litter(p) = new_patch%leaf_litter(p) - burned_litter - currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/dat + currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day enddo !************************************/ @@ -615,6 +620,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) SF_val_CWD_frac(c) * bstem * currentCohort%cfa currentSite%flux_out = currentSite%flux_out + dead_tree_density * & AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + dead_tree_density * & + AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa enddo @@ -625,6 +632,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) dead_tree_density * currentCohort%bl * currentCohort%cfa currentSite%flux_out = currentSite%flux_out + & dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + & + dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa enddo @@ -655,6 +664,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) !KgC/gridcell/day currentSite%flux_out = currentSite%flux_out + burned_leaves * currentCohort%n * & patch_site_areadis/currentPatch%area * AREA + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm+ burned_leaves * currentCohort%n * & + patch_site_areadis/currentPatch%area * AREA endif currentCohort%cfa = 0.0_r8 @@ -1113,7 +1124,7 @@ subroutine fuse_patches( csite ) if(nopatches > maxpatch)then iterate = 1 profiletol = profiletol * 1.1_r8 - write(iulog,*) 'maxpatch exceeded, triggering patch fusion iteration.',profiletol,nopatches + !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! !---------------------------------------------------------------------! diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index 8de0ecc9ad..846f8c11fc 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -125,15 +125,14 @@ subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, wa currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - currentPatch%root_litter_out(p) enddo - currentPatch%leaf_litter_in(:) = 0.0_r8 - currentPatch%root_litter_in(:) = 0.0_r8 - currentPatch%leaf_litter_out(:) = 0.0_r8 - currentPatch%root_litter_out(:) = 0.0_r8 - - currentPatch%CWD_AG_in(:) = 0.0_r8 - currentPatch%cwd_bg_in(:) = 0.0_r8 - currentPatch%CWD_AG_out(:) = 0.0_r8 - currentPatch%cwd_bg_out(:) = 0.0_r8 + ! currentPatch%leaf_litter_in(:) = 0.0_r8 + ! currentPatch%root_litter_in(:) = 0.0_r8 + ! currentPatch%leaf_litter_out(:) = 0.0_r8 + ! currentPatch%root_litter_out(:) = 0.0_r8 + ! currentPatch%CWD_AG_in(:) = 0.0_r8 + ! currentPatch%cwd_bg_in(:) = 0.0_r8 + ! currentPatch%CWD_AG_out(:) = 0.0_r8 + ! currentPatch%cwd_bg_out(:) = 0.0_r8 end subroutine non_canopy_derivs @@ -636,6 +635,8 @@ subroutine seeds_in( cp_pnt ) currentSite => currentPatch%siteptr currentPatch%seeds_in(:) = 0.0_r8 + currentPatch%seed_rain_flux(:) = 0.0_r8 + currentCohort => currentPatch%tallest do while (associated(currentCohort)) p = currentCohort%pft @@ -649,6 +650,7 @@ subroutine seeds_in( cp_pnt ) if (EXTERNAL_RECRUITMENT == 1) then !external seed rain - needed to prevent extinction do p = 1,numpft_ed currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + EDecophyscon%seed_rain(p) !KgC/m2/year + currentPatch%seed_rain_flux(p) = currentPatch%seed_rain_flux(p) + EDecophyscon%seed_rain(p) !KgC/m2/year enddo endif currentPatch => currentPatch%younger @@ -1251,4 +1253,6 @@ subroutine cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_i end subroutine cwd_out + + end module EDPhysiologyMod diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 02fd161a06..d601d5bf49 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -857,10 +857,10 @@ subroutine ED_Norman_Radiation (bounds, & enddo enddo if (lai_change(1,2,1).gt.0.0.and.lai_change(1,2,2).gt.0.0)then - write(iulog,*) 'lai_change(1,2,12)',lai_change(1,2,1:4) +! write(iulog,*) 'lai_change(1,2,12)',lai_change(1,2,1:4) endif if (lai_change(1,2,2).gt.0.0.and.lai_change(1,2,3).gt.0.0)then - write(iulog,*) ' lai_change (1,2,23)',lai_change(1,2,1:4) +! write(iulog,*) ' lai_change (1,2,23)',lai_change(1,2,1:4) endif if (lai_change(1,1,3).gt.0.0.and.lai_change(1,1,2).gt.0.0)then ! NO-OP diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 5257751e2d..53b79e1f8a 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -10,15 +10,26 @@ module EDCLMLinkMod use decompMod , only : bounds_type use clm_varpar , only : nclmax, nlevcan_ed, numpft, numcft, mxpft use clm_varctl , only : iulog - use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, nlevsclass_ed - use EDtypesMod , only : sclass_ed,nlevsclass_ed,AREA + + use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, ncwd + use EDtypesMod , only : sclass_ed, nlevsclass_ed, AREA + use CanopyStateType , only : canopystate_type + use clm_varctl , only : use_vertsoilc use EDParamsMod , only : ED_val_ag_biomass + use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type + use SoilBiogeochemCarbonStatetype , only : soilbiogeochem_carbonstate_type + use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep + use shr_const_mod, only: SHR_CONST_CDAY ! implicit none private ! logical :: DEBUG = .false. ! for debugging this module (EDCLMLinkMod.F90) + + ! !PUBLIC DATA MEMBERS + real(r8), public :: cwd_fcel_ed + real(r8), public :: cwd_flig_ed type, public :: ed_clm_type @@ -66,14 +77,14 @@ module EDCLMLinkMod real(r8), pointer, private :: ED_bleaf_patch (:) ! kGC/m2 Total leaf biomass. real(r8), pointer, private :: ED_biomass_patch (:) ! kGC/m2 Total biomass. - real(r8), pointer, private :: storvegc_patch (:) ! (gC/m2) stored vegetation carbon, excluding cpool - real(r8), pointer, private :: dispvegc_patch (:) ! (gC/m2) displayed veg carbon, excluding storage and cpool - real(r8), pointer, private :: leafc_patch (:) ! (gC/m2) leaf C - real(r8), pointer, private :: livestemc_patch (:) ! (gC/m2) live stem C - real(r8), pointer, private :: deadstemc_patch (:) ! (gC/m2) dead stem C - real(r8), pointer, private :: livestemn_patch (:) ! (gN/m2) live stem N - real(r8), pointer, private :: npp_patch (:) ! (gC/m2/s) patch net primary production + ! real(r8), pointer, private :: storvegc_patch (:) ! (gC/m2) stored vegetation carbon, excluding cpool + ! real(r8), pointer, private :: dispvegc_patch (:) ! (gC/m2) displayed veg carbon, excluding storage and cpool + ! real(r8), pointer, private :: leafc_patch (:) ! (gC/m2) leaf C + ! real(r8), pointer, private :: livestemc_patch (:) ! (gC/m2) live stem C + ! real(r8), pointer, private :: deadstemc_patch (:) ! (gC/m2) dead stem C + ! real(r8), pointer, private :: livestemn_patch (:) ! (gN/m2) live stem N + real(r8), pointer, private :: npp_patch (:) ! (gC/m2/s) patch net primary production real(r8), pointer, private :: gpp_patch (:) ! (gC/m2/s) patch gross primary production real(r8), pointer :: ed_gpp_gd_scpf (:,:) ! [kg/m2/yr] gross primary production @@ -106,6 +117,45 @@ module EDCLMLinkMod real(r8), pointer :: ed_m4_gd_scpf (:,:) ! [Stems/ha/yr] Mean Impact Mortality real(r8), pointer :: ed_m5_gd_scpf (:,:) ! [Stems/ha/yr] Mean Fire Mortality + ! litterfall fluxes of C from ED patches to BGC columns + + real(r8), pointer, public :: ED_c_to_litr_lab_c_col(:,:) !total labile litter coming from ED. gC/m3/s + real(r8), pointer, public :: ED_c_to_litr_cel_c_col(:,:) !total cellulose litter coming from ED. gC/m3/s + real(r8), pointer, public :: ED_c_to_litr_lig_c_col(:,:) !total lignin litter coming from ED. gC/m3/s + real(r8), pointer, private :: leaf_prof_col(:,:) !(1/m) profile of leaves + real(r8), pointer, private :: froot_prof_col(:,:,:) !(1/m) profile of fine roots + real(r8), pointer, private :: croot_prof_col(:,:) !(1/m) profile of coarse roots + real(r8), pointer, private :: stem_prof_col(:,:) !(1/m) profile of leaves + + ! summary carbon fluxes at the column level + real(r8), pointer, public :: nep_col(:) ! [gC/m2/s] Net ecosystem production, i.e. fast-timescale carbon balance that does not include disturbance + real(r8), pointer, public :: nep_timeintegrated_col(:) ! [gC/m2/s] Net ecosystem production, i.e. fast-timescale carbon balance that does not include disturbance + real(r8), pointer, public :: npp_timeintegrated_col(:) ! [gC/m2/s] Net primary production, time integrated at column level for carbon balance checking + real(r8), pointer, public :: hr_timeintegrated_col(:) ! [gC/m2/s] heterotrophic respiration, time integrated for carbon balance checking + real(r8), pointer, public :: nbp_col(:) ! [gC/m2/s] Net biosphere production, i.e. slow-timescale carbon balance that integrates to total carbon change + real(r8), pointer, public :: npp_hifreq_col(:) ! [gC/m2/s] Net primary production at the fast timescale, aggregated to the column level + real(r8), pointer, public :: fire_c_to_atm_col(:) ! [gC/m2/s] total fire carbon loss to atmosphere + real(r8), pointer, public :: ed_to_bgc_this_edts_col(:) ! [gC/m2/s] total flux of carbon from ED to BGC models on current ED timestep + real(r8), pointer, public :: ed_to_bgc_last_edts_col(:) ! [gC/m2/s] total flux of carbon from ED to BGC models on prior ED timestep + real(r8), pointer, public :: seed_rain_flux_col(:) ! [gC/m2/s] total flux of carbon from seed rain + + ! summary carbon states at the column level + real(r8), pointer, public :: totecosysc_col(:) ! [gC/m2] Total ecosystem carbon at the column level, including vegetation, CWD, litter, and soil pools + real(r8), pointer, public :: totecosysc_old_col(:) ! [gC/m2] Total ecosystem C at the column level from last call to balance check + real(r8), pointer, public :: totedc_col(:) ! [gC/m2] Total ED carbon at the column level, including vegetation, CWD, seeds, and ED litter + real(r8), pointer, public :: totedc_old_col(:) ! [gC/m2] Total ED C at the column level from last call to balance check + real(r8), pointer, public :: totbgcc_col(:) ! [gC/m2] Total BGC carbon at the column level, including litter, and soil pools + real(r8), pointer, public :: totbgcc_old_col(:) ! [gC/m2] Total BGC C at the column level from last call to balance check + real(r8), pointer, public :: biomass_stock_col(:) ! [gC/m2] total biomass at the column level in gC / m2 + real(r8), pointer, public :: ed_litter_stock_col(:) ! [gC/m2] ED litter at the column level in gC / m2 + real(r8), pointer, public :: cwd_stock_col(:) ! [gC/m2] ED CWD at the column level in gC / m2 + real(r8), pointer, public :: seed_stock_col(:) ! [gC/m2] ED seed mass carbon at the column level in gC / m2 + + ! carbon balance errors. at some point we'll reduce these to close to zero and delete, but for now we'll just keep[ track of them + real(r8), pointer, public :: cbalance_error_ed_col(:) ! [gC/m2/s] total carbon balance error for the ED side + real(r8), pointer, public :: cbalance_error_bgc_col(:) ! [gC/m2/s] total carbon balance error for the BGC side + real(r8), pointer, public :: cbalance_error_total_col(:) ! [gC/m2/s] total carbon balance error for the whole thing + contains ! Public routines @@ -113,13 +163,16 @@ module EDCLMLinkMod procedure , public :: Restart procedure , public :: SetValues procedure , public :: ed_clm_link + procedure , public :: Summary + procedure , public :: ED_BGC_Carbon_Balancecheck ! Private routines procedure , private :: ed_clm_leaf_area_profile procedure , private :: ed_update_history_variables procedure , private :: InitAllocate procedure , private :: InitHistory - procedure , private :: InitCold +! procedure , private :: InitCold + procedure , private :: flux_into_litter_pools end type ed_clm_type @@ -141,7 +194,7 @@ subroutine Init(this, bounds) call this%InitAllocate(bounds) call this%InitHistory(bounds) - call this%InitCold(bounds) + !call this%InitCold(bounds) end subroutine Init @@ -150,7 +203,8 @@ subroutine InitAllocate(this, bounds) ! ! !USES: use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : nlevgrnd + use clm_varpar , only : nlevgrnd, nlevdecomp_full + use EDtypesMod , only : numpft_ed ! ! !ARGUMENTS: class (ed_clm_type) :: this @@ -164,6 +218,7 @@ subroutine InitAllocate(this, bounds) begp = bounds%begp; endp = bounds%endp begc = bounds%begc; endc = bounds%endc + begg = bounds%begg; endg = bounds%endg allocate(this%trimming_patch (begp:endp)) ; this%trimming_patch (:) = 0.0_r8 @@ -200,17 +255,53 @@ subroutine InitAllocate(this, bounds) allocate(this%ED_bleaf_patch (begp:endp)) ; this%ED_bleaf_patch (:) = 0.0_r8 allocate(this%ED_biomass_patch (begp:endp)) ; this%ED_biomass_patch (:) = 0.0_r8 - allocate(this%storvegc_patch (begp:endp)) ; this%storvegc_patch (:) = nan - allocate(this%dispvegc_patch (begp:endp)) ; this%dispvegc_patch (:) = nan - allocate(this%leafc_patch (begp:endp)) ; this%leafc_patch (:) = nan - allocate(this%livestemc_patch (begp:endp)) ; this%livestemc_patch (:) = nan - allocate(this%deadstemc_patch (begp:endp)) ; this%deadstemc_patch (:) = nan - allocate(this%livestemn_patch (begp:endp)) ; this%livestemn_patch (:) = nan + ! allocate(this%storvegc_patch (begp:endp)) ; this%storvegc_patch (:) = nan + ! allocate(this%dispvegc_patch (begp:endp)) ; this%dispvegc_patch (:) = nan + ! allocate(this%leafc_patch (begp:endp)) ; this%leafc_patch (:) = nan + ! allocate(this%livestemc_patch (begp:endp)) ; this%livestemc_patch (:) = nan + ! allocate(this%deadstemc_patch (begp:endp)) ; this%deadstemc_patch (:) = nan + ! allocate(this%livestemn_patch (begp:endp)) ; this%livestemn_patch (:) = nan allocate(this%gpp_patch (begp:endp)) ; this%gpp_patch (:) = nan allocate(this%npp_patch (begp:endp)) ; this%npp_patch (:) = nan - begg = bounds%begg; endg = bounds%endg + allocate(this%ED_c_to_litr_lab_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_lab_c_col (:,:) = nan + allocate(this%ED_c_to_litr_cel_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_cel_c_col (:,:) = nan + allocate(this%ED_c_to_litr_lig_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_lig_c_col (:,:) = nan + + allocate(this%leaf_prof_col (begc:endc,1:nlevdecomp_full)) ; this%leaf_prof_col (:,:) = nan + allocate(this%froot_prof_col (begc:endc,1:numpft_ed,1:nlevdecomp_full)); this%froot_prof_col (:,:,:) = nan + allocate(this%croot_prof_col (begc:endc,1:nlevdecomp_full)) ; this%croot_prof_col (:,:) = nan + allocate(this%stem_prof_col (begc:endc,1:nlevdecomp_full)) ; this%stem_prof_col (:,:) = nan + + allocate(this%ed_to_bgc_this_edts_col (begc:endc)) ; this%ed_to_bgc_this_edts_col (:) = nan + allocate(this%ed_to_bgc_last_edts_col (begc:endc)) ; this%ed_to_bgc_last_edts_col (:) = nan + allocate(this%seed_rain_flux_col (begc:endc)) ; this%seed_rain_flux_col (:) = nan + + allocate(this%nep_col (begc:endc)) ; this%nep_col (:) = nan + allocate(this%nep_timeintegrated_col (begc:endc)) ; this%nep_timeintegrated_col (:) = nan + allocate(this%npp_timeintegrated_col (begc:endc)) ; this%npp_timeintegrated_col (:) = nan + allocate(this%hr_timeintegrated_col (begc:endc)) ; this%hr_timeintegrated_col (:) = nan + + allocate(this%nbp_col (begc:endc)) ; this%nbp_col (:) = nan + allocate(this%npp_hifreq_col (begc:endc)) ; this%npp_hifreq_col (:) = nan + allocate(this%fire_c_to_atm_col (begc:endc)) ; this%fire_c_to_atm_col (:) = nan + + allocate(this%totecosysc_col (begc:endc)) ; this%totecosysc_col (:) = nan + allocate(this%totecosysc_old_col (begc:endc)) ; this%totecosysc_old_col (:) = nan + allocate(this%totedc_col (begc:endc)) ; this%totedc_col (:) = nan + allocate(this%totedc_old_col (begc:endc)) ; this%totedc_old_col (:) = nan + allocate(this%totbgcc_col (begc:endc)) ; this%totbgcc_col (:) = nan + allocate(this%totbgcc_old_col (begc:endc)) ; this%totbgcc_old_col (:) = nan + allocate(this%biomass_stock_col (begc:endc)) ; this%biomass_stock_col (:) = nan + allocate(this%ed_litter_stock_col (begc:endc)) ; this%ed_litter_stock_col (:) = nan + allocate(this%cwd_stock_col (begc:endc)) ; this%cwd_stock_col (:) = nan + allocate(this%seed_stock_col (begc:endc)) ; this%seed_stock_col (:) = nan + + allocate(this%cbalance_error_ed_col (begc:endc)) ; this%cbalance_error_ed_col (:) = nan + allocate(this%cbalance_error_bgc_col (begc:endc)) ; this%cbalance_error_bgc_col (:) = nan + allocate(this%cbalance_error_total_col (begc:endc)) ; this%cbalance_error_total_col (:) = nan + allocate(this%ed_gpp_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_gpp_gd_scpf (:,:) = 0.0_r8 allocate(this%ed_npp_totl_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_totl_gd_scpf (:,:) = 0.0_r8 allocate(this%ed_npp_leaf_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_leaf_gd_scpf (:,:) = 0.0_r8 @@ -288,20 +379,20 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='Scaling factor between tree basal area and canopy area', & ptr_patch=this%canopy_spread_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTbiomass', units='kgC/m2', type2d='levgrnd', & + call hist_addfld2d (fname='PFTbiomass', units='gC/m2', type2d='levgrnd', & avgflag='A', long_name='total PFT level biomass', & ptr_patch=this%PFTbiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTleafbiomass', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTleafbiomass', units='gC/m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level leaf biomass', & ptr_patch=this%PFTleafbiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTstorebiomass', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTstorebiomass', units='gC/m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level stored biomass', & ptr_patch=this%PFTstorebiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTnindivs', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTnindivs', units='indiv / m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level number of individuals', & ptr_patch=this%PFTnindivs_patch, set_lake=0._r8, set_urb=0._r8) call hist_addfld1d (fname='FIRE_NESTEROV_INDEX', units='none', & @@ -348,59 +439,55 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='spitfire fuel surface/volume ', & ptr_patch=this%fire_fuel_sav_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='TFC_ROS', units='m', & - avgflag='A', long_name='spitfire fuel surface/volume ', & - ptr_patch=this%TFC_ROS_patch, set_lake=0._r8, set_urb=0._r8) - - call hist_addfld1d (fname='SUM_FUEL', units=' KgC m-2 y-1', & - avgflag='A', long_name='Litter flux in leaves', & + call hist_addfld1d (fname='SUM_FUEL', units='gC m-2', & + avgflag='A', long_name='total ground fuel related to ros (omits 1000hr fuels)', & ptr_patch=this%sum_fuel_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='LITTER_IN', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='LITTER_IN', units='gC m-2 s-1', & avgflag='A', long_name='Litter flux in leaves', & ptr_patch=this%litter_in_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='LITTER_OUT', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='LITTER_OUT', units='gC m-2 s-1', & avgflag='A', long_name='Litter flux out leaves', & ptr_patch=this%litter_out_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_BANK', units=' KgC m-2', & + call hist_addfld1d (fname='SEED_BANK', units='gC m-2', & avgflag='A', long_name='Total Seed Mass of all PFTs', & ptr_patch=this%seed_bank_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEEDS_IN', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEEDS_IN', units='gC m-2 s-1', & avgflag='A', long_name='Seed Production Rate', & ptr_patch=this%seeds_in_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_GERMINATION', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEED_GERMINATION', units='gC m-2 s-1', & avgflag='A', long_name='Seed mass converted into new cohorts', & ptr_patch=this%seed_germination_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_DECAY', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEED_DECAY', units='gC m-2 s-1', & avgflag='A', long_name='Seed mass decay', & ptr_patch=this%seed_decay_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bstore', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bstore', units='gC m-2', & avgflag='A', long_name='ED stored biomass', & ptr_patch=this%ED_bstore_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bdead', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bdead', units='gC m-2', & avgflag='A', long_name='ED dead biomass', & ptr_patch=this%ED_bdead_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_balive', units=' KgC m-2', & + call hist_addfld1d (fname='ED_balive', units='gC m-2', & avgflag='A', long_name='ED live biomass', & ptr_patch=this%ED_balive_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bleaf', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bleaf', units='gC m-2', & avgflag='A', long_name='ED leaf biomass', & ptr_patch=this%ED_bleaf_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_biomass', units=' KgC m-2', & + call hist_addfld1d (fname='ED_biomass', units='gC m-2', & avgflag='A', long_name='ED total biomass', & ptr_patch=this%ED_biomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='RB', units=' s m-1', & + call hist_addfld1d (fname='RB', units='s m-1', & avgflag='A', long_name='leaf boundary resistance', & ptr_patch=this%rb_patch, set_lake=0._r8, set_urb=0._r8) @@ -408,35 +495,35 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='potential evap', & ptr_patch=this%efpot_patch, set_lake=0._r8, set_urb=0._r8) - this%dispvegc_patch(begp:endp) = spval - call hist_addfld1d (fname='DISPVEGC', units='gC/m^2', & - avgflag='A', long_name='displayed veg carbon, excluding storage and cpool', & - ptr_patch=this%dispvegc_patch) + ! this%dispvegc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='DISPVEGC', units='gC/m^2', & + ! avgflag='A', long_name='displayed veg carbon, excluding storage and cpool', & + ! ptr_patch=this%dispvegc_patch) - this%storvegc_patch(begp:endp) = spval - call hist_addfld1d (fname='STORVEGC', units='gC/m^2', & - avgflag='A', long_name='stored vegetation carbon, excluding cpool', & - ptr_patch=this%storvegc_patch) + ! this%storvegc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='STORVEGC', units='gC/m^2', & + ! avgflag='A', long_name='stored vegetation carbon, excluding cpool', & + ! ptr_patch=this%storvegc_patch) - this%leafc_patch(begp:endp) = spval - call hist_addfld1d (fname='LEAFC', units='gC/m^2', & - avgflag='A', long_name='leaf C', & - ptr_patch=this%leafc_patch) + ! this%leafc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LEAFC', units='gC/m^2', & + ! avgflag='A', long_name='leaf C', & + ! ptr_patch=this%leafc_patch) - this%livestemc_patch(begp:endp) = spval - call hist_addfld1d (fname='LIVESTEMC', units='gC/m^2', & - avgflag='A', long_name='live stem C', & - ptr_patch=this%livestemc_patch) + ! this%livestemc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LIVESTEMC', units='gC/m^2', & + ! avgflag='A', long_name='live stem C', & + ! ptr_patch=this%livestemc_patch) - this%deadstemc_patch(begp:endp) = spval - call hist_addfld1d (fname='DEADSTEMC', units='gC/m^2', & - avgflag='A', long_name='dead stem C', & - ptr_patch=this%deadstemc_patch) + ! this%deadstemc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='DEADSTEMC', units='gC/m^2', & + ! avgflag='A', long_name='dead stem C', & + ! ptr_patch=this%deadstemc_patch) - this%livestemn_patch(begp:endp) = spval - call hist_addfld1d (fname='LIVESTEMN', units='gN/m^2', & - avgflag='A', long_name='live stem N', & - ptr_patch=this%livestemn_patch) + ! this%livestemn_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LIVESTEMN', units='gN/m^2', & + ! avgflag='A', long_name='live stem N', & + ! ptr_patch=this%livestemn_patch) this%gpp_patch(begp:endp) = spval call hist_addfld1d (fname='GPP', units='gC/m^2/s', & @@ -448,6 +535,97 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='net primary production', & ptr_patch=this%npp_patch) + this%nep_col(begc:endc) = spval + call hist_addfld1d (fname='NEP', units='gC/m^2/s', & + avgflag='A', long_name='net ecosystem production', & + ptr_col=this%nep_col) + + this%fire_c_to_atm_col(begc:endc) = spval + call hist_addfld1d (fname='Fire_Closs', units='gC/m^2/s', & + avgflag='A', long_name='ED/SPitfire Carbon loss to atmosphere', & + ptr_col=this%fire_c_to_atm_col) + + this%nbp_col(begc:endc) = spval + call hist_addfld1d (fname='NBP', units='gC/m^2/s', & + avgflag='A', long_name='net biosphere production', & + ptr_col=this%nbp_col) + + this%npp_hifreq_col(begc:endc) = spval + call hist_addfld1d (fname='NPP_hifreq', units='gC/m^2/s', & + avgflag='A', long_name='net primary production at high frequency', & + ptr_col=this%npp_hifreq_col,default='inactive') + + this%totecosysc_col(begc:endc) = spval + call hist_addfld1d (fname='TOTECOSYSC', units='gC/m^2', & + avgflag='A', long_name='total ecosystem carbon', & + ptr_col=this%totecosysc_col) + + this%cbalance_error_ed_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_ED', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error on ED side', & + ptr_col=this%cbalance_error_ed_col) + + this%cbalance_error_bgc_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_BGC', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error on BGC side', & + ptr_col=this%cbalance_error_bgc_col) + + this%cbalance_error_total_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_TOTAL', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error total', & + ptr_col=this%cbalance_error_total_col) + + this%biomass_stock_col(begc:endc) = spval + call hist_addfld1d (fname='BIOMASS_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total ED biomass carbon at the column level', & + ptr_col=this%biomass_stock_col) + + this%ed_litter_stock_col(begc:endc) = spval + call hist_addfld1d (fname='ED_LITTER_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total ED litter carbon at the column level', & + ptr_col=this%ed_litter_stock_col) + + this%cwd_stock_col(begc:endc) = spval + call hist_addfld1d (fname='CWD_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total CWD carbon at the column level', & + ptr_col=this%cwd_stock_col) + + this%seed_stock_col(begc:endc) = spval + call hist_addfld1d (fname='SEED_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total seed carbon at the column level', & + ptr_col=this%seed_stock_col) + + !!! carbon fluxes into soil grid (dimensioned depth x column) + this%ED_c_to_litr_lab_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_lab_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_lab_c', & + ptr_col=this%ED_c_to_litr_lab_c_col) + + this%ED_c_to_litr_cel_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_cel_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_cel_c', & + ptr_col=this%ED_c_to_litr_cel_c_col) + + this%ED_c_to_litr_lig_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_lig_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_lig_c', & + ptr_col=this%ED_c_to_litr_lig_c_col) + + this%leaf_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='leaf_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='leaf_prof', & + ptr_col=this%leaf_prof_col,default='inactive') + + this%croot_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='croot_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='croot_prof', & + ptr_col=this%croot_prof_col,default='inactive') + + this%stem_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='stem_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='stem_prof', & + ptr_col=this%stem_prof_col,default='inactive') + ! Carbon Flux (grid dimension x scpf) ! ============================================================== @@ -523,26 +701,25 @@ subroutine InitHistory(this, bounds) end subroutine InitHistory !----------------------------------------------------------------------- - subroutine InitCold(this, bounds) - ! - ! !DESCRIPTION: - ! Initialize relevant time varying variables - ! - ! !ARGUMENTS: - class (ed_clm_type) :: this - type(bounds_type), intent(in) :: bounds - ! - ! !LOCAL VARIABLES: - integer :: p - !----------------------------------------------------------------------- - - do p = bounds%begp,bounds%endp - this%dispvegc_patch(p) = 0._r8 - this%storvegc_patch(p) = 0._r8 - end do - - end subroutine InitCold - + ! subroutine InitCold(this, bounds) + ! ! + ! ! !DESCRIPTION: + ! ! Initialize relevant time varying variables + ! ! + ! ! !ARGUMENTS: + ! class (ed_clm_type) :: this + ! type(bounds_type), intent(in) :: bounds + ! ! + ! ! !LOCAL VARIABLES: + ! integer :: p + ! !----------------------------------------------------------------------- + + ! ! do p = bounds%begp,bounds%endp + ! ! this%dispvegc_patch(p) = 0._r8 + ! ! this%storvegc_patch(p) = 0._r8 + ! ! end do + + ! end subroutine InitCold !----------------------------------------------------------------------- subroutine Restart ( this, bounds, ncid, flag ) ! @@ -552,6 +729,8 @@ subroutine Restart ( this, bounds, ncid, flag ) ! !USES: use restUtilMod use ncdio_pio + ! use EDtypesMod , only : numpft_ed + ! ! !ARGUMENTS: class (ed_clm_type) :: this @@ -561,23 +740,173 @@ subroutine Restart ( this, bounds, ncid, flag ) ! ! !LOCAL VARIABLES: logical :: readvar + real(r8), pointer :: ptr2d(:,:) ! temp. pointers for slicing larger arrays + real(r8), pointer :: ptr1d(:) ! temp. pointers for slicing larger arrays + ! character(LEN=3) :: istr1 + ! integer :: k !------------------------------------------------------------------------ - call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%leafc_patch) - - call restartvar(ncid=ncid, flag=flag, varname='livestemc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%livestemc_patch) - - call restartvar(ncid=ncid, flag=flag, varname='deadstemc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%deadstemc_patch) + ! call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%leafc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='livestemc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%livestemc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='deadstemc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%deadstemc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='livestemn', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%livestemn_patch) + + ptr1d => this%nep_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='nep_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%npp_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='npp_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%hr_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='hr_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totecosysc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totecosysc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_ed_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_ed_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_bgc_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_bgc_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_total_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_total_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totedc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totedc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totbgcc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totbgcc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ed_to_bgc_this_edts_col(:) + call restartvar(ncid=ncid, flag=flag, varname='ed_to_bgc_this_edts_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ed_to_bgc_last_edts_col(:) + call restartvar(ncid=ncid, flag=flag, varname='ed_to_bgc_last_edts_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%seed_rain_flux_col(:) + call restartvar(ncid=ncid, flag=flag, varname='seed_rain_flux_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + if (use_vertsoilc) then + ptr2d => this%ED_c_to_litr_lab_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lab_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ptr2d => this%ED_c_to_litr_cel_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_cel_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ptr2d => this%ED_c_to_litr_lig_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lig_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%leaf_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='leaf_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%croot_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='croot_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%stem_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='stem_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! do k = 1, numpft_ed + ! write(istr1,"(I3.3)") k + ! ptr2d => this%froot_prof_col(:,k,:) + ! call restartvar(ncid=ncid, flag=flag, varname='froot_prof_col_PFT'//istr1, xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + ! end do + else + ptr1d => this%ED_c_to_litr_lab_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lab_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ED_c_to_litr_cel_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_cel_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ED_c_to_litr_lig_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lig_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%leaf_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='leaf_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%croot_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='croot_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%stem_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='stem_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! do k = 1, numpft_ed + ! write(istr1,"(I3.3)") k + ! ptr1d => this%froot_prof_col(:,k,1) + ! call restartvar(ncid=ncid, flag=flag, varname='froot_prof_col_PFT'//istr1, xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + ! end do + end if - call restartvar(ncid=ncid, flag=flag, varname='livestemn', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%livestemn_patch) end subroutine Restart @@ -870,7 +1199,10 @@ subroutine ed_clm_link( this, bounds, ed_allsites_inst, ed_phenology_inst, & end do !grid loop - call this%ed_update_history_variables( bounds, ed_allsites_inst(begg:endg), & + call this%flux_into_litter_pools(bounds, ed_allsites_inst(begg:endg), firstsoilpatch, & + canopystate_inst) + + call this%ed_update_history_variables(bounds, ed_allsites_inst(begg:endg), & firstsoilpatch, ed_Phenology_inst, canopystate_inst) end associate @@ -1104,16 +1436,16 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & write(iulog,*) 'EDCLMLinkMod II ',p,ED_bstore(p) endif - ED_bleaf(p) = ED_bleaf(p) + n_density * currentCohort%bl - ED_bstore(p) = ED_bstore(p) + n_density * currentCohort%bstore - ED_biomass(p) = ED_biomass(p) + n_density * currentCohort%b - ED_bdead(p) = ED_bdead(p) + n_density * currentCohort%bdead - ED_balive(p) = ED_balive(p) + n_density * currentCohort%balive - npp(p) = npp(p) + n_density * currentCohort%npp - gpp(p) = gpp(p) + n_density * currentCohort%gpp - PFTbiomass(p,ft) = PFTbiomass(p,ft) + n_density * currentCohort%b - PFTleafbiomass(p,ft) = PFTleafbiomass(p,ft) + n_density * currentCohort%bl - PFTstorebiomass(p,ft) = PFTstorebiomass(p,ft) + n_density * currentCohort%bstore + ED_bleaf(p) = ED_bleaf(p) + n_density * currentCohort%bl * 1.e3_r8 + ED_bstore(p) = ED_bstore(p) + n_density * currentCohort%bstore * 1.e3_r8 + ED_biomass(p) = ED_biomass(p) + n_density * currentCohort%b * 1.e3_r8 + ED_bdead(p) = ED_bdead(p) + n_density * currentCohort%bdead * 1.e3_r8 + ED_balive(p) = ED_balive(p) + n_density * currentCohort%balive * 1.e3_r8 + npp(p) = npp(p) + n_density * currentCohort%npp * 1.e3_r8 / (365. * SHR_CONST_CDAY) + gpp(p) = gpp(p) + n_density * currentCohort%gpp * 1.e3_r8 / (365. * SHR_CONST_CDAY) + PFTbiomass(p,ft) = PFTbiomass(p,ft) + n_density * currentCohort%b * 1.e3_r8 + PFTleafbiomass(p,ft) = PFTleafbiomass(p,ft) + n_density * currentCohort%bl * 1.e3_r8 + PFTstorebiomass(p,ft) = PFTstorebiomass(p,ft) + n_density * currentCohort%bstore * 1.e3_r8 PFTnindivs(p,ft) = PFTnindivs(p,ft) + currentCohort%n dbh = currentCohort%dbh !-0.5*(1./365.25)*currentCohort%ddbhdt @@ -1188,13 +1520,13 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & fire_fuel_eff_moist(p) = currentPatch%fuel_eff_moist fire_fuel_sav(p) = currentPatch%fuel_sav fire_fuel_mef(p) = currentPatch%fuel_mef - sum_fuel(p) = currentPatch%sum_fuel - litter_in(p) = sum(currentPatch%CWD_AG_in) +sum(currentPatch%leaf_litter_in) - litter_out(p) = sum(currentPatch%CWD_AG_out)+sum(currentPatch%leaf_litter_out) - seed_bank(p) = sum(currentPatch%seed_bank) - seeds_in(p) = sum(currentPatch%seeds_in) - seed_decay(p) = sum(currentPatch%seed_decay) - seed_germination(p) = sum(currentPatch%seed_germination) + sum_fuel(p) = currentPatch%sum_fuel * 1.e3_r8 + litter_in(p) = (sum(currentPatch%CWD_AG_in) +sum(currentPatch%leaf_litter_in)) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + litter_out(p) = (sum(currentPatch%CWD_AG_out)+sum(currentPatch%leaf_litter_out)) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_bank(p) = sum(currentPatch%seed_bank) * 1.e3_r8 + seeds_in(p) = sum(currentPatch%seeds_in) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_decay(p) = sum(currentPatch%seed_decay) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_germination(p) = sum(currentPatch%seed_germination) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY canopy_spread(p) = currentPatch%spread(1) area_plant(p) = currentPatch%total_canopy_area /currentPatch%area area_trees(p) = currentPatch%total_tree_area /currentPatch%area @@ -1701,4 +2033,674 @@ subroutine ed_clm_leaf_area_profile( this, currentSite, waterstate_inst, canopys end subroutine ed_clm_leaf_area_profile + + subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch, canopystate_inst) + ! Created by Charlie Koven and Rosie Fisher, 2014-2015 + ! take the flux out of the fragmenting litter pools and port into the decomposing litter pools. + ! in this implementation, decomposing pools are assumed to be humus and non-flammable, whereas fragmenting pools + ! are assumed to be physically fragmenting but not respiring. This is a simplification, but allows us to + ! a) reconcile the need to track both chemical fractions (lignin, cellulose, labile) and size fractions (trunk, branch, etc.) + ! b) to impose a realistic delay on the surge of nutrients into the litter pools when large CWD is added to the system via mortality + + ! because of the different subgrid structure, this subroutine includes the functionality that in the big-leaf BGC model, is calculated in SoilBiogeochemVerticalProfileMod + + ! The ED code is resolved at a daily timestep, but all of the CN-BGC fluxes are passed in as derivatives per second, + ! and then accumulated in the CNStateUpdate routines. One way of doing this is to pass back the CN fluxes per second, + ! and keep them constant for the whole day (making sure they are not overwritten. + ! This means that the carbon gets passed back and forth between the photosynthesis code (fast timestepping) to the ED code (slow timestepping), back to the BGC code (fast timestepping). + ! This means that the state update for the litter pools and for the CWD pools occurs at different timescales. + + use SFParamsMod, only: SF_val_max_decomp + use clm_varpar, only : mxpft,nlevdecomp, nlevdecomp_full + use EDTypesMod, only : AREA, numpft_ed + use SoilBiogeochemVerticalProfileMod, only: exponential_rooting_profile, pftspecific_rootingprofile, rootprof_exp, surfprof_exp + use pftconMod, only : pftcon + + use clm_varcon, only : zisoi, dzsoi_decomp, zsoi + use ColumnType , only : col + use abortutils , only : endrun + use shr_log_mod , only : errMsg => shr_log_errMsg + use EDParamsMod, only : ED_val_ag_biomass + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + integer :: firstsoilpatch(bounds%begg:bounds%endg) ! the first patch in this gridcell that is soil and thus bare... + type(canopystate_type) , intent(inout) :: canopystate_inst + ! + ! !LOCAL VARIABLES: + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort + type(ed_site_type), pointer :: cs + integer c,p,cc,j,g + real(r8) time_convert ! from year to seconds + real(r8) mass_convert ! ED uses kg, CLM uses g + integer :: begp,endp + integer :: begc,endc !bounds + !------------------------------------------------------------------------ + real(r8) :: cinput_rootfr(bounds%begc:bounds%endc, 1:numpft_ed, 1:nlevdecomp_full) ! column by pft root fraction used for calculating inputs + real(r8) :: croot_prof_perpatch(1:nlevdecomp_full) + real(r8) :: surface_prof(1:nlevdecomp_full) + integer :: ft, lev + real(r8) :: rootfr_tot(1:numpft_ed), biomass_bg_ft(1:numpft_ed) + real(r8) :: surface_prof_tot, leaf_prof_sum, stem_prof_sum, froot_prof_sum, biomass_bg_tot + real(r8) :: delta + + begp = bounds%begp; endp = bounds%endp + begc = bounds%begc; endc = bounds%endc + delta = 0.001_r8 + !no of seconds in a year. + time_convert = 365.0_r8*SHR_CONST_CDAY + + ! number of grams in a kilogram + mass_convert = 1000._r8 + + associate( & + ED_c_to_litr_lab_c => this%ED_c_to_litr_lab_c_col , & ! Output: total labile litter coming from ED. gC/m3/s + ED_c_to_litr_cel_c => this%ED_c_to_litr_cel_c_col , & ! Output: total cellulose litter coming from ED. gC/m3/s + ED_c_to_litr_lig_c => this%ED_c_to_litr_lig_c_col , & ! Output: total lignin litter coming from ED. gC/m3/s + leaf_prof => this%leaf_prof_col , & ! Output: (1/m) profile of leaves + froot_prof => this%froot_prof_col , & ! Output: (1/m) profile of fine roots + croot_prof => this%croot_prof_col , & ! Output: (1/m) profile of coarse roots + stem_prof => this%stem_prof_col , & ! Output: (1/m) profile of leaves + altmax_lastyear_indx => canopystate_inst%altmax_lastyear_indx_col & ! Input: [integer (:) ]frost table depth (m) + ) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! first calculate vertical profiles + ! define two types of profiles: + ! (1) a surface profile, for leaves and stem inputs, which is the same for each pft but differs from one column to the next to avoid inputting any C into permafrost + ! (2) a fine root profile, which is indexed by both column and pft, differs for each pft and also from one column to the next to avoid inputting any C into permafrost + ! (3) a coarse root profile, which is the root-biomass=weighted average of the fine root profiles + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (use_vertsoilc) then + + ! define a single shallow surface profile for surface additions (leaves, stems, and N deposition) + surface_prof(:) = 0._r8 + do j = 1, nlevdecomp + surface_prof(j) = exp(-surfprof_exp * zsoi(j)) / dzsoi_decomp(j) + end do + + ! initialize profiles to zero + leaf_prof(begc:endc, :) = 0._r8 + froot_prof(begc:endc, 1:numpft_ed, :) = 0._r8 + croot_prof(begc:endc, :) = 0._r8 + stem_prof(begc:endc, :) = 0._r8 + + cinput_rootfr(begc:endc, 1:numpft_ed, :) = 0._r8 + + do c = bounds%begc,bounds%endc + + ! calculate pft-specific rooting profiles in the absence of permafrost limitations + if ( exponential_rooting_profile ) then + if ( .not. pftspecific_rootingprofile ) then + ! define rooting profile from exponential parameters + do ft = 1, numpft_ed + do j = 1, nlevdecomp + cinput_rootfr(c,ft,j) = exp(-rootprof_exp * zsoi(j)) / dzsoi_decomp(j) + end do + end do + else + ! use beta distribution parameter from Jackson et al., 1996 + do ft = 1, numpft_ed + do j = 1, nlevdecomp + cinput_rootfr(c,ft,j) = ( pftcon%rootprof_beta(ft) ** (zisoi(j-1)*100._r8) - & + pftcon%rootprof_beta(ft) ** (zisoi(j)*100._r8) ) & + / dzsoi_decomp(j) + end do + end do + endif + else + do ft = 1,numpft_ed + do j = 1, nlevdecomp + ! use standard CLM root fraction profiles; + cinput_rootfr(c,ft,j) = ( .5_r8*( & + exp(-pftcon%roota_par(ft) * col%zi(c,lev-1)) & + + exp(-pftcon%rootb_par(ft) * col%zi(c,lev-1)) & + - exp(-pftcon%roota_par(ft) * col%zi(c,lev)) & + - exp(-pftcon%rootb_par(ft) * col%zi(c,lev)))) / dzsoi_decomp(j) + end do + end do + endif + ! + ! + ! now add permafrost constraint: integrate rootfr over active layer of soil column, + ! truncate below permafrost table where present, and rescale so that integral = 1 + do ft = 1,numpft_ed + rootfr_tot(ft) = 0._r8 + end do + surface_prof_tot = 0._r8 + ! + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + surface_prof_tot = surface_prof_tot + surface_prof(j) * dzsoi_decomp(j) + end do + do ft = 1,numpft_ed + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(c,ft,j) * dzsoi_decomp(j) + end do + end do + ! + ! rescale the fine root profile + do ft = 1,numpft_ed + if ( (altmax_lastyear_indx(c) > 0) .and. (rootfr_tot(ft) > 0._r8) ) then + ! where there is not permafrost extending to the surface, integrate the profiles over the active layer + ! this is equivalent to integrating over all soil layers outside of permafrost regions + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + froot_prof(c,ft,j) = cinput_rootfr(c,ft,j) / rootfr_tot(ft) + end do + else + ! if fully frozen, or no roots, put everything in the top layer + froot_prof(c,ft,1) = 1._r8/dzsoi_decomp(1) + endif + end do + ! + ! rescale the shallow profiles + if ( (altmax_lastyear_indx(c) > 0) .and. (surface_prof_tot > 0._r8) ) then + ! where there is not permafrost extending to the surface, integrate the profiles over the active layer + ! this is equivalent to integrating over all soil layers outside of permafrost regions + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + ! set all surface processes to shallower profile + leaf_prof(c,j) = surface_prof(j)/ surface_prof_tot + stem_prof(c,j) = surface_prof(j)/ surface_prof_tot + end do + else + ! if fully frozen, or no roots, put everything in the top layer + leaf_prof(c,1) = 1._r8/dzsoi_decomp(1) + stem_prof(c,1) = 1._r8/dzsoi_decomp(1) + do j = 2, nlevdecomp + leaf_prof(c,j) = 0._r8 + stem_prof(c,j) = 0._r8 + end do + endif + end do + + else + + ! for one layer decomposition model, set profiles to unity + leaf_prof(bounds%begc:bounds%endc, :) = 1._r8 + froot_prof(bounds%begc:bounds%endc, 1:numpft_ed, :) = 1._r8 + stem_prof(bounds%begc:bounds%endc, :) = 1._r8 + + end if + + ! sanity check to ensure they integrate to 1 + do c = bounds%begc,bounds%endc + ! check the leaf and stem profiles + leaf_prof_sum = 0._r8 + stem_prof_sum = 0._r8 + do j = 1, nlevdecomp + leaf_prof_sum = leaf_prof_sum + leaf_prof(c,j) * dzsoi_decomp(j) + stem_prof_sum = stem_prof_sum + stem_prof(c,j) * dzsoi_decomp(j) + end do + if ( ( abs(stem_prof_sum - 1._r8) > delta ) .or. ( abs(leaf_prof_sum - 1._r8) > delta ) ) then + write(iulog, *) 'profile sums: ', leaf_prof_sum, stem_prof_sum + write(iulog, *) 'surface_prof: ', surface_prof + write(iulog, *) 'surface_prof_tot: ', surface_prof_tot + write(iulog, *) 'leaf_prof: ', leaf_prof(c,:) + write(iulog, *) 'stem_prof: ', stem_prof(c,:) + write(iulog, *) 'altmax_lastyear_indx: ', altmax_lastyear_indx(c) + write(iulog, *) 'dzsoi_decomp: ', dzsoi_decomp + call endrun(msg=' ERROR: sum-1 > delta'//errMsg(__FILE__, __LINE__)) + endif + ! now check each fine root profile + do ft = 1,numpft_ed + froot_prof_sum = 0._r8 + do j = 1, nlevdecomp + froot_prof_sum = froot_prof_sum + froot_prof(c,ft,j) * dzsoi_decomp(j) + end do + if ( ( abs(froot_prof_sum - 1._r8) > delta ) ) then + write(iulog, *) 'profile sums: ', froot_prof_sum + call endrun(msg=' ERROR: sum-1 > delta'//errMsg(__FILE__, __LINE__)) + endif + end do + end do + + ! zero the column-level C input variables + do c = bounds%begc,bounds%endc + do j = 1, nlevdecomp + ED_c_to_litr_lab_c(c,j) = 0._r8 + ED_c_to_litr_cel_c(c,j) = 0._r8 + ED_c_to_litr_lig_c(c,j) = 0._r8 + croot_prof(c,j) = 0._r8 + end do + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! now disaggregate the inputs vertically, using the vertical profiles + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + currentPatch => ed_allsites_inst(g)%oldest_patch + + do while(associated(currentPatch)) + + cs => currentpatch%siteptr + cc = cs%clmcolumn + + ! the CWD pools lose information about which PFT they came from; for the stems this doesn't matter as they all have the same profile, + ! however for the coarse roots they may have different profiles. to approximately recover this information, loop over all cohorts in patch + ! to calculate the total root biomass in that patch of each pft, and then rescale the croot_prof as the weighted average of the froot_prof + biomass_bg_ft(1:numpft_ed) = 0._r8 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + biomass_bg_ft(currentCohort%pft) = biomass_bg_ft(currentCohort%pft) + currentCohort%b * (currentCohort%n / currentPatch%area) * (1.0_r8-ED_val_ag_biomass) + currentCohort => currentCohort%shorter + enddo !currentCohort + ! + biomass_bg_tot = 0._r8 + do ft = 1,numpft_ed + biomass_bg_tot = biomass_bg_tot + biomass_bg_ft(ft) + end do + ! + do j = 1, nlevdecomp + ! zero this for each patch + croot_prof_perpatch(j) = 0._r8 + end do + ! + if ( biomass_bg_tot .gt. 0._r8) then + do ft = 1,numpft_ed + do j = 1, nlevdecomp + croot_prof_perpatch(j) = croot_prof_perpatch(j) + froot_prof(cc,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot + end do + end do + else ! no biomass + croot_prof_perpatch(1) = 1./dzsoi_decomp(1) + end if + ! + ! add croot_prof as weighted average (weighted by patch area) of croot_prof_perpatch + do j = 1, nlevdecomp + croot_prof(cc, j) = croot_prof(cc, j) + croot_prof_perpatch(j) * currentPatch%area / AREA + end do + ! + ! now disaggregate, vertically and by decomposition substrate type, the actual fluxes from CWD and litter pools + ! + ! do c = 1, ncwd + ! write(iulog,*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), cwd_fcel_ed, currentpatch%area/AREA + ! write(iulog,*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), cwd_fcel_ed, currentpatch%area/AREA + ! end do + ! do ft = 1,numpft_ed + ! write(iulog,*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), cwd_fcel_ed, currentpatch%area/AREA + ! write(iulog,*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), cwd_fcel_ed, currentpatch%area/AREA + ! end do + ! ! + ! CWD pools fragmenting into decomposing litter pools. + do c = 1, ncwd + do j = 1, nlevdecomp + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%CWD_AG_out(c) * cwd_fcel_ed * currentpatch%area/AREA * stem_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%CWD_AG_out(c) * cwd_flig_ed * currentpatch%area/AREA * stem_prof(cc,j) + ! + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%CWD_BG_out(c) * cwd_fcel_ed * currentpatch%area/AREA * croot_prof_perpatch(j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%CWD_BG_out(c) * cwd_flig_ed * currentpatch%area/AREA * croot_prof_perpatch(j) + end do + end do + + ! leaf and fine root pools. + do ft = 1,numpft_ed + do j = 1, nlevdecomp + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ! + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_flab(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_fcel(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_flig(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ! + !! and seed_decay too. for now, use the same lability fractions as for leaf litter + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ! + enddo + end do + + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + do cc = bounds%begc,bounds%endc + do j = 1, nlevdecomp + ! time unit conversion + ED_c_to_litr_lab_c(cc,j)=ED_c_to_litr_lab_c(cc,j) * mass_convert / time_convert + ED_c_to_litr_cel_c(cc,j)=ED_c_to_litr_cel_c(cc,j) * mass_convert / time_convert + ED_c_to_litr_lig_c(cc,j)=ED_c_to_litr_lig_c(cc,j) * mass_convert / time_convert + + end do + end do + + ! write(iulog,*)'cdk ED_c_to_litr_lab_c: ', ED_c_to_litr_lab_c + ! write(iulog,*)'cdk ED_c_to_litr_cel_c: ', ED_c_to_litr_cel_c + ! write(iulog,*)'cdk ED_c_to_litr_lig_c: ', ED_c_to_litr_lig_c + ! write(iulog,*)'cdk nlevdecomp_full, bounds%begc, bounds%endc: ', nlevdecomp_full, bounds%begc, bounds%endc + ! write(iulog,*)'cdk leaf_prof: ', leaf_prof + ! write(iulog,*)'cdk stem_prof: ', stem_prof + ! write(iulog,*)'cdk froot_prof: ', froot_prof + ! write(iulog,*)'cdk croot_prof_perpatch: ', croot_prof_perpatch + ! write(iulog,*)'cdk croot_prof: ', croot_prof + + end associate + end subroutine flux_into_litter_pools + + !------------------------------------------------------------------------ + + subroutine Summary(this, bounds, num_soilc, filter_soilc, & + ed_allsites_inst, soilbiogeochem_carbonflux_inst, & + soilbiogeochem_carbonstate_inst) + + ! Summarize the combined production and decomposition fluxes into net fluxes + ! Written by Charlie Koven, Feb 2016 + ! + ! !USES: + use ColumnType , only : col + use LandunitType , only : lun + use landunit_varcon , only : istsoil + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_soilc ! number of soil columns in filter + integer , intent(in) :: filter_soilc(:) ! filter for soil columns + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst + ! + ! !LOCAL VARIABLES: + real(r8) :: npp_hifreq_col(bounds%begc:bounds%endc) ! column-level, high frequency NPP + real(r8) :: dt ! radiation time step (seconds) + integer :: c, g, cc, fc, l + type(ed_site_type), pointer :: cs + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort + integer :: firstsoilpatch(bounds%begg:bounds%endg) ! the first patch in this gridcell that is soil and thus bare... + + associate(& + hr => soilbiogeochem_carbonflux_inst%hr_col, & ! (gC/m2/s) total heterotrophic respiration + totsomc => soilbiogeochem_carbonstate_inst%totsomc_col, & ! (gC/m2) total soil organic matter carbon + totlitc => soilbiogeochem_carbonstate_inst%totlitc_col, & ! (gC/m2) total litter carbon in BGC pools + npp_hifreq => this%npp_hifreq_col, & + nep => this%nep_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & + nbp => this%nbp_col, & + totecosysc => this%totecosysc_col, & + totedc => this%totedc_col, & + totbgcc => this%totbgcc_col, & + biomass_stock => this%biomass_stock_col, & ! total biomass in gC / m2 + ed_litter_stock => this%ed_litter_stock_col, & ! ED litter in gC / m2 + cwd_stock => this%cwd_stock_col, & ! total CWD in gC / m2 + seed_stock => this%seed_stock_col, & ! total seed mass in gC / m2 + ed_to_bgc_this_edts => this%ed_to_bgc_this_edts_col, & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col & + ) + + ! set time steps + dt = real( get_step_size(), r8 ) + + ! zero variables first + do c = bounds%begc,bounds%endc + ! summary flux variables + npp_hifreq(c) = 0._r8 + fire_c_to_atm(c) = 0._r8 + + ! summary stock variables + ed_litter_stock(c) = 0._r8 + cwd_stock(c) = 0._r8 + seed_stock(c) = 0._r8 + biomass_stock(c) = 0._r8 + end do + + ! retrieve the first soil patch associated with each gridcell. + ! make sure we only get the first patch value for places which have soil. + firstsoilpatch(bounds%begg:bounds%endg) = -999 + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + l = col%landunit(c) + if (lun%itype(l) == istsoil .and. col%itype(c) == istsoil) then + firstsoilpatch(g) = col%patchi(c) + endif + enddo + + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + + ! map ed site-level fire fluxes to clm column fluxes + fire_c_to_atm(cc) = ed_allsites_inst(g)%total_burn_flux_to_atm / ( AREA * SHR_CONST_CDAY * 1.e3_r8) + + currentPatch => ed_allsites_inst(g)%oldest_patch + do while(associated(currentPatch)) + + ! map litter, CWD, and seed pools to column level + cwd_stock(cc) = cwd_stock(cc) + (currentPatch%area / AREA) * (sum(currentPatch%cwd_ag)+ & + sum(currentPatch%cwd_bg)) * 1.e3_r8 + ed_litter_stock(cc) = ed_litter_stock(cc) + (currentPatch%area / AREA) * & + (sum(currentPatch%leaf_litter)+sum(currentPatch%root_litter)) * 1.e3_r8 + seed_stock(cc) = seed_stock(cc) + (currentPatch%area / AREA) * sum(currentPatch%seed_bank) * 1.e3_r8 + + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + + ! map ed cohort-level npp fluxes to clm column fluxes + npp_hifreq(cc) = npp_hifreq(cc) + currentCohort%npp_clm * 1e3 * currentCohort%n / ( AREA * dt) + + ! map biomass pools to column level + biomass_stock(cc) = biomass_stock(cc) + (currentCohort%bdead + currentCohort%balive + & + currentCohort%bstore) * currentCohort%n * 1.e3_r8 / AREA + + currentCohort => currentCohort%shorter + enddo !currentCohort + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + ! calculate NEP and NBP fluxes. + do fc = 1,num_soilc + c = filter_soilc(fc) + nep(c) = npp_hifreq(c) - hr(c) + nbp(c) = npp_hifreq(c) - ( hr(c) + fire_c_to_atm(c) ) + end do + + ! calculate total stocks + do fc = 1,num_soilc + c = filter_soilc(fc) + + totedc(c) = ed_litter_stock(c) + cwd_stock(c) + seed_stock(c) + biomass_stock(c) ! ED stocks + totbgcc(c) = totsomc(c) + totlitc(c) ! BGC stocks + totecosysc(c) = totedc(c) + totbgcc(c) + + end do + + ! in ED timesteps, because of offset between when ED and BGC reconcile the gain and loss of litterfall carbon, + ! (i.e. ED reconciles it instantly, while BGC reconciles it incrementally over the subsequent day) + ! calculate the total ED -> BGC flux and keep track of the last day's info for balance checking purposes + if ( is_beg_curr_day() ) then + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + ed_to_bgc_last_edts(cc) = ed_to_bgc_this_edts(cc) + endif + end do + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + ed_to_bgc_this_edts(cc) = 0._r8 + seed_rain_flux(cc) = 0._r8 + endif + end do + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + currentPatch => ed_allsites_inst(g)%oldest_patch + do while(associated(currentPatch)) + ! + ed_to_bgc_this_edts(cc) = ed_to_bgc_this_edts(cc) + (sum(currentPatch%CWD_AG_out) + sum(currentPatch%CWD_BG_out) + & + + sum(currentPatch%seed_decay) + sum(currentPatch%leaf_litter_out) + sum(currentPatch%root_litter_out)) & + * ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + ! + seed_rain_flux(cc) = seed_rain_flux(cc) + sum(currentPatch%seed_rain_flux) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + ! + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + endif + + end associate + + end subroutine Summary + + + subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soilbiogeochem_carbonflux_inst) + + ! Integrate in time the fluxes into and out of the ecosystem, and compare these on a daily timestep + ! to the chagne in carbon stocks of the ecosystem + ! + ! Written by Charlie Koven, Feb 2016 + ! + ! !USES: + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_soilc ! number of soil columns in filter + integer , intent(in) :: filter_soilc(:) ! filter for soil columns + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + ! + ! !LOCAL VARIABLES: + real(r8) :: dtime ! land model time step (sec) + integer :: nstep ! model timestep + real(r8) :: nbp_integrated(bounds%begc:bounds%endc) ! total net biome production integrated + real(r8) :: error_total(bounds%begc:bounds%endc) + real(r8) :: error_ed(bounds%begc:bounds%endc) + real(r8) :: error_bgc(bounds%begc:bounds%endc) + real(r8) :: error_tolerance = 1.e-6_r8 + real(r8) :: max_error_ed + real(r8) :: max_error_bgc + real(r8) :: max_error_total + integer :: fc,c + + associate(& + nep => this%nep_col, & + nep_timeintegrated => this%nep_timeintegrated_col, & + hr => soilbiogeochem_carbonflux_inst%hr_col, & + hr_timeintegrated => this%hr_timeintegrated_col, & + npp_hifreq => this%npp_hifreq_col, & + npp_timeintegrated => this%npp_timeintegrated_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & + totecosysc_old => this%totecosysc_old_col, & + totecosysc => this%totecosysc_col, & + totedc_old => this%totedc_old_col, & + totedc => this%totedc_col, & + totbgcc_old => this%totbgcc_old_col, & + totbgcc => this%totbgcc_col, & + ed_to_bgc_this_edts => this%ed_to_bgc_this_edts_col, & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col, & + cbalance_error_ed => this%cbalance_error_ed_col, & + cbalance_error_bgc => this%cbalance_error_bgc_col, & + cbalance_error_total=> this%cbalance_error_total_col & + ) + + dtime = get_step_size() + nstep = get_nstep() + + if (nstep .le. 1) then + ! when starting up the model, initialize the integrator variables + do fc = 1,num_soilc + c = filter_soilc(fc) + totecosysc_old(c) = totecosysc(c) + totedc_old(c) = totedc(c) + totbgcc_old(c) = totbgcc(c) + nep_timeintegrated(c) = 0._r8 + hr_timeintegrated(c) = 0._r8 + npp_timeintegrated(c) = 0._r8 + ! + ! also initialize the ed-BGC flux variables + ed_to_bgc_this_edts(c) = 0._r8 + ed_to_bgc_last_edts(c) = 0._r8 + ! + cbalance_error_ed(c) = 0._r8 + cbalance_error_bgc(c) = 0._r8 + cbalance_error_total(c) = 0._r8 + end do + endif + + if ( .not. is_beg_curr_day() ) then + ! on CLM (half-hourly) timesteps, integrate the NEP fluxes + do fc = 1,num_soilc + c = filter_soilc(fc) + nep_timeintegrated(c) = nep_timeintegrated(c) + nep(c) * dtime + hr_timeintegrated(c) = hr_timeintegrated(c) + hr(c) * dtime + npp_timeintegrated(c) = npp_timeintegrated(c) + npp_hifreq(c) * dtime + end do + else + ! on ED (daily) timesteps, first integrate the NEP fluxes and add in the daily disturbance flux + do fc = 1,num_soilc + c = filter_soilc(fc) + nep_timeintegrated(c) = nep_timeintegrated(c) + nep(c) * dtime + hr_timeintegrated(c) = hr_timeintegrated(c) + hr(c) * dtime + npp_timeintegrated(c) = npp_timeintegrated(c) + npp_hifreq(c) * dtime + nbp_integrated(c) = nep_timeintegrated(c) - fire_c_to_atm(c) * SHR_CONST_CDAY + seed_rain_flux(c)* SHR_CONST_CDAY + end do + + ! next compare the change in carbon and calculate the error + do fc = 1,num_soilc + c = filter_soilc(fc) + error_ed(c) = totedc(c) - totedc_old(c) - (npp_timeintegrated(c) + seed_rain_flux(c)* SHR_CONST_CDAY - ed_to_bgc_this_edts(c)* SHR_CONST_CDAY - fire_c_to_atm(c) * SHR_CONST_CDAY) + error_bgc(c) = totbgcc(c) - totbgcc_old(c) - (ed_to_bgc_last_edts(c)* SHR_CONST_CDAY - hr_timeintegrated(c)) + error_total(c) = totecosysc(c) - totecosysc_old(c) - (nbp_integrated(c) + ed_to_bgc_last_edts(c)* SHR_CONST_CDAY - ed_to_bgc_this_edts(c)* SHR_CONST_CDAY) + end do + ! + ! put in consistent flux units and send to history so we can keep track of the errors + do fc = 1,num_soilc + c = filter_soilc(fc) + cbalance_error_ed(c) = error_ed(c) / SHR_CONST_CDAY + cbalance_error_bgc(c) = error_bgc(c) / SHR_CONST_CDAY + cbalance_error_total(c) = error_total(c) / SHR_CONST_CDAY + end do + + ! for now, rather than crashing the model, lets just report the largest error to see what we're up against + ! + ! RETURN TO THIS LATER AND ADD A CRASHER IF BALANCE EXCEEDS THRESHOLD + ! + ! max_error_total = 0._r8 + ! do fc = 1,num_soilc + ! c = filter_soilc(fc) + ! if (abs(error_total(c)) .gt. max_error_total) then + ! max_error_ed = abs(error_ed(c)) + ! max_error_bgc = abs(error_bgc(c)) + ! max_error_total = abs(error_total(c)) + ! endif + ! end do + ! write(iulog,*) 'ED_BGC_Carbon_Balancecheck: max_error_ed, max_error_bgc, max_error_total (gC / m2 / day): ', max_error_ed, max_error_bgc, max_error_total + + ! reset the C stock and flux integrators + do fc = 1,num_soilc + c = filter_soilc(fc) + totecosysc_old(c) = totecosysc(c) + totedc_old(c) = totedc(c) + totbgcc_old(c) = totbgcc(c) + nep_timeintegrated(c) = 0._r8 + npp_timeintegrated(c) = 0._r8 + hr_timeintegrated(c) = 0._r8 + end do + + endif + + end associate + + end subroutine ED_BGC_Carbon_Balancecheck + end module EDCLMLinkMod diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index 2aedbc276b..47934da698 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -298,6 +298,7 @@ module EDTypesMod real(r8) :: seed_decay(numpft_ed) ! seed decay in KgC/m2/year real(r8) :: seed_germination(numpft_ed) ! germination rate of seed pool in KgC/m2/year real(r8) :: dseed_dt(numpft_ed) + real(r8) :: seed_rain_flux(numpft_ed) ! flux of seeds from exterior KgC/m2/year (needed for C balance purposes) ! PHOTOSYNTHESIS real(r8) :: psn_z(nclmax,numpft_ed,nlevcan_ed) ! carbon assimilation in each canopy layer, pft, and leaf layer. umolC/m2/s @@ -421,13 +422,14 @@ module EDTypesMod integer :: dleafondate ! doy of leaf on drought:- integer :: dleafoffdate ! doy of leaf on drought:- real(r8) :: water_memory(10) ! last 10 days of soil moisture memory... - real(r8) :: cwd_ag_burned(ncwd) - real(r8) :: leaf_litter_burned(numpft_ed) ! FIRE real(r8) :: acc_ni ! daily nesterov index accumulating over time. real(r8) :: ab ! daily burnt area: m2 real(r8) :: frac_burnt ! fraction of soil burnt in this day. + real(r8) :: total_burn_flux_to_atm ! total carbon burnt to the atmosphere in this day. KgC/site + real(r8) :: cwd_ag_burned(ncwd) + real(r8) :: leaf_litter_burned(numpft_ed) end type ed_site_type