From 2b650b77bba10d020b68c45c547ac1aefdd54ac5 Mon Sep 17 00:00:00 2001 From: Ben Andre Date: Mon, 14 Dec 2015 19:02:06 -0700 Subject: [PATCH 01/22] Import ed_v0.1.0 unstable science changes. Created a patch file of changes Rosie considered unstable science from the ed_v0.1.0 branch: svn diff -r68491:r72748 \ https://svn-ccsm-models.cgd.ucar.edu/clm2/branches/ed_v0.1.0/ > \ ed_v0.1.0-buggy-science.patch Patch applied cleanly. Testing: Branch has not been compiled or tested. --- biogeochem/EDBGCDynMod.F90 | 355 ++++++++++++++++++++++++++ biogeochem/EDPatchDynamicsMod.F90 | 2 +- biogeochem/EDPhysiologyMod.F90 | 18 +- biogeophys/EDSurfaceAlbedoMod.F90 | 4 +- main/EDCLMLinkMod.F90 | 406 +++++++++++++++++++++++++++++- main/EDMainMod.F90 | 2 +- 6 files changed, 771 insertions(+), 16 deletions(-) create mode 100644 biogeochem/EDBGCDynMod.F90 diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 new file mode 100644 index 0000000000..6ba670fbae --- /dev/null +++ b/biogeochem/EDBGCDynMod.F90 @@ -0,0 +1,355 @@ +module EDBGCDynMod + +! Interface from ED calls to CLM belowground biogeochemistry module + + use shr_kind_mod , only : r8 => shr_kind_r8 + use clm_varctl , only : use_c13, use_c14, use_ed + use decompMod , only : bounds_type + use perf_mod , only : t_startf, t_stopf + use clm_varctl , only : use_century_decomp, use_nitrif_denitrif + use CNVegStateType , only : cnveg_state_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonFluxType , only : cnveg_carbonflux_type + use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type + use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type + use SoilBiogeochemStateType , only : soilbiogeochem_state_type + use SoilBiogeochemCarbonStateType , only : soilbiogeochem_carbonstate_type + use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type + use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type + use SoilBiogeochemNitrogenFluxType , only : soilbiogeochem_nitrogenflux_type + use EDCLMLinkMod , only : ed_clm_type + use CanopyStateType , only : canopystate_type + use SoilStateType , only : soilstate_type + use SoilHydrologyType , only : soilhydrology_type + use TemperatureType , only : temperature_type + use WaterstateType , only : waterstate_type + use WaterfluxType , only : waterflux_type + use atm2lndType , only : atm2lnd_type + use SoilStateType , only : soilstate_type + use ch4Mod , only : ch4_type + + + ! public :: EDBGCDynInit ! BGC dynamics: initialization + public :: EDBGCDyn ! BGC Dynamics + public :: EDBGCDynSummary ! BGC dynamics: summary + +contains + + + !----------------------------------------------------------------------- + subroutine EDBGCDyn(bounds, & + num_soilc, filter_soilc, num_soilp, filter_soilp, num_pcropp, filter_pcropp, doalb, & + cnveg_state_inst, & + cnveg_carbonflux_inst, cnveg_carbonstate_inst, & + ed_clm_inst, & + soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & + soilbiogeochem_state_inst, & + soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, & + c13_soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonflux_inst, & + c14_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonflux_inst, & + atm2lnd_inst, waterstate_inst, waterflux_inst, & + canopystate_inst, soilstate_inst, temperature_inst, crop_inst, ch4_inst) + ! + ! !DESCRIPTION: + + ! + ! !USES: + use clm_varpar , only: crop_prog, nlevgrnd, nlevdecomp_full + use clm_varpar , only: nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools + use subgridAveMod , only: p2c + use CropType , only: crop_type + use CNNDynamicsMod , only: CNNDeposition,CNNFixation, CNNFert, CNSoyfix + use CNMRespMod , only: CNMResp + use CNPhenologyMod , only: CNPhenology + use CNGRespMod , only: CNGResp + use CNFireMod , only: CNFireArea, CNFireFluxes + use CNCIsoFluxMod , only: CIsoFlux1, CIsoFlux2, CIsoFlux2h, CIsoFlux3 + use CNC14DecayMod , only: C14Decay + use CNWoodProductsMod , only: CNWoodProducts + use CNCStateUpdate1Mod , only: CStateUpdate1,CStateUpdate0 + use CNCStateUpdate2Mod , only: CStateUpdate2, CStateUpdate2h + use CNCStateUpdate3Mod , only: CStateUpdate3 + use CNNStateUpdate1Mod , only: NStateUpdate1 + use CNNStateUpdate2Mod , only: NStateUpdate2, NStateUpdate2h + use CNGapMortalityMod , only: CNGapMortality + use dynHarvestMod , only: CNHarvest + use SoilBiogeochemDecompCascadeBGCMod , only: decomp_rate_constants_bgc + use SoilBiogeochemDecompCascadeCNMod , only: decomp_rate_constants_cn + use SoilBiogeochemCompetitionMod , only: SoilBiogeochemCompetition + use SoilBiogeochemDecompMod , only: SoilBiogeochemDecomp + use SoilBiogeochemLittVertTranspMod , only: SoilBiogeochemLittVertTransp + use SoilBiogeochemPotentialMod , only: SoilBiogeochemPotential + use SoilBiogeochemVerticalProfileMod , only: SoilBiogeochemVerticalProfile + use SoilBiogeochemNitrifDenitrifMod , only: SoilBiogeochemNitrifDenitrif + use SoilBiogeochemNStateUpdate1Mod , only: SoilBiogeochemNStateUpdate1 + ! + ! !ARGUMENTS: + 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 + integer , intent(in) :: num_soilp ! number of soil patches in filter + integer , intent(in) :: filter_soilp(:) ! filter for soil patches + integer , intent(in) :: num_pcropp ! number of prog. crop patches in filter + integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches + logical , intent(in) :: doalb ! true = surface albedo calculation time step + type(cnveg_state_type) , intent(inout) :: cnveg_state_inst + type(cnveg_carbonflux_type) , intent(inout) :: cnveg_carbonflux_inst + type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst + type(ed_clm_type) , intent(inout) :: ed_clm_inst + type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: c13_soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: c13_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst + type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst + type(atm2lnd_type) , intent(in) :: atm2lnd_inst + type(waterstate_type) , intent(in) :: waterstate_inst + type(waterflux_type) , intent(in) :: waterflux_inst + type(canopystate_type) , intent(in) :: canopystate_inst + type(soilstate_type) , intent(in) :: soilstate_inst + type(temperature_type) , intent(inout) :: temperature_inst + type(crop_type) , intent(in) :: crop_inst + type(ch4_type) , intent(in) :: ch4_inst + ! + ! !LOCAL VARIABLES: + real(r8):: cn_decomp_pools(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools) + real(r8):: p_decomp_cpool_loss(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential C loss from one pool to another + real(r8):: pmnf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential mineral N flux, from one pool to another + real(r8):: arepr(bounds%begp:bounds%endp) ! reproduction allocation coefficient (only used for crop_prog) + real(r8):: aroot(bounds%begp:bounds%endp) ! root allocation coefficient (only used for crop_prog) + integer :: begp,endp + integer :: begc,endc + !----------------------------------------------------------------------- + + begp = bounds%begp; endp = bounds%endp + begc = bounds%begc; endc = bounds%endc + + !real(r8) , intent(in) :: rootfr_patch(bounds%begp:, 1:) + !integer , intent(in) :: altmax_lastyear_indx_col(bounds%begc:) ! frost table depth (m) + + associate( & + rootfr_patch => soilstate_inst%rootfr_patch , & ! fraction of roots in each soil layer (nlevgrnd) + altmax_lastyear_indx_col => canopystate_inst%altmax_lastyear_indx_col , & ! frost table depth (m) + laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit projected leaf area index + laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded projected leaf area index + frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] + frac_veg_nosno_alb => canopystate_inst%frac_veg_nosno_alb_patch , & ! Output: [integer (:) ] frac of vegetation not covered by snow [-] + tlai => canopystate_inst%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow + tsai => canopystate_inst%tsai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index, no burying by snow + elai => canopystate_inst%elai_patch , & ! Output: [real(r8) (:) ] one-sided leaf area index with burying by snow + esai => canopystate_inst%esai_patch , & ! Output: [real(r8) (:) ] one-sided stem area index with burying by snow + htop => canopystate_inst%htop_patch , & ! Output: [real(r8) (:) ] canopy top (m) + hbot => canopystate_inst%hbot_patch & ! Output: [real(r8) (:) ] canopy bottom (m) + ) + + ! -------------------------------------------------- + ! zero the column-level C and N fluxes + ! -------------------------------------------------- + + call t_startf('BGCZero') + + call soilbiogeochem_carbonflux_inst%SetValues( & + num_soilc, filter_soilc, 0._r8) + if ( use_c13 ) then + call c13_soilbiogeochem_carbonflux_inst%SetValues( & + num_soilc, filter_soilc, 0._r8) + end if + if ( use_c14 ) then + call c14_soilbiogeochem_carbonflux_inst%SetValues( & + num_soilc, filter_soilc, 0._r8) + end if + + call t_stopf('BGCZero') + + ! -------------------------------------------------- + ! Nitrogen Deposition, Fixation and Respiration + ! -------------------------------------------------- + + ! call t_startf('CNDeposition') + ! call CNNDeposition(bounds, & + ! atm2lnd_inst, soilbiogeochem_nitrogenflux_inst) + ! call t_stopf('CNDeposition') + + + ! if (crop_prog) then + ! call CNNFert(bounds, num_soilc,filter_soilc, & + ! cnveg_nitrogenflux_inst, soilbiogeochem_nitrogenflux_inst) + + ! call CNSoyfix (bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & + ! waterstate_inst, crop_inst, cnveg_state_inst, cnveg_nitrogenflux_inst , & + ! soilbiogeochem_state_inst, soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst) + ! end if + + !-------------------------------------------- + ! Soil Biogeochemistry + !-------------------------------------------- + + if (use_century_decomp) then + call decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & + canopystate_inst, soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) + else + call decomp_rate_constants_cn(bounds, num_soilc, filter_soilc, & + canopystate_inst, soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) + end if + + ! calculate potential decomp rates and total immobilization demand (previously inlined in CNDecompAlloc) + call SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & + soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & + soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & + cn_decomp_pools=cn_decomp_pools(begc:endc,1:nlevdecomp,1:ndecomp_pools), & + p_decomp_cpool_loss=p_decomp_cpool_loss(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) + + + !-------------------------------------------- + ! Resolve the competition between plants and soil heterotrophs + ! for available soil mineral N resource + !-------------------------------------------- + ! will add this back in when integrtating hte nutirent cycles + + + !-------------------------------------------- + ! Calculate litter and soil decomposition rate + !-------------------------------------------- + + ! Calculation of actual immobilization and decomp rates, following + ! resolution of plant/heterotroph competition for mineral N (previously inlined in CNDecompAllocation in CNDecompMod) + + call t_startf('SoilBiogeochemDecomp') + + call SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, & + soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & + soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & + cn_decomp_pools=cn_decomp_pools(begc:endc,1:nlevdecomp,1:ndecomp_pools), & + p_decomp_cpool_loss=p_decomp_cpool_loss(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) + + call t_stopf('SoilBiogeochemDecomp') + + + !-------------------------------------------- + ! Update1 + !-------------------------------------------- + + call t_startf('BNGCUpdate1') + + + ! Update all prognostic carbon state variables (except for gap-phase mortality and fire fluxes) + call CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & + cnveg_state_inst, cnveg_carbonflux_inst, cnveg_carbonstate_inst, & + ed_clm_inst, & + soilbiogeochem_carbonflux_inst) + + call t_stopf('BNGCUpdate1') + + !-------------------------------------------- + ! Calculate vertical mixing of soil and litter pools + !-------------------------------------------- + + call t_startf('SoilBiogeochemLittVertTransp') + + call SoilBiogeochemLittVertTransp(bounds, num_soilc, filter_soilc, & + canopystate_inst, soilbiogeochem_state_inst, & + soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & + c13_soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonflux_inst, & + c14_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonflux_inst, & + soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst) + + call t_stopf('SoilBiogeochemLittVertTransp') + + end associate + + end subroutine EDBGCDyn + + + !----------------------------------------------------------------------- + subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & + soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & + c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & + c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & + soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) + ! + ! !DESCRIPTION: + ! Call to all CN and SoilBiogeochem summary routines + ! + ! !USES: + use clm_varpar , only: ndecomp_cascade_transitions + use CNPrecisionControlMod , only: CNPrecisionControl + use SoilBiogeochemPrecisionControlMod , only: SoilBiogeochemPrecisionControl + ! + ! !ARGUMENTS: + 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 + integer , intent(in) :: num_soilp ! number of soil patches in filter + integer , intent(in) :: filter_soilp(:) ! filter for soil patches + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: c13_soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: c13_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst + type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst + ! + ! !LOCAL VARIABLES: + integer :: begc,endc + !----------------------------------------------------------------------- + + begc = bounds%begc; endc= bounds%endc + + ! Call to all summary routines + + call t_startf('BGCsum') + + ! Set controls on very low values in critical state variables + + call SoilBiogeochemPrecisionControl(num_soilc, filter_soilc, & + soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonstate_inst, & + c14_soilbiogeochem_carbonstate_inst,soilbiogeochem_nitrogenstate_inst) + + ! Note - all summary updates to cnveg_carbonstate_inst and cnveg_carbonflux_inst are done in + ! soilbiogeochem_carbonstate_inst%summary and CNVeg_carbonstate_inst%summary + + ! ---------------------------------------------- + ! soilbiogeochem carbon/nitrogen state summary + ! ---------------------------------------------- + + call soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) + if ( use_c13 ) then + call c13_soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) + end if + if ( use_c14 ) then + call c14_soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) + end if + ! call soilbiogeochem_nitrogenstate_inst%summary(bounds, num_soilc, filter_soilc) + + ! ---------------------------------------------- + ! soilbiogeochem carbon/nitrogen flux summary + ! ---------------------------------------------- + + call soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) + if ( use_c13 ) then + call c13_soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) + end if + if ( use_c14 ) then + call c14_soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) + end if + ! call soilbiogeochem_nitrogenflux_inst%Summary(bounds, num_soilc, filter_soilc) + + ! ---------------------------------------------- + ! ed veg carbon state summary + ! ---------------------------------------------- + + ! ---------------------------------------------- + ! ed veg carbon/nitrogen flux summary + ! ---------------------------------------------- + + call t_stopf('BGCsum') + + end subroutine EDBGCDynSummary + +end module EDBGCDynMod diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index 397606ced1..d5ed6bc157 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -1020,7 +1020,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 969f8481b1..4277ab17cf 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -120,14 +120,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 @@ -1174,4 +1174,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 a061c5c275..e6c007eaba 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -849,10 +849,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 ! write(iulog,*) 'first layer of lai_change 2 3',lai_change(1,1,1:3) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index b0b559e9ca..bc0c3e0f99 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -9,12 +9,19 @@ module EDCLMLinkMod use decompMod , only : bounds_type use clm_varpar , only : nclmax, nlevcan_ed, numpft, numcft use clm_varctl , only : iulog - use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type + use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, ncwd + use CanopyStateType , only : canopystate_type + use clm_varctl , only : use_vertsoilc + ! 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 @@ -71,6 +78,16 @@ module EDCLMLinkMod 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 + ! 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 + contains ! Public routines @@ -84,7 +101,8 @@ module EDCLMLinkMod 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 @@ -115,7 +133,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 @@ -123,9 +142,12 @@ subroutine InitAllocate(this, bounds) ! ! !LOCAL VARIABLES: integer :: begp,endp + integer :: begc,endc !bounds !------------------------------------------------------------------------ begp = bounds%begp; endp = bounds%endp + begc = bounds%begc; endc = bounds%endc + allocate(this%trimming_patch (begp:endp)) ; this%trimming_patch (:) = 0.0_r8 allocate(this%canopy_spread_patch (begp:endp)) ; this%canopy_spread_patch (:) = 0.0_r8 @@ -171,6 +193,15 @@ subroutine InitAllocate(this, bounds) allocate(this%gpp_patch (begp:endp)) ; this%gpp_patch (:) = nan allocate(this%npp_patch (begp:endp)) ; this%npp_patch (:) = nan + 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 + end subroutine InitAllocate !------------------------------------------------------------------------ @@ -379,6 +410,32 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='net primary production', & ptr_patch=this%npp_patch) + 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) + + 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) + + 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) + + call hist_addfld_decomp (fname='leaf_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='leaf_prof', & + ptr_col=this%leaf_prof_col) + + call hist_addfld_decomp (fname='croot_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='croot_prof', & + ptr_col=this%croot_prof_col) + + call hist_addfld_decomp (fname='stem_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='stem_prof', & + ptr_col=this%stem_prof_col) + + + end subroutine InitHistory !----------------------------------------------------------------------- @@ -717,7 +774,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 @@ -1426,4 +1486,342 @@ 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 shr_const_mod, only: SHR_CONST_CDAY + 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./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./dzsoi_decomp(1) + stem_prof(c,1) = 1./dzsoi_decomp(1) + 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) + 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 + + end associate + end subroutine flux_into_litter_pools + + + end module EDCLMLinkMod diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 66e929a386..ace3cfd73c 100755 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -285,7 +285,7 @@ subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature enddo - write(6,*)'DEBUG18: calling non_canopy_derivs with pno= ',currentPatch%clm_pno + call non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, waterstate_inst ) !update state variables simultaneously according to derivatives for this time period. From d3f63293ca8a422ec7ef6c5065f1fd361d6b33e2 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 17 Feb 2016 16:24:50 -0800 Subject: [PATCH 02/22] added restart variables to edclmlinkmod --- main/EDCLMLinkMod.F90 | 47 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index df6eb3a48e..11f2b63c4a 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -614,6 +614,8 @@ subroutine Restart ( this, bounds, ncid, flag ) ! ! !LOCAL VARIABLES: logical :: readvar + character(LEN=3) :: istr1 + integer :: k !------------------------------------------------------------------------ call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & @@ -626,12 +628,55 @@ subroutine Restart ( this, bounds, ncid, flag ) 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) + interpinic_flag='interp', readvar=readvar, data=this%deadstemc_patch) 4 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) + 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=this%ED_c_to_litr_lab_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=this%ED_c_to_litr_cel_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=this%ED_c_to_litr_lig_c_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=this%leaf_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=this%croot_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=this%stem_prof_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=this%ED_c_to_litr_lab_c_col) + + do k = 1, numpft_ed + write(istr1,"(I3.3)") 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=this%froot_prof_col(:,k,:)) + end do + end subroutine Restart !----------------------------------------------------------------------- From 5581373d20fa799eee94cc53c2788fd58ed3a7d1 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 17 Feb 2016 16:32:39 -0800 Subject: [PATCH 03/22] typo bugfix --- main/EDCLMLinkMod.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 11f2b63c4a..6170c3434c 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -628,7 +628,7 @@ subroutine Restart ( this, bounds, ncid, flag ) 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) 4 + 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='', & From b5a3d347e8df0eca914830611a6c231851efb9a4 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Thu, 18 Feb 2016 14:16:51 -0800 Subject: [PATCH 04/22] removed profile variables from restart file as I don't believe they are needed --- main/EDCLMLinkMod.F90 | 60 ++++++++++++++++++++++--------------------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 6170c3434c..0a2ec3876a 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -605,6 +605,8 @@ subroutine Restart ( this, bounds, ncid, flag ) ! !USES: use restUtilMod use ncdio_pio + ! use EDtypesMod , only : numpft_ed + ! ! !ARGUMENTS: class (ed_clm_type) :: this @@ -614,8 +616,8 @@ subroutine Restart ( this, bounds, ncid, flag ) ! ! !LOCAL VARIABLES: logical :: readvar - character(LEN=3) :: istr1 - integer :: k + ! character(LEN=3) :: istr1 + ! integer :: k !------------------------------------------------------------------------ call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & @@ -649,33 +651,33 @@ subroutine Restart ( this, bounds, ncid, flag ) long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%ED_c_to_litr_lig_c_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=this%leaf_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=this%croot_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=this%stem_prof_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=this%ED_c_to_litr_lab_c_col) - - do k = 1, numpft_ed - write(istr1,"(I3.3)") 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=this%froot_prof_col(:,k,:)) - end do + ! 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=this%leaf_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=this%croot_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=this%stem_prof_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=this%ED_c_to_litr_lab_c_col) + + ! do k = 1, numpft_ed + ! write(istr1,"(I3.3)") 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=this%froot_prof_col(:,k,:)) + ! end do end subroutine Restart From fd4a9e8708a3e775be35d7ad6d1cadf2abf28c14 Mon Sep 17 00:00:00 2001 From: ckoven Date: Thu, 18 Feb 2016 22:59:12 -0800 Subject: [PATCH 05/22] working on restarts still.. --- main/EDCLMLinkMod.F90 | 87 ++++++++++++++++++++++++++----------------- 1 file changed, 52 insertions(+), 35 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 0a2ec3876a..d75c8208be 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -616,40 +616,62 @@ 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='livestemn', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%livestemn_patch) - - 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=this%ED_c_to_litr_lab_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=this%ED_c_to_litr_cel_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=this%ED_c_to_litr_lig_c_col) + ! 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) + + 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) + 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) + end if ! call restartvar(ncid=ncid, flag=flag, varname='leaf_prof_col', xtype=ncd_double, & ! dim1name='column', dim2name='levgrnd', switchdim=.true., & @@ -666,11 +688,6 @@ subroutine Restart ( this, bounds, ncid, flag ) ! long_name='', units='', & ! interpinic_flag='interp', readvar=readvar, data=this%stem_prof_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=this%ED_c_to_litr_lab_c_col) - ! do k = 1, numpft_ed ! write(istr1,"(I3.3)") k ! call restartvar(ncid=ncid, flag=flag, varname='froot_prof_col_PFT'//istr1, xtype=ncd_double, & From 71b347c0c4e99611837cc6ef59b1c52044a547f6 Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 19 Feb 2016 17:04:24 -0800 Subject: [PATCH 06/22] more fixes on ED/BGC intersection. model should now pass tests --- main/EDCLMLinkMod.F90 | 72 ++++++++++++++++++++++++++++++------------- 1 file changed, 50 insertions(+), 22 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index d75c8208be..ee17c99e2c 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -656,6 +656,33 @@ subroutine Restart ( this, bounds, ncid, flag ) 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, & @@ -671,30 +698,31 @@ subroutine Restart ( this, bounds, ncid, flag ) 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='leaf_prof_col', xtype=ncd_double, & - ! dim1name='column', dim2name='levgrnd', switchdim=.true., & - ! long_name='', units='', & - ! interpinic_flag='interp', readvar=readvar, data=this%leaf_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=this%croot_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=this%stem_prof_col) - - ! do k = 1, numpft_ed - ! write(istr1,"(I3.3)") 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=this%froot_prof_col(:,k,:)) - ! end do end subroutine Restart From 0d3ce80353e23801c136646d61be4d45c645689e Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Fri, 26 Feb 2016 15:42:40 -0800 Subject: [PATCH 07/22] added first instance of NEP and NBP fluxes --- biogeochem/EDBGCDynMod.F90 | 13 ++++-- main/EDCLMLinkMod.F90 | 84 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+), 3 deletions(-) diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 index 6ba670fbae..d90b36c34a 100644 --- a/biogeochem/EDBGCDynMod.F90 +++ b/biogeochem/EDBGCDynMod.F90 @@ -27,7 +27,7 @@ module EDBGCDynMod use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type use ch4Mod , only : ch4_type - + use EDtypesMod , only : ed_site_type ! public :: EDBGCDynInit ! BGC dynamics: initialization public :: EDBGCDyn ! BGC Dynamics @@ -270,10 +270,12 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & - soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) + soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, & + ed_clm_inst, ed_allsites_inst) ! ! !DESCRIPTION: ! Call to all CN and SoilBiogeochem summary routines + ! also aggregate production and decomposition fluxes to whole-ecosystem balance fluxes ! ! !USES: use clm_varpar , only: ndecomp_cascade_transitions @@ -294,6 +296,8 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst + type(ed_clm_type) , intent(inout) :: ed_clm_inst + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) ! ! !LOCAL VARIABLES: integer :: begc,endc @@ -341,8 +345,11 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so ! call soilbiogeochem_nitrogenflux_inst%Summary(bounds, num_soilc, filter_soilc) ! ---------------------------------------------- - ! ed veg carbon state summary + ! ed veg carbon flux summary ! ---------------------------------------------- + + call ed_clm_inst%Summary(bounds, numsoilc, filter_soilc, num_soilp, filter_soilp, & + ed_allsites_inst(bounds%begg:bounds%endg), soilbiogeochem_carbonflux_inst) ! ---------------------------------------------- ! ed veg carbon/nitrogen flux summary diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index ee17c99e2c..e2ceb437b1 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -16,6 +16,8 @@ module EDCLMLinkMod 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 clm_time_manager , only : get_step_size ! implicit none @@ -123,6 +125,11 @@ module EDCLMLinkMod 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 :: 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 + contains ! Public routines @@ -130,6 +137,7 @@ module EDCLMLinkMod procedure , public :: Restart procedure , public :: SetValues procedure , public :: ed_clm_link + procedure , public :: Summary ! Private routines procedure , private :: ed_clm_leaf_area_profile @@ -239,6 +247,10 @@ subroutine InitAllocate(this, bounds) 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%nep_col (begc:endc)) ; this%nep_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%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 @@ -476,6 +488,21 @@ 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%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 at high frequency', units='gC/m^2/s', & + avgflag='A', long_name='net primary production at high frequency', & + ptr_col=this%npp_hifreq_col) + !!! carbon fluxes into soil grid (dimensioned depth x column) 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', & @@ -2184,7 +2211,64 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch end associate end subroutine flux_into_litter_pools + + !------------------------------------------------------------------------ + subroutine Summary(this, bounds, numsoilc, filter_soilc, num_soilp, filter_soilp, & + ed_allsites_inst, soilbiogeochem_carbonflux_inst) + + ! Summarize the combined production and decomposition fluxes into net fluxes + ! Written by Charlie Koven, Feb 2016 + + 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 + integer , intent(in) :: num_soilp ! number of soil patches in filter + integer , intent(in) :: filter_soilp(:) ! filter for soil patches + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + + real(r8) :: npp_hifreq_col(bounds%begc:bounds%endc ! column-level, high frequency NPP + real(r8) :: dt ! radiation time step (seconds) + + associate(& + hr => soilbiogeochem_carbonflux_inst%hr_col & ! Output: (gC/m2/s) total heterotrophic respiration + npp_hifreq => this%npp_hifreq_col + nep => this%nep_col + nbp => this%nbp_col + ) + + ! set time steps + dt = real( get_step_size(), r8 ) + + do c = bounds%begc,bounds%endc + npp_hifreq(c) = 0._r8 + end do + + 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 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + npp_hifreq(cc) = npp_hifreq(cc) + currentpatch%npp_clm * 1e3 / ( AREA * dt) + enddo !currentCohort + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + ! calculate NEP and NBP fluxes. + !!!! CDK FEB/26/2016: THIS IS IMPORTANT AND NEEDS TO CHANGE AS IT IS ONLY A PLACEHOLDER. + !!!! NEP AND NBP ARE BOTH THE SAME RIGHT NOW BECAUSE I DON'T KNOW HOW TO ADD IN THE FIRE, DISTURBANCE, ETC FLUXES INTO THE NBP FLUX YET + do fc = 1,num_soilc + c = filter_soilc(fc) + nep(c) = npp_hifreq(c) - hr(c) + nbp(c) = npp_hifreq(c) - hr(c) + end do + end module EDCLMLinkMod From 935b2ad825ceee03ae69de43cb74108d66292316 Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 26 Feb 2016 16:15:50 -0800 Subject: [PATCH 08/22] first set of bugfixes on NEP/NBP code --- main/EDCLMLinkMod.F90 | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index e2ceb437b1..8fb5600712 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -2214,12 +2214,13 @@ end subroutine flux_into_litter_pools !------------------------------------------------------------------------ - subroutine Summary(this, bounds, numsoilc, filter_soilc, num_soilp, filter_soilp, & + subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & ed_allsites_inst, soilbiogeochem_carbonflux_inst) ! Summarize the combined production and decomposition fluxes into net fluxes ! Written by Charlie Koven, Feb 2016 + 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 @@ -2228,14 +2229,18 @@ subroutine Summary(this, bounds, numsoilc, filter_soilc, num_soilp, filter_soilp type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst - real(r8) :: npp_hifreq_col(bounds%begc:bounds%endc ! column-level, high frequency NPP + 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 + type(ed_site_type), pointer :: cs + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort associate(& - hr => soilbiogeochem_carbonflux_inst%hr_col & ! Output: (gC/m2/s) total heterotrophic respiration - npp_hifreq => this%npp_hifreq_col - nep => this%nep_col - nbp => this%nbp_col + hr => soilbiogeochem_carbonflux_inst%hr_col, & ! Output: (gC/m2/s) total heterotrophic respiration + npp_hifreq => this%npp_hifreq_col, & + nep => this%nep_col, & + nbp => this%nbp_col & ) ! set time steps @@ -2246,14 +2251,14 @@ subroutine Summary(this, bounds, numsoilc, filter_soilc, num_soilp, filter_soilp end do do g = bounds%begg,bounds%endg - if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + if (ed_allsites_inst(g)%istheresoil) then currentPatch => ed_allsites_inst(g)%oldest_patch do while(associated(currentPatch)) - cs => currentpatch%siteptr + cs => currentPatch%siteptr cc = cs%clmcolumn currentCohort => currentPatch%tallest do while(associated(currentCohort)) - npp_hifreq(cc) = npp_hifreq(cc) + currentpatch%npp_clm * 1e3 / ( AREA * dt) + npp_hifreq(cc) = npp_hifreq(cc) + currentCohort%npp_clm * 1e3 * currentCohort%n / ( AREA * dt) enddo !currentCohort currentPatch => currentPatch%younger end do !currentPatch @@ -2269,6 +2274,9 @@ subroutine Summary(this, bounds, numsoilc, filter_soilc, num_soilp, filter_soilp nbp(c) = npp_hifreq(c) - hr(c) end do + end associate + +end subroutine Summary end module EDCLMLinkMod From a779c53b94085a2f8bb66deb4049ce113e20faa9 Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 26 Feb 2016 18:21:13 -0800 Subject: [PATCH 09/22] second set of bugfixes on NEP/NBP code --- biogeochem/EDBGCDynMod.F90 | 2 +- main/EDCLMLinkMod.F90 | 146 +++++++++++++++++++++---------------- 2 files changed, 83 insertions(+), 65 deletions(-) diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 index d90b36c34a..9e23408432 100644 --- a/biogeochem/EDBGCDynMod.F90 +++ b/biogeochem/EDBGCDynMod.F90 @@ -348,7 +348,7 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so ! ed veg carbon flux summary ! ---------------------------------------------- - call ed_clm_inst%Summary(bounds, numsoilc, filter_soilc, num_soilp, filter_soilp, & + call ed_clm_inst%Summary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & ed_allsites_inst(bounds%begg:bounds%endg), soilbiogeochem_carbonflux_inst) ! ---------------------------------------------- diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 8fb5600712..697e625461 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -499,9 +499,9 @@ subroutine InitHistory(this, bounds) ptr_col=this%nbp_col) this%npp_hifreq_col(begc:endc) = spval - call hist_addfld1d (fname='NPP at high frequency', units='gC/m^2/s', & + 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) + ptr_col=this%npp_hifreq_col,default='inactive') !!! carbon fluxes into soil grid (dimensioned depth x column) call hist_addfld_decomp (fname='ED_c_to_litr_lab_c', units='gC/m^2/s', type2d='levdcmp', & @@ -518,15 +518,15 @@ subroutine InitHistory(this, bounds) call hist_addfld_decomp (fname='leaf_prof', units='1/m', type2d='levdcmp', & avgflag='A', long_name='leaf_prof', & - ptr_col=this%leaf_prof_col) + ptr_col=this%leaf_prof_col,default='inactive') call hist_addfld_decomp (fname='croot_prof', units='1/m', type2d='levdcmp', & avgflag='A', long_name='croot_prof', & - ptr_col=this%croot_prof_col) + ptr_col=this%croot_prof_col,default='inactive') call hist_addfld_decomp (fname='stem_prof', units='1/m', type2d='levdcmp', & avgflag='A', long_name='stem_prof', & - ptr_col=this%stem_prof_col) + ptr_col=this%stem_prof_col,default='inactive') ! Carbon Flux (grid dimension x scpf) @@ -2219,64 +2219,82 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil ! Summarize the combined production and decomposition fluxes into net fluxes ! Written by Charlie Koven, Feb 2016 - - 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 - integer , intent(in) :: num_soilp ! number of soil patches in filter - integer , intent(in) :: filter_soilp(:) ! filter for soil patches - type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) - type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst - - 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 - type(ed_site_type), pointer :: cs - type (ed_patch_type) , pointer :: currentPatch - type (ed_cohort_type) , pointer :: currentCohort - - associate(& - hr => soilbiogeochem_carbonflux_inst%hr_col, & ! Output: (gC/m2/s) total heterotrophic respiration - npp_hifreq => this%npp_hifreq_col, & - nep => this%nep_col, & - nbp => this%nbp_col & - ) - - ! set time steps - dt = real( get_step_size(), r8 ) - - do c = bounds%begc,bounds%endc - npp_hifreq(c) = 0._r8 - end do - - do g = bounds%begg,bounds%endg - if (ed_allsites_inst(g)%istheresoil) then - currentPatch => ed_allsites_inst(g)%oldest_patch - do while(associated(currentPatch)) - cs => currentPatch%siteptr - cc = cs%clmcolumn - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - npp_hifreq(cc) = npp_hifreq(cc) + currentCohort%npp_clm * 1e3 * currentCohort%n / ( AREA * dt) - enddo !currentCohort - currentPatch => currentPatch%younger - end do !currentPatch - end if - end do - - ! calculate NEP and NBP fluxes. - !!!! CDK FEB/26/2016: THIS IS IMPORTANT AND NEEDS TO CHANGE AS IT IS ONLY A PLACEHOLDER. - !!!! NEP AND NBP ARE BOTH THE SAME RIGHT NOW BECAUSE I DON'T KNOW HOW TO ADD IN THE FIRE, DISTURBANCE, ETC FLUXES INTO THE NBP FLUX YET - do fc = 1,num_soilc - c = filter_soilc(fc) - nep(c) = npp_hifreq(c) - hr(c) - nbp(c) = npp_hifreq(c) - hr(c) - end do - - end associate - -end subroutine Summary - + ! + ! !USES: + use ColumnType , only : col + use LandunitType , only : lun + use landunit_varcon , only : istsoil + + 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 + integer , intent(in) :: num_soilp ! number of soil patches in filter + integer , intent(in) :: filter_soilp(:) ! filter for soil patches + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + + 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, & ! Output: (gC/m2/s) total heterotrophic respiration + npp_hifreq => this%npp_hifreq_col, & + nep => this%nep_col, & + nbp => this%nbp_col & + ) + + ! set time steps + dt = real( get_step_size(), r8 ) + + do c = bounds%begc,bounds%endc + npp_hifreq(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 + currentPatch => ed_allsites_inst(g)%oldest_patch + do while(associated(currentPatch)) + cs => currentPatch%siteptr + cc = cs%clmcolumn + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + npp_hifreq(cc) = npp_hifreq(cc) + currentCohort%npp_clm * 1e3 * currentCohort%n / ( AREA * dt) + currentCohort => currentCohort%shorter + enddo !currentCohort + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + ! calculate NEP and NBP fluxes. + !!!! CDK FEB/26/2016: THIS IS IMPORTANT AND NEEDS TO CHANGE AS IT IS ONLY A PLACEHOLDER. + !!!! NEP AND NBP ARE BOTH THE SAME RIGHT NOW BECAUSE I DON'T KNOW HOW TO ADD IN THE FIRE, DISTURBANCE, ETC FLUXES INTO THE NBP FLUX YET + do fc = 1,num_soilc + c = filter_soilc(fc) + nep(c) = npp_hifreq(c) - hr(c) + nbp(c) = npp_hifreq(c) - hr(c) + end do + + end associate + + end subroutine Summary + + end module EDCLMLinkMod From 339f6c47a093bd3227b793a2f75c9c9f32413dfd Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Sat, 27 Feb 2016 17:51:57 -0800 Subject: [PATCH 10/22] added code to aggregate column-level fire carbon fluxes and use these as the difference between NEP and NBP --- biogeochem/EDPatchDynamicsMod.F90 | 13 ++++++++++++- main/EDCLMLinkMod.F90 | 26 ++++++++++++++++++-------- main/EDTypesMod.F90 | 5 +++-- 3 files changed, 33 insertions(+), 11 deletions(-) diff --git a/biogeochem/EDPatchDynamicsMod.F90 b/biogeochem/EDPatchDynamicsMod.F90 index dab3b7fe5e..1f3f1c618c 100755 --- a/biogeochem/EDPatchDynamicsMod.F90 +++ b/biogeochem/EDPatchDynamicsMod.F90 @@ -198,8 +198,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)) @@ -546,12 +549,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 !************************************/ @@ -613,6 +618,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 @@ -623,6 +630,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 @@ -653,6 +662,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 diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 697e625461..c53da614cf 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -128,7 +128,9 @@ module EDCLMLinkMod ! 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 :: 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 :: 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 + contains @@ -250,6 +252,7 @@ subroutine InitAllocate(this, bounds) allocate(this%nep_col (begc:endc)) ; this%nep_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%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 @@ -388,10 +391,6 @@ 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', & ptr_patch=this%sum_fuel_patch, set_lake=0._r8, set_urb=0._r8) @@ -493,6 +492,11 @@ subroutine InitHistory(this, bounds) 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', & @@ -2246,14 +2250,17 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil hr => soilbiogeochem_carbonflux_inst%hr_col, & ! Output: (gC/m2/s) total heterotrophic respiration npp_hifreq => this%npp_hifreq_col, & nep => this%nep_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & nbp => this%nbp_col & ) ! set time steps dt = real( get_step_size(), r8 ) + ! zero npp first do c = bounds%begc,bounds%endc npp_hifreq(c) = 0._r8 + fire_c_to_atm(c) = 0._r8 end do ! retrieve the first soil patch associated with each gridcell. @@ -2269,6 +2276,11 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil do g = bounds%begg,bounds%endg if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + ! first map ed site-level fire fluxes to clm column fluxes + cc = ed_allsites_inst(g)%clmcolumn + fire_c_to_atm(cc) = ed_allsites_inst(g)%total_burn_flux_to_atm / ( AREA * SHR_CONST_CDAY * 1.e3_r8) + + ! second map ed cohort-level npp fluxes to clm column fluxes currentPatch => ed_allsites_inst(g)%oldest_patch do while(associated(currentPatch)) cs => currentPatch%siteptr @@ -2284,12 +2296,10 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil end do ! calculate NEP and NBP fluxes. - !!!! CDK FEB/26/2016: THIS IS IMPORTANT AND NEEDS TO CHANGE AS IT IS ONLY A PLACEHOLDER. - !!!! NEP AND NBP ARE BOTH THE SAME RIGHT NOW BECAUSE I DON'T KNOW HOW TO ADD IN THE FIRE, DISTURBANCE, ETC FLUXES INTO THE NBP FLUX YET do fc = 1,num_soilc c = filter_soilc(fc) nep(c) = npp_hifreq(c) - hr(c) - nbp(c) = npp_hifreq(c) - hr(c) + nbp(c) = npp_hifreq(c) - ( hr(c) + fire_c_to_atm(c) ) end do end associate diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index ab816a9d7c..c4cb994dc5 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -413,13 +413,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 From f9571fd4421e861a0352e4478efdfbd30dff4023 Mon Sep 17 00:00:00 2001 From: ckoven Date: Sun, 28 Feb 2016 10:00:45 -0800 Subject: [PATCH 11/22] bugfixes on fire NEP/NBP code --- main/EDCLMLinkMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index c53da614cf..94cb3ef1b4 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -2228,6 +2228,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil use ColumnType , only : col use LandunitType , only : lun use landunit_varcon , only : istsoil + use shr_const_mod, only: SHR_CONST_CDAY class(ed_clm_type) :: this type(bounds_type) , intent(in) :: bounds From 3430b2af5b1ad4dec799fded5a831701393dfa83 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Mon, 29 Feb 2016 11:56:37 -0800 Subject: [PATCH 12/22] added combined ED/BGC C stock balance check versus daily-integrated NBP infrastructure --- biogeochem/EDBGCDynMod.F90 | 19 +++- main/EDCLMLinkMod.F90 | 217 +++++++++++++++++++++++++++++++++---- 2 files changed, 212 insertions(+), 24 deletions(-) diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 index 9e23408432..05443a48c7 100644 --- a/biogeochem/EDBGCDynMod.F90 +++ b/biogeochem/EDBGCDynMod.F90 @@ -345,16 +345,27 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so ! call soilbiogeochem_nitrogenflux_inst%Summary(bounds, num_soilc, filter_soilc) ! ---------------------------------------------- - ! ed veg carbon flux summary + ! ed veg carbon state and flux summary ! ---------------------------------------------- - call ed_clm_inst%Summary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & - ed_allsites_inst(bounds%begg:bounds%endg), soilbiogeochem_carbonflux_inst) + call ed_clm_inst%Summary(bounds, num_soilc, filter_soilc, & + ed_allsites_inst(bounds%begg:bounds%endg), & + soilbiogeochem_carbonflux_inst, & + soilbiogeochem_carbonstate_inst) ! ---------------------------------------------- - ! ed veg carbon/nitrogen flux summary + ! ed veg nitrogen flux summary ! ---------------------------------------------- + ! TBD... + + ! ---------------------------------------------- + ! calculate balance checks on entire carbon cycle (ED + BGC) + ! ---------------------------------------------- + + call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc, & + ed_allsites_inst(bounds%begg:bounds%endg)) + call t_stopf('BGCsum') end subroutine EDBGCDynSummary diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 94cb3ef1b4..3abd5ef311 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -17,6 +17,7 @@ module EDCLMLinkMod 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 : get_step_size ! @@ -127,10 +128,18 @@ module EDCLMLinkMod ! 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 :: 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 + ! 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 :: 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 contains @@ -140,6 +149,8 @@ module EDCLMLinkMod procedure , public :: SetValues procedure , public :: ed_clm_link procedure , public :: Summary + procedure , public :: ED_BGC_Carbon_Balancecheck + procedure , public :: Summary ! Private routines procedure , private :: ed_clm_leaf_area_profile @@ -250,10 +261,18 @@ subroutine InitAllocate(this, bounds) allocate(this%stem_prof_col (begc:endc,1:nlevdecomp_full)) ; this%stem_prof_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%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%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%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 @@ -336,15 +355,15 @@ subroutine InitHistory(this, bounds) 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', & + 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', & + 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', & @@ -507,6 +526,31 @@ subroutine InitHistory(this, bounds) 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%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) 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', & @@ -669,6 +713,16 @@ subroutine Restart ( this, bounds, ncid, flag ) ! 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%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) + 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, & @@ -2218,8 +2272,9 @@ end subroutine flux_into_litter_pools !------------------------------------------------------------------------ - subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & - ed_allsites_inst, soilbiogeochem_carbonflux_inst) + 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 @@ -2229,16 +2284,19 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil use LandunitType , only : lun use landunit_varcon , only : istsoil use shr_const_mod, only: SHR_CONST_CDAY - - class(ed_clm_type) :: this + ! + 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 - integer , intent(in) :: num_soilp ! number of soil patches in filter - integer , intent(in) :: filter_soilp(:) ! filter for soil patches 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 @@ -2248,20 +2306,34 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil 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, & ! Output: (gC/m2/s) total heterotrophic respiration + 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 & + nbp => this%nbp_col, & + totecosysc => this%totecosysc_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 ) ! set time steps dt = real( get_step_size(), r8 ) - ! zero npp first + ! 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. @@ -2277,18 +2349,31 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil do g = bounds%begg,bounds%endg if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then - ! first map ed site-level fire fluxes to clm column fluxes 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) - - ! second map ed cohort-level npp fluxes to clm column fluxes + currentPatch => ed_allsites_inst(g)%oldest_patch do while(associated(currentPatch)) - cs => currentPatch%siteptr - cc = cs%clmcolumn + + ! 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) + ed_litter_stock(cc) = ed_litter_stock(cc) + (currentPatch%area / AREA) * & + (sum(currentPatch%leaf_litter)+sum(currentPatch%root_litter)) + seed_stock(cc) = seed_stock(cc) + (currentPatch%area / AREA) * sum(currentPatch%seed_bank) + 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 / AREA + currentCohort => currentCohort%shorter enddo !currentCohort currentPatch => currentPatch%younger @@ -2303,9 +2388,101 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, num_soilp, filter_soil 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) + + totecosysc(c) = totsomc(c) + totlitc(c) + & ! BGC stocks + ed_litter_stock(c) + cwd_stock(c) + seed_stock(c) + biomass_stock(c) ! ED stocks + end do + end associate end subroutine Summary - + + + subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, & + ed_allsites_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: + use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep + ! + 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 + ! + ! !LOCAL VARIABLES: + real(r8) :: dtime ! land model time step (sec) + real(r8) :: nstep ! model timestep + real(r8) :: nbp_integrated(bounds%begc:bounds%endc) ! total net biome production integrated + real(r8) :: error_tolerance = 1.e-6_r8 + + associate(& + nep => this%nep_col, & + nep_timeintegrated => this%nep_timeintegrated_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & + totecosysc_old => this%totecosysc_old_col, & + totecosysc => this%totecosysc_col & + ) + + dtime = get_step_size() + nstep = get_nstep() + + if (nstep .eq. 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) + nep_timeintegrated(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 + 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 + nbp_integrated(c) = nep_timeintegrated(c) + fire_c_to_atm(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(c) = totecosysc_old(c) - totecosysc(c) - nbp_integrated(c) + 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 + ! + write(iulog,*) 'ED_BGC_Carbon_Balancecheck: max carbon balance error (gC / m2 / day): ', max(abs(error(:))) + + ! reset the C stock and flux integrators + do fc = 1,num_soilc + c = filter_soilc(fc) + totecosysc_old(c) = totecosysc(c) + nep_timeintegrated(c) = 0._r8 + end do + + endif + + end associate + + end subroutine ED_BGC_Carbon_Balancecheck end module EDCLMLinkMod From c12e399bcb8c3f9306855374218bf11d3afd321b Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 29 Feb 2016 15:05:40 -0800 Subject: [PATCH 13/22] compile bugfixes. still crashing in some tests though --- biogeochem/EDBGCDynMod.F90 | 3 +-- main/EDCLMLinkMod.F90 | 22 ++++++++++++++++------ 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 index 05443a48c7..9a59030519 100644 --- a/biogeochem/EDBGCDynMod.F90 +++ b/biogeochem/EDBGCDynMod.F90 @@ -363,8 +363,7 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so ! calculate balance checks on entire carbon cycle (ED + BGC) ! ---------------------------------------------- - call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc, & - ed_allsites_inst(bounds%begg:bounds%endg)) + call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc) call t_stopf('BGCsum') diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 3abd5ef311..3de0c8165e 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -150,7 +150,6 @@ module EDCLMLinkMod procedure , public :: ed_clm_link procedure , public :: Summary procedure , public :: ED_BGC_Carbon_Balancecheck - procedure , public :: Summary ! Private routines procedure , private :: ed_clm_leaf_area_profile @@ -1956,12 +1955,13 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch use EDTypesMod, only : AREA, numpft_ed use SoilBiogeochemVerticalProfileMod, only: exponential_rooting_profile, pftspecific_rootingprofile, rootprof_exp, surfprof_exp use pftconMod, only : pftcon - use shr_const_mod, only: SHR_CONST_CDAY + 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 + use shr_const_mod, only: SHR_CONST_CDAY ! implicit none ! @@ -2359,7 +2359,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ! 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) + sum(currentPatch%cwd_bg)) ed_litter_stock(cc) = ed_litter_stock(cc) + (currentPatch%area / AREA) * & (sum(currentPatch%leaf_litter)+sum(currentPatch%root_litter)) seed_stock(cc) = seed_stock(cc) + (currentPatch%area / AREA) * sum(currentPatch%seed_bank) @@ -2401,8 +2401,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & end subroutine Summary - subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, & - ed_allsites_inst) + subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) ! 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 @@ -2411,6 +2410,7 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, & ! ! !USES: use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep + use shr_const_mod, only: SHR_CONST_CDAY ! implicit none ! @@ -2424,7 +2424,10 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, & real(r8) :: dtime ! land model time step (sec) real(r8) :: nstep ! model timestep real(r8) :: nbp_integrated(bounds%begc:bounds%endc) ! total net biome production integrated + real(r8) :: error(bounds%begc:bounds%endc) real(r8) :: error_tolerance = 1.e-6_r8 + real(r8) :: max_error + integer :: fc,c associate(& nep => this%nep_col, & @@ -2470,7 +2473,14 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, & ! ! RETURN TO THIS LATER AND ADD A CRASHER IF BALANCE EXCEEDS THRESHOLD ! - write(iulog,*) 'ED_BGC_Carbon_Balancecheck: max carbon balance error (gC / m2 / day): ', max(abs(error(:))) + max_error = 0._r8 + do fc = 1,num_soilc + c = filter_soilc(fc) + if (abs(error(c)) .gt. max_error) then + max_error = abs(error(c)) + endif + end do + write(iulog,*) 'ED_BGC_Carbon_Balancecheck: max carbon balance error (gC / m2 / day): ', max_error ! reset the C stock and flux integrators do fc = 1,num_soilc From dc2bb37f87e9be07053808777aa4c5e2f4348ced Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 29 Feb 2016 16:48:29 -0800 Subject: [PATCH 14/22] runtime bugfix --- main/EDCLMLinkMod.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 3de0c8165e..4fc885c3eb 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -2422,7 +2422,7 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) ! ! !LOCAL VARIABLES: real(r8) :: dtime ! land model time step (sec) - real(r8) :: nstep ! model timestep + integer :: nstep ! model timestep real(r8) :: nbp_integrated(bounds%begc:bounds%endc) ! total net biome production integrated real(r8) :: error(bounds%begc:bounds%endc) real(r8) :: error_tolerance = 1.e-6_r8 @@ -2440,13 +2440,13 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) dtime = get_step_size() nstep = get_nstep() - if (nstep .eq. 1) then + 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) nep_timeintegrated(c) = 0._r8 - end do + end do endif if ( .not. is_beg_curr_day() ) then From 3054567f597ca3a396f2b3490d6cebc88b401696 Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 2 Mar 2016 15:15:44 -0800 Subject: [PATCH 15/22] runtime bugfix to prevent history file crashes on bgc profiles --- main/EDCLMLinkMod.F90 | 22 +++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 4fc885c3eb..e3c8f64cdc 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -551,26 +551,32 @@ subroutine InitHistory(this, bounds) 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') @@ -2094,7 +2100,7 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch end do else ! if fully frozen, or no roots, put everything in the top layer - froot_prof(c,ft,1) = 1./dzsoi_decomp(1) + froot_prof(c,ft,1) = 1._r8/dzsoi_decomp(1) endif end do ! @@ -2109,8 +2115,12 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch end do else ! if fully frozen, or no roots, put everything in the top layer - leaf_prof(c,1) = 1./dzsoi_decomp(1) - stem_prof(c,1) = 1./dzsoi_decomp(1) + 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 @@ -2266,6 +2276,12 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch ! 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 From c8cfa2060bca31a559fef48f92df88639222c207 Mon Sep 17 00:00:00 2001 From: ckoven Date: Fri, 4 Mar 2016 17:08:05 -0800 Subject: [PATCH 16/22] regularized the units of many of the history outputs to be consistent with the rest of CLM --- biogeochem/EDCohortDynamicsMod.F90 | 2 +- main/EDCLMLinkMod.F90 | 196 ++++++++++++++--------------- 2 files changed, 98 insertions(+), 100 deletions(-) diff --git a/biogeochem/EDCohortDynamicsMod.F90 b/biogeochem/EDCohortDynamicsMod.F90 index d921e42533..111a4a4d17 100755 --- a/biogeochem/EDCohortDynamicsMod.F90 +++ b/biogeochem/EDCohortDynamicsMod.F90 @@ -723,7 +723,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/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index e3c8f64cdc..4d44d8dee5 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -19,6 +19,7 @@ module EDCLMLinkMod use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type use SoilBiogeochemCarbonStatetype , only : soilbiogeochem_carbonstate_type use clm_time_manager , only : get_step_size + use shr_const_mod, only: SHR_CONST_CDAY ! implicit none @@ -76,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 @@ -156,7 +157,7 @@ module EDCLMLinkMod 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 @@ -179,7 +180,7 @@ subroutine Init(this, bounds) call this%InitAllocate(bounds) call this%InitHistory(bounds) - call this%InitCold(bounds) + !call this%InitCold(bounds) end subroutine Init @@ -240,12 +241,12 @@ 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 @@ -349,15 +350,15 @@ 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', & + 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', & + 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) @@ -409,55 +410,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='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) @@ -465,35 +466,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', & @@ -656,26 +657,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 ) ! @@ -1342,16 +1342,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 @@ -1426,13 +1426,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 @@ -1967,7 +1967,6 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch use abortutils , only : endrun use shr_log_mod , only : errMsg => shr_log_errMsg use EDParamsMod, only : ED_val_ag_biomass - use shr_const_mod, only: SHR_CONST_CDAY ! implicit none ! @@ -2299,7 +2298,6 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & use ColumnType , only : col use LandunitType , only : lun use landunit_varcon , only : istsoil - use shr_const_mod, only: SHR_CONST_CDAY ! implicit none ! @@ -2375,10 +2373,10 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ! 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)) + 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)) - seed_stock(cc) = seed_stock(cc) + (currentPatch%area / AREA) * sum(currentPatch%seed_bank) + (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)) @@ -2388,7 +2386,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ! map biomass pools to column level biomass_stock(cc) = biomass_stock(cc) + (currentCohort%bdead + currentCohort%balive + & - currentCohort%bstore) * currentCohort%n / AREA + currentCohort%bstore) * currentCohort%n * 1.e3_r8 / AREA currentCohort => currentCohort%shorter enddo !currentCohort From 57c6fca3ac45ae737a98039711038dee175ce634 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Mon, 7 Mar 2016 11:10:02 -0800 Subject: [PATCH 17/22] added term to take into acocunt time offset of reconciling ED and CLM litter fluxes for balance check purposes --- main/EDCLMLinkMod.F90 | 65 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 62 insertions(+), 3 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 4d44d8dee5..83d994ded3 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -133,6 +133,8 @@ module EDCLMLinkMod 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 ! 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 @@ -260,6 +262,9 @@ subroutine InitAllocate(this, bounds) 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%nep_col (begc:endc)) ; this%nep_col (:) = nan allocate(this%nep_timeintegrated_col (begc:endc)) ; this%nep_timeintegrated_col (:) = nan allocate(this%nbp_col (begc:endc)) ; this%nbp_col (:) = nan @@ -728,6 +733,16 @@ subroutine Restart ( this, bounds, ncid, flag ) 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) + 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, & @@ -2331,7 +2346,9 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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 + 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 & ) ! set time steps @@ -2410,6 +2427,36 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ed_litter_stock(c) + cwd_stock(c) + seed_stock(c) + biomass_stock(c) ! ED stocks 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 fc = 1,num_soilc + c = filter_soilc(fc) + ed_to_bgc_last_edts(c) = ed_to_bgc_this_edts(c) + end do + ! + do fc = 1,num_soilc + c = filter_soilc(fc) + ed_to_bgc_this_edts(c) = 0._r8 + 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(c) = ed_to_bgc_this_edts(c) + (currentpatch%CWD_AG_out(c) + currentpatch%CWD_BG_out(c) + & + sum(currentpatch%leaf_litter_out(:)) + sum(currentpatch%root_litter_out(:))) & + * currentpatch%area/AREA * 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 @@ -2448,7 +2495,9 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) nep_timeintegrated => this%nep_timeintegrated_col, & fire_c_to_atm => this%fire_c_to_atm_col, & totecosysc_old => this%totecosysc_old_col, & - totecosysc => this%totecosysc_col & + totecosysc => this%totecosysc_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 & ) dtime = get_step_size() @@ -2460,6 +2509,10 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) c = filter_soilc(fc) totecosysc_old(c) = totecosysc(c) nep_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 end do endif @@ -2477,10 +2530,16 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) nbp_integrated(c) = nep_timeintegrated(c) + fire_c_to_atm(c) * SHR_CONST_CDAY end do + ! adjust the nbp to take into account the fact that ED litter stocks reflect today's changes while BGC reflects yesterday's + do fc = 1,num_soilc + c = filter_soilc(fc) + nbp_integrated(c) = nbp_integrated(c) - (ed_to_bgc_this_edts(c) - ed_to_bgc_last_edts(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(c) = totecosysc_old(c) - totecosysc(c) - nbp_integrated(c) + error(c) = totecosysc(c) - totecosysc_old(c) - nbp_integrated(c) end do ! for now, rather than crashing the model, lets just report the largest error to see what we're up against From b772df7363f7b5a442430dda8b7a56c285428145 Mon Sep 17 00:00:00 2001 From: ckoven Date: Mon, 7 Mar 2016 17:30:07 -0800 Subject: [PATCH 18/22] still chasing carbon balance errors... --- biogeochem/EDBGCDynMod.F90 | 3 +- main/EDCLMLinkMod.F90 | 108 +++++++++++++++++++++++++++++-------- 2 files changed, 88 insertions(+), 23 deletions(-) diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 index 9a59030519..fe8acf0959 100644 --- a/biogeochem/EDBGCDynMod.F90 +++ b/biogeochem/EDBGCDynMod.F90 @@ -363,7 +363,8 @@ subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_so ! calculate balance checks on entire carbon cycle (ED + BGC) ! ---------------------------------------------- - call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc) + call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc, & + soilbiogeochem_carbonflux_inst) call t_stopf('BGCsum') diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 83d994ded3..3df957e3b9 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -18,7 +18,7 @@ module EDCLMLinkMod use EDParamsMod , only : ED_val_ag_biomass use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type use SoilBiogeochemCarbonStatetype , only : soilbiogeochem_carbonstate_type - use clm_time_manager , only : get_step_size + use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep use shr_const_mod, only: SHR_CONST_CDAY ! @@ -130,6 +130,8 @@ module EDCLMLinkMod ! 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 @@ -139,6 +141,10 @@ module EDCLMLinkMod ! 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 @@ -267,12 +273,19 @@ subroutine InitAllocate(this, bounds) 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 @@ -728,11 +741,31 @@ subroutine Restart ( this, bounds, ncid, flag ) 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%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='', & @@ -2269,6 +2302,12 @@ subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch 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 @@ -2343,6 +2382,8 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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 @@ -2422,9 +2463,11 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ! 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) - totecosysc(c) = totsomc(c) + totlitc(c) + & ! BGC stocks - ed_litter_stock(c) + cwd_stock(c) + seed_stock(c) + biomass_stock(c) ! ED stocks end do ! in ED timesteps, because of offset between when ED and BGC reconcile the gain and loss of litterfall carbon, @@ -2447,8 +2490,8 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & currentPatch => ed_allsites_inst(g)%oldest_patch do while(associated(currentPatch)) ! - ed_to_bgc_this_edts(c) = ed_to_bgc_this_edts(c) + (currentpatch%CWD_AG_out(c) + currentpatch%CWD_BG_out(c) + & - sum(currentpatch%leaf_litter_out(:)) + sum(currentpatch%root_litter_out(:))) & + 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) ! currentPatch => currentPatch%younger @@ -2462,7 +2505,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & end subroutine Summary - subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) + 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 @@ -2470,8 +2513,6 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) ! Written by Charlie Koven, Feb 2016 ! ! !USES: - use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep - use shr_const_mod, only: SHR_CONST_CDAY ! implicit none ! @@ -2480,22 +2521,35 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) 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(bounds%begc:bounds%endc) + 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 + 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 & ) @@ -2508,7 +2562,11 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) 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 @@ -2521,45 +2579,51 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc) 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 - nbp_integrated(c) = nep_timeintegrated(c) + fire_c_to_atm(c) * SHR_CONST_CDAY - end do - - ! adjust the nbp to take into account the fact that ED litter stocks reflect today's changes while BGC reflects yesterday's - do fc = 1,num_soilc - c = filter_soilc(fc) - nbp_integrated(c) = nbp_integrated(c) - (ed_to_bgc_this_edts(c) - ed_to_bgc_last_edts(c)) * SHR_CONST_CDAY + 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 end do ! next compare the change in carbon and calculate the error do fc = 1,num_soilc c = filter_soilc(fc) - error(c) = totecosysc(c) - totecosysc_old(c) - nbp_integrated(c) + error_ed(c) = totedc(c) - totedc_old(c) - (npp_timeintegrated(c) - 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 ! 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 = 0._r8 + max_error_total = 0._r8 do fc = 1,num_soilc c = filter_soilc(fc) - if (abs(error(c)) .gt. max_error) then - max_error = abs(error(c)) + 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 carbon balance error (gC / m2 / day): ', max_error + 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 From 71ebb915b4e1b28a6d5097e20cb8e74f75b7a689 Mon Sep 17 00:00:00 2001 From: ckoven Date: Tue, 8 Mar 2016 13:33:27 -0800 Subject: [PATCH 19/22] runtime bugfix to prevent crash due to filter mismatch; though filter mismatch is likely still happening. --- main/EDCLMLinkMod.F90 | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 3df957e3b9..32366e9ae4 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -2474,14 +2474,19 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & ! (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 fc = 1,num_soilc - c = filter_soilc(fc) - ed_to_bgc_last_edts(c) = ed_to_bgc_this_edts(c) + ! + 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 fc = 1,num_soilc - c = filter_soilc(fc) - ed_to_bgc_this_edts(c) = 0._r8 + 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 + endif end do ! do g = bounds%begg,bounds%endg @@ -2490,9 +2495,9 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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) + 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 ) ! currentPatch => currentPatch%younger end do !currentPatch From 00a104d49c7d841e0c96abec8e48b5c89fbb0977 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Wed, 9 Mar 2016 15:57:39 -0800 Subject: [PATCH 20/22] added seed rain flux to carbon balance logic --- biogeochem/EDPhysiologyMod.F90 | 3 +++ main/EDCLMLinkMod.F90 | 20 ++++++++++++++++---- main/EDTypesMod.F90 | 1 + 3 files changed, 20 insertions(+), 4 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index f8e3464fb9..f773adb199 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -634,6 +634,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 @@ -647,6 +649,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 diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 32366e9ae4..4908cdb75c 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -137,6 +137,7 @@ module EDCLMLinkMod 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 @@ -270,6 +271,7 @@ subroutine InitAllocate(this, bounds) 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 @@ -775,6 +777,11 @@ subroutine Restart ( this, bounds, ncid, flag ) 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 @@ -2389,7 +2396,8 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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 & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col ) ! set time steps @@ -2486,6 +2494,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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 ! @@ -2499,6 +2508,8 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & + 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 @@ -2556,7 +2567,8 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soi 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 & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col & ) dtime = get_step_size() @@ -2594,13 +2606,13 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soi 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 + 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) - ed_to_bgc_this_edts(c)* SHR_CONST_CDAY - fire_c_to_atm(c) * SHR_CONST_CDAY) + 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 diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index c4cb994dc5..b06b26c670 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -290,6 +290,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 From de3fee30d5a84b75e64fc23ec1c734458d72992e Mon Sep 17 00:00:00 2001 From: ckoven Date: Wed, 9 Mar 2016 17:20:17 -0800 Subject: [PATCH 21/22] sending C balance errors to history file to keep track of --- main/EDCLMLinkMod.F90 | 80 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 67 insertions(+), 13 deletions(-) diff --git a/main/EDCLMLinkMod.F90 b/main/EDCLMLinkMod.F90 index 4908cdb75c..53b79e1f8a 100755 --- a/main/EDCLMLinkMod.F90 +++ b/main/EDCLMLinkMod.F90 @@ -151,6 +151,11 @@ module EDCLMLinkMod 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 @@ -291,7 +296,11 @@ subroutine InitAllocate(this, bounds) 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%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 @@ -551,6 +560,21 @@ subroutine InitHistory(this, bounds) 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', & @@ -758,6 +782,21 @@ subroutine Restart ( this, bounds, ncid, flag ) 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='', & @@ -2397,7 +2436,7 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc, & 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 + seed_rain_flux => this%seed_rain_flux_col & ) ! set time steps @@ -2568,7 +2607,10 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soi 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 & + 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() @@ -2588,6 +2630,10 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soi ! 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 @@ -2616,21 +2662,29 @@ subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soi 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 + ! 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 From eba7a1791de1f0f5ba062e3f00f55033a4bbad19 Mon Sep 17 00:00:00 2001 From: Charlie Koven Date: Tue, 22 Mar 2016 11:22:53 -0700 Subject: [PATCH 22/22] changed inline description in EDBGCDynMod and moved it outside the ED subdirectory hierarchy since it properly sits outside the area defined by the interface --- biogeochem/EDBGCDynMod.F90 | 373 ------------------------------------- 1 file changed, 373 deletions(-) delete mode 100644 biogeochem/EDBGCDynMod.F90 diff --git a/biogeochem/EDBGCDynMod.F90 b/biogeochem/EDBGCDynMod.F90 deleted file mode 100644 index fe8acf0959..0000000000 --- a/biogeochem/EDBGCDynMod.F90 +++ /dev/null @@ -1,373 +0,0 @@ -module EDBGCDynMod - -! Interface from ED calls to CLM belowground biogeochemistry module - - use shr_kind_mod , only : r8 => shr_kind_r8 - use clm_varctl , only : use_c13, use_c14, use_ed - use decompMod , only : bounds_type - use perf_mod , only : t_startf, t_stopf - use clm_varctl , only : use_century_decomp, use_nitrif_denitrif - use CNVegStateType , only : cnveg_state_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type - use CNVegCarbonFluxType , only : cnveg_carbonflux_type - use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type - use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type - use SoilBiogeochemStateType , only : soilbiogeochem_state_type - use SoilBiogeochemCarbonStateType , only : soilbiogeochem_carbonstate_type - use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type - use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type - use SoilBiogeochemNitrogenFluxType , only : soilbiogeochem_nitrogenflux_type - use EDCLMLinkMod , only : ed_clm_type - use CanopyStateType , only : canopystate_type - use SoilStateType , only : soilstate_type - use SoilHydrologyType , only : soilhydrology_type - use TemperatureType , only : temperature_type - use WaterstateType , only : waterstate_type - use WaterfluxType , only : waterflux_type - use atm2lndType , only : atm2lnd_type - use SoilStateType , only : soilstate_type - use ch4Mod , only : ch4_type - use EDtypesMod , only : ed_site_type - - ! public :: EDBGCDynInit ! BGC dynamics: initialization - public :: EDBGCDyn ! BGC Dynamics - public :: EDBGCDynSummary ! BGC dynamics: summary - -contains - - - !----------------------------------------------------------------------- - subroutine EDBGCDyn(bounds, & - num_soilc, filter_soilc, num_soilp, filter_soilp, num_pcropp, filter_pcropp, doalb, & - cnveg_state_inst, & - cnveg_carbonflux_inst, cnveg_carbonstate_inst, & - ed_clm_inst, & - soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & - soilbiogeochem_state_inst, & - soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, & - c13_soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonflux_inst, & - c14_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonflux_inst, & - atm2lnd_inst, waterstate_inst, waterflux_inst, & - canopystate_inst, soilstate_inst, temperature_inst, crop_inst, ch4_inst) - ! - ! !DESCRIPTION: - - ! - ! !USES: - use clm_varpar , only: crop_prog, nlevgrnd, nlevdecomp_full - use clm_varpar , only: nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools - use subgridAveMod , only: p2c - use CropType , only: crop_type - use CNNDynamicsMod , only: CNNDeposition,CNNFixation, CNNFert, CNSoyfix - use CNMRespMod , only: CNMResp - use CNPhenologyMod , only: CNPhenology - use CNGRespMod , only: CNGResp - use CNFireMod , only: CNFireArea, CNFireFluxes - use CNCIsoFluxMod , only: CIsoFlux1, CIsoFlux2, CIsoFlux2h, CIsoFlux3 - use CNC14DecayMod , only: C14Decay - use CNWoodProductsMod , only: CNWoodProducts - use CNCStateUpdate1Mod , only: CStateUpdate1,CStateUpdate0 - use CNCStateUpdate2Mod , only: CStateUpdate2, CStateUpdate2h - use CNCStateUpdate3Mod , only: CStateUpdate3 - use CNNStateUpdate1Mod , only: NStateUpdate1 - use CNNStateUpdate2Mod , only: NStateUpdate2, NStateUpdate2h - use CNGapMortalityMod , only: CNGapMortality - use dynHarvestMod , only: CNHarvest - use SoilBiogeochemDecompCascadeBGCMod , only: decomp_rate_constants_bgc - use SoilBiogeochemDecompCascadeCNMod , only: decomp_rate_constants_cn - use SoilBiogeochemCompetitionMod , only: SoilBiogeochemCompetition - use SoilBiogeochemDecompMod , only: SoilBiogeochemDecomp - use SoilBiogeochemLittVertTranspMod , only: SoilBiogeochemLittVertTransp - use SoilBiogeochemPotentialMod , only: SoilBiogeochemPotential - use SoilBiogeochemVerticalProfileMod , only: SoilBiogeochemVerticalProfile - use SoilBiogeochemNitrifDenitrifMod , only: SoilBiogeochemNitrifDenitrif - use SoilBiogeochemNStateUpdate1Mod , only: SoilBiogeochemNStateUpdate1 - ! - ! !ARGUMENTS: - 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 - integer , intent(in) :: num_soilp ! number of soil patches in filter - integer , intent(in) :: filter_soilp(:) ! filter for soil patches - integer , intent(in) :: num_pcropp ! number of prog. crop patches in filter - integer , intent(in) :: filter_pcropp(:) ! filter for prognostic crop patches - logical , intent(in) :: doalb ! true = surface albedo calculation time step - type(cnveg_state_type) , intent(inout) :: cnveg_state_inst - type(cnveg_carbonflux_type) , intent(inout) :: cnveg_carbonflux_inst - type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst - type(ed_clm_type) , intent(inout) :: ed_clm_inst - type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst - type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst - type(soilbiogeochem_carbonflux_type) , intent(inout) :: c13_soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: c13_soilbiogeochem_carbonstate_inst - type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst - type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst - type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst - type(atm2lnd_type) , intent(in) :: atm2lnd_inst - type(waterstate_type) , intent(in) :: waterstate_inst - type(waterflux_type) , intent(in) :: waterflux_inst - type(canopystate_type) , intent(in) :: canopystate_inst - type(soilstate_type) , intent(in) :: soilstate_inst - type(temperature_type) , intent(inout) :: temperature_inst - type(crop_type) , intent(in) :: crop_inst - type(ch4_type) , intent(in) :: ch4_inst - ! - ! !LOCAL VARIABLES: - real(r8):: cn_decomp_pools(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools) - real(r8):: p_decomp_cpool_loss(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential C loss from one pool to another - real(r8):: pmnf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) !potential mineral N flux, from one pool to another - real(r8):: arepr(bounds%begp:bounds%endp) ! reproduction allocation coefficient (only used for crop_prog) - real(r8):: aroot(bounds%begp:bounds%endp) ! root allocation coefficient (only used for crop_prog) - integer :: begp,endp - integer :: begc,endc - !----------------------------------------------------------------------- - - begp = bounds%begp; endp = bounds%endp - begc = bounds%begc; endc = bounds%endc - - !real(r8) , intent(in) :: rootfr_patch(bounds%begp:, 1:) - !integer , intent(in) :: altmax_lastyear_indx_col(bounds%begc:) ! frost table depth (m) - - associate( & - rootfr_patch => soilstate_inst%rootfr_patch , & ! fraction of roots in each soil layer (nlevgrnd) - altmax_lastyear_indx_col => canopystate_inst%altmax_lastyear_indx_col , & ! frost table depth (m) - laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit projected leaf area index - laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded projected leaf area index - frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] - frac_veg_nosno_alb => canopystate_inst%frac_veg_nosno_alb_patch , & ! Output: [integer (:) ] frac of vegetation not covered by snow [-] - tlai => canopystate_inst%tlai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index, no burying by snow - tsai => canopystate_inst%tsai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index, no burying by snow - elai => canopystate_inst%elai_patch , & ! Output: [real(r8) (:) ] one-sided leaf area index with burying by snow - esai => canopystate_inst%esai_patch , & ! Output: [real(r8) (:) ] one-sided stem area index with burying by snow - htop => canopystate_inst%htop_patch , & ! Output: [real(r8) (:) ] canopy top (m) - hbot => canopystate_inst%hbot_patch & ! Output: [real(r8) (:) ] canopy bottom (m) - ) - - ! -------------------------------------------------- - ! zero the column-level C and N fluxes - ! -------------------------------------------------- - - call t_startf('BGCZero') - - call soilbiogeochem_carbonflux_inst%SetValues( & - num_soilc, filter_soilc, 0._r8) - if ( use_c13 ) then - call c13_soilbiogeochem_carbonflux_inst%SetValues( & - num_soilc, filter_soilc, 0._r8) - end if - if ( use_c14 ) then - call c14_soilbiogeochem_carbonflux_inst%SetValues( & - num_soilc, filter_soilc, 0._r8) - end if - - call t_stopf('BGCZero') - - ! -------------------------------------------------- - ! Nitrogen Deposition, Fixation and Respiration - ! -------------------------------------------------- - - ! call t_startf('CNDeposition') - ! call CNNDeposition(bounds, & - ! atm2lnd_inst, soilbiogeochem_nitrogenflux_inst) - ! call t_stopf('CNDeposition') - - - ! if (crop_prog) then - ! call CNNFert(bounds, num_soilc,filter_soilc, & - ! cnveg_nitrogenflux_inst, soilbiogeochem_nitrogenflux_inst) - - ! call CNSoyfix (bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & - ! waterstate_inst, crop_inst, cnveg_state_inst, cnveg_nitrogenflux_inst , & - ! soilbiogeochem_state_inst, soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst) - ! end if - - !-------------------------------------------- - ! Soil Biogeochemistry - !-------------------------------------------- - - if (use_century_decomp) then - call decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & - canopystate_inst, soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) - else - call decomp_rate_constants_cn(bounds, num_soilc, filter_soilc, & - canopystate_inst, soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) - end if - - ! calculate potential decomp rates and total immobilization demand (previously inlined in CNDecompAlloc) - call SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & - soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & - soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & - cn_decomp_pools=cn_decomp_pools(begc:endc,1:nlevdecomp,1:ndecomp_pools), & - p_decomp_cpool_loss=p_decomp_cpool_loss(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & - pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) - - - !-------------------------------------------- - ! Resolve the competition between plants and soil heterotrophs - ! for available soil mineral N resource - !-------------------------------------------- - ! will add this back in when integrtating hte nutirent cycles - - - !-------------------------------------------- - ! Calculate litter and soil decomposition rate - !-------------------------------------------- - - ! Calculation of actual immobilization and decomp rates, following - ! resolution of plant/heterotroph competition for mineral N (previously inlined in CNDecompAllocation in CNDecompMod) - - call t_startf('SoilBiogeochemDecomp') - - call SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, & - soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & - soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & - cn_decomp_pools=cn_decomp_pools(begc:endc,1:nlevdecomp,1:ndecomp_pools), & - p_decomp_cpool_loss=p_decomp_cpool_loss(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & - pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) - - call t_stopf('SoilBiogeochemDecomp') - - - !-------------------------------------------- - ! Update1 - !-------------------------------------------- - - call t_startf('BNGCUpdate1') - - - ! Update all prognostic carbon state variables (except for gap-phase mortality and fire fluxes) - call CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & - cnveg_state_inst, cnveg_carbonflux_inst, cnveg_carbonstate_inst, & - ed_clm_inst, & - soilbiogeochem_carbonflux_inst) - - call t_stopf('BNGCUpdate1') - - !-------------------------------------------- - ! Calculate vertical mixing of soil and litter pools - !-------------------------------------------- - - call t_startf('SoilBiogeochemLittVertTransp') - - call SoilBiogeochemLittVertTransp(bounds, num_soilc, filter_soilc, & - canopystate_inst, soilbiogeochem_state_inst, & - soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & - c13_soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonflux_inst, & - c14_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonflux_inst, & - soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst) - - call t_stopf('SoilBiogeochemLittVertTransp') - - end associate - - end subroutine EDBGCDyn - - - !----------------------------------------------------------------------- - subroutine EDBGCDynSummary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & - soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & - c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & - c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & - soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst, & - ed_clm_inst, ed_allsites_inst) - ! - ! !DESCRIPTION: - ! Call to all CN and SoilBiogeochem summary routines - ! also aggregate production and decomposition fluxes to whole-ecosystem balance fluxes - ! - ! !USES: - use clm_varpar , only: ndecomp_cascade_transitions - use CNPrecisionControlMod , only: CNPrecisionControl - use SoilBiogeochemPrecisionControlMod , only: SoilBiogeochemPrecisionControl - ! - ! !ARGUMENTS: - 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 - integer , intent(in) :: num_soilp ! number of soil patches in filter - integer , intent(in) :: filter_soilp(:) ! filter for soil patches - type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst - type(soilbiogeochem_carbonflux_type) , intent(inout) :: c13_soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: c13_soilbiogeochem_carbonstate_inst - type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst - type(soilbiogeochem_carbonstate_type) , intent(inout) :: c14_soilbiogeochem_carbonstate_inst - type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst - type(soilbiogeochem_nitrogenstate_type) , intent(inout) :: soilbiogeochem_nitrogenstate_inst - type(ed_clm_type) , intent(inout) :: ed_clm_inst - type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) - ! - ! !LOCAL VARIABLES: - integer :: begc,endc - !----------------------------------------------------------------------- - - begc = bounds%begc; endc= bounds%endc - - ! Call to all summary routines - - call t_startf('BGCsum') - - ! Set controls on very low values in critical state variables - - call SoilBiogeochemPrecisionControl(num_soilc, filter_soilc, & - soilbiogeochem_carbonstate_inst, c13_soilbiogeochem_carbonstate_inst, & - c14_soilbiogeochem_carbonstate_inst,soilbiogeochem_nitrogenstate_inst) - - ! Note - all summary updates to cnveg_carbonstate_inst and cnveg_carbonflux_inst are done in - ! soilbiogeochem_carbonstate_inst%summary and CNVeg_carbonstate_inst%summary - - ! ---------------------------------------------- - ! soilbiogeochem carbon/nitrogen state summary - ! ---------------------------------------------- - - call soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) - if ( use_c13 ) then - call c13_soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) - end if - if ( use_c14 ) then - call c14_soilbiogeochem_carbonstate_inst%summary(bounds, num_soilc, filter_soilc) - end if - ! call soilbiogeochem_nitrogenstate_inst%summary(bounds, num_soilc, filter_soilc) - - ! ---------------------------------------------- - ! soilbiogeochem carbon/nitrogen flux summary - ! ---------------------------------------------- - - call soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) - if ( use_c13 ) then - call c13_soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) - end if - if ( use_c14 ) then - call c14_soilbiogeochem_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc) - end if - ! call soilbiogeochem_nitrogenflux_inst%Summary(bounds, num_soilc, filter_soilc) - - ! ---------------------------------------------- - ! ed veg carbon state and flux summary - ! ---------------------------------------------- - - call ed_clm_inst%Summary(bounds, num_soilc, filter_soilc, & - ed_allsites_inst(bounds%begg:bounds%endg), & - soilbiogeochem_carbonflux_inst, & - soilbiogeochem_carbonstate_inst) - - ! ---------------------------------------------- - ! ed veg nitrogen flux summary - ! ---------------------------------------------- - - ! TBD... - - ! ---------------------------------------------- - ! calculate balance checks on entire carbon cycle (ED + BGC) - ! ---------------------------------------------- - - call ed_clm_inst%ED_BGC_Carbon_Balancecheck(bounds, num_soilc, filter_soilc, & - soilbiogeochem_carbonflux_inst) - - call t_stopf('BGCsum') - - end subroutine EDBGCDynSummary - -end module EDBGCDynMod