diff --git a/components/clm/bld/CLMBuildNamelist.pm b/components/clm/bld/CLMBuildNamelist.pm index 7e12d75746..f58e81cb5a 100755 --- a/components/clm/bld/CLMBuildNamelist.pm +++ b/components/clm/bld/CLMBuildNamelist.pm @@ -1105,9 +1105,9 @@ sub setup_cmdl_bgc_spinup { my @valid_values = $definition->get_valid_values( $var ); fatal_error("$var has an invalid value ($val). Valid values are: @valid_values\n"); } - if ( $nl_flags->{'bgc_spinup'} eq "on" && $nl_flags->{'use_cn'} ne ".true.") { - fatal_error("$var can not be '$nl_flags->{'bgc_spinup'}' if CN is turned off (use_cn=$nl_flags->{'use_cn'})."); - } + # if ( $nl_flags->{'bgc_spinup'} eq "on" && $nl_flags->{'use_cn'} ne ".true.") { + # fatal_error("$var can not be '$nl_flags->{'bgc_spinup'}' if CN is turned off (use_cn=$nl_flags->{'use_cn'})."); + # } if ( $nl->get_value("spinup_state") eq 0 && $nl_flags->{'bgc_spinup'} eq "on" ) { fatal_error("Namelist spinup_state contradicts the command line option bgc_spinup" ); } @@ -2399,9 +2399,14 @@ sub setup_logic_bgc_shared { if ( $physv->as_long() >= $physv->as_long("clm4_5")) { if ( $nl_flags->{'bgc_mode'} ne "sp" ) { - add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'decomp_depth_efolding', 'phys'=>$physv->as_string() ); add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'constrain_stress_deciduous_onset', 'phys'=>$physv->as_string() ); } + # FIXME(bja, 201606) the logic around ed / bgc_mode / + # use_century_decomp is confusing and messed up. This is a hack + # workaround. + if ($nl_flags->{'use_century_decomp'} == '.true.') { + add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'decomp_depth_efolding', 'phys'=>$physv->as_string() ); + } } } diff --git a/components/clm/cime_config/testdefs/testmods_dirs/clm/edTest/user_nl_clm b/components/clm/cime_config/testdefs/testmods_dirs/clm/edTest/user_nl_clm index 12415896e2..e32e6b3cb8 100644 --- a/components/clm/cime_config/testdefs/testmods_dirs/clm/edTest/user_nl_clm +++ b/components/clm/cime_config/testdefs/testmods_dirs/clm/edTest/user_nl_clm @@ -1,11 +1,12 @@ use_cn = .false. -use_century_decomp = .false. -use_vertsoilc = .false. +use_century_decomp = .true. +use_vertsoilc = .true. finidat = '' hist_mfilt = 365 hist_nhtfrq = -24 hist_empty_htapes = .true. hist_fincl1 = 'NPP','GPP','BTRAN','H2OSOI','TLAI','LITTER_IN','LITTER_OUT', - 'STORVEGC','FIRE_AREA','SCORCH_HEIGHT','FIRE_INTENSITY','FIRE_TFC_ROS','fire_fuel_mef', + 'FIRE_AREA','SCORCH_HEIGHT','FIRE_INTENSITY','FIRE_TFC_ROS','fire_fuel_mef', 'fire_fuel_bulkd','fire_fuel_sav','FIRE_NESTEROV_INDEX','PFTbiomass', - 'PFTleafbiomass','FIRE_ROS','WIND','TFC_ROS','DISPVEGC','AREA_TREES','AREA_PLANT' \ No newline at end of file + 'PFTleafbiomass','FIRE_ROS','WIND','AREA_TREES','AREA_PLANT', + 'TOTSOMC','TOTLITC','T_SCALAR','NEP','NBP','TOTECOSYSC' \ No newline at end of file diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 288c8c9dd9..5b383735ea 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -733,7 +733,7 @@ subroutine fuse_cohorts(patchptr) nextc%n*nextc%year_net_uptake(i))/newn endif enddo - + currentCohort%n = newn !remove fused cohort from the list nextc%taller%shorter => nextnextc diff --git a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 index 5bfdf1c71a..58c4656c06 100755 --- a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 @@ -200,8 +200,11 @@ subroutine spawn_patches( currentSite ) ! calculate area of disturbed land, in this timestep, by summing contributions from each existing patch. currentPatch => currentSite%youngest_patch + + ! zero site-level fire fluxes currentSite%cwd_ag_burned = 0.0_r8 currentSite%leaf_litter_burned = 0.0_r8 + currentSite%total_burn_flux_to_atm = 0.0_r8 site_areadis = 0.0_r8 do while(associated(currentPatch)) @@ -548,12 +551,14 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) burned_litter = new_patch%cwd_ag(c) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(c+1) !kG/m2/day new_patch%cwd_ag(c) = new_patch%cwd_ag(c) - burned_litter currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day enddo do p = 1,numpft_ed burned_litter = new_patch%leaf_litter(p) * patch_site_areadis/new_patch%area * currentPatch%burnt_frac_litter(dg_sf) new_patch%leaf_litter(p) = new_patch%leaf_litter(p) - burned_litter - currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/dat + currentSite%flux_out = currentSite%flux_out + burned_litter * new_patch%area !kG/site/day + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + burned_litter * new_patch%area !kG/site/day enddo !************************************/ @@ -615,6 +620,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) SF_val_CWD_frac(c) * bstem * currentCohort%cfa currentSite%flux_out = currentSite%flux_out + dead_tree_density * & AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + dead_tree_density * & + AREA * SF_val_CWD_frac(c) * bstem * currentCohort%cfa enddo @@ -625,6 +632,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) dead_tree_density * currentCohort%bl * currentCohort%cfa currentSite%flux_out = currentSite%flux_out + & dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm + & + dead_tree_density * AREA * currentCohort%bl * currentCohort%cfa enddo @@ -655,6 +664,8 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis) !KgC/gridcell/day currentSite%flux_out = currentSite%flux_out + burned_leaves * currentCohort%n * & patch_site_areadis/currentPatch%area * AREA + currentSite%total_burn_flux_to_atm = currentSite%total_burn_flux_to_atm+ burned_leaves * currentCohort%n * & + patch_site_areadis/currentPatch%area * AREA endif currentCohort%cfa = 0.0_r8 @@ -1113,7 +1124,7 @@ subroutine fuse_patches( csite ) if(nopatches > maxpatch)then iterate = 1 profiletol = profiletol * 1.1_r8 - write(iulog,*) 'maxpatch exceeded, triggering patch fusion iteration.',profiletol,nopatches + !---------------------------------------------------------------------! ! Making profile tolerance larger means that more fusion will happen ! !---------------------------------------------------------------------! diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index 8de0ecc9ad..846f8c11fc 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -125,15 +125,14 @@ subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, wa currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - currentPatch%root_litter_out(p) enddo - currentPatch%leaf_litter_in(:) = 0.0_r8 - currentPatch%root_litter_in(:) = 0.0_r8 - currentPatch%leaf_litter_out(:) = 0.0_r8 - currentPatch%root_litter_out(:) = 0.0_r8 - - currentPatch%CWD_AG_in(:) = 0.0_r8 - currentPatch%cwd_bg_in(:) = 0.0_r8 - currentPatch%CWD_AG_out(:) = 0.0_r8 - currentPatch%cwd_bg_out(:) = 0.0_r8 + ! currentPatch%leaf_litter_in(:) = 0.0_r8 + ! currentPatch%root_litter_in(:) = 0.0_r8 + ! currentPatch%leaf_litter_out(:) = 0.0_r8 + ! currentPatch%root_litter_out(:) = 0.0_r8 + ! currentPatch%CWD_AG_in(:) = 0.0_r8 + ! currentPatch%cwd_bg_in(:) = 0.0_r8 + ! currentPatch%CWD_AG_out(:) = 0.0_r8 + ! currentPatch%cwd_bg_out(:) = 0.0_r8 end subroutine non_canopy_derivs @@ -636,6 +635,8 @@ subroutine seeds_in( cp_pnt ) currentSite => currentPatch%siteptr currentPatch%seeds_in(:) = 0.0_r8 + currentPatch%seed_rain_flux(:) = 0.0_r8 + currentCohort => currentPatch%tallest do while (associated(currentCohort)) p = currentCohort%pft @@ -649,6 +650,7 @@ subroutine seeds_in( cp_pnt ) if (EXTERNAL_RECRUITMENT == 1) then !external seed rain - needed to prevent extinction do p = 1,numpft_ed currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + EDecophyscon%seed_rain(p) !KgC/m2/year + currentPatch%seed_rain_flux(p) = currentPatch%seed_rain_flux(p) + EDecophyscon%seed_rain(p) !KgC/m2/year enddo endif currentPatch => currentPatch%younger @@ -1251,4 +1253,6 @@ subroutine cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_i end subroutine cwd_out + + end module EDPhysiologyMod diff --git a/components/clm/src/ED/biogeophys/EDSurfaceAlbedoMod.F90 b/components/clm/src/ED/biogeophys/EDSurfaceAlbedoMod.F90 index 02fd161a06..d601d5bf49 100644 --- a/components/clm/src/ED/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/components/clm/src/ED/biogeophys/EDSurfaceAlbedoMod.F90 @@ -857,10 +857,10 @@ subroutine ED_Norman_Radiation (bounds, & enddo enddo if (lai_change(1,2,1).gt.0.0.and.lai_change(1,2,2).gt.0.0)then - write(iulog,*) 'lai_change(1,2,12)',lai_change(1,2,1:4) +! write(iulog,*) 'lai_change(1,2,12)',lai_change(1,2,1:4) endif if (lai_change(1,2,2).gt.0.0.and.lai_change(1,2,3).gt.0.0)then - write(iulog,*) ' lai_change (1,2,23)',lai_change(1,2,1:4) +! write(iulog,*) ' lai_change (1,2,23)',lai_change(1,2,1:4) endif if (lai_change(1,1,3).gt.0.0.and.lai_change(1,1,2).gt.0.0)then ! NO-OP diff --git a/components/clm/src/ED/main/EDCLMLinkMod.F90 b/components/clm/src/ED/main/EDCLMLinkMod.F90 index e405fc21b4..aced5963c3 100755 --- a/components/clm/src/ED/main/EDCLMLinkMod.F90 +++ b/components/clm/src/ED/main/EDCLMLinkMod.F90 @@ -10,15 +10,26 @@ module EDCLMLinkMod use decompMod , only : bounds_type use clm_varpar , only : nclmax, nlevcan_ed, numpft, numcft, mxpft use clm_varctl , only : iulog - use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, nlevsclass_ed - use EDtypesMod , only : sclass_ed,nlevsclass_ed,AREA + + use EDtypesMod , only : ed_site_type, ed_cohort_type, ed_patch_type, ncwd + use EDtypesMod , only : sclass_ed, nlevsclass_ed, AREA + use CanopyStateType , only : canopystate_type + use clm_varctl , only : use_vertsoilc use EDParamsMod , only : ED_val_ag_biomass + use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type + use SoilBiogeochemCarbonStatetype , only : soilbiogeochem_carbonstate_type + use clm_time_manager , only : is_beg_curr_day, get_step_size, get_nstep + use shr_const_mod, only: SHR_CONST_CDAY ! implicit none private ! logical :: DEBUG = .false. ! for debugging this module (EDCLMLinkMod.F90) + + ! !PUBLIC DATA MEMBERS + real(r8), public :: cwd_fcel_ed + real(r8), public :: cwd_flig_ed type, public :: ed_clm_type @@ -66,14 +77,14 @@ module EDCLMLinkMod real(r8), pointer, private :: ED_bleaf_patch (:) ! kGC/m2 Total leaf biomass. real(r8), pointer, private :: ED_biomass_patch (:) ! kGC/m2 Total biomass. - real(r8), pointer, private :: storvegc_patch (:) ! (gC/m2) stored vegetation carbon, excluding cpool - real(r8), pointer, private :: dispvegc_patch (:) ! (gC/m2) displayed veg carbon, excluding storage and cpool - real(r8), pointer, private :: leafc_patch (:) ! (gC/m2) leaf C - real(r8), pointer, private :: livestemc_patch (:) ! (gC/m2) live stem C - real(r8), pointer, private :: deadstemc_patch (:) ! (gC/m2) dead stem C - real(r8), pointer, private :: livestemn_patch (:) ! (gN/m2) live stem N - real(r8), pointer, private :: npp_patch (:) ! (gC/m2/s) patch net primary production + ! real(r8), pointer, private :: storvegc_patch (:) ! (gC/m2) stored vegetation carbon, excluding cpool + ! real(r8), pointer, private :: dispvegc_patch (:) ! (gC/m2) displayed veg carbon, excluding storage and cpool + ! real(r8), pointer, private :: leafc_patch (:) ! (gC/m2) leaf C + ! real(r8), pointer, private :: livestemc_patch (:) ! (gC/m2) live stem C + ! real(r8), pointer, private :: deadstemc_patch (:) ! (gC/m2) dead stem C + ! real(r8), pointer, private :: livestemn_patch (:) ! (gN/m2) live stem N + real(r8), pointer, private :: npp_patch (:) ! (gC/m2/s) patch net primary production real(r8), pointer, private :: gpp_patch (:) ! (gC/m2/s) patch gross primary production real(r8), pointer :: ed_gpp_gd_scpf (:,:) ! [kg/m2/yr] gross primary production @@ -106,6 +117,45 @@ module EDCLMLinkMod real(r8), pointer :: ed_m4_gd_scpf (:,:) ! [Stems/ha/yr] Mean Impact Mortality real(r8), pointer :: ed_m5_gd_scpf (:,:) ! [Stems/ha/yr] Mean Fire Mortality + ! litterfall fluxes of C from ED patches to BGC columns + + real(r8), pointer, public :: ED_c_to_litr_lab_c_col(:,:) !total labile litter coming from ED. gC/m3/s + real(r8), pointer, public :: ED_c_to_litr_cel_c_col(:,:) !total cellulose litter coming from ED. gC/m3/s + real(r8), pointer, public :: ED_c_to_litr_lig_c_col(:,:) !total lignin litter coming from ED. gC/m3/s + real(r8), pointer, private :: leaf_prof_col(:,:) !(1/m) profile of leaves + real(r8), pointer, private :: froot_prof_col(:,:,:) !(1/m) profile of fine roots + real(r8), pointer, private :: croot_prof_col(:,:) !(1/m) profile of coarse roots + real(r8), pointer, private :: stem_prof_col(:,:) !(1/m) profile of leaves + + ! summary carbon fluxes at the column level + real(r8), pointer, public :: nep_col(:) ! [gC/m2/s] Net ecosystem production, i.e. fast-timescale carbon balance that does not include disturbance + real(r8), pointer, public :: nep_timeintegrated_col(:) ! [gC/m2/s] Net ecosystem production, i.e. fast-timescale carbon balance that does not include disturbance + real(r8), pointer, public :: npp_timeintegrated_col(:) ! [gC/m2/s] Net primary production, time integrated at column level for carbon balance checking + real(r8), pointer, public :: hr_timeintegrated_col(:) ! [gC/m2/s] heterotrophic respiration, time integrated for carbon balance checking + real(r8), pointer, public :: nbp_col(:) ! [gC/m2/s] Net biosphere production, i.e. slow-timescale carbon balance that integrates to total carbon change + real(r8), pointer, public :: npp_hifreq_col(:) ! [gC/m2/s] Net primary production at the fast timescale, aggregated to the column level + real(r8), pointer, public :: fire_c_to_atm_col(:) ! [gC/m2/s] total fire carbon loss to atmosphere + real(r8), pointer, public :: ed_to_bgc_this_edts_col(:) ! [gC/m2/s] total flux of carbon from ED to BGC models on current ED timestep + real(r8), pointer, public :: ed_to_bgc_last_edts_col(:) ! [gC/m2/s] total flux of carbon from ED to BGC models on prior ED timestep + real(r8), pointer, public :: seed_rain_flux_col(:) ! [gC/m2/s] total flux of carbon from seed rain + + ! summary carbon states at the column level + real(r8), pointer, public :: totecosysc_col(:) ! [gC/m2] Total ecosystem carbon at the column level, including vegetation, CWD, litter, and soil pools + real(r8), pointer, public :: totecosysc_old_col(:) ! [gC/m2] Total ecosystem C at the column level from last call to balance check + real(r8), pointer, public :: totedc_col(:) ! [gC/m2] Total ED carbon at the column level, including vegetation, CWD, seeds, and ED litter + real(r8), pointer, public :: totedc_old_col(:) ! [gC/m2] Total ED C at the column level from last call to balance check + real(r8), pointer, public :: totbgcc_col(:) ! [gC/m2] Total BGC carbon at the column level, including litter, and soil pools + real(r8), pointer, public :: totbgcc_old_col(:) ! [gC/m2] Total BGC C at the column level from last call to balance check + real(r8), pointer, public :: biomass_stock_col(:) ! [gC/m2] total biomass at the column level in gC / m2 + real(r8), pointer, public :: ed_litter_stock_col(:) ! [gC/m2] ED litter at the column level in gC / m2 + real(r8), pointer, public :: cwd_stock_col(:) ! [gC/m2] ED CWD at the column level in gC / m2 + real(r8), pointer, public :: seed_stock_col(:) ! [gC/m2] ED seed mass carbon at the column level in gC / m2 + + ! carbon balance errors. at some point we'll reduce these to close to zero and delete, but for now we'll just keep[ track of them + real(r8), pointer, public :: cbalance_error_ed_col(:) ! [gC/m2/s] total carbon balance error for the ED side + real(r8), pointer, public :: cbalance_error_bgc_col(:) ! [gC/m2/s] total carbon balance error for the BGC side + real(r8), pointer, public :: cbalance_error_total_col(:) ! [gC/m2/s] total carbon balance error for the whole thing + contains ! Public routines @@ -113,13 +163,16 @@ module EDCLMLinkMod procedure , public :: Restart procedure , public :: SetValues procedure , public :: ed_clm_link + procedure , public :: Summary + procedure , public :: ED_BGC_Carbon_Balancecheck ! Private routines procedure , private :: ed_clm_leaf_area_profile procedure , private :: ed_update_history_variables procedure , private :: InitAllocate procedure , private :: InitHistory - procedure , private :: InitCold +! procedure , private :: InitCold + procedure , private :: flux_into_litter_pools end type ed_clm_type @@ -141,7 +194,7 @@ subroutine Init(this, bounds) call this%InitAllocate(bounds) call this%InitHistory(bounds) - call this%InitCold(bounds) + !call this%InitCold(bounds) end subroutine Init @@ -150,7 +203,8 @@ subroutine InitAllocate(this, bounds) ! ! !USES: use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : nlevgrnd + use clm_varpar , only : nlevgrnd, nlevdecomp_full + use EDtypesMod , only : numpft_ed ! ! !ARGUMENTS: class (ed_clm_type) :: this @@ -164,6 +218,7 @@ subroutine InitAllocate(this, bounds) begp = bounds%begp; endp = bounds%endp begc = bounds%begc; endc = bounds%endc + begg = bounds%begg; endg = bounds%endg allocate(this%trimming_patch (begp:endp)) ; this%trimming_patch (:) = 0.0_r8 @@ -200,17 +255,53 @@ subroutine InitAllocate(this, bounds) allocate(this%ED_bleaf_patch (begp:endp)) ; this%ED_bleaf_patch (:) = 0.0_r8 allocate(this%ED_biomass_patch (begp:endp)) ; this%ED_biomass_patch (:) = 0.0_r8 - allocate(this%storvegc_patch (begp:endp)) ; this%storvegc_patch (:) = nan - allocate(this%dispvegc_patch (begp:endp)) ; this%dispvegc_patch (:) = nan - allocate(this%leafc_patch (begp:endp)) ; this%leafc_patch (:) = nan - allocate(this%livestemc_patch (begp:endp)) ; this%livestemc_patch (:) = nan - allocate(this%deadstemc_patch (begp:endp)) ; this%deadstemc_patch (:) = nan - allocate(this%livestemn_patch (begp:endp)) ; this%livestemn_patch (:) = nan + ! allocate(this%storvegc_patch (begp:endp)) ; this%storvegc_patch (:) = nan + ! allocate(this%dispvegc_patch (begp:endp)) ; this%dispvegc_patch (:) = nan + ! allocate(this%leafc_patch (begp:endp)) ; this%leafc_patch (:) = nan + ! allocate(this%livestemc_patch (begp:endp)) ; this%livestemc_patch (:) = nan + ! allocate(this%deadstemc_patch (begp:endp)) ; this%deadstemc_patch (:) = nan + ! allocate(this%livestemn_patch (begp:endp)) ; this%livestemn_patch (:) = nan allocate(this%gpp_patch (begp:endp)) ; this%gpp_patch (:) = nan allocate(this%npp_patch (begp:endp)) ; this%npp_patch (:) = nan - begg = bounds%begg; endg = bounds%endg + allocate(this%ED_c_to_litr_lab_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_lab_c_col (:,:) = nan + allocate(this%ED_c_to_litr_cel_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_cel_c_col (:,:) = nan + allocate(this%ED_c_to_litr_lig_c_col (begc:endc,1:nlevdecomp_full)) ; this%ED_c_to_litr_lig_c_col (:,:) = nan + + allocate(this%leaf_prof_col (begc:endc,1:nlevdecomp_full)) ; this%leaf_prof_col (:,:) = nan + allocate(this%froot_prof_col (begc:endc,1:numpft_ed,1:nlevdecomp_full)); this%froot_prof_col (:,:,:) = nan + allocate(this%croot_prof_col (begc:endc,1:nlevdecomp_full)) ; this%croot_prof_col (:,:) = nan + allocate(this%stem_prof_col (begc:endc,1:nlevdecomp_full)) ; this%stem_prof_col (:,:) = nan + + allocate(this%ed_to_bgc_this_edts_col (begc:endc)) ; this%ed_to_bgc_this_edts_col (:) = nan + allocate(this%ed_to_bgc_last_edts_col (begc:endc)) ; this%ed_to_bgc_last_edts_col (:) = nan + allocate(this%seed_rain_flux_col (begc:endc)) ; this%seed_rain_flux_col (:) = nan + + allocate(this%nep_col (begc:endc)) ; this%nep_col (:) = nan + allocate(this%nep_timeintegrated_col (begc:endc)) ; this%nep_timeintegrated_col (:) = nan + allocate(this%npp_timeintegrated_col (begc:endc)) ; this%npp_timeintegrated_col (:) = nan + allocate(this%hr_timeintegrated_col (begc:endc)) ; this%hr_timeintegrated_col (:) = nan + + allocate(this%nbp_col (begc:endc)) ; this%nbp_col (:) = nan + allocate(this%npp_hifreq_col (begc:endc)) ; this%npp_hifreq_col (:) = nan + allocate(this%fire_c_to_atm_col (begc:endc)) ; this%fire_c_to_atm_col (:) = nan + + allocate(this%totecosysc_col (begc:endc)) ; this%totecosysc_col (:) = nan + allocate(this%totecosysc_old_col (begc:endc)) ; this%totecosysc_old_col (:) = nan + allocate(this%totedc_col (begc:endc)) ; this%totedc_col (:) = nan + allocate(this%totedc_old_col (begc:endc)) ; this%totedc_old_col (:) = nan + allocate(this%totbgcc_col (begc:endc)) ; this%totbgcc_col (:) = nan + allocate(this%totbgcc_old_col (begc:endc)) ; this%totbgcc_old_col (:) = nan + allocate(this%biomass_stock_col (begc:endc)) ; this%biomass_stock_col (:) = nan + allocate(this%ed_litter_stock_col (begc:endc)) ; this%ed_litter_stock_col (:) = nan + allocate(this%cwd_stock_col (begc:endc)) ; this%cwd_stock_col (:) = nan + allocate(this%seed_stock_col (begc:endc)) ; this%seed_stock_col (:) = nan + + allocate(this%cbalance_error_ed_col (begc:endc)) ; this%cbalance_error_ed_col (:) = nan + allocate(this%cbalance_error_bgc_col (begc:endc)) ; this%cbalance_error_bgc_col (:) = nan + allocate(this%cbalance_error_total_col (begc:endc)) ; this%cbalance_error_total_col (:) = nan + allocate(this%ed_gpp_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_gpp_gd_scpf (:,:) = 0.0_r8 allocate(this%ed_npp_totl_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_totl_gd_scpf (:,:) = 0.0_r8 allocate(this%ed_npp_leaf_gd_scpf (begg:endg,1:nlevsclass_ed*mxpft)); this%ed_npp_leaf_gd_scpf (:,:) = 0.0_r8 @@ -288,20 +379,20 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='Scaling factor between tree basal area and canopy area', & ptr_patch=this%canopy_spread_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTbiomass', units='kgC/m2', type2d='levgrnd', & + call hist_addfld2d (fname='PFTbiomass', units='gC/m2', type2d='levgrnd', & avgflag='A', long_name='total PFT level biomass', & ptr_patch=this%PFTbiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTleafbiomass', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTleafbiomass', units='gC/m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level leaf biomass', & ptr_patch=this%PFTleafbiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTstorebiomass', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTstorebiomass', units='gC/m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level stored biomass', & ptr_patch=this%PFTstorebiomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld2d (fname='PFTnindivs', units='kgC/m2', type2d='levgrnd', & - avgflag='A', long_name='total PFT level biomass', & + call hist_addfld2d (fname='PFTnindivs', units='indiv / m2', type2d='levgrnd', & + avgflag='A', long_name='total PFT level number of individuals', & ptr_patch=this%PFTnindivs_patch, set_lake=0._r8, set_urb=0._r8) call hist_addfld1d (fname='FIRE_NESTEROV_INDEX', units='none', & @@ -348,59 +439,55 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='spitfire fuel surface/volume ', & ptr_patch=this%fire_fuel_sav_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='TFC_ROS', units='m', & - avgflag='A', long_name='spitfire fuel surface/volume ', & - ptr_patch=this%TFC_ROS_patch, set_lake=0._r8, set_urb=0._r8) - - call hist_addfld1d (fname='SUM_FUEL', units=' KgC m-2 y-1', & - avgflag='A', long_name='Litter flux in leaves', & + call hist_addfld1d (fname='SUM_FUEL', units='gC m-2', & + avgflag='A', long_name='total ground fuel related to ros (omits 1000hr fuels)', & ptr_patch=this%sum_fuel_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='LITTER_IN', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='LITTER_IN', units='gC m-2 s-1', & avgflag='A', long_name='Litter flux in leaves', & ptr_patch=this%litter_in_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='LITTER_OUT', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='LITTER_OUT', units='gC m-2 s-1', & avgflag='A', long_name='Litter flux out leaves', & ptr_patch=this%litter_out_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_BANK', units=' KgC m-2', & + call hist_addfld1d (fname='SEED_BANK', units='gC m-2', & avgflag='A', long_name='Total Seed Mass of all PFTs', & ptr_patch=this%seed_bank_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEEDS_IN', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEEDS_IN', units='gC m-2 s-1', & avgflag='A', long_name='Seed Production Rate', & ptr_patch=this%seeds_in_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_GERMINATION', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEED_GERMINATION', units='gC m-2 s-1', & avgflag='A', long_name='Seed mass converted into new cohorts', & ptr_patch=this%seed_germination_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='SEED_DECAY', units=' KgC m-2 y-1', & + call hist_addfld1d (fname='SEED_DECAY', units='gC m-2 s-1', & avgflag='A', long_name='Seed mass decay', & ptr_patch=this%seed_decay_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bstore', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bstore', units='gC m-2', & avgflag='A', long_name='ED stored biomass', & ptr_patch=this%ED_bstore_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bdead', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bdead', units='gC m-2', & avgflag='A', long_name='ED dead biomass', & ptr_patch=this%ED_bdead_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_balive', units=' KgC m-2', & + call hist_addfld1d (fname='ED_balive', units='gC m-2', & avgflag='A', long_name='ED live biomass', & ptr_patch=this%ED_balive_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_bleaf', units=' KgC m-2', & + call hist_addfld1d (fname='ED_bleaf', units='gC m-2', & avgflag='A', long_name='ED leaf biomass', & ptr_patch=this%ED_bleaf_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='ED_biomass', units=' KgC m-2', & + call hist_addfld1d (fname='ED_biomass', units='gC m-2', & avgflag='A', long_name='ED total biomass', & ptr_patch=this%ED_biomass_patch, set_lake=0._r8, set_urb=0._r8) - call hist_addfld1d (fname='RB', units=' s m-1', & + call hist_addfld1d (fname='RB', units='s m-1', & avgflag='A', long_name='leaf boundary resistance', & ptr_patch=this%rb_patch, set_lake=0._r8, set_urb=0._r8) @@ -408,35 +495,35 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='potential evap', & ptr_patch=this%efpot_patch, set_lake=0._r8, set_urb=0._r8) - this%dispvegc_patch(begp:endp) = spval - call hist_addfld1d (fname='DISPVEGC', units='gC/m^2', & - avgflag='A', long_name='displayed veg carbon, excluding storage and cpool', & - ptr_patch=this%dispvegc_patch) + ! this%dispvegc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='DISPVEGC', units='gC/m^2', & + ! avgflag='A', long_name='displayed veg carbon, excluding storage and cpool', & + ! ptr_patch=this%dispvegc_patch) - this%storvegc_patch(begp:endp) = spval - call hist_addfld1d (fname='STORVEGC', units='gC/m^2', & - avgflag='A', long_name='stored vegetation carbon, excluding cpool', & - ptr_patch=this%storvegc_patch) + ! this%storvegc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='STORVEGC', units='gC/m^2', & + ! avgflag='A', long_name='stored vegetation carbon, excluding cpool', & + ! ptr_patch=this%storvegc_patch) - this%leafc_patch(begp:endp) = spval - call hist_addfld1d (fname='LEAFC', units='gC/m^2', & - avgflag='A', long_name='leaf C', & - ptr_patch=this%leafc_patch) + ! this%leafc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LEAFC', units='gC/m^2', & + ! avgflag='A', long_name='leaf C', & + ! ptr_patch=this%leafc_patch) - this%livestemc_patch(begp:endp) = spval - call hist_addfld1d (fname='LIVESTEMC', units='gC/m^2', & - avgflag='A', long_name='live stem C', & - ptr_patch=this%livestemc_patch) + ! this%livestemc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LIVESTEMC', units='gC/m^2', & + ! avgflag='A', long_name='live stem C', & + ! ptr_patch=this%livestemc_patch) - this%deadstemc_patch(begp:endp) = spval - call hist_addfld1d (fname='DEADSTEMC', units='gC/m^2', & - avgflag='A', long_name='dead stem C', & - ptr_patch=this%deadstemc_patch) + ! this%deadstemc_patch(begp:endp) = spval + ! call hist_addfld1d (fname='DEADSTEMC', units='gC/m^2', & + ! avgflag='A', long_name='dead stem C', & + ! ptr_patch=this%deadstemc_patch) - this%livestemn_patch(begp:endp) = spval - call hist_addfld1d (fname='LIVESTEMN', units='gN/m^2', & - avgflag='A', long_name='live stem N', & - ptr_patch=this%livestemn_patch) + ! this%livestemn_patch(begp:endp) = spval + ! call hist_addfld1d (fname='LIVESTEMN', units='gN/m^2', & + ! avgflag='A', long_name='live stem N', & + ! ptr_patch=this%livestemn_patch) this%gpp_patch(begp:endp) = spval call hist_addfld1d (fname='GPP', units='gC/m^2/s', & @@ -448,6 +535,97 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='net primary production', & ptr_patch=this%npp_patch) + this%nep_col(begc:endc) = spval + call hist_addfld1d (fname='NEP', units='gC/m^2/s', & + avgflag='A', long_name='net ecosystem production', & + ptr_col=this%nep_col) + + this%fire_c_to_atm_col(begc:endc) = spval + call hist_addfld1d (fname='Fire_Closs', units='gC/m^2/s', & + avgflag='A', long_name='ED/SPitfire Carbon loss to atmosphere', & + ptr_col=this%fire_c_to_atm_col) + + this%nbp_col(begc:endc) = spval + call hist_addfld1d (fname='NBP', units='gC/m^2/s', & + avgflag='A', long_name='net biosphere production', & + ptr_col=this%nbp_col) + + this%npp_hifreq_col(begc:endc) = spval + call hist_addfld1d (fname='NPP_hifreq', units='gC/m^2/s', & + avgflag='A', long_name='net primary production at high frequency', & + ptr_col=this%npp_hifreq_col,default='inactive') + + this%totecosysc_col(begc:endc) = spval + call hist_addfld1d (fname='TOTECOSYSC', units='gC/m^2', & + avgflag='A', long_name='total ecosystem carbon', & + ptr_col=this%totecosysc_col) + + this%cbalance_error_ed_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_ED', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error on ED side', & + ptr_col=this%cbalance_error_ed_col) + + this%cbalance_error_bgc_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_BGC', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error on BGC side', & + ptr_col=this%cbalance_error_bgc_col) + + this%cbalance_error_total_col(begc:endc) = spval + call hist_addfld1d (fname='CBALANCE_ERROR_TOTAL', units='gC/m^2/s', & + avgflag='A', long_name='total carbon balance error total', & + ptr_col=this%cbalance_error_total_col) + + this%biomass_stock_col(begc:endc) = spval + call hist_addfld1d (fname='BIOMASS_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total ED biomass carbon at the column level', & + ptr_col=this%biomass_stock_col) + + this%ed_litter_stock_col(begc:endc) = spval + call hist_addfld1d (fname='ED_LITTER_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total ED litter carbon at the column level', & + ptr_col=this%ed_litter_stock_col) + + this%cwd_stock_col(begc:endc) = spval + call hist_addfld1d (fname='CWD_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total CWD carbon at the column level', & + ptr_col=this%cwd_stock_col) + + this%seed_stock_col(begc:endc) = spval + call hist_addfld1d (fname='SEED_STOCK_COL', units='gC/m^2', & + avgflag='A', long_name='total seed carbon at the column level', & + ptr_col=this%seed_stock_col) + + !!! carbon fluxes into soil grid (dimensioned depth x column) + this%ED_c_to_litr_lab_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_lab_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_lab_c', & + ptr_col=this%ED_c_to_litr_lab_c_col) + + this%ED_c_to_litr_cel_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_cel_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_cel_c', & + ptr_col=this%ED_c_to_litr_cel_c_col) + + this%ED_c_to_litr_lig_c_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='ED_c_to_litr_lig_c', units='gC/m^2/s', type2d='levdcmp', & + avgflag='A', long_name='ED_c_to_litr_lig_c', & + ptr_col=this%ED_c_to_litr_lig_c_col) + + this%leaf_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='leaf_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='leaf_prof', & + ptr_col=this%leaf_prof_col,default='inactive') + + this%croot_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='croot_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='croot_prof', & + ptr_col=this%croot_prof_col,default='inactive') + + this%stem_prof_col(begc:endc,1:nlevdecomp_full) = spval + call hist_addfld_decomp (fname='stem_prof', units='1/m', type2d='levdcmp', & + avgflag='A', long_name='stem_prof', & + ptr_col=this%stem_prof_col,default='inactive') + ! Carbon Flux (grid dimension x scpf) ! ============================================================== @@ -523,26 +701,25 @@ subroutine InitHistory(this, bounds) end subroutine InitHistory !----------------------------------------------------------------------- - subroutine InitCold(this, bounds) - ! - ! !DESCRIPTION: - ! Initialize relevant time varying variables - ! - ! !ARGUMENTS: - class (ed_clm_type) :: this - type(bounds_type), intent(in) :: bounds - ! - ! !LOCAL VARIABLES: - integer :: p - !----------------------------------------------------------------------- - - do p = bounds%begp,bounds%endp - this%dispvegc_patch(p) = 0._r8 - this%storvegc_patch(p) = 0._r8 - end do - - end subroutine InitCold - + ! subroutine InitCold(this, bounds) + ! ! + ! ! !DESCRIPTION: + ! ! Initialize relevant time varying variables + ! ! + ! ! !ARGUMENTS: + ! class (ed_clm_type) :: this + ! type(bounds_type), intent(in) :: bounds + ! ! + ! ! !LOCAL VARIABLES: + ! integer :: p + ! !----------------------------------------------------------------------- + + ! ! do p = bounds%begp,bounds%endp + ! ! this%dispvegc_patch(p) = 0._r8 + ! ! this%storvegc_patch(p) = 0._r8 + ! ! end do + + ! end subroutine InitCold !----------------------------------------------------------------------- subroutine Restart ( this, bounds, ncid, flag ) ! @@ -552,6 +729,8 @@ subroutine Restart ( this, bounds, ncid, flag ) ! !USES: use restUtilMod use ncdio_pio + ! use EDtypesMod , only : numpft_ed + ! ! !ARGUMENTS: class (ed_clm_type) :: this @@ -561,23 +740,173 @@ subroutine Restart ( this, bounds, ncid, flag ) ! ! !LOCAL VARIABLES: logical :: readvar + real(r8), pointer :: ptr2d(:,:) ! temp. pointers for slicing larger arrays + real(r8), pointer :: ptr1d(:) ! temp. pointers for slicing larger arrays + ! character(LEN=3) :: istr1 + ! integer :: k !------------------------------------------------------------------------ - call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%leafc_patch) - - call restartvar(ncid=ncid, flag=flag, varname='livestemc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%livestemc_patch) - - call restartvar(ncid=ncid, flag=flag, varname='deadstemc', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%deadstemc_patch) + ! call restartvar(ncid=ncid, flag=flag, varname='leafc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%leafc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='livestemc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%livestemc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='deadstemc', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%deadstemc_patch) + + ! call restartvar(ncid=ncid, flag=flag, varname='livestemn', xtype=ncd_double, & + ! dim1name='pft', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=this%livestemn_patch) + + ptr1d => this%nep_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='nep_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%npp_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='npp_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%hr_timeintegrated_col(:) + call restartvar(ncid=ncid, flag=flag, varname='hr_timeintegrated_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totecosysc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totecosysc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_ed_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_ed_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_bgc_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_bgc_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%cbalance_error_total_col(:) + call restartvar(ncid=ncid, flag=flag, varname='cbalance_error_total_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totedc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totedc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%totbgcc_old_col(:) + call restartvar(ncid=ncid, flag=flag, varname='totbgcc_old_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ed_to_bgc_this_edts_col(:) + call restartvar(ncid=ncid, flag=flag, varname='ed_to_bgc_this_edts_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ed_to_bgc_last_edts_col(:) + call restartvar(ncid=ncid, flag=flag, varname='ed_to_bgc_last_edts_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%seed_rain_flux_col(:) + call restartvar(ncid=ncid, flag=flag, varname='seed_rain_flux_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + if (use_vertsoilc) then + ptr2d => this%ED_c_to_litr_lab_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lab_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ptr2d => this%ED_c_to_litr_cel_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_cel_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ptr2d => this%ED_c_to_litr_lig_c_col + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lig_c_col', xtype=ncd_double, & + dim1name='column', dim2name='levgrnd', switchdim=.true., & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%leaf_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='leaf_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%croot_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='croot_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! ptr2d => this%stem_prof_col + ! call restartvar(ncid=ncid, flag=flag, varname='stem_prof_col', xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + + ! do k = 1, numpft_ed + ! write(istr1,"(I3.3)") k + ! ptr2d => this%froot_prof_col(:,k,:) + ! call restartvar(ncid=ncid, flag=flag, varname='froot_prof_col_PFT'//istr1, xtype=ncd_double, & + ! dim1name='column', dim2name='levgrnd', switchdim=.true., & + ! long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr2d) + ! end do + else + ptr1d => this%ED_c_to_litr_lab_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lab_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ED_c_to_litr_cel_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_cel_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ptr1d => this%ED_c_to_litr_lig_c_col(:,1) + call restartvar(ncid=ncid, flag=flag, varname='ED_c_to_litr_lig_c_col', xtype=ncd_double, & + dim1name='column', long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%leaf_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='leaf_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%croot_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='croot_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! ptr1d => this%stem_prof_col(:,1) + ! call restartvar(ncid=ncid, flag=flag, varname='stem_prof_col', xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + + ! do k = 1, numpft_ed + ! write(istr1,"(I3.3)") k + ! ptr1d => this%froot_prof_col(:,k,1) + ! call restartvar(ncid=ncid, flag=flag, varname='froot_prof_col_PFT'//istr1, xtype=ncd_double, & + ! dim1name='column', long_name='', units='', & + ! interpinic_flag='interp', readvar=readvar, data=ptr1d) + ! end do + end if - call restartvar(ncid=ncid, flag=flag, varname='livestemn', xtype=ncd_double, & - dim1name='pft', long_name='', units='', & - interpinic_flag='interp', readvar=readvar, data=this%livestemn_patch) end subroutine Restart @@ -870,7 +1199,10 @@ subroutine ed_clm_link( this, bounds, ed_allsites_inst, ed_phenology_inst, & end do !grid loop - call this%ed_update_history_variables( bounds, ed_allsites_inst(begg:endg), & + call this%flux_into_litter_pools(bounds, ed_allsites_inst(begg:endg), firstsoilpatch, & + canopystate_inst) + + call this%ed_update_history_variables(bounds, ed_allsites_inst(begg:endg), & firstsoilpatch, ed_Phenology_inst, canopystate_inst) end associate @@ -1104,16 +1436,16 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & write(iulog,*) 'EDCLMLinkMod II ',p,ED_bstore(p) endif - ED_bleaf(p) = ED_bleaf(p) + n_density * currentCohort%bl - ED_bstore(p) = ED_bstore(p) + n_density * currentCohort%bstore - ED_biomass(p) = ED_biomass(p) + n_density * currentCohort%b - ED_bdead(p) = ED_bdead(p) + n_density * currentCohort%bdead - ED_balive(p) = ED_balive(p) + n_density * currentCohort%balive - npp(p) = npp(p) + n_density * currentCohort%npp - gpp(p) = gpp(p) + n_density * currentCohort%gpp - PFTbiomass(p,ft) = PFTbiomass(p,ft) + n_density * currentCohort%b - PFTleafbiomass(p,ft) = PFTleafbiomass(p,ft) + n_density * currentCohort%bl - PFTstorebiomass(p,ft) = PFTstorebiomass(p,ft) + n_density * currentCohort%bstore + ED_bleaf(p) = ED_bleaf(p) + n_density * currentCohort%bl * 1.e3_r8 + ED_bstore(p) = ED_bstore(p) + n_density * currentCohort%bstore * 1.e3_r8 + ED_biomass(p) = ED_biomass(p) + n_density * currentCohort%b * 1.e3_r8 + ED_bdead(p) = ED_bdead(p) + n_density * currentCohort%bdead * 1.e3_r8 + ED_balive(p) = ED_balive(p) + n_density * currentCohort%balive * 1.e3_r8 + npp(p) = npp(p) + n_density * currentCohort%npp * 1.e3_r8 / (365. * SHR_CONST_CDAY) + gpp(p) = gpp(p) + n_density * currentCohort%gpp * 1.e3_r8 / (365. * SHR_CONST_CDAY) + PFTbiomass(p,ft) = PFTbiomass(p,ft) + n_density * currentCohort%b * 1.e3_r8 + PFTleafbiomass(p,ft) = PFTleafbiomass(p,ft) + n_density * currentCohort%bl * 1.e3_r8 + PFTstorebiomass(p,ft) = PFTstorebiomass(p,ft) + n_density * currentCohort%bstore * 1.e3_r8 PFTnindivs(p,ft) = PFTnindivs(p,ft) + currentCohort%n dbh = currentCohort%dbh !-0.5*(1./365.25)*currentCohort%ddbhdt @@ -1188,13 +1520,13 @@ subroutine ed_update_history_variables( this, bounds, ed_allsites_inst, & fire_fuel_eff_moist(p) = currentPatch%fuel_eff_moist fire_fuel_sav(p) = currentPatch%fuel_sav fire_fuel_mef(p) = currentPatch%fuel_mef - sum_fuel(p) = currentPatch%sum_fuel - litter_in(p) = sum(currentPatch%CWD_AG_in) +sum(currentPatch%leaf_litter_in) - litter_out(p) = sum(currentPatch%CWD_AG_out)+sum(currentPatch%leaf_litter_out) - seed_bank(p) = sum(currentPatch%seed_bank) - seeds_in(p) = sum(currentPatch%seeds_in) - seed_decay(p) = sum(currentPatch%seed_decay) - seed_germination(p) = sum(currentPatch%seed_germination) + sum_fuel(p) = currentPatch%sum_fuel * 1.e3_r8 + litter_in(p) = (sum(currentPatch%CWD_AG_in) +sum(currentPatch%leaf_litter_in)) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + litter_out(p) = (sum(currentPatch%CWD_AG_out)+sum(currentPatch%leaf_litter_out)) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_bank(p) = sum(currentPatch%seed_bank) * 1.e3_r8 + seeds_in(p) = sum(currentPatch%seeds_in) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_decay(p) = sum(currentPatch%seed_decay) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY + seed_germination(p) = sum(currentPatch%seed_germination) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY canopy_spread(p) = currentPatch%spread(1) area_plant(p) = currentPatch%total_canopy_area /currentPatch%area area_trees(p) = currentPatch%total_tree_area /currentPatch%area @@ -1718,4 +2050,674 @@ subroutine ed_clm_leaf_area_profile( this, currentSite, waterstate_inst, canopys end subroutine ed_clm_leaf_area_profile + + subroutine flux_into_litter_pools(this, bounds, ed_allsites_inst, firstsoilpatch, canopystate_inst) + ! Created by Charlie Koven and Rosie Fisher, 2014-2015 + ! take the flux out of the fragmenting litter pools and port into the decomposing litter pools. + ! in this implementation, decomposing pools are assumed to be humus and non-flammable, whereas fragmenting pools + ! are assumed to be physically fragmenting but not respiring. This is a simplification, but allows us to + ! a) reconcile the need to track both chemical fractions (lignin, cellulose, labile) and size fractions (trunk, branch, etc.) + ! b) to impose a realistic delay on the surge of nutrients into the litter pools when large CWD is added to the system via mortality + + ! because of the different subgrid structure, this subroutine includes the functionality that in the big-leaf BGC model, is calculated in SoilBiogeochemVerticalProfileMod + + ! The ED code is resolved at a daily timestep, but all of the CN-BGC fluxes are passed in as derivatives per second, + ! and then accumulated in the CNStateUpdate routines. One way of doing this is to pass back the CN fluxes per second, + ! and keep them constant for the whole day (making sure they are not overwritten. + ! This means that the carbon gets passed back and forth between the photosynthesis code (fast timestepping) to the ED code (slow timestepping), back to the BGC code (fast timestepping). + ! This means that the state update for the litter pools and for the CWD pools occurs at different timescales. + + use SFParamsMod, only: SF_val_max_decomp + use clm_varpar, only : mxpft,nlevdecomp, nlevdecomp_full + use EDTypesMod, only : AREA, numpft_ed + use SoilBiogeochemVerticalProfileMod, only: exponential_rooting_profile, pftspecific_rootingprofile, rootprof_exp, surfprof_exp + use pftconMod, only : pftcon + + use clm_varcon, only : zisoi, dzsoi_decomp, zsoi + use ColumnType , only : col + use abortutils , only : endrun + use shr_log_mod , only : errMsg => shr_log_errMsg + use EDParamsMod, only : ED_val_ag_biomass + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + integer :: firstsoilpatch(bounds%begg:bounds%endg) ! the first patch in this gridcell that is soil and thus bare... + type(canopystate_type) , intent(inout) :: canopystate_inst + ! + ! !LOCAL VARIABLES: + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort + type(ed_site_type), pointer :: cs + integer c,p,cc,j,g + real(r8) time_convert ! from year to seconds + real(r8) mass_convert ! ED uses kg, CLM uses g + integer :: begp,endp + integer :: begc,endc !bounds + !------------------------------------------------------------------------ + real(r8) :: cinput_rootfr(bounds%begc:bounds%endc, 1:numpft_ed, 1:nlevdecomp_full) ! column by pft root fraction used for calculating inputs + real(r8) :: croot_prof_perpatch(1:nlevdecomp_full) + real(r8) :: surface_prof(1:nlevdecomp_full) + integer :: ft, lev + real(r8) :: rootfr_tot(1:numpft_ed), biomass_bg_ft(1:numpft_ed) + real(r8) :: surface_prof_tot, leaf_prof_sum, stem_prof_sum, froot_prof_sum, biomass_bg_tot + real(r8) :: delta + + begp = bounds%begp; endp = bounds%endp + begc = bounds%begc; endc = bounds%endc + delta = 0.001_r8 + !no of seconds in a year. + time_convert = 365.0_r8*SHR_CONST_CDAY + + ! number of grams in a kilogram + mass_convert = 1000._r8 + + associate( & + ED_c_to_litr_lab_c => this%ED_c_to_litr_lab_c_col , & ! Output: total labile litter coming from ED. gC/m3/s + ED_c_to_litr_cel_c => this%ED_c_to_litr_cel_c_col , & ! Output: total cellulose litter coming from ED. gC/m3/s + ED_c_to_litr_lig_c => this%ED_c_to_litr_lig_c_col , & ! Output: total lignin litter coming from ED. gC/m3/s + leaf_prof => this%leaf_prof_col , & ! Output: (1/m) profile of leaves + froot_prof => this%froot_prof_col , & ! Output: (1/m) profile of fine roots + croot_prof => this%croot_prof_col , & ! Output: (1/m) profile of coarse roots + stem_prof => this%stem_prof_col , & ! Output: (1/m) profile of leaves + altmax_lastyear_indx => canopystate_inst%altmax_lastyear_indx_col & ! Input: [integer (:) ]frost table depth (m) + ) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! first calculate vertical profiles + ! define two types of profiles: + ! (1) a surface profile, for leaves and stem inputs, which is the same for each pft but differs from one column to the next to avoid inputting any C into permafrost + ! (2) a fine root profile, which is indexed by both column and pft, differs for each pft and also from one column to the next to avoid inputting any C into permafrost + ! (3) a coarse root profile, which is the root-biomass=weighted average of the fine root profiles + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + if (use_vertsoilc) then + + ! define a single shallow surface profile for surface additions (leaves, stems, and N deposition) + surface_prof(:) = 0._r8 + do j = 1, nlevdecomp + surface_prof(j) = exp(-surfprof_exp * zsoi(j)) / dzsoi_decomp(j) + end do + + ! initialize profiles to zero + leaf_prof(begc:endc, :) = 0._r8 + froot_prof(begc:endc, 1:numpft_ed, :) = 0._r8 + croot_prof(begc:endc, :) = 0._r8 + stem_prof(begc:endc, :) = 0._r8 + + cinput_rootfr(begc:endc, 1:numpft_ed, :) = 0._r8 + + do c = bounds%begc,bounds%endc + + ! calculate pft-specific rooting profiles in the absence of permafrost limitations + if ( exponential_rooting_profile ) then + if ( .not. pftspecific_rootingprofile ) then + ! define rooting profile from exponential parameters + do ft = 1, numpft_ed + do j = 1, nlevdecomp + cinput_rootfr(c,ft,j) = exp(-rootprof_exp * zsoi(j)) / dzsoi_decomp(j) + end do + end do + else + ! use beta distribution parameter from Jackson et al., 1996 + do ft = 1, numpft_ed + do j = 1, nlevdecomp + cinput_rootfr(c,ft,j) = ( pftcon%rootprof_beta(ft) ** (zisoi(j-1)*100._r8) - & + pftcon%rootprof_beta(ft) ** (zisoi(j)*100._r8) ) & + / dzsoi_decomp(j) + end do + end do + endif + else + do ft = 1,numpft_ed + do j = 1, nlevdecomp + ! use standard CLM root fraction profiles; + cinput_rootfr(c,ft,j) = ( .5_r8*( & + exp(-pftcon%roota_par(ft) * col%zi(c,lev-1)) & + + exp(-pftcon%rootb_par(ft) * col%zi(c,lev-1)) & + - exp(-pftcon%roota_par(ft) * col%zi(c,lev)) & + - exp(-pftcon%rootb_par(ft) * col%zi(c,lev)))) / dzsoi_decomp(j) + end do + end do + endif + ! + ! + ! now add permafrost constraint: integrate rootfr over active layer of soil column, + ! truncate below permafrost table where present, and rescale so that integral = 1 + do ft = 1,numpft_ed + rootfr_tot(ft) = 0._r8 + end do + surface_prof_tot = 0._r8 + ! + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + surface_prof_tot = surface_prof_tot + surface_prof(j) * dzsoi_decomp(j) + end do + do ft = 1,numpft_ed + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + rootfr_tot(ft) = rootfr_tot(ft) + cinput_rootfr(c,ft,j) * dzsoi_decomp(j) + end do + end do + ! + ! rescale the fine root profile + do ft = 1,numpft_ed + if ( (altmax_lastyear_indx(c) > 0) .and. (rootfr_tot(ft) > 0._r8) ) then + ! where there is not permafrost extending to the surface, integrate the profiles over the active layer + ! this is equivalent to integrating over all soil layers outside of permafrost regions + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + froot_prof(c,ft,j) = cinput_rootfr(c,ft,j) / rootfr_tot(ft) + end do + else + ! if fully frozen, or no roots, put everything in the top layer + froot_prof(c,ft,1) = 1._r8/dzsoi_decomp(1) + endif + end do + ! + ! rescale the shallow profiles + if ( (altmax_lastyear_indx(c) > 0) .and. (surface_prof_tot > 0._r8) ) then + ! where there is not permafrost extending to the surface, integrate the profiles over the active layer + ! this is equivalent to integrating over all soil layers outside of permafrost regions + do j = 1, min(max(altmax_lastyear_indx(c), 1), nlevdecomp) + ! set all surface processes to shallower profile + leaf_prof(c,j) = surface_prof(j)/ surface_prof_tot + stem_prof(c,j) = surface_prof(j)/ surface_prof_tot + end do + else + ! if fully frozen, or no roots, put everything in the top layer + leaf_prof(c,1) = 1._r8/dzsoi_decomp(1) + stem_prof(c,1) = 1._r8/dzsoi_decomp(1) + do j = 2, nlevdecomp + leaf_prof(c,j) = 0._r8 + stem_prof(c,j) = 0._r8 + end do + endif + end do + + else + + ! for one layer decomposition model, set profiles to unity + leaf_prof(bounds%begc:bounds%endc, :) = 1._r8 + froot_prof(bounds%begc:bounds%endc, 1:numpft_ed, :) = 1._r8 + stem_prof(bounds%begc:bounds%endc, :) = 1._r8 + + end if + + ! sanity check to ensure they integrate to 1 + do c = bounds%begc,bounds%endc + ! check the leaf and stem profiles + leaf_prof_sum = 0._r8 + stem_prof_sum = 0._r8 + do j = 1, nlevdecomp + leaf_prof_sum = leaf_prof_sum + leaf_prof(c,j) * dzsoi_decomp(j) + stem_prof_sum = stem_prof_sum + stem_prof(c,j) * dzsoi_decomp(j) + end do + if ( ( abs(stem_prof_sum - 1._r8) > delta ) .or. ( abs(leaf_prof_sum - 1._r8) > delta ) ) then + write(iulog, *) 'profile sums: ', leaf_prof_sum, stem_prof_sum + write(iulog, *) 'surface_prof: ', surface_prof + write(iulog, *) 'surface_prof_tot: ', surface_prof_tot + write(iulog, *) 'leaf_prof: ', leaf_prof(c,:) + write(iulog, *) 'stem_prof: ', stem_prof(c,:) + write(iulog, *) 'altmax_lastyear_indx: ', altmax_lastyear_indx(c) + write(iulog, *) 'dzsoi_decomp: ', dzsoi_decomp + call endrun(msg=' ERROR: sum-1 > delta'//errMsg(__FILE__, __LINE__)) + endif + ! now check each fine root profile + do ft = 1,numpft_ed + froot_prof_sum = 0._r8 + do j = 1, nlevdecomp + froot_prof_sum = froot_prof_sum + froot_prof(c,ft,j) * dzsoi_decomp(j) + end do + if ( ( abs(froot_prof_sum - 1._r8) > delta ) ) then + write(iulog, *) 'profile sums: ', froot_prof_sum + call endrun(msg=' ERROR: sum-1 > delta'//errMsg(__FILE__, __LINE__)) + endif + end do + end do + + ! zero the column-level C input variables + do c = bounds%begc,bounds%endc + do j = 1, nlevdecomp + ED_c_to_litr_lab_c(c,j) = 0._r8 + ED_c_to_litr_cel_c(c,j) = 0._r8 + ED_c_to_litr_lig_c(c,j) = 0._r8 + croot_prof(c,j) = 0._r8 + end do + end do + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! now disaggregate the inputs vertically, using the vertical profiles + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + currentPatch => ed_allsites_inst(g)%oldest_patch + + do while(associated(currentPatch)) + + cs => currentpatch%siteptr + cc = cs%clmcolumn + + ! the CWD pools lose information about which PFT they came from; for the stems this doesn't matter as they all have the same profile, + ! however for the coarse roots they may have different profiles. to approximately recover this information, loop over all cohorts in patch + ! to calculate the total root biomass in that patch of each pft, and then rescale the croot_prof as the weighted average of the froot_prof + biomass_bg_ft(1:numpft_ed) = 0._r8 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + biomass_bg_ft(currentCohort%pft) = biomass_bg_ft(currentCohort%pft) + currentCohort%b * (currentCohort%n / currentPatch%area) * (1.0_r8-ED_val_ag_biomass) + currentCohort => currentCohort%shorter + enddo !currentCohort + ! + biomass_bg_tot = 0._r8 + do ft = 1,numpft_ed + biomass_bg_tot = biomass_bg_tot + biomass_bg_ft(ft) + end do + ! + do j = 1, nlevdecomp + ! zero this for each patch + croot_prof_perpatch(j) = 0._r8 + end do + ! + if ( biomass_bg_tot .gt. 0._r8) then + do ft = 1,numpft_ed + do j = 1, nlevdecomp + croot_prof_perpatch(j) = croot_prof_perpatch(j) + froot_prof(cc,ft,j) * biomass_bg_ft(ft) / biomass_bg_tot + end do + end do + else ! no biomass + croot_prof_perpatch(1) = 1./dzsoi_decomp(1) + end if + ! + ! add croot_prof as weighted average (weighted by patch area) of croot_prof_perpatch + do j = 1, nlevdecomp + croot_prof(cc, j) = croot_prof(cc, j) + croot_prof_perpatch(j) * currentPatch%area / AREA + end do + ! + ! now disaggregate, vertically and by decomposition substrate type, the actual fluxes from CWD and litter pools + ! + ! do c = 1, ncwd + ! write(iulog,*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), cwd_fcel_ed, currentpatch%area/AREA + ! write(iulog,*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), cwd_fcel_ed, currentpatch%area/AREA + ! end do + ! do ft = 1,numpft_ed + ! write(iulog,*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), cwd_fcel_ed, currentpatch%area/AREA + ! write(iulog,*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), cwd_fcel_ed, currentpatch%area/AREA + ! end do + ! ! + ! CWD pools fragmenting into decomposing litter pools. + do c = 1, ncwd + do j = 1, nlevdecomp + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%CWD_AG_out(c) * cwd_fcel_ed * currentpatch%area/AREA * stem_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%CWD_AG_out(c) * cwd_flig_ed * currentpatch%area/AREA * stem_prof(cc,j) + ! + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%CWD_BG_out(c) * cwd_fcel_ed * currentpatch%area/AREA * croot_prof_perpatch(j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%CWD_BG_out(c) * cwd_flig_ed * currentpatch%area/AREA * croot_prof_perpatch(j) + end do + end do + + ! leaf and fine root pools. + do ft = 1,numpft_ed + do j = 1, nlevdecomp + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%leaf_litter_out(ft) * pftcon%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ! + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_flab(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_fcel(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%root_litter_out(ft) * pftcon%fr_flig(ft) * currentpatch%area/AREA * froot_prof(cc,ft,j) + ! + !! and seed_decay too. for now, use the same lability fractions as for leaf litter + ED_c_to_litr_lab_c(cc,j) = ED_c_to_litr_lab_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_flab(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_cel_c(cc,j) = ED_c_to_litr_cel_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_fcel(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ED_c_to_litr_lig_c(cc,j) = ED_c_to_litr_lig_c(cc,j) + currentpatch%seed_decay(ft) * pftcon%lf_flig(ft) * currentpatch%area/AREA * leaf_prof(cc,j) + ! + enddo + end do + + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + do cc = bounds%begc,bounds%endc + do j = 1, nlevdecomp + ! time unit conversion + ED_c_to_litr_lab_c(cc,j)=ED_c_to_litr_lab_c(cc,j) * mass_convert / time_convert + ED_c_to_litr_cel_c(cc,j)=ED_c_to_litr_cel_c(cc,j) * mass_convert / time_convert + ED_c_to_litr_lig_c(cc,j)=ED_c_to_litr_lig_c(cc,j) * mass_convert / time_convert + + end do + end do + + ! write(iulog,*)'cdk ED_c_to_litr_lab_c: ', ED_c_to_litr_lab_c + ! write(iulog,*)'cdk ED_c_to_litr_cel_c: ', ED_c_to_litr_cel_c + ! write(iulog,*)'cdk ED_c_to_litr_lig_c: ', ED_c_to_litr_lig_c + ! write(iulog,*)'cdk nlevdecomp_full, bounds%begc, bounds%endc: ', nlevdecomp_full, bounds%begc, bounds%endc + ! write(iulog,*)'cdk leaf_prof: ', leaf_prof + ! write(iulog,*)'cdk stem_prof: ', stem_prof + ! write(iulog,*)'cdk froot_prof: ', froot_prof + ! write(iulog,*)'cdk croot_prof_perpatch: ', croot_prof_perpatch + ! write(iulog,*)'cdk croot_prof: ', croot_prof + + end associate + end subroutine flux_into_litter_pools + + !------------------------------------------------------------------------ + + subroutine Summary(this, bounds, num_soilc, filter_soilc, & + ed_allsites_inst, soilbiogeochem_carbonflux_inst, & + soilbiogeochem_carbonstate_inst) + + ! Summarize the combined production and decomposition fluxes into net fluxes + ! Written by Charlie Koven, Feb 2016 + ! + ! !USES: + use ColumnType , only : col + use LandunitType , only : lun + use landunit_varcon , only : istsoil + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_soilc ! number of soil columns in filter + integer , intent(in) :: filter_soilc(:) ! filter for soil columns + type(ed_site_type) , intent(inout), target :: ed_allsites_inst( bounds%begg: ) + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(inout) :: soilbiogeochem_carbonstate_inst + ! + ! !LOCAL VARIABLES: + real(r8) :: npp_hifreq_col(bounds%begc:bounds%endc) ! column-level, high frequency NPP + real(r8) :: dt ! radiation time step (seconds) + integer :: c, g, cc, fc, l + type(ed_site_type), pointer :: cs + type (ed_patch_type) , pointer :: currentPatch + type (ed_cohort_type) , pointer :: currentCohort + integer :: firstsoilpatch(bounds%begg:bounds%endg) ! the first patch in this gridcell that is soil and thus bare... + + associate(& + hr => soilbiogeochem_carbonflux_inst%hr_col, & ! (gC/m2/s) total heterotrophic respiration + totsomc => soilbiogeochem_carbonstate_inst%totsomc_col, & ! (gC/m2) total soil organic matter carbon + totlitc => soilbiogeochem_carbonstate_inst%totlitc_col, & ! (gC/m2) total litter carbon in BGC pools + npp_hifreq => this%npp_hifreq_col, & + nep => this%nep_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & + nbp => this%nbp_col, & + totecosysc => this%totecosysc_col, & + totedc => this%totedc_col, & + totbgcc => this%totbgcc_col, & + biomass_stock => this%biomass_stock_col, & ! total biomass in gC / m2 + ed_litter_stock => this%ed_litter_stock_col, & ! ED litter in gC / m2 + cwd_stock => this%cwd_stock_col, & ! total CWD in gC / m2 + seed_stock => this%seed_stock_col, & ! total seed mass in gC / m2 + ed_to_bgc_this_edts => this%ed_to_bgc_this_edts_col, & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col & + ) + + ! set time steps + dt = real( get_step_size(), r8 ) + + ! zero variables first + do c = bounds%begc,bounds%endc + ! summary flux variables + npp_hifreq(c) = 0._r8 + fire_c_to_atm(c) = 0._r8 + + ! summary stock variables + ed_litter_stock(c) = 0._r8 + cwd_stock(c) = 0._r8 + seed_stock(c) = 0._r8 + biomass_stock(c) = 0._r8 + end do + + ! retrieve the first soil patch associated with each gridcell. + ! make sure we only get the first patch value for places which have soil. + firstsoilpatch(bounds%begg:bounds%endg) = -999 + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + l = col%landunit(c) + if (lun%itype(l) == istsoil .and. col%itype(c) == istsoil) then + firstsoilpatch(g) = col%patchi(c) + endif + enddo + + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + + ! map ed site-level fire fluxes to clm column fluxes + fire_c_to_atm(cc) = ed_allsites_inst(g)%total_burn_flux_to_atm / ( AREA * SHR_CONST_CDAY * 1.e3_r8) + + currentPatch => ed_allsites_inst(g)%oldest_patch + do while(associated(currentPatch)) + + ! map litter, CWD, and seed pools to column level + cwd_stock(cc) = cwd_stock(cc) + (currentPatch%area / AREA) * (sum(currentPatch%cwd_ag)+ & + sum(currentPatch%cwd_bg)) * 1.e3_r8 + ed_litter_stock(cc) = ed_litter_stock(cc) + (currentPatch%area / AREA) * & + (sum(currentPatch%leaf_litter)+sum(currentPatch%root_litter)) * 1.e3_r8 + seed_stock(cc) = seed_stock(cc) + (currentPatch%area / AREA) * sum(currentPatch%seed_bank) * 1.e3_r8 + + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + + ! map ed cohort-level npp fluxes to clm column fluxes + npp_hifreq(cc) = npp_hifreq(cc) + currentCohort%npp_clm * 1e3 * currentCohort%n / ( AREA * dt) + + ! map biomass pools to column level + biomass_stock(cc) = biomass_stock(cc) + (currentCohort%bdead + currentCohort%balive + & + currentCohort%bstore) * currentCohort%n * 1.e3_r8 / AREA + + currentCohort => currentCohort%shorter + enddo !currentCohort + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + + ! calculate NEP and NBP fluxes. + do fc = 1,num_soilc + c = filter_soilc(fc) + nep(c) = npp_hifreq(c) - hr(c) + nbp(c) = npp_hifreq(c) - ( hr(c) + fire_c_to_atm(c) ) + end do + + ! calculate total stocks + do fc = 1,num_soilc + c = filter_soilc(fc) + + totedc(c) = ed_litter_stock(c) + cwd_stock(c) + seed_stock(c) + biomass_stock(c) ! ED stocks + totbgcc(c) = totsomc(c) + totlitc(c) ! BGC stocks + totecosysc(c) = totedc(c) + totbgcc(c) + + end do + + ! in ED timesteps, because of offset between when ED and BGC reconcile the gain and loss of litterfall carbon, + ! (i.e. ED reconciles it instantly, while BGC reconciles it incrementally over the subsequent day) + ! calculate the total ED -> BGC flux and keep track of the last day's info for balance checking purposes + if ( is_beg_curr_day() ) then + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + ed_to_bgc_last_edts(cc) = ed_to_bgc_this_edts(cc) + endif + end do + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + ed_to_bgc_this_edts(cc) = 0._r8 + seed_rain_flux(cc) = 0._r8 + endif + end do + ! + do g = bounds%begg,bounds%endg + if (firstsoilpatch(g) >= 0 .and. ed_allsites_inst(g)%istheresoil) then + cc = ed_allsites_inst(g)%clmcolumn + currentPatch => ed_allsites_inst(g)%oldest_patch + do while(associated(currentPatch)) + ! + ed_to_bgc_this_edts(cc) = ed_to_bgc_this_edts(cc) + (sum(currentPatch%CWD_AG_out) + sum(currentPatch%CWD_BG_out) + & + + sum(currentPatch%seed_decay) + sum(currentPatch%leaf_litter_out) + sum(currentPatch%root_litter_out)) & + * ( currentPatch%area/AREA ) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + ! + seed_rain_flux(cc) = seed_rain_flux(cc) + sum(currentPatch%seed_rain_flux) * 1.e3_r8 / ( 365.0_r8*SHR_CONST_CDAY ) + ! + currentPatch => currentPatch%younger + end do !currentPatch + end if + end do + endif + + end associate + + end subroutine Summary + + + subroutine ED_BGC_Carbon_Balancecheck(this, bounds, num_soilc, filter_soilc, soilbiogeochem_carbonflux_inst) + + ! Integrate in time the fluxes into and out of the ecosystem, and compare these on a daily timestep + ! to the chagne in carbon stocks of the ecosystem + ! + ! Written by Charlie Koven, Feb 2016 + ! + ! !USES: + ! + implicit none + ! + ! !ARGUMENTS + class(ed_clm_type) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_soilc ! number of soil columns in filter + integer , intent(in) :: filter_soilc(:) ! filter for soil columns + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + ! + ! !LOCAL VARIABLES: + real(r8) :: dtime ! land model time step (sec) + integer :: nstep ! model timestep + real(r8) :: nbp_integrated(bounds%begc:bounds%endc) ! total net biome production integrated + real(r8) :: error_total(bounds%begc:bounds%endc) + real(r8) :: error_ed(bounds%begc:bounds%endc) + real(r8) :: error_bgc(bounds%begc:bounds%endc) + real(r8) :: error_tolerance = 1.e-6_r8 + real(r8) :: max_error_ed + real(r8) :: max_error_bgc + real(r8) :: max_error_total + integer :: fc,c + + associate(& + nep => this%nep_col, & + nep_timeintegrated => this%nep_timeintegrated_col, & + hr => soilbiogeochem_carbonflux_inst%hr_col, & + hr_timeintegrated => this%hr_timeintegrated_col, & + npp_hifreq => this%npp_hifreq_col, & + npp_timeintegrated => this%npp_timeintegrated_col, & + fire_c_to_atm => this%fire_c_to_atm_col, & + totecosysc_old => this%totecosysc_old_col, & + totecosysc => this%totecosysc_col, & + totedc_old => this%totedc_old_col, & + totedc => this%totedc_col, & + totbgcc_old => this%totbgcc_old_col, & + totbgcc => this%totbgcc_col, & + ed_to_bgc_this_edts => this%ed_to_bgc_this_edts_col, & + ed_to_bgc_last_edts => this%ed_to_bgc_last_edts_col, & + seed_rain_flux => this%seed_rain_flux_col, & + cbalance_error_ed => this%cbalance_error_ed_col, & + cbalance_error_bgc => this%cbalance_error_bgc_col, & + cbalance_error_total=> this%cbalance_error_total_col & + ) + + dtime = get_step_size() + nstep = get_nstep() + + if (nstep .le. 1) then + ! when starting up the model, initialize the integrator variables + do fc = 1,num_soilc + c = filter_soilc(fc) + totecosysc_old(c) = totecosysc(c) + totedc_old(c) = totedc(c) + totbgcc_old(c) = totbgcc(c) + nep_timeintegrated(c) = 0._r8 + hr_timeintegrated(c) = 0._r8 + npp_timeintegrated(c) = 0._r8 + ! + ! also initialize the ed-BGC flux variables + ed_to_bgc_this_edts(c) = 0._r8 + ed_to_bgc_last_edts(c) = 0._r8 + ! + cbalance_error_ed(c) = 0._r8 + cbalance_error_bgc(c) = 0._r8 + cbalance_error_total(c) = 0._r8 + end do + endif + + if ( .not. is_beg_curr_day() ) then + ! on CLM (half-hourly) timesteps, integrate the NEP fluxes + do fc = 1,num_soilc + c = filter_soilc(fc) + nep_timeintegrated(c) = nep_timeintegrated(c) + nep(c) * dtime + hr_timeintegrated(c) = hr_timeintegrated(c) + hr(c) * dtime + npp_timeintegrated(c) = npp_timeintegrated(c) + npp_hifreq(c) * dtime + end do + else + ! on ED (daily) timesteps, first integrate the NEP fluxes and add in the daily disturbance flux + do fc = 1,num_soilc + c = filter_soilc(fc) + nep_timeintegrated(c) = nep_timeintegrated(c) + nep(c) * dtime + hr_timeintegrated(c) = hr_timeintegrated(c) + hr(c) * dtime + npp_timeintegrated(c) = npp_timeintegrated(c) + npp_hifreq(c) * dtime + nbp_integrated(c) = nep_timeintegrated(c) - fire_c_to_atm(c) * SHR_CONST_CDAY + seed_rain_flux(c)* SHR_CONST_CDAY + end do + + ! next compare the change in carbon and calculate the error + do fc = 1,num_soilc + c = filter_soilc(fc) + error_ed(c) = totedc(c) - totedc_old(c) - (npp_timeintegrated(c) + seed_rain_flux(c)* SHR_CONST_CDAY - ed_to_bgc_this_edts(c)* SHR_CONST_CDAY - fire_c_to_atm(c) * SHR_CONST_CDAY) + error_bgc(c) = totbgcc(c) - totbgcc_old(c) - (ed_to_bgc_last_edts(c)* SHR_CONST_CDAY - hr_timeintegrated(c)) + error_total(c) = totecosysc(c) - totecosysc_old(c) - (nbp_integrated(c) + ed_to_bgc_last_edts(c)* SHR_CONST_CDAY - ed_to_bgc_this_edts(c)* SHR_CONST_CDAY) + end do + ! + ! put in consistent flux units and send to history so we can keep track of the errors + do fc = 1,num_soilc + c = filter_soilc(fc) + cbalance_error_ed(c) = error_ed(c) / SHR_CONST_CDAY + cbalance_error_bgc(c) = error_bgc(c) / SHR_CONST_CDAY + cbalance_error_total(c) = error_total(c) / SHR_CONST_CDAY + end do + + ! for now, rather than crashing the model, lets just report the largest error to see what we're up against + ! + ! RETURN TO THIS LATER AND ADD A CRASHER IF BALANCE EXCEEDS THRESHOLD + ! + ! max_error_total = 0._r8 + ! do fc = 1,num_soilc + ! c = filter_soilc(fc) + ! if (abs(error_total(c)) .gt. max_error_total) then + ! max_error_ed = abs(error_ed(c)) + ! max_error_bgc = abs(error_bgc(c)) + ! max_error_total = abs(error_total(c)) + ! endif + ! end do + ! write(iulog,*) 'ED_BGC_Carbon_Balancecheck: max_error_ed, max_error_bgc, max_error_total (gC / m2 / day): ', max_error_ed, max_error_bgc, max_error_total + + ! reset the C stock and flux integrators + do fc = 1,num_soilc + c = filter_soilc(fc) + totecosysc_old(c) = totecosysc(c) + totedc_old(c) = totedc(c) + totbgcc_old(c) = totbgcc(c) + nep_timeintegrated(c) = 0._r8 + npp_timeintegrated(c) = 0._r8 + hr_timeintegrated(c) = 0._r8 + end do + + endif + + end associate + + end subroutine ED_BGC_Carbon_Balancecheck + end module EDCLMLinkMod diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index ae20235171..cc31ec4ded 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -298,6 +298,7 @@ module EDTypesMod real(r8) :: seed_decay(numpft_ed) ! seed decay in KgC/m2/year real(r8) :: seed_germination(numpft_ed) ! germination rate of seed pool in KgC/m2/year real(r8) :: dseed_dt(numpft_ed) + real(r8) :: seed_rain_flux(numpft_ed) ! flux of seeds from exterior KgC/m2/year (needed for C balance purposes) ! PHOTOSYNTHESIS real(r8) :: psn_z(nclmax,numpft_ed,nlevcan_ed) ! carbon assimilation in each canopy layer, pft, and leaf layer. umolC/m2/s @@ -421,13 +422,14 @@ module EDTypesMod integer :: dleafondate ! doy of leaf on drought:- integer :: dleafoffdate ! doy of leaf on drought:- real(r8) :: water_memory(10) ! last 10 days of soil moisture memory... - real(r8) :: cwd_ag_burned(ncwd) - real(r8) :: leaf_litter_burned(numpft_ed) ! FIRE real(r8) :: acc_ni ! daily nesterov index accumulating over time. real(r8) :: ab ! daily burnt area: m2 real(r8) :: frac_burnt ! fraction of soil burnt in this day. + real(r8) :: total_burn_flux_to_atm ! total carbon burnt to the atmosphere in this day. KgC/site + real(r8) :: cwd_ag_burned(ncwd) + real(r8) :: leaf_litter_burned(numpft_ed) end type ed_site_type diff --git a/components/clm/src/biogeochem/CNCStateUpdate1Mod.F90 b/components/clm/src/biogeochem/CNCStateUpdate1Mod.F90 index 0b6fca3e37..94eb7913d9 100644 --- a/components/clm/src/biogeochem/CNCStateUpdate1Mod.F90 +++ b/components/clm/src/biogeochem/CNCStateUpdate1Mod.F90 @@ -17,6 +17,8 @@ module CNCStateUpdate1Mod use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type use PatchType , only : patch + use EDCLMLinkMod , only : ed_clm_type + use clm_varctl , only : use_ed, use_cn, iulog ! implicit none private @@ -73,6 +75,7 @@ end subroutine CStateUpdate0 !----------------------------------------------------------------------- subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & crop_inst, cnveg_carbonflux_inst, cnveg_carbonstate_inst, & + ed_clm_inst, & soilbiogeochem_carbonflux_inst) ! ! !DESCRIPTION: @@ -88,6 +91,7 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & type(crop_type) , intent(in) :: crop_inst type(cnveg_carbonflux_type) , intent(inout) :: cnveg_carbonflux_inst ! See note below for xsmrpool_to_atm_patch type(cnveg_carbonstate_type) , intent(inout) :: cnveg_carbonstate_inst + type(ed_clm_type) , intent(in) :: ed_clm_inst type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst ! ! !LOCAL VARIABLES: @@ -110,7 +114,8 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & cf_veg => cnveg_carbonflux_inst , & ! Output: cs_veg => cnveg_carbonstate_inst , & ! Output: - cf_soil => soilbiogeochem_carbonflux_inst & ! Output: + cf_soil => soilbiogeochem_carbonflux_inst , & ! Output: + edclm => ed_clm_inst & ) ! set time steps @@ -119,6 +124,7 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & ! Below is the input into the soil biogeochemistry model ! plant to litter fluxes + if (.not. use_ed) then do j = 1,nlevdecomp do fc = 1,num_soilc c = filter_soilc(fc) @@ -133,7 +139,18 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & ( cf_veg%dwt_livecrootc_to_cwdc_col(c,j) + cf_veg%dwt_deadcrootc_to_cwdc_col(c,j) ) *dt end do end do - + else !use_ed + ! here add all ed litterfall and CWD breakdown to litter fluxes + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + cf_soil%decomp_cpools_sourcesink_col(c,j,i_met_lit) = edclm%ED_c_to_litr_lab_c_col(c,j) * dt + cf_soil%decomp_cpools_sourcesink_col(c,j,i_cel_lit) = edclm%ED_c_to_litr_cel_c_col(c,j) * dt + cf_soil%decomp_cpools_sourcesink_col(c,j,i_lig_lit) = edclm%ED_c_to_litr_lig_c_col(c,j) * dt + end do + end do + endif + ! litter and SOM HR fluxes do k = 1, ndecomp_cascade_transitions do j = 1,nlevdecomp @@ -158,6 +175,7 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & end if end do + if (.not. use_ed) then ! seeding fluxes, from dynamic landcover do fc = 1,num_soilc c = filter_soilc(fc) @@ -384,9 +402,10 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & end if end do ! end of patch loop - - end associate - + end if + + end associate + end subroutine CStateUpdate1 end module CNCStateUpdate1Mod diff --git a/components/clm/src/biogeochem/CNDriverMod.F90 b/components/clm/src/biogeochem/CNDriverMod.F90 index e5a684f9e5..3997b208c6 100644 --- a/components/clm/src/biogeochem/CNDriverMod.F90 +++ b/components/clm/src/biogeochem/CNDriverMod.F90 @@ -31,12 +31,12 @@ module CNDriverMod use WaterfluxType , only : waterflux_type use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type - use CanopyStateType , only : canopystate_type use TemperatureType , only : temperature_type use PhotosynthesisMod , only : photosyns_type use ch4Mod , only : ch4_type use EnergyFluxType , only : energyflux_type use SoilHydrologyType , only : soilhydrology_type + use EDCLMLinkMod , only : ed_clm_type ! ! !PUBLIC TYPES: implicit none @@ -83,6 +83,7 @@ subroutine CNDriverNoLeaching(bounds, c14_cnveg_carbonflux_inst, c14_cnveg_carbonstate_inst, & cnveg_nitrogenflux_inst, cnveg_nitrogenstate_inst, & c_products_inst, c13_products_inst, c14_products_inst, n_products_inst, & + ed_clm_inst, & soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -156,6 +157,7 @@ subroutine CNDriverNoLeaching(bounds, type(cn_products_type) , intent(inout) :: c14_products_inst type(cn_products_type) , intent(inout) :: n_products_inst type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst + type(ed_clm_type) , intent(in) :: ed_clm_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 @@ -539,15 +541,18 @@ subroutine CNDriverNoLeaching(bounds, ! Update all prognostic carbon state variables (except for gap-phase mortality and fire fluxes) call CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & crop_inst, cnveg_carbonflux_inst, cnveg_carbonstate_inst, & + ed_clm_inst, & soilbiogeochem_carbonflux_inst) if ( use_c13 ) then call CStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp, & crop_inst, c13_cnveg_carbonflux_inst, c13_cnveg_carbonstate_inst, & + ed_clm_inst, & c13_soilbiogeochem_carbonflux_inst) end if if ( use_c14 ) then call CStateUpdate1(num_soilc, filter_soilc, num_soilp, filter_soilp, & crop_inst, c14_cnveg_carbonflux_inst, c14_cnveg_carbonstate_inst, & + ed_clm_inst, & c14_soilbiogeochem_carbonflux_inst) end if diff --git a/components/clm/src/biogeochem/CNVegetationFacade.F90 b/components/clm/src/biogeochem/CNVegetationFacade.F90 index 740a2ede59..c0b81612c6 100644 --- a/components/clm/src/biogeochem/CNVegetationFacade.F90 +++ b/components/clm/src/biogeochem/CNVegetationFacade.F90 @@ -92,13 +92,17 @@ module CNVegetationFacade ! !PUBLIC TYPES: type, public :: cn_vegetation_type - private - + ! FIXME(bja, 2016-06) These need to be public for use when ED is + ! turned on. Should either be moved out of here or create some ED + ! version of the facade.... type(cnveg_state_type) :: cnveg_state_inst type(cnveg_carbonstate_type) :: cnveg_carbonstate_inst + type(cnveg_carbonflux_type) :: cnveg_carbonflux_inst + + !X!private + type(cnveg_carbonstate_type) :: c13_cnveg_carbonstate_inst type(cnveg_carbonstate_type) :: c14_cnveg_carbonstate_inst - type(cnveg_carbonflux_type) :: cnveg_carbonflux_inst type(cnveg_carbonflux_type) :: c13_cnveg_carbonflux_inst type(cnveg_carbonflux_type) :: c14_cnveg_carbonflux_inst type(cnveg_nitrogenstate_type) :: cnveg_nitrogenstate_inst @@ -573,7 +577,7 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, & num_soilc, filter_soilc, & num_soilp, filter_soilp, & num_pcropp, filter_pcropp, & - doalb, & + doalb, ed_clm_inst, & soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -590,6 +594,8 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, & ! Should only be called if use_cn is true ! ! !USES: + use EDCLMLinkMod , only : ed_clm_type + ! ! !ARGUMENTS: class(cn_vegetation_type) , intent(inout) :: this @@ -601,6 +607,7 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, & 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(ed_clm_type) , intent(in) :: 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 @@ -640,6 +647,7 @@ subroutine EcosystemDynamicsPreDrainage(this, bounds, & this%cnveg_nitrogenflux_inst, this%cnveg_nitrogenstate_inst, & this%c_products_inst, this%c13_products_inst, this%c14_products_inst, & this%n_products_inst, & + ed_clm_inst, & soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & diff --git a/components/clm/src/biogeochem/EDBGCDynMod.F90 b/components/clm/src/biogeochem/EDBGCDynMod.F90 new file mode 100644 index 0000000000..02f04a9b0d --- /dev/null +++ b/components/clm/src/biogeochem/EDBGCDynMod.F90 @@ -0,0 +1,371 @@ +module EDBGCDynMod + +! This module creates a pathway to call the belowground biogeochemistry code as driven by the ED vegetation model +! but bypassing the aboveground CN vegetation code. It is modeled after the CNDriverMod in its call sequence and +! functionality. + + 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_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: 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 CNCIsoFluxMod , only: CIsoFlux1, CIsoFlux2, CIsoFlux2h, CIsoFlux3 + use CNC14DecayMod , only: C14Decay + 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_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, & + crop_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 diff --git a/components/clm/src/biogeophys/BalanceCheckMod.F90 b/components/clm/src/biogeophys/BalanceCheckMod.F90 index 58767d4312..09140a6c1d 100644 --- a/components/clm/src/biogeophys/BalanceCheckMod.F90 +++ b/components/clm/src/biogeophys/BalanceCheckMod.F90 @@ -568,7 +568,7 @@ subroutine BalanceCheck( bounds, num_do_smb_c, filter_do_smb_c, & found = .false. do p = bounds%begp, bounds%endp if (patch%active(p)) then - if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > 1.e-7_r8) ) then + if ( (errsol(p) /= spval) .and. (abs(errsol(p)) > 1.e-6_r8) ) then found = .true. indexp = p indexg = patch%gridcell(indexp) diff --git a/components/clm/src/biogeophys/CanopyStateType.F90 b/components/clm/src/biogeophys/CanopyStateType.F90 index 4a295010ce..91abfa9d9b 100644 --- a/components/clm/src/biogeophys/CanopyStateType.F90 +++ b/components/clm/src/biogeophys/CanopyStateType.F90 @@ -593,7 +593,7 @@ subroutine Restart(this, bounds, ncid, flag) end do end if - if (use_cn) then + if (use_cn .or. use_ed) then call restartvar(ncid=ncid, flag=flag, varname='altmax', xtype=ncd_double, & dim1name='column', long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%altmax_col) diff --git a/components/clm/src/biogeophys/HydrologyNoDrainageMod.F90 b/components/clm/src/biogeophys/HydrologyNoDrainageMod.F90 index 46925e181a..3006f05782 100644 --- a/components/clm/src/biogeophys/HydrologyNoDrainageMod.F90 +++ b/components/clm/src/biogeophys/HydrologyNoDrainageMod.F90 @@ -370,7 +370,7 @@ subroutine HydrologyNoDrainage(bounds, & end do end do - if (use_cn) then +! if (use_cn) then ! Update soilpsi. ! ZMS: Note this could be merged with the following loop updating smp_l in the future. do j = 1, nlevgrnd @@ -394,7 +394,7 @@ subroutine HydrologyNoDrainage(bounds, & end if end do end do - end if +! end if ! Update smp_l for history and for ch4Mod. ! ZMS: Note, this form, which seems to be the same as used in SoilWater, DOES NOT distinguish between @@ -412,7 +412,7 @@ subroutine HydrologyNoDrainage(bounds, & end do end do - if (use_cn) then + ! if (use_cn) then ! Available soil water up to a depth of 0.05 m. ! Potentially available soil water (=whc) up to a depth of 0.05 m. ! Water content as fraction of whc up to a depth of 0.05 m. @@ -474,7 +474,7 @@ subroutine HydrologyNoDrainage(bounds, & end if wf2(c) = tsw/stsw end do - end if + ! end if ! top-layer diagnostics do fc = 1, num_snowc diff --git a/components/clm/src/cpl/lnd_comp_mct.F90 b/components/clm/src/cpl/lnd_comp_mct.F90 index 38b6aa8ca7..2207f0276d 100644 --- a/components/clm/src/cpl/lnd_comp_mct.F90 +++ b/components/clm/src/cpl/lnd_comp_mct.F90 @@ -347,7 +347,6 @@ subroutine lnd_run_mct(EClock, cdata_l, x2l_l, l2x_l) !--------------------------------------------------------------------------- ! Determine processor bounds - write(iulog,*) 'in run_mct' call get_proc_bounds(bounds) #if (defined _MEMTRACE) diff --git a/components/clm/src/main/clm_driver.F90 b/components/clm/src/main/clm_driver.F90 index c789fedbfd..5d63a4803c 100644 --- a/components/clm/src/main/clm_driver.F90 +++ b/components/clm/src/main/clm_driver.F90 @@ -77,6 +77,7 @@ module clm_driver use PatchType , only : patch use clm_instMod use clm_initializeMod , only : soil_water_retention_curve + use EDBGCDynMod , only : EDBGCDyn, EDBGCDynSummary ! ! !PUBLIC TYPES: implicit none @@ -687,6 +688,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) filter(nc)%num_soilc, filter(nc)%soilc, & filter(nc)%num_soilp, filter(nc)%soilp, & filter(nc)%num_pcropp, filter(nc)%pcropp, doalb, & + ed_clm_inst, & soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -699,6 +701,32 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) call t_stopf('ecosysdyn') + elseif (use_ed) then + + ! call EDBGCDyn(bounds_clump, & + ! filter(nc)%num_soilc, filter(nc)%soilc, & + ! ed_clm_inst, & + ! soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst, & + ! c13_soilbiogeochem_carbonflux_inst, c13_soilbiogeochem_carbonstate_inst, & + ! c14_soilbiogeochem_carbonflux_inst, c14_soilbiogeochem_carbonstate_inst, & + ! soilbiogeochem_state_inst, & + ! soilbiogeochem_nitrogenflux_inst, soilbiogeochem_nitrogenstate_inst) + + call EDBGCDyn(bounds_clump, & + filter(nc)%num_soilc, filter(nc)%soilc, & + filter(nc)%num_soilp, filter(nc)%soilp, & + filter(nc)%num_pcropp, filter(nc)%pcropp, doalb, & + bgc_vegetation_inst%cnveg_carbonflux_inst, & + bgc_vegetation_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) + end if ! Prescribed biogeography - prescribed canopy structure, some prognostic carbon fluxes @@ -758,6 +786,17 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate) end if + if ( use_ed ) then + call EDBGCDynSummary(bounds_clump, & + filter(nc)%num_soilc, filter(nc)%soilc, & + filter(nc)%num_soilp, filter(nc)%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(bounds_clump%begg:bounds_clump%endg)) + end if + ! ============================================================================ ! Check the energy and water balance and also carbon and nitrogen balance ! ============================================================================ diff --git a/components/clm/src/main/clm_instMod.F90 b/components/clm/src/main/clm_instMod.F90 index 32df96eccb..22f7b00168 100644 --- a/components/clm/src/main/clm_instMod.F90 +++ b/components/clm/src/main/clm_instMod.F90 @@ -325,7 +325,7 @@ subroutine clm_instInit(bounds) call drydepvel_inst%Init(bounds) - if (use_cn) then + if (use_cn .or. use_ed ) then ! Initialize soilbiogeochem_state_inst @@ -342,7 +342,7 @@ subroutine clm_instInit(bounds) call init_decompcascade_cn(bounds, soilbiogeochem_state_inst) end if - ! Initalize soilbiogeochem carbon and nitrogen types + ! Initalize soilbiogeochem carbon types call soilbiogeochem_carbonstate_inst%Init(bounds, carbon_type='c12', ratio=1._r8) if (use_c13) then @@ -353,10 +353,6 @@ subroutine clm_instInit(bounds) call c14_soilbiogeochem_carbonstate_inst%Init(bounds, carbon_type='c14', ratio=c14ratio, & c12_soilbiogeochem_carbonstate_inst=soilbiogeochem_carbonstate_inst) end if - call soilbiogeochem_nitrogenstate_inst%Init(bounds, & - soilbiogeochem_carbonstate_inst%decomp_cpools_vr_col(begc:endc,1:nlevdecomp_full,1:ndecomp_pools), & - soilbiogeochem_carbonstate_inst%decomp_cpools_col(begc:endc,1:ndecomp_pools), & - soilbiogeochem_carbonstate_inst%decomp_cpools_1m_col(begc:endc, 1:ndecomp_pools)) call soilbiogeochem_carbonflux_inst%Init(bounds, carbon_type='c12') if (use_c13) then @@ -365,6 +361,18 @@ subroutine clm_instInit(bounds) if (use_c14) then call c14_soilbiogeochem_carbonflux_inst%Init(bounds, carbon_type='c14') end if + + end if + + if ( use_cn ) then + + ! Initalize soilbiogeochem nitrogen types + + call soilbiogeochem_nitrogenstate_inst%Init(bounds, & + soilbiogeochem_carbonstate_inst%decomp_cpools_vr_col(begc:endc,1:nlevdecomp_full,1:ndecomp_pools), & + soilbiogeochem_carbonstate_inst%decomp_cpools_col(begc:endc,1:ndecomp_pools), & + soilbiogeochem_carbonstate_inst%decomp_cpools_1m_col(begc:endc, 1:ndecomp_pools)) + call soilbiogeochem_nitrogenflux_inst%Init(bounds) end if ! end of if use_cn @@ -488,7 +496,7 @@ subroutine clm_instRest(bounds, ncid, flag) call ch4_inst%restart(bounds, ncid, flag=flag) end if - if (use_cn) then + if (use_cn .or. use_ed) then call soilbiogeochem_state_inst%restart(bounds, ncid, flag=flag) call soilbiogeochem_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c12') @@ -501,6 +509,9 @@ subroutine clm_instRest(bounds, ncid, flag) c12_soilbiogeochem_carbonstate_inst=soilbiogeochem_carbonstate_inst) end if call soilbiogeochem_carbonflux_inst%restart(bounds, ncid, flag=flag) + endif + if ( use_cn ) then + call soilbiogeochem_nitrogenstate_inst%restart(bounds, ncid, flag=flag) call soilbiogeochem_nitrogenflux_inst%restart(bounds, ncid, flag=flag) @@ -513,6 +524,7 @@ subroutine clm_instRest(bounds, ncid, flag) call ED_Phenology_inst%restart(bounds, ncid, flag=flag) call EDRest ( bounds, ncid, flag, ed_allsites_inst(bounds%begg:bounds%endg), & ed_clm_inst, ed_phenology_inst, waterstate_inst, canopystate_inst ) + call ed_clm_inst%Restart(bounds, ncid, flag=flag) end if end subroutine clm_instRest diff --git a/components/clm/src/main/clm_varpar.F90 b/components/clm/src/main/clm_varpar.F90 index 863c515f4b..af9b810923 100644 --- a/components/clm/src/main/clm_varpar.F90 +++ b/components/clm/src/main/clm_varpar.F90 @@ -11,6 +11,8 @@ module clm_varpar use clm_varctl , only: use_century_decomp, use_c13, use_c14 use clm_varctl , only: iulog, use_crop, create_crop_landunit, irrigate use clm_varctl , only: use_vichydro, soil_layerstruct + use clm_varctl , only: use_ed + ! ! !PUBLIC TYPES: implicit none @@ -188,21 +190,39 @@ subroutine clm_varpar_init() write(iulog, *) end if - if (use_century_decomp) then - ndecomp_pools = 7 - ndecomp_cascade_transitions = 10 - i_met_lit = 1 - i_cel_lit = 2 - i_lig_lit = 3 - i_cwd = 4 + if ( use_ed ) then + if (use_century_decomp) then + ndecomp_pools = 6 + ndecomp_cascade_transitions = 8 + i_met_lit = 1 + i_cel_lit = 2 + i_lig_lit = 3 + i_cwd = 0 + else + ndecomp_pools = 7 + ndecomp_cascade_transitions = 7 + i_met_lit = 1 + i_cel_lit = 2 + i_lig_lit = 3 + i_cwd = 0 + end if else - ndecomp_pools = 8 - ndecomp_cascade_transitions = 9 - i_met_lit = 1 - i_cel_lit = 2 - i_lig_lit = 3 - i_cwd = 4 - end if + if (use_century_decomp) then + ndecomp_pools = 7 + ndecomp_cascade_transitions = 10 + i_met_lit = 1 + i_cel_lit = 2 + i_lig_lit = 3 + i_cwd = 4 + else + ndecomp_pools = 8 + ndecomp_cascade_transitions = 9 + i_met_lit = 1 + i_cel_lit = 2 + i_lig_lit = 3 + i_cwd = 4 + end if + endif end subroutine clm_varpar_init diff --git a/components/clm/src/main/readParamsMod.F90 b/components/clm/src/main/readParamsMod.F90 index 5deffa74af..fa62a49be7 100644 --- a/components/clm/src/main/readParamsMod.F90 +++ b/components/clm/src/main/readParamsMod.F90 @@ -71,8 +71,11 @@ subroutine readParameters (nutrient_competition_method) call SFParamsRead(ncid) end if + if (use_ed .or. use_cn) then + call CNParamsReadShared(ncid, NLFilename_in) ! this is called CN params but really is for the soil biogeochem parameters + end if + if (use_cn) then - call CNParamsReadShared(ncid, NLFilename_in) call nutrient_competition_method%readParams(ncid) call readCNGapMortParams(ncid) call readCNMRespParams(ncid) diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 index 17625b0e97..0ab2b7e7b0 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 @@ -7,7 +7,7 @@ module SoilBiogeochemCarbonStateType use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools, nlevcan use clm_varpar , only : nlevdecomp_full, nlevdecomp use clm_varcon , only : spval, ispval, dzsoi_decomp, zisoi, zsoi, c3_r2 - use clm_varctl , only : iulog, use_vertsoilc, spinup_state + use clm_varctl , only : iulog, use_vertsoilc, spinup_state, use_ed use landunit_varcon , only : istcrop, istsoil use abortutils , only : endrun use spmdMod , only : masterproc @@ -96,7 +96,9 @@ subroutine InitAllocate(this, bounds) this%decomp_cpools_vr_col(:,:,:)= nan allocate(this%ctrunc_col (begc :endc)) ; this%ctrunc_col (:) = nan - allocate(this%cwdc_col (begc :endc)) ; this%cwdc_col (:) = nan + if ( .not. use_ed ) then + allocate(this%cwdc_col (begc :endc)) ; this%cwdc_col (:) = nan + endif allocate(this%totlitc_col (begc :endc)) ; this%totlitc_col (:) = nan allocate(this%totsomc_col (begc :endc)) ; this%totsomc_col (:) = nan allocate(this%totlitc_1m_col (begc :endc)) ; this%totlitc_1m_col (:) = nan @@ -409,19 +411,20 @@ subroutine InitCold(this, bounds, ratio, c12_soilbiogeochem_carbonstate_inst) endif end if - if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then - if (present(c12_soilbiogeochem_carbonstate_inst)) then - this%cwdc_col(c) = c12_soilbiogeochem_carbonstate_inst%cwdc_col(c) * ratio - else - this%cwdc_col(c) = 0._r8 + if ( .not. use_ed ) then + if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then + if (present(c12_soilbiogeochem_carbonstate_inst)) then + this%cwdc_col(c) = c12_soilbiogeochem_carbonstate_inst%cwdc_col(c) * ratio + else + this%cwdc_col(c) = 0._r8 + end if + this%ctrunc_col(c) = 0._r8 + this%totlitc_col(c) = 0._r8 + this%totsomc_col(c) = 0._r8 + this%totlitc_1m_col(c) = 0._r8 + this%totsomc_1m_col(c) = 0._r8 end if - this%ctrunc_col(c) = 0._r8 - this%totlitc_col(c) = 0._r8 - this%totsomc_col(c) = 0._r8 - this%totlitc_1m_col(c) = 0._r8 - this%totsomc_1m_col(c) = 0._r8 end if - end do ! now loop through special filters and explicitly set the variables that @@ -740,7 +743,9 @@ subroutine SetValues ( this, num_column, filter_column, value_column) do fi = 1,num_column i = filter_column(fi) - this%cwdc_col(i) = value_column + if ( .not. use_ed ) then + this%cwdc_col(i) = value_column + end if this%ctrunc_col(i) = value_column this%totlitc_col(i) = value_column this%totlitc_1m_col(i) = value_column @@ -918,19 +923,22 @@ subroutine Summary(this, bounds, num_allc, filter_allc) end do ! coarse woody debris carbon - do fc = 1,num_allc - c = filter_allc(fc) - this%cwdc_col(c) = 0._r8 - end do - do l = 1, ndecomp_pools - if ( decomp_cascade_con%is_cwd(l) ) then - do fc = 1,num_allc - c = filter_allc(fc) - this%cwdc_col(c) = this%cwdc_col(c) + this%decomp_cpools_col(c,l) - end do - end if - end do + if (.not. use_ed ) then + do fc = 1,num_allc + c = filter_allc(fc) + this%cwdc_col(c) = 0._r8 + end do + do l = 1, ndecomp_pools + if ( decomp_cascade_con%is_cwd(l) ) then + do fc = 1,num_allc + c = filter_allc(fc) + this%cwdc_col(c) = this%cwdc_col(c) + this%decomp_cpools_col(c,l) + end do + end if + end do + end if + end subroutine Summary !----------------------------------------------------------------------- diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 index e65fc0127d..555c858877 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 @@ -11,7 +11,7 @@ module SoilBiogeochemDecompCascadeBGCMod use shr_log_mod , only : errMsg => shr_log_errMsg use clm_varpar , only : nlevsoi, nlevgrnd, nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools use clm_varpar , only : i_met_lit, i_cel_lit, i_lig_lit, i_cwd - use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_vertsoilc + use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_vertsoilc, use_ed use clm_varcon , only : zsoi use decompMod , only : bounds_type use spmdMod , only : masterproc @@ -27,6 +27,7 @@ module SoilBiogeochemDecompCascadeBGCMod use ColumnType , only : col use GridcellType , only : grc use SoilBiogeochemStateType , only : get_spinup_latitude_term + use EDCLMLinkMod , only : cwd_fcel_ed, cwd_flig_ed ! implicit none private @@ -204,6 +205,11 @@ subroutine readParams ( ncid ) call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) params_inst%cwd_flig_bgc=tempr + + if ( use_ed ) then + cwd_fcel_ed = params_inst%cwd_fcel_bgc + cwd_flig_ed = params_inst%cwd_flig_bgc + endif end subroutine readParams @@ -383,22 +389,28 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_cellulose(i_litr3) = .false. is_lignin(i_litr3) = .true. - ! CWD - floating_cn_ratio_decomp_pools(i_cwd) = .true. - decomp_pool_name_restart(i_cwd) = 'cwd' - decomp_pool_name_history(i_cwd) = 'CWD' - decomp_pool_name_long(i_cwd) = 'coarse woody debris' - decomp_pool_name_short(i_cwd) = 'CWD' - is_litter(i_cwd) = .false. - is_soil(i_cwd) = .false. - is_cwd(i_cwd) = .true. - initial_cn_ratio(i_cwd) = 90._r8 - initial_stock(i_cwd) = 0._r8 - is_metabolic(i_cwd) = .false. - is_cellulose(i_cwd) = .false. - is_lignin(i_cwd) = .false. - - i_soil1 = 5 + if (.not. use_ed) then + ! CWD + floating_cn_ratio_decomp_pools(i_cwd) = .true. + decomp_pool_name_restart(i_cwd) = 'cwd' + decomp_pool_name_history(i_cwd) = 'CWD' + decomp_pool_name_long(i_cwd) = 'coarse woody debris' + decomp_pool_name_short(i_cwd) = 'CWD' + is_litter(i_cwd) = .false. + is_soil(i_cwd) = .false. + is_cwd(i_cwd) = .true. + initial_cn_ratio(i_cwd) = 90._r8 + initial_stock(i_cwd) = 0._r8 + is_metabolic(i_cwd) = .false. + is_cellulose(i_cwd) = .false. + is_lignin(i_cwd) = .false. + endif + + if (.not. use_ed) then + i_soil1 = 5 + else + i_soil1 = 4 + endif floating_cn_ratio_decomp_pools(i_soil1) = .false. decomp_pool_name_restart(i_soil1) = 'soil1' decomp_pool_name_history(i_soil1) = 'SOIL1' @@ -413,7 +425,11 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_cellulose(i_soil1) = .false. is_lignin(i_soil1) = .false. - i_soil2 = 6 + if (.not. use_ed) then + i_soil2 = 6 + else + i_soil2 = 5 + endif floating_cn_ratio_decomp_pools(i_soil2) = .false. decomp_pool_name_restart(i_soil2) = 'soil2' decomp_pool_name_history(i_soil2) = 'SOIL2' @@ -428,7 +444,11 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_cellulose(i_soil2) = .false. is_lignin(i_soil2) = .false. - i_soil3 = 7 + if (.not. use_ed) then + i_soil3 = 7 + else + i_soil3 = 6 + endif floating_cn_ratio_decomp_pools(i_soil3) = .false. decomp_pool_name_restart(i_soil3) = 'soil3' decomp_pool_name_history(i_soil3) = 'SOIL3' @@ -444,12 +464,6 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_lignin(i_soil3) = .false. - i_litr1 = i_met_lit - i_litr2 = i_cel_lit - i_litr3 = i_lig_lit - i_soil1 = 5 - i_soil2 = 6 - i_soil3 = 7 speedup_fac = 1._r8 !lit1 @@ -458,7 +472,9 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i spinup_factor(i_litr2) = 1._r8 spinup_factor(i_litr3) = 1._r8 !CWD - spinup_factor(i_cwd) = max(1._r8, (speedup_fac * params_inst%tau_cwd_bgc / 2._r8 )) + if (.not. use_ed) then + spinup_factor(i_cwd) = max(1._r8, (speedup_fac * params_inst%tau_cwd_bgc / 2._r8 )) + end if !som1 spinup_factor(i_soil1) = 1._r8 !som2,3 @@ -470,7 +486,6 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i write(iulog,*) 'Spinup factors ',spinup_factor end if - !---------------- list of transitions and their time-independent coefficients ---------------! i_l1s1 = 1 cascade_step_name(i_l1s1) = 'L1S1' @@ -528,19 +543,21 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i cascade_receiver_pool(i_s3s1) = i_soil1 pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = 1.0_r8 - i_cwdl2 = 9 - cascade_step_name(i_cwdl2) = 'CWDL2' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = rf_cwdl2 - cascade_donor_pool(i_cwdl2) = i_cwd - cascade_receiver_pool(i_cwdl2) = i_litr2 - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel - - i_cwdl3 = 10 - cascade_step_name(i_cwdl3) = 'CWDL3' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = rf_cwdl3 - cascade_donor_pool(i_cwdl3) = i_cwd - cascade_receiver_pool(i_cwdl3) = i_litr3 - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig + if (.not. use_ed) then + i_cwdl2 = 9 + cascade_step_name(i_cwdl2) = 'CWDL2' + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = rf_cwdl2 + cascade_donor_pool(i_cwdl2) = i_cwd + cascade_receiver_pool(i_cwdl2) = i_litr2 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel + + i_cwdl3 = 10 + cascade_step_name(i_cwdl3) = 'CWDL3' + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = rf_cwdl3 + cascade_donor_pool(i_cwdl3) = i_cwd + cascade_receiver_pool(i_cwdl3) = i_litr3 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig + end if deallocate(rf_s1s2) deallocate(rf_s1s3) @@ -695,9 +712,18 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & i_litr1 = i_met_lit i_litr2 = i_cel_lit i_litr3 = i_lig_lit - i_soil1 = 5 - i_soil2 = 6 - i_soil3 = 7 + if ( use_ed ) then + i_soil1 = 4 + i_soil2 = 5 + i_soil3 = 6 + else + i_soil1 = 5 + i_soil2 = 6 + i_soil3 = 7 + endif + + ! calc ref rate + catanf_30 = catanf(30._r8) if ( spinup_state >= 1 ) then do fc = 1,num_soilc @@ -752,9 +778,6 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & end do endif - ! calc ref rate - catanf_30 = catanf(30._r8) - !--- time dependent coefficients-----! if ( nlevdecomp .eq. 1 ) then @@ -989,6 +1012,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & end do end if + ! calculate rate constants for all litter and som pools if (use_vertsoilc) then do j = 1,nlevdecomp do fc = 1,num_soilc @@ -999,8 +1023,6 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & * spinup_geogterm_l23(c) decomp_k(c,j,i_litr3) = k_l2_l3 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) & * spinup_geogterm_l23(c) - decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) & - * spinup_geogterm_cwd(c) decomp_k(c,j,i_soil1) = k_s1 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) & * spinup_geogterm_s1(c) decomp_k(c,j,i_soil2) = k_s2 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) & @@ -1016,7 +1038,6 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & decomp_k(c,j,i_litr1) = k_l1 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_l1(c) decomp_k(c,j,i_litr2) = k_l2_l3 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_l23(c) decomp_k(c,j,i_litr3) = k_l2_l3 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_l23(c) - decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_s1(c) decomp_k(c,j,i_soil1) = k_s1 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_s1(c) decomp_k(c,j,i_soil2) = k_s2 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_s2(c) decomp_k(c,j,i_soil3) = k_s3 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) * spinup_geogterm_s3(c) @@ -1024,6 +1045,25 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & end do end if + ! do the same for cwd, but only if ed is not enabled (because ED handles CWD on its own structure + if (.not. use_ed) then + if (use_vertsoilc) then + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) + end do + end do + else + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) + end do + end do + end if + end if + end associate end subroutine decomp_rate_constants_bgc diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 index 26fb3149d5..a109f05295 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 @@ -11,7 +11,7 @@ module SoilBiogeochemDecompCascadeCNMod use shr_log_mod , only : errMsg => shr_log_errMsg use clm_varpar , only : nlevsoi, nlevgrnd, nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools use clm_varpar , only : i_met_lit, i_cel_lit, i_lig_lit, i_cwd - use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_vertsoilc + use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_vertsoilc, use_ed use clm_varcon , only : zsoi use decompMod , only : bounds_type use abortutils , only : endrun @@ -24,6 +24,7 @@ module SoilBiogeochemDecompCascadeCNMod use TemperatureType , only : temperature_type use ch4Mod , only : ch4_type use ColumnType , only : col + use EDCLMLinkMod , only : cwd_fcel_ed, cwd_flig_ed ! implicit none private @@ -205,6 +206,11 @@ subroutine readParams ( ncid ) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(__FILE__, __LINE__)) params_inst%cwd_flig_cn=tempr + if ( use_ed ) then + cwd_fcel_ed = params_inst%cwd_fcel_cn + cwd_flig_ed = params_inst%cwd_flig_cn + endif + end subroutine readParams !----------------------------------------------------------------------- @@ -342,21 +348,27 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) is_cellulose(i_litr3) = .false. is_lignin(i_litr3) = .true. - floating_cn_ratio_decomp_pools(i_cwd) = .true. - decomp_pool_name_restart(i_cwd) = 'cwd' - decomp_pool_name_history(i_cwd) = 'CWD' - decomp_pool_name_long(i_cwd) = 'coarse woody debris' - decomp_pool_name_short(i_cwd) = 'CWD' - is_litter(i_cwd) = .false. - is_soil(i_cwd) = .false. - is_cwd(i_cwd) = .true. - initial_cn_ratio(i_cwd) = 500._r8 - initial_stock(i_cwd) = 0._r8 - is_metabolic(i_cwd) = .false. - is_cellulose(i_cwd) = .false. - is_lignin(i_cwd) = .false. - - i_soil1 = 5 + if (.not. use_ed) then + floating_cn_ratio_decomp_pools(i_cwd) = .true. + decomp_pool_name_restart(i_cwd) = 'cwd' + decomp_pool_name_history(i_cwd) = 'CWD' + decomp_pool_name_long(i_cwd) = 'coarse woody debris' + decomp_pool_name_short(i_cwd) = 'CWD' + is_litter(i_cwd) = .false. + is_soil(i_cwd) = .false. + is_cwd(i_cwd) = .true. + initial_cn_ratio(i_cwd) = 500._r8 + initial_stock(i_cwd) = 0._r8 + is_metabolic(i_cwd) = .false. + is_cellulose(i_cwd) = .false. + is_lignin(i_cwd) = .false. + end if + + if ( .not. use_ed ) then + i_soil1 = 5 + else + i_soil1 = 4 + endif floating_cn_ratio_decomp_pools(i_soil1) = .false. decomp_pool_name_restart(i_soil1) = 'soil1' decomp_pool_name_history(i_soil1) = 'SOIL1' @@ -371,7 +383,11 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) is_cellulose(i_soil1) = .false. is_lignin(i_soil1) = .false. - i_soil2 = 6 + if ( .not. use_ed ) then + i_soil2 = 6 + else + i_soil2 = 5 + endif floating_cn_ratio_decomp_pools(i_soil2) = .false. decomp_pool_name_restart(i_soil2) = 'soil2' decomp_pool_name_history(i_soil2) = 'SOIL2' @@ -386,7 +402,11 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) is_cellulose(i_soil2) = .false. is_lignin(i_soil2) = .false. - i_soil3 = 7 + if ( .not. use_ed ) then + i_soil3 = 7 + else + i_soil3 = 6 + endif floating_cn_ratio_decomp_pools(i_soil3) = .false. decomp_pool_name_restart(i_soil3) = 'soil3' decomp_pool_name_history(i_soil3) = 'SOIL3' @@ -401,7 +421,11 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) is_cellulose(i_soil3) = .false. is_lignin(i_soil3) = .false. - i_soil4 = 8 + if ( .not. use_ed ) then + i_soil4 = 8 + else + i_soil4 = 7 + endif floating_cn_ratio_decomp_pools(i_soil4) = .false. decomp_pool_name_restart(i_soil4) = 'soil4' decomp_pool_name_history(i_soil4) = 'SOIL4' @@ -435,7 +459,9 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) spinup_factor(i_litr1) = 1._r8 spinup_factor(i_litr2) = 1._r8 spinup_factor(i_litr3) = 1._r8 - spinup_factor(i_cwd) = 1._r8 + if (.not. use_ed) then + spinup_factor(i_cwd) = 1._r8 + end if spinup_factor(i_soil1) = params_inst%spinup_vector(1) spinup_factor(i_soil2) = params_inst%spinup_vector(2) spinup_factor(i_soil3) = params_inst%spinup_vector(3) @@ -492,20 +518,21 @@ subroutine init_decompcascade_cn(bounds, soilbiogeochem_state_inst) cascade_receiver_pool(i_s4atm) = i_atm pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s4atm) = 1.0_r8 - i_cwdl2 = 8 - cascade_step_name(i_cwdl2) = 'CWDL2' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = 0._r8 - cascade_donor_pool(i_cwdl2) = i_cwd - cascade_receiver_pool(i_cwdl2) = i_litr2 - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel - - i_cwdl3 = 9 - cascade_step_name(i_cwdl3) = 'CWDL3' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = 0._r8 - cascade_donor_pool(i_cwdl3) = i_cwd - cascade_receiver_pool(i_cwdl3) = i_litr3 - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig - + if (.not. use_ed) then + i_cwdl2 = 8 + cascade_step_name(i_cwdl2) = 'CWDL2' + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = 0._r8 + cascade_donor_pool(i_cwdl2) = i_cwd + cascade_receiver_pool(i_cwdl2) = i_litr2 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel + + i_cwdl3 = 9 + cascade_step_name(i_cwdl3) = 'CWDL3' + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = 0._r8 + cascade_donor_pool(i_cwdl3) = i_cwd + cascade_receiver_pool(i_cwdl3) = i_litr3 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig + end if end associate @@ -658,11 +685,17 @@ subroutine decomp_rate_constants_cn(bounds, & i_litr1 = 1 i_litr2 = 2 i_litr3 = 3 - i_soil1 = 5 - i_soil2 = 6 - i_soil3 = 7 - i_soil4 = 8 - + if (use_ed) then + i_soil1 = 4 + i_soil2 = 5 + i_soil3 = 6 + i_soil4 = 7 + else + i_soil1 = 5 + i_soil2 = 6 + i_soil3 = 7 + i_soil4 = 8 + endif !--- time dependent coefficients-----! if ( nlevdecomp .eq. 1 ) then @@ -867,6 +900,7 @@ subroutine decomp_rate_constants_cn(bounds, & end do end if + ! calculate rate constants for all litter and som pools if (use_vertsoilc) then do j = 1,nlevdecomp do fc = 1,num_soilc @@ -874,7 +908,6 @@ subroutine decomp_rate_constants_cn(bounds, & decomp_k(c,j,i_litr1) = k_l1 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_litr2) = k_l2 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_litr3) = k_l3 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt - decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil1) = k_s1 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil2) = k_s2 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil3) = k_s3 * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt @@ -888,7 +921,6 @@ subroutine decomp_rate_constants_cn(bounds, & decomp_k(c,j,i_litr1) = k_l1 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_litr2) = k_l2 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_litr3) = k_l3 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt - decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil1) = k_s1 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil2) = k_s2 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt decomp_k(c,j,i_soil3) = k_s3 * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt @@ -897,6 +929,25 @@ subroutine decomp_rate_constants_cn(bounds, & end do end if + ! do the same for cwd, but only if ed is not enabled (because ED handles CWD on its own structure + if (.not. use_ed) then + if (use_vertsoilc) then + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) / dt + end do + end do + else + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + decomp_k(c,j,i_cwd) = k_frag * t_scalar(c,j) * w_scalar(c,j) * o_scalar(c,j) / dt + end do + end do + end if + end if + end associate end subroutine decomp_rate_constants_cn diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 index 3d98c29fc2..3b6b2fe82e 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 @@ -11,7 +11,7 @@ module SoilBiogeochemDecompMod use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use clm_varpar , only : nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools - use clm_varctl , only : use_nitrif_denitrif, use_lch4 + use clm_varctl , only : use_nitrif_denitrif, use_lch4, use_ed use clm_varcon , only : dzsoi_decomp use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use SoilBiogeochemStateType , only : soilbiogeochem_state_type @@ -132,6 +132,7 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, ! column loop to calculate actual immobilization and decomp rates, following ! resolution of plant/heterotroph competition for mineral N + if ( .not. use_ed) then ! calculate c:n ratios of applicable pools do l = 1, ndecomp_pools if ( floating_cn_ratio_decomp_pools(l) ) then @@ -203,6 +204,20 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, end do end do end do + else + do k = 1, ndecomp_cascade_transitions + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + ! + decomp_cascade_hr_vr(c,j,k) = rf_decomp_cascade(c,j,k) * p_decomp_cpool_loss(c,j,k) + ! + decomp_cascade_ctransfer_vr(c,j,k) = (1._r8 - rf_decomp_cascade(c,j,k)) * p_decomp_cpool_loss(c,j,k) + ! + end do + end do + end do + end if if (use_lch4) then ! Calculate total fraction of potential HR, for methane code @@ -221,13 +236,14 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, end do end do - ! Nitrogen limitation / (low)-moisture limitation + + ! Nitrogen limitation / (low)-moisture limitation do j = 1,nlevdecomp do fc = 1,num_soilc c = filter_soilc(fc) if (phr_vr(c,j) > 0._r8) then fphr(c,j) = hrsum(c,j) / phr_vr(c,j) * w_scalar(c,j) - fphr(c,j) = max(fphr(c,j), 0.01_r8) ! Prevent overflow errors for 0 respiration + fphr(c,j) = max(fphr(c,j), 0.01_r8) ! Prevent overflow errors for 0 respiration else fphr(c,j) = 1._r8 end if @@ -235,12 +251,19 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, end do end if - ! vertically integrate net and gross mineralization fluxes for diagnostic output - do j = 1,nlevdecomp - do fc = 1,num_soilc - c = filter_soilc(fc) - net_nmin(c) = net_nmin(c) + net_nmin_vr(c,j) * dzsoi_decomp(j) - gross_nmin(c) = gross_nmin(c) + gross_nmin_vr(c,j) * dzsoi_decomp(j) + + ! vertically integrate net and gross mineralization fluxes for diagnostic output + + do fc = 1,num_soilc + c = filter_soilc(fc) + do j = 1,nlevdecomp + if(.not.use_ed)then + net_nmin(c) = net_nmin(c) + net_nmin_vr(c,j) * dzsoi_decomp(j) + gross_nmin(c) = gross_nmin(c) + gross_nmin_vr(c,j) * dzsoi_decomp(j) + ! else + ! net_nmin(c) = 0.0_r8 + ! gross_nmin(c) = 0.0_r8 + endif end do end do diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemLittVertTranspMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemLittVertTranspMod.F90 index 24c993cb22..430c9d00a0 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemLittVertTranspMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemLittVertTranspMod.F90 @@ -5,7 +5,7 @@ module SoilBiogeochemLittVertTranspMod ! use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg - use clm_varctl , only : iulog, use_c13, use_c14, spinup_state, use_vertsoilc + use clm_varctl , only : iulog, use_c13, use_c14, spinup_state, use_vertsoilc, use_ed, use_cn use clm_varcon , only : secspday use decompMod , only : bounds_type use abortutils , only : endrun @@ -183,6 +183,9 @@ subroutine SoilBiogeochemLittVertTransp(bounds, num_soilc, filter_soilc, & if ( use_c14 ) then ntype = ntype+1 endif + if ( use_ed ) then + ntype = 1 + endif spinup_term = 1._r8 epsilon = 1.e-30 @@ -247,9 +250,11 @@ subroutine SoilBiogeochemLittVertTransp(bounds, num_soilc, filter_soilc, & source => soilbiogeochem_carbonflux_inst%decomp_cpools_sourcesink_col trcr_tendency_ptr => soilbiogeochem_carbonflux_inst%decomp_cpools_transport_tendency_col case (2) ! N - conc_ptr => soilbiogeochem_nitrogenstate_inst%decomp_npools_vr_col - source => soilbiogeochem_nitrogenflux_inst%decomp_npools_sourcesink_col - trcr_tendency_ptr => soilbiogeochem_nitrogenflux_inst%decomp_npools_transport_tendency_col + if (use_cn ) then + conc_ptr => soilbiogeochem_nitrogenstate_inst%decomp_npools_vr_col + source => soilbiogeochem_nitrogenflux_inst%decomp_npools_sourcesink_col + trcr_tendency_ptr => soilbiogeochem_nitrogenflux_inst%decomp_npools_transport_tendency_col + endif case (3) if ( use_c13 ) then ! C13 diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 index e5be13f4fd..44794ca50d 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 @@ -16,6 +16,7 @@ module SoilBiogeochemPotentialMod use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type use SoilBiogeochemNitrogenStateType , only : soilbiogeochem_nitrogenstate_type use SoilBiogeochemNitrogenFluxType , only : soilbiogeochem_nitrogenflux_type + use clm_varctl , only : use_ed, iulog ! implicit none private @@ -130,6 +131,7 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & fphr => soilbiogeochem_carbonflux_inst%fphr_col & ! Output: [real(r8) (:,:) ] fraction of potential SOM + LITTER heterotrophic ) + if ( .not. use_ed ) then ! set initial values for potential C and N fluxes p_decomp_cpool_loss(begc:endc, :, :) = 0._r8 pmnf_decomp_cascade(begc:endc, :, :) = 0._r8 @@ -223,6 +225,20 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & potential_immob_vr(c,j) = immob(c,j) end do end do + else ! use_ed + ! As a first step we are making this a C-only model, so no N downregulation of fluxes. + do k = 1, ndecomp_cascade_transitions + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + ! + p_decomp_cpool_loss(c,j,k) = decomp_cpools_vr(c,j,cascade_donor_pool(k)) & + * decomp_k(c,j,cascade_donor_pool(k)) * pathfrac_decomp_cascade(c,j,k) + ! + end do + end do + end do + end if ! Add up potential hr for methane calculations do j = 1,nlevdecomp diff --git a/components/clm/src/soilbiogeochem/SoilBiogeochemPrecisionControlMod.F90 b/components/clm/src/soilbiogeochem/SoilBiogeochemPrecisionControlMod.F90 index e99e745631..da3789c183 100644 --- a/components/clm/src/soilbiogeochem/SoilBiogeochemPrecisionControlMod.F90 +++ b/components/clm/src/soilbiogeochem/SoilBiogeochemPrecisionControlMod.F90 @@ -31,7 +31,7 @@ subroutine SoilBiogeochemPrecisionControl(num_soilc, filter_soilc, & ! they get too small. ! ! !USES: - use clm_varctl , only : iulog, use_c13, use_c14, use_nitrif_denitrif + use clm_varctl , only : iulog, use_c13, use_c14, use_nitrif_denitrif, use_cn use clm_varpar , only : nlevdecomp use CNSharedParamsMod, only: use_fun ! @@ -98,8 +98,10 @@ subroutine SoilBiogeochemPrecisionControl(num_soilc, filter_soilc, & cc = cc + cs%decomp_cpools_vr_col(c,j,k) cs%decomp_cpools_vr_col(c,j,k) = 0._r8 - cn = cn + ns%decomp_npools_vr_col(c,j,k) - ns%decomp_npools_vr_col(c,j,k) = 0._r8 + if (use_cn) then + cn = cn + ns%decomp_npools_vr_col(c,j,k) + ns%decomp_npools_vr_col(c,j,k) = 0._r8 + endif if ( use_c13 ) then cc13 = cc13 + c13cs%decomp_cpools_vr_col(c,j,k) @@ -117,7 +119,10 @@ subroutine SoilBiogeochemPrecisionControl(num_soilc, filter_soilc, & ! be getting the N truncation flux anyway. cs%ctrunc_vr_col(c,j) = cs%ctrunc_vr_col(c,j) + cc - ns%ntrunc_vr_col(c,j) = ns%ntrunc_vr_col(c,j) + cn + + if (use_cn) then + ns%ntrunc_vr_col(c,j) = ns%ntrunc_vr_col(c,j) + cn + endif if ( use_c13 ) then c13cs%ctrunc_vr_col(c,j) = c13cs%ctrunc_vr_col(c,j) + cc13 endif