diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index a172ac445b..30809918ee 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -90,7 +90,7 @@ OPTIONS sp = Satellite Phenology (SP) This toggles off the namelist variable: use_cn bgc = Carbon Nitrogen with methane, nitrification, vertical soil C, - CENTURY decomposition + CENTURY or MIMICS decomposition This toggles on the namelist variables: use_cn, use_lch4, use_nitrif_denitrif fates = FATES/Ecosystem Demography with below ground BGC @@ -852,7 +852,6 @@ sub setup_cmdl_bgc { $log->fatal_error("The namelist variable use_fates is inconsistent with the -bgc option"); } - # Now set use_cn and use_fates foreach $var ( "use_cn", "use_fates" ) { $val = $nl_flags->{$var}; @@ -876,7 +875,7 @@ sub setup_cmdl_bgc { $log->fatal_error("$var must NOT be None if use_cn or use_fates are on"); } } elsif ( $soil_decomp_method ne "None" ) { - $log->fatal_error("$var must be None if use_cn or use_fates are not"); + $log->fatal_error("$var must be None if use_cn and use_fates are off"); } # # Soil decomposition control variables, methane and Nitrification-Denitrification @@ -2825,7 +2824,6 @@ sub setup_logic_bgc_shared { if ( $nl_flags->{'bgc_mode'} ne "sp" ) { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'constrain_stress_deciduous_onset', 'phys'=>$physv->as_string() ); } - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'decomp_depth_efolding', 'phys'=>$physv->as_string() ); } #------------------------------------------------------------------------------- diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 1ee46abecb..a9535ffd33 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -481,9 +481,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case). -lnd/clm2/paramdata/ctsm51_params.c210923.nc -lnd/clm2/paramdata/clm50_params.c210803.nc -lnd/clm2/paramdata/clm45_params.c210803.nc +lnd/clm2/paramdata/ctsm51_params.c211112.nc +lnd/clm2/paramdata/clm50_params.c211112.nc +lnd/clm2/paramdata/clm45_params.c211112.nc @@ -566,11 +566,6 @@ attributes from the config_cache.xml file (with keys converted to upper-case). .false. Constant - -0.5 -10.0 -10.0 - .false. .true. @@ -2535,7 +2530,7 @@ lnd/clm2/surfdata_map/release-clm5.0.30/surfdata_ne0np4.CONUS.ne30x8_hist_78pfts >lnd/clm2/paramdata/finundated_inversiondata_0.9x1_ESMFmesh_cdf5_130621.nc - + CENTURYKoven2013 diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index a79991efd8..84564c6aa3 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -2127,12 +2127,6 @@ Base advective flux (downwards) for SOM. Maximum depth to mix soils to by croturbation, in permafrost soils. - - -E-folding depth over which decomposition is slowed with depth in all soils. - - If TRUE, reduce heterotrophic respiration according to available oxygen predicted by CH4 submodel. @@ -2169,10 +2163,11 @@ Number of days over which to use exponential relaxation of NPP in N fixation cal + group="soilbgc_decomp" valid_values="None,CENTURYKoven2013,MIMICSWieder2015" > Soil decomposition method None -- No soil decomposition is done CENTURYKoven2013 -- CENTURY model in CTSM from Koven et. al. 2013 + MIMICSWieder2015 -- MIMICS model in CTSM from Wieder et. al. 2015 An active soil decomposition method requires the BGC or FATES model to work And both BGC and FATES models require an active soil decomposition model diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index 77fc2053a8..34b4c7363d 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -522,6 +522,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/README b/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/README index e8eee75e53..39ac318b86 100644 --- a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/README +++ b/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/README @@ -1,5 +1,5 @@ This points to an alternate params file with -rf_cwdl2_bgc = 0.5 +rf_cwdl2 = 0.5 rf_cwdl3_bgc = 0.5 while by default these parameters equal zero. diff --git a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm index b7672df8ca..7ae4a69aad 100644 --- a/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm +++ b/cime_config/testdefs/testmods_dirs/clm/ciso_cwd_hr/user_nl_clm @@ -1,2 +1,2 @@ -paramfile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_params.c210827.nc' +paramfile = '/glade/p/cesm/cseg/inputdata/lnd/clm2/paramdata/ctsm51_ciso_cwd_hr_params.c211112.nc' hist_fincl1 = 'CWDC_HR','C13_CWDC_HR','C14_CWDC_HR','CWD_HR_L2','CWD_HR_L2_vr','CWD_HR_L3','CWD_HR_L3_vr' diff --git a/cime_config/testdefs/testmods_dirs/clm/mimics/README b/cime_config/testdefs/testmods_dirs/clm/mimics/README new file mode 100644 index 0000000000..082e324f4d --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/mimics/README @@ -0,0 +1,6 @@ +This test mod turns on the MIMICS instead of the CENTURY below-ground +biogeochemistry. + +As of 2021/11/15 this test inherits mods from the ../coldStart directory; +however the plan is to change that to ../default when we have an acceptable +spun up finidat file to initialize MIMICS with. diff --git a/cime_config/testdefs/testmods_dirs/clm/mimics/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/mimics/include_user_mods new file mode 100644 index 0000000000..3dbc936366 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/mimics/include_user_mods @@ -0,0 +1 @@ +../coldStart diff --git a/cime_config/testdefs/testmods_dirs/clm/mimics/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/mimics/user_nl_clm new file mode 100644 index 0000000000..152d91b21e --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/mimics/user_nl_clm @@ -0,0 +1 @@ +soil_decomp_method = 'MIMICSWieder2015' diff --git a/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/README b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/README new file mode 100644 index 0000000000..1363ea696b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/README @@ -0,0 +1,8 @@ +slevis 2022/2/1: +This testmod tests the newton_krylov_spinup usermod found in +cime_config/usermods_dirs/newton_krylov_spinup +combined with the mimics testmod. + +By default this usermod runs for 20 yrs and writes a history file at that +time. To save time running the test-suite, we're changing +hist_nhtfrq = 0,-43800 (from 0,-175200) diff --git a/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/include_user_mods new file mode 100644 index 0000000000..192aa023a8 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/include_user_mods @@ -0,0 +1,2 @@ +../mimics +../../../../usermods_dirs/newton_krylov_spinup diff --git a/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/user_nl_clm new file mode 100644 index 0000000000..bce7009c36 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/newton_krylov_spinup/user_nl_clm @@ -0,0 +1 @@ +hist_nhtfrq = 0,-43800 diff --git a/cime_config/usermods_dirs/newton_krylov_spinup/README b/cime_config/usermods_dirs/newton_krylov_spinup/README new file mode 100644 index 0000000000..037709ee7d --- /dev/null +++ b/cime_config/usermods_dirs/newton_krylov_spinup/README @@ -0,0 +1,4 @@ +user_nl_clm contains this line: +hist_nhtfrq = 0,-175200 +0 means monthly history frequency and +-175200 means 20 years history frequency (indicated in hours) diff --git a/cime_config/usermods_dirs/newton_krylov_spinup/shell_commands b/cime_config/usermods_dirs/newton_krylov_spinup/shell_commands new file mode 100644 index 0000000000..a66f52f6fd --- /dev/null +++ b/cime_config/usermods_dirs/newton_krylov_spinup/shell_commands @@ -0,0 +1,3 @@ +#!/bin/bash +./xmlchange CLM_FORCE_COLDSTART="on" + diff --git a/cime_config/usermods_dirs/newton_krylov_spinup/user_nl_clm b/cime_config/usermods_dirs/newton_krylov_spinup/user_nl_clm new file mode 100644 index 0000000000..318105a043 --- /dev/null +++ b/cime_config/usermods_dirs/newton_krylov_spinup/user_nl_clm @@ -0,0 +1,15 @@ +hist_dov2xy = .true.,.false. +hist_nhtfrq = 0,-175200 +hist_mfilt = 1,1 +hist_fincl2 = 'FPI_vr', 'K_PAS_SOM', 'K_SLO_SOM', 'K_ACT_SOM', + 'K_CWD', 'K_CEL_LIT', 'K_LIG_LIT', 'K_MET_LIT', + 'CWD_PATHFRAC_L2_vr', 'CWD_RESP_FRAC_L2_vr', + 'CWD_PATHFRAC_L3_vr', 'CWD_RESP_FRAC_L3_vr', + 'L1_PATHFRAC_S1_vr', 'L1_RESP_FRAC_S1_vr', + 'L2_PATHFRAC_S1_vr', 'L2_RESP_FRAC_S1_vr', + 'L3_PATHFRAC_S2_vr', 'L3_RESP_FRAC_S2_vr', + 'S1_PATHFRAC_S2_vr', 'S1_RESP_FRAC_S2_vr' + 'S1_PATHFRAC_S3_vr', 'S1_RESP_FRAC_S3_vr' + 'S2_PATHFRAC_S1_vr', 'S2_RESP_FRAC_S1_vr' + 'S2_PATHFRAC_S3_vr', 'S2_RESP_FRAC_S3_vr' + 'S3_PATHFRAC_S1_vr', 'S3_RESP_FRAC_S1_vr' diff --git a/doc/ChangeLog b/doc/ChangeLog index 027ec028e9..a93c31092f 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,118 @@ =============================================================== +Tag name: ctsm5.1.dev074 +Originator(s): slevis (Samuel Levis,SLevis Consulting,303-665-1310) +Date: Wed Feb 2 00:44:27 MST 2022 +One-line Summary: Introduce vert. resolved MIMICS as new method to solve below ground decomp. + +Purpose and description of changes +---------------------------------- + + Introducing new option to solve below ground decomposition: MIMICSWieder2015. + The old option (CENTURYKoven2013) remains available. User will select one or + the other via the namelist variable soil_decomp_method. + + Elin's MIMICS+ github issue relates: #1260. + MIMICS spinup issues relate and are listed in caveats for users. + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + + +Bugs fixed or introduced +------------------------ +Issues fixed (include CTSM Issue #): + #1473 : decomp_depth_efolding appears (with the same value) in the namelist + and the params files; the model ends up using the one from the namelist + + +Notes of particular relevance for users +--------------------------------------- +Caveats for users (e.g., need to interpolate initial conditions): + MIMICS spinup method is in development + Issues related to spinup + #1451 -- MIMICS spinup + #1455 -- Additional history output for Newton-Krylov spinup + #1457 -- Adding history variables for Newton-Krylov spinup + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): + New option MIMICSWieder2015 for use with namelist variable soil_decomp_method + decomp_depth_efolding removed as values on param files was already being used + +Changes to the datasets (e.g., parameter, surface or initial files): + CLM parameter files now include new MIMICS-related parameters prefixed with + mimics_. Existing CENTURYKoven-specific parameters were renamed to include the + prefix bgc_. Parameters used by both CENTURYKoven2013 and MIMICSWieder2015 have + neither prefix. + +Substantial timing or memory changes: + Timing comparison between CENTURYKoven (BGC for short) and MIMICSWieder: + In a 1x1_brazil test writing annual output and restarts every 10 yrs, + BGC took ~2.8 hrs/100 yrs, MIMICS took ~3.6 hrs/100 yrs. + + More timing comparisons shown here: + https://github.com/ESCOMP/CTSM/pull/1318#issuecomment-1008211485 + +Notes of particular relevance for developers: +--------------------------------------------- +Changes to tests or testing: + New test-suite tests: + ERP_D_P36x2_Ld3.f10_f10_mg37.I1850Clm51BgcCrop.cheyenne_gnu.clm-mimics + SMS_Ld5_Mmpi-serial.1x1_brazil.IHistClm50BgcQianRs.izumi_gnu.clm-mimics + SMS_Ly5_Mmpi-serial.1x1_brazil.IHistClm50BgcQianRs.cheyenne_intel.clm-newton_krylov_spinup + +Testing summary: +---------------- + + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + cheyenne - PASS (831 tests are different from baseline because of new param files) + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + cheyenne ---- OK (black test fails) + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + cheyenne ---- PASS + izumi ------- OK (see one answer change below) + + +Answer changes +-------------- + +Changes answers relative to baseline: NO (almost) + + One test in the izumi test-suite fails + SMS_Vmct.f10_f10_mg37.I2000Clm50BgcCrop.izumi_pgi.clm-crop + and this is due to diffs from baseline (dev073). + cprnc.out shows diffs only in one variable: + RMS TOTECOSYSC 1.9853E-12 NORMALIZED 6.5087E-17 + + +Other details +------------- + +Pull Requests that document the changes (include PR ids): + https://github.com/ESCOMP/ctsm/pull/1318 -- MIMICS into vertically resolved CTSM + +=============================================================== +=============================================================== Tag name: ctsm5.1.dev073 Originator(s): sacks (Bill Sacks) Date: Tue Jan 25 16:33:06 MST 2022 diff --git a/doc/ChangeSum b/doc/ChangeSum index bcf6687654..93e87f9495 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.1.dev074 slevis 02/02/2022 Introduce vert. resolved MIMICS as new method to solve below ground decomp. ctsm5.1.dev073 sacks 01/25/2022 Some fixes for Gregorian calendar ctsm5.1.dev072 negins 01/17/2022 mksurfdat toolchain part 1: gen_mksurf_namelist ctsm5.1.dev071 glemieux 01/16/2022 Small changes to enable new fates dimension and update fates tag diff --git a/src/biogeochem/CNBalanceCheckMod.F90 b/src/biogeochem/CNBalanceCheckMod.F90 index 67e83713c3..925c3479d0 100644 --- a/src/biogeochem/CNBalanceCheckMod.F90 +++ b/src/biogeochem/CNBalanceCheckMod.F90 @@ -325,7 +325,7 @@ subroutine CBalanceCheck(this, bounds, num_soilc, filter_soilc, & do g = bounds%begg, bounds%endg ! calculate gridcell-level carbon storage for mass conservation check ! Notes: - ! totgrcc = totcolc = totc_p2c_col(c) + soilbiogeochem_cwdc_col(c) + soilbiogeochem_totlitc_col(c) + soilbiogeochem_totsomc_col(c) + soilbiogeochem_ctrunc_col(c) + ! totgrcc = totcolc = totc_p2c_col(c) + soilbiogeochem_cwdc_col(c) + soilbiogeochem_totlitc_col(c) + soilbiogeochem_totmicc_col(c) + soilbiogeochem_totsomc_col(c) + soilbiogeochem_ctrunc_col(c) ! totc_p2c_col = totc_patch = totvegc_patch(p) + xsmrpool_patch(p) + ctrunc_patch(p) + cropseedc_deficit_patch(p) ! Not including seedc_grc in grc_begcb and grc_endcb because ! seedc_grc forms out of thin air, for now, and equals @@ -540,7 +540,7 @@ subroutine NBalanceCheck(this, bounds, num_soilc, filter_soilc, & write(iulog,*)'output mass = ',col_noutputs(c)*dt write(iulog,*)'net flux = ',(col_ninputs(c)-col_noutputs(c))*dt write(iulog,*)'inputs,ffix,nfix,ndep = ',ffix_to_sminn(c)*dt,nfix_to_sminn(c)*dt,ndep_to_sminn(c)*dt - write(iulog,*)'outputs,ffix,nfix,ndep = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt + write(iulog,*)'outputs,lch,roff,dnit = ',smin_no3_leached(c)*dt, smin_no3_runoff(c)*dt,f_n2o_nit(c)*dt call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, msg=errMsg(sourcefile, __LINE__)) end if diff --git a/src/biogeochem/CNDriverMod.F90 b/src/biogeochem/CNDriverMod.F90 index 2f51846606..45119a452a 100644 --- a/src/biogeochem/CNDriverMod.F90 +++ b/src/biogeochem/CNDriverMod.F90 @@ -11,8 +11,8 @@ module CNDriverMod use decompMod , only : bounds_type use perf_mod , only : t_startf, t_stopf use clm_varctl , only : use_nitrif_denitrif, use_nguardrail - use SoilBiogeochemDecompCascadeConType, only : century_decomp, decomp_method - use clm_varctl , only : use_crop + use clm_varctl , only : iulog, use_crop + use SoilBiogeochemDecompCascadeConType, only : mimics_decomp, century_decomp, decomp_method use CNSharedParamsMod , only : use_fun use CNVegStateType , only : cnveg_state_type use CNVegCarbonStateType , only : cnveg_carbonstate_type @@ -129,6 +129,7 @@ subroutine CNDriverNoLeaching(bounds, use CNGapMortalityMod , only: CNGapMortality use CNSharedParamsMod , only: use_fun use dynHarvestMod , only: CNHarvest + use SoilBiogeochemDecompCascadeMIMICSMod, only: decomp_rates_mimics use SoilBiogeochemDecompCascadeBGCMod , only: decomp_rate_constants_bgc use SoilBiogeochemCompetitionMod , only: SoilBiogeochemCompetition use SoilBiogeochemDecompMod , only: SoilBiogeochemDecomp @@ -200,6 +201,8 @@ subroutine CNDriverNoLeaching(bounds, 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):: p_decomp_npool_to_din(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) ! potential flux to dissolved inorganic N + real(r8):: p_decomp_cn_gain(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools) ! C:N ratio of the flux gained by the receiver pool real(r8):: arepr(bounds%begp:bounds%endp) ! reproduction allocation coefficient (only used for use_crop) real(r8):: aroot(bounds%begp:bounds%endp) ! root allocation coefficient (only used for use_crop) integer :: begp,endp @@ -311,18 +314,28 @@ subroutine CNDriverNoLeaching(bounds, !-------------------------------------------- call t_startf('SoilBiogeochem') + call t_startf('DecompRate') if (decomp_method == century_decomp) then call decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) + else if (decomp_method == mimics_decomp) then + call decomp_rates_mimics(bounds, num_soilc, filter_soilc, & + soilstate_inst, temperature_inst, cnveg_carbonflux_inst, ch4_inst, & + soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst) end if + call t_stopf('DecompRate') + call t_startf('SoilBiogeochemPotential') ! 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)) + p_decomp_cn_gain=p_decomp_cn_gain(begc:endc,1:nlevdecomp,1:ndecomp_pools), & + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + p_decomp_npool_to_din=p_decomp_npool_to_din(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) + call t_stopf('SoilBiogeochemPotential') ! calculate vertical profiles for distributing soil and litter C and N (previously subroutine decomp_vertprofiles called from CNDecompAlloc) call SoilBiogeochemVerticalProfile(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & @@ -395,7 +408,8 @@ subroutine CNDriverNoLeaching(bounds, call t_startf('soilbiogeochemcompetition') - call SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, filter_soilp, waterstatebulk_inst, & + call SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, filter_soilp, & + p_decomp_cn_gain, pmnf_decomp_cascade, waterstatebulk_inst, & waterfluxbulk_inst,temperature_inst,soilstate_inst,cnveg_state_inst, & cnveg_carbonstate_inst ,& cnveg_carbonflux_inst,cnveg_nitrogenstate_inst,cnveg_nitrogenflux_inst, & @@ -436,7 +450,8 @@ subroutine CNDriverNoLeaching(bounds, 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)) + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + p_decomp_npool_to_din=p_decomp_npool_to_din(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) call t_stopf('SoilBiogeochemDecomp') @@ -937,6 +952,7 @@ subroutine CNDriverSummarizeStates(bounds, num_allc, filter_allc, & num_soilc, filter_soilc, num_soilp, filter_soilp, & soilbiogeochem_cwdc_col=soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & soilbiogeochem_totlitc_col=soilbiogeochem_carbonstate_inst%totlitc_col(begc:endc), & + soilbiogeochem_totmicc_col=soilbiogeochem_carbonstate_inst%totmicc_col(begc:endc), & soilbiogeochem_totsomc_col=soilbiogeochem_carbonstate_inst%totsomc_col(begc:endc), & soilbiogeochem_ctrunc_col=soilbiogeochem_carbonstate_inst%ctrunc_col(begc:endc)) @@ -945,6 +961,7 @@ subroutine CNDriverSummarizeStates(bounds, num_allc, filter_allc, & num_soilc, filter_soilc, num_soilp, filter_soilp, & soilbiogeochem_cwdc_col=c13_soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & soilbiogeochem_totlitc_col=c13_soilbiogeochem_carbonstate_inst%totlitc_col(begc:endc), & + soilbiogeochem_totmicc_col=c13_soilbiogeochem_carbonstate_inst%totmicc_col(begc:endc), & soilbiogeochem_totsomc_col=c13_soilbiogeochem_carbonstate_inst%totsomc_col(begc:endc), & soilbiogeochem_ctrunc_col=c13_soilbiogeochem_carbonstate_inst%ctrunc_col(begc:endc)) end if @@ -954,6 +971,7 @@ subroutine CNDriverSummarizeStates(bounds, num_allc, filter_allc, & num_soilc, filter_soilc, num_soilp, filter_soilp, & soilbiogeochem_cwdc_col=c14_soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & soilbiogeochem_totlitc_col=c14_soilbiogeochem_carbonstate_inst%totlitc_col(begc:endc), & + soilbiogeochem_totmicc_col=c14_soilbiogeochem_carbonstate_inst%totmicc_col(begc:endc), & soilbiogeochem_totsomc_col=c14_soilbiogeochem_carbonstate_inst%totsomc_col(begc:endc), & soilbiogeochem_ctrunc_col=c14_soilbiogeochem_carbonstate_inst%ctrunc_col(begc:endc)) end if @@ -975,6 +993,10 @@ subroutine CNDriverSummarizeFluxes(bounds, & soilbiogeochem_carbonflux_inst, & c13_soilbiogeochem_carbonflux_inst, & c14_soilbiogeochem_carbonflux_inst, & + soilbiogeochem_carbonstate_inst, & + c13_soilbiogeochem_carbonstate_inst, & + c14_soilbiogeochem_carbonstate_inst, & + soilbiogeochem_nitrogenstate_inst, & soilbiogeochem_nitrogenflux_inst) ! ! !DESCRIPTION: @@ -999,6 +1021,10 @@ subroutine CNDriverSummarizeFluxes(bounds, & type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst type(soilbiogeochem_carbonflux_type) , intent(inout) :: c13_soilbiogeochem_carbonflux_inst type(soilbiogeochem_carbonflux_type) , intent(inout) :: c14_soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type) , intent(in) :: soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonstate_type) , intent(in) :: c13_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_carbonstate_type) , intent(in) :: c14_soilbiogeochem_carbonstate_inst + type(soilbiogeochem_nitrogenstate_type) , intent(in) :: soilbiogeochem_nitrogenstate_inst type(soilbiogeochem_nitrogenflux_type) , intent(inout) :: soilbiogeochem_nitrogenflux_inst ! ! !LOCAL VARIABLES: @@ -1030,6 +1056,7 @@ subroutine CNDriverSummarizeFluxes(bounds, & ! cnveg carbon/nitrogen flux summary ! ---------------------------------------------- + call t_startf('CNvegCflux_summary') call cnveg_carbonflux_inst%Summary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & isotope='bulk', & soilbiogeochem_hr_col=soilbiogeochem_carbonflux_inst%hr_col(begc:endc), & @@ -1037,6 +1064,8 @@ subroutine CNDriverSummarizeFluxes(bounds, & soilbiogeochem_lithr_col=soilbiogeochem_carbonflux_inst%lithr_col(begc:endc), & soilbiogeochem_decomp_cascade_ctransfer_col=& soilbiogeochem_carbonflux_inst%decomp_cascade_ctransfer_col(begc:endc,1:ndecomp_cascade_transitions), & + soilbiogeochem_cwdc_col=soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & + soilbiogeochem_cwdn_col=soilbiogeochem_nitrogenstate_inst%cwdn_col(begc:endc), & product_closs_grc=c_products_inst%product_loss_grc(begg:endg)) if ( use_c13 ) then @@ -1047,6 +1076,8 @@ subroutine CNDriverSummarizeFluxes(bounds, & soilbiogeochem_lithr_col=c13_soilbiogeochem_carbonflux_inst%lithr_col(begc:endc), & soilbiogeochem_decomp_cascade_ctransfer_col=& c13_soilbiogeochem_carbonflux_inst%decomp_cascade_ctransfer_col(begc:endc,1:ndecomp_cascade_transitions), & + soilbiogeochem_cwdc_col=c13_soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & + soilbiogeochem_cwdn_col=soilbiogeochem_nitrogenstate_inst%cwdn_col(begc:endc), & product_closs_grc=c13_products_inst%product_loss_grc(begg:endg)) end if @@ -1058,8 +1089,11 @@ subroutine CNDriverSummarizeFluxes(bounds, & soilbiogeochem_lithr_col=c14_soilbiogeochem_carbonflux_inst%lithr_col(begc:endc), & soilbiogeochem_decomp_cascade_ctransfer_col=& c14_soilbiogeochem_carbonflux_inst%decomp_cascade_ctransfer_col(begc:endc,1:ndecomp_cascade_transitions), & + soilbiogeochem_cwdc_col=c14_soilbiogeochem_carbonstate_inst%cwdc_col(begc:endc), & + soilbiogeochem_cwdn_col=soilbiogeochem_nitrogenstate_inst%cwdn_col(begc:endc), & product_closs_grc=c14_products_inst%product_loss_grc(begg:endg)) end if + call t_stopf('CNvegCflux_summary') call cnveg_nitrogenflux_inst%Summary(bounds, num_soilc, filter_soilc, num_soilp, filter_soilp) diff --git a/src/biogeochem/CNFUNMod.F90 b/src/biogeochem/CNFUNMod.F90 index 57465f8fba..37593a2b57 100644 --- a/src/biogeochem/CNFUNMod.F90 +++ b/src/biogeochem/CNFUNMod.F90 @@ -505,7 +505,6 @@ subroutine CNFUN(bounds,num_soilc, filter_soilc,num_soilp& !--------------------------------- associate(ivt => patch%itype , & ! Input: [integer (:) ] p leafcn => pftcon%leafcn , & ! Input: leaf C:N (gC/gN) - lflitcn => pftcon%lflitcn , & ! Input: leaf litter C:N (gC/gN) season_decid => pftcon%season_decid , & ! Input: binary flag for seasonal ! -deciduous leaf habit (0 or 1) stress_decid => pftcon%stress_decid , & ! Input: binary flag for stress diff --git a/src/biogeochem/CNGRespMod.F90 b/src/biogeochem/CNGRespMod.F90 index d95761e61f..44a2225548 100644 --- a/src/biogeochem/CNGRespMod.F90 +++ b/src/biogeochem/CNGRespMod.F90 @@ -60,7 +60,6 @@ subroutine CNGResp(num_soilp, filter_soilp, cnveg_carbonflux_inst, canopystate_i grperc => pftcon%grperc , & ! Input: growth respiration parameter grpnow => pftcon%grpnow , & ! Input: growth respiration parameter leafcn => pftcon%leafcn , & ! Input: leaf C:N (gC/gN) - frootcn => pftcon%frootcn , & ! Input: fine root C:N (gC/gN) livewdcn => pftcon%livewdcn , & ! Input: live wood (phloem and ray parenchyma) C:N (gC/gN) laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:)] sunlit projected leaf area index diff --git a/src/biogeochem/CNGapMortalityMod.F90 b/src/biogeochem/CNGapMortalityMod.F90 index a41542f3a1..bd6867b999 100644 --- a/src/biogeochem/CNGapMortalityMod.F90 +++ b/src/biogeochem/CNGapMortalityMod.F90 @@ -134,7 +134,6 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so heatstress => dgvs_inst%heatstress_patch , & ! Input: [real(r8) (:) ] leafcn => pftcon%leafcn , & ! Input: [real(r8) (:)] leaf C:N (gC/gN) - frootcn => pftcon%frootcn , & ! Input: [real(r8) (:)] fine root C:N (gC/gN) livewdcn => pftcon%livewdcn , & ! Input: [real(r8) (:)] live wood (phloem and ray parenchyma) C:N (gC/gN) 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 diff --git a/src/biogeochem/CNNDynamicsMod.F90 b/src/biogeochem/CNNDynamicsMod.F90 index fc99e25a69..7a80d701a8 100644 --- a/src/biogeochem/CNNDynamicsMod.F90 +++ b/src/biogeochem/CNNDynamicsMod.F90 @@ -125,7 +125,7 @@ subroutine CNNDeposition( bounds, & ! directly into the canopy and mineral N entering the soil pool. ! ! !USES: - use CNSharedParamsMod , only: use_fun + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds type(atm2lnd_type) , intent(in) :: atm2lnd_inst @@ -218,7 +218,7 @@ subroutine CNNFixation(num_soilc, filter_soilc, & !----------------------------------------------------------------------- associate( & - cannsum_npp => cnveg_carbonflux_inst%annsum_npp_col , & ! Input: [real(r8) (:)] nitrogen deposition rate (gN/m2/s) + cannsum_npp => cnveg_carbonflux_inst%annsum_npp_col , & ! Input: [real(r8) (:)] (gC/m2/yr) annual sum of NPP, averaged from patch-level col_lag_npp => cnveg_carbonflux_inst%lag_npp_col , & ! Input: [real(r8) (:)] (gC/m2/s) lagged net primary production nfix_to_sminn => soilbiogeochem_nitrogenflux_inst%nfix_to_sminn_col & ! Output: [real(r8) (:)] symbiotic/asymbiotic N fixation to soil mineral N (gN/m2/s) diff --git a/src/biogeochem/CNSharedParamsMod.F90 b/src/biogeochem/CNSharedParamsMod.F90 index f38a7debb5..93e464eda1 100644 --- a/src/biogeochem/CNSharedParamsMod.F90 +++ b/src/biogeochem/CNSharedParamsMod.F90 @@ -15,7 +15,9 @@ module CNSharedParamsMod type, public :: CNParamsShareType real(r8) :: Q10 ! temperature dependence real(r8) :: minpsi ! minimum soil water potential for heterotrophic resp - real(r8) :: cwd_fcel ! cellulose fraction of coarse woody debris + real(r8) :: maxpsi ! maximum soil water potential for heterotrophic resp + real(r8) :: rf_cwdl2 ! respiration fraction in CWD to litter2 transition (frac) + real(r8) :: tau_cwd ! corrected fragmentation rate constant CWD, century leaves wood decomposition rates open, within range of 0 - 0.5 yr^-1 (1/0.3) (1/yr) real(r8) :: cwd_flig ! lignin fraction of coarse woody debris real(r8) :: froz_q10 ! separate q10 for frozen soil respiration rates real(r8) :: decomp_depth_efolding ! e-folding depth for reduction in decomposition (m) @@ -96,16 +98,31 @@ subroutine CNParamsReadShared_netcdf(ncid) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) CNParamsShareInst%minpsi=tempr - tString='cwd_fcel' + tString='maxpsi_hr' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - CNParamsShareInst%cwd_fcel=tempr + CNParamsShareInst%maxpsi=tempr + + tString='rf_cwdl2' + call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + CNParamsShareInst%rf_cwdl2=tempr + + tString='tau_cwd' + call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + CNParamsShareInst%tau_cwd=tempr tString='cwd_flig' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) CNParamsShareInst%cwd_flig=tempr + tString='decomp_depth_efolding' + call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + CNParamsShareInst%decomp_depth_efolding=tempr + tString='froz_q10' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) @@ -149,7 +166,6 @@ subroutine CNParamsReadShared_namelist(namelist_file) integer :: ierr ! error code integer :: unitn ! unit for namelist file - real(r8) :: decomp_depth_efolding = 0.0_r8 logical :: constrain_stress_deciduous_onset = .false. character(len=32) :: subroutine_name = 'CNParamsReadNamelist' @@ -162,7 +178,6 @@ subroutine CNParamsReadShared_namelist(namelist_file) ! ---------------------------------------------------------------------- namelist /bgc_shared/ & - decomp_depth_efolding, & constrain_stress_deciduous_onset @@ -189,18 +204,15 @@ subroutine CNParamsReadShared_namelist(namelist_file) end if ! masterproc ! Broadcast the parameters from master - call shr_mpi_bcast ( decomp_depth_efolding, mpicom ) call shr_mpi_bcast ( constrain_stress_deciduous_onset, mpicom ) ! Save the parameter to the instance - CNParamsShareInst%decomp_depth_efolding = decomp_depth_efolding CNParamsShareInst%constrain_stress_deciduous_onset = constrain_stress_deciduous_onset ! Output read parameters to the lnd.log if (masterproc) then write(iulog,*) 'CN/BGC shared namelist parameters:' write(iulog,*)' ' - write(iulog,*)' decomp_depth_efolding = ', decomp_depth_efolding write(iulog,*)' constrain_stress_deciduous_onset = ',constrain_stress_deciduous_onset write(iulog,*) diff --git a/src/biogeochem/CNVegCarbonFluxType.F90 b/src/biogeochem/CNVegCarbonFluxType.F90 index 9eb226d67f..718f7b0330 100644 --- a/src/biogeochem/CNVegCarbonFluxType.F90 +++ b/src/biogeochem/CNVegCarbonFluxType.F90 @@ -11,13 +11,13 @@ module CNVegCarbonFluxType use decompMod , only : bounds_type use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use clm_varpar , only : ndecomp_cascade_transitions, ndecomp_pools - use clm_varpar , only : nlevdecomp_full, nlevdecomp, i_litr_min, i_litr_max + use clm_varpar , only : nlevdecomp_full, nlevdecomp, i_litr_min, i_litr_max, i_cwdl2 use clm_varcon , only : spval, dzsoi_decomp use clm_varctl , only : use_cndv, use_c13, use_nitrif_denitrif, use_crop use clm_varctl , only : use_grainproduct use clm_varctl , only : iulog use landunit_varcon , only : istsoil, istcrop, istdlak - use pftconMod , only : npcropmin + use pftconMod , only : npcropmin, pftcon use LandunitType , only : lun use ColumnType , only : col use PatchType , only : patch @@ -334,6 +334,9 @@ module CNVegCarbonFluxType real(r8), pointer :: nbp_grc (:) ! (gC/m2/s) net biome production, includes fire, landuse, harvest and hrv_xsmrpool flux, positive for sink (same as net carbon exchange between land and atmosphere) real(r8), pointer :: nee_grc (:) ! (gC/m2/s) net ecosystem exchange of carbon, includes fire and hrv_xsmrpool, excludes landuse and harvest flux, positive for source + ! Plant C to N ratios + real(r8), pointer :: ligninNratioAvg_col(:) ! Average of leaf, fine root, and CWD lignin to N ratio + ! Dynamic landcover fluxnes real(r8), pointer :: landuseflux_grc(:) ! (gC/m2/s) dwt_conv_cflux+product_closs real(r8), pointer :: npp_Nactive_patch (:) ! C used by mycorrhizal uptake (gC/m2/s) @@ -697,6 +700,8 @@ subroutine InitAllocate(this, bounds, carbon_type) allocate(this%annsum_npp_col (begc:endc)) ; this%annsum_npp_col (:) = nan allocate(this%lag_npp_col (begc:endc)) ; this%lag_npp_col (:) = spval + allocate(this%ligninNratioAvg_col (begc:endc)) ; this%ligninNratioAvg_col (:) = nan + allocate(this%nep_col (begc:endc)) ; this%nep_col (:) = nan allocate(this%nbp_grc (begg:endg)) ; this%nbp_grc (:) = nan allocate(this%nee_grc (begg:endg)) ; this%nee_grc (:) = nan @@ -3411,12 +3416,12 @@ subroutine InitCold(this, bounds) if (lun%ifspecial(l)) then this%annsum_npp_col(c) = spval + this%ligninNratioAvg_col(c) = spval end if - ! also initialize dynamic landcover fluxes so that they have - ! real values on first timestep, prior to calling pftdyn_cnbal if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then this%annsum_npp_col(c) = 0._r8 + this%ligninNratioAvg_col(c) = 0._r8 end if end do @@ -3584,6 +3589,11 @@ subroutine RestartBulkOnly ( this, bounds, ncid, flag ) long_name='', units='', & interpinic_flag='interp', readvar=readvar, data=this%annsum_npp_col) + call restartvar(ncid=ncid, flag=flag, varname='ligninNratioAvg', xtype=ncd_double, & + dim1name='column', & + long_name='', units='', & + interpinic_flag='interp', readvar=readvar, data=this%ligninNratioAvg_col) + call restartvar(ncid=ncid, flag=flag, varname='tempsum_litfall', xtype=ncd_double, & dim1name='pft', & long_name='', units='', & @@ -3993,6 +4003,7 @@ subroutine Summary_carbonflux(this, & bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, isotope, & soilbiogeochem_hr_col, soilbiogeochem_cwdhr_col, soilbiogeochem_lithr_col, & soilbiogeochem_decomp_cascade_ctransfer_col, & + soilbiogeochem_cwdc_col, soilbiogeochem_cwdn_col, & product_closs_grc) ! ! !DESCRIPTION: @@ -4003,8 +4014,8 @@ subroutine Summary_carbonflux(this, & use clm_varcon , only: secspday use clm_varctl , only: nfix_timeconst, carbon_resp_opt use subgridAveMod , only: p2c, c2g - use SoilBiogeochemDecompCascadeConType , only: decomp_cascade_con - use CNSharedParamsMod , only: use_fun + use SoilBiogeochemDecompCascadeConType , only: decomp_cascade_con, mimics_decomp, decomp_method + use CNSharedParamsMod , only: use_fun, CNParamsShareInst ! ! !ARGUMENTS: class(cnveg_carbonflux_type) :: this @@ -4018,6 +4029,8 @@ subroutine Summary_carbonflux(this, & real(r8) , intent(in) :: soilbiogeochem_cwdhr_col(bounds%begc:) real(r8) , intent(in) :: soilbiogeochem_lithr_col(bounds%begc:) real(r8) , intent(in) :: soilbiogeochem_decomp_cascade_ctransfer_col(bounds%begc:,1:) + real(r8) , intent(in) :: soilbiogeochem_cwdc_col(bounds%begc:) + real(r8) , intent(in) :: soilbiogeochem_cwdn_col(bounds%begc:) real(r8) , intent(in) :: product_closs_grc(bounds%begg:) ! ! !LOCAL VARIABLES: @@ -4025,6 +4038,13 @@ subroutine Summary_carbonflux(this, & integer :: fp,fc ! lake filter indices real(r8) :: nfixlags, dtime ! temp variables for making lagged npp real(r8) :: maxdepth ! depth to integrate soil variables + real(r8) :: ligninNratio_cwd ! lignin to N ratio of CWD + real(r8) :: ligninNratio_leaf_patch(bounds%begp:bounds%endp) ! lignin to N ratio of leaves, patch level + real(r8) :: ligninNratio_froot_patch(bounds%begp:bounds%endp) ! lignin to N ratio of fine roots, patch level + real(r8) :: ligninNratio_leaf_col(bounds%begc:bounds%endc) ! lignin to N ratio of leaves, column level + real(r8) :: ligninNratio_froot_col(bounds%begc:bounds%endc) ! lignin to N ratio of fine roots, column level + real(r8) :: leafc_to_litter_col(bounds%begc:bounds%endc) ! leaf C to litter C, column level + real(r8) :: frootc_to_litter_col(bounds%begc:bounds%endc) ! fine root C to litter C, column level real(r8) :: nep_grc(bounds%begg:bounds%endg) ! nep_col averaged to gridcell real(r8) :: fire_closs_grc(bounds%begg:bounds%endg) ! fire_closs_col averaged to gridcell real(r8) :: hrv_xsmrpool_to_atm_grc(bounds%begg:bounds%endg) ! hrv_xsmrpool_to_atm_col averaged to gridcell (gC/m2/s) @@ -4396,6 +4416,17 @@ subroutine Summary_carbonflux(this, & this%hrv_gresp_storage_to_litter_patch(p) + & this%hrv_gresp_xfer_to_litter_patch(p) + if (decomp_method == mimics_decomp) then + ! Calculate ligninNratio for leaves and fine roots + associate(ivt => patch%itype) ! Input: [integer (:)] patch plant type + ligninNratio_leaf_patch(p) = pftcon%lf_flig(ivt(p)) * & + pftcon%lflitcn(ivt(p)) * & + this%leafc_to_litter_patch(p) + ligninNratio_froot_patch(p) = pftcon%fr_flig(ivt(p)) * & + pftcon%frootcn(ivt(p)) * & + this%frootc_to_litter_patch(p) + end associate + end if end do ! end of patches loop !------------------------------------------------ @@ -4440,6 +4471,39 @@ subroutine Summary_carbonflux(this, & this%gpp_patch(bounds%begp:bounds%endp), & this%gpp_col(bounds%begc:bounds%endc)) + if (decomp_method == mimics_decomp) then + call p2c(bounds, num_soilc, filter_soilc, & + ligninNratio_leaf_patch(bounds%begp:bounds%endp), & + ligninNratio_leaf_col(bounds%begc:bounds%endc)) + call p2c(bounds, num_soilc, filter_soilc, & + ligninNratio_froot_patch(bounds%begp:bounds%endp), & + ligninNratio_froot_col(bounds%begc:bounds%endc)) + call p2c(bounds, num_soilc, filter_soilc, & + this%leafc_to_litter_patch(bounds%begp:bounds%endp), & + leafc_to_litter_col(bounds%begc:bounds%endc)) + call p2c(bounds, num_soilc, filter_soilc, & + this%frootc_to_litter_patch(bounds%begp:bounds%endp), & + frootc_to_litter_col(bounds%begc:bounds%endc)) + + ! Calculate ligninNratioAve + do fc = 1,num_soilc + c = filter_soilc(fc) + if (soilbiogeochem_cwdn_col(c) > 0._r8) then + ligninNratio_cwd = CNParamsShareInst%cwd_flig * & + (soilbiogeochem_cwdc_col(c) / soilbiogeochem_cwdn_col(c)) * & + soilbiogeochem_decomp_cascade_ctransfer_col(c,i_cwdl2) + else + ligninNratio_cwd = 0._r8 + end if + this%ligninNratioAvg_col(c) = & + (ligninNratio_leaf_col(c) + ligninNratio_froot_col(c) + & + ligninNratio_cwd) / & + max(1.0e-3_r8, leafc_to_litter_col(c) + & + frootc_to_litter_col(c) + & + soilbiogeochem_decomp_cascade_ctransfer_col(c,i_cwdl2)) + end do + end if + ! this code is to calculate an exponentially-relaxed npp value for use in NDynamics code if ( trim(isotope) == 'bulk') then diff --git a/src/biogeochem/CNVegCarbonStateType.F90 b/src/biogeochem/CNVegCarbonStateType.F90 index 28cf07869b..318f00eac2 100644 --- a/src/biogeochem/CNVegCarbonStateType.F90 +++ b/src/biogeochem/CNVegCarbonStateType.F90 @@ -1077,8 +1077,6 @@ subroutine Restart ( this, bounds, ncid, flag, carbon_type, reseed_dead_plants, integer :: idata logical :: exit_spinup = .false. logical :: enter_spinup = .false. - ! flags for comparing the model and restart decomposition cascades - integer :: decomp_cascade_state, restart_file_decomp_cascade_state ! spinup state as read from restart file, for determining whether to enter or exit spinup mode. integer :: restart_file_spinup_state integer :: total_num_reseed_patch ! Total number of patches to reseed across all processors @@ -2394,7 +2392,8 @@ end subroutine ZeroDwt !----------------------------------------------------------------------- subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & num_soilc, filter_soilc, num_soilp, filter_soilp, & - soilbiogeochem_cwdc_col, soilbiogeochem_totlitc_col, soilbiogeochem_totsomc_col, & + soilbiogeochem_cwdc_col, soilbiogeochem_totlitc_col, & + soilbiogeochem_totmicc_col, soilbiogeochem_totsomc_col, & soilbiogeochem_ctrunc_col) ! ! !USES: @@ -2415,6 +2414,7 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & integer , intent(in) :: num_soilp ! number of soil patches in filter integer , intent(in) :: filter_soilp(:) ! filter for soil patches real(r8) , intent(in) :: soilbiogeochem_cwdc_col(bounds%begc:) + real(r8) , intent(in) :: soilbiogeochem_totmicc_col(bounds%begc:) real(r8) , intent(in) :: soilbiogeochem_totlitc_col(bounds%begc:) real(r8) , intent(in) :: soilbiogeochem_totsomc_col(bounds%begc:) real(r8) , intent(in) :: soilbiogeochem_ctrunc_col(bounds%begc:) @@ -2425,6 +2425,7 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & !----------------------------------------------------------------------- SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_cwdc_col) == (/bounds%endc/)), sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_totmicc_col) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_totlitc_col) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_totsomc_col) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(soilbiogeochem_ctrunc_col) == (/bounds%endc/)), sourcefile, __LINE__) @@ -2515,6 +2516,7 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & ! total ecosystem carbon, including veg but excluding cpool (TOTECOSYSC) this%totecosysc_col(c) = & soilbiogeochem_cwdc_col(c) + & + soilbiogeochem_totmicc_col(c) + & soilbiogeochem_totlitc_col(c) + & soilbiogeochem_totsomc_col(c) + & this%totvegc_col(c) @@ -2522,6 +2524,7 @@ subroutine Summary_carbonstate(this, bounds, num_allc, filter_allc, & ! total column carbon, including veg and cpool (TOTCOLC) this%totc_col(c) = this%totc_p2c_col(c) + & soilbiogeochem_cwdc_col(c) + & + soilbiogeochem_totmicc_col(c) + & soilbiogeochem_totlitc_col(c) + & soilbiogeochem_totsomc_col(c) + & soilbiogeochem_ctrunc_col(c) diff --git a/src/biogeochem/CNVegNitrogenStateType.F90 b/src/biogeochem/CNVegNitrogenStateType.F90 index a0cc51e654..f8d491faf4 100644 --- a/src/biogeochem/CNVegNitrogenStateType.F90 +++ b/src/biogeochem/CNVegNitrogenStateType.F90 @@ -4,16 +4,13 @@ module CNVegNitrogenStateType use shr_kind_mod , only : r8 => shr_kind_r8 use shr_infnan_mod , only : isnan => shr_infnan_isnan, nan => shr_infnan_nan, assignment(=) - 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 + use clm_varcon , only : spval use landunit_varcon , only : istcrop, istsoil - use clm_varctl , only : iulog, override_bgc_restart_mismatch_dump + use clm_varctl , only : iulog use clm_varctl , only : use_crop use CNSharedParamsMod , only : use_fun use decompMod , only : bounds_type use pftconMod , only : npcropmin, noveg, pftcon - use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con use abortutils , only : endrun use spmdMod , only : masterproc use LandunitType , only : lun @@ -1041,6 +1038,7 @@ subroutine Summary_nitrogenstate(this, bounds, num_allc, filter_allc, & this%totecosysn_col(c) = & soilbiogeochem_nitrogenstate_inst%cwdn_col(c) + & soilbiogeochem_nitrogenstate_inst%totlitn_col(c) + & + soilbiogeochem_nitrogenstate_inst%totmicn_col(c) + & soilbiogeochem_nitrogenstate_inst%totsomn_col(c) + & soilbiogeochem_nitrogenstate_inst%sminn_col(c) + & this%totvegn_col(c) @@ -1050,6 +1048,7 @@ subroutine Summary_nitrogenstate(this, bounds, num_allc, filter_allc, & this%totn_col(c) = this%totn_p2c_col(c) + & soilbiogeochem_nitrogenstate_inst%cwdn_col(c) + & soilbiogeochem_nitrogenstate_inst%totlitn_col(c) + & + soilbiogeochem_nitrogenstate_inst%totmicn_col(c) + & soilbiogeochem_nitrogenstate_inst%totsomn_col(c) + & soilbiogeochem_nitrogenstate_inst%sminn_col(c) + & soilbiogeochem_nitrogenstate_inst%ntrunc_col(c) diff --git a/src/biogeochem/CNVegetationFacade.F90 b/src/biogeochem/CNVegetationFacade.F90 index 98995626b0..457f5818ec 100644 --- a/src/biogeochem/CNVegetationFacade.F90 +++ b/src/biogeochem/CNVegetationFacade.F90 @@ -1083,6 +1083,10 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & soilbiogeochem_carbonflux_inst, & c13_soilbiogeochem_carbonflux_inst, & c14_soilbiogeochem_carbonflux_inst, & + soilbiogeochem_carbonstate_inst, & + c13_soilbiogeochem_carbonstate_inst, & + c14_soilbiogeochem_carbonstate_inst, & + soilbiogeochem_nitrogenstate_inst, & soilbiogeochem_nitrogenflux_inst) ! On the radiation time step, use C state variables to calculate diff --git a/src/biogeochem/EDBGCDynMod.F90 b/src/biogeochem/EDBGCDynMod.F90 index 7f30d9ef13..754f7ff4a0 100644 --- a/src/biogeochem/EDBGCDynMod.F90 +++ b/src/biogeochem/EDBGCDynMod.F90 @@ -10,7 +10,7 @@ module EDBGCDynMod use perf_mod , only : t_startf, t_stopf use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun - use SoilBiogeochemDecompCascadeConType , only : century_decomp, decomp_method + use SoilBiogeochemDecompCascadeConType , only : mimics_decomp, century_decomp, decomp_method use CNVegCarbonStateType , only : cnveg_carbonstate_type use CNVegCarbonFluxType , only : cnveg_carbonflux_type use SoilBiogeochemStateType , only : soilbiogeochem_state_type @@ -71,8 +71,8 @@ subroutine EDBGCDyn(bounds, & use CNNStateUpdate1Mod , only: NStateUpdate1 use CNNStateUpdate2Mod , only: NStateUpdate2, NStateUpdate2h use CNGapMortalityMod , only: CNGapMortality + use SoilBiogeochemDecompCascadeMIMICSMod, only: decomp_rates_mimics use SoilBiogeochemDecompCascadeBGCMod , only: decomp_rate_constants_bgc - use SoilBiogeochemCompetitionMod , only: SoilBiogeochemCompetition use SoilBiogeochemDecompMod , only: SoilBiogeochemDecomp use SoilBiogeochemLittVertTranspMod , only: SoilBiogeochemLittVertTransp use SoilBiogeochemPotentialMod , only: SoilBiogeochemPotential @@ -113,6 +113,8 @@ subroutine EDBGCDyn(bounds, & 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):: p_decomp_npool_to_din(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) ! potential flux to dissolved inorganic N + real(r8):: p_decomp_cn_gain(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_pools) ! C:N ratio of the flux gained by the receiver pool 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 @@ -180,6 +182,10 @@ subroutine EDBGCDyn(bounds, & if (decomp_method == century_decomp) then call decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & soilstate_inst, temperature_inst, ch4_inst, soilbiogeochem_carbonflux_inst) + else if (decomp_method == mimics_decomp) then + call decomp_rates_mimics(bounds, num_soilc, filter_soilc, & + soilstate_inst, temperature_inst, cnveg_carbonflux_inst, ch4_inst, & + soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst) end if ! calculate potential decomp rates and total immobilization demand (previously inlined in CNDecompAlloc) @@ -188,7 +194,9 @@ subroutine EDBGCDyn(bounds, & 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)) + p_decomp_cn_gain=p_decomp_cn_gain(begc:endc,1:nlevdecomp,1:ndecomp_pools), & + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + p_decomp_npool_to_din=p_decomp_npool_to_din(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) !-------------------------------------------- @@ -212,7 +220,8 @@ subroutine EDBGCDyn(bounds, & 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)) + pmnf_decomp_cascade=pmnf_decomp_cascade(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions), & + p_decomp_npool_to_din=p_decomp_npool_to_din(begc:endc,1:nlevdecomp,1:ndecomp_cascade_transitions)) call t_stopf('SoilBiogeochemDecomp') diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index 44ab2f8ef2..a593664217 100644 --- a/src/main/clm_instMod.F90 +++ b/src/main/clm_instMod.F90 @@ -9,8 +9,9 @@ module clm_instMod use decompMod , only : bounds_type use clm_varpar , only : ndecomp_pools, nlevdecomp_full use clm_varctl , only : use_cn, use_c13, use_c14, use_lch4, use_cndv, use_fates + use clm_varctl , only : iulog use clm_varctl , only : use_crop, snow_cover_fraction_method, paramfile - use SoilBiogeochemDecompCascadeConType , only : century_decomp, decomp_method + use SoilBiogeochemDecompCascadeConType , only : mimics_decomp, century_decomp, decomp_method use clm_varcon , only : bdsno, c13ratio, c14ratio use landunit_varcon , only : istice, istsoil use perf_mod , only : t_startf, t_stopf @@ -188,6 +189,7 @@ subroutine clm_instInit(bounds) use clm_varpar , only : nlevsno use controlMod , only : nlfilename, fsurdat use domainMod , only : ldomain + use SoilBiogeochemDecompCascadeMIMICSMod, only : init_decompcascade_mimics use SoilBiogeochemDecompCascadeBGCMod , only : init_decompcascade_bgc use SoilBiogeochemDecompCascadeContype , only : init_decomp_cascade_constants use SoilBiogeochemCompetitionMod , only : SoilBiogeochemCompetitionInit @@ -383,6 +385,9 @@ subroutine clm_instInit(bounds) if (decomp_method == century_decomp ) then call init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, & soilstate_inst ) + else if (decomp_method == mimics_decomp ) then + call init_decompcascade_mimics(bounds, soilbiogeochem_state_inst, & + soilstate_inst) end if ! Initalize soilbiogeochem carbon types diff --git a/src/main/clm_varcon.F90 b/src/main/clm_varcon.F90 index c8d84a3b07..340115fe90 100644 --- a/src/main/clm_varcon.F90 +++ b/src/main/clm_varcon.F90 @@ -122,6 +122,10 @@ module clm_varcon real(r8), public, parameter :: aquifer_water_baseline = 5000._r8 ! baseline value for water in the unconfined aquifer [mm] real(r8), public, parameter :: c_to_b = 2.0_r8 ! conversion between mass carbon and total biomass (g biomass /g C) + ! Some non-tunable conversions (may need to place elsewhere) + real(r8), public, parameter :: g_to_mg = 1.0e3_r8 ! coefficient to convert g to mg + real(r8), public, parameter :: cm3_to_m3 = 1.0e-6_r8 ! coefficient to convert cm3 to m3 + real(r8), public, parameter :: pct_to_frac = 1.0e-2_r8 ! coefficient to convert % to fraction !!! C13 real(r8), public, parameter :: preind_atm_del13c = -6.0_r8 ! preindustrial value for atmospheric del13C diff --git a/src/main/clm_varpar.F90 b/src/main/clm_varpar.F90 index e8d62453d3..6e73555af1 100644 --- a/src/main/clm_varpar.F90 +++ b/src/main/clm_varpar.F90 @@ -71,7 +71,10 @@ module clm_varpar integer, public :: i_litr_min = -9 ! min index of litter pools; overwritten in SoilBiogeochemDecompCascade*Mod integer, public :: i_litr_max = -9 ! max index of litter pools; overwritten in SoilBiogeochemDecompCascade*Mod integer, public :: i_met_lit = -9 ! index of metabolic litter pool; overwritten in SoilBiogeochemDecompCascade*Mod + integer, public :: i_cop_mic = -9 ! index of copiotrophic microbial pool; overwritten in SoilBiogeochemDecompCascade*Mod + integer, public :: i_oli_mic = -9 ! index of oligotrophic microbial pool; overwritten in SoilBiogeochemDecompCascade*Mod integer, public :: i_cwd = -9 ! index of cwd pool; overwritten in SoilBiogeochemDecompCascade*Mod + integer, public :: i_cwdl2 = -9 ! index of cwd to l2 transition; overwritten in SoilBiogeochemDecompCascade*Mod integer, public :: ndecomp_pools_max integer, public :: ndecomp_pools diff --git a/src/main/pftconMod.F90 b/src/main/pftconMod.F90 index 2ac733505f..16dd97240c 100644 --- a/src/main/pftconMod.F90 +++ b/src/main/pftconMod.F90 @@ -508,6 +508,7 @@ subroutine InitRead(this) use clm_varctl , only : paramfile, use_fates, use_flexibleCN, use_dynroot, use_biomass_heat_storage use spmdMod , only : masterproc use CLMFatesParamInterfaceMod, only : FatesReadPFTs + use SoilBiogeochemDecompCascadeConType, only : mimics_decomp, decomp_method ! ! !ARGUMENTS: class(pftcon_type) :: this @@ -765,11 +766,18 @@ subroutine InitRead(this) ! in do-loops. While executing the next few lines, we do not yet have access ! to i_litr_min, i_litr_max. this%fr_f(:,1) = this%fr_flab - this%fr_f(:,2) = this%fr_fcel - this%fr_f(:,3) = this%fr_flig this%lf_f(:,1) = this%lf_flab - this%lf_f(:,2) = this%lf_fcel - this%lf_f(:,3) = this%lf_flig + if (decomp_method == mimics_decomp) then + this%fr_f(:,2) = this%fr_fcel + this%fr_flig + this%fr_f(:,3) = 0.0_r8 + this%lf_f(:,2) = this%lf_fcel + this%lf_flig + this%lf_f(:,3) = 0.0_r8 + else + this%fr_f(:,2) = this%fr_fcel + this%fr_f(:,3) = this%fr_flig + this%lf_f(:,2) = this%lf_fcel + this%lf_f(:,3) = this%lf_flig + end if call ncd_io('leaf_long', this%leaf_long, 'read', ncid, readvar=readv, posNOTonfile=.true.) if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) diff --git a/src/main/readParamsMod.F90 b/src/main/readParamsMod.F90 index 6cad023781..31a116ab28 100644 --- a/src/main/readParamsMod.F90 +++ b/src/main/readParamsMod.F90 @@ -8,6 +8,7 @@ module readParamsMod ! ! ! USES: use clm_varctl , only : paramfile, iulog, use_fates, use_cn + use SoilBiogeochemDecompCascadeConType, only : mimics_decomp, century_decomp, decomp_method use spmdMod , only : masterproc use fileutils , only : getfil use ncdio_pio , only : ncd_pio_closefile, ncd_pio_openfile @@ -37,6 +38,7 @@ subroutine readParameters (nutrient_competition_method, photosyns_inst) use SoilBiogeochemLittVertTranspMod , only : readSoilBiogeochemLittVertTranspParams => readParams use SoilBiogeochemPotentialMod , only : readSoilBiogeochemPotentialParams => readParams use SoilBiogeochemDecompMod , only : readSoilBiogeochemDecompParams => readParams + use SoilBiogeochemDecompCascadeMIMICSMod, only : readSoilBiogeochemDecompMimicsParams => readParams use SoilBiogeochemDecompCascadeBGCMod , only : readSoilBiogeochemDecompBgcParams => readParams use ch4Mod , only : readCH4Params => readParams use LunaMod , only : readParams_Luna => readParams @@ -98,7 +100,11 @@ subroutine readParameters (nutrient_competition_method, photosyns_inst) ! if (use_cn .or. use_fates) then call readSoilBiogeochemCompetitionParams(ncid) - call readSoilBiogeochemDecompBgcParams(ncid) + if (decomp_method == mimics_decomp) then + call readSoilBiogeochemDecompMimicsParams(ncid) + else if (decomp_method == century_decomp) then + call readSoilBiogeochemDecompBgcParams(ncid) + end if call readSoilBiogeochemDecompParams(ncid) call readSoilBiogeochemLittVertTranspParams(ncid) call readSoilBiogeochemNitrifDenitrifParams(ncid) diff --git a/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 b/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 index cb8f3beef4..343c30eb42 100644 --- a/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemCarbonFluxType.F90 @@ -25,11 +25,15 @@ module SoilBiogeochemCarbonFluxType ! decomposition fluxes real(r8), pointer :: decomp_cpools_sourcesink_col (:,:,:) ! change in decomposing c pools. Used to update concentrations concurrently with vertical transport (gC/m3/timestep) + real(r8), pointer :: c_overflow_vr (:,:,:) ! vertically-resolved C rejected by microbes that cannot process it (gC/m3/s) real(r8), pointer :: decomp_cascade_hr_vr_col (:,:,:) ! vertically-resolved het. resp. from decomposing C pools (gC/m3/s) real(r8), pointer :: decomp_cascade_hr_col (:,:) ! vertically-integrated (diagnostic) het. resp. from decomposing C pools (gC/m2/s) real(r8), pointer :: decomp_cascade_ctransfer_vr_col (:,:,:) ! vertically-resolved C transferred along deomposition cascade (gC/m3/s) real(r8), pointer :: decomp_cascade_ctransfer_col (:,:) ! vertically-integrated (diagnostic) C transferred along decomposition cascade (gC/m2/s) - real(r8), pointer :: decomp_k_col (:,:,:) ! rate constant for decomposition (1./sec) + real(r8), pointer :: cn_col (:,:) ! (gC/gN) C:N ratio by pool + real(r8), pointer :: rf_decomp_cascade_col (:,:,:) ! (frac) respired fraction in decomposition step + real(r8), pointer :: pathfrac_decomp_cascade_col (:,:,:) ! (frac) what fraction of C passes from donor to receiver pool through a given transition + real(r8), pointer :: decomp_k_col (:,:,:) ! rate coefficient for decomposition (1./sec) real(r8), pointer :: hr_vr_col (:,:) ! (gC/m3/s) total vertically-resolved het. resp. from decomposing C pools real(r8), pointer :: o_scalar_col (:,:) ! fraction by which decomposition is limited by anoxia real(r8), pointer :: w_scalar_col (:,:) ! fraction by which decomposition is limited by moisture availability @@ -43,6 +47,7 @@ module SoilBiogeochemCarbonFluxType real(r8), pointer :: fphr_col (:,:) ! fraction of potential heterotrophic respiration real(r8), pointer :: hr_col (:) ! (gC/m2/s) total heterotrophic respiration + real(r8), pointer :: michr_col (:) ! (gC/m2/s) microbial heterotrophic respiration real(r8), pointer :: cwdhr_col (:) ! (gC/m2/s) coarse woody debris heterotrophic respiration real(r8), pointer :: lithr_col (:) ! (gC/m2/s) litter heterotrophic respiration real(r8), pointer :: somhr_col (:) ! (gC/m2/s) soil organic matter heterotrophic res @@ -109,6 +114,9 @@ subroutine InitAllocate(this, bounds) allocate(this%decomp_cpools_sourcesink_col(begc:endc,1:nlevdecomp_full,1:ndecomp_pools)) this%decomp_cpools_sourcesink_col(:,:,:)= nan + allocate(this%c_overflow_vr(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)) + this%c_overflow_vr(:,:,:) = nan + allocate(this%decomp_cascade_hr_vr_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)) this%decomp_cascade_hr_vr_col(:,:,:)= spval @@ -121,7 +129,16 @@ subroutine InitAllocate(this, bounds) allocate(this%decomp_cascade_ctransfer_col(begc:endc,1:ndecomp_cascade_transitions)) this%decomp_cascade_ctransfer_col(:,:)= nan - allocate(this%decomp_k_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)) + allocate(this%cn_col(begc:endc,1:ndecomp_pools)) + this%cn_col(:,:)= spval + + allocate(this%rf_decomp_cascade_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)) + this%rf_decomp_cascade_col(:,:,:) = nan + + allocate(this%pathfrac_decomp_cascade_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)) + this%pathfrac_decomp_cascade_col(:,:,:) = nan + + allocate(this%decomp_k_col(begc:endc,1:nlevdecomp_full,1:ndecomp_pools)) this%decomp_k_col(:,:,:)= spval allocate(this%decomp_cpools_leached_col(begc:endc,1:ndecomp_pools)) @@ -131,6 +148,7 @@ subroutine InitAllocate(this, bounds) this%decomp_cpools_transport_tendency_col(:,:,:)= nan allocate(this%hr_col (begc:endc)) ; this%hr_col (:) = nan + allocate(this%michr_col (begc:endc)) ; this%michr_col (:) = nan allocate(this%cwdhr_col (begc:endc)) ; this%cwdhr_col (:) = nan allocate(this%lithr_col (begc:endc)) ; this%lithr_col (:) = nan allocate(this%somhr_col (begc:endc)) ; this%somhr_col (:) = nan @@ -206,6 +224,11 @@ subroutine InitHistory(this, bounds, carbon_type) avgflag='A', long_name='total heterotrophic respiration', & ptr_col=this%hr_col) + this%michr_col(begc:endc) = spval + call hist_addfld1d (fname='MICC_HR', units='gC/m^2/s', & + avgflag='A', long_name='microbial C heterotrophic respiration', & + ptr_col=this%michr_col) + this%cwdhr_col(begc:endc) = spval call hist_addfld1d (fname='CWDC_HR', units='gC/m^2/s', & avgflag='A', long_name='cwd C heterotrophic respiration', & @@ -247,6 +270,8 @@ subroutine InitHistory(this, bounds, carbon_type) this%decomp_cascade_hr_vr_col(begc:endc,:,:) = spval this%decomp_cascade_ctransfer_col(begc:endc,:) = spval this%decomp_cascade_ctransfer_vr_col(begc:endc,:,:) = spval + this%pathfrac_decomp_cascade_col(begc:endc,:,:) = spval + this%rf_decomp_cascade_col(begc:endc,:,:) = spval do l = 1, ndecomp_cascade_transitions ! output the vertically integrated fluxes only as default @@ -323,6 +348,32 @@ subroutine InitHistory(this, bounds, carbon_type) call hist_addfld_decomp (fname=fieldname, units='gC/m^3/s', type2d='levdcmp', & avgflag='A', long_name=longname, & ptr_col=data2dptr, default='inactive') + + ! pathfrac_decomp_cascade_col and rf_decomp_cascade_col needed + ! when using Newton-Krylov to spin up the decomposition + data2dptr => this%rf_decomp_cascade_col(:,:,l) + fieldname = & + trim(decomp_cascade_con%decomp_pool_name_short(decomp_cascade_con%cascade_donor_pool(l)))//'_RESP_FRAC_'//& + trim(decomp_cascade_con%decomp_pool_name_short(decomp_cascade_con%cascade_receiver_pool(l)))& + //trim(vr_suffix) + longname = 'respired from '//& + trim(decomp_cascade_con%decomp_pool_name_long(decomp_cascade_con%cascade_donor_pool(l)))//' to '// & + trim(decomp_cascade_con%decomp_pool_name_long(decomp_cascade_con%cascade_receiver_pool(l))) + call hist_addfld_decomp (fname=fieldname, units='fraction', type2d='levdcmp', & + avgflag='A', long_name=longname, & + ptr_col=data2dptr, default='inactive') + + data2dptr => this%pathfrac_decomp_cascade_col(:,:,l) + fieldname = & + trim(decomp_cascade_con%decomp_pool_name_short(decomp_cascade_con%cascade_donor_pool(l)))//'_PATHFRAC_'//& + trim(decomp_cascade_con%decomp_pool_name_short(decomp_cascade_con%cascade_receiver_pool(l)))& + //trim(vr_suffix) + longname = 'PATHFRAC from '//& + trim(decomp_cascade_con%decomp_pool_name_long(decomp_cascade_con%cascade_donor_pool(l)))//' to '// & + trim(decomp_cascade_con%decomp_pool_name_long(decomp_cascade_con%cascade_receiver_pool(l))) + call hist_addfld_decomp (fname=fieldname, units='fraction', type2d='levdcmp', & + avgflag='A', long_name=longname, & + ptr_col=data2dptr, default='inactive') endif end if @@ -390,6 +441,11 @@ subroutine InitHistory(this, bounds, carbon_type) avgflag='A', long_name='C13 total heterotrophic respiration', & ptr_col=this%hr_col) + this%michr_col(begc:endc) = spval + call hist_addfld1d (fname='C13_MICC_HR', units='gC13/m^2/s', & + avgflag='A', long_name='C13 microbial heterotrophic respiration', & + ptr_col=this%michr_col, default='inactive') + this%cwdhr_col(begc:endc) = spval call hist_addfld1d (fname='C13_CWDC_HR', units='gC/m^2/s', & avgflag='A', long_name='C13 cwd C heterotrophic respiration', & @@ -464,6 +520,11 @@ subroutine InitHistory(this, bounds, carbon_type) avgflag='A', long_name='C14 total heterotrophic respiration', & ptr_col=this%hr_col) + this%michr_col(begc:endc) = spval + call hist_addfld1d (fname='C14_MICC_HR', units='gC13/m^2/s', & + avgflag='A', long_name='C14 microbial heterotrophic respiration', & + ptr_col=this%michr_col, default='inactive') + this%cwdhr_col(begc:endc) = spval call hist_addfld1d (fname='C14_CWDC_HR', units='gC/m^2/s', & avgflag='A', long_name='C14 cwd C heterotrophic respiration', & @@ -670,10 +731,12 @@ subroutine SetValues ( this, num_column, filter_column, value_column) do fi = 1,num_column i = filter_column(fi) this%decomp_cascade_hr_col(i,l) = value_column + this%c_overflow_vr(i,j,l) = value_column this%decomp_cascade_hr_vr_col(i,j,l) = value_column this%decomp_cascade_ctransfer_col(i,l) = value_column this%decomp_cascade_ctransfer_vr_col(i,j,l) = value_column - this%decomp_k_col(i,j,l) = value_column + this%pathfrac_decomp_cascade_col(i,j,l) = value_column + this%rf_decomp_cascade_col(i,j,l) = value_column end do end do end do @@ -682,12 +745,14 @@ subroutine SetValues ( this, num_column, filter_column, value_column) do fi = 1,num_column i = filter_column(fi) this%decomp_cpools_leached_col(i,k) = value_column + this%cn_col(i,k) = value_column end do do j = 1, nlevdecomp_full do fi = 1,num_column i = filter_column(fi) this%decomp_cpools_transport_tendency_col(i,j,k) = value_column this%decomp_cpools_sourcesink_col(i,j,k) = value_column + this%decomp_k_col(i,j,k) = value_column end do end do end do @@ -707,6 +772,7 @@ subroutine SetValues ( this, num_column, filter_column, value_column) this%somhr_col(i) = value_column this%lithr_col(i) = value_column this%cwdhr_col(i) = value_column + this%michr_col(i) = value_column this%soilc_change_col(i) = value_column end do @@ -826,11 +892,24 @@ subroutine Summary(this, bounds, num_soilc, filter_soilc) end do end associate + ! microbial heterotrophic respiration (MICHR) + associate(is_microbe => decomp_cascade_con%is_microbe) ! TRUE => pool is a microbial pool + do k = 1, ndecomp_cascade_transitions + if ( is_microbe(decomp_cascade_con%cascade_donor_pool(k)) ) then + do fc = 1,num_soilc + c = filter_soilc(fc) + this%michr_col(c) = this%michr_col(c) + this%decomp_cascade_hr_col(c,k) + end do + end if + end do + end associate + ! total heterotrophic respiration (HR) do fc = 1,num_soilc c = filter_soilc(fc) this%hr_col(c) = & + this%michr_col(c) + & this%cwdhr_col(c) + & this%lithr_col(c) + & this%somhr_col(c) diff --git a/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 b/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 index a8f312a55c..7ef4265a4c 100644 --- a/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemCarbonStateType.F90 @@ -30,6 +30,7 @@ module SoilBiogeochemCarbonStateType ! summary (diagnostic) state variables, not involved in mass balance real(r8), pointer :: ctrunc_col (:) ! (gC/m2) column-level sink for C truncation + real(r8), pointer :: totmicc_col (:) ! (gC/m2) total microbial carbon real(r8), pointer :: totlitc_col (:) ! (gC/m2) total litter carbon real(r8), pointer :: totlitc_1m_col (:) ! (gC/m2) total litter carbon to 1 meter real(r8), pointer :: totsomc_col (:) ! (gC/m2) total soil organic matter carbon @@ -110,6 +111,7 @@ subroutine InitAllocate(this, bounds) if ( .not. use_fates ) then allocate(this%cwdc_col (begc :endc)) ; this%cwdc_col (:) = nan endif + allocate(this%totmicc_col (begc :endc)) ; this%totmicc_col (:) = nan 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 @@ -183,6 +185,11 @@ subroutine InitHistory(this, bounds, carbon_type) endif end do + this%totmicc_col(begc:endc) = spval + call hist_addfld1d (fname='TOTMICC', units='gC/m^2', & + avgflag='A', long_name='total microbial carbon', & + ptr_col=this%totmicc_col) + this%totlitc_col(begc:endc) = spval call hist_addfld1d (fname='TOTLITC', units='gC/m^2', & avgflag='A', long_name='total litter carbon', & @@ -253,6 +260,11 @@ subroutine InitHistory(this, bounds, carbon_type) ptr_col=data1dptr, default='inactive') end do + this%totmicc_col(begc:endc) = spval + call hist_addfld1d (fname='C13_TOTMICC', units='gC/m^2', & + avgflag='A', long_name='C13 total microbial carbon', & + ptr_col=this%totmicc_col) + this%totlitc_col(begc:endc) = spval call hist_addfld1d (fname='C13_TOTLITC', units='gC13/m^2', & avgflag='A', long_name='C13 total litter carbon', & @@ -328,6 +340,11 @@ subroutine InitHistory(this, bounds, carbon_type) endif end do + this%totmicc_col(begc:endc) = spval + call hist_addfld1d (fname='C14_TOTMICC', units='gC/m^2', & + avgflag='A', long_name='C14 total microbial carbon', & + ptr_col=this%totmicc_col) + this%totlitc_col(begc:endc) = spval call hist_addfld1d (fname='C14_TOTLITC', units='gC14/m^2', & avgflag='A', long_name='C14 total litter carbon', & @@ -447,6 +464,7 @@ subroutine InitCold(this, bounds, ratio, c12_soilbiogeochem_carbonstate_inst) this%cwdc_col(c) = 0._r8 end if this%ctrunc_col(c) = 0._r8 + this%totmicc_col(c) = 0._r8 this%totlitc_col(c) = 0._r8 this%totsomc_col(c) = 0._r8 this%totlitc_1m_col(c) = 0._r8 @@ -510,8 +528,6 @@ subroutine Restart ( this, bounds, ncid, flag, carbon_type, totvegc_col, c12_so integer :: idata logical :: exit_spinup = .false. logical :: enter_spinup = .false. - ! flags for comparing the model and restart decomposition cascades - integer :: decomp_cascade_state, restart_file_decomp_cascade_state !------------------------------------------------------------------------ if (carbon_type == 'c12') then @@ -719,6 +735,7 @@ subroutine SetValues ( this, num_column, filter_column, value_column) this%cwdc_col(i) = value_column end if this%ctrunc_col(i) = value_column + this%totmicc_col(i) = value_column this%totlitc_col(i) = value_column this%totlitc_1m_col(i) = value_column this%totsomc_col(i) = value_column @@ -887,6 +904,20 @@ subroutine Summary(this, bounds, num_allc, filter_allc) end do end if + ! total microbial carbon (TOTMICC) + do fc = 1,num_allc + c = filter_allc(fc) + this%totmicc_col(c) = 0._r8 + end do + do l = 1, ndecomp_pools + if ( decomp_cascade_con%is_microbe(l) ) then + do fc = 1,num_allc + c = filter_allc(fc) + this%totmicc_col(c) = this%totmicc_col(c) + this%decomp_cpools_col(c,l) + end do + endif + end do + ! total litter carbon (TOTLITC) do fc = 1,num_allc c = filter_allc(fc) diff --git a/src/soilbiogeochem/SoilBiogeochemCompetitionMod.F90 b/src/soilbiogeochem/SoilBiogeochemCompetitionMod.F90 index 6f5580bb8b..eec9e00f80 100644 --- a/src/soilbiogeochem/SoilBiogeochemCompetitionMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemCompetitionMod.F90 @@ -165,7 +165,8 @@ subroutine SoilBiogeochemCompetitionInit ( bounds) end subroutine SoilBiogeochemCompetitionInit !----------------------------------------------------------------------- - subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, filter_soilp, waterstatebulk_inst, & + subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, filter_soilp, & + p_decomp_cn_gain, pmnf_decomp_cascade, waterstatebulk_inst, & waterfluxbulk_inst, temperature_inst,soilstate_inst, & cnveg_state_inst,cnveg_carbonstate_inst, & cnveg_carbonflux_inst,cnveg_nitrogenstate_inst,cnveg_nitrogenflux_inst, & @@ -176,11 +177,13 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, ! !USES: use clm_varctl , only: cnallocate_carbon_only, iulog use clm_varpar , only: nlevdecomp, ndecomp_cascade_transitions + use clm_varpar , only: i_cop_mic, i_oli_mic use clm_varcon , only: nitrif_n2o_loss_frac use CNSharedParamsMod, only: use_fun use CNFUNMod , only: CNFUN use subgridAveMod , only: p2c use perf_mod , only : t_startf, t_stopf + use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, mimics_decomp, decomp_method ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -188,6 +191,8 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, 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 + real(r8) , intent(in) :: pmnf_decomp_cascade(bounds%begc:,1:,1:) ! potential mineral N flux from one pool to another (gN/m3/s) + real(r8) , intent(in) :: p_decomp_cn_gain(bounds%begc:,1:,1:) ! C:N ratio of the flux gained by the receiver pool type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(in) :: waterfluxbulk_inst type(temperature_type) , intent(in) :: temperature_inst @@ -205,9 +210,11 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, ! ! ! !LOCAL VARIABLES: - integer :: c,p,l,pi,j ! indices + integer :: c,p,l,pi,j,k ! indices integer :: fc ! filter column index logical :: local_use_fun ! local version of use_fun + real(r8) :: amnf_immob_vr ! actual mineral N flux from immobilization (gN/m3/s) + real(r8) :: n_deficit_vr ! microbial N deficit, vertically resolved (gN/m3/s) real(r8) :: compet_plant_no3 ! (unitless) relative compettiveness of plants for NO3 real(r8) :: compet_plant_nh4 ! (unitless) relative compettiveness of plants for NH4 real(r8) :: compet_decomp_no3 ! (unitless) relative competitiveness of immobilizers for NO3 @@ -247,6 +254,9 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, smin_nh4_vr => soilbiogeochem_nitrogenstate_inst%smin_nh4_vr_col , & ! Input: [real(r8) (:,:) ] (gN/m3) soil mineral NH4 smin_no3_vr => soilbiogeochem_nitrogenstate_inst%smin_no3_vr_col , & ! Input: [real(r8) (:,:) ] (gN/m3) soil mineral NO3 + c_overflow_vr => soilbiogeochem_carbonflux_inst%c_overflow_vr , & ! Output: [real(r8) (:,:,:)] (gC/m3/s) vertically-resolved C rejected by microbes that cannot process it + cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Input: [integer (:) ] which pool is C added to for a given decomposition step + pot_f_nit_vr => soilbiogeochem_nitrogenflux_inst%pot_f_nit_vr_col , & ! Input: [real(r8) (:,:) ] (gN/m3/s) potential soil nitrification flux pot_f_denit_vr => soilbiogeochem_nitrogenflux_inst%pot_f_denit_vr_col , & ! Input: [real(r8) (:,:) ] (gN/m3/s) potential soil denitrification flux f_nit_vr => soilbiogeochem_nitrogenflux_inst%f_nit_vr_col , & ! Output: [real(r8) (:,:) ] (gN/m3/s) soil nitrification flux @@ -600,7 +610,13 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, end if - + if (decomp_method == mimics_decomp) then + ! turn off fpi for MIMICS and only lets plants + ! take up available mineral nitrogen. + ! TODO slevis: -ve or tiny sminn_vr could cause problems + fpi_nh4_vr(c,j) = 1.0_r8 + actual_immob_nh4_vr(c,j) = potential_immob_vr(c,j) + end if if(.not.local_use_fun)then sum_no3_demand(c,j) = (plant_ndemand(c)*nuptake_prof(c,j)-smin_nh4_to_plant_vr(c,j)) + & @@ -693,9 +709,15 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, end if end if - - + if (decomp_method == mimics_decomp) then + ! turn off fpi for MIMICS and only lets plants + ! take up available mineral nitrogen. + ! TODO slevis: -ve or tiny sminn_vr could cause problems + fpi_no3_vr(c,j) = 1.0_r8 - fpi_nh4_vr(c,j) ! => 0 + actual_immob_no3_vr(c,j) = potential_immob_vr(c,j) - & + actual_immob_nh4_vr(c,j) ! => 0 + end if ! n2o emissions: n2o from nitr is const fraction, n2o from denitr is calculated in nitrif_denitrif f_n2o_nit_vr(c,j) = f_nit_vr(c,j) * nitrif_n2o_loss_frac @@ -791,6 +813,38 @@ subroutine SoilBiogeochemCompetition (bounds, num_soilc, filter_soilc,num_soilp, end if + if (decomp_method == mimics_decomp) then + do j = 1, nlevdecomp + do fc=1,num_soilc + c = filter_soilc(fc) + + do k = 1, ndecomp_cascade_transitions + if (cascade_receiver_pool(k) == i_cop_mic .or. & + cascade_receiver_pool(k) == i_oli_mic) then + sum_ndemand_vr(c,j) = sum_no3_demand_scaled(c,j) + & + sum_nh4_demand_scaled(c,j) + if (pmnf_decomp_cascade(c,j,k) > 0.0_r8 .and. & + sum_ndemand_vr(c,j) > 0.0_r8) then + amnf_immob_vr = (sminn_vr(c,j) / dt) * & + (pmnf_decomp_cascade(c,j,k) / & + sum_ndemand_vr(c,j)) + n_deficit_vr = pmnf_decomp_cascade(c,j,k) - & + amnf_immob_vr + c_overflow_vr(c,j,k) = & + n_deficit_vr * p_decomp_cn_gain(c,j,cascade_receiver_pool(k)) + else ! not pmnf and sum_ndemand > 0 + c_overflow_vr(c,j,k) = 0.0_r8 + end if + else ! not microbes receiving + c_overflow_vr(c,j,k) = 0.0_r8 + end if + end do + end do + end do + else ! not mimics_decomp + c_overflow_vr(:,:,:) = 0.0_r8 + end if + if(.not.local_use_fun)then ! give plants a second pass to see if there is any mineral N left over with which to satisfy residual N demand. ! first take frm nh4 pool; then take from no3 pool diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 index 6749b8ebc4..99a736a382 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemDecompCascadeBGCMod.F90 @@ -10,7 +10,7 @@ module SoilBiogeochemDecompCascadeBGCMod use shr_const_mod , only : SHR_CONST_TKFRZ use shr_log_mod , only : errMsg => shr_log_errMsg use clm_varpar , only : nlevdecomp, ndecomp_pools_max - use clm_varpar , only : i_litr_min, i_litr_max, i_met_lit, i_cwd + use clm_varpar , only : i_litr_min, i_litr_max, i_met_lit, i_cwd, i_cwdl2 use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_fates use clm_varcon , only : zsoi use decompMod , only : bounds_type @@ -48,6 +48,31 @@ module SoilBiogeochemDecompCascadeBGCMod integer, private :: i_cel_lit ! index of cellulose litter pool integer, private :: i_lig_lit ! index of lignin litter pool + real(r8), private :: cwd_fcel + real(r8), private :: rf_l1s1 + real(r8), private :: rf_l2s1 + real(r8), private :: rf_l3s2 + real(r8), private :: rf_s2s1 + real(r8), private :: rf_s2s3 + real(r8), private :: rf_s3s1 + real(r8), private :: rf_cwdl3 + real(r8), private, allocatable :: rf_s1s2(:,:) + real(r8), private, allocatable :: rf_s1s3(:,:) + real(r8), private, allocatable :: f_s1s2(:,:) + real(r8), private, allocatable :: f_s1s3(:,:) + real(r8), private :: f_s2s1 + real(r8), private :: f_s2s3 + + integer, private :: i_l1s1 + integer, private :: i_l2s1 + integer, private :: i_l3s2 + integer, private :: i_s1s2 + integer, private :: i_s1s3 + integer, private :: i_s2s1 + integer, private :: i_s2s3 + integer, private :: i_s3s1 + integer, private :: i_cwdl3 + type, private :: params_type real(r8):: cn_s1_bgc !C:N for SOM 1 real(r8):: cn_s2_bgc !C:N for SOM 2 @@ -61,7 +86,6 @@ module SoilBiogeochemDecompCascadeBGCMod real(r8):: rf_s2s3_bgc real(r8):: rf_s3s1_bgc - real(r8):: rf_cwdl2_bgc real(r8):: rf_cwdl3_bgc real(r8):: tau_l1_bgc ! 1/turnover time of litter 1 from Century (l/18.5) (1/yr) @@ -69,16 +93,11 @@ module SoilBiogeochemDecompCascadeBGCMod real(r8):: tau_s1_bgc ! 1/turnover time of SOM 1 from Century (1/7.3) (1/yr) real(r8):: tau_s2_bgc ! 1/turnover time of SOM 2 from Century (1/0.2) (1/yr) real(r8):: tau_s3_bgc ! 1/turnover time of SOM 3 from Century (1/0.0045) (1/yr) - real(r8):: tau_cwd_bgc ! corrected fragmentation rate constant CWD, century leaves wood decomposition rates open, within range of 0 - 0.5 yr^-1 (1/0.3) (1/yr) real(r8) :: cwd_fcel_bgc !cellulose fraction for CWD - real(r8) :: cwd_flig - - real(r8) :: minpsi_bgc !minimum soil water potential for heterotrophic resp - real(r8) :: maxpsi_bgc !maximum soil water potential for heterotrophic resp - real(r8), allocatable :: initial_Cstocks(:) ! Initial Carbon stocks for a cold-start - real(r8) :: initial_Cstocks_depth ! Soil depth for initial Carbon stocks for a cold-start + real(r8), allocatable :: bgc_initial_Cstocks(:) ! Initial Carbon stocks for a cold-start + real(r8) :: bgc_initial_Cstocks_depth ! Soil depth for initial Carbon stocks for a cold-start end type params_type ! @@ -111,120 +130,95 @@ subroutine readParams ( ncid ) !----------------------------------------------------------------------- ! Read off of netcdf file - tString='tau_l1' + tString='bgc_tau_l1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%tau_l1_bgc=tempr - tString='tau_l2_l3' + tString='bgc_tau_l2_l3' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%tau_l2_l3_bgc=tempr - tString='tau_s1' + tString='bgc_tau_s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%tau_s1_bgc=tempr - tString='tau_s2' + tString='bgc_tau_s2' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%tau_s2_bgc=tempr - tString='tau_s3' + tString='bgc_tau_s3' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%tau_s3_bgc=tempr - tString='tau_cwd_bgc' - call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) - if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%tau_cwd_bgc=tempr - - tString='cn_s1_bgc' + tString='bgc_cn_s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%cn_s1_bgc=tempr - tString='cn_s2_bgc' + tString='bgc_cn_s2' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%cn_s2_bgc=tempr - tString='cn_s3_bgc' + tString='bgc_cn_s3' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%cn_s3_bgc=tempr - tString='rf_l1s1_bgc' + tString='bgc_rf_l1s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_l1s1_bgc=tempr - tString='rf_l2s1_bgc' + tString='bgc_rf_l2s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_l2s1_bgc=tempr - tString='rf_l3s2_bgc' + tString='bgc_rf_l3s2' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_l3s2_bgc=tempr - tString='rf_s2s1_bgc' + tString='bgc_rf_s2s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_s2s1_bgc=tempr - tString='rf_s2s3_bgc' + tString='bgc_rf_s2s3' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_s2s3_bgc=tempr - tString='rf_s3s1_bgc' + tString='bgc_rf_s3s1' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_s3s1_bgc=tempr - tString='rf_cwdl2_bgc' - call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) - if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%rf_cwdl2_bgc=tempr - - tString='rf_cwdl3_bgc' + tString='bgc_rf_cwdl3' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%rf_cwdl3_bgc=tempr - tString='cwd_fcel' + tString='bgc_cwd_fcel' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) params_inst%cwd_fcel_bgc=tempr - tString='minpsi_hr' - call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) + allocate(params_inst%bgc_initial_Cstocks(ndecomp_pools_max)) + tString='bgc_initial_Cstocks' + call ncd_io(trim(tString), params_inst%bgc_initial_Cstocks(:), 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%minpsi_bgc=tempr - tString='maxpsi_hr' + tString='bgc_initial_Cstocks_depth' call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%maxpsi_bgc=tempr - - tString='cwd_flig' - call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) - if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%cwd_flig=tempr - - allocate(params_inst%initial_Cstocks(ndecomp_pools_max)) - tString='initial_Cstocks_bgc' - call ncd_io(trim(tString), params_inst%initial_Cstocks(:), 'read', ncid, readvar=readv) - if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - - tString='initial_Cstocks_depth_bgc' - call ncd_io(trim(tString),tempr, 'read', ncid, readvar=readv) - if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) - params_inst%initial_Cstocks_depth=tempr + params_inst%bgc_initial_Cstocks_depth=tempr end subroutine readParams @@ -244,40 +238,9 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i ! ! !LOCAL VARIABLES !-- properties of each decomposing pool - real(r8) :: rf_l1s1 - real(r8) :: rf_l2s1 - real(r8) :: rf_l3s2 - !real(r8) :: rf_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) - !real(r8) :: rf_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) - real(r8), allocatable :: rf_s1s2(:,:) - real(r8), allocatable :: rf_s1s3(:,:) - real(r8) :: rf_s2s1 - real(r8) :: rf_s2s3 - real(r8) :: rf_s3s1 - real(r8) :: rf_cwdl2 - real(r8) :: rf_cwdl3 - real(r8) :: cwd_fcel - real(r8) :: cwd_flig real(r8) :: cn_s1 real(r8) :: cn_s2 real(r8) :: cn_s3 - !real(r8) :: f_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) - !real(r8) :: f_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) - real(r8), allocatable :: f_s1s2(:,:) - real(r8), allocatable :: f_s1s3(:,:) - real(r8) :: f_s2s1 - real(r8) :: f_s2s3 - - integer :: i_l1s1 - integer :: i_l2s1 - integer :: i_l3s2 - integer :: i_s1s2 - integer :: i_s1s3 - integer :: i_s2s1 - integer :: i_s2s3 - integer :: i_s3s1 - integer :: i_cwdl2 - integer :: i_cwdl3 real(r8):: speedup_fac ! acceleration factor, higher when vertsoilc = .true. integer :: c, j ! indices @@ -285,9 +248,6 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i !----------------------------------------------------------------------- associate( & - rf_decomp_cascade => soilbiogeochem_state_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) - pathfrac_decomp_cascade => soilbiogeochem_state_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] what fraction of C leaving a given pool passes through a given transition (frac) - cellsand => soilstate_inst%cellsand_col , & ! Input: [real(r8) (:,:) ] column 3D sand cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Output: [integer (:) ] which pool is C taken from for a given decomposition step @@ -325,12 +285,10 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i rf_s2s3 = params_inst%rf_s2s3_bgc rf_s3s1 = params_inst%rf_s3s1_bgc - rf_cwdl2 = params_inst%rf_cwdl2_bgc rf_cwdl3 = params_inst%rf_cwdl3_bgc ! set the cellulose and lignin fractions for coarse woody debris cwd_fcel = params_inst%cwd_fcel_bgc - cwd_flig = params_inst%cwd_flig ! set path fractions f_s2s1 = 0.42_r8/(0.45_r8) @@ -346,7 +304,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i rf_s1s3(c,j) = t end do end do - initial_stock_soildepth = params_inst%initial_Cstocks_depth + initial_stock_soildepth = params_inst%bgc_initial_Cstocks_depth !------------------- list of pools and their attributes ------------ i_litr_min = 1 @@ -360,7 +318,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_met_lit) = .false. is_cwd(i_met_lit) = .false. initial_cn_ratio(i_met_lit) = 90._r8 - initial_stock(i_met_lit) = params_inst%initial_Cstocks(i_met_lit) + initial_stock(i_met_lit) = params_inst%bgc_initial_Cstocks(i_met_lit) is_metabolic(i_met_lit) = .true. is_cellulose(i_met_lit) = .false. is_lignin(i_met_lit) = .false. @@ -375,7 +333,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_cel_lit) = .false. is_cwd(i_cel_lit) = .false. initial_cn_ratio(i_cel_lit) = 90._r8 - initial_stock(i_cel_lit) = params_inst%initial_Cstocks(i_cel_lit) + initial_stock(i_cel_lit) = params_inst%bgc_initial_Cstocks(i_cel_lit) is_metabolic(i_cel_lit) = .false. is_cellulose(i_cel_lit) = .true. is_lignin(i_cel_lit) = .false. @@ -390,7 +348,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_lig_lit) = .false. is_cwd(i_lig_lit) = .false. initial_cn_ratio(i_lig_lit) = 90._r8 - initial_stock(i_lig_lit) = params_inst%initial_Cstocks(i_lig_lit) + initial_stock(i_lig_lit) = params_inst%bgc_initial_Cstocks(i_lig_lit) is_metabolic(i_lig_lit) = .false. is_cellulose(i_lig_lit) = .false. is_lignin(i_lig_lit) = .true. @@ -415,7 +373,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_act_som) = .true. is_cwd(i_act_som) = .false. initial_cn_ratio(i_act_som) = cn_s1 - initial_stock(i_act_som) = params_inst%initial_Cstocks(i_act_som) + initial_stock(i_act_som) = params_inst%bgc_initial_Cstocks(i_act_som) is_metabolic(i_act_som) = .false. is_cellulose(i_act_som) = .false. is_lignin(i_act_som) = .false. @@ -430,7 +388,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_slo_som) = .true. is_cwd(i_slo_som) = .false. initial_cn_ratio(i_slo_som) = cn_s2 - initial_stock(i_slo_som) = params_inst%initial_Cstocks(i_slo_som) + initial_stock(i_slo_som) = params_inst%bgc_initial_Cstocks(i_slo_som) is_metabolic(i_slo_som) = .false. is_cellulose(i_slo_som) = .false. is_lignin(i_slo_som) = .false. @@ -445,7 +403,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_pas_som) = .true. is_cwd(i_pas_som) = .false. initial_cn_ratio(i_pas_som) = cn_s3 - initial_stock(i_pas_som) = params_inst%initial_Cstocks(i_pas_som) + initial_stock(i_pas_som) = params_inst%bgc_initial_Cstocks(i_pas_som) is_metabolic(i_pas_som) = .false. is_cellulose(i_pas_som) = .false. is_lignin(i_pas_som) = .false. @@ -462,7 +420,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i is_soil(i_cwd) = .false. is_cwd(i_cwd) = .true. initial_cn_ratio(i_cwd) = 90._r8 - initial_stock(i_cwd) = params_inst%initial_Cstocks(i_cwd) + initial_stock(i_cwd) = params_inst%bgc_initial_Cstocks(i_cwd) is_metabolic(i_cwd) = .false. is_cellulose(i_cwd) = .false. is_lignin(i_cwd) = .false. @@ -477,7 +435,7 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i spinup_factor(i_lig_lit) = 1._r8 !CWD if (.not. use_fates) then - spinup_factor(i_cwd) = max(1._r8, (speedup_fac * params_inst%tau_cwd_bgc / 2._r8 )) + spinup_factor(i_cwd) = max(1._r8, (speedup_fac * CNParamsShareInst%tau_cwd / 2._r8 )) end if !som1 spinup_factor(i_act_som) = 1._r8 @@ -493,81 +451,57 @@ subroutine init_decompcascade_bgc(bounds, soilbiogeochem_state_inst, soilstate_i !---------------- list of transitions and their time-independent coefficients ---------------! i_l1s1 = 1 decomp_cascade_con%cascade_step_name(i_l1s1) = 'L1S1' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1s1) = rf_l1s1 cascade_donor_pool(i_l1s1) = i_met_lit cascade_receiver_pool(i_l1s1) = i_act_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1s1) = 1.0_r8 i_l2s1 = 2 decomp_cascade_con%cascade_step_name(i_l2s1) = 'L2S1' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2s1) = rf_l2s1 cascade_donor_pool(i_l2s1) = i_cel_lit cascade_receiver_pool(i_l2s1) = i_act_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2s1)= 1.0_r8 i_l3s2 = 3 decomp_cascade_con%cascade_step_name(i_l3s2) = 'L3S2' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l3s2) = rf_l3s2 cascade_donor_pool(i_l3s2) = i_lig_lit cascade_receiver_pool(i_l3s2) = i_slo_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l3s2) = 1.0_r8 i_s1s2 = 4 decomp_cascade_con%cascade_step_name(i_s1s2) = 'S1S2' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s2) = rf_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) cascade_donor_pool(i_s1s2) = i_act_som cascade_receiver_pool(i_s1s2) = i_slo_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s2) = f_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) i_s1s3 = 5 decomp_cascade_con%cascade_step_name(i_s1s3) = 'S1S3' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s3) = rf_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) cascade_donor_pool(i_s1s3) = i_act_som cascade_receiver_pool(i_s1s3) = i_pas_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s3) = f_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) i_s2s1 = 6 decomp_cascade_con%cascade_step_name(i_s2s1) = 'S2S1' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = rf_s2s1 cascade_donor_pool(i_s2s1) = i_slo_som cascade_receiver_pool(i_s2s1) = i_act_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = f_s2s1 i_s2s3 = 7 decomp_cascade_con%cascade_step_name(i_s2s3) = 'S2S3' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s3) = rf_s2s3 cascade_donor_pool(i_s2s3) = i_slo_som cascade_receiver_pool(i_s2s3) = i_pas_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s3) = f_s2s3 i_s3s1 = 8 decomp_cascade_con%cascade_step_name(i_s3s1) = 'S3S1' - rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = rf_s3s1 cascade_donor_pool(i_s3s1) = i_pas_som cascade_receiver_pool(i_s3s1) = i_act_som - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = 1.0_r8 if (.not. use_fates) then i_cwdl2 = 9 decomp_cascade_con%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_cel_lit - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel i_cwdl3 = 10 decomp_cascade_con%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_lig_lit - pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig end if - deallocate(rf_s1s2) - deallocate(rf_s1s3) - deallocate(f_s1s2) - deallocate(f_s1s3) - deallocate(params_inst%initial_Cstocks) + deallocate(params_inst%bgc_initial_Cstocks) end associate @@ -596,6 +530,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst ! ! !LOCAL VARIABLES: + real(r8), parameter :: eps = 1.e-6_r8 real(r8):: frw(bounds%begc:bounds%endc) ! rooting fraction weight real(r8), allocatable:: fr(:,:) ! column-level rooting fraction by soil depth real(r8):: psi ! temporary soilpsi for water scalar @@ -630,8 +565,10 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & catanf(t1) = 11.75_r8 +(29.7_r8 / SHR_CONST_PI) * atan( SHR_CONST_PI * 0.031_r8 * ( t1 - 15.4_r8 )) associate( & - minpsi => params_inst%minpsi_bgc , & ! Input: [real(r8) ] minimum soil suction (mm) - maxpsi => params_inst%maxpsi_bgc , & ! Input: [real(r8) ] maximum soil suction (mm) + cwd_flig => CNParamsShareInst%cwd_flig , & ! Input: [real(r8) ] lignin fraction of coarse woody debris (frac) + rf_cwdl2 => CNParamsShareInst%rf_cwdl2 , & ! Input: [real(r8) ] respiration fraction in CWD to litter2 transition (frac) + minpsi => CNParamsShareInst%minpsi , & ! Input: [real(r8) ] minimum soil suction (mm) + maxpsi => CNParamsShareInst%maxpsi , & ! Input: [real(r8) ] maximum soil suction (mm) soilpsi => soilstate_inst%soilpsi_col , & ! Input: [real(r8) (:,:) ] soil water potential in each soil layer (MPa) t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) @@ -639,7 +576,8 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & o2stress_sat => ch4_inst%o2stress_sat_col , & ! Input: [real(r8) (:,:) ] Ratio of oxygen available to that demanded by roots, aerobes, & methanotrophs (nlevsoi) o2stress_unsat => ch4_inst%o2stress_unsat_col , & ! Input: [real(r8) (:,:) ] Ratio of oxygen available to that demanded by roots, aerobes, & methanotrophs (nlevsoi) finundated => ch4_inst%finundated_col , & ! Input: [real(r8) (:) ] fractional inundated area - + rf_decomp_cascade => soilbiogeochem_carbonflux_inst%rf_decomp_cascade_col , & ! Output: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) + pathfrac_decomp_cascade => soilbiogeochem_carbonflux_inst%pathfrac_decomp_cascade_col , & ! Output: [real(r8) (:,:,:) ] what fraction of C passes from donor to receiver pool through a given transition (frac) t_scalar => soilbiogeochem_carbonflux_inst%t_scalar_col , & ! Output: [real(r8) (:,:) ] soil temperature scalar for decomp w_scalar => soilbiogeochem_carbonflux_inst%w_scalar_col , & ! Output: [real(r8) (:,:) ] soil water scalar for decomp o_scalar => soilbiogeochem_carbonflux_inst%o_scalar_col , & ! Output: [real(r8) (:,:) ] fraction by which decomposition is limited by anoxia @@ -671,7 +609,7 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & k_s1 = 1._r8 / (secspday * days_per_year * params_inst%tau_s1_bgc) k_s2 = 1._r8 / (secspday * days_per_year * params_inst%tau_s2_bgc) k_s3 = 1._r8 / (secspday * days_per_year * params_inst%tau_s3_bgc) - k_frag = 1._r8 / (secspday * days_per_year * params_inst%tau_cwd_bgc) + k_frag = 1._r8 / (secspday * days_per_year * CNParamsShareInst%tau_cwd) ! calc ref rate catanf_30 = catanf(30._r8) @@ -680,39 +618,39 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & do fc = 1,num_soilc c = filter_soilc(fc) ! - if ( abs(spinup_factor(i_met_lit) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_met_lit) - 1._r8) .gt. eps) then spinup_geogterm_l1(c) = spinup_factor(i_met_lit) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_l1(c) = 1._r8 endif ! - if ( abs(spinup_factor(i_cel_lit) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_cel_lit) - 1._r8) .gt. eps) then spinup_geogterm_l23(c) = spinup_factor(i_cel_lit) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_l23(c) = 1._r8 endif ! if ( .not. use_fates ) then - if ( abs(spinup_factor(i_cwd) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_cwd) - 1._r8) .gt. eps) then spinup_geogterm_cwd(c) = spinup_factor(i_cwd) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_cwd(c) = 1._r8 endif endif ! - if ( abs(spinup_factor(i_act_som) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_act_som) - 1._r8) .gt. eps) then spinup_geogterm_s1(c) = spinup_factor(i_act_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_s1(c) = 1._r8 endif ! - if ( abs(spinup_factor(i_slo_som) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_slo_som) - 1._r8) .gt. eps) then spinup_geogterm_s2(c) = spinup_factor(i_slo_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_s2(c) = 1._r8 endif ! - if ( abs(spinup_factor(i_pas_som) - 1._r8) .gt. .000001_r8) then + if ( abs(spinup_factor(i_pas_som) - 1._r8) .gt. eps) then spinup_geogterm_s3(c) = spinup_factor(i_pas_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) else spinup_geogterm_s3(c) = 1._r8 @@ -955,6 +893,28 @@ subroutine decomp_rate_constants_bgc(bounds, num_soilc, filter_soilc, & end if end do end do + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1s1) = 1.0_r8 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2s1) = 1.0_r8 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l3s2) = 1.0_r8 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s2) = f_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s3) = f_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = f_s2s1 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s3) = f_s2s3 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = 1.0_r8 + if (.not. use_fates) then + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = cwd_fcel + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = cwd_flig + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = rf_cwdl2 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl3) = rf_cwdl3 + end if + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1s1) = rf_l1s1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2s1) = rf_l2s1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l3s2) = rf_l3s2 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s2) = rf_s1s2(bounds%begc:bounds%endc,1:nlevdecomp) + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1s3) = rf_s1s3(bounds%begc:bounds%endc,1:nlevdecomp) + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = rf_s2s1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s3) = rf_s2s3 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = rf_s3s1 end associate diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeConType.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeConType.F90 index e62971ce59..9f8357de37 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompCascadeConType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemDecompCascadeConType.F90 @@ -30,6 +30,7 @@ module SoilBiogeochemDecompCascadeConType character(len=8) , pointer :: decomp_pool_name_history(:) ! name of pool for history files character(len=20) , pointer :: decomp_pool_name_long(:) ! name of pool for netcdf long names character(len=8) , pointer :: decomp_pool_name_short(:) ! name of pool for netcdf short names + logical , pointer :: is_microbe(:) ! TRUE => pool is a microbe pool logical , pointer :: is_litter(:) ! TRUE => pool is a litter pool logical , pointer :: is_soil(:) ! TRUE => pool is a soil pool logical , pointer :: is_cwd(:) ! TRUE => pool is a cwd pool @@ -45,6 +46,7 @@ module SoilBiogeochemDecompCascadeConType integer, public, parameter :: i_atm = 0 ! for terminal pools (i.e. 100% respiration) (only used for CN not for BGC) integer, public, parameter :: no_soil_decomp = 0 ! No soil decomposition is done integer, public, parameter :: century_decomp = 1 ! CENTURY decomposition method type + integer, public, parameter :: mimics_decomp = 2 ! MIMICS decomposition method type integer, public :: decomp_method = ispval ! Type of decomposition to use type(decomp_cascade_type), public :: decomp_cascade_con !------------------------------------------------------------------------ @@ -87,6 +89,8 @@ subroutine decomp_cascade_par_init( NLFilename ) decomp_method = no_soil_decomp case( 'CENTURYKoven2013' ) decomp_method = century_decomp + case( 'MIMICSWieder2015' ) + decomp_method = mimics_decomp case default call endrun('Bad soil_decomp_method = '//soil_decomp_method ) end select @@ -109,17 +113,17 @@ subroutine decomp_cascade_par_init( NLFilename ) if (decomp_method == century_decomp) then ndecomp_pools = 6 ndecomp_cascade_transitions = 8 - else ! TODO slevis: Currently for CN. MIMICS will get its own. + else if (decomp_method == mimics_decomp) then ndecomp_pools = 7 - ndecomp_cascade_transitions = 7 + ndecomp_cascade_transitions = 14 end if else if (decomp_method == century_decomp) then ndecomp_pools = 7 ndecomp_cascade_transitions = 10 - else ! TODO slevis: Currently for CN. MIMICS will get its own. + else if (decomp_method == mimics_decomp) then ndecomp_pools = 8 - ndecomp_cascade_transitions = 9 + ndecomp_cascade_transitions = 15 end if endif ! The next param also appears as a dimension in the params files dated @@ -166,6 +170,7 @@ subroutine init_decomp_cascade_constants( ) allocate(decomp_cascade_con%decomp_pool_name_history(ibeg:ndecomp_pools)) allocate(decomp_cascade_con%decomp_pool_name_long(ibeg:ndecomp_pools)) allocate(decomp_cascade_con%decomp_pool_name_short(ibeg:ndecomp_pools)) + allocate(decomp_cascade_con%is_microbe(ibeg:ndecomp_pools)) allocate(decomp_cascade_con%is_litter(ibeg:ndecomp_pools)) allocate(decomp_cascade_con%is_soil(ibeg:ndecomp_pools)) allocate(decomp_cascade_con%is_cwd(ibeg:ndecomp_pools)) @@ -187,6 +192,7 @@ subroutine init_decomp_cascade_constants( ) decomp_cascade_con%decomp_pool_name_restart(ibeg:ndecomp_pools) = '' decomp_cascade_con%decomp_pool_name_long(ibeg:ndecomp_pools) = '' decomp_cascade_con%decomp_pool_name_short(ibeg:ndecomp_pools) = '' + decomp_cascade_con%is_microbe(ibeg:ndecomp_pools) = .false. decomp_cascade_con%is_litter(ibeg:ndecomp_pools) = .false. decomp_cascade_con%is_soil(ibeg:ndecomp_pools) = .false. decomp_cascade_con%is_cwd(ibeg:ndecomp_pools) = .false. diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 new file mode 100644 index 0000000000..677aa4b04d --- /dev/null +++ b/src/soilbiogeochem/SoilBiogeochemDecompCascadeMIMICSMod.F90 @@ -0,0 +1,1284 @@ +module SoilBiogeochemDecompCascadeMIMICSMod + + !----------------------------------------------------------------------- + ! !DESCRIPTION: + ! Sets the coeffiecients used in the decomposition cascade submodel. + ! This uses the MIMICS parameters. + ! + ! !USES: + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_const_mod , only : SHR_CONST_TKFRZ + use shr_log_mod , only : errMsg => shr_log_errMsg + use clm_varpar , only : nlevdecomp, ndecomp_pools_max + use clm_varpar , only : i_met_lit, i_cop_mic, i_oli_mic, i_cwd + use clm_varpar , only : i_litr_min, i_litr_max, i_cwdl2 + use clm_varctl , only : iulog, spinup_state, anoxia, use_lch4, use_fates + use clm_varcon , only : zsoi + use decompMod , only : bounds_type + use spmdMod , only : masterproc + use abortutils , only : endrun + use CNSharedParamsMod , only : CNParamsShareInst, nlev_soildecomp_standard + use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con + use SoilBiogeochemStateType , only : soilbiogeochem_state_type + use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type + use SoilBiogeochemCarbonStateType , only : soilbiogeochem_carbonstate_type + use SoilStateType , only : soilstate_type + use TemperatureType , only : temperature_type + use CNVegCarbonFluxType , only : cnveg_carbonflux_type + use ch4Mod , only : ch4_type + use ColumnType , only : col + use GridcellType , only : grc + use SoilBiogeochemStateType , only : get_spinup_latitude_term + + ! + implicit none + private + ! + ! !PUBLIC MEMBER FUNCTIONS: + public :: readParams ! Read in parameters from params file + public :: init_decompcascade_mimics ! Initialization + public :: decomp_rates_mimics ! Figure out decomposition rates + ! + ! !PUBLIC DATA MEMBERS + ! + ! !PRIVATE DATA MEMBERS + ! next four 2d vars are dimensioned (columns,nlevdecomp) + real(r8), private, allocatable :: desorp(:,:) + real(r8), private, allocatable :: fphys_m1(:,:) + real(r8), private, allocatable :: fphys_m2(:,:) + real(r8), private, allocatable :: p_scalar(:,:) + integer, private :: i_phys_som ! index of physically protected Soil Organic Matter (SOM) + integer, private :: i_chem_som ! index of chemically protected SOM + integer, private :: i_avl_som ! index of available (aka active) SOM + integer, private :: i_str_lit ! index of structural litter pool + integer, private :: i_l1m1 ! indices of transitions, eg l1m1: litter 1 -> first microbial pool + integer, private :: i_l1m2 + integer, private :: i_l2m1 + integer, private :: i_l2m2 + integer, private :: i_s1m1 + integer, private :: i_s1m2 + integer, private :: i_m1s1 + integer, private :: i_m1s2 + integer, private :: i_m2s1 + integer, private :: i_m2s2 + integer, private :: i_s2s1 + integer, private :: i_s3s1 + integer, private :: i_m1s3 + integer, private :: i_m2s3 + real(r8), private :: rf_l1m1 ! respiration fractions by transition + real(r8), private :: rf_l1m2 + real(r8), private :: rf_l2m1 + real(r8), private :: rf_l2m2 + real(r8), private :: rf_s1m1 + real(r8), private :: rf_s1m2 + real(r8), private :: vint_l1_m1 ! regression intercepts by transition + real(r8), private :: vint_l2_m1 + real(r8), private :: vint_s1_m1 + real(r8), private :: vint_l1_m2 + real(r8), private :: vint_l2_m2 + real(r8), private :: vint_s1_m2 + real(r8), private :: kint_l1_m1 ! regression intercepts by transition + real(r8), private :: kint_l2_m1 + real(r8), private :: kint_s1_m1 + real(r8), private :: kint_l1_m2 + real(r8), private :: kint_l2_m2 + real(r8), private :: kint_s1_m2 + real(r8), private :: vmod_l1_m1 ! vmod = vmod * av from Wieder et al 2015 + real(r8), private :: vmod_l2_m1 + real(r8), private :: vmod_s1_m1 + real(r8), private :: vmod_l1_m2 + real(r8), private :: vmod_l2_m2 + real(r8), private :: vmod_s1_m2 + real(r8), private :: kmod_l1_m1 ! kmod = ak / kmod from Wieder et al 2015 + real(r8), private :: kmod_l2_m1 + real(r8), private :: kmod_s1_m1 + real(r8), private :: kmod_l1_m2 + real(r8), private :: kmod_l2_m2 + real(r8), private :: kmod_s1_m2 + real(r8), private :: vslope_l1_m1 ! regression coefficients by transition + real(r8), private :: vslope_l2_m1 + real(r8), private :: vslope_s1_m1 + real(r8), private :: vslope_l1_m2 + real(r8), private :: vslope_l2_m2 + real(r8), private :: vslope_s1_m2 + real(r8), private :: kslope_l1_m1 ! regression coefficients by transition + real(r8), private :: kslope_l2_m1 + real(r8), private :: kslope_s1_m1 + real(r8), private :: kslope_l1_m2 + real(r8), private :: kslope_l2_m2 + real(r8), private :: kslope_s1_m2 + + type, private :: params_type + real(r8) :: mimics_nue_into_mic ! microbial N use efficiency for N fluxes + real(r8) :: mimics_desorpQ10 + real(r8) :: mimics_densdep ! exponent controling the density dependence of microbial turnover + real(r8) :: mimics_tau_mod_factor ! (1 / tauModDenom) from testbed code + real(r8) :: mimics_tau_mod_min + real(r8) :: mimics_tau_mod_max + real(r8) :: mimics_ko_r + real(r8) :: mimics_ko_k + real(r8) :: mimics_cn_r ! C:N of MICr + real(r8) :: mimics_cn_k ! C:N of MICk + real(r8) :: mimics_cn_mod_num ! adjusts microbial CN based on fmet + real(r8) :: mimics_t_soi_ref ! reference soil temperature (degC) + real(r8) :: mimics_initial_Cstocks_depth ! Soil depth for initial C stocks for a cold-start (m) + real(r8), allocatable :: mimics_initial_Cstocks(:) ! Initial C stocks for a cold-start (gC/m3) + ! The next few vectors are dimensioned by the number of decomposition + ! transitions that make use of the corresponding parameters, currently + ! six. The transitions are represented in this order: + ! l1m1 l2m1 s1m1 l1m2 l2m2 s1m2 + real(r8), allocatable :: mimics_mge(:) ! Microbial growth efficiency (mg/mg) + real(r8), allocatable :: mimics_vmod(:) ! vmod = vmod * av from Wieder et al 2015 + real(r8), allocatable :: mimics_vint(:) ! regression intercepts (5.47 ln(mg Cs (mg MIC)-1 h-1) ) + real(r8), allocatable :: mimics_vslope(:) ! regression coeffs (ln(mg Cs (mg MIC)-1 h-1) ¡C-1) + real(r8), allocatable :: mimics_kmod(:) ! kmod = ak / kmod from Wieder et al 2015 + real(r8), allocatable :: mimics_kint(:) ! regression intercepts + real(r8), allocatable :: mimics_kslope(:) ! regression coeffs + ! The next few vectors are dimensioned by the number of parameters with + ! the same name (eg 2 for mimics_tau_r_p1, mimics_tau_r_p2) used in the + ! respective formula. In the formulas, we use scalar copies of each + ! parameter with suffixes _p1, p2, p3, ... to distinguish among them. + ! See allocate statements below for the size of each of the following + ! vectors. + real(r8), allocatable :: mimics_fmet(:) + real(r8), allocatable :: mimics_p_scalar(:) + real(r8), allocatable :: mimics_fphys_r(:) + real(r8), allocatable :: mimics_fphys_k(:) + real(r8), allocatable :: mimics_fchem_r(:) + real(r8), allocatable :: mimics_fchem_k(:) + real(r8), allocatable :: mimics_desorp(:) + real(r8), allocatable :: mimics_tau_r(:) + real(r8), allocatable :: mimics_tau_k(:) + end type params_type + ! + type(params_type), private :: params_inst + + character(len=*), parameter, private :: sourcefile = & + __FILE__ + + !----------------------------------------------------------------------- + +contains + + !----------------------------------------------------------------------- + subroutine readParams ( ncid ) + ! + ! !DESCRIPTION: + ! + ! !USES: + use ncdio_pio , only: file_desc_t,ncd_io + ! + ! !ARGUMENTS: + type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id + ! + ! !LOCAL VARIABLES: + character(len=32) :: subname = 'readMimicsParams' + character(len=100) :: errCode = 'Error reading MIMICS params ' + logical :: readv ! has variable been read in or not + real(r8) :: tempr ! temporary to read in constant + character(len=100) :: tString ! temp. var for reading + !----------------------------------------------------------------------- + + ! Read off of netcdf file + tString='mimics_initial_Cstocks_depth' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_initial_Cstocks_depth=tempr + + allocate(params_inst%mimics_initial_Cstocks(ndecomp_pools_max)) + tString='mimics_initial_Cstocks' + call ncd_io(trim(tString), params_inst%mimics_initial_Cstocks(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_mge(ndecomp_pools_max)) + tString='mimics_mge' + call ncd_io(trim(tString), params_inst%mimics_mge(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_vmod(ndecomp_pools_max)) + tString='mimics_vmod' + call ncd_io(trim(tString), params_inst%mimics_vmod(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_vslope(ndecomp_pools_max)) + tString='mimics_vslope' + call ncd_io(trim(tString), params_inst%mimics_vslope(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_vint(ndecomp_pools_max)) + tString='mimics_vint' + call ncd_io(trim(tString), params_inst%mimics_vint(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_kmod(ndecomp_pools_max)) + tString='mimics_kmod' + call ncd_io(trim(tString), params_inst%mimics_kmod(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_kslope(ndecomp_pools_max)) + tString='mimics_kslope' + call ncd_io(trim(tString), params_inst%mimics_kslope(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_kint(ndecomp_pools_max)) + tString='mimics_kint' + call ncd_io(trim(tString), params_inst%mimics_kint(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_p_scalar(2)) + tString='mimics_p_scalar' + call ncd_io(trim(tString), params_inst%mimics_p_scalar(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_desorp(2)) + tString='mimics_desorp' + call ncd_io(trim(tString), params_inst%mimics_desorp(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_fphys_r(2)) + tString='mimics_fphys_r' + call ncd_io(trim(tString), params_inst%mimics_fphys_r(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_fphys_k(2)) + tString='mimics_fphys_k' + call ncd_io(trim(tString), params_inst%mimics_fphys_k(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_fmet(4)) + tString='mimics_fmet' + call ncd_io(trim(tString), params_inst%mimics_fmet(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_fchem_r(4)) + tString='mimics_fchem_r' + call ncd_io(trim(tString), params_inst%mimics_fchem_r(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_fchem_k(4)) + tString='mimics_fchem_k' + call ncd_io(trim(tString), params_inst%mimics_fchem_k(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_tau_r(2)) + tString='mimics_tau_r' + call ncd_io(trim(tString), params_inst%mimics_tau_r(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + allocate(params_inst%mimics_tau_k(2)) + tString='mimics_tau_k' + call ncd_io(trim(tString), params_inst%mimics_tau_k(:), 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + + tString='mimics_nue_into_mic' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_nue_into_mic = tempr + + tString='mimics_tau_mod_factor' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_tau_mod_factor = tempr + + tString='mimics_tau_mod_min' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_tau_mod_min = tempr + + tString='mimics_tau_mod_max' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_tau_mod_max = tempr + + tString='mimics_ko_r' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_ko_r = tempr + + tString='mimics_ko_k' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_ko_k = tempr + + tString='mimics_densdep' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_densdep = tempr + + tString='mimics_desorpQ10' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_desorpQ10 = tempr + + tString='mimics_t_soi_ref' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_t_soi_ref = tempr + + tString='mimics_cn_mod_num' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_cn_mod_num = tempr + + tString='mimics_cn_r' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_cn_r = tempr + + tString='mimics_cn_k' + call ncd_io(trim(tString), tempr, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=trim(errCode)//trim(tString)//errMsg(sourcefile, __LINE__)) + params_inst%mimics_cn_k = tempr + + end subroutine readParams + + !----------------------------------------------------------------------- + subroutine init_decompcascade_mimics(bounds, soilbiogeochem_state_inst, soilstate_inst ) + ! + ! !DESCRIPTION: + ! initialize rate constants and decomposition pathways following the + ! decomposition cascade of the MIMICS model. + ! + ! !USES: + use clm_varcon, only: pct_to_frac + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + type(soilbiogeochem_state_type) , intent(inout) :: soilbiogeochem_state_inst + type(soilstate_type) , intent(in) :: soilstate_inst + ! + ! !LOCAL VARIABLES + !-- properties of each decomposing pool + real(r8) :: mimics_nue_into_mic + real(r8) :: mimics_p_scalar_p1 + real(r8) :: mimics_p_scalar_p2 + real(r8) :: mimics_fphys_r_p1 + real(r8) :: mimics_fphys_r_p2 + real(r8) :: mimics_fphys_k_p1 + real(r8) :: mimics_fphys_k_p2 + real(r8) :: mimics_desorp_p1 + real(r8) :: mimics_desorp_p2 + + real(r8):: speedup_fac ! acceleration factor, higher when vertsoilc = .true. + + real(r8) :: clay_frac ! local copy of cellclay converted from % (fraction) + integer :: c, j ! indices + !----------------------------------------------------------------------- + + associate( & + nue_decomp_cascade => soilbiogeochem_state_inst%nue_decomp_cascade_col , & ! Output: [real(r8) (:) ] N use efficiency for a given transition (gN going into microbe / gN decomposed) + + cellclay => soilstate_inst%cellclay_col , & ! Input: [real(r8) (:,:) ] column 3D clay (%) + + cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Output: [integer (:) ] which pool is C taken from for a given decomposition step + cascade_receiver_pool => decomp_cascade_con%cascade_receiver_pool , & ! Output: [integer (:) ] which pool is C added to for a given decomposition step + floating_cn_ratio_decomp_pools => decomp_cascade_con%floating_cn_ratio_decomp_pools , & ! Output: [logical (:) ] TRUE => pool has fixed C:N ratio + is_microbe => decomp_cascade_con%is_microbe , & ! Output: [logical (:) ] TRUE => pool is microbial + is_litter => decomp_cascade_con%is_litter , & ! Output: [logical (:) ] TRUE => pool is a litter pool + is_soil => decomp_cascade_con%is_soil , & ! Output: [logical (:) ] TRUE => pool is a soil pool + is_cwd => decomp_cascade_con%is_cwd , & ! Output: [logical (:) ] TRUE => pool is a cwd pool + initial_cn_ratio => decomp_cascade_con%initial_cn_ratio , & ! Output: [real(r8) (:) ] c:n ratio for initialization of pools + initial_stock => decomp_cascade_con%initial_stock , & ! Output: [real(r8) (:) ] initial concentration for seeding at spinup + initial_stock_soildepth => decomp_cascade_con%initial_stock_soildepth , & ! Output: [real(r8) (:) ] soil depth for initial concentration for seeding at spinup + is_metabolic => decomp_cascade_con%is_metabolic , & ! Output: [logical (:) ] TRUE => pool is metabolic material + is_cellulose => decomp_cascade_con%is_cellulose , & ! Output: [logical (:) ] TRUE => pool is cellulose + is_lignin => decomp_cascade_con%is_lignin , & ! Output: [logical (:) ] TRUE => pool is lignin + spinup_factor => decomp_cascade_con%spinup_factor & ! Output: [real(r8) (:) ] factor for AD spinup associated with each pool + + ) + + allocate(desorp(bounds%begc:bounds%endc,1:nlevdecomp)) + allocate(fphys_m1(bounds%begc:bounds%endc,1:nlevdecomp)) + allocate(fphys_m2(bounds%begc:bounds%endc,1:nlevdecomp)) + allocate(p_scalar(bounds%begc:bounds%endc,1:nlevdecomp)) + + !------- time-constant coefficients ---------- ! + mimics_nue_into_mic = params_inst%mimics_nue_into_mic + mimics_p_scalar_p1 = params_inst%mimics_p_scalar(1) + mimics_p_scalar_p2 = params_inst%mimics_p_scalar(2) + mimics_fphys_r_p1 = params_inst%mimics_fphys_r(1) + mimics_fphys_r_p2 = params_inst%mimics_fphys_r(2) + mimics_fphys_k_p1 = params_inst%mimics_fphys_k(1) + mimics_fphys_k_p2 = params_inst%mimics_fphys_k(2) + mimics_desorp_p1 = params_inst%mimics_desorp(1) + mimics_desorp_p2 = params_inst%mimics_desorp(2) + + ! set respiration fractions for fluxes between compartments + rf_l1m1 = 1.0_r8 - params_inst%mimics_mge(1) + rf_l2m1 = 1.0_r8 - params_inst%mimics_mge(2) + rf_s1m1 = 1.0_r8 - params_inst%mimics_mge(3) + rf_l1m2 = 1.0_r8 - params_inst%mimics_mge(4) + rf_l2m2 = 1.0_r8 - params_inst%mimics_mge(5) + rf_s1m2 = 1.0_r8 - params_inst%mimics_mge(6) + + ! vmod = "old" vmod * av AND kmod = ak / "old" kmod + ! Table B1 Wieder et al. (2015) and MIMICS params file give diff + ! ak and av values. I used the params file values. + vmod_l1_m1 = params_inst%mimics_vmod(1) + vmod_l2_m1 = params_inst%mimics_vmod(2) + vmod_s1_m1 = params_inst%mimics_vmod(3) + vmod_l1_m2 = params_inst%mimics_vmod(4) + vmod_l2_m2 = params_inst%mimics_vmod(5) + vmod_s1_m2 = params_inst%mimics_vmod(6) + kmod_l1_m1 = params_inst%mimics_kmod(1) + kmod_l2_m1 = params_inst%mimics_kmod(2) + kmod_s1_m1 = params_inst%mimics_kmod(3) + kmod_l1_m2 = params_inst%mimics_kmod(4) + kmod_l2_m2 = params_inst%mimics_kmod(5) + kmod_s1_m2 = params_inst%mimics_kmod(6) + vslope_l1_m1 = params_inst%mimics_vslope(1) + vslope_l2_m1 = params_inst%mimics_vslope(2) + vslope_s1_m1 = params_inst%mimics_vslope(3) + vslope_l1_m2 = params_inst%mimics_vslope(4) + vslope_l2_m2 = params_inst%mimics_vslope(5) + vslope_s1_m2 = params_inst%mimics_vslope(6) + kslope_l1_m1 = params_inst%mimics_kslope(1) + kslope_l2_m1 = params_inst%mimics_kslope(2) + kslope_s1_m1 = params_inst%mimics_kslope(3) + kslope_l1_m2 = params_inst%mimics_kslope(4) + kslope_l2_m2 = params_inst%mimics_kslope(5) + kslope_s1_m2 = params_inst%mimics_kslope(6) + vint_l1_m1 = params_inst%mimics_vint(1) + vint_l2_m1 = params_inst%mimics_vint(2) + vint_s1_m1 = params_inst%mimics_vint(3) + vint_l1_m2 = params_inst%mimics_vint(4) + vint_l2_m2 = params_inst%mimics_vint(5) + vint_s1_m2 = params_inst%mimics_vint(6) + kint_l1_m1 = params_inst%mimics_kint(1) + kint_l2_m1 = params_inst%mimics_kint(2) + kint_s1_m1 = params_inst%mimics_kint(3) + kint_l1_m2 = params_inst%mimics_kint(4) + kint_l2_m2 = params_inst%mimics_kint(5) + kint_s1_m2 = params_inst%mimics_kint(6) + + ! some of these are dependent on the soil texture properties + ! One-time initializations here. + ! Time-dep params in subr. decomp_rates_mimics. + + do c = bounds%begc, bounds%endc + do j = 1, nlevdecomp + ! The parameter values currently in the params files always lead to + ! positive values for the expressions below, so we do not + ! need to use the max function to limit these expressions. + ! We apply the min function on cellclay because we are looping over + ! some non-soil columns here that contain cellclay = 1e36. + clay_frac = pct_to_frac * & + dmin1(100.0_r8, cellclay(c,j)) ! conv. % to fraction + desorp(c,j) = mimics_desorp_p1 * dexp(mimics_desorp_p2 * clay_frac) + fphys_m1(c,j) = dmin1(1.0_r8, mimics_fphys_r_p1 * & + dexp(mimics_fphys_r_p2 * clay_frac)) + fphys_m2(c,j) = dmin1(1.0_r8, mimics_fphys_k_p1 * & + dexp(mimics_fphys_k_p2 * clay_frac)) + p_scalar(c,j) = 1.0_r8 / (mimics_p_scalar_p1 * & + dexp(mimics_p_scalar_p2 * dsqrt(clay_frac))) + end do + end do + initial_stock_soildepth = params_inst%mimics_initial_Cstocks_depth + + !------------------- list of pools and their attributes ------------ + i_litr_min = 1 + i_met_lit = i_litr_min + floating_cn_ratio_decomp_pools(i_met_lit) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_met_lit) = 'litr1' + decomp_cascade_con%decomp_pool_name_history(i_met_lit) = 'MET_LIT' + decomp_cascade_con%decomp_pool_name_long(i_met_lit) = 'metabolic litter' + decomp_cascade_con%decomp_pool_name_short(i_met_lit) = 'L1' + is_microbe(i_met_lit) = .false. + is_litter(i_met_lit) = .true. + is_soil(i_met_lit) = .false. + is_cwd(i_met_lit) = .false. + initial_cn_ratio(i_met_lit) = 10._r8 ! 90 in BGC; not used in MIMICS + initial_stock(i_met_lit) = params_inst%mimics_initial_Cstocks(i_met_lit) + is_metabolic(i_met_lit) = .true. + is_cellulose(i_met_lit) = .false. + is_lignin(i_met_lit) = .false. + + i_str_lit = i_met_lit + 1 + floating_cn_ratio_decomp_pools(i_str_lit) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_str_lit) = 'litr2' + decomp_cascade_con%decomp_pool_name_history(i_str_lit) = 'STR_LIT' + decomp_cascade_con%decomp_pool_name_long(i_str_lit) = 'structural litter' + decomp_cascade_con%decomp_pool_name_short(i_str_lit) = 'L2' + is_microbe(i_str_lit) = .false. + is_litter(i_str_lit) = .true. + is_soil(i_str_lit) = .false. + is_cwd(i_str_lit) = .false. + initial_cn_ratio(i_str_lit) = 10._r8 ! 90 in BGC; not used in MIMICS + initial_stock(i_str_lit) = params_inst%mimics_initial_Cstocks(i_str_lit) + is_metabolic(i_str_lit) = .false. + is_cellulose(i_str_lit) = .true. + is_lignin(i_str_lit) = .true. + + i_litr_max = i_str_lit + if (i_litr_min /= 1 .or. i_litr_max < 2 .or. i_litr_max > 3) then + write(iulog,*) 'Expecting i_litr_min = 1 and i_litr_max = 2 or 3.' + write(iulog,*) 'See pftconMod, SoilBiogeochemCarbonFluxType, and' + write(iulog,*) 'clmfates_interfaceMod for ramifications of changing' + write(iulog,*) 'this assumption.' + call endrun(msg='ERROR: i_litr_min and/or i_litr_max out of range '// & + errMsg(sourcefile, __LINE__)) + end if + + i_avl_som = i_str_lit + 1 + floating_cn_ratio_decomp_pools(i_avl_som) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_avl_som) = 'soil1' + decomp_cascade_con%decomp_pool_name_history(i_avl_som) = 'AVL_SOM' + decomp_cascade_con%decomp_pool_name_long(i_avl_som) = 'available soil organic matter' + decomp_cascade_con%decomp_pool_name_short(i_avl_som) = 'S1' + is_microbe(i_avl_som) = .false. + is_litter(i_avl_som) = .false. + is_soil(i_avl_som) = .true. + is_cwd(i_avl_som) = .false. + initial_cn_ratio(i_avl_som) = 10._r8 ! cn_s1 in BGC; not used in MIMICS + initial_stock(i_avl_som) = params_inst%mimics_initial_Cstocks(i_avl_som) + is_metabolic(i_avl_som) = .false. + is_cellulose(i_avl_som) = .false. + is_lignin(i_avl_som) = .false. + + i_chem_som = i_avl_som + 1 + floating_cn_ratio_decomp_pools(i_chem_som) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_chem_som) = 'soil2' + decomp_cascade_con%decomp_pool_name_history(i_chem_som) = 'CHEM_SOM' + decomp_cascade_con%decomp_pool_name_long(i_chem_som) = 'chemically protected soil organic matter' + decomp_cascade_con%decomp_pool_name_short(i_chem_som) = 'S2' + is_microbe(i_chem_som) = .false. + is_litter(i_chem_som) = .false. + is_soil(i_chem_som) = .true. + is_cwd(i_chem_som) = .false. + initial_cn_ratio(i_chem_som) = 10._r8 ! cn_s2 in BGC; not used in MIMICS + initial_stock(i_chem_som) = params_inst%mimics_initial_Cstocks(i_chem_som) + is_metabolic(i_chem_som) = .false. + is_cellulose(i_chem_som) = .false. + is_lignin(i_chem_som) = .false. + + i_phys_som = i_chem_som + 1 + floating_cn_ratio_decomp_pools(i_phys_som) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_phys_som) = 'soil3' + decomp_cascade_con%decomp_pool_name_history(i_phys_som) = 'PHYS_SOM' + decomp_cascade_con%decomp_pool_name_long(i_phys_som) = 'physically protected soil organic matter' + decomp_cascade_con%decomp_pool_name_short(i_phys_som) = 'S3' + is_microbe(i_phys_som) = .false. + is_litter(i_phys_som) = .false. + is_soil(i_phys_som) = .true. + is_cwd(i_phys_som) = .false. + initial_cn_ratio(i_phys_som) = 10._r8 ! cn_s3 in BGC; not used in MIMICS + initial_stock(i_phys_som) = params_inst%mimics_initial_Cstocks(i_phys_som) + is_metabolic(i_phys_som) = .false. + is_cellulose(i_phys_som) = .false. + is_lignin(i_phys_som) = .false. + + i_cop_mic = i_phys_som + 1 + floating_cn_ratio_decomp_pools(i_cop_mic) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_cop_mic) = 'micr1' + decomp_cascade_con%decomp_pool_name_history(i_cop_mic) = 'COP_MIC' + decomp_cascade_con%decomp_pool_name_long(i_cop_mic) = 'copiotrophic microbes' + decomp_cascade_con%decomp_pool_name_short(i_cop_mic) = 'M1' + is_microbe(i_cop_mic) = .true. + is_litter(i_cop_mic) = .false. + is_soil(i_cop_mic) = .false. + is_cwd(i_cop_mic) = .false. + initial_cn_ratio(i_cop_mic) = 10._r8 ! MIMICS may use this + initial_stock(i_cop_mic) = params_inst%mimics_initial_Cstocks(i_cop_mic) + is_metabolic(i_cop_mic) = .false. + is_cellulose(i_cop_mic) = .false. + is_lignin(i_cop_mic) = .false. + + i_oli_mic = i_cop_mic + 1 + floating_cn_ratio_decomp_pools(i_oli_mic) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_oli_mic) = 'micr2' + decomp_cascade_con%decomp_pool_name_history(i_oli_mic) = 'OLI_MIC' + decomp_cascade_con%decomp_pool_name_long(i_oli_mic) = 'oligotrophic microbes' + decomp_cascade_con%decomp_pool_name_short(i_oli_mic) = 'M2' + is_microbe(i_oli_mic) = .true. + is_litter(i_oli_mic) = .false. + is_soil(i_oli_mic) = .false. + is_cwd(i_oli_mic) = .false. + initial_cn_ratio(i_oli_mic) = 10._r8 ! MIMICS may use this + initial_stock(i_oli_mic) = params_inst%mimics_initial_Cstocks(i_oli_mic) + is_metabolic(i_oli_mic) = .false. + is_cellulose(i_oli_mic) = .false. + is_lignin(i_oli_mic) = .false. + + if (.not. use_fates) then + ! CWD + i_cwd = i_oli_mic + 1 + floating_cn_ratio_decomp_pools(i_cwd) = .true. + decomp_cascade_con%decomp_pool_name_restart(i_cwd) = 'cwd' + decomp_cascade_con%decomp_pool_name_history(i_cwd) = 'CWD' + decomp_cascade_con%decomp_pool_name_long(i_cwd) = 'coarse woody debris' + decomp_cascade_con%decomp_pool_name_short(i_cwd) = 'CWD' + is_microbe(i_cwd) = .false. + is_litter(i_cwd) = .false. + is_soil(i_cwd) = .false. + is_cwd(i_cwd) = .true. + initial_cn_ratio(i_cwd) = 10._r8 ! 90 in BGC; not used in MIMICS + initial_stock(i_cwd) = params_inst%mimics_initial_Cstocks(i_cwd) + is_metabolic(i_cwd) = .false. + is_cellulose(i_cwd) = .false. + is_lignin(i_cwd) = .false. + endif + + speedup_fac = 1._r8 + + !lit1,2 + spinup_factor(i_met_lit) = 1._r8 + spinup_factor(i_str_lit) = 1._r8 + !CWD + if (.not. use_fates) then + spinup_factor(i_cwd) = max(1._r8, (speedup_fac * CNParamsShareInst%tau_cwd * 0.5_r8 )) + end if + !som1,2,3 + spinup_factor(i_avl_som) = 1._r8 + spinup_factor(i_chem_som) = 1._r8 ! BGC used cwd formula above but + spinup_factor(i_phys_som) = 1._r8 ! ...w the respective tau_s values + ! micr1,2 + spinup_factor(i_cop_mic) = 1._r8 + spinup_factor(i_oli_mic) = 1._r8 + + if ( masterproc ) then + write(iulog,*) 'Spinup_state ',spinup_state + write(iulog,*) 'Spinup factors ',spinup_factor + end if + + !---------------- list of transitions and their time-independent coefficients ---------------! + i_l1m1 = 1 + decomp_cascade_con%cascade_step_name(i_l1m1) = 'L1M1' + cascade_donor_pool(i_l1m1) = i_met_lit + cascade_receiver_pool(i_l1m1) = i_cop_mic + nue_decomp_cascade(i_l1m1) = mimics_nue_into_mic + + i_l1m2 = 2 + decomp_cascade_con%cascade_step_name(i_l1m2) = 'L1M2' + cascade_donor_pool(i_l1m2) = i_met_lit + cascade_receiver_pool(i_l1m2) = i_oli_mic + nue_decomp_cascade(i_l1m2) = mimics_nue_into_mic + + i_l2m1 = 3 + decomp_cascade_con%cascade_step_name(i_l2m1) = 'L2M1' + cascade_donor_pool(i_l2m1) = i_str_lit + cascade_receiver_pool(i_l2m1) = i_cop_mic + nue_decomp_cascade(i_l2m1) = mimics_nue_into_mic + + i_l2m2 = 4 + decomp_cascade_con%cascade_step_name(i_l2m2) = 'L2M2' + cascade_donor_pool(i_l2m2) = i_str_lit + cascade_receiver_pool(i_l2m2) = i_oli_mic + nue_decomp_cascade(i_l2m2) = mimics_nue_into_mic + + i_s1m1 = 5 + decomp_cascade_con%cascade_step_name(i_s1m1) = 'S1M1' + cascade_donor_pool(i_s1m1) = i_avl_som + cascade_receiver_pool(i_s1m1) = i_cop_mic + nue_decomp_cascade(i_s1m1) = mimics_nue_into_mic + + i_s1m2 = 6 + decomp_cascade_con%cascade_step_name(i_s1m2) = 'S1M2' + cascade_donor_pool(i_s1m2) = i_avl_som + cascade_receiver_pool(i_s1m2) = i_oli_mic + nue_decomp_cascade(i_s1m2) = mimics_nue_into_mic + + i_s2s1 = 7 + decomp_cascade_con%cascade_step_name(i_s2s1) = 'S2S1' + cascade_donor_pool(i_s2s1) = i_chem_som + cascade_receiver_pool(i_s2s1) = i_avl_som + nue_decomp_cascade(i_s2s1) = 1.0_r8 + + i_s3s1 = 8 + decomp_cascade_con%cascade_step_name(i_s3s1) = 'S3S1' + cascade_donor_pool(i_s3s1) = i_phys_som + cascade_receiver_pool(i_s3s1) = i_avl_som + nue_decomp_cascade(i_s3s1) = 1.0_r8 + + i_m1s1 = 9 + decomp_cascade_con%cascade_step_name(i_m1s1) = 'M1S1' + cascade_donor_pool(i_m1s1) = i_cop_mic + cascade_receiver_pool(i_m1s1) = i_avl_som + nue_decomp_cascade(i_m1s1) = 1.0_r8 + + i_m1s2 = 10 + decomp_cascade_con%cascade_step_name(i_m1s2) = 'M1S2' + cascade_donor_pool(i_m1s2) = i_cop_mic + cascade_receiver_pool(i_m1s2) = i_chem_som + nue_decomp_cascade(i_m1s2) = 1.0_r8 + + i_m1s3 = 11 + decomp_cascade_con%cascade_step_name(i_m1s3) = 'M1S3' + cascade_donor_pool(i_m1s3) = i_cop_mic + cascade_receiver_pool(i_m1s3) = i_phys_som + nue_decomp_cascade(i_m1s3) = 1.0_r8 + + i_m2s1 = 12 + decomp_cascade_con%cascade_step_name(i_m2s1) = 'M2S1' + cascade_donor_pool(i_m2s1) = i_oli_mic + cascade_receiver_pool(i_m2s1) = i_avl_som + nue_decomp_cascade(i_m2s1) = 1.0_r8 + + i_m2s2 = 13 + decomp_cascade_con%cascade_step_name(i_m2s2) = 'M2S2' + cascade_donor_pool(i_m2s2) = i_oli_mic + cascade_receiver_pool(i_m2s2) = i_chem_som + nue_decomp_cascade(i_m2s2) = 1.0_r8 + + i_m2s3 = 14 + decomp_cascade_con%cascade_step_name(i_m2s3) = 'M2S3' + cascade_donor_pool(i_m2s3) = i_oli_mic + cascade_receiver_pool(i_m2s3) = i_phys_som + nue_decomp_cascade(i_m2s3) = 1.0_r8 + + if (.not. use_fates) then + i_cwdl2 = 15 + decomp_cascade_con%cascade_step_name(i_cwdl2) = 'CWDL2' + cascade_donor_pool(i_cwdl2) = i_cwd + cascade_receiver_pool(i_cwdl2) = i_str_lit + nue_decomp_cascade(i_cwdl2) = 1.0_r8 + end if + + deallocate(params_inst%mimics_mge) + deallocate(params_inst%mimics_vmod) + deallocate(params_inst%mimics_vint) + deallocate(params_inst%mimics_vslope) + deallocate(params_inst%mimics_kmod) + deallocate(params_inst%mimics_kint) + deallocate(params_inst%mimics_kslope) + deallocate(params_inst%mimics_p_scalar) + deallocate(params_inst%mimics_desorp) + deallocate(params_inst%mimics_fphys_r) + deallocate(params_inst%mimics_fphys_k) + deallocate(params_inst%mimics_initial_Cstocks) + + end associate + + end subroutine init_decompcascade_mimics + + !----------------------------------------------------------------------- + subroutine decomp_rates_mimics(bounds, num_soilc, filter_soilc, & + soilstate_inst, temperature_inst, cnveg_carbonflux_inst, & + ch4_inst, soilbiogeochem_carbonflux_inst, soilbiogeochem_carbonstate_inst) + ! + ! !DESCRIPTION: + ! Calculate rates and decomposition pathways for the MIMICS + ! decomposition cascade model + ! + ! !USES: + use clm_time_manager , only : get_curr_days_per_year + use clm_varcon , only : secspday, secsphr, tfrz + use clm_varcon , only : g_to_mg, cm3_to_m3 + ! + ! !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 + type(soilstate_type) , intent(in) :: soilstate_inst + type(temperature_type) , intent(in) :: temperature_inst + type(cnveg_carbonflux_type) , intent(in) :: cnveg_carbonflux_inst + type(ch4_type) , intent(in) :: ch4_inst + type(soilbiogeochem_carbonflux_type) , intent(inout) :: soilbiogeochem_carbonflux_inst + type(soilbiogeochem_carbonstate_type), intent(in) :: soilbiogeochem_carbonstate_inst + ! + ! !LOCAL VARIABLES: + real(r8), parameter :: eps = 1.e-6_r8 + real(r8):: frw(bounds%begc:bounds%endc) ! rooting fraction weight + real(r8), allocatable:: fr(:,:) ! column-level rooting fraction by soil depth + real(r8):: psi ! temporary soilpsi for water scalar + real(r8):: k_frag ! fragmentation rate constant CWD (1/sec) + real(r8):: fmet + real(r8):: favl + real(r8):: fchem_m1 + real(r8):: fchem_m2 + real(r8):: desorption + real(r8):: vmax_l1_m1 ! + real(r8):: vmax_l2_m1 ! + real(r8):: vmax_s1_m1 ! + real(r8):: vmax_l1_m2 ! + real(r8):: vmax_l2_m2 ! + real(r8):: vmax_s1_m2 ! + real(r8):: km_l1_m1 ! + real(r8):: km_l2_m1 ! + real(r8):: km_s1_m1 ! + real(r8):: km_l1_m2 ! + real(r8):: km_l2_m2 ! + real(r8):: km_s1_m2 ! + real(r8):: tau_m1 ! + real(r8):: tau_m2 ! + real(r8):: tau_mod + real(r8):: m1_conc ! + real(r8):: m2_conc ! + real(r8):: term_1 ! + real(r8):: term_2 ! + real(r8):: t_soi_degC +! real(r8):: decomp_depth_efolding ! (meters) e-folding depth for reduction in decomposition [ + integer :: c, fc, j, k, l + real(r8):: days_per_year ! days per year + real(r8):: depth_scalar(bounds%begc:bounds%endc,1:nlevdecomp) + real(r8):: w_d_o_scalars ! product of w_scalar * depth_scalar * o_scalar + real(r8):: mino2lim !minimum anaerobic decomposition rate + real(r8):: mimics_fmet_p1 + real(r8):: mimics_fmet_p2 + real(r8):: mimics_fmet_p3 + real(r8):: mimics_fmet_p4 + real(r8):: mimics_fchem_r_p1 + real(r8):: mimics_fchem_r_p2 + real(r8):: mimics_fchem_r_p3 + real(r8):: mimics_fchem_k_p1 + real(r8):: mimics_fchem_k_p2 + real(r8):: mimics_fchem_k_p3 + real(r8):: mimics_tau_mod_min + real(r8):: mimics_tau_mod_max + real(r8):: mimics_tau_mod_factor + real(r8):: mimics_tau_r_p1 + real(r8):: mimics_tau_r_p2 + real(r8):: mimics_tau_k_p1 + real(r8):: mimics_tau_k_p2 + real(r8):: mimics_ko_r + real(r8):: mimics_ko_k + real(r8):: mimics_densdep + real(r8):: mimics_desorpQ10 + real(r8):: mimics_t_soi_ref + real(r8):: mimics_cn_mod_num + real(r8):: mimics_cn_r + real(r8):: mimics_cn_k + + real(r8):: spinup_geogterm_l1(bounds%begc:bounds%endc) ! geographically-varying spinup term for l1 + real(r8):: spinup_geogterm_l2(bounds%begc:bounds%endc) ! geographically-varying spinup term for l2 + real(r8):: spinup_geogterm_cwd(bounds%begc:bounds%endc) ! geographically-varying spinup term for cwd + real(r8):: spinup_geogterm_s1(bounds%begc:bounds%endc) ! geographically-varying spinup term for s1 + real(r8):: spinup_geogterm_s2(bounds%begc:bounds%endc) ! geographically-varying spinup term for s2 + real(r8):: spinup_geogterm_s3(bounds%begc:bounds%endc) ! geographically-varying spinup term for s3 + real(r8):: spinup_geogterm_m1(bounds%begc:bounds%endc) ! geographically-varying spinup term for m1 + real(r8):: spinup_geogterm_m2(bounds%begc:bounds%endc) ! geographically-varying spinup term for m2 + + !----------------------------------------------------------------------- + + associate( & + ligninNratioAvg => cnveg_carbonflux_inst%ligninNratioAvg_col , & ! Input: [real(r8) (:) ] column-level lignin to nitrogen ratio + annsum_npp_col => cnveg_carbonflux_inst%annsum_npp_col , & ! Input: [real(r8) (:) ] annual sum of NPP at the column level (gC/m2/yr) + + rf_cwdl2 => CNParamsShareInst%rf_cwdl2 , & ! Input: [real(r8) ] respiration fraction in CWD to litter2 transition (frac) + minpsi => CNParamsShareInst%minpsi , & ! Input: [real(r8) ] minimum soil suction (mm) + maxpsi => CNParamsShareInst%maxpsi , & ! Input: [real(r8) ] maximum soil suction (mm) + soilpsi => soilstate_inst%soilpsi_col , & ! Input: [real(r8) (:,:) ] soil water potential in each soil layer (MPa) + t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevgrnd) + o2stress_sat => ch4_inst%o2stress_sat_col , & ! Input: [real(r8) (:,:) ] Ratio of oxygen available to that demanded by roots, aerobes, & methanotrophs (nlevsoi) + o2stress_unsat => ch4_inst%o2stress_unsat_col , & ! Input: [real(r8) (:,:) ] Ratio of oxygen available to that demanded by roots, aerobes, & methanotrophs (nlevsoi) + finundated => ch4_inst%finundated_col , & ! Input: [real(r8) (:) ] fractional inundated area + decomp_cpools_vr => soilbiogeochem_carbonstate_inst%decomp_cpools_vr_col , & ! Input: [real(r8) (:,:,:) ] (gC/m3) vertically-resolved decomposing (litter, cwd, soil) C pools + pathfrac_decomp_cascade => soilbiogeochem_carbonflux_inst%pathfrac_decomp_cascade_col , & ! Output: [real(r8) (:,:,:) ] what fraction of C leaving a given pool passes through a given transition (frac) + rf_decomp_cascade => soilbiogeochem_carbonflux_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) + w_scalar => soilbiogeochem_carbonflux_inst%w_scalar_col , & ! Output: [real(r8) (:,:) ] soil water scalar for decomp + o_scalar => soilbiogeochem_carbonflux_inst%o_scalar_col , & ! Output: [real(r8) (:,:) ] fraction by which decomposition is limited by anoxia + cn_col => soilbiogeochem_carbonflux_inst%cn_col , & ! Output: [real(r8) (:,:) ] C:N ratio + decomp_k => soilbiogeochem_carbonflux_inst%decomp_k_col , & ! Output: [real(r8) (:,:,:) ] rate for decomposition (1./sec) + spinup_factor => decomp_cascade_con%spinup_factor & ! Input: [real(r8) (:) ] factor for AD spinup associated with each pool + ) + + mino2lim = CNParamsShareInst%mino2lim + + days_per_year = get_curr_days_per_year() + +! ! Set "decomp_depth_efolding" parameter +! decomp_depth_efolding = CNParamsShareInst%decomp_depth_efolding + + ! translate to per-second time constant + k_frag = 1._r8 / (secspday * days_per_year * CNParamsShareInst%tau_cwd) + + ! calc ref rate + if ( spinup_state >= 1 ) then + do fc = 1,num_soilc + c = filter_soilc(fc) + ! + if ( abs(spinup_factor(i_met_lit) - 1._r8) .gt. eps) then + spinup_geogterm_l1(c) = spinup_factor(i_met_lit) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_l1(c) = 1._r8 + endif + ! + if ( abs(spinup_factor(i_str_lit) - 1._r8) .gt. eps) then + spinup_geogterm_l2(c) = spinup_factor(i_str_lit) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_l2(c) = 1._r8 + endif + ! + if ( .not. use_fates ) then + if ( abs(spinup_factor(i_cwd) - 1._r8) .gt. eps) then + spinup_geogterm_cwd(c) = spinup_factor(i_cwd) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_cwd(c) = 1._r8 + endif + endif + ! + if ( abs(spinup_factor(i_avl_som) - 1._r8) .gt. eps) then + spinup_geogterm_s1(c) = spinup_factor(i_avl_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_s1(c) = 1._r8 + endif + ! + if ( abs(spinup_factor(i_chem_som) - 1._r8) .gt. eps) then + spinup_geogterm_s2(c) = spinup_factor(i_chem_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_s2(c) = 1._r8 + endif + ! + if ( abs(spinup_factor(i_phys_som) - 1._r8) .gt. eps) then + spinup_geogterm_s3(c) = spinup_factor(i_phys_som) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_s3(c) = 1._r8 + endif + ! + if ( abs(spinup_factor(i_cop_mic) - 1._r8) .gt. eps) then + spinup_geogterm_m1(c) = spinup_factor(i_cop_mic) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_m1(c) = 1._r8 + endif + ! + if ( abs(spinup_factor(i_oli_mic) - 1._r8) .gt. eps) then + spinup_geogterm_m2(c) = spinup_factor(i_oli_mic) * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) + else + spinup_geogterm_m2(c) = 1._r8 + endif + ! + end do + else + do fc = 1,num_soilc + c = filter_soilc(fc) + spinup_geogterm_l1(c) = 1._r8 + spinup_geogterm_l2(c) = 1._r8 + spinup_geogterm_cwd(c) = 1._r8 + spinup_geogterm_s1(c) = 1._r8 + spinup_geogterm_s2(c) = 1._r8 + spinup_geogterm_s3(c) = 1._r8 + spinup_geogterm_m1(c) = 1._r8 + spinup_geogterm_m2(c) = 1._r8 + end do + endif + + !--- time dependent coefficients-----! + if ( nlevdecomp .eq. 1 ) then + + ! calculate function to weight the temperature and water potential scalars + ! for decomposition control. + + + ! the following normalizes values in fr so that they + ! sum to 1.0 across top nlevdecomp levels on a column + frw(bounds%begc:bounds%endc) = 0._r8 + allocate(fr(bounds%begc:bounds%endc,nlev_soildecomp_standard)) + do j=1,nlev_soildecomp_standard + do fc = 1,num_soilc + c = filter_soilc(fc) + frw(c) = frw(c) + col%dz(c,j) + end do + end do + do j = 1,nlev_soildecomp_standard + do fc = 1,num_soilc + c = filter_soilc(fc) + if (frw(c) /= 0._r8) then + fr(c,j) = col%dz(c,j) / frw(c) + else + fr(c,j) = 0._r8 + end if + end do + end do + + ! calculate the rate constant scalar for soil water content. + ! Uses the log relationship with water potential given in + ! Andren, O., and K. Paustian, 1987. Barley straw decomposition in the field: + ! a comparison of models. Ecology, 68(5):1190-1200. + ! and supported by data in + ! Orchard, V.A., and F.J. Cook, 1983. Relationship between soil respiration + ! and soil moisture. Soil Biol. Biochem., 15(4):447-453. + + do j = 1,nlev_soildecomp_standard + do fc = 1,num_soilc + c = filter_soilc(fc) + if (j==1) w_scalar(c,:) = 0._r8 + psi = min(soilpsi(c,j),maxpsi) + ! decomp only if soilpsi is higher than minpsi + if (psi > minpsi) then + w_scalar(c,1) = w_scalar(c,1) + (log(minpsi/psi)/log(minpsi/maxpsi))*fr(c,j) + end if + end do + end do + + ! Calculate ANOXIA + ! anoxia = .true. when (use_lch4) + + if (anoxia) then + + do j = 1,nlev_soildecomp_standard + do fc = 1,num_soilc + c = filter_soilc(fc) + + if (j==1) o_scalar(c,:) = 0._r8 + + o_scalar(c,1) = o_scalar(c,1) + fr(c,j) * max(o2stress_unsat(c,j), mino2lim) + end do + end do + else + o_scalar(bounds%begc:bounds%endc,1:nlevdecomp) = 1._r8 + end if + + deallocate(fr) + + else + + ! calculate the rate constant scalar for soil water content. + ! Uses the log relationship with water potential given in + ! Andren, O., and K. Paustian, 1987. Barley straw decomposition in the field: + ! a comparison of models. Ecology, 68(5):1190-1200. + ! and supported by data in + ! Orchard, V.A., and F.J. Cook, 1983. Relationship between soil respiration + ! and soil moisture. Soil Biol. Biochem., 15(4):447-453. + + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + psi = min(soilpsi(c,j),maxpsi) + ! decomp only if soilpsi is higher than minpsi + if (psi > minpsi) then + w_scalar(c,j) = (log(minpsi/psi)/log(minpsi/maxpsi)) + else + w_scalar(c,j) = 0._r8 + end if + end do + end do + + ! Calculate ANOXIA + ! anoxia = .true. when (use_lch4) + + if (anoxia) then + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + + o_scalar(c,j) = max(o2stress_unsat(c,j), mino2lim) + end do + end do + else + o_scalar(bounds%begc:bounds%endc,1:nlevdecomp) = 1._r8 + end if + + end if + + ! Term that reduces decomposition rate at depth + ! Placeholder. For now depth_scalar = 1. + do j = 1, nlevdecomp + do fc = 1, num_soilc + c = filter_soilc(fc) + ! Using fixed e-folding depth as in + ! SoilBiogeochemDecompCascadeBGCMod.F90 +! depth_scalar(c,j) = exp(-zsoi(j) / decomp_depth_efolding) + depth_scalar(c,j) = 1.0_r8 + end do + end do + + ! TODO @ekluzek suggested possibly making the Left Hand Sides into arrays + ! and I wonder in that case whether to skip these assignments altogether + ! and use the Right Hand Sides directly + mimics_fmet_p1 = params_inst%mimics_fmet(1) + mimics_fmet_p2 = params_inst%mimics_fmet(2) + mimics_fmet_p3 = params_inst%mimics_fmet(3) + mimics_fmet_p4 = params_inst%mimics_fmet(4) + mimics_fchem_r_p1 = params_inst%mimics_fchem_r(1) + mimics_fchem_r_p2 = params_inst%mimics_fchem_r(2) + mimics_fchem_r_p3 = params_inst%mimics_fchem_r(3) + mimics_fchem_k_p1 = params_inst%mimics_fchem_k(1) + mimics_fchem_k_p2 = params_inst%mimics_fchem_k(2) + mimics_fchem_k_p3 = params_inst%mimics_fchem_k(3) + mimics_tau_mod_min = params_inst%mimics_tau_mod_min + mimics_tau_mod_max = params_inst%mimics_tau_mod_max + mimics_tau_mod_factor = params_inst%mimics_tau_mod_factor + mimics_tau_r_p1 = params_inst%mimics_tau_r(1) + mimics_tau_r_p2 = params_inst%mimics_tau_r(2) + mimics_tau_k_p1 = params_inst%mimics_tau_k(1) + mimics_tau_k_p2 = params_inst%mimics_tau_k(2) + mimics_ko_r = params_inst%mimics_ko_r + mimics_ko_k = params_inst%mimics_ko_k + mimics_densdep = params_inst%mimics_densdep + mimics_desorpQ10 = params_inst%mimics_desorpQ10 + mimics_t_soi_ref = params_inst%mimics_t_soi_ref + mimics_cn_mod_num = params_inst%mimics_cn_mod_num + mimics_cn_r = params_inst%mimics_cn_r + mimics_cn_k = params_inst%mimics_cn_k + + ! calculate rates for all litter and som pools + do fc = 1,num_soilc + c = filter_soilc(fc) + + ! Time-dependent params from Wieder et al. 2015 & testbed code + + ! Set limits on Lignin:N to keep fmet > 0.2 + ! Necessary for litter quality in boreal forests with high cwd flux + ! TODO Check for high-freq variations in ligninNratioAvg. To avoid, + ! replace pool_to_litter terms with ann or other long term mean + ! in CNVegCarbonFluxType. + fmet = mimics_fmet_p1 * (mimics_fmet_p2 - mimics_fmet_p3 * min(mimics_fmet_p4, ligninNratioAvg(c))) + tau_mod = min(mimics_tau_mod_max, max(mimics_tau_mod_min, & + sqrt(mimics_tau_mod_factor * & + annsum_npp_col(c)))) + + ! tau_m1 is tauR and tau_m2 is tauK in Wieder et al. 2015 + ! tau ends up in units of per hour but is expected + ! in units of per second, so convert here; alternatively + ! place the conversion once in w_d_o_scalars + tau_m1 = mimics_tau_r_p1 * exp(mimics_tau_r_p2 * fmet) * tau_mod / & + secsphr + tau_m2 = mimics_tau_k_p1 * exp(mimics_tau_k_p2 * fmet) * tau_mod / & + secsphr + + ! These two get used in SoilBiogeochemPotentialMod.F90 + ! cn(c,i_cop_mic), cn(c,i_oli_mic) are CN_r, CN_k in the testbed code + ! cn_r, cn_k are CNr, CNk in the testbed code + ! https://github.com/wwieder/biogeochem_testbed_1.1/blob/Testbed_CN/SOURCE_CODE/mimics_cycle_CN.f90#L1613 + cn_col(c,i_cop_mic) = mimics_cn_r * sqrt(mimics_cn_mod_num / fmet) + cn_col(c,i_oli_mic) = mimics_cn_k * sqrt(mimics_cn_mod_num / fmet) + + ! Used in the update of certain pathfrac terms that vary with time + ! in the next loop + fchem_m1 = min(1._r8, max(0._r8, mimics_fchem_r_p1 * & + exp(mimics_fchem_r_p2 * fmet) * mimics_fchem_r_p3)) + fchem_m2 = min(1._r8, max(0._r8, mimics_fchem_k_p1 * & + exp(mimics_fchem_k_p2 * fmet) * mimics_fchem_k_p3)) + + do j = 1,nlevdecomp + ! vmax ends up in units of per hour but is expected + ! in units of per second, so convert here; alternatively + ! place the conversion once in w_d_o_scalars + ! Table B1 Wieder et al. 2015 & MIMICS params file give diff + ! kslope. I used the params file value(s). + t_soi_degC = t_soisno(c,j) - tfrz + + vmax_l1_m1 = exp(vslope_l1_m1 * t_soi_degC + vint_l1_m1) * & + vmod_l1_m1 / secsphr + vmax_l1_m2 = exp(vslope_l1_m2 * t_soi_degC + vint_l1_m2) * & + vmod_l1_m2 / secsphr + vmax_l2_m1 = exp(vslope_l2_m1 * t_soi_degC + vint_l2_m1) * & + vmod_l2_m1 / secsphr + vmax_l2_m2 = exp(vslope_l2_m2 * t_soi_degC + vint_l2_m2) * & + vmod_l2_m2 / secsphr + vmax_s1_m1 = exp(vslope_s1_m1 * t_soi_degC + vint_s1_m1) * & + vmod_s1_m1 / secsphr + vmax_s1_m2 = exp(vslope_s1_m2 * t_soi_degC + vint_s1_m2) * & + vmod_s1_m2 / secsphr + + km_l1_m1 = exp(kslope_l1_m1 * t_soi_degC + kint_l1_m1) * & + kmod_l1_m1 + km_l1_m2 = exp(kslope_l1_m2 * t_soi_degC + kint_l1_m2) * & + kmod_l1_m2 + km_l2_m1 = exp(kslope_l2_m1 * t_soi_degC + kint_l2_m1) * & + kmod_l2_m1 + km_l2_m2 = exp(kslope_l2_m2 * t_soi_degC + kint_l2_m2) * & + kmod_l2_m2 + km_s1_m1 = exp(kslope_s1_m1 * t_soi_degC + kint_s1_m1) * & + kmod_s1_m1 * p_scalar(c,j) + km_s1_m2 = exp(kslope_s1_m2 * t_soi_degC + kint_s1_m2) * & + kmod_s1_m2 * p_scalar(c,j) + + ! Desorption a function of soil temperature and + ! Q10 = 1.1 w/ reference temperature of 25C. + ! Expected in units of per second, so convert; alternatively + ! place the conversion once in w_d_o_scalars + desorption = (desorp(c,j) / secsphr) * mimics_desorpQ10 * & + exp((t_soi_degC - mimics_t_soi_ref) / 10.0_r8) + + ! Microbial concentration with necessary unit conversions + ! mgC/cm3 = gC/m3 * (1e3 mg/g) / (1e6 cm3/m3) + m1_conc = (decomp_cpools_vr(c,j,i_cop_mic) / col%dz(c,j)) * & + g_to_mg * cm3_to_m3 + m2_conc = (decomp_cpools_vr(c,j,i_oli_mic) / col%dz(c,j)) * & + g_to_mg * cm3_to_m3 + + ! Product of w_scalar * depth_scalar * o_scalar + w_d_o_scalars = w_scalar(c,j) * depth_scalar(c,j) * o_scalar(c,j) + + ! decomp_k used in SoilBiogeochemPotentialMod.F90 + ! also updating pathfrac terms that vary with time + term_1 = vmax_l1_m1 * m1_conc / (km_l1_m1 + m1_conc) + term_2 = vmax_l1_m2 * m2_conc / (km_l1_m2 + m2_conc) + decomp_k(c,j,i_met_lit) = (term_1 + term_2) * w_d_o_scalars + if (term_1 + term_2 /= 0._r8) then + pathfrac_decomp_cascade(c,j,i_l1m1) = term_1 / (term_1 + term_2) + pathfrac_decomp_cascade(c,j,i_l1m2) = term_2 / (term_1 + term_2) + else + pathfrac_decomp_cascade(c,j,i_l1m1) = 0._r8 + pathfrac_decomp_cascade(c,j,i_l1m2) = 0._r8 + end if + + term_1 = vmax_l2_m1 * m1_conc / (km_l2_m1 + m1_conc) + term_2 = vmax_l2_m2 * m2_conc / (km_l2_m2 + m2_conc) + decomp_k(c,j,i_str_lit) = (term_1 + term_2) * w_d_o_scalars + if (term_1 + term_2 /= 0._r8) then + pathfrac_decomp_cascade(c,j,i_l2m1) = term_1 / (term_1 + term_2) + pathfrac_decomp_cascade(c,j,i_l2m2) = term_2 / (term_1 + term_2) + else + pathfrac_decomp_cascade(c,j,i_l2m1) = 0._r8 + pathfrac_decomp_cascade(c,j,i_l2m2) = 0._r8 + end if + + term_1 = vmax_s1_m1 * m1_conc / (km_s1_m1 + m1_conc) + term_2 = vmax_s1_m2 * m2_conc / (km_s1_m2 + m2_conc) + decomp_k(c,j,i_avl_som) = (term_1 + term_2) * w_d_o_scalars + if (term_1 + term_2 /= 0._r8) then + pathfrac_decomp_cascade(c,j,i_s1m1) = term_1 / (term_1 + term_2) + pathfrac_decomp_cascade(c,j,i_s1m2) = term_2 / (term_1 + term_2) + else + pathfrac_decomp_cascade(c,j,i_s1m1) = 0._r8 + pathfrac_decomp_cascade(c,j,i_s1m2) = 0._r8 + end if + + decomp_k(c,j,i_phys_som) = desorption * depth_scalar(c,j) + + term_1 = vmax_l2_m1 * m1_conc / (mimics_ko_r * km_l2_m1 + m1_conc) + term_2 = vmax_l2_m2 * m2_conc / (mimics_ko_k * km_l2_m2 + m2_conc) + ! The right hand side is OXIDAT in the testbed (line 1145) + decomp_k(c,j,i_chem_som) = (term_1 + term_2) * w_d_o_scalars + + decomp_k(c,j,i_cop_mic) = tau_m1 * & + m1_conc**(mimics_densdep - 1.0_r8) * w_d_o_scalars + favl = min(1.0_r8, max(0.0_r8, 1.0_r8 - fphys_m1(c,j) - fchem_m1)) + pathfrac_decomp_cascade(c,j,i_m1s1) = favl + pathfrac_decomp_cascade(c,j,i_m1s2) = fchem_m1 + + decomp_k(c,j,i_oli_mic) = tau_m2 * & + m2_conc**(mimics_densdep - 1.0_r8) * w_d_o_scalars + favl = min(1.0_r8, max(0.0_r8, 1.0_r8 - fphys_m2(c,j) - fchem_m2)) + pathfrac_decomp_cascade(c,j,i_m2s1) = favl + pathfrac_decomp_cascade(c,j,i_m2s2) = fchem_m2 + + ! Same for cwd but only if fates not enabled; fates handles cwd on + ! its own structure + ! TODO This shows how BGC applies the spinup coefficients + if (.not. use_fates) then + decomp_k(c,j,i_cwd) = k_frag * w_d_o_scalars ! * spinup_geogterm_cwd(c) + end if + end do + end do + + ! pathfrac terms not calculated in the previous loop + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = 1.0_r8 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = 1.0_r8 + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m1s3) = fphys_m1(bounds%begc:bounds%endc,1:nlevdecomp) + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m2s3) = fphys_m2(bounds%begc:bounds%endc,1:nlevdecomp) + if (.not. use_fates) then + pathfrac_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = 1.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_cwdl2) = rf_cwdl2 + end if + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1m1) = rf_l1m1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l1m2) = rf_l1m2 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2m1) = rf_l2m1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_l2m2) = rf_l2m2 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1m1) = rf_s1m1 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s1m2) = rf_s1m2 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s2s1) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_s3s1) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m1s1) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m1s2) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m1s3) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m2s1) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m2s2) = 0.0_r8 + rf_decomp_cascade(bounds%begc:bounds%endc,1:nlevdecomp,i_m2s3) = 0.0_r8 + + end associate + + end subroutine decomp_rates_mimics + +end module SoilBiogeochemDecompCascadeMIMICSMod diff --git a/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 index 43b817c65b..b59db9c5ce 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemDecompMod.F90 @@ -11,9 +11,9 @@ 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_fates + use clm_varctl , only : use_nitrif_denitrif, use_lch4, use_fates, iulog use clm_varcon , only : dzsoi_decomp - use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con + use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, mimics_decomp, decomp_method use SoilBiogeochemStateType , only : soilbiogeochem_state_type use SoilBiogeochemCarbonStateType , only : soilbiogeochem_carbonstate_type use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type @@ -70,7 +70,8 @@ end subroutine readParams subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & - cn_decomp_pools, p_decomp_cpool_loss, pmnf_decomp_cascade) + cn_decomp_pools, p_decomp_cpool_loss, pmnf_decomp_cascade, & + p_decomp_npool_to_din) ! ! !USES: use SoilBiogeochemDecompCascadeConType, only : i_atm @@ -87,6 +88,7 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, real(r8) , intent(inout) :: cn_decomp_pools(bounds%begc:,1:,1:) ! c:n ratios of applicable pools real(r8) , intent(inout) :: p_decomp_cpool_loss(bounds%begc:,1:,1:) ! potential C loss from one pool to another real(r8) , intent(inout) :: pmnf_decomp_cascade(bounds%begc:,1:,1:) ! potential mineral N flux from one pool to another + real(r8) , intent(in) :: p_decomp_npool_to_din(bounds%begc:,1:,1:) ! potential flux to dissolved inorganic N ! ! !LOCAL VARIABLES: integer :: c,j,k,l,m ! indices @@ -101,6 +103,7 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, SHR_ASSERT_ALL_FL((ubound(cn_decomp_pools) == (/endc,nlevdecomp,ndecomp_pools/)) , sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(p_decomp_cpool_loss) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(pmnf_decomp_cascade) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(p_decomp_npool_to_din) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) associate( & cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step @@ -109,8 +112,7 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, initial_cn_ratio => decomp_cascade_con%initial_cn_ratio , & ! Input: [real(r8) (:) ] c:n ratio for initialization of pools fpi_vr => soilbiogeochem_state_inst%fpi_vr_col , & ! Input: [real(r8) (:,:) ] fraction of potential immobilization (no units) - rf_decomp_cascade => soilbiogeochem_state_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) - pathfrac_decomp_cascade => soilbiogeochem_state_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] what fraction of C leaving a given pool passes through a given transition (frac) + rf_decomp_cascade => soilbiogeochem_carbonflux_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) decomp_npools_vr => soilbiogeochem_nitrogenstate_inst%decomp_npools_vr_col , & ! Input: [real(r8) (:,:,:) ] (gC/m3) vertically-resolved decomposing (litter, cwd, soil) N pools decomp_cpools_vr => soilbiogeochem_carbonstate_inst%decomp_cpools_vr_col , & ! Input: [real(r8) (:,:,:) ] (gC/m3) vertically-resolved decomposing (litter, cwd, soil) c pools @@ -119,12 +121,13 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, decomp_cascade_sminn_flux_vr => soilbiogeochem_nitrogenflux_inst%decomp_cascade_sminn_flux_vr_col , & ! Output: [real(r8) (:,:,:) ] vert-res mineral N flux for transition along decomposition cascade (gN/m3/s) potential_immob_vr => soilbiogeochem_nitrogenflux_inst%potential_immob_vr_col , & ! Output: [real(r8) (:,:) ] sminn_to_denit_decomp_cascade_vr => soilbiogeochem_nitrogenflux_inst%sminn_to_denit_decomp_cascade_vr_col , & ! Output: [real(r8) (:,:,:) ] - gross_nmin_vr => soilbiogeochem_nitrogenflux_inst%gross_nmin_vr_col , & ! Output: [real(r8) (:,:) ] + gross_nmin_vr => soilbiogeochem_nitrogenflux_inst%gross_nmin_vr_col , & ! Input: [real(r8) (:,:) ] net_nmin_vr => soilbiogeochem_nitrogenflux_inst%net_nmin_vr_col , & ! Output: [real(r8) (:,:) ] gross_nmin => soilbiogeochem_nitrogenflux_inst%gross_nmin_col , & ! Output: [real(r8) (:) ] gross rate of N mineralization (gN/m2/s) net_nmin => soilbiogeochem_nitrogenflux_inst%net_nmin_col , & ! Output: [real(r8) (:) ] net rate of N mineralization (gN/m2/s) w_scalar => soilbiogeochem_carbonflux_inst%w_scalar_col , & ! Input: [real(r8) (:,:) ] fraction by which decomposition is limited by moisture availability + c_overflow_vr => soilbiogeochem_carbonflux_inst%c_overflow_vr , & ! Input: [real(r8) (:,:,:) ] vertically-resolved C rejected by microbes that cannot process it (gC/m3/s) decomp_cascade_hr_vr => soilbiogeochem_carbonflux_inst%decomp_cascade_hr_vr_col , & ! Output: [real(r8) (:,:,:) ] vertically-resolved het. resp. from decomposing C pools (gC/m3/s) decomp_cascade_ctransfer_vr => soilbiogeochem_carbonflux_inst%decomp_cascade_ctransfer_vr_col , & ! Output: [real(r8) (:,:,:) ] vertically-resolved het. resp. from decomposing C pools (gC/m3/s) phr_vr => soilbiogeochem_carbonflux_inst%phr_vr_col , & ! Input: [real(r8) (:,:) ] potential HR (gC/m3/s) @@ -184,6 +187,12 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, end if 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) + if (decomp_method == mimics_decomp) then + decomp_cascade_hr_vr(c,j,k) = min( & + p_decomp_cpool_loss(c,j,k), & + decomp_cascade_hr_vr(c,j,k) + c_overflow_vr(c,j,k)) + decomp_cascade_ctransfer_vr(c,j,k) = max(0.0_r8, p_decomp_cpool_loss(c,j,k) - decomp_cascade_hr_vr(c,j,k)) + end if if (decomp_npools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. cascade_receiver_pool(k) /= i_atm) then decomp_cascade_ntransfer_vr(c,j,k) = p_decomp_cpool_loss(c,j,k) / cn_decomp_pools(c,j,cascade_donor_pool(k)) else @@ -195,6 +204,10 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_soilc, decomp_cascade_sminn_flux_vr(c,j,k) = - pmnf_decomp_cascade(c,j,k) endif net_nmin_vr(c,j) = net_nmin_vr(c,j) - pmnf_decomp_cascade(c,j,k) + if (decomp_method == mimics_decomp) then + decomp_cascade_sminn_flux_vr(c,j,k) = decomp_cascade_sminn_flux_vr(c,j,k) - p_decomp_npool_to_din(c,j,k) + net_nmin_vr(c,j) = net_nmin_vr(c,j) + p_decomp_npool_to_din(c,j,k) + end if else decomp_cascade_ntransfer_vr(c,j,k) = 0._r8 if (.not. use_nitrif_denitrif) then @@ -213,8 +226,13 @@ subroutine SoilBiogeochemDecomp (bounds, num_soilc, filter_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) + if (decomp_method == mimics_decomp) then + decomp_cascade_hr_vr(c,j,k) = min( & + p_decomp_cpool_loss(c,j,k), & + decomp_cascade_hr_vr(c,j,k) + c_overflow_vr(c,j,k)) + decomp_cascade_ctransfer_vr(c,j,k) = max(0.0_r8, p_decomp_cpool_loss(c,j,k) - decomp_cascade_hr_vr(c,j,k)) + end if ! end do end do diff --git a/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90 b/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90 index 613d02c608..da5e2d4240 100644 --- a/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNitrogenStateType.F90 @@ -12,7 +12,7 @@ module SoilBiogeochemNitrogenStateType use clm_varpar , only : nlevdecomp_full, nlevdecomp, nlevsoi use clm_varcon , only : spval, dzsoi_decomp, zisoi use clm_varctl , only : use_nitrif_denitrif - use SoilBiogeochemDecompCascadeConType , only : century_decomp, decomp_method + use SoilBiogeochemDecompCascadeConType , only : mimics_decomp, century_decomp, decomp_method use clm_varctl , only : iulog, override_bgc_restart_mismatch_dump, spinup_state use landunit_varcon , only : istcrop, istsoil use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con @@ -46,6 +46,7 @@ module SoilBiogeochemNitrogenStateType real(r8), pointer :: ntrunc_col (:) ! col (gN/m2) column-level sink for N truncation real(r8), pointer :: cwdn_col (:) ! col (gN/m2) Diagnostic: coarse woody debris N real(r8), pointer :: totlitn_col (:) ! col (gN/m2) total litter nitrogen + real(r8), pointer :: totmicn_col (:) ! col (gN/m2) total microbial nitrogen real(r8), pointer :: totsomn_col (:) ! col (gN/m2) total soil organic matter nitrogen real(r8), pointer :: totlitn_1m_col (:) ! col (gN/m2) total litter nitrogen to 1 meter real(r8), pointer :: totsomn_1m_col (:) ! col (gN/m2) total soil organic matter nitrogen to 1 meter @@ -120,6 +121,7 @@ subroutine InitAllocate(this, bounds) allocate(this%sminn_col (begc:endc)) ; this%sminn_col (:) = nan allocate(this%ntrunc_col (begc:endc)) ; this%ntrunc_col (:) = nan allocate(this%totlitn_col (begc:endc)) ; this%totlitn_col (:) = nan + allocate(this%totmicn_col (begc:endc)) ; this%totmicn_col (:) = nan allocate(this%totsomn_col (begc:endc)) ; this%totsomn_col (:) = nan allocate(this%totlitn_1m_col (begc:endc)) ; this%totlitn_1m_col (:) = nan allocate(this%totsomn_1m_col (begc:endc)) ; this%totsomn_1m_col (:) = nan @@ -277,6 +279,11 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='total litter N', & ptr_col=this%totlitn_col) + this%totmicn_col(begc:endc) = spval + call hist_addfld1d (fname='TOTMICN', units='gN/m^2', & + avgflag='A', long_name='total microbial N', & + ptr_col=this%totmicn_col) + this%totsomn_col(begc:endc) = spval call hist_addfld1d (fname='TOTSOMN', units='gN/m^2', & avgflag='A', long_name='total soil organic matter N', & @@ -368,6 +375,7 @@ subroutine InitCold(this, bounds, & this%smin_no3_col(c) = 0._r8 end if this%totlitn_col(c) = 0._r8 + this%totmicn_col(c) = 0._r8 this%totsomn_col(c) = 0._r8 this%totlitn_1m_col(c) = 0._r8 this%totsomn_1m_col(c) = 0._r8 @@ -495,6 +503,8 @@ subroutine Restart ( this, bounds, ncid, flag, totvegc_col ) if (decomp_method == century_decomp ) then decomp_cascade_state = 1 + else if (decomp_method == mimics_decomp ) then + decomp_cascade_state = 2 else decomp_cascade_state = 0 end if @@ -656,6 +666,7 @@ subroutine SetValues ( this, num_column, filter_column, value_column ) this%smin_nh4_col(i) = value_column end if this%totlitn_col(i) = value_column + this%totmicn_col(i) = value_column this%totsomn_col(i) = value_column this%totsomn_1m_col(i) = value_column this%totlitn_1m_col(i) = value_column @@ -849,6 +860,22 @@ subroutine Summary(this, bounds, num_allc, filter_allc) end if end do + ! total microbial nitrogen (TOTMICN) + do fc = 1,num_allc + c = filter_allc(fc) + this%totmicn_col(c) = 0._r8 + end do + do l = 1, ndecomp_pools + if ( decomp_cascade_con%is_microbe(l) ) then + do fc = 1,num_allc + c = filter_allc(fc) + this%totmicn_col(c) = & + this%totmicn_col(c) + & + this%decomp_npools_col(c,l) + end do + end if + end do + ! total soil organic matter nitrogen (TOTSOMN) do fc = 1,num_allc c = filter_allc(fc) diff --git a/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 b/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 index d7d6ed727f..da46e178b7 100644 --- a/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemPotentialMod.F90 @@ -10,7 +10,8 @@ module SoilBiogeochemPotentialMod use shr_kind_mod , only : r8 => shr_kind_r8 use decompMod , only : bounds_type use clm_varpar , only : nlevdecomp, ndecomp_cascade_transitions, ndecomp_pools - use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con + use clm_varpar , only : i_cop_mic, i_oli_mic + use SoilBiogeochemDecompCascadeConType , only : decomp_cascade_con, mimics_decomp, decomp_method use SoilBiogeochemStateType , only : soilbiogeochem_state_type use SoilBiogeochemCarbonStateType , only : soilbiogeochem_carbonstate_type use SoilBiogeochemCarbonFluxType , only : soilbiogeochem_carbonflux_type @@ -70,7 +71,8 @@ end subroutine readParams subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, soilbiogeochem_carbonflux_inst, & soilbiogeochem_nitrogenstate_inst, soilbiogeochem_nitrogenflux_inst, & - cn_decomp_pools, p_decomp_cpool_loss, pmnf_decomp_cascade) + cn_decomp_pools, p_decomp_cpool_loss, p_decomp_cn_gain, & + pmnf_decomp_cascade, p_decomp_npool_to_din) ! ! !USES: use shr_log_mod , only : errMsg => shr_log_errMsg @@ -88,20 +90,31 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & real(r8) , intent(out) :: cn_decomp_pools(bounds%begc:,1:,1:) ! c:n ratios of applicable pools real(r8) , intent(out) :: p_decomp_cpool_loss(bounds%begc:,1:,1:) ! potential C loss from one pool to another real(r8) , intent(out) :: pmnf_decomp_cascade(bounds%begc:,1:,1:) ! potential mineral N flux, from one pool to another + real(r8) , intent(out) :: p_decomp_npool_to_din(bounds%begc:,1:,1:) ! potential flux to dissolved inorganic N + real(r8) , intent(out) :: p_decomp_cn_gain(bounds%begc:,1:,1:) ! C:N ratio of the flux gained by the receiver pool ! ! !LOCAL VARIABLES: integer :: c,j,k,l,m !indices integer :: fc !filter column index integer :: begc,endc !bounds real(r8):: immob(bounds%begc:bounds%endc,1:nlevdecomp) !potential N immobilization + real(r8):: p_decomp_cpool_gain(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) ! potential C gain by receiver pool + real(r8):: p_decomp_npool_gain(bounds%begc:bounds%endc,1:nlevdecomp,1:ndecomp_cascade_transitions) ! potential N gain by receiver pool + real(r8):: p_decomp_cpool_gain_sum(1:ndecomp_pools) ! total potential C gain by receiver pool (only microbial pools) + real(r8):: p_decomp_npool_gain_sum(1:ndecomp_pools) ! total potential N gain by receiver pool (only microbial pools) + real(r8):: decomp_nc_loss_donor ! N:C ratio of donor pool + real(r8):: p_decomp_cn_diff_ratio ! relative change in receiver pool C:N + real(r8):: p_decomp_npool_loss ! potential N flux out of donor pool real(r8):: ratio !temporary variable !----------------------------------------------------------------------- begc = bounds%begc; endc = bounds%endc SHR_ASSERT_ALL_FL((ubound(cn_decomp_pools) == (/endc,nlevdecomp,ndecomp_pools/)) , sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(p_decomp_cn_gain) == (/endc,nlevdecomp,ndecomp_pools/)) , sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(p_decomp_cpool_loss) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(pmnf_decomp_cascade) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(p_decomp_npool_to_din) == (/endc,nlevdecomp,ndecomp_cascade_transitions/)) , sourcefile, __LINE__) associate( & cascade_donor_pool => decomp_cascade_con%cascade_donor_pool , & ! Input: [integer (:) ] which pool is C taken from for a given decomposition step @@ -109,8 +122,9 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & floating_cn_ratio_decomp_pools => decomp_cascade_con%floating_cn_ratio_decomp_pools , & ! Input: [logical (:) ] TRUE => pool has fixed C:N ratio initial_cn_ratio => decomp_cascade_con%initial_cn_ratio , & ! Input: [real(r8) (:) ] c:n ratio for initialization of pools - rf_decomp_cascade => soilbiogeochem_state_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) - pathfrac_decomp_cascade => soilbiogeochem_state_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] what fraction of C leaving a given pool passes through a given transition (frac) + nue_decomp_cascade => soilbiogeochem_state_inst%nue_decomp_cascade_col , & ! Input: [real(r8) (:) ] N use efficiency for a given transition (gN going into microbe / gN decomposed) + rf_decomp_cascade => soilbiogeochem_carbonflux_inst%rf_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] respired fraction in decomposition step (frac) + pathfrac_decomp_cascade => soilbiogeochem_carbonflux_inst%pathfrac_decomp_cascade_col , & ! Input: [real(r8) (:,:,:) ] what fraction of C passes from donor to receiver pool through a given transition (frac) decomp_npools_vr => soilbiogeochem_nitrogenstate_inst%decomp_npools_vr_col , & ! Input: [real(r8) (:,:,:) ] (gC/m3) vertically-resolved decomposing (litter, cwd, soil) N pools @@ -119,7 +133,8 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & potential_immob_vr => soilbiogeochem_nitrogenflux_inst%potential_immob_vr_col , & ! Output: [real(r8) (:,:) ] gross_nmin_vr => soilbiogeochem_nitrogenflux_inst%gross_nmin_vr_col , & ! Output: [real(r8) (:,:) ] - decomp_k => soilbiogeochem_carbonflux_inst%decomp_k_col , & ! Output: [real(r8) (:,:,:) ] rate constant for decomposition (1./sec) + cn_col => soilbiogeochem_carbonflux_inst%cn_col , & ! Input: [real(r8) (:,:) ] C:N ratio + decomp_k => soilbiogeochem_carbonflux_inst%decomp_k_col , & ! Input: [real(r8) (:,:,:) ] decomposition rate coefficient (1./sec) phr_vr => soilbiogeochem_carbonflux_inst%phr_vr_col & ! Output: [real(r8) (:,:) ] potential HR (gC/m3/s) ) @@ -165,6 +180,7 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & decomp_k(c,j,cascade_donor_pool(k)) > 0._r8 ) then 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) + if ( .not. floating_cn_ratio_decomp_pools(cascade_receiver_pool(k)) ) then !! not transition of cwd to litter if (cascade_receiver_pool(k) /= i_atm ) then ! not 100% respiration @@ -180,15 +196,99 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & else ! 100% respiration pmnf_decomp_cascade(c,j,k) = - p_decomp_cpool_loss(c,j,k) / cn_decomp_pools(c,j,cascade_donor_pool(k)) endif - - else ! CWD -> litter + else ! CWD -> litter OR mimics_decomp is true pmnf_decomp_cascade(c,j,k) = 0._r8 - end if - end if - end do - end do - end do + if (decomp_method == mimics_decomp) then + ! N:C ratio of donor pools (N:C instead of C:N because + ! already checked that we're not dividing by zero) + decomp_nc_loss_donor = & + decomp_npools_vr(c,j,cascade_donor_pool(k)) / & + decomp_cpools_vr(c,j,cascade_donor_pool(k)) + ! calculate N fluxes from donor pools + p_decomp_npool_loss = & + p_decomp_cpool_loss(c,j,k) * decomp_nc_loss_donor + ! Track N lost to DIN from messy eating, ie the + ! diff between p_decomp_npool_loss (pertains to + ! donor) and p_decomp_npool_gain (receiver) + p_decomp_cpool_gain(c,j,k) = & + p_decomp_cpool_loss(c,j,k) * (1.0_r8 - rf_decomp_cascade(c,j,k)) + p_decomp_npool_gain(c,j,k) = & + p_decomp_npool_loss * nue_decomp_cascade(k) + p_decomp_npool_to_din(c,j,k) = & + p_decomp_npool_loss - p_decomp_npool_gain(c,j,k) + end if ! using mimics + end if ! floating_cn_ratio_decomp_pools = .true. + else ! donors not donating + p_decomp_cpool_gain(c,j,k) = 0.0_r8 + p_decomp_npool_gain(c,j,k) = 0.0_r8 + p_decomp_npool_to_din(c,j,k) = 0.0_r8 + end if ! donors donating? (ie decomp_cpools_vr & decomp_k > 0) + end do ! column loop + + end do ! soil level loop + end do ! transitions loop + + ! Calculate cn_gain into microbial biomass (in first do k + ! transitions loop). + ! Compare cn_gain to target C:N ratio of microbial biomass pools + ! to determine immobilization vs. mineralization (in second do k + ! transitions loop). + if (decomp_method == mimics_decomp) then + do j = 1,nlevdecomp + do fc = 1,num_soilc + c = filter_soilc(fc) + ! Sum C & N fluxes from all transitions into m1 & m2 pools. + ! Had to form a new loop for the summation due to the order + ! necessary, ie do k as the innermost loop. + p_decomp_cpool_gain_sum(:) = 0.0_r8 + p_decomp_npool_gain_sum(:) = 0.0_r8 + do k = 1, ndecomp_cascade_transitions + if (cascade_receiver_pool(k) == i_cop_mic .or. & + cascade_receiver_pool(k) == i_oli_mic) then + if (decomp_cpools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. & + decomp_k(c,j,cascade_donor_pool(k)) > 0._r8 ) then + p_decomp_cpool_gain_sum(cascade_receiver_pool(k)) = & + p_decomp_cpool_gain_sum(cascade_receiver_pool(k)) + & + p_decomp_cpool_gain(c,j,k) + p_decomp_npool_gain_sum(cascade_receiver_pool(k)) = & + p_decomp_npool_gain_sum(cascade_receiver_pool(k)) + & + p_decomp_npool_gain(c,j,k) + ! Once k-loop completes, the left hand side should end + ! up with the correct cn ratio + if (p_decomp_npool_gain_sum(cascade_receiver_pool(k)) > 0.0_r8) then + p_decomp_cn_gain(c,j,cascade_receiver_pool(k)) = & + p_decomp_cpool_gain_sum(cascade_receiver_pool(k)) / & + p_decomp_npool_gain_sum(cascade_receiver_pool(k)) + else + p_decomp_cn_gain(c,j,cascade_receiver_pool(k)) = 0.0_r8 + end if ! denominator check + end if ! donors donating (decomp_cpools_vr & decomp_k > 0) + end if ! microbes receiving + end do ! transitions loop + do k = 1, ndecomp_cascade_transitions + if (cascade_receiver_pool(k) == i_cop_mic .or. & + cascade_receiver_pool(k) == i_oli_mic) then + if (decomp_cpools_vr(c,j,cascade_donor_pool(k)) > 0._r8 .and. & + decomp_k(c,j,cascade_donor_pool(k)) > 0._r8 ) then + ! if p_decomp_cn_diff < 0 N mineralization + ! > 0 immobilization + ! "min" in next line turns off immobilization flux + p_decomp_cn_diff_ratio = min(0.0_r8, & + (p_decomp_cn_gain(c,j,cascade_receiver_pool(k)) - & + cn_col(c,cascade_receiver_pool(k))) / cn_col(c,cascade_receiver_pool(k))) + ! Actual amount of N that's mineralized or that would + ! need to be immobilized + ! negative=mineralization: add to the DIN pool + ! positive=immobilizaiton: compete for N with plants to + ! see how much we get + pmnf_decomp_cascade(c,j,k) = p_decomp_cn_diff_ratio * p_decomp_npool_gain(c,j,k) + end if ! donors donating (decomp_cpools_vr & decomp_k > 0) + end if ! microbes receiving + end do ! transitions loop + end do ! column loop + end do ! soil levels loop + end if ! using mimics ! Sum up all the potential immobilization fluxes (positive pmnf flux) ! and all the mineralization fluxes (negative pmnf flux) @@ -207,6 +307,9 @@ subroutine SoilBiogeochemPotential (bounds, num_soilc, filter_soilc, & else gross_nmin_vr(c,j) = gross_nmin_vr(c,j) - pmnf_decomp_cascade(c,j,k) end if + if (decomp_method == mimics_decomp) then + gross_nmin_vr(c,j) = gross_nmin_vr(c,j) + p_decomp_npool_to_din(c,j,k) + end if end do end do end do diff --git a/src/soilbiogeochem/SoilBiogeochemStateType.F90 b/src/soilbiogeochem/SoilBiogeochemStateType.F90 index 308db43eae..fcdced386d 100644 --- a/src/soilbiogeochem/SoilBiogeochemStateType.F90 +++ b/src/soilbiogeochem/SoilBiogeochemStateType.F90 @@ -7,7 +7,7 @@ module SoilBiogeochemStateType use abortutils , only : endrun use spmdMod , only : masterproc use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevsoifl, nlevsoi - use clm_varpar , only : ndecomp_cascade_transitions, nlevdecomp, nlevdecomp_full + use clm_varpar , only : ndecomp_pools, ndecomp_cascade_transitions, nlevdecomp, nlevdecomp_full use clm_varcon , only : spval, ispval, c14ratio, grlnd use landunit_varcon, only : istsoil, istcrop use clm_varpar , only : nlevsno, nlevgrnd, nlevlak @@ -33,8 +33,7 @@ module SoilBiogeochemStateType real(r8) , pointer :: fpi_vr_col (:,:) ! (no units) fraction of potential immobilization real(r8) , pointer :: fpi_col (:) ! (no units) fraction of potential immobilization real(r8), pointer :: fpg_col (:) ! (no units) fraction of potential gpp - real(r8) , pointer :: rf_decomp_cascade_col (:,:,:) ! (frac) respired fraction in decomposition step - real(r8) , pointer :: pathfrac_decomp_cascade_col (:,:,:) ! (frac) what fraction of C leaving a given pool passes through a given transition + real(r8) , pointer :: nue_decomp_cascade_col (:) ! (gN going into microbe / gN decomposed) N use efficiency for a given transition real(r8) , pointer :: nfixation_prof_col (:,:) ! (1/m) profile for N fixation additions real(r8) , pointer :: ndep_prof_col (:,:) ! (1/m) profile for N fixation additions real(r8) , pointer :: som_adv_coef_col (:,:) ! (m2/s) SOM advective flux @@ -102,11 +101,8 @@ subroutine InitAllocate(this, bounds) allocate(this%som_diffus_coef_col (begc:endc,1:nlevdecomp_full)) ; this%som_diffus_coef_col (:,:) = spval allocate(this%plant_ndemand_col (begc:endc)) ; this%plant_ndemand_col (:) = nan - allocate(this%rf_decomp_cascade_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)); - this%rf_decomp_cascade_col(:,:,:) = nan - - allocate(this%pathfrac_decomp_cascade_col(begc:endc,1:nlevdecomp_full,1:ndecomp_cascade_transitions)); - this%pathfrac_decomp_cascade_col(:,:,:) = nan + allocate(this%nue_decomp_cascade_col(1:ndecomp_cascade_transitions)); + this%nue_decomp_cascade_col(:) = nan end subroutine InitAllocate