diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index c80ccab573..3326c4b105 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -2063,7 +2063,6 @@ sub error_if_set { } } - #------------------------------------------------------------------------------- sub setup_logic_soilstate { @@ -3676,6 +3675,11 @@ sub setup_logic_canopyfluxes { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_undercanopy_stability' ); add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'itmax_canopy_fluxes', 'structure'=>$nl_flags->{'structure'}); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_biomass_heat_storage', + 'use_fates'=>$nl_flags->{'use_fates'}, 'phys'=>$nl_flags->{'phys'} ); + if ( &value_is_true($nl->get_value('use_biomass_heat_storage') ) && &value_is_true( $nl_flags->{'use_fates'}) ) { + $log->fatal_error('use_biomass_heat_storage can NOT be set to true when fates is on'); + } } #------------------------------------------------------------------------------- diff --git a/bld/namelist_files/namelist_defaults_ctsm.xml b/bld/namelist_files/namelist_defaults_ctsm.xml index 57793cc83d..b80f529904 100644 --- a/bld/namelist_files/namelist_defaults_ctsm.xml +++ b/bld/namelist_files/namelist_defaults_ctsm.xml @@ -182,7 +182,6 @@ attributes from the config_cache.xml file (with keys converted to upper-case). .true. .false. - 1 1 @@ -312,6 +311,10 @@ attributes from the config_cache.xml file (with keys converted to upper-case). .false. .true. +.true. +.false. +.false. + 40 3 @@ -470,9 +473,9 @@ attributes from the config_cache.xml file (with keys converted to upper-case). -lnd/clm2/paramdata/ctsm51_params.c200905.nc -lnd/clm2/paramdata/clm50_params.c200905.nc -lnd/clm2/paramdata/clm45_params.c200905.nc +lnd/clm2/paramdata/ctsm51_params.c210112.nc +lnd/clm2/paramdata/clm50_params.c210112.nc +lnd/clm2/paramdata/clm45_params.c210112.nc diff --git a/bld/namelist_files/namelist_definition_ctsm.xml b/bld/namelist_files/namelist_definition_ctsm.xml index 23612b2c88..d0db01ddc4 100644 --- a/bld/namelist_files/namelist_definition_ctsm.xml +++ b/bld/namelist_files/namelist_definition_ctsm.xml @@ -331,6 +331,10 @@ Intercept of free living Nitrogen fixation with zero annual ET If TRUE use the undercanopy stability term used with CLM4.5 (Sakaguchi&Zeng, 2008) + +If TRUE, include biomass heat storage in canopy energy balance. + Max number of iterations used in subr. CanopyFluxes. For many years, 40 was the hardwired default value. diff --git a/bld/unit_testers/build-namelist_test.pl b/bld/unit_testers/build-namelist_test.pl index 1a86b5faf6..24c7a0add3 100755 --- a/bld/unit_testers/build-namelist_test.pl +++ b/bld/unit_testers/build-namelist_test.pl @@ -138,9 +138,9 @@ sub make_config_cache { # # Figure out number of tests that will run # -my $ntests = 1513; +my $ntests = 1550; if ( defined($opts{'compare'}) ) { - $ntests += 1017; + $ntests += 1044; } plan( tests=>$ntests ); @@ -837,6 +837,11 @@ sub make_config_cache { GLC_TWO_WAY_COUPLING=>"FALSE", phys=>"clm5_0", }, + "useFATESWbMH" =>{ options=>"-bgc fates -envxml_dir . -no-megan", + namelst=>"use_biomass_heat_storage=.true.", + GLC_TWO_WAY_COUPLING=>"FALSE", + phys=>"clm5_1", + }, "createcropFalse" =>{ options=>"-bgc bgc -envxml_dir . -no-megan", namelst=>"create_crop_landunit=.false.", GLC_TWO_WAY_COUPLING=>"FALSE", @@ -1413,22 +1418,30 @@ sub make_config_cache { &cleanup(); # Run FATES mode for several resolutions and configurations my $clmoptions = "-bgc fates -envxml_dir . -no-megan"; - my @clmres = ( "1x1_brazil", "5x5_amazon", "10x15", "1.9x2.5" ); + my @clmres = ( "1x1_brazil", "5x5_amazon", "4x5", "1.9x2.5" ); foreach my $res ( @clmres ) { - $options = "-res $res"; - my @edoptions = ( "-use_case 2000_control", "", "-namelist \"&a use_lch4=.true.,use_nitrif_denitrif=.true./\"", "-clm_accelerated_spinup on" ); + $options = "-res $res -clm_start_type cold"; + my @edoptions = ( "-use_case 2000_control", + "-use_case 1850_control", + "", + "-namelist \"&a use_lch4=.true.,use_nitrif_denitrif=.true./\"", + "-clm_accelerated_spinup on" + ); foreach my $edop (@edoptions ) { + if ( $res eq "5x5_amazon" && ($edop =~ /1850_control/) ) { + next; + } &make_env_run( ); eval{ system( "$bldnml $options $clmoptions $edop > $tempfile 2>&1 " ); }; is( $@, '', "$options $edop" ); - $cfiles->checkfilesexist( "$options $edop", $mode ); + $cfiles->checkfilesexist( "$options $clmoptions $edop", $mode ); $cfiles->shownmldiff( "default", "standard" ); if ( defined($opts{'compare'}) ) { - $cfiles->doNOTdodiffonfile( "$tempfile", "$options $edop", $mode ); - $cfiles->comparefiles( "$options $edop", $mode, $opts{'compare'} ); + $cfiles->doNOTdodiffonfile( "$tempfile", "$options $clmoptions $edop", $mode ); + $cfiles->comparefiles( "$options $clmoptions $edop", $mode, $opts{'compare'} ); } if ( defined($opts{'generate'}) ) { - $cfiles->copyfiles( "$options $edop", $mode ); + $cfiles->copyfiles( "$options $clmoptions $edop", $mode ); } &cleanup(); } diff --git a/cime_config/config_compsets.xml b/cime_config/config_compsets.xml index de0139fc63..4567e755e7 100644 --- a/cime_config/config_compsets.xml +++ b/cime_config/config_compsets.xml @@ -144,6 +144,11 @@ + + I1850Clm51BgcCrop + 1850_DATM%GSWP3v1_CLM50%BGC-CROP_SICE_SOCN_MOSART_SGLC_SWAV + + I1850Clm51Sp 1850_DATM%GSWP3v1_CLM51%SP_SICE_SOCN_MOSART_SGLC_SWAV @@ -238,6 +243,12 @@ + + I1850Clm51Bgc + 1850_DATM%GSWP3v1_CLM50%BGC_SICE_SOCN_MOSART_SGLC_SWAV + + + I1850Clm51SpNoAnthro 1850_DATM%GSWP3v1_CLM51%SP-NOANTHRO_SICE_SOCN_MOSART_SGLC_SWAV diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index cebb54f7f5..c673d4e544 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1973,12 +1973,13 @@ for ERS test as otherwise it won't work for a sub-day test (no need to run this - + + diff --git a/doc/ChangeLog b/doc/ChangeLog index 4eac6c847f..2ca737aa91 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,133 @@ =============================================================== +Tag name: ctsm5.1.dev021 +Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) +Date: Tue Jan 12 20:21:52 MST 2021 +One-line Summary: Add option for biomass heat storage (BHS) to clm5_1 physics + +Purpose of changes +------------------ + +Add heat stored in biomass (for trees and shrubs) to the surface energy balance calculation. Add +a switch for it and turn it on by default for clm5_1 physics. It's turned off for clm4_5, clm5_0 +physics and when FATES is turned on. Those cases are identical to before, answers only change +when it's turned on. + +Papers describing BHS simulations: +R. Meier, Davin, E., Swenson, S., Lawrence, D., and Schwaab, Jo. (2019). Biomass heat +storage dampens diurnal temperature variations +in forests. Environmental Research Letters. 14. 084026. 10.1088/1748-9326/ab2b4e. + +S.C. Swenson, Burns, S. P., and Lawrence, D. M. ( 2019). The impact of biomass heat storage +on the canopy energy balance and atmospheric stability in the community land model, Journal +of Advances in Modeling Earth Systems, 11, 83– 98. +https://doi.org/10.1029/2018MS001476 + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): + #342 --- set medlynslope to C4 appropriate value for millet and sorghum + (had already been done for miscanthus and switchgrass) + #1246 -- Make spinup_factor_deadwood a variable rather than hardcoded constant + +Known bugs introduced in this tag (include github issue ID): + #1247 -- FATES doesn't work with BHS + +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.] + +[X] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): + New variables are added on the restart files (stem/leaf biomass and stem Temp) + Testing didn't show this as a problem, but theoretically could require + updating initial conditions + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): New compsets, new hist fields + I1850Clm51BgcCrop + I1850Clm51Bgc + History fields: AGSB, AGLB, FSH_STEM, DHDT_CANOPY, RAH1, RAH2, RAW1, RAW2, USTAR, UM, UAF, UM, UAF, + TAF, QAF, OBU, ZETA, VPD, num_iter, RB, TSTEM + DHSDT_CANOPY, TSTEM are default active + +Changes made to namelist defaults (e.g., changed parameter values): New use_biomass_heat_storage + New namelist item: use_biomass_heat_storage to turn on or off (by default only on for clm5_1 + physics for both SP and BGC modes) + +Changes to the datasets (e.g., parameter, surface or initial files): New parameter datasets + All of the params files were updated. New terms were BHS parameters and taper + stocking is now set as "nstem" on the parameter file + taper is now on the parameter file + +Substantial timing or memory changes: Doesn't seem to + Only one short test failed the TPUTCOMP test + Longer tests were not significantly different + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + Setting of leaf/stem biomass should be refactored and removed from CanopyFluxes (as described in #1247) + There's a change in CNGapMortality that could change for DV when AD mode is on + Stability cap (zera) is different for BHS on than off + +Changes to tests or testing: New CLM51 FATES test + +CTSM testing: regular + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - OK (462 comparisons fail because of new namelist and params files) + + python testing (see instructions in python/README.md; document testing done): + + cheyenne -- PASS + + regular tests (aux_clm): + + cheyenne ---- OK + izumi ------- OK + +If the tag used for baseline comparisons was NOT the previous tag, note that here: + + +Answer changes +-------------- + +Changes answers relative to baseline: Yes! + + Summarize any changes to answers, i.e., + - what code configurations: clm5_1 + - what platforms/compilers: all + - nature of change: new climate for tree and shrubs + +Detailed list of changes +------------------------ + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + #1016 -- Heat Storage biomass + +=============================================================== +=============================================================== Tag name: ctsm5.1.dev020 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) Date: Wed Dec 30 00:42:16 MST 2020 diff --git a/doc/ChangeSum b/doc/ChangeSum index 6dc5b4fc0d..0ca1854990 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.1.dev021 erik 01/12/2021 Add option for biomass heat storage (BHS) to clm5_1 physics ctsm5.1.dev020 erik 12/30/2020 Potential roundoff changes in preparation for bio-mass heat storage option ctsm5.1.dev019 sacks 12/19/2020 Fix ndep from coupler ctsm5.1.dev018 slevis 12/08/2020 Add ACTIVE (T/F) column to master hist fields table and alphabetize diff --git a/src/biogeochem/CNDVLightMod.F90 b/src/biogeochem/CNDVLightMod.F90 index 3c498742b9..aa429890ba 100644 --- a/src/biogeochem/CNDVLightMod.F90 +++ b/src/biogeochem/CNDVLightMod.F90 @@ -72,8 +72,9 @@ subroutine Light(bounds, num_natvegp, filter_natvegp, & slatop => pftcon%slatop , & ! Input: specific leaf area at top of canopy, projected area basis (m2/gC) dsladlai => pftcon%dsladlai , & ! Input: dSLA/dLAI, projected area basis (m2/gC) woody => pftcon%woody , & ! Input: woody patch or not - tree => pftcon%tree , & ! Input: tree patch or not - + is_tree => pftcon%is_tree , & ! Input: tree patch or not + is_shrub => pftcon%is_shrub , & ! Input: shrub patch or not + deadstemc => cnveg_carbonstate_inst%deadstemc_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead stem C leafcmax => cnveg_carbonstate_inst%leafcmax_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C storage @@ -132,11 +133,11 @@ subroutine Light(bounds, num_natvegp, filter_natvegp, & fpc_inc(p) = max(0._r8, fpcgrid(p) - fpcgrid_old) if (woody(ivt(p)) == 1._r8) then - if (tree(ivt(p)) == 1) then + if (is_tree(ivt(p))) then numtrees(g) = numtrees(g) + 1 fpc_tree_total(g) = fpc_tree_total(g) + fpcgrid(p) fpc_inc_tree(g) = fpc_inc_tree(g) + fpc_inc(p) - else ! if shrubs + else if (is_shrub(ivt(p))) then fpc_shrub_total(g) = fpc_shrub_total(g) + fpcgrid(p) end if else ! if grass @@ -162,7 +163,7 @@ subroutine Light(bounds, num_natvegp, filter_natvegp, & ! light competition - if (woody(ivt(p))==1._r8 .and. tree(ivt(p))==1._r8) then + if (woody(ivt(p))==1._r8 .and. is_tree(ivt(p))) then if (fpc_tree_total(g) > fpc_tree_max) then @@ -200,7 +201,7 @@ subroutine Light(bounds, num_natvegp, filter_natvegp, & end if - else if (woody(ivt(p))==1._r8 .and. tree(ivt(p))==0._r8) then ! shrub + else if (woody(ivt(p))==1._r8 .and. is_shrub(ivt(p))) then ! shrub if (fpc_shrub_total(g) > fpc_shrub_max(g)) then diff --git a/src/biogeochem/CNDVType.F90 b/src/biogeochem/CNDVType.F90 index 6bc87d8f4b..065e972a15 100644 --- a/src/biogeochem/CNDVType.F90 +++ b/src/biogeochem/CNDVType.F90 @@ -97,7 +97,7 @@ subroutine InitAllocate(this, bounds) use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use clm_varpar , only : maxveg use pftconMod , only : allom1s, allom2s, allom1, allom2, allom3, reinickerp - use pftconMod , only : ntree, nbrdlf_dcd_brl_shrub + use pftconMod , only : nbrdlf_dcd_brl_shrub use pftconMod , only : pftcon ! ! !ARGUMENTS: @@ -148,7 +148,7 @@ subroutine InitAllocate(this, bounds) dgv_ecophyscon%allom2(m) = allom2 dgv_ecophyscon%allom3(m) = allom3 ! modification for shrubs by X.D.Z - if (m > ntree .and. m <= nbrdlf_dcd_brl_shrub ) then + if (pftcon%is_shrub(m)) then dgv_ecophyscon%allom1(m) = allom1s dgv_ecophyscon%allom2(m) = allom2s end if diff --git a/src/biogeochem/CNFireBaseMod.F90 b/src/biogeochem/CNFireBaseMod.F90 index 1ba0efc4d4..b9c125716e 100644 --- a/src/biogeochem/CNFireBaseMod.F90 +++ b/src/biogeochem/CNFireBaseMod.F90 @@ -24,7 +24,7 @@ module CNFireBaseMod use atm2lndType , only : atm2lnd_type use CNDVType , only : dgvs_type use CNVegStateType , only : cnveg_state_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CNVegCarbonFluxType , only : cnveg_carbonflux_type use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type @@ -441,7 +441,7 @@ subroutine CNFireFluxes (this, bounds, num_soilc, filter_soilc, num_soilp, filte ! ! !USES: use clm_time_manager , only: get_step_size_real,get_days_per_year,get_curr_date - use clm_varctl , only: use_cndv, spinup_state + use clm_varctl , only: use_cndv use clm_varcon , only: secspday use pftconMod , only: nc3crop use dynSubgridControlMod , only: run_has_transient_landcover @@ -718,10 +718,7 @@ subroutine CNFireFluxes (this, bounds, num_soilc, filter_soilc, num_soilp, filte ! apply this rate to the patch state variables to get flux rates ! biomass burning ! carbon fluxes - m = 1._r8 - if (spinup_state == 2) then - m = 10._r8 - end if + m = spinup_factor_deadwood m_leafc_to_fire(p) = leafc(p) * f * cc_leaf(patch%itype(p)) m_leafc_storage_to_fire(p) = leafc_storage(p) * f * cc_other(patch%itype(p)) diff --git a/src/biogeochem/CNFireLi2014Mod.F90 b/src/biogeochem/CNFireLi2014Mod.F90 index a69efbfeae..e8fd78230e 100644 --- a/src/biogeochem/CNFireLi2014Mod.F90 +++ b/src/biogeochem/CNFireLi2014Mod.F90 @@ -17,7 +17,7 @@ module CNFireLi2014Mod use shr_kind_mod , only : r8 => shr_kind_r8, CL => shr_kind_CL use shr_const_mod , only : SHR_CONST_PI,SHR_CONST_TKFRZ use shr_infnan_mod , only : shr_infnan_isnan - use clm_varctl , only : iulog, spinup_state + use clm_varctl , only : iulog use clm_varpar , only : nlevdecomp, ndecomp_pools, nlevdecomp_full use clm_varcon , only : dzsoi_decomp use pftconMod , only : noveg, pftcon @@ -27,7 +27,7 @@ module CNFireLi2014Mod use atm2lndType , only : atm2lnd_type use CNDVType , only : dgvs_type use CNVegStateType , only : cnveg_state_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CNVegCarbonFluxType , only : cnveg_carbonflux_type use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type @@ -920,10 +920,7 @@ subroutine CNFireFluxes (this, bounds, num_soilc, filter_soilc, num_soilp, filte ! apply this rate to the patch state variables to get flux rates ! biomass burning ! carbon fluxes - m = 1._r8 - if (spinup_state == 2) then - m = 10._r8 - end if + m = spinup_factor_deadwood m_leafc_to_fire(p) = leafc(p) * f * cc_leaf(patch%itype(p)) m_leafc_storage_to_fire(p) = leafc_storage(p) * f * cc_other(patch%itype(p)) diff --git a/src/biogeochem/CNFireLi2016Mod.F90 b/src/biogeochem/CNFireLi2016Mod.F90 index b1cdafb4ee..afd661cd28 100644 --- a/src/biogeochem/CNFireLi2016Mod.F90 +++ b/src/biogeochem/CNFireLi2016Mod.F90 @@ -27,7 +27,7 @@ module CNFireLi2016Mod use atm2lndType , only : atm2lnd_type use CNDVType , only : dgvs_type use CNVegStateType , only : cnveg_state_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CNVegCarbonFluxType , only : cnveg_carbonflux_type use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type @@ -420,19 +420,11 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ end if end if end if - if (spinup_state == 2) then - rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & - frootc_xfer(p) + deadcrootc(p) * 10._r8 + & + rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & + frootc_xfer(p) + deadcrootc(p) * spinup_factor_deadwood + & deadcrootc_storage(p) + deadcrootc_xfer(p) + & livecrootc(p)+livecrootc_storage(p) + & livecrootc_xfer(p))*patch%wtcol(p) - else - rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & - frootc_xfer(p) + deadcrootc(p) + & - deadcrootc_storage(p) + deadcrootc_xfer(p) + & - livecrootc(p)+livecrootc_storage(p) + & - livecrootc_xfer(p))*patch%wtcol(p) - endif fsr_col(c) = fsr_col(c) + fsr_pft(patch%itype(p))*patch%wtcol(p)/(1.0_r8-cropf_col(c)) @@ -588,7 +580,7 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ if( cropf_col(c) < 1._r8 )then fuelc(c) = totlitc(c)+totvegc(c)-rootc_col(c)-fuelc_crop(c)*cropf_col(c) if (spinup_state == 2) then - fuelc(c) = fuelc(c) + ((10._r8 - 1._r8)*deadstemc_col(c)) + fuelc(c) = fuelc(c) + ((spinup_factor_deadwood - 1._r8)*deadstemc_col(c)) do j = 1, nlevdecomp fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j) * spinup_factor(i_cwd) & * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) diff --git a/src/biogeochem/CNFireLi2021Mod.F90 b/src/biogeochem/CNFireLi2021Mod.F90 index 77b0693fca..973bbf46ed 100644 --- a/src/biogeochem/CNFireLi2021Mod.F90 +++ b/src/biogeochem/CNFireLi2021Mod.F90 @@ -27,7 +27,7 @@ module CNFireLi2021Mod use atm2lndType , only : atm2lnd_type use CNDVType , only : dgvs_type use CNVegStateType , only : cnveg_state_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CNVegCarbonFluxType , only : cnveg_carbonflux_type use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type @@ -421,19 +421,11 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ end if end if end if - if (spinup_state == 2) then - rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & - frootc_xfer(p) + deadcrootc(p) * 10._r8 + & + rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & + frootc_xfer(p) + deadcrootc(p) * spinup_factor_deadwood + & deadcrootc_storage(p) + deadcrootc_xfer(p) + & livecrootc(p)+livecrootc_storage(p) + & livecrootc_xfer(p))*patch%wtcol(p) - else - rootc_col(c) = rootc_col(c) + (frootc(p) + frootc_storage(p) + & - frootc_xfer(p) + deadcrootc(p) + & - deadcrootc_storage(p) + deadcrootc_xfer(p) + & - livecrootc(p)+livecrootc_storage(p) + & - livecrootc_xfer(p))*patch%wtcol(p) - endif fsr_col(c) = fsr_col(c) + fsr_pft(patch%itype(p))*patch%wtcol(p)/(1.0_r8-cropf_col(c)) @@ -589,7 +581,7 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ if( cropf_col(c) < 1._r8 )then fuelc(c) = totlitc(c)+totvegc(c)-rootc_col(c)-fuelc_crop(c)*cropf_col(c) if (spinup_state == 2) then - fuelc(c) = fuelc(c) + ((10._r8 - 1._r8)*deadstemc_col(c)) + fuelc(c) = fuelc(c) + ((spinup_factor_deadwood - 1._r8)*deadstemc_col(c)) do j = 1, nlevdecomp fuelc(c) = fuelc(c)+decomp_cpools_vr(c,j,i_cwd) * dzsoi_decomp(j) * spinup_factor(i_cwd) & * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) diff --git a/src/biogeochem/CNGapMortalityMod.F90 b/src/biogeochem/CNGapMortalityMod.F90 index e313cdf213..cd02221de4 100644 --- a/src/biogeochem/CNGapMortalityMod.F90 +++ b/src/biogeochem/CNGapMortalityMod.F90 @@ -14,7 +14,7 @@ module CNGapMortalityMod use shr_log_mod , only : errMsg => shr_log_errMsg use pftconMod , only : pftcon use CNDVType , only : dgvs_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CNVegCarbonFluxType , only : cnveg_carbonflux_type use CNVegNitrogenStateType , only : cnveg_nitrogenstate_type use CNVegNitrogenFluxType , only : cnveg_nitrogenflux_type @@ -191,13 +191,8 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so cnveg_carbonflux_inst%m_frootc_to_litter_patch(p) = cnveg_carbonstate_inst%frootc_patch(p) * m cnveg_carbonflux_inst%m_livestemc_to_litter_patch(p) = cnveg_carbonstate_inst%livestemc_patch(p) * m cnveg_carbonflux_inst%m_livecrootc_to_litter_patch(p) = cnveg_carbonstate_inst%livecrootc_patch(p) * m - if (spinup_state == 2 .and. .not. use_cndv) then !accelerate mortality of dead woody pools - cnveg_carbonflux_inst%m_deadstemc_to_litter_patch(p) = cnveg_carbonstate_inst%deadstemc_patch(p) * m * 10._r8 - cnveg_carbonflux_inst%m_deadcrootc_to_litter_patch(p) = cnveg_carbonstate_inst%deadcrootc_patch(p) * m * 10._r8 - else - cnveg_carbonflux_inst%m_deadstemc_to_litter_patch(p) = cnveg_carbonstate_inst%deadstemc_patch(p) * m - cnveg_carbonflux_inst%m_deadcrootc_to_litter_patch(p) = cnveg_carbonstate_inst%deadcrootc_patch(p) * m - end if + cnveg_carbonflux_inst%m_deadstemc_to_litter_patch(p) = cnveg_carbonstate_inst%deadstemc_patch(p) * m * spinup_factor_deadwood + cnveg_carbonflux_inst%m_deadcrootc_to_litter_patch(p) = cnveg_carbonstate_inst%deadcrootc_patch(p) * m * spinup_factor_deadwood ! storage pools cnveg_carbonflux_inst%m_leafc_storage_to_litter_patch(p) = cnveg_carbonstate_inst%leafc_storage_patch(p) * m @@ -228,8 +223,8 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so cnveg_nitrogenflux_inst%m_livecrootn_to_litter_patch(p) = cnveg_nitrogenstate_inst%livecrootn_patch(p) * m if (spinup_state == 2 .and. .not. use_cndv) then !accelerate mortality of dead woody pools - cnveg_nitrogenflux_inst%m_deadstemn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadstemn_patch(p) * m * 10._r8 - cnveg_nitrogenflux_inst%m_deadcrootn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadcrootn_patch(p) * m * 10._r8 + cnveg_nitrogenflux_inst%m_deadstemn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadstemn_patch(p) * m * spinup_factor_deadwood + cnveg_nitrogenflux_inst%m_deadcrootn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadcrootn_patch(p) * m * spinup_factor_deadwood else cnveg_nitrogenflux_inst%m_deadstemn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadstemn_patch(p) * m cnveg_nitrogenflux_inst%m_deadcrootn_to_litter_patch(p) = cnveg_nitrogenstate_inst%deadcrootn_patch(p) * m diff --git a/src/biogeochem/CNVegCarbonStateType.F90 b/src/biogeochem/CNVegCarbonStateType.F90 index 0a8ea2e03c..f7c9178453 100644 --- a/src/biogeochem/CNVegCarbonStateType.F90 +++ b/src/biogeochem/CNVegCarbonStateType.F90 @@ -102,6 +102,9 @@ module CNVegCarbonStateType end type cnveg_carbonstate_type + real(r8), public :: spinup_factor_deadwood = 1.0_r8 ! Spinup factor used for this simulation + real(r8), public :: spinup_factor_AD = 10.0_r8 ! Spinup factor used when in Accelerated Decomposition mode + ! !PRIVATE DATA: type, private :: cnvegcarbonstate_const_type @@ -1036,7 +1039,7 @@ end subroutine InitCold !----------------------------------------------------------------------- subroutine Restart ( this, bounds, ncid, flag, carbon_type, reseed_dead_plants, & c12_cnveg_carbonstate_inst, filter_reseed_patch, & - num_reseed_patch) + num_reseed_patch, spinup_factor4deadwood ) ! ! !DESCRIPTION: ! Read/write CN restart data for carbon state @@ -1062,6 +1065,7 @@ subroutine Restart ( this, bounds, ncid, flag, carbon_type, reseed_dead_plants, type (cnveg_carbonstate_type) , intent(in), optional :: c12_cnveg_carbonstate_inst integer , intent(out), optional :: filter_reseed_patch(:) integer , intent(out), optional :: num_reseed_patch + real(r8) , intent(out), optional :: spinup_factor4deadwood ! ! !LOCAL VARIABLES: integer :: i,j,k,l,c,p @@ -1243,22 +1247,23 @@ subroutine Restart ( this, bounds, ncid, flag, carbon_type, reseed_dead_plants, if (flag == 'read' .and. spinup_state /= restart_file_spinup_state .and. .not. use_cndv) then if ( masterproc ) write(iulog, *) 'exit_spinup ',exit_spinup,' restart_file_spinup_state ',restart_file_spinup_state + if ( spinup_state == 2 ) spinup_factor_deadwood = spinup_factor_AD if (spinup_state <= 1 .and. restart_file_spinup_state == 2 ) then if ( masterproc ) write(iulog,*) ' CNRest: taking Dead wood C pools out of AD spinup mode' exit_spinup = .true. - if ( masterproc ) write(iulog, *) 'Multiplying stemc and crootc by 10 for exit spinup' + if ( masterproc ) write(iulog, *) 'Multiplying stemc and crootc by ', spinup_factor_AD, ' for exit spinup' do i = bounds%begp,bounds%endp - this%deadstemc_patch(i) = this%deadstemc_patch(i) * 10._r8 - this%deadcrootc_patch(i) = this%deadcrootc_patch(i) * 10._r8 + this%deadstemc_patch(i) = this%deadstemc_patch(i) * spinup_factor_AD + this%deadcrootc_patch(i) = this%deadcrootc_patch(i) * spinup_factor_AD end do else if (spinup_state == 2 .and. restart_file_spinup_state <= 1 )then if (spinup_state == 2 .and. restart_file_spinup_state <= 1 )then if ( masterproc ) write(iulog,*) ' CNRest: taking Dead wood C pools into AD spinup mode' enter_spinup = .true. - if ( masterproc ) write(iulog, *) 'Dividing stemc and crootc by 10 for enter spinup ' + if ( masterproc ) write(iulog, *) 'Dividing stemc and crootc by ', spinup_factor_AD, 'for enter spinup ' do i = bounds%begp,bounds%endp - this%deadstemc_patch(i) = this%deadstemc_patch(i) / 10._r8 - this%deadcrootc_patch(i) = this%deadcrootc_patch(i) / 10._r8 + this%deadstemc_patch(i) = this%deadstemc_patch(i) / spinup_factor_AD + this%deadcrootc_patch(i) = this%deadcrootc_patch(i) / spinup_factor_AD end do end if end if @@ -2281,6 +2286,9 @@ subroutine Restart ( this, bounds, ncid, flag, carbon_type, reseed_dead_plants, end if end if + ! Output spinup factor for deadwood (dead stem and dead course root) + if ( present(spinup_factor4deadwood) ) spinup_factor4deadwood = spinup_factor_AD + end subroutine Restart !----------------------------------------------------------------------- diff --git a/src/biogeochem/CNVegNitrogenStateType.F90 b/src/biogeochem/CNVegNitrogenStateType.F90 index 1b06cd3fc3..10a191ded9 100644 --- a/src/biogeochem/CNVegNitrogenStateType.F90 +++ b/src/biogeochem/CNVegNitrogenStateType.F90 @@ -529,7 +529,8 @@ end subroutine InitCold !----------------------------------------------------------------------- subroutine Restart ( this, bounds, ncid, flag, leafc_patch, & leafc_storage_patch, frootc_patch, frootc_storage_patch, & - deadstemc_patch, filter_reseed_patch, num_reseed_patch ) + deadstemc_patch, filter_reseed_patch, num_reseed_patch, & + spinup_factor_deadwood ) ! ! !DESCRIPTION: ! Read/write restart data @@ -554,6 +555,7 @@ subroutine Restart ( this, bounds, ncid, flag, leafc_patch, & real(r8) , intent(in) :: deadstemc_patch(bounds%begp:) integer , intent(in) :: filter_reseed_patch(:) integer , intent(in) :: num_reseed_patch + real(r8) , intent(in) :: spinup_factor_deadwood ! ! !LOCAL VARIABLES: integer :: i, p, l @@ -719,18 +721,18 @@ subroutine Restart ( this, bounds, ncid, flag, leafc_patch, & if (spinup_state <= 1 .and. restart_file_spinup_state == 2 ) then if ( masterproc ) write(iulog,*) ' CNRest: taking Dead wood N pools out of AD spinup mode' exit_spinup = .true. - if ( masterproc ) write(iulog, *) 'Multiplying stemn and crootn by 10 for exit spinup ' + if ( masterproc ) write(iulog, *) 'Multiplying stemn and crootn by ', spinup_factor_deadwood, 'for exit spinup ' do i = bounds%begp,bounds%endp - this%deadstemn_patch(i) = this%deadstemn_patch(i) * 10._r8 - this%deadcrootn_patch(i) = this%deadcrootn_patch(i) * 10._r8 + this%deadstemn_patch(i) = this%deadstemn_patch(i) * spinup_factor_deadwood + this%deadcrootn_patch(i) = this%deadcrootn_patch(i) * spinup_factor_deadwood end do else if (spinup_state == 2 .and. restart_file_spinup_state <= 1 ) then if ( masterproc ) write(iulog,*) ' CNRest: taking Dead wood N pools into AD spinup mode' enter_spinup = .true. - if ( masterproc ) write(iulog, *) 'Dividing stemn and crootn by 10 for enter spinup ' + if ( masterproc ) write(iulog, *) 'Dividing stemn and crootn by ', spinup_factor_deadwood, 'for enter spinup ' do i = bounds%begp,bounds%endp - this%deadstemn_patch(i) = this%deadstemn_patch(i) / 10._r8 - this%deadcrootn_patch(i) = this%deadcrootn_patch(i) / 10._r8 + this%deadstemn_patch(i) = this%deadstemn_patch(i) / spinup_factor_deadwood + this%deadcrootn_patch(i) = this%deadcrootn_patch(i) / spinup_factor_deadwood end do endif diff --git a/src/biogeochem/CNVegStructUpdateMod.F90 b/src/biogeochem/CNVegStructUpdateMod.F90 index 22bbea8614..1c0ec70fa0 100644 --- a/src/biogeochem/CNVegStructUpdateMod.F90 +++ b/src/biogeochem/CNVegStructUpdateMod.F90 @@ -13,9 +13,10 @@ module CNVegStructUpdateMod use CNDVType , only : dgvs_type use CNVegStateType , only : cnveg_state_type use CropType , only : crop_type - use CNVegCarbonStateType , only : cnveg_carbonstate_type + use CNVegCarbonStateType , only : cnveg_carbonstate_type, spinup_factor_deadwood use CanopyStateType , only : canopystate_type use PatchType , only : patch + use decompMod , only : bounds_type ! implicit none private @@ -27,7 +28,7 @@ module CNVegStructUpdateMod contains !----------------------------------------------------------------------- - subroutine CNVegStructUpdate(num_soilp, filter_soilp, & + subroutine CNVegStructUpdate(bounds,num_soilp, filter_soilp, & waterdiagnosticbulk_inst, frictionvel_inst, dgvs_inst, cnveg_state_inst, crop_inst, & cnveg_carbonstate_inst, canopystate_inst) ! @@ -44,10 +45,12 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & use pftconMod , only : nmiscanthus, nirrig_miscanthus, nswitchgrass, nirrig_switchgrass use pftconMod , only : pftcon - use clm_varctl , only : spinup_state + use clm_varctl , only : spinup_state, use_biomass_heat_storage + use clm_varcon , only : c_to_b use clm_time_manager , only : get_rad_step_size ! ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_soilp ! number of column soil points in patch filter integer , intent(in) :: filter_soilp(:) ! patch filter for soil points type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst @@ -65,8 +68,7 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & ! !LOCAL VARIABLES: integer :: p,c,g ! indices integer :: fp ! lake filter indices - real(r8) :: taper ! ratio of height:radius_breast_height (tree allometry) - real(r8) :: stocking ! #stems / ha (stocking density) + real(r8) :: stocking ! #stems / ha (stocking density) real(r8) :: ol ! thickness of canopy layer covered by snow (m) real(r8) :: fb ! fraction of canopy layer covered by snow real(r8) :: tlai_old ! for use in Zeng tsai formula @@ -101,7 +103,10 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & dwood => pftcon%dwood , & ! Input: density of wood (gC/m^3) ztopmx => pftcon%ztopmx , & ! Input: laimx => pftcon%laimx , & ! Input: - + nstem => pftcon%nstem , & ! Input: Tree number density (#ind/m2) + taper => pftcon%taper , & ! Input: ratio of height:radius_breast_height (tree allometry) + fbw => pftcon%fbw , & ! Input: Fraction of fresh biomass that is water + allom2 => dgv_ecophyscon%allom2 , & ! Input: [real(r8) (:) ] ecophys const allom3 => dgv_ecophyscon%allom3 , & ! Input: [real(r8) (:) ] ecophys const @@ -115,6 +120,7 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & leafc => cnveg_carbonstate_inst%leafc_patch , & ! Input: [real(r8) (:) ] (gC/m2) leaf C deadstemc => cnveg_carbonstate_inst%deadstemc_patch , & ! Input: [real(r8) (:) ] (gC/m2) dead stem C + livestemc => cnveg_carbonstate_inst%livestemc_patch , & ! Input: [real(r8) (:) ] (gC/m2) live stem C farea_burned => cnveg_state_inst%farea_burned_col , & ! Input: [real(r8) (:) ] F. Li and S. Levis htmx => cnveg_state_inst%htmx_patch , & ! Output: [real(r8) (:) ] max hgt attained by a crop during yr (m) @@ -125,6 +131,8 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & ! *** Key Output from CN*** tlai => canopystate_inst%tlai_patch , & ! Output: [real(r8) (:) ] one-sided leaf area index, no burying by snow tsai => canopystate_inst%tsai_patch , & ! Output: [real(r8) (:) ] one-sided stem area index, no burying by snow + stem_biomass => canopystate_inst%stem_biomass_patch , & ! Output: [real(r8) (:) ] Aboveground stem biomass (kg/m**2) + leaf_biomass => canopystate_inst%leaf_biomass_patch , & ! Output: [real(r8) (:) ] Aboveground leave biomass (kg/m**2) htop => canopystate_inst%htop_patch , & ! Output: [real(r8) (:) ] canopy top (m) hbot => canopystate_inst%hbot_patch , & ! Output: [real(r8) (:) ] canopy bottom (m) elai => canopystate_inst%elai_patch , & ! Output: [real(r8) (:) ] one-sided leaf area index with burying by snow @@ -134,13 +142,6 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & dt = real( get_rad_step_size(), r8 ) - ! constant allometric parameters - taper = 200._r8 - stocking = 1000._r8 - - ! convert from stems/ha -> stems/m^2 - stocking = stocking / 10000._r8 - ! patch loop do fp = 1,num_soilp p = filter_soilp(fp) @@ -182,25 +183,15 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & if (woody(ivt(p)) == 1._r8) then - ! trees and shrubs - - ! if shrubs have a squat taper - if (ivt(p) >= nbrdlf_evr_shrub .and. ivt(p) <= nbrdlf_dcd_brl_shrub) then - taper = 10._r8 - ! otherwise have a tall taper - else - taper = 200._r8 - end if - ! trees and shrubs for now have a very simple allometry, with hard-wired - ! stem taper (height:radius) and hard-wired stocking density (#individuals/area) + ! stem taper (height:radius) and nstem from PFT parameter file if (use_cndv) then if (fpcgrid(p) > 0._r8 .and. nind(p) > 0._r8) then stocking = nind(p)/fpcgrid(p) !#ind/m2 nat veg area -> #ind/m2 patch area htop(p) = allom2(ivt(p)) * ( (24._r8 * deadstemc(p) / & - (SHR_CONST_PI * stocking * dwood(ivt(p)) * taper))**(1._r8/3._r8) )**allom3(ivt(p)) ! lpj's htop w/ cn's stemdiam + (SHR_CONST_PI * stocking * dwood(ivt(p)) * taper(ivt(p))))**(1._r8/3._r8) )**allom3(ivt(p)) ! lpj's htop w/ cn's stemdiam else htop(p) = 0._r8 @@ -208,16 +199,27 @@ subroutine CNVegStructUpdate(num_soilp, filter_soilp, & else !correct height calculation if doing accelerated spinup - if (spinup_state == 2) then - htop(p) = ((3._r8 * deadstemc(p) * 10._r8 * taper * taper)/ & - (SHR_CONST_PI * stocking * dwood(ivt(p))))**(1._r8/3._r8) - else - htop(p) = ((3._r8 * deadstemc(p) * taper * taper)/ & - (SHR_CONST_PI * stocking * dwood(ivt(p))))**(1._r8/3._r8) - end if + htop(p) = ((3._r8 * deadstemc(p) * spinup_factor_deadwood * taper(ivt(p)) * taper(ivt(p)))/ & + (SHR_CONST_PI * nstem(ivt(p)) * dwood(ivt(p))))**(1._r8/3._r8) endif + ! + ! calculate vegetation physiological parameters used in biomass heat storage + ! + if (use_biomass_heat_storage) then + ! Assumes fbw (fraction of biomass that is water) is the same for leaves and stems + leaf_biomass(p) = max(0.0025_r8,leafc(p)) & + * c_to_b * 1.e-3_r8 / (1._r8 - fbw(ivt(p))) + + stem_biomass(p) = (spinup_factor_deadwood*deadstemc(p) + livestemc(p)) & + * c_to_b * 1.e-3_r8 / (1._r8 - fbw(ivt(p))) + + else + leaf_biomass(p) = 0_r8 + stem_biomass(p) = 0_r8 + end if + ! Peter Thornton, 5/3/2004 ! Adding test to keep htop from getting too close to forcing height for windspeed ! Also added for grass, below, although it is not likely to ever be an issue. diff --git a/src/biogeochem/CNVegetationFacade.F90 b/src/biogeochem/CNVegetationFacade.F90 index 44e904bf2c..9ef32b4563 100644 --- a/src/biogeochem/CNVegetationFacade.F90 +++ b/src/biogeochem/CNVegetationFacade.F90 @@ -451,6 +451,7 @@ subroutine Restart(this, bounds, ncid, flag) ! !LOCAL VARIABLES: integer :: begp, endp + real(r8) :: spinup_factor4deadwood ! Spinup factor used for deadwood (dead-stem and dead course root) character(len=*), parameter :: subname = 'Restart' !----------------------------------------------------------------------- @@ -460,10 +461,13 @@ subroutine Restart(this, bounds, ncid, flag) endp = bounds%endp call this%cnveg_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c12', & reseed_dead_plants=this%reseed_dead_plants, filter_reseed_patch=reseed_patch, & - num_reseed_patch=num_reseed_patch ) + num_reseed_patch=num_reseed_patch, spinup_factor4deadwood=spinup_factor4deadwood ) if ( flag /= 'read' .and. num_reseed_patch /= 0 )then call endrun(msg="ERROR num_reseed should be zero and is not"//errmsg(sourcefile, __LINE__)) end if + if ( flag /= 'read' .and. spinup_factor4deadwood /= 10_r8 )then + call endrun(msg="ERROR spinup_factor4deadwood should be 10 and is not"//errmsg(sourcefile, __LINE__)) + end if if (use_c13) then call this%c13_cnveg_carbonstate_inst%restart(bounds, ncid, flag=flag, carbon_type='c13', & reseed_dead_plants=this%reseed_dead_plants, c12_cnveg_carbonstate_inst=this%cnveg_carbonstate_inst) @@ -487,7 +491,8 @@ subroutine Restart(this, bounds, ncid, flag) frootc_patch=this%cnveg_carbonstate_inst%frootc_patch(begp:endp), & frootc_storage_patch=this%cnveg_carbonstate_inst%frootc_storage_patch(begp:endp), & deadstemc_patch=this%cnveg_carbonstate_inst%deadstemc_patch(begp:endp), & - filter_reseed_patch=reseed_patch, num_reseed_patch=num_reseed_patch) + filter_reseed_patch=reseed_patch, num_reseed_patch=num_reseed_patch, & + spinup_factor_deadwood=spinup_factor4deadwood ) call this%cnveg_nitrogenflux_inst%restart(bounds, ncid, flag=flag) call this%cnveg_state_inst%restart(bounds, ncid, flag=flag, & cnveg_carbonstate=this%cnveg_carbonstate_inst, & @@ -1084,7 +1089,7 @@ subroutine EcosystemDynamicsPostDrainage(this, bounds, num_allc, filter_allc, & ! vegetation structure (LAI, SAI, height) if (doalb) then - call CNVegStructUpdate(num_soilp, filter_soilp, & + call CNVegStructUpdate(bounds,num_soilp, filter_soilp, & waterdiagnosticbulk_inst, frictionvel_inst, this%dgvs_inst, this%cnveg_state_inst, & crop_inst, this%cnveg_carbonstate_inst, canopystate_inst) end if diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 69f7cdb167..279e6cc13a 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -340,6 +340,7 @@ subroutine BalanceCheck( bounds, & qflx_sfc_irrig => waterflux_inst%qflx_sfc_irrig_col , & ! Input: [real(r8) (:) ] irrigation flux (mm H2O /s) qflx_glcice_dyn_water_flux => waterflux_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system) + dhsdt_canopy => energyflux_inst%dhsdt_canopy_patch , & ! Input: [real(r8) (:) ] change in heat content of canopy (W/m**2) [+ to atm] eflx_lwrad_out => energyflux_inst%eflx_lwrad_out_patch , & ! Input: [real(r8) (:) ] emitted infrared (longwave) radiation (W/m**2) eflx_lwrad_net => energyflux_inst%eflx_lwrad_net_patch , & ! Input: [real(r8) (:) ] net infrared (longwave) rad (W/m**2) [+ = to atm] eflx_sh_tot => energyflux_inst%eflx_sh_tot_patch , & ! Input: [real(r8) (:) ] total sensible heat flux (W/m**2) [+ to atm] @@ -616,7 +617,7 @@ subroutine BalanceCheck( bounds, & if (.not. lun%urbpoi(l)) then errseb(p) = sabv(p) + sabg_chk(p) + forc_lwrad(c) - eflx_lwrad_out(p) & - - eflx_sh_tot(p) - eflx_lh_tot(p) - eflx_soil_grnd(p) + - eflx_sh_tot(p) - eflx_lh_tot(p) - eflx_soil_grnd(p) - dhsdt_canopy(p) else errseb(p) = sabv(p) + sabg(p) & - eflx_lwrad_net(p) & @@ -702,6 +703,7 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'eflx_sh_tot = ' ,eflx_sh_tot(indexp) write(iulog,*)'eflx_lh_tot = ' ,eflx_lh_tot(indexp) write(iulog,*)'eflx_soil_grnd = ' ,eflx_soil_grnd(indexp) + write(iulog,*)'dhsdt_canopy = ' ,dhsdt_canopy(indexp) write(iulog,*)'fsa fsr = ' ,fsa(indexp), fsr(indexp) write(iulog,*)'fabd fabi = ' ,fabd(indexp,:), fabi(indexp,:) write(iulog,*)'albd albi = ' ,albd(indexp,:), albi(indexp,:) diff --git a/src/biogeophys/BareGroundFluxesMod.F90 b/src/biogeophys/BareGroundFluxesMod.F90 index 6751459246..89893e3aa4 100644 --- a/src/biogeophys/BareGroundFluxesMod.F90 +++ b/src/biogeophys/BareGroundFluxesMod.F90 @@ -145,7 +145,9 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & real(r8) :: www ! surface soil wetness [-] !------------------------------------------------------------------------------ - associate( & + associate( & + dhsdt_canopy => energyflux_inst%dhsdt_canopy_patch , & ! Output: [real(r8) (:) ] change in heat storage of stem (W/m**2) [+ to atm] + eflx_sh_stem => energyflux_inst%eflx_sh_stem_patch , & ! Output: [real(r8) (:) ] sensible heat flux from stems (W/m**2) [+ to atm] soilresis => soilstate_inst%soilresis_col , & ! Input: [real(r8) (:,:) ] evaporative soil resistance (s/m) snl => col%snl , & ! Input: [integer (:) ] number of snow layers dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) @@ -236,19 +238,19 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & forc_hgt_u_patch => frictionvel_inst%forc_hgt_u_patch , & ! Input: u10_clm => frictionvel_inst%u10_clm_patch , & ! Input: [real(r8) (:) ] 10 m height winds (m/s) zetamax => frictionvel_inst%zetamaxstable , & ! Input: [real(r8) ] max zeta value under stable conditions - z0mg_col => frictionvel_inst%z0mg_col , & ! Output: [real(r8) (:) ] roughness length, momentum [m] - z0hg_col => frictionvel_inst%z0hg_col , & ! Output: [real(r8) (:) ] roughness length, sensible heat [m] - z0qg_col => frictionvel_inst%z0qg_col , & ! Output: [real(r8) (:) ] roughness length, latent heat [m] - ram1 => frictionvel_inst%ram1_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance (s/m) - + z0mg_col => frictionvel_inst%z0mg_col , & ! Output: [real(r8) (:) ] roughness length, momentum [m] + z0hg_col => frictionvel_inst%z0hg_col , & ! Output: [real(r8) (:) ] roughness length, sensible heat [m] + z0qg_col => frictionvel_inst%z0qg_col , & ! Output: [real(r8) (:) ] roughness length, latent heat [m] + ram1 => frictionvel_inst%ram1_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance (s/m) + num_iter => frictionvel_inst%num_iter_patch , & ! Output: [real(r8) (:) ] number of iterations htvp => energyflux_inst%htvp_col , & ! Input: [real(r8) (:) ] latent heat of evaporation (/sublimation) [J/kg] - qflx_ev_snow => waterfluxbulk_inst%qflx_ev_snow_patch , & ! Output: [real(r8) (:) ] evaporation flux from snow (mm H2O/s) [+ to atm] - qflx_ev_soil => waterfluxbulk_inst%qflx_ev_soil_patch , & ! Output: [real(r8) (:) ] evaporation flux from soil (mm H2O/s) [+ to atm] - qflx_ev_h2osfc => waterfluxbulk_inst%qflx_ev_h2osfc_patch , & ! Output: [real(r8) (:) ] evaporation flux from h2osfc (mm H2O/s) [+ to atm] - qflx_evap_soi => waterfluxbulk_inst%qflx_evap_soi_patch , & ! Output: [real(r8) (:) ] soil evaporation (mm H2O/s) (+ = to atm) - qflx_evap_tot => waterfluxbulk_inst%qflx_evap_tot_patch , & ! Output: [real(r8) (:) ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg - qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) - qflx_evap_veg => waterfluxbulk_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm) + qflx_ev_snow => waterfluxbulk_inst%qflx_ev_snow_patch , & ! Output: [real(r8) (:) ] evaporation flux from snow (mm H2O/s) [+ to atm] + qflx_ev_soil => waterfluxbulk_inst%qflx_ev_soil_patch , & ! Output: [real(r8) (:) ] evaporation flux from soil (mm H2O/s) [+ to atm] + qflx_ev_h2osfc => waterfluxbulk_inst%qflx_ev_h2osfc_patch , & ! Output: [real(r8) (:) ] evaporation flux from h2osfc (mm H2O/s) [+ to atm] + qflx_evap_soi => waterfluxbulk_inst%qflx_evap_soi_patch , & ! Output: [real(r8) (:) ] soil evaporation (mm H2O/s) (+ = to atm) + qflx_evap_tot => waterfluxbulk_inst%qflx_evap_tot_patch , & ! Output: [real(r8) (:) ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg + qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) + qflx_evap_veg => waterfluxbulk_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm) rssun => photosyns_inst%rssun_patch , & ! Output: [real(r8) (:) ] leaf sunlit stomatal resistance (s/m) (output from Photosynthesis) rssha => photosyns_inst%rssha_patch , & ! Output: [real(r8) (:) ] leaf shaded stomatal resistance (s/m) (output from Photosynthesis) @@ -287,6 +289,8 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & displa(p) = 0._r8 dlrad(p) = 0._r8 ulrad(p) = 0._r8 + dhsdt_canopy(p) = 0._r8 + eflx_sh_stem(p) = 0._r8 ur(p) = max(params_inst%wind_min,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g))) dth(p) = thm(p)-t_grnd(c) @@ -304,6 +308,8 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & call frictionvel_inst%MoninObukIni(ur(p), thv(c), dthv, zldis(p), z0mg_patch(p), um(p), obu(p)) + ! Initialize iteration counter + num_iter(p) = 0 end do ! Perform stability iteration @@ -338,8 +344,9 @@ subroutine BareGroundFluxes(bounds, num_noexposedvegp, filter_noexposedvegp, & um(p) = sqrt(ur(p)*ur(p) + wc*wc) end if obu(p) = zldis(p)/zeta - end do + num_iter(p) = iter + end do end do ! end stability iteration do f = 1, num_noexposedvegp diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index b96b43187c..348af3a613 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -14,8 +14,8 @@ module CanopyFluxesMod use shr_log_mod , only : errMsg => shr_log_errMsg use abortutils , only : endrun use clm_varctl , only : iulog, use_cn, use_lch4, use_c13, use_c14, use_cndv, use_fates, & - use_luna, use_hydrstress - use clm_varpar , only : nlevgrnd, nlevsno + use_luna, use_hydrstress, use_biomass_heat_storage + use clm_varpar , only : nlevgrnd, nlevsno, mxpft use clm_varcon , only : namep use pftconMod , only : pftcon use decompMod , only : bounds_type @@ -67,6 +67,7 @@ module CanopyFluxesMod real(r8) :: wind_min ! Minimum wind speed at the atmospheric forcing height (m/s) end type params_type type(params_type), private :: params_inst + ! ! !PUBLIC DATA MEMBERS: ! true => btran is based only on unfrozen soil levels @@ -77,7 +78,7 @@ module CanopyFluxesMod logical, public :: perchroot_alt = .false. ! ! !PRIVATE DATA MEMBERS: - logical, private :: use_undercanopy_stability = .true. ! use undercanopy stability term or not + logical, private :: use_undercanopy_stability = .false. ! use undercanopy stability term or not integer, private :: itmax_canopy_fluxes = -1 ! max # of iterations used in subroutine CanopyFluxes character(len=*), parameter, private :: sourcefile = & @@ -111,6 +112,7 @@ subroutine CanopyFluxesReadNML(NLFilename) !----------------------------------------------------------------------- namelist /canopyfluxes_inparm/ use_undercanopy_stability + namelist /canopyfluxes_inparm/ use_biomass_heat_storage namelist /canopyfluxes_inparm/ itmax_canopy_fluxes ! Initialize options to default values, in case they are not specified in @@ -139,6 +141,7 @@ subroutine CanopyFluxesReadNML(NLFilename) end if call shr_mpi_bcast (use_undercanopy_stability, mpicom) + call shr_mpi_bcast (use_biomass_heat_storage, mpicom) call shr_mpi_bcast (itmax_canopy_fluxes, mpicom) if (masterproc) then @@ -220,11 +223,12 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! less than 0.1 W/m2; or the iterative steps over 40. ! ! !USES: - use shr_const_mod , only : SHR_CONST_RGAS + use shr_const_mod , only : SHR_CONST_RGAS, shr_const_pi use clm_time_manager , only : get_step_size_real, get_prev_date,is_end_curr_day - use clm_varcon , only : sb, cpair, hvap, vkc, grav, denice + use clm_varcon , only : sb, cpair, hvap, vkc, grav, denice, c_to_b use clm_varcon , only : denh2o, tfrz, tlsai_crit, alpha_aero use clm_varcon , only : c14ratio + use clm_varcon , only : c_water, c_dry_biomass use perf_mod , only : t_startf, t_stopf use QSatMod , only : QSat use CLMFatesInterfaceMod, only : hlm_fates_interface_type @@ -279,39 +283,34 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, real(r8) :: dtime ! land model time step (sec) real(r8) :: zldis(bounds%begp:bounds%endp) ! reference height "minus" zero displacement height [m] - real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory real(r8) :: wc ! convective velocity [m/s] real(r8) :: dth(bounds%begp:bounds%endp) ! diff of virtual temp. between ref. height and surface real(r8) :: dthv(bounds%begp:bounds%endp) ! diff of vir. poten. temp. between ref. height and surface real(r8) :: dqh(bounds%begp:bounds%endp) ! diff of humidity between ref. height and surface - real(r8) :: obu(bounds%begp:bounds%endp) ! Monin-Obukhov length (m) - real(r8) :: um(bounds%begp:bounds%endp) ! wind speed including the stablity effect [m/s] real(r8) :: ur(bounds%begp:bounds%endp) ! wind speed at reference height [m/s] - real(r8) :: uaf(bounds%begp:bounds%endp) ! velocity of air within foliage [m/s] real(r8) :: temp1(bounds%begp:bounds%endp) ! relation for potential temperature profile real(r8) :: temp12m(bounds%begp:bounds%endp) ! relation for potential temperature profile applied at 2-m real(r8) :: temp2(bounds%begp:bounds%endp) ! relation for specific humidity profile real(r8) :: temp22m(bounds%begp:bounds%endp) ! relation for specific humidity profile applied at 2-m - real(r8) :: ustar(bounds%begp:bounds%endp) ! friction velocity [m/s] real(r8) :: tstar ! temperature scaling parameter real(r8) :: qstar ! moisture scaling parameter real(r8) :: thvstar ! virtual potential temperature scaling parameter - real(r8) :: taf(bounds%begp:bounds%endp) ! air temperature within canopy space [K] - real(r8) :: qaf(bounds%begp:bounds%endp) ! humidity of canopy air [kg/kg] real(r8) :: rpp ! fraction of potential evaporation from leaf [-] real(r8) :: rppdry ! fraction of potential evaporation through transp [-] real(r8) :: cf ! heat transfer coefficient from leaves [-] real(r8) :: rb(bounds%begp:bounds%endp) ! leaf boundary layer resistance [s/m] - real(r8) :: rah(bounds%begp:bounds%endp,2) ! thermal resistance [s/m] - real(r8) :: raw(bounds%begp:bounds%endp,2) ! moisture resistance [s/m] + real(r8) :: rah(bounds%begp:bounds%endp,2) ! thermal resistance [s/m] (air, ground) + real(r8) :: raw(bounds%begp:bounds%endp,2) ! moisture resistance [s/m] (air, ground) real(r8) :: wta ! heat conductance for air [m/s] real(r8) :: wtg(bounds%begp:bounds%endp) ! heat conductance for ground [m/s] real(r8) :: wtl ! heat conductance for leaf [m/s] + real(r8) :: wtstem ! heat conductance for stem [m/s] real(r8) :: wta0(bounds%begp:bounds%endp) ! normalized heat conductance for air [-] real(r8) :: wtl0(bounds%begp:bounds%endp) ! normalized heat conductance for leaf [-] real(r8) :: wtg0 ! normalized heat conductance for ground [-] + real(r8) :: wtstem0(bounds%begp:bounds%endp) ! normalized heat conductance for stem [-] real(r8) :: wtal(bounds%begp:bounds%endp) ! normalized heat conductance for air and leaf [-] - real(r8) :: wtga ! normalized heat cond. for air and ground [-] + real(r8) :: wtga(bounds%begp:bounds%endp) ! normalized heat cond. for air and ground [-] real(r8) :: wtaq ! latent heat conductance for air [m/s] real(r8) :: wtlq ! latent heat conductance for leaf [m/s] real(r8) :: wtgq(bounds%begp:bounds%endp) ! latent heat conductance for ground [m/s] @@ -343,6 +342,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, real(r8) :: efsh ! sensible heat from leaf [mm/s] real(r8) :: obuold(bounds%begp:bounds%endp) ! monin-obukhov length from previous iteration real(r8) :: tlbef(bounds%begp:bounds%endp) ! leaf temperature from previous iteration [K] + real(r8) :: tl_ini(bounds%begp:bounds%endp) ! leaf temperature from beginning of time step [K] + real(r8) :: ts_ini(bounds%begp:bounds%endp) ! stem temperature from beginning of time step [K] real(r8) :: ecidif ! excess energies [W/m2] real(r8) :: err(bounds%begp:bounds%endp) ! balance error real(r8) :: erre ! balance error @@ -402,14 +403,32 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, real(r8) :: dt_veg_temp(bounds%begp:bounds%endp) integer :: iv logical :: is_end_day ! is end of current day - real(r8) :: tl_ini ! leaf temperature from beginning of time step [K] - real(r8) :: ts_ini ! stem temperature from beginning of time step [K] - real(r8) :: cp_leaf !heat capacity of leaves - real(r8) :: dt_stem !change in stem temperature - real(r8) :: frac_rad_abs_by_stem !fraction of incoming radiation absorbed by stems - real(r8) :: lw_stem !internal longwave stem - real(r8) :: lw_leaf !internal longwave leaf - real(r8) :: sa_internal !min(sa_stem,sa_leaf) + real(r8) :: dbh(bounds%begp:bounds%endp) ! diameter at breast height of vegetation + real(r8) :: cp_leaf(bounds%begp:bounds%endp) ! heat capacity of leaves + real(r8) :: cp_stem(bounds%begp:bounds%endp) ! heat capacity of stems + real(r8) :: rstem(bounds%begp:bounds%endp) ! stem resistance to heat transfer + real(r8) :: dt_stem(bounds%begp:bounds%endp) ! change in stem temperature + real(r8) :: frac_rad_abs_by_stem(bounds%begp:bounds%endp) ! fraction of incoming radiation absorbed by stems + real(r8) :: lw_stem(bounds%begp:bounds%endp) ! internal longwave stem + real(r8) :: lw_leaf(bounds%begp:bounds%endp) ! internal longwave leaf + real(r8) :: sa_stem(bounds%begp:bounds%endp) ! surface area stem m2/m2_ground + real(r8) :: sa_leaf(bounds%begp:bounds%endp) ! surface area leaf m2/m2_ground + real(r8) :: sa_internal(bounds%begp:bounds%endp) ! min(sa_stem,sa_leaf) + real(r8) :: uuc(bounds%begp:bounds%endp) ! undercanopy windspeed + real(r8) :: carea_stem ! cross-sectional area of stem + real(r8) :: dlrad_leaf ! Downward longwave radition from leaf + + ! Indices for raw and rah + integer, parameter :: above_canopy = 1 ! Above canopy + integer, parameter :: below_canopy = 2 ! Below canopy + ! Biomass heat storage tuning parameters + ! These parameters can be used to account for differences + ! in vegetation shape. + real(r8), parameter :: k_vert = 0.1_r8 !vertical distribution of stem + real(r8), parameter :: k_cyl_vol = 1.0_r8 !departure from cylindrical volume + real(r8), parameter :: k_cyl_area = 1.0_r8 !departure from cylindrical area + real(r8), parameter :: k_internal = 0.0_r8 !self-absorbtion of leaf/stem longwave + real(r8), parameter :: min_stem_diameter = 0.05_r8 !minimum stem diameter for which to calculate stem interactions integer :: dummy_to_make_pgi_happy !------------------------------------------------------------------------------ @@ -417,16 +436,24 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, SHR_ASSERT_ALL_FL((ubound(downreg_patch) == (/bounds%endp/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(leafn_patch) == (/bounds%endp/)), sourcefile, __LINE__) - associate( & - soilresis => soilstate_inst%soilresis_col , & ! Input: [real(r8) (:) ] soil evaporative resistance - snl => col%snl , & ! Input: [integer (:) ] number of snow layers - dayl => grc%dayl , & ! Input: [real(r8) (:) ] daylength (s) - max_dayl => grc%max_dayl , & ! Input: [real(r8) (:) ] maximum daylength for this grid cell (s) - - dleaf => pftcon%dleaf , & ! Input: characteristic leaf dimension (m) + associate( & + t_stem => temperature_inst%t_stem_patch , & ! Output: [real(r8) (:) ] stem temperature (Kelvin) + dhsdt_canopy => energyflux_inst%dhsdt_canopy_patch , & ! Output: [real(r8) (:) ] change in heat storage of stem (W/m**2) [+ to atm] + soilresis => soilstate_inst%soilresis_col , & ! Input: [real(r8) (:) ] soil evaporative resistance + snl => col%snl , & ! Input: [integer (:) ] number of snow layers + dayl => grc%dayl , & ! Input: [real(r8) (:) ] daylength (s) + max_dayl => grc%max_dayl , & ! Input: [real(r8) (:) ] maximum daylength for this grid cell (s) + is_tree => pftcon%is_tree , & ! Input: tree patch or not + is_shrub => pftcon%is_shrub , & ! Input: shrub patch or not + dleaf => pftcon%dleaf , & ! Input: characteristic leaf dimension (m) + dbh_param => pftcon%dbh , & ! Input: diameter at brest height (m) + fbw => pftcon%fbw , & ! Input: fraction of biomass that is water + nstem => pftcon%nstem , & ! Input: stem number density (#ind/m2) + rstem_per_dbh => pftcon%rstem_per_dbh , & ! Input: stem resistance per stem diameter (s/m**2) + wood_density => pftcon%wood_density , & ! Input: dry wood density (kg/m3) forc_lwrad => atm2lnd_inst%forc_lwrad_downscaled_col , & ! Input: [real(r8) (:) ] downward infrared (longwave) radiation (W/m**2) - forc_q => wateratm2lndbulk_inst%forc_q_downscaled_col , & ! Input: [real(r8) (:) ] atmospheric specific humidity (kg/kg) + forc_q => wateratm2lndbulk_inst%forc_q_downscaled_col , & ! Input: [real(r8) (:) ] atmospheric specific humidity (kg/kg) forc_pbot => atm2lnd_inst%forc_pbot_downscaled_col , & ! Input: [real(r8) (:) ] atmospheric pressure (Pa) forc_th => atm2lnd_inst%forc_th_downscaled_col , & ! Input: [real(r8) (:) ] atmospheric potential temperature (Kelvin) forc_rho => atm2lnd_inst%forc_rho_downscaled_col , & ! Input: [real(r8) (:) ] density (kg/m**3) @@ -476,6 +503,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit leaf area laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded leaf area displa => canopystate_inst%displa_patch , & ! Input: [real(r8) (:) ] displacement height (m) + stem_biomass => canopystate_inst%stem_biomass_patch , & ! Output: [real(r8) (:) ] Aboveground stem biomass (kg/m**2) + leaf_biomass => canopystate_inst%leaf_biomass_patch , & ! Output: [real(r8) (:) ] Aboveground leaf biomass (kg/m**2) htop => canopystate_inst%htop_patch , & ! Input: [real(r8) (:) ] canopy top(m) dleaf_patch => canopystate_inst%dleaf_patch , & ! Output: [real(r8) (:) ] mean leaf diameter for this patch/pft @@ -555,8 +584,23 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, eflx_sh_snow => energyflux_inst%eflx_sh_snow_patch , & ! Output: [real(r8) (:) ] sensible heat flux from snow (W/m**2) [+ to atm] eflx_sh_h2osfc => energyflux_inst%eflx_sh_h2osfc_patch , & ! Output: [real(r8) (:) ] sensible heat flux from soil (W/m**2) [+ to atm] eflx_sh_soil => energyflux_inst%eflx_sh_soil_patch , & ! Output: [real(r8) (:) ] sensible heat flux from soil (W/m**2) [+ to atm] + eflx_sh_stem => energyflux_inst%eflx_sh_stem_patch , & ! Output: [real(r8) (:) ] sensible heat flux from stems (W/m**2) [+ to atm] eflx_sh_veg => energyflux_inst%eflx_sh_veg_patch , & ! Output: [real(r8) (:) ] sensible heat flux from leaves (W/m**2) [+ to atm] eflx_sh_grnd => energyflux_inst%eflx_sh_grnd_patch , & ! Output: [real(r8) (:) ] sensible heat flux from ground (W/m**2) [+ to atm] + rah1 => frictionvel_inst%rah1_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance [s/m] + rah2 => frictionvel_inst%rah2_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance [s/m] + raw1 => frictionvel_inst%raw1_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance [s/m] + raw2 => frictionvel_inst%raw2_patch , & ! Output: [real(r8) (:) ] aerodynamical resistance [s/m] + ustar => frictionvel_inst%ustar_patch , & ! Output: [real(r8) (:) ] friction velocity [m/s] + um => frictionvel_inst%um_patch , & ! Output: [real(r8) (:) ] wind speed including the stablity effect [m/s] + uaf => frictionvel_inst%uaf_patch , & ! Output: [real(r8) (:) ] canopy air speed [m/s] + taf => frictionvel_inst%taf_patch , & ! Output: [real(r8) (:) ] canopy air temperature [K] + qaf => frictionvel_inst%qaf_patch , & ! Output: [real(r8) (:) ] canopy air humidity [kg/kg] + obu => frictionvel_inst%obu_patch , & ! Output: [real(r8) (:) ] Monin-Obukhov length [m] + zeta => frictionvel_inst%zeta_patch , & ! Output: [real(r8) (:) ] dimensionless stability parameter + vpd => frictionvel_inst%vpd_patch , & ! Output: [real(r8) (:) ] vapor pressure deficit [Pa] + num_iter => frictionvel_inst%num_iter_patch , & ! Output: [real(r8) (:) ] number of iterations + begp => bounds%begp , & endp => bounds%endp , & begg => bounds%begg , & @@ -634,15 +678,104 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, wtaq0(p) = 0._r8 obuold(p) = 0._r8 btran(p) = btran0 + dhsdt_canopy(p) = 0._r8 + eflx_sh_stem(p) = 0._r8 end do - frac_rad_abs_by_stem = 0._r8 - lw_leaf = 0._r8 - cp_leaf = 0._r8 - lw_stem = 0._r8 - sa_internal = 0._r8 - dt_stem = 0._r8 - tl_ini = 0._r8 - ts_ini = 0._r8 + ! + ! Calculate biomass heat capacities + ! + if(use_biomass_heat_storage) then +bioms: do f = 1, fn + p = filterp(f) + + ! fraction of stem receiving incoming radiation + frac_rad_abs_by_stem(p) = (esai(p))/(elai(p)+esai(p)) + + ! when elai = 0, do not multiply by k_vert (i.e. frac_rad_abs_by_stem = 1) + if(elai(p) > 0._r8) frac_rad_abs_by_stem(p) = k_vert * frac_rad_abs_by_stem(p) + + ! if using Satellite Phenology mode, use values in parameter file + ! otherwise calculate dbh from stem biomass + if(use_cn) then + if(stem_biomass(p) > 0._r8) then + dbh(p) = 2._r8 * sqrt(stem_biomass(p) * (1._r8 - fbw(patch%itype(p))) & + / ( shr_const_pi * htop(p) * k_cyl_vol & + * nstem(patch%itype(p)) * wood_density(patch%itype(p)))) + else + dbh(p) = 0._r8 + endif + else + dbh(p) = dbh_param(patch%itype(p)) + endif + + ! leaf and stem surface area + sa_leaf(p) = elai(p) + ! double in spirit of full surface area for sensible heat + sa_leaf(p) = 2.*sa_leaf(p) + + ! Surface area for stem + sa_stem(p) = nstem(patch%itype(p))*(htop(p)*shr_const_pi*dbh(p)) + ! adjust for departure of cylindrical stem model + sa_stem(p) = k_cyl_area * sa_stem(p) + + ! + ! only calculate separate leaf/stem heat capacity for trees + ! and shrubs if dbh is greater than some minimum value + ! (set surface area for stem, and fraction absorbed by stem to zero) + if(.not.(is_tree(patch%itype(p)) .or. is_shrub(patch%itype(p))) & + .or. dbh(p) < min_stem_diameter) then + frac_rad_abs_by_stem(p) = 0.0 + sa_stem(p) = 0.0 + endif + + ! cross-sectional area of stems + carea_stem = shr_const_pi * (dbh(p)*0.5)**2 + + ! if using Satellite Phenology mode, calculate leaf and stem biomass + if(.not. use_cn) then + ! boreal needleleaf lma*c2b ~ 0.25 kg dry mass/m2(leaf) + leaf_biomass(p) = 0.25_r8 * max(0.01_r8, sa_leaf(p)) & + / (1.-fbw(patch%itype(p))) + stem_biomass(p) = carea_stem * htop(p) * k_cyl_vol & + * nstem(patch%itype(p)) * wood_density(patch%itype(p)) & + /(1.-fbw(patch%itype(p))) + endif + + ! internal longwave fluxes between leaf and stem + ! (use same area of interaction i.e. ignore leaf <-> leaf) + sa_internal(p) = min(sa_leaf(p),sa_stem(p)) + sa_internal(p) = k_internal * sa_internal(p) + + ! calculate specify heat capacity of vegetation + ! as weighted averaged of dry biomass and water + ! lma_dry has units of kg dry mass/m2 here + ! (Appendix B of Bonan et al., GMD, 2018) + + cp_leaf(p) = leaf_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + + ! cp-stem will have units J/k/ground_area + cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + ! adjust for departure from cylindrical stem model + cp_stem(p) = k_cyl_vol * cp_stem(p) + + ! resistance between internal stem temperature and canopy air + rstem(p) = rstem_per_dbh(patch%itype(p))*dbh(p) + + enddo bioms + else + ! Otherwise set biomass heat storage terms to zero + do f = 1, fn + p = filterp(f) + sa_leaf(p) = (elai(p)+esai(p)) + frac_rad_abs_by_stem(p) = 0._r8 + sa_stem(p) = 0._r8 + sa_internal(p) = 0._r8 + cp_leaf(p) = 0._r8 + cp_stem(p) = 0._r8 + rstem(p) = 0._r8 + end do + end if + ! calculate daylength control for Vcmax do f = 1, fn @@ -721,7 +854,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, end if - ! Modify aerodynamic parameters for sparse/dense canopy (X. Zeng) + !! Modify aerodynamic parameters for sparse/dense canopy (X. Zeng) + do f = 1, fn p = filterp(f) c = patch%column(p) @@ -799,6 +933,11 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Initialize Monin-Obukhov length and wind speed call frictionvel_inst%MoninObukIni(ur(p), thv(c), dthv(p), zldis(p), z0mv(p), um(p), obu(p)) + num_iter(p) = 0 + + ! Record initial veg/stem temperatures + tl_ini(p) = t_veg(p) + ts_ini(p) = t_stem(p) end do @@ -832,13 +971,16 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Determine aerodynamic resistances ram1(p) = 1._r8/(ustar(p)*ustar(p)/um(p)) - rah(p,1) = 1._r8/(temp1(p)*ustar(p)) - raw(p,1) = 1._r8/(temp2(p)*ustar(p)) + rah(p,above_canopy) = 1._r8/(temp1(p)*ustar(p)) + raw(p,above_canopy) = 1._r8/(temp2(p)*ustar(p)) ! Bulk boundary layer resistance of leaves uaf(p) = um(p)*sqrt( 1._r8/(ram1(p)*um(p)) ) + ! empirical undercanopy wind speed + uuc(p) = min(0.4_r8,(0.03_r8*um(p)/ustar(p))) + ! Use pft parameter for leaf characteristic width ! dleaf_patch if this is not an fates patch. ! Otherwise, the value has already been loaded @@ -865,8 +1007,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ri = ( grav*htop(p) * (taf(p) - t_grnd(c)) ) / (taf(p) * uaf(p) **2.00_r8) - !! modify csoilc value (0.004) if the under-canopy is in stable condition - + ! modify csoilc value (0.004) if the under-canopy is in stable condition if (use_undercanopy_stability .and. (taf(p) - t_grnd(c) ) > 0._r8) then ! decrease the value of csoilc by dividing it with (1+gamma*min(S, 10.0)) ! ria ("gmanna" in Sakaguchi&Zeng, 2008) is a constant (=0.5) @@ -878,10 +1019,16 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, !! Sakaguchi changes for stability formulation ends here - rah(p,2) = 1._r8/(csoilcn*uaf(p)) - raw(p,2) = rah(p,2) + if (use_biomass_heat_storage) then + ! use uuc for ground fluxes (keep uaf for canopy terms) + rah(p,below_canopy) = 1._r8/(csoilcn*uuc(p)) + else + rah(p,below_canopy) = 1._r8/(csoilcn*uaf(p)) + endif + + raw(p,below_canopy) = rah(p,below_canopy) if (use_lch4) then - grnd_ch4_cond(p) = 1._r8/(raw(p,1)+raw(p,2)) + grnd_ch4_cond(p) = 1._r8/(raw(p,above_canopy)+raw(p,below_canopy)) end if ! Stomatal resistances for sunlit and shaded fractions of canopy. @@ -890,6 +1037,13 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, svpts(p) = el(p) ! pa eah(p) = forc_pbot(c) * qaf(p) / 0.622_r8 ! pa rhaf(p) = eah(p)/svpts(p) + ! variables for history fields + rah1(p) = rah(p,above_canopy) + raw1(p) = raw(p,above_canopy) + rah2(p) = rah(p,below_canopy) + raw2(p) = raw(p,below_canopy) + vpd(p) = max((svpts(p) - eah(p)), 50._r8) * 0.001_r8 + end do if ( use_fates ) then @@ -947,17 +1101,24 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Sensible heat conductance for air, leaf and ground ! Moved the original subroutine in-line... - wta = 1._r8/rah(p,1) ! air - wtl = (elai(p)+esai(p))/rb(p) ! leaf - wtg(p) = 1._r8/rah(p,2) ! ground - wtshi = 1._r8/(wta+wtl+wtg(p)) + wta = 1._r8/rah(p,above_canopy) ! air + wtl = sa_leaf(p)/rb(p) ! leaf + wtg(p) = 1._r8/rah(p,below_canopy) ! ground + wtstem = sa_stem(p)/(rstem(p) + rb(p)) ! stem + + wtshi = 1._r8/(wta+wtl+wtstem+wtg(p)) ! air, leaf, stem and ground wtl0(p) = wtl*wtshi ! leaf wtg0 = wtg(p)*wtshi ! ground wta0(p) = wta*wtshi ! air - wtga = wta0(p)+wtg0 ! ground + air - wtal(p) = wta0(p)+wtl0(p) ! air + leaf + wtstem0(p) = wtstem*wtshi ! stem + wtga(p) = wta0(p)+wtg0+wtstem0(p) ! ground + air + stem + wtal(p) = wta0(p)+wtl0(p)+wtstem0(p) ! air + leaf + stem + + ! internal longwave fluxes between leaf and stem + lw_stem(p) = sa_internal(p) * emv(p) * sb * t_stem(p)**4 + lw_leaf(p) = sa_internal(p) * emv(p) * sb * t_veg(p)**4 ! Fraction of potential evaporation from leaf @@ -972,7 +1133,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, canopy_cond(p) = (laisun(p)/(rb(p)+rssun(p)) + laisha(p)/(rb(p)+rssha(p)))/max(elai(p), 0.01_r8) end if - efpot = forc_rho(c)*wtl*(qsatl(p)-qaf(p)) + efpot = forc_rho(c)*((elai(p)+esai(p))/rb(p))*(qsatl(p)-qaf(p)) h2ocan = liqcan(p) + snocan(p) ! When the hydraulic stress parameterization is active calculate rpp @@ -1013,7 +1174,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Air has same conductance for both sensible and latent heat. ! Moved the original subroutine in-line... - wtaq = frac_veg_nosno(p)/raw(p,1) ! air + wtaq = frac_veg_nosno(p)/raw(p,above_canopy) ! air wtlq = frac_veg_nosno(p)*(elai(p)+esai(p))/rb(p) * rpp ! leaf !Litter layer resistance. Added by K.Sakaguchi @@ -1024,13 +1185,13 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! add litter resistance and Lee and Pielke 1992 beta if (delq(p) < 0._r8) then !dew. Do not apply beta for negative flux (follow old rsoil) - wtgq(p) = frac_veg_nosno(p)/(raw(p,2)+rdl) + wtgq(p) = frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) else if (do_soilevap_beta()) then - wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,2)+rdl) + wtgq(p) = soilbeta(c)*frac_veg_nosno(p)/(raw(p,below_canopy)+rdl) endif if (do_soil_resistance_sl14()) then - wtgq(p) = frac_veg_nosno(p)/(raw(p,2)+soilresis(c)) + wtgq(p) = frac_veg_nosno(p)/(raw(p,below_canopy)+soilresis(c)) endif end if @@ -1046,7 +1207,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, dc1 = forc_rho(c)*cpair*wtl dc2 = hvap*forc_rho(c)*wtlq - efsh = dc1*(wtga*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)) + efsh = dc1*(wtga(p)*t_veg(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)-wtstem0(p)*t_stem(p)) + eflx_sh_stem(p) = forc_rho(c)*cpair*wtstem*((wta0(p)+wtg0+wtl0(p))*t_stem(p)-wtg0*t_grnd(c)-wta0(p)*thm(p)-wtl0(p)*t_veg(p)) efe(p) = dc2*(wtgaq*qsatl(p)-wtgq0*qg(c)-wtaq0(p)*forc_q(c)) ! Evaporation flux from foliage @@ -1062,28 +1224,30 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, +(1._r8-frac_sno(c)-frac_h2osfc(c))*t_soisno(c,1)**4 & +frac_h2osfc(c)*t_h2osfc(c)**4) - dt_veg(p) = ((1._r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) & + dt_veg(p) = ((1._r8-frac_rad_abs_by_stem(p))*(sabv(p) + air(p) & + bir(p)*t_veg(p)**4 + cir(p)*lw_grnd) & - - efsh - efe(p) - lw_leaf + lw_stem & - - (cp_leaf/dtime)*(t_veg(p) - tl_ini)) & - / ((1._r8-frac_rad_abs_by_stem)*(- 4._r8*bir(p)*t_veg(p)**3 & - + 4._r8*sa_internal*emv(p)*sb*t_veg(p)**3 & - +dc1*wtga+dc2*wtgaq*qsatldT(p))+ cp_leaf/dtime) + - efsh - efe(p) - lw_leaf(p) + lw_stem(p) & + - (cp_leaf(p)/dtime)*(t_veg(p) - tl_ini(p))) & + / ((1._r8-frac_rad_abs_by_stem(p))*(- 4._r8*bir(p)*t_veg(p)**3) & + + 4._r8*sa_internal(p)*emv(p)*sb*t_veg(p)**3 & + +dc1*wtga(p) +dc2*wtgaq*qsatldT(p)+ cp_leaf(p)/dtime) + t_veg(p) = tlbef(p) + dt_veg(p) + dels = dt_veg(p) del(p) = abs(dels) err(p) = 0._r8 if (del(p) > delmax) then dt_veg(p) = delmax*dels/del(p) t_veg(p) = tlbef(p) + dt_veg(p) - err(p) = (1._r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) & + err(p) = (1._r8-frac_rad_abs_by_stem(p))*(sabv(p) + air(p) & + bir(p)*tlbef(p)**3*(tlbef(p) + & 4._r8*dt_veg(p)) + cir(p)*lw_grnd) & - -sa_internal*emv(p)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) & - + lw_stem & - - (efsh + dc1*wtga*dt_veg(p)) - (efe(p) + & + -sa_internal(p)*emv(p)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) & + + lw_stem(p) & + - (efsh + dc1*wtga(p)*dt_veg(p)) - (efe(p) + & dc2*wtgaq*qsatldT(p)*dt_veg(p)) & - - (cp_leaf/dtime)*(t_veg(p) - tl_ini) + - (cp_leaf(p)/dtime)*(t_veg(p) - tl_ini(p)) end if ! Fluxes from leaves to canopy space @@ -1091,7 +1255,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! result in an imbalance in "hvap*qflx_evap_veg" and ! "efe + dc2*wtgaq*qsatdt_veg" - efpot = forc_rho(c)*wtl*(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) & + efpot = forc_rho(c)*((elai(p)+esai(p))/rb(p)) & + *(wtgaq*(qsatl(p)+qsatldT(p)*dt_veg(p)) & -wtgq0*qg(c)-wtaq0(p)*forc_q(c)) qflx_evap_veg(p) = rpp*efpot @@ -1124,7 +1289,11 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! The energy loss due to above two limits is added to ! the sensible heat flux. - eflx_sh_veg(p) = efsh + dc1*wtga*dt_veg(p) + err(p) + erre + hvap*ecidif + eflx_sh_veg(p) = efsh + dc1*wtga(p)*dt_veg(p) + err(p) + erre + hvap*ecidif + + ! Update SH and lw_leaf for changes in t_veg + eflx_sh_stem(p) = eflx_sh_stem(p) + forc_rho(c)*cpair*wtstem*(-wtl0(p)*dt_veg(p)) + lw_leaf(p) = sa_internal(p)*emv(p)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) ! Re-calculate saturated vapor pressure, specific humidity, and their ! derivatives at the leaf surface @@ -1137,7 +1306,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! temperature, canopy vapor pressure, aerodynamic temperature, and ! Monin-Obukhov stability parameter for next iteration. - taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p) + taf(p) = wtg0*t_grnd(c) + wta0(p)*thm(p) + wtl0(p)*t_veg(p) + wtstem0(p)*t_stem(p) qaf(p) = wtlq0(p)*qsatl(p) + wtgq0*qg(c) + forc_q(c)*wtaq0(p) ! Update Monin-Obukhov length and wind speed including the @@ -1151,17 +1320,22 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, qstar = temp2(p)*dqh(p) thvstar = tstar*(1._r8+0.61_r8*forc_q(c)) + 0.61_r8*forc_th(c)*qstar - zeta = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) + zeta(p) = zldis(p)*vkc*grav*thvstar/(ustar(p)**2*thv(c)) - if (zeta >= 0._r8) then !stable - zeta = min(zetamax,max(zeta,0.01_r8)) + if (zeta(p) >= 0._r8) then !stable + ! remove stability cap when biomass heat storage is active + if(use_biomass_heat_storage) then + zeta(p) = min(100._r8,max(zeta(p),0.01_r8)) + else + zeta(p) = min(zetamax,max(zeta(p),0.01_r8)) + endif um(p) = max(ur(p),0.1_r8) else !unstable - zeta = max(-100._r8,min(zeta,-0.01_r8)) + zeta(p) = max(-100._r8,min(zeta(p),-0.01_r8)) wc = beta*(-grav*ustar(p)*thvstar*zii/thv(c))**0.333_r8 um(p) = sqrt(ur(p)*ur(p)+wc*wc) end if - obu(p) = zldis(p)/zeta + obu(p) = zldis(p)/zeta(p) if (obuold(p)*obu(p) < 0._r8) nmozsgn(p) = nmozsgn(p)+1 if (nmozsgn(p) >= 4) obu(p) = zldis(p)/(-0.01_r8) @@ -1172,6 +1346,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Test for convergence itlef = itlef+1 + num_iter(p) = itlef if (itlef > itmin) then do f = 1, fn p = filterp(f) @@ -1189,7 +1364,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, end if end do end if - end do ITERATION ! End stability iteration call t_stopf('can_iter') @@ -1207,26 +1381,47 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, +(1._r8-frac_sno(c)-frac_h2osfc(c))*t_soisno(c,1)**4 & +frac_h2osfc(c)*t_h2osfc(c)**4) - err(p) = (1.0_r8-frac_rad_abs_by_stem)*(sabv(p) + air(p) + bir(p)*tlbef(p)**3 & + err(p) = (1.0_r8-frac_rad_abs_by_stem(p))*(sabv(p) + air(p) + bir(p)*tlbef(p)**3 & *(tlbef(p) + 4._r8*dt_veg(p)) + cir(p)*lw_grnd) & - - lw_leaf + lw_stem - eflx_sh_veg(p) - hvap*qflx_evap_veg(p) & - - ((t_veg(p)-tl_ini)*cp_leaf/dtime) + - lw_leaf(p) + lw_stem(p) - eflx_sh_veg(p) - hvap*qflx_evap_veg(p) & + - ((t_veg(p)-tl_ini(p))*cp_leaf(p)/dtime) + + ! Update stem temperature; adjust outgoing longwave + ! does not account for changes in SH or internal LW, + ! as that would change result for t_veg above + if (use_biomass_heat_storage) then + if (stem_biomass(p) > 0._r8) then + dt_stem(p) = (frac_rad_abs_by_stem(p)*(sabv(p) + air(p) + bir(p)*ts_ini(p)**4 & + + cir(p)*lw_grnd) - eflx_sh_stem(p) & + + lw_leaf(p)- lw_stem(p))/(cp_stem(p)/dtime & + - frac_rad_abs_by_stem(p)*bir(p)*4.*ts_ini(p)**3) + else + dt_stem(p) = 0._r8 + endif + + dhsdt_canopy(p) = dt_stem(p)*cp_stem(p)/dtime & + +(t_veg(p)-tl_ini(p))*cp_leaf(p)/dtime + + t_stem(p) = t_stem(p) + dt_stem(p) + else + dt_stem(p) = 0._r8 + endif + + delt = wtal(p)*t_grnd(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p)-wtstem0(p)*t_stem(p) ! Fluxes from ground to canopy space - delt = wtal(p)*t_grnd(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p) taux(p) = -forc_rho(c)*forc_u(g)/ram1(p) tauy(p) = -forc_rho(c)*forc_v(g)/ram1(p) eflx_sh_grnd(p) = cpair*forc_rho(c)*wtg(p)*delt ! compute individual sensible heat fluxes - delt_snow = wtal(p)*t_soisno(c,snl(c)+1)-wtl0(p)*t_veg(p)-wta0(p)*thm(p) - eflx_sh_snow(p) = cpair*forc_rho(c)*wtg(p)*delt_snow + delt_snow = wtal(p)*t_soisno(c,snl(c)+1)-wtl0(p)*t_veg(p)-wta0(p)*thm(p)-wtstem0(p)*t_stem(p) + delt_soil = wtal(p)*t_soisno(c,1)-wtl0(p)*t_veg(p)-wta0(p)*thm(p)-wtstem0(p)*t_stem(p) + delt_h2osfc = wtal(p)*t_h2osfc(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p)-wtstem0(p)*t_stem(p) - delt_soil = wtal(p)*t_soisno(c,1)-wtl0(p)*t_veg(p)-wta0(p)*thm(p) + eflx_sh_snow(p) = cpair*forc_rho(c)*wtg(p)*delt_snow eflx_sh_soil(p) = cpair*forc_rho(c)*wtg(p)*delt_soil - - delt_h2osfc = wtal(p)*t_h2osfc(c)-wtl0(p)*t_veg(p)-wta0(p)*thm(p) eflx_sh_h2osfc(p) = cpair*forc_rho(c)*wtg(p)*delt_h2osfc qflx_evap_soi(p) = forc_rho(c)*wtgq(p)*delq(p) @@ -1296,17 +1491,17 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(c) & + emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p)) & - *(1.0_r8-frac_rad_abs_by_stem) & - + emv(p)*emg(c)*sb*ts_ini**3*(ts_ini + 4._r8*dt_stem) & - *frac_rad_abs_by_stem + *(1.0_r8-frac_rad_abs_by_stem(p)) & + + emv(p)*emg(c)*sb*ts_ini(p)**3*(ts_ini(p) + 4._r8*dt_stem(p)) & + *frac_rad_abs_by_stem(p) ! Upward longwave radiation above the canopy ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(c) & + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb & - *tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))*(1._r8-frac_rad_abs_by_stem) & + *tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))*(1._r8-frac_rad_abs_by_stem(p)) & + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb & - *ts_ini**3*(ts_ini+ 4._r8*dt_stem)*frac_rad_abs_by_stem & + *ts_ini(p)**3*(ts_ini(p)+ 4._r8*dt_stem(p))*frac_rad_abs_by_stem(p) & + emg(c)*(1._r8-emv(p))*sb*lw_grnd) ! Calculate the skin temperature as a weighted sum of all the ground and vegetated fraction @@ -1431,7 +1626,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, write(iulog,*) 'energy balance in canopy ',p,', err=',err(p) end do - end associate diff --git a/src/biogeophys/CanopyStateType.F90 b/src/biogeophys/CanopyStateType.F90 index 8aa8a24f68..abb32c9acc 100644 --- a/src/biogeophys/CanopyStateType.F90 +++ b/src/biogeophys/CanopyStateType.F90 @@ -35,6 +35,8 @@ module CanopyStateType real(r8) , pointer :: laisha_z_patch (:,:) ! patch patch shaded leaf area for canopy layer real(r8) , pointer :: mlaidiff_patch (:) ! patch difference between lai month one and month two (for dry deposition of chemical tracers) real(r8) , pointer :: annlai_patch (:,:) ! patch 12 months of monthly lai from input data set (for dry deposition of chemical tracers) + real(r8) , pointer :: stem_biomass_patch (:) ! Aboveground stem biomass (kg/m**2) + real(r8) , pointer :: leaf_biomass_patch (:) ! Aboveground leaf biomass (kg/m**2) real(r8) , pointer :: htop_patch (:) ! patch canopy top (m) real(r8) , pointer :: hbot_patch (:) ! patch canopy bottom (m) real(r8) , pointer :: z0m_patch (:) ! patch momentum roughness length (m) @@ -121,6 +123,8 @@ subroutine InitAllocate(this, bounds) allocate(this%laisha_z_patch (begp:endp,1:nlevcan)) ; this%laisha_z_patch (:,:) = nan allocate(this%mlaidiff_patch (begp:endp)) ; this%mlaidiff_patch (:) = nan allocate(this%annlai_patch (12,begp:endp)) ; this%annlai_patch (:,:) = nan + allocate(this%stem_biomass_patch (begp:endp)) ; this%stem_biomass_patch (:) = nan + allocate(this%leaf_biomass_patch (begp:endp)) ; this%leaf_biomass_patch (:) = nan allocate(this%htop_patch (begp:endp)) ; this%htop_patch (:) = nan allocate(this%hbot_patch (begp:endp)) ; this%hbot_patch (:) = nan allocate(this%z0m_patch (begp:endp)) ; this%z0m_patch (:) = nan @@ -186,6 +190,16 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='shaded projected leaf area index', & ptr_patch=this%laisha_patch, set_urb=0._r8) + this%stem_biomass_patch(begp:endp) = spval + call hist_addfld1d (fname='AGSB', units='kg/m^2', & + avgflag='A', long_name='Aboveground stem biomass', & + ptr_patch=this%stem_biomass_patch, default='inactive') + + this%leaf_biomass_patch(begp:endp) = spval + call hist_addfld1d (fname='AGLB', units='kg/m^2', & + avgflag='A', long_name='Aboveground leaf biomass', & + ptr_patch=this%leaf_biomass_patch, default='inactive') + if (use_cn .or. use_fates) then this%fsun_patch(begp:endp) = spval call hist_addfld1d (fname='FSUN', units='proportion', & @@ -466,14 +480,15 @@ subroutine InitCold(this, bounds) do p = bounds%begp, bounds%endp l = patch%landunit(p) - this%frac_veg_nosno_patch(p) = 0._r8 - this%tlai_patch(p) = 0._r8 - this%tsai_patch(p) = 0._r8 - this%elai_patch(p) = 0._r8 - this%esai_patch(p) = 0._r8 - this%htop_patch(p) = 0._r8 - this%hbot_patch(p) = 0._r8 - this%vegwp_patch(p,:) = -2.5e4_r8 + this%tlai_patch(p) = 0._r8 + this%tsai_patch(p) = 0._r8 + this%elai_patch(p) = 0._r8 + this%esai_patch(p) = 0._r8 + this%stem_biomass_patch(p)= 0._r8 + this%leaf_biomass_patch(p)= 0._r8 + this%htop_patch(p) = 0._r8 + this%hbot_patch(p) = 0._r8 + this%vegwp_patch(p,:) = -2.5e4_r8 if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then this%laisun_patch(p) = 0._r8 @@ -528,6 +543,14 @@ subroutine Restart(this, bounds, ncid, flag) dim1name='pft', long_name='one-sided stem area index, with burying by snow', units='', & interpinic_flag='interp', readvar=readvar, data=this%esai_patch) + call restartvar(ncid=ncid, flag=flag, varname='stem_biomass', xtype=ncd_double, & + dim1name='pft', long_name='stem biomass', units='kg/m^2', & + interpinic_flag='interp', readvar=readvar, data=this%stem_biomass_patch) + + call restartvar(ncid=ncid, flag=flag, varname='leaf_biomass', xtype=ncd_double, & + dim1name='pft', long_name='leaf biomass', units='kg/m^2', & + interpinic_flag='interp', readvar=readvar, data=this%leaf_biomass_patch) + call restartvar(ncid=ncid, flag=flag, varname='htop', xtype=ncd_double, & dim1name='pft', long_name='canopy top', units='m', & interpinic_flag='interp', readvar=readvar, data=this%htop_patch) diff --git a/src/biogeophys/EnergyFluxType.F90 b/src/biogeophys/EnergyFluxType.F90 index 9b1a37a235..5634d26e50 100644 --- a/src/biogeophys/EnergyFluxType.F90 +++ b/src/biogeophys/EnergyFluxType.F90 @@ -8,6 +8,7 @@ module EnergyFluxType use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg use clm_varcon , only : spval + use clm_varctl , only : use_biomass_heat_storage use decompMod , only : bounds_type use LandunitType , only : lun use ColumnType , only : col @@ -21,6 +22,7 @@ module EnergyFluxType type, public :: energyflux_type ! Fluxes + real(r8), pointer :: eflx_sh_stem_patch (:) ! patch sensible heat flux from stem (W/m**2) [+ to atm] real(r8), pointer :: eflx_h2osfc_to_snow_col (:) ! col snow melt to h2osfc heat flux (W/m**2) real(r8), pointer :: eflx_sh_grnd_patch (:) ! patch sensible heat flux from ground (W/m**2) [+ to atm] real(r8), pointer :: eflx_sh_veg_patch (:) ! patch sensible heat flux from leaves (W/m**2) [+ to atm] @@ -100,6 +102,9 @@ module EnergyFluxType ! Latent heat real(r8), pointer :: htvp_col (:) ! latent heat of vapor of water (or sublimation) [j/kg] + ! Canopy heat + real(r8), pointer :: dhsdt_canopy_patch (:) ! patch change in heat content of canopy (leaf+stem) (W/m**2) [+ to atm] + ! Balance Checks real(r8), pointer :: errsoi_patch (:) ! soil/lake energy conservation error (W/m**2) real(r8), pointer :: errsoi_col (:) ! soil/lake energy conservation error (W/m**2) @@ -190,6 +195,7 @@ subroutine InitAllocate(this, bounds) allocate( this%eflx_sh_tot_u_patch (begp:endp)) ; this%eflx_sh_tot_u_patch (:) = nan allocate( this%eflx_sh_tot_r_patch (begp:endp)) ; this%eflx_sh_tot_r_patch (:) = nan allocate( this%eflx_sh_grnd_patch (begp:endp)) ; this%eflx_sh_grnd_patch (:) = nan + allocate( this%eflx_sh_stem_patch (begp:endp)) ; this%eflx_sh_stem_patch (:) = nan allocate( this%eflx_sh_veg_patch (begp:endp)) ; this%eflx_sh_veg_patch (:) = nan allocate( this%eflx_sh_precip_conversion_col(begc:endc)) ; this%eflx_sh_precip_conversion_col(:) = nan allocate( this%eflx_lh_tot_u_patch (begp:endp)) ; this%eflx_lh_tot_u_patch (:) = nan @@ -245,6 +251,8 @@ subroutine InitAllocate(this, bounds) allocate( this%htvp_col (begc:endc)) ; this%htvp_col (:) = nan + allocate( this%dhsdt_canopy_patch (begp:endp)) ; this%dhsdt_canopy_patch (:) = nan + allocate(this%rresis_patch (begp:endp,1:nlevgrnd)) ; this%rresis_patch (:,:) = nan allocate(this%btran_patch (begp:endp)) ; this%btran_patch (:) = nan allocate(this%btran_min_patch (begp:endp)) ; this%btran_min_patch (:) = nan @@ -431,6 +439,18 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp) avgflag='A', long_name='sensible heat from veg', & ptr_patch=this%eflx_sh_veg_patch, set_lake=0._r8, c2l_scale_type='urbanf') + if (use_biomass_heat_storage) then + this%eflx_sh_stem_patch(begp:endp) = spval + call hist_addfld1d (fname='FSH_STEM', units='W/m^2', & + avgflag='A', long_name='sensible heat from stem', & + ptr_patch=this%eflx_sh_stem_patch, c2l_scale_type='urbanf',default = 'inactive') + + this%dhsdt_canopy_patch(begp:endp) = spval + call hist_addfld1d (fname='DHSDT_CANOPY', units='W/m^2', & + avgflag='A', long_name='change in canopy heat storage', & + ptr_patch=this%dhsdt_canopy_patch, set_lake=0._r8, c2l_scale_type='urbanf',default='active') + endif + this%eflx_sh_grnd_patch(begp:endp) = spval call hist_addfld1d (fname='FSH_G', units='W/m^2', & avgflag='A', long_name='sensible heat from ground', & @@ -528,12 +548,10 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp) ptr_patch=this%cgrnds_patch, default='inactive', c2l_scale_type='urbanf') end if - if (use_cn) then - this%eflx_gnet_patch(begp:endp) = spval - call hist_addfld1d (fname='EFLX_GNET', units='W/m^2', & - avgflag='A', long_name='net heat flux into ground', & - ptr_patch=this%eflx_gnet_patch, default='inactive', c2l_scale_type='urbanf') - end if + this%eflx_gnet_patch(begp:endp) = spval + call hist_addfld1d (fname='EFLX_GNET', units='W/m^2', & + avgflag='A', long_name='net heat flux into ground', & + ptr_patch=this%eflx_gnet_patch, default='inactive', c2l_scale_type='urbanf') this%eflx_grnd_lake_patch(begp:endp) = spval call hist_addfld1d (fname='EFLX_GRND_LAKE', units='W/m^2', & diff --git a/src/biogeophys/FrictionVelocityMod.F90 b/src/biogeophys/FrictionVelocityMod.F90 index 19af9c9343..9efd167879 100644 --- a/src/biogeophys/FrictionVelocityMod.F90 +++ b/src/biogeophys/FrictionVelocityMod.F90 @@ -55,6 +55,20 @@ module FrictionVelocityMod real(r8), pointer, public :: z0mg_col (:) ! col roughness length over ground, momentum [m] real(r8), pointer, public :: z0hg_col (:) ! col roughness length over ground, sensible heat [m] real(r8), pointer, public :: z0qg_col (:) ! col roughness length over ground, latent heat [m] + ! variables to add history output from CanopyFluxesMod + real(r8), pointer, public :: rah1_patch (:) ! patch sensible heat flux resistance [s/m] + real(r8), pointer, public :: rah2_patch (:) ! patch below-canopy sensible heat flux resistance [s/m] + real(r8), pointer, public :: raw1_patch (:) ! patch moisture flux resistance [s/m] + real(r8), pointer, public :: raw2_patch (:) ! patch below-canopy moisture flux resistance [s/m] + real(r8), pointer, public :: ustar_patch (:) ! patch friction velocity [m/s] + real(r8), pointer, public :: um_patch (:) ! patch wind speed including the stablity effect [m/s] + real(r8), pointer, public :: uaf_patch (:) ! patch canopy air speed [m/s] + real(r8), pointer, public :: taf_patch (:) ! patch canopy air temperature [K] + real(r8), pointer, public :: qaf_patch (:) ! patch canopy humidity [kg/kg] + real(r8), pointer, public :: obu_patch (:) ! patch Monin-Obukhov length [m] + real(r8), pointer, public :: zeta_patch (:) ! patch dimensionless stability parameter + real(r8), pointer, public :: vpd_patch (:) ! patch vapor pressure deficit [Pa] + real(r8), pointer, public :: num_iter_patch (:) ! patch number of iterations real(r8), pointer, public :: z0m_actual_patch (:) ! patch roughness length actually used in flux calculations, momentum [m] contains @@ -138,6 +152,19 @@ subroutine InitAllocate(this, bounds) allocate(this%z0mg_col (begc:endc)) ; this%z0mg_col (:) = nan allocate(this%z0qg_col (begc:endc)) ; this%z0qg_col (:) = nan allocate(this%z0hg_col (begc:endc)) ; this%z0hg_col (:) = nan + allocate(this%rah1_patch (begp:endp)) ; this%rah1_patch (:) = nan + allocate(this%rah2_patch (begp:endp)) ; this%rah2_patch (:) = nan + allocate(this%raw1_patch (begp:endp)) ; this%raw1_patch (:) = nan + allocate(this%raw2_patch (begp:endp)) ; this%raw2_patch (:) = nan + allocate(this%um_patch (begp:endp)) ; this%um_patch (:) = nan + allocate(this%uaf_patch (begp:endp)) ; this%uaf_patch (:) = nan + allocate(this%taf_patch (begp:endp)) ; this%taf_patch (:) = nan + allocate(this%qaf_patch (begp:endp)) ; this%qaf_patch (:) = nan + allocate(this%ustar_patch (begp:endp)) ; this%ustar_patch (:) = nan + allocate(this%obu_patch (begp:endp)) ; this%obu_patch (:) = nan + allocate(this%zeta_patch (begp:endp)) ; this%zeta_patch (:) = nan + allocate(this%vpd_patch (begp:endp)) ; this%vpd_patch (:) = nan + allocate(this%num_iter_patch (begp:endp)) ; this%num_iter_patch (:) = nan allocate(this%z0m_actual_patch (begp:endp)) ; this%z0m_actual_patch (:) = nan end subroutine InitAllocate @@ -211,7 +238,63 @@ subroutine InitHistory(this, bounds) ptr_patch=this%fv_patch, default='inactive') end if - if (use_cn) then + call hist_addfld1d (fname='RAH1', units='s/m', & + avgflag='A', long_name='aerodynamical resistance ', & + ptr_patch=this%rah1_patch, default='inactive') + this%rah2_patch(begp:endp) = spval + call hist_addfld1d (fname='RAH2', units='s/m', & + avgflag='A', long_name='aerodynamical resistance ', & + ptr_patch=this%rah2_patch, default='inactive') + this%raw1_patch(begp:endp) = spval + call hist_addfld1d (fname='RAW1', units='s/m', & + avgflag='A', long_name='aerodynamical resistance ', & + ptr_patch=this%raw1_patch, default='inactive') + this%raw2_patch(begp:endp) = spval + call hist_addfld1d (fname='RAW2', units='s/m', & + avgflag='A', long_name='aerodynamical resistance ', & + ptr_patch=this%raw2_patch, default='inactive') + this%ustar_patch(begp:endp) = spval + call hist_addfld1d (fname='USTAR', units='s/m', & + avgflag='A', long_name='aerodynamical resistance ', & + ptr_patch=this%ustar_patch, default='inactive') + this%um_patch(begp:endp) = spval + call hist_addfld1d (fname='UM', units='m/s', & + avgflag='A', long_name='wind speed plus stability effect', & + ptr_patch=this%um_patch, default='inactive') + this%uaf_patch(begp:endp) = spval + call hist_addfld1d (fname='UAF', units='m/s', & + avgflag='A', long_name='canopy air speed ', & + ptr_patch=this%uaf_patch, default='inactive') + this%taf_patch(begp:endp) = spval + call hist_addfld1d (fname='TAF', units='K', & + avgflag='A', long_name='canopy air temperature', & + ptr_patch=this%taf_patch, default='inactive') + this%qaf_patch(begp:endp) = spval + call hist_addfld1d (fname='QAF', units='kg/kg', & + avgflag='A', long_name='canopy air humidity', & + ptr_patch=this%qaf_patch, default='inactive') + this%obu_patch(begp:endp) = spval + call hist_addfld1d (fname='OBU', units='m', & + avgflag='A', long_name='Monin-Obukhov length', & + ptr_patch=this%obu_patch, default='inactive') + this%zeta_patch(begp:endp) = spval + call hist_addfld1d (fname='ZETA', units='unitless', & + avgflag='A', long_name='dimensionless stability parameter', & + ptr_patch=this%zeta_patch, default='inactive') + this%vpd_patch(begp:endp) = spval + call hist_addfld1d (fname='VPD', units='Pa', & + avgflag='A', long_name='vpd', & + ptr_patch=this%vpd_patch, default='inactive') + this%num_iter_patch(begp:endp) = spval + call hist_addfld1d (fname='num_iter', units='unitless', & + avgflag='A', long_name='number of iterations', & + ptr_patch=this%num_iter_patch, default='inactive') + this%rb1_patch(begp:endp) = spval + call hist_addfld1d (fname='RB', units='s/m', & + avgflag='A', long_name='leaf boundary resistance', & + ptr_patch=this%rb1_patch, default='inactive') + + if (use_cn) then this%z0hv_patch(begp:endp) = spval call hist_addfld1d (fname='Z0HV', units='m', & avgflag='A', long_name='roughness length over vegetation, sensible heat', & diff --git a/src/biogeophys/LakeFluxesMod.F90 b/src/biogeophys/LakeFluxesMod.F90 index 3477c565d6..212d0ca7d1 100644 --- a/src/biogeophys/LakeFluxesMod.F90 +++ b/src/biogeophys/LakeFluxesMod.F90 @@ -269,6 +269,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, eflx_gnet => energyflux_inst%eflx_gnet_patch , & ! Output: [real(r8) (:) ] net heat flux into ground (W/m**2) taux => energyflux_inst%taux_patch , & ! Output: [real(r8) (:) ] wind (shear) stress: e-w (kg/m/s**2) tauy => energyflux_inst%tauy_patch , & ! Output: [real(r8) (:) ] wind (shear) stress: n-s (kg/m/s**2) + dhsdt_canopy => energyflux_inst%dhsdt_canopy_patch , & ! Output: [real(r8) (:) ] change in heat storage of stem (W/m**2) [+ to atm] ks => lakestate_inst%ks_col , & ! Output: [real(r8) (:) ] coefficient passed to LakeTemperature ws => lakestate_inst%ws_col , & ! Output: [real(r8) (:) ] surface friction velocity (m/s) @@ -377,6 +378,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, c = patch%column(p) g = patch%gridcell(p) + dhsdt_canopy(p) = 0.0_r8 nmozsgn(p) = 0 obuold(p) = 0._r8 displa(p) = 0._r8 diff --git a/src/biogeophys/SoilFluxesMod.F90 b/src/biogeophys/SoilFluxesMod.F90 index 303f4402d6..10082db373 100644 --- a/src/biogeophys/SoilFluxesMod.F90 +++ b/src/biogeophys/SoilFluxesMod.F90 @@ -88,6 +88,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, & !----------------------------------------------------------------------- associate( & + eflx_sh_stem => energyflux_inst%eflx_sh_stem_patch , & ! Output: [real(r8) (:) ] sensible heat flux from stems (W/m**2) [+ to atm] eflx_h2osfc_to_snow_col => energyflux_inst%eflx_h2osfc_to_snow_col , & ! Input: [real(r8) (:) ] col snow melt to h2osfc heat flux (W/m**2) forc_lwrad => atm2lnd_inst%forc_lwrad_downscaled_col , & ! Input: [real(r8) (:) ] downward infrared (longwave) radiation (W/m**2) @@ -315,6 +316,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, & ! Total fluxes (vegetation + ground) eflx_sh_tot(p) = eflx_sh_veg(p) + eflx_sh_grnd(p) + if (.not. lun%urbpoi(l)) eflx_sh_tot(p) = eflx_sh_tot(p) + eflx_sh_stem(p) qflx_evap_tot(p) = qflx_evap_veg(p) + qflx_evap_soi(p) eflx_lh_tot(p)= hvap*qflx_evap_veg(p) + htvp(c)*qflx_evap_soi(p) if (lun%itype(l) == istsoil .or. lun%itype(l) == istcrop) then diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index bb25b2bd53..050c6513b0 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -7,8 +7,8 @@ module TemperatureType use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use abortutils , only : endrun - use clm_varctl , only : use_cndv, iulog, use_luna, use_crop - use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevlak, nlevurb, nlevmaxurbgrnd + use clm_varctl , only : use_cndv, iulog, use_luna, use_crop, use_biomass_heat_storage + use clm_varpar , only : nlevsno, nlevgrnd, nlevlak, nlevurb, nlevmaxurbgrnd use clm_varcon , only : spval, ispval use GridcellType , only : grc use LandunitType , only : lun @@ -22,6 +22,7 @@ module TemperatureType type, public :: temperature_type ! Temperatures + real(r8), pointer :: t_stem_patch (:) ! patch stem temperatu\re (Kelvin) real(r8), pointer :: t_veg_patch (:) ! patch vegetation temperature (Kelvin) real(r8), pointer :: t_skin_patch (:) ! patch skin temperature (Kelvin) real(r8), pointer :: t_veg_day_patch (:) ! patch daytime accumulative vegetation temperature (Kelvinx*nsteps), LUNA specific, from midnight to current step @@ -192,6 +193,7 @@ subroutine InitAllocate(this, bounds) begg = bounds%begg; endg= bounds%endg ! Temperatures + allocate(this%t_stem_patch (begp:endp)) ; this%t_stem_patch (:) = nan allocate(this%t_veg_patch (begp:endp)) ; this%t_veg_patch (:) = nan allocate(this%t_skin_patch (begp:endp)) ; this%t_skin_patch (:) = nan if(use_luna) then @@ -388,6 +390,13 @@ subroutine InitHistory(this, bounds, is_simple_buildtemp, is_prog_buildtemp ) avgflag='A', long_name='Urban daily maximum of average 2-m temperature', & ptr_patch=this%t_ref2m_max_u_patch, set_nourb=spval, default='inactive') + if (use_biomass_heat_storage) then + this%t_stem_patch(begp:endp) = spval + call hist_addfld1d (fname='TSTEM', units='K', & + avgflag='A', long_name='stem temperature', & + ptr_patch=this%t_stem_patch, default='active') + endif + this%t_veg_patch(begp:endp) = spval call hist_addfld1d (fname='TV', units='K', & avgflag='A', long_name='vegetation temperature', & @@ -790,6 +799,8 @@ subroutine InitCold(this, bounds, & this%t_veg_patch(p) = 283._r8 end if + this%t_stem_patch(p) = this%t_veg_patch(p) + if (use_vancouver) then this%t_ref2m_patch(p) = 297.56 else if (use_mexicocity) then @@ -889,6 +900,11 @@ subroutine Restart(this, bounds, ncid, flag, is_simple_buildtemp, is_prog_buildt long_name='vegetation temperature', units='K', & interpinic_flag='interp', readvar=readvar, data=this%t_veg_patch) + call restartvar(ncid=ncid, flag=flag, varname='T_STEM', xtype=ncd_double, & + dim1name='pft', & + long_name='stem temperature', units='K', & + interpinic_flag='interp', readvar=readvar, data=this%t_stem_patch) + call restartvar(ncid=ncid, flag=flag, varname='TH2OSFC', xtype=ncd_double, & dim1name='column', & long_name='surface water temperature', units='K', & diff --git a/src/main/clm_varcon.F90 b/src/main/clm_varcon.F90 index bf4beaec7d..9f66d335ad 100644 --- a/src/main/clm_varcon.F90 +++ b/src/main/clm_varcon.F90 @@ -13,7 +13,7 @@ module clm_varcon SHR_CONST_RHOICE,SHR_CONST_TKFRZ,SHR_CONST_REARTH, & SHR_CONST_PDB, SHR_CONST_PI, SHR_CONST_CDAY, & SHR_CONST_RGAS, SHR_CONST_PSTD, & - SHR_CONST_MWDAIR, SHR_CONST_MWWV + SHR_CONST_MWDAIR, SHR_CONST_MWWV, SHR_CONST_CPFW use clm_varpar , only: numrad, nlevgrnd, nlevlak, nlevdecomp_full use clm_varpar , only: ngases use clm_varpar , only: nlayer @@ -81,6 +81,8 @@ module clm_varcon real(r8), public :: alpha_aero = 1.0_r8 ! constant for aerodynamic parameter weighting real(r8), public :: tlsai_crit = 2.0_r8 ! critical value of elai+esai for which aerodynamic parameters are maximum real(r8), public :: watmin = 0.01_r8 ! minimum soil moisture (mm) + real(r8), public :: c_water = SHR_CONST_CPFW ! specific heat of water [J/kg/K] + real(r8), public :: c_dry_biomass = 1400_r8 ! specific heat of dry biomass real(r8), public :: re = SHR_CONST_REARTH*0.001_r8 ! radius of earth (km) @@ -117,11 +119,12 @@ module clm_varcon real(r8), public :: thk_bedrock = 3.0_r8 ! thermal conductivity of 'typical' saturated granitic rock ! (Clauser and Huenges, 1995)(W/m/K) - real(r8), public :: csol_bedrock = 2.0e6_r8 ! vol. heat capacity of granite/sandstone J/(m3 K)(Shabbir, 2000) !scs + real(r8), public :: csol_bedrock = 2.0e6_r8 ! vol. heat capacity of granite/sandstone J/(m3 K)(Shabbir, 2000) real(r8), public, parameter :: zmin_bedrock = 0.4_r8 ! minimum soil depth [m] 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) + !!! C13 real(r8), public, parameter :: preind_atm_del13c = -6.0 ! preindustrial value for atmospheric del13C real(r8), public, parameter :: preind_atm_ratio = SHR_CONST_PDB + (preind_atm_del13c * SHR_CONST_PDB)/1000.0 ! 13C/12C diff --git a/src/main/clm_varctl.F90 b/src/main/clm_varctl.F90 index b70827a084..aadb32ce09 100644 --- a/src/main/clm_varctl.F90 +++ b/src/main/clm_varctl.F90 @@ -278,6 +278,12 @@ module clm_varctl logical, public :: use_lai_streams = .false. ! true => use lai streams in SatellitePhenologyMod.F90 + !---------------------------------------------------------- + ! biomass heat storage switch + !---------------------------------------------------------- + + logical, public :: use_biomass_heat_storage = .false. ! true => include biomass heat storage in canopy energy budget + !---------------------------------------------------------- ! bedrock / soil depth switch !---------------------------------------------------------- diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index ef06a19285..55b4cc338d 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -241,6 +241,8 @@ subroutine control_init(dtime) namelist /clm_inparm/ use_bedrock + namelist /clm_inparm/ use_biomass_heat_storage + namelist /clm_inparm/ use_hydrstress namelist /clm_inparm/ use_dynroot @@ -736,6 +738,8 @@ subroutine control_spmd() call mpi_bcast (use_bedrock, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (use_biomass_heat_storage, 1, MPI_LOGICAL, 0, mpicom, ier) + call mpi_bcast (use_hydrstress, 1, MPI_LOGICAL, 0, mpicom, ier) call mpi_bcast (use_dynroot, 1, MPI_LOGICAL, 0, mpicom, ier) diff --git a/src/main/pftconMod.F90 b/src/main/pftconMod.F90 index 702a03b511..88e5965051 100644 --- a/src/main/pftconMod.F90 +++ b/src/main/pftconMod.F90 @@ -26,7 +26,6 @@ module pftconMod integer, public :: nbrdlf_dcd_trp_tree ! value for Broadleaf deciduous tropical tree integer, public :: nbrdlf_dcd_tmp_tree ! value for Broadleaf deciduous temperate tree integer, public :: nbrdlf_dcd_brl_tree ! value for Broadleaf deciduous boreal tree - integer, public :: ntree ! value for last type of tree integer, public :: nbrdlf_evr_shrub ! value for Broadleaf evergreen shrub integer, public :: nbrdlf_dcd_tmp_shrub ! value for Broadleaf deciduous temperate shrub integer, public :: nbrdlf_dcd_brl_shrub ! value for Broadleaf deciduous boreal shrub @@ -110,7 +109,9 @@ module pftconMod type, public :: pftcon_type integer , allocatable :: noveg (:) ! value for not vegetated - integer , allocatable :: tree (:) ! tree or not? + logical , allocatable :: is_tree (:) ! tree or not? + logical , allocatable :: is_shrub (:) ! shrub or not? + logical , allocatable :: is_grass (:) ! grass or not? real(r8), allocatable :: dleaf (:) ! characteristic leaf dimension (m) real(r8), allocatable :: c3psn (:) ! photosynthetic pathway: 0. = c4, 1. = c3 @@ -147,6 +148,13 @@ module pftconMod real(r8), allocatable :: root_radius (:) ! root radius (m) real(r8), allocatable :: root_density (:) ! root density (gC/m3) + real(r8), allocatable :: dbh (:) ! diameter at breast height (m) + real(r8), allocatable :: fbw (:) ! fraction of biomass that is water + real(r8), allocatable :: nstem (:) ! stem density (#/m2) + real(r8), allocatable :: taper (:) ! tapering ratio of height:radius_breast_height + real(r8), allocatable :: rstem_per_dbh (:) ! stem resistance per dbh (s/m/m) + real(r8), allocatable :: wood_density (:) ! wood density (kg/m3) + ! crop ! These arrays give information about the merge of unused crop types to the types CLM @@ -336,8 +344,10 @@ subroutine InitAllocate (this) class(pftcon_type) :: this !----------------------------------------------------------------------- - allocate( this%noveg (0:mxpft)); this%noveg (:) =huge(1) - allocate( this%tree (0:mxpft)); this%tree (:) =huge(1) + allocate( this%noveg (0:mxpft)); this%noveg (:) = huge(1) + allocate( this%is_tree (0:mxpft)); this%is_tree (:) = .false. + allocate( this%is_shrub (0:mxpft)); this%is_shrub (:) = .false. + allocate( this%is_grass (0:mxpft)); this%is_grass (:) = .false. allocate( this%dleaf (0:mxpft) ) allocate( this%c3psn (0:mxpft) ) @@ -467,6 +477,12 @@ subroutine InitAllocate (this) allocate( this%fun_cn_flex_c (0:mxpft) ) allocate( this%FUN_fracfixers(0:mxpft) ) + allocate( this%dbh (0:mxpft) ) + allocate( this%fbw (0:mxpft) ) + allocate( this%nstem (0:mxpft) ) + allocate( this%taper (0:mxpft) ) + allocate( this%rstem_per_dbh (0:mxpft) ) + allocate( this%wood_density (0:mxpft) ) end subroutine InitAllocate @@ -481,7 +497,7 @@ subroutine InitRead(this) use fileutils , only : getfil use ncdio_pio , only : ncd_io, ncd_pio_closefile, ncd_pio_openfile, file_desc_t use ncdio_pio , only : ncd_inqdid, ncd_inqdlen - use clm_varctl , only : paramfile, use_fates, use_flexibleCN, use_dynroot + use clm_varctl , only : paramfile, use_fates, use_flexibleCN, use_dynroot, use_biomass_heat_storage use spmdMod , only : masterproc use CLMFatesParamInterfaceMod, only : FatesReadPFTs ! @@ -986,13 +1002,8 @@ subroutine InitRead(this) this%dwood(m) = dwood this%root_radius(m) = root_radius this%root_density(m) = root_density - - if (m <= ntree) then - this%tree(m) = 1 - else - this%tree(m) = 0 - end if end do + ! ! clm 5 nitrogen variables ! @@ -1018,6 +1029,33 @@ subroutine InitRead(this) if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) end if + call ncd_io('nstem',this%nstem, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + call ncd_io('taper',this%taper, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + ! + ! Biomass heat storage variables + ! + if (use_biomass_heat_storage ) then + ! + ! These variables are used for stem biomass and only for tree and shrub + ! (They are effectively unused for other veg types) + ! + call ncd_io('dbh',this%dbh, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + call ncd_io('fbw',this%fbw, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + call ncd_io('rstem',this%rstem_per_dbh, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + call ncd_io('wood_density',this%wood_density, 'read', ncid, readvar=readv) + if ( .not. readv ) call endrun(msg=' ERROR: error in reading in pft data'//errMsg(sourcefile, __LINE__)) + else + this%dbh = 0.0_r8 + this%fbw = 0.0_r8 + this%rstem_per_dbh = 0.0_r8 + this%wood_density = 0.0_r8 + end if + call ncd_pio_closefile(ncid) call FatesReadPFTs() @@ -1112,13 +1150,35 @@ subroutine InitRead(this) if ( trim(pftname(i)) == 'irrigated_tropical_soybean' ) nirrig_trp_soybean = i end do - ntree = nbrdlf_dcd_brl_tree ! value for last type of tree npcropmin = ntmp_corn ! first prognostic crop npcropmax = mxpft ! last prognostic crop in list call this%set_is_pft_known_to_model() call this%set_num_cfts_known_to_model() + ! Set vegetation family identifier (tree/shrub/grass) + do m = 0,mxpft + if (m == ndllf_evr_tmp_tree .or. m == ndllf_evr_brl_tree & + .or. m == ndllf_dcd_brl_tree .or. m == nbrdlf_evr_trp_tree & + .or. m == nbrdlf_evr_tmp_tree .or. m == nbrdlf_dcd_trp_tree & + .or. m == nbrdlf_dcd_tmp_tree .or. m == nbrdlf_dcd_brl_tree) then + this%is_tree(m) = .true. + else + this%is_tree(m) = .false. + endif + if(m == nbrdlf_evr_shrub .or. m == nbrdlf_dcd_tmp_shrub .or. m == nbrdlf_dcd_brl_shrub) then + this%is_shrub(m) = .true. + else + this%is_shrub(m) = .false. + endif + if(m == nc3_arctic_grass .or. m == nc3_nonarctic_grass .or. m == nc4_grass) then + this%is_grass(m) = .true. + else + this%is_grass(m) = .false. + endif + + end do + if (use_cndv) then this%fcur(:) = this%fcurdv(:) end if @@ -1277,7 +1337,9 @@ subroutine Clean(this) !----------------------------------------------------------------------- deallocate( this%noveg) - deallocate( this%tree) + deallocate( this%is_tree) + deallocate( this%is_shrub) + deallocate( this%is_grass) deallocate( this%dleaf) deallocate( this%c3psn) @@ -1407,6 +1469,12 @@ subroutine Clean(this) deallocate( this%fun_cn_flex_c) deallocate( this%FUN_fracfixers) + deallocate( this%dbh) + deallocate( this%fbw) + deallocate( this%nstem) + deallocate( this%rstem_per_dbh) + deallocate( this%wood_density) + deallocate( this%taper) end subroutine Clean end module pftconMod diff --git a/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 b/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 index 98951f9b56..8305fcefe6 100644 --- a/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemDecompCascadeCNMod.F90 @@ -668,7 +668,9 @@ subroutine decomp_rate_constants_cn(bounds, & end if ! The following code implements the acceleration part of the AD spinup - ! algorithm, by multiplying all of the SOM decomposition base rates by 10.0. + ! algorithm, by multiplying all of the SOM decomposition base rates by + ! spinup_vector, scalar between 1 and 70X, defined as a constant for each + ! pool here if ( spinup_state .eq. 1 ) then k_s1 = k_s1 * params_inst%spinup_vector(1)