diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 3ce5fe24a5..698f2396ab 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -3532,11 +3532,12 @@ sub setup_logic_methane { add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'finundation_method', 'use_cn'=>$nl_flags->{'use_cn'}, 'use_fates'=>$nl_flags->{'use_fates'} ); my $finundation_method = remove_leading_and_trailing_quotes($nl->get_value('finundation_method' )); - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_ch4finundated', - 'finundation_method'=>$finundation_method); - if ($opts->{'driver'} eq "nuopc" ) { - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_ch4finundated', - 'finundation_method'=>$finundation_method); + # prognostic inundation does not require an input stream; other methods do + if($finundation_method ne 'h2osfc') { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_fldfilename_ch4finundated', + 'finundation_method'=>$finundation_method); + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'stream_meshfile_ch4finundated', + 'finundation_method'=>$finundation_method); } add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_aereoxid_prog', 'use_cn'=>$nl_flags->{'use_cn'}, 'use_fates'=>$nl_flags->{'use_fates'} ); diff --git a/doc/ChangeLog b/doc/ChangeLog index 6504fe8948..9d3c8474d7 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,131 @@ =============================================================== +Tag name: ctsm5.4.056 +Originator(s): slevis (Samuel Levis,UCAR/TSS,slevis@ucar.edu) +Date: Wed Sep 9 04:48:03 PM MDT 2026 +One-line Summary: Merge the wetlands branch to master + +Purpose and description of changes +---------------------------------- + + The wetlands branch has these tags... + branch_tags/wetlands.n01.ctsm5.4.029 + branch_tags/wetlands.n02.ctsm5.4.029 + branch_tags/wetlands.n03.ctsm5.4.029 + branch_tags/wetlands.n03.ctsm5.4.044 + branch_tags/wetlands.n04.ctsm5.4.044 + branch_tags/wetlands.n04.ctsm5.4.054 <-- PR #4204 started with this tag, I merged ctsm5.4.055 to it, then I merged to master (explained next) + branch_tags/wetlands.n05.ctsm5.4.054 + We wanted wetlands.n04.ctsm5.4.055 merged to master, and we accomplished this with PR #4204, without merging or tagging anything new in the wetlands branch. + + Main updates coming in with the wetlands branch: + 1. Changes to snow-covered fraction and inundated fraction where the two overlap. + 2. Changes to the methane parameterization. + Best evaluation of updates (1) and (2): https://webext.cgd.ucar.edu/I20TR/ctsm5.4.024_f19_wetlandCH4_151_HIST/ctsm5.4.024_f19_wetlandCH4_151_HIST_1850_2023_vs_ctsm5.4.024_f19_150_HIST_1850_2023/website + 3. Changes to phenology to improve high-lat survivability by allowing plants to green up after the summer solstice. + Best evaluation of (3): https://github.com/NCAR/LMWG_dev/issues/174 + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- +[Remove entire section if none of the boxes are checked.] + +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] clm6_0 + +[X] clm5_0 + +[X] ctsm5_0-nwp + +Bugs fixed +---------- +[Remove any lines that don't apply. Remove entire section if nothing applies.] + +List of CTSM issues fixed (include CTSM Issue # and description) [one per line]: + List of PRs available below + +Notes of particular relevance for users +--------------------------------------- +[Remove any lines that don't apply. Remove entire section if nothing applies.] + +Caveats for users (e.g., need to interpolate initial conditions): + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): + +Changes made to namelist defaults (e.g., changed parameter values): + +Changes to the datasets (e.g., streams, surface or initial files): + +Changes to the parameter file (output of tools/param_utils/compare_paramfiles): + +Changes to documentation: + Updating doc/source/tech_note/Vegetation_Phenology_Turnover/CLM50_Tech_Note_Vegetation_Phenology_Turnover.rst + Opened issue for further updates: https://github.com/ESCOMP/CTSM/issues/4205 + +Contributors: + swensosc, wwieder, lifang0209 + +Notes of particular relevance for developers: +--------------------------------------------- +NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide +[Remove any lines that don't apply. Remove entire section if nothing applies.] + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + +Changes to tests or testing: + Updating src/dyn_subgrid/test/dynColumnStateUpdater_test/test_column_state_updater.pf + +Testing summary: +---------------- + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + derecho - OK (2 expected fails) + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + derecho ----- PEND (rerunning 3 tests) + izumi ------- OK + + any other testing (give details below): + + ctsm_sci + derecho ---- OK + + +Answer changes +-------------- + +Changes answers relative to baseline: Yes + + Summarize any changes to answers, i.e., + - what code configurations: All + - what platforms/compilers: All + - nature of change: larger than roundoff, subtly new climate + + If this tag changes climate list the run(s) done to evaluate the new + climate. Preferably in https://github.com/NCAR/LMWG_dev (or give details below): + - Best evaluation of updates (1) and (2) listed in the changes section above: + https://webext.cgd.ucar.edu/I20TR/ctsm5.4.024_f19_wetlandCH4_151_HIST/ctsm5.4.024_f19_wetlandCH4_151_HIST_1850_2023_vs_ctsm5.4.024_f19_150_HIST_1850_2023/website + - Best evaluation of update (3): + https://github.com/NCAR/LMWG_dev/issues/174 + +Other details +------------- +Pull Requests that document the changes (include PR ids): + https://github.com/ESCOMP/ctsm/pull/3893 from swensosc/wetland_update + https://github.com/ESCOMP/ctsm/pull/3934 from wwieder/wetlands_fire + https://github.com/ESCOMP/ctsm/pull/3933 from swensosc/wetland_update + https://github.com/ESCOMP/ctsm/pull/4091 merge master to wetlands + https://github.com/ESCOMP/ctsm/pull/4108 from swensosc/seasonal_decid_onset + https://github.com/ESCOMP/ctsm/pull/4191 merge master to wetlands + https://github.com/ESCOMP/ctsm/pull/4204 this master tag ctsm5.4.056 + +=============================================================== +=============================================================== Tag name: ctsm5.4.055 Originator(s): sacks (Bill Sacks) Date: Tue Sep 8 16:01:56 MDT 2026 diff --git a/doc/ChangeSum b/doc/ChangeSum index 520580bd66..9f7548d1ca 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.4.056 multiple 09/10/2026 Merge the wetlands branch to master ctsm5.4.055 sacks 09/08/2026 Add temporal smoothing of DYNBAL conservation correction fluxes ctsm5.4.054 erik 08/31/2026 Update submodules to the latest needed for cesm3_0_alpha10a ctsm5.4.053 oleson 08/18/2026 Resolve problems with running PLUMBER2 sites diff --git a/doc/source/tech_note/Vegetation_Phenology_Turnover/CLM50_Tech_Note_Vegetation_Phenology_Turnover.rst b/doc/source/tech_note/Vegetation_Phenology_Turnover/CLM50_Tech_Note_Vegetation_Phenology_Turnover.rst index 2178baa32f..b5babb2e08 100644 --- a/doc/source/tech_note/Vegetation_Phenology_Turnover/CLM50_Tech_Note_Vegetation_Phenology_Turnover.rst +++ b/doc/source/tech_note/Vegetation_Phenology_Turnover/CLM50_Tech_Note_Vegetation_Phenology_Turnover.rst @@ -412,7 +412,9 @@ and the associated nitrogen fluxes are: where :math:`{f}_{stor,xfer}` is the fraction of current storage pool moved into the transfer pool for display over the incipient onset period. This fraction is set to 0.5, based on the observation that seasonal deciduous trees are capable of replacing their canopies from storage reserves in the event of a severe early-season disturbance such as frost damage or defoliation due to insect herbivory. -If the onset criteria are not met before the summer solstice, then :math:`{GDD}_{sum}` is set to 0.0 and the growing-degree-day accumulation will not start again until the following winter solstice. This mechanism prevents the initiation of very short growing seasons late in the summer in cold climates. The onset counter is decremented on each time step after initiation of the onset period, until it reaches zero, signaling the end of the onset period: +The *min_critical_daylength_onset* parameter specifies a daylength below which onset may not occur. If the onset criterion (:math:`{GDD}_{sum} > {GDD}_{sum\_crit}`) is not attained before the daylength becomes less than this parameter in autumn, then :math:`{GDD}_{sum}` is set to 0 and the growing-degree-day accumulation will not start again until the following winter solstice. In cold climates, this mechanism ensures that if the onset criteria are not met, the accumulation does not continue into the following year. + +The onset counter is decremented on each time step after initiation of the onset period, until it reaches zero, signaling the end of the onset period: .. math:: :label: t_onset_decrement diff --git a/src/biogeochem/CNFireLi2024Mod.F90 b/src/biogeochem/CNFireLi2024Mod.F90 index d6c15f25ca..1aa5f7b1f3 100644 --- a/src/biogeochem/CNFireLi2024Mod.F90 +++ b/src/biogeochem/CNFireLi2024Mod.F90 @@ -217,7 +217,7 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ lgdp_col => cnveg_state_inst%lgdp_col , & ! Output: [real(r8) (:) ] gdp limitation factor for nfire lgdp1_col => cnveg_state_inst%lgdp1_col , & ! Output: [real(r8) (:) ] gdp limitation factor for baf per fire lpop_col => cnveg_state_inst%lpop_col , & ! Output: [real(r8) (:) ] pop limitation factor for baf per fire - lfwt => cnveg_state_inst%lfwt_col , & ! Output: [real(r8) (:) ] fractional coverage of non-crop and non-bare-soil Patches + tgswt => cnveg_state_inst%tgswt_grc , & ! Output: [real(r8) (:) ] gridcell fractional coverage of tree, grass and shrub PFTs (0-1) trotr1_col => cnveg_state_inst%trotr1_col , & ! Output: [real(r8) (:) ] patch weight of BET on the column (0-1) trotr2_col => cnveg_state_inst%trotr2_col , & ! Output: [real(r8) (:) ] patch weight of BDT on the column (0-1) dtrotr_col => cnveg_state_inst%dtrotr_col , & ! Output: [real(r8) (:) ] decreased frac. coverage of BET+BDT on grid for dt @@ -288,24 +288,26 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ return end if ! - ! Calculate fraction of crop (cropf_col) and non-crop and non-bare-soil - ! vegetation (lfwt) in vegetated column + ! Calculate fraction of crop in vegetation column (cropf_col) + ! and fraction of tree, grass, and shrub in grid cell (tgswt) ! do fc = 1,num_soilc c = filter_soilc(fc) + g = col%gridcell(c) cropf_col(c) = 0._r8 - lfwt(c) = 0._r8 + tgswt(g) = 0._r8 end do do fp = 1, num_soilp p = filter_soilp(fp) c = patch%column(p) + g = col%gridcell(c) ! For crop veg types if( patch%itype(p) > nc4_grass )then cropf_col(c) = cropf_col(c) + patch%wtcol(p) end if - ! Exclude crops and bare soil + ! For tree, grass, and shrub if (patch%itype(p) <= nc4_grass .and. patch%itype(p) >= ndllf_evr_tmp_tree) then - lfwt(c) = lfwt(c) + patch%wtgcell(p) + tgswt(g) = tgswt(g) + patch%wtgcell(p) end if end do ! @@ -543,8 +545,9 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ do fc = 1, num_soilc c = filter_soilc(fc) g= col%gridcell(c) + tgswt(g) = max(0._r8, min(1._r8, tgswt(g))) if(grc%latdeg(g) < cnfire_const%borealat )then - if ((trotr1_col(c)+trotr2_col(c))*col%wtgcell(c)<=0.8_r8.and.trotr1_col(c)+trotr2_col(c)>0.0_r8) then + if ((trotr1_col(c)+trotr2_col(c))*tgswt(g)<=0.8_r8.and.trotr1_col(c)+trotr2_col(c)>0.0_r8) then baf_peatf(c) = non_boreal_peatfire_c/secsphr*max(0._r8, & min(1._r8,(1._r8-prec30_col(c)*secspday/nonborpeat_fire_precip_denom)))*peatf_lf(c) else @@ -609,12 +612,12 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ ig = (lh+this%forc_lnfm(g)/(5.16_r8+2.16_r8* & cos(SHR_CONST_PI/180._r8*3*min(60._r8,abs(grc%latdeg(g)))))* & cnfire_params%ignition_efficiency)*(1._r8-fs)* & - (lfwt(c)**0.5) + (tgswt(g)**0.5) else ig = this%forc_lnfm(g)/(5.16_r8+2.16_r8* & cos(SHR_CONST_PI/180._r8*3*min(60._r8,abs(grc%latdeg(g)))))* & cnfire_params%ignition_efficiency*(1._r8-fs)* & - (lfwt(c)**0.5) + (tgswt(g)**0.5) end if ! Reduce burnability at high elevations @@ -627,7 +630,7 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ nfire(c) = ig/secsphr*fb*fire_m*lgdp_col(c) * topoi !fire counts/km2/sec Lb_lf = 1._r8+10._r8*(1._r8-EXP(-0.06_r8*forc_wind(g))) spread_m = fire_m**0.5_r8 - fd_col(c) = (lfwt(c)*lgdp1_col(c)*lpop_col(c))**0.5_r8 * fd_col(c) + fd_col(c) = (tgswt(g)*lgdp1_col(c)*lpop_col(c))**0.5_r8 * fd_col(c) farea_burned(c) = min(1._r8,(cnfire_const%g0*spread_m*fsr_col(c)* & fd_col(c)/1000._r8)**2*nfire(c)*SHR_CONST_PI*Lb_lf+ & baf_crop(c)+baf_peatf(c)) ! fraction (0-1) per sec diff --git a/src/biogeochem/CNPhenologyMod.F90 b/src/biogeochem/CNPhenologyMod.F90 index 71ec677961..28693464be 100644 --- a/src/biogeochem/CNPhenologyMod.F90 +++ b/src/biogeochem/CNPhenologyMod.F90 @@ -1289,10 +1289,17 @@ function SeasonalDecidOnset( onset_gdd, onset_gddflag, soilt, soila10, t_a5min, ! Test to turn off growing degree-day sum, if on. ! This test resets the growing degree day sum if it gets past ! the summer solstice without reaching the threshold value. - ! In that case, it will take until the next winter solstice - ! before the growing degree-day summation starts again. - - if (onset_gddflag == 1._r8 .and. ws_flag == 0._r8) then + ! After the summer solstice, a daylength criterion is used + ! to trigger the reset: if daylength is less than the + ! parameter value, gdd is reset. This criterion was added to + ! allow for situations in the high latitudes when using only + ! the summer solstice as a trigger caused some vegetation to never + ! achieve onset. + ! If the conditional is met, it will take until the next winter + ! solstice before the growing degree-day summation starts again. + + if (onset_gddflag == 1._r8 .and. ws_flag == 0._r8.and. & + dayl crit_onset_gdd .and. season_decid_temperate == 1) then do_onset = .true. - ! Note: The check "dayl>min_critical_daylength_onset" in the if - ! statement was added because for some coastal - ! points the other triggers could allow onset in January/February - ! which isn't sustainable and is a degenerate case. To prevent this - ! condition was added, but now the other conditions aren't triggered - ! until much later so it's value just needs to be high enough to prevent - ! the degenerate case of happening too early, and low enough that it - ! doesn't restrict onset. As such the value of this parameter shouldn't - ! matter for reasonable values between the two degenerate cases. + ! Note: The "dayl>min_critical_daylength_onset" criterion + ! was added because for some coastal points, onset was triggered + ! too early (e.g January/February). This criterion ensures that + ! onset does not occur too early, and at the same time (because the + ! criterion is symmetric about the solstices), it does not allow + ! onset when the day length becomes less than + ! min_critical_daylength_onset after the summer solstice. + ! The value of min_critical_daylength_onset needs to be high + ! enough to prevent onset happening too early, and low enough + ! that it doesn't restrict onset. else if (season_decid_temperate == 0 .and. onset_gddflag == 1.0_r8 .and. & soila10 > SHR_CONST_TKFRZ .and. & - t_a5min > SHR_CONST_TKFRZ .and. ws_flag==1.0_r8 .and. & + t_a5min > SHR_CONST_TKFRZ .and. & dayl>min_critical_daylength_onset .and. & snow_5day ch4_inst%ch4_oxid_depth_lake_col , & ! Output: [real(r8) (:,:) ] CH4 consumption rate via oxidation in each soil layer (mol/m3/s) (nlevsoi) ch4_prod_depth_sat => ch4_inst%ch4_prod_depth_sat_col , & ! Output: [real(r8) (:,:) ] production of CH4 in each soil layer (nlevsoi) (mol/m3/s) ch4_prod_depth_unsat => ch4_inst%ch4_prod_depth_unsat_col , & ! Output: [real(r8) (:,:) ] production of CH4 in each soil layer (nlevsoi) (mol/m3/s) + ch4_prod_tot => ch4_inst%ch4_prod_tot_col , & ! Output: [real(r8) (:,:) ] col CH4 production flux (mol/m2/s) + ch4_oxid_tot => ch4_inst%ch4_oxid_tot_col , & ! Output: [real(r8) (:,:) ] col CH4 oxidation flux (mol/m2/s) ch4_prod_depth_lake => ch4_inst%ch4_prod_depth_lake_col , & ! Output: [real(r8) (:,:) ] production of CH4 in each soil layer (nlevsoi) (mol/m3/s) lake_soilc => ch4_inst%lake_soilc_col , & ! Output: [real(r8) (:,:) ] total soil organic matter found in level (g C / m^3) (nlevsoi) conc_ch4_sat => ch4_inst%conc_ch4_sat_col , & ! Output: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) conc_ch4_unsat => ch4_inst%conc_ch4_unsat_col , & ! Output: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) + + surface_layer_conc_ch4_sat => ch4_inst%conc_ch4_surface_layer_sat_col , & ! Input: [real(r8) (:) ] CH4 conc surface water layer (mol/m3) + surface_layer_conc_ch4_unsat => ch4_inst%conc_ch4_surface_layer_unsat_col , & ! Input: [real(r8) (:) ] CH4 conc surface soil layer (mol/m3) + surface_layer_thickness_sat => ch4_inst%surface_layer_thickness_sat_col , & ! Input: [real(r8) (:) ] surface water layer thickness (m) + surface_layer_thickness_unsat => ch4_inst%surface_layer_thickness_unsat_col , & ! Input: [real(r8) (:) ] surface soil layer thickness (m) + conc_ch4_lake => ch4_inst%conc_ch4_lake_col , & ! Output: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) conc_o2_lake => ch4_inst%conc_o2_lake_col , & ! Output: [real(r8) (:,:) ] O2 conc in each soil layer (mol/m3) (nlevsoi) ch4_dfsat_flux => ch4_inst%ch4_dfsat_flux_col , & ! Output: [real(r8) (:) ] CH4 flux to atm due to decreasing finundated (kg C/m^2/s) [+] @@ -1881,9 +2004,7 @@ subroutine ch4 (bounds, num_soilc, filter_soilc, num_lakec, filter_lakec, & end do - ! Check to see if finundated changed since the last timestep. If it increased, then reduce conc_ch4_sat - ! proportionally. If it decreased, then add flux to atm. - + ! Check to see if finundated changed since the last timestep and adjust concentrations to conserve vertically integrated methane do j=1,nlevsoi do fc = 1, num_soilc c = filter_soilc(fc) @@ -1894,18 +2015,39 @@ subroutine ch4 (bounds, num_soilc, filter_soilc, num_lakec, filter_lakec, & g = col%gridcell(c) if (.not. ch4_first_time_grc(g)) then - if (finundated(c) > fsat_bef(c)) then !Reduce conc_ch4_sat - dfsat = finundated(c) - fsat_bef(c) + dfsat = finundated(c) - fsat_bef(c) + if (dfsat > 0._r8) then ! adjust conc_ch4_sat conc_ch4_sat(c,j) = (fsat_bef(c)*conc_ch4_sat(c,j) + dfsat*conc_ch4_unsat(c,j)) / finundated(c) - else if (finundated(c) < fsat_bef(c)) then - ch4_dfsat_flux(c) = ch4_dfsat_flux(c) + & - (fsat_bef(c) - finundated(c))*(conc_ch4_sat(c,j) - conc_ch4_unsat(c,j)) * & - dz(c,j) / dtime * catomw / 1000._r8 ! mol --> kg + else if (dfsat < 0._r8) then ! adjust conc_ch4_unsat + if (finundated(c) < 1._r8) then + conc_ch4_unsat(c,j) = ((1._r8-fsat_bef(c))*conc_ch4_unsat(c,j) - dfsat*conc_ch4_sat(c,j)) / (1._r8 - finundated(c)) + endif end if end if end do end do - + ! Account for surface layer + do fc = 1, num_soilc + c = filter_soilc(fc) + g = col%gridcell(c) + if (.not. ch4_first_time_grc(g)) then + dfsat = finundated(c) - fsat_bef(c) + if (dfsat > 0._r8) then ! adjust conc_ch4_sat + surface_layer_conc_ch4_sat(c) = & + (fsat_bef(c)*surface_layer_conc_ch4_sat(c) & + + dfsat*surface_layer_conc_ch4_unsat(c)*(surface_layer_thickness_unsat(c)/surface_layer_thickness_sat(c))) & + / finundated(c) + else if (dfsat < 0._r8) then ! adjust conc_ch4_unsat + if (finundated(c) < 1._r8) then + surface_layer_conc_ch4_unsat(c) = & + ((1._r8-fsat_bef(c))*surface_layer_conc_ch4_unsat(c) & + - dfsat*surface_layer_conc_ch4_sat(c)*(surface_layer_thickness_sat(c)/surface_layer_thickness_unsat(c))) & + / (1._r8 - finundated(c)) + endif + end if + end if + end do + !!!! Begin biochemistry ! First for soil @@ -1989,7 +2131,7 @@ subroutine ch4 (bounds, num_soilc, filter_soilc, num_lakec, filter_lakec, & call get_jwt (bounds, num_soilc, filter_soilc, jwt(begc:endc), & soilstate_inst, waterstatebulk_inst, temperature_inst) - + do fc = 1, num_soilc c = filter_soilc(fc) zwt_ch4_unsat(c) = zi(c,jwt(c)) @@ -2558,6 +2700,7 @@ subroutine ch4_prod (bounds, num_methc, filter_methc, num_methp, & else if ( (somhr(c) + lithr(c)) > 0._r8) then partition_z = hr_vr(c,j) * dz(c,j) / (somhr(c) + lithr(c)) + else partition_z = 1._r8 end if @@ -2628,6 +2771,7 @@ subroutine ch4_prod (bounds, num_methc, filter_methc, num_methp, & ! Competition will be done in ch4_oxid o2_decomp_depth(c,j) = base_decomp * partition_z / dz (c,j) + if (anoxia) then ! Divide off o_scalar to use potential O2-unlimited HR to represent aerobe demand for oxygen competition if (.not. lake .and. j > nlevdecomp) then @@ -2641,6 +2785,7 @@ subroutine ch4_prod (bounds, num_methc, filter_methc, num_methp, & end if end if ! anoxia + ! Add root respiration if (.not. lake) then o2_decomp_depth(c,j) = o2_decomp_depth(c,j) + rr_vr(c,j)/catomw/dz(c,j) ! mol/m^3/s @@ -2788,7 +2933,7 @@ subroutine ch4_oxid (bounds, & k_m_eff = k_m_unsat vmax_eff = vmax_oxid_unsat end if - + porevol = max(watsat(c,j) - h2osoi_vol(c,j), 0._r8) h2osoi_vol_min = min(watsat(c,j), h2osoi_vol(c,j)) if (j <= jwt(c) .and. smp_l(c,j) < 0._r8) then @@ -2821,7 +2966,7 @@ subroutine ch4_oxid (bounds, & ch4_oxid_depth(c,j) = oxid_a o2_oxid_depth(c,j) = ch4_oxid_depth(c,j) * 2._r8 - + end do end do @@ -3117,6 +3262,7 @@ subroutine SiteOxAere(is_vegetated, & if (transpirationloss .and. is_vegetated) then ! Calculate water concentration h2osoi_vol_min = min(watsat(j), h2osoi_vol(j)) + k_h_inv = exp(-c_h_inv(1) * (1._r8 / t_soisno(j) - 1._r8 / kh_tbase) + log (kh_theta(1))) k_h_cc = t_soisno(j) / k_h_inv * rgasLatm conc_ch4_wat = conc_ch4(j) / ( (watsat(j)-h2osoi_vol_min)/k_h_cc + h2osoi_vol_min) @@ -3170,7 +3316,7 @@ subroutine SiteOxAere(is_vegetated, & aerecond = area_tiller * rootfr(j) * diffus_aere / (z(j)*params_inst%rob) ! Add in boundary layer resistance aerecond = 1._r8 / (1._r8/(aerecond+smallnumber) + 1._r8/(grnd_ch4_cond+smallnumber)) - + aere(j) = aerecond * (conc_ch4(j)/watsat(j)/k_h_cc - c_atm(1)) / dz(j) ![mol/m3-total/s] !ZS: Added watsat & Henry's const. aere(j) = max(aere(j), 0._r8) ! prevent backwards diffusion @@ -3182,7 +3328,9 @@ subroutine SiteOxAere(is_vegetated, & aerecond = area_tiller * rootfr(j) * oxdiffus / (z(j)*params_inst%rob) aerecond = 1._r8 / (1._r8/(aerecond+smallnumber) + 1._r8/(grnd_ch4_cond+smallnumber)) oxaere(j) = -aerecond *(conc_o2(j)/watsat(j)/k_h_cc - c_atm(2)) / dz(j) ![mol/m3-total/s] + oxaere(j) = max(oxaere(j), 0._r8) + ! Diffusion in is positive; prevent backwards diffusion if ( .not. use_aereoxid_prog ) then ! fixed aere oxid proportion; will be done in ch4_tran oxaere(j) = 0._r8 @@ -3197,7 +3345,6 @@ subroutine SiteOxAere(is_vegetated, & return end subroutine SiteOxAere - !----------------------------------------------------------------------- subroutine ch4_ebul (bounds, & num_methc, filter_methc, & @@ -3341,6 +3488,104 @@ subroutine ch4_ebul (bounds, & end subroutine ch4_ebul + !----------------------------------------------------------------------- + subroutine surface_layer_oxid(watsat, sucsat, bsw, h2osoi_vol, & + smp_l, t_soisno, conc_ch4, conc_o2, ch4_oxid, o2_oxid) + ! + ! !DESCRIPTION: + ! Calculate oxidation in virtual unsaturated layer + ! Oxidation is based on double Michaelis-Mentin kinetics, and is adjusted for low soil moisture. + ! Oxidation will be limited by available oxygen in ch4_tran. + ! + + ! !USES: + + ! + ! Arguments (in) + real(r8), intent(in) :: watsat ! volumetric soil water at saturation (porosity) + real(r8), intent(in) :: sucsat ! minimum soil suction (mm) + real(r8), intent(in) :: bsw ! Clapp and Hornberger "b" + real(r8), intent(in) :: h2osoi_vol ! volumetric soil water (m3/m3) + real(r8), intent(in) :: smp_l ! soil matric potential (mm) + real(r8), intent(in) :: t_soisno ! soil temperature (Kelvin) + real(r8), intent(in) :: conc_ch4 ! CH4 conc in layer (mol/m3) + real(r8), intent(in) :: conc_o2 ! O2 conc in layer (mol/m3) + + ! Arguments (out) + real(r8), intent(out) :: ch4_oxid ! Oxidation predicted by method A (temperature & enzyme limited) (mol CH4/m3/s) + real(r8), intent(out) :: o2_oxid ! Oxidation predicted by method A (temperature & enzyme limited) (mol O2/m3/s) + + ! + ! !LOCAL VARIABLES: + real(r8):: t0 ! Base temperature for Q10 + real(r8):: porevol ! air-filled volume ratio to total soil volume + real(r8):: h2osoi_vol_min ! h2osoi_vol restricted to be below watsat + real(r8):: conc_ch4_rel ! concentration with respect to water volume (mol/m^3 water) + real(r8):: conc_o2_rel ! concentration with respect to water volume (mol/m^3 water) + real(r8):: smp_fact ! factor for reduction based on soil moisture (unitless) + real(r8):: k_h_cc, k_h_inv ! see functions below for description + real(r8):: k_m_eff ! effective k_m + real(r8):: vmax_eff ! effective vmax + ! ch4 oxidation parameters + real(r8) :: vmax_ch4_oxid ! oxidation rate constant (= 45.e-6_r8 * 1000._r8 / 3600._r8) [mol/m3-w/s]; + real(r8) :: k_m ! Michaelis-Menten oxidation rate constant for CH4 concentration + real(r8) :: q10_ch4oxid ! Q10 oxidation constant + real(r8) :: smp_crit ! Critical soil moisture potential + real(r8) :: k_m_o2 ! Michaelis-Menten oxidation rate constant for O2 concentration + real(r8) :: k_m_unsat ! Michaelis-Menten oxidation rate constant for CH4 concentration + real(r8) :: vmax_oxid_unsat ! (= 45.e-6_r8 * 1000._r8 / 3600._r8 / 10._r8) [mol/m3-w/s] + ! + !----------------------------------------------------------------------- + + + ! Set oxidation parameters + vmax_ch4_oxid = params_inst%vmax_ch4_oxid + k_m = params_inst%k_m + q10_ch4oxid = params_inst%q10_ch4oxid + smp_crit = params_inst%smp_crit + k_m_o2 = params_inst%k_m_o2 + k_m_unsat = params_inst%k_m_unsat + vmax_oxid_unsat = params_inst%vmax_oxid_unsat + + t0 = tfrz + 12._r8 ! Walter, for Michigan site where the 45 M/h comes from + + ! only unsaturated layer considered here + k_m_eff = k_m_unsat + vmax_eff = vmax_oxid_unsat + + ! maximum water content in unsaturated layer is field capacity + h2osoi_vol_min = min(watsat*(sucsat/psi_fc)**(1._r8/bsw), h2osoi_vol) + + porevol = max(watsat - h2osoi_vol_min, 0._r8) + + ! moisture limitation + if (smp_l < 0._r8) then + ! Schnell & King, 1996, Figure 3 + smp_fact = exp(-min(smp_l,psi_fc)/smp_crit) + else + smp_fact = 1._r8 + end if + + k_h_inv = exp(-c_h_inv(1) * (1._r8 / t_soisno - 1._r8 / kh_tbase) + log (kh_theta(1))) + k_h_cc = t_soisno / k_h_inv * rgasLatm ! (4.21) Wania [(mol/m3w) / (mol/m3g)] + conc_ch4_rel = conc_ch4 / (h2osoi_vol_min + porevol/k_h_cc) + + k_h_inv = exp(-c_h_inv(2) * (1._r8 / t_soisno - 1._r8 / kh_tbase) + log (kh_theta(2))) + k_h_cc = t_soisno / k_h_inv * rgasLatm ! (4.21) Wania [(mol/m3w) / (mol/m3g)] + conc_o2_rel = conc_o2 / (h2osoi_vol_min + porevol/k_h_cc) + + ch4_oxid = vmax_eff * h2osoi_vol_min* conc_ch4_rel / (k_m_eff + conc_ch4_rel) & + ![mol/m3-t/s] [mol/m3-w/s] [m3-w/m3-t] [mol/m3-w] [mol/m3-w] [mol/m3-w] + * conc_o2_rel / (k_m_o2 + conc_o2_rel) & + * q10_ch4oxid ** ((t_soisno - t0) / 10._r8) * smp_fact + + ! For all landunits / levels, prevent oxidation if at or below freezing + if (t_soisno <= tfrz) ch4_oxid = 0._r8 + + o2_oxid = ch4_oxid * 2._r8 + + end subroutine surface_layer_oxid + !----------------------------------------------------------------------- subroutine ch4_tran (bounds, & num_methc, filter_methc, & @@ -3355,6 +3600,7 @@ subroutine ch4_tran (bounds, & ! are not tracked and oxidation is not allowed inside snow and lake water. ! Diffusivity is set based on soil texture and organic matter fraction. A Crank-Nicholson solution is used. ! Then CH4 diffusive flux is calculated and consistency is checked. + ! 3.27.2026: Surface layer added above resolved soil layers ! !USES: use clm_time_manager , only : get_step_size_real, get_nstep @@ -3383,27 +3629,27 @@ subroutine ch4_tran (bounds, & integer :: jtop(bounds%begc:bounds%endc) ! top level at each column integer :: iter ! iteration counter when dtime_ch4 < dtime real(r8) :: dtime ! land model time step (sec) - real(r8) :: at (bounds%begc:bounds%endc,0:nlevsoi) ! "a" vector for tridiagonal matrix - real(r8) :: bt (bounds%begc:bounds%endc,0:nlevsoi) ! "b" vector for tridiagonal matrix - real(r8) :: ct (bounds%begc:bounds%endc,0:nlevsoi) ! "c" vector for tridiagonal matrix - real(r8) :: rt (bounds%begc:bounds%endc,0:nlevsoi) ! "r" vector for tridiagonal solution + real(r8) :: at (bounds%begc:bounds%endc,-1:nlevsoi) ! "a" vector for tridiagonal matrix + real(r8) :: bt (bounds%begc:bounds%endc,-1:nlevsoi) ! "b" vector for tridiagonal matrix + real(r8) :: ct (bounds%begc:bounds%endc,-1:nlevsoi) ! "c" vector for tridiagonal matrix + real(r8) :: rt (bounds%begc:bounds%endc,-1:nlevsoi) ! "r" vector for tridiagonal solution real(r8) :: f_a ! air-filled fraction of available pore space - real(r8) :: diffus (bounds%begc:bounds%endc,0:nlevsoi) ! diffusivity (m2/s) + real(r8) :: diffus (bounds%begc:bounds%endc,-1:nlevsoi) ! diffusivity (m2/s) real(r8) :: k_h_inv ! 1/Henry's Law Constant in Latm/mol - real(r8) :: k_h_cc(bounds%begc:bounds%endc,0:nlevsoi,ngases) ! ratio of mol/m3 in liquid to mol/m3 in gas + real(r8) :: k_h_cc(bounds%begc:bounds%endc,-1:nlevsoi,ngases) ! ratio of mol/m3 in liquid to mol/m3 in gas real(r8) :: dzj ! - real(r8) :: dp1_zp1 (bounds%begc:bounds%endc,0:nlevsoi) ! diffusivity/delta_z for next j - real(r8) :: dm1_zm1 (bounds%begc:bounds%endc,0:nlevsoi) ! diffusivity/delta_z for previous j + real(r8) :: dp1_zp1 (bounds%begc:bounds%endc,-1:nlevsoi) ! diffusivity/delta_z for next j + real(r8) :: dm1_zm1 (bounds%begc:bounds%endc,-1:nlevsoi) ! diffusivity/delta_z for previous j real(r8) :: t_soisno_c ! soil temperature (C) (-nlevsno+1:nlevsoi) real(r8) :: eps ! either epsilon_a or epsilon_w, depending on where in soil, wrt WT real(r8) :: deficit ! mol CH4 /m^2 that must be subtracted from diffusive flux to atm. to make up ! for keeping concentrations always above zero - real(r8) :: conc_ch4_bef(bounds%begc:bounds%endc,1:nlevsoi) ! concentration at the beginning of the timestep + real(r8) :: conc_ch4_bef(bounds%begc:bounds%endc,0:nlevsoi) ! concentration at the beginning of the timestep real(r8) :: errch4(bounds%begc:bounds%endc) ! Error (Mol CH4 /m^2) [+ = too much CH4] - real(r8) :: conc_ch4_rel(bounds%begc:bounds%endc,0:nlevsoi) ! Concentration per volume of air or water - real(r8) :: conc_o2_rel(bounds%begc:bounds%endc,0:nlevsoi) ! Concentration per volume of air or water - real(r8) :: conc_ch4_rel_old(bounds%begc:bounds%endc,0:nlevsoi) ! Concentration during last Crank-Nich. loop - real(r8) :: h2osoi_vol_min(bounds%begc:bounds%endc,1:nlevsoi) ! h2osoi_vol restricted to be <= watsat + real(r8) :: conc_ch4_rel(bounds%begc:bounds%endc,-1:nlevsoi) ! Concentration per volume of air or water + real(r8) :: conc_o2_rel(bounds%begc:bounds%endc,-1:nlevsoi) ! Concentration per volume of air or water + real(r8) :: conc_ch4_rel_old(bounds%begc:bounds%endc,-1:nlevsoi) ! Concentration during last Crank-Nich. loop + real(r8) :: h2osoi_vol_min(bounds%begc:bounds%endc,0:nlevsoi) ! h2osoi_vol restricted to be <= watsat real(r8), parameter :: smallnumber = 1.e-12_r8 real(r8) :: snowdiff ! snow diffusivity (m^2/s) real(r8) :: snowres(bounds%begc:bounds%endc) ! Cumulative Snow resistance (s/m). Also includes @@ -3414,13 +3660,13 @@ subroutine ch4_tran (bounds, & real(r8) :: airfrac ! air fraction in snow real(r8) :: waterfrac ! water fraction in snow real(r8) :: icefrac ! ice fraction in snow - real(r8) :: epsilon_t (bounds%begc:bounds%endc,1:nlevsoi,1:ngases) ! - real(r8) :: epsilon_t_old (bounds%begc:bounds%endc,1:nlevsoi,1:ngases) ! epsilon_t from last time step !Currently deprecated - real(r8) :: source (bounds%begc:bounds%endc,1:nlevsoi,1:ngases) ! source - real(r8) :: source_old (bounds%begc:bounds%endc,1:nlevsoi,1:ngases) ! source from last time step !Currently deprecated + real(r8) :: epsilon_t (bounds%begc:bounds%endc,0:nlevsoi,1:ngases) ! + real(r8) :: epsilon_t_old (bounds%begc:bounds%endc,0:nlevsoi,1:ngases) ! epsilon_t from last time step !Currently deprecated + real(r8) :: source (bounds%begc:bounds%endc,0:nlevsoi,1:ngases) ! source + real(r8) :: source_old (bounds%begc:bounds%endc,0:nlevsoi,1:ngases) ! source from last time step !Currently deprecated real(r8) :: om_frac ! organic matter fraction real(r8) :: o2demand, ch4demand ! mol/m^3/s - real(r8) :: liqfrac(bounds%begc:bounds%endc, 1:nlevsoi) + real(r8) :: liqfrac(bounds%begc:bounds%endc, 0:nlevsoi) real(r8) :: capthick ! (mm) min thickness before assuming h2osfc is impermeable real(r8) :: satpow ! exponent on watsat for saturated soil solute diffusion real(r8) :: scale_factor_gasdiff ! For sensitivity tests; convection would allow this to be > 1 @@ -3428,6 +3674,12 @@ subroutine ch4_tran (bounds, & real(r8) :: organic_max ! organic matter content (kg/m3) where soil is assumed to act like peat real(r8) :: aereoxid ! fraction of methane flux entering aerenchyma rhizosphere + real(r8) :: surface_layer_thickness_old(bounds%begc:bounds%endc) ! + real(r8) :: surface_layer_ch4_oxid(bounds%begc:bounds%endc) ! + real(r8) :: surface_layer_o2_oxid(bounds%begc:bounds%endc) ! + real(r8) :: surface_layer_ch4stress + real(r8) :: surface_layer_o2stress + real(r8) :: dzm1 real(r8), pointer :: ch4_prod_depth (:,:) real(r8), pointer :: ch4_oxid_depth (:,:) real(r8), pointer :: ch4_aere_depth (:,:) @@ -3444,9 +3696,13 @@ subroutine ch4_tran (bounds, & real(r8), pointer :: co2_decomp_depth (:,:) real(r8), pointer :: conc_o2 (:,:) real(r8), pointer :: conc_ch4 (:,:) + real(r8), pointer :: surface_layer_conc_ch4(:) + real(r8), pointer :: surface_layer_conc_o2(:) + real(r8), pointer :: surface_layer_thickness(:) - integer :: nstep ! time step number - character(len=32) :: subname='ch4_tran' ! subroutine name + integer :: nstep ! time step number + character(len=128) :: err_msg='unset' ! error message + character(len=32) :: subname='ch4_tran' ! subroutine name !----------------------------------------------------------------------- SHR_ASSERT_ALL_FL((ubound(jwt) == (/bounds%endc/)), sourcefile, __LINE__) @@ -3459,6 +3715,8 @@ subroutine ch4_tran (bounds, & bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" (nlevgrnd) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) + smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] + sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) cellorg => soilstate_inst%cellorg_col , & ! Input: [real(r8) (:,:) ] column 3D org (kg/m^3 organic matter) (nlevgrnd) t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) (-nlevsno+1:nlevsoi) @@ -3494,6 +3752,9 @@ subroutine ch4_tran (bounds, & co2_decomp_depth => ch4_inst%co2_decomp_depth_unsat_col ! Output: [real(r8) (:,:) ] CO2 production during decomposition in each soil layer (nlevsoi) (mol/m3/s) conc_ch4 => ch4_inst%conc_ch4_unsat_col ! Output: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) conc_o2 => ch4_inst%conc_o2_unsat_col ! Output: [real(r8) (:,:) ] O2 conc in each soil layer (mol/m3) (nlevsoi) + surface_layer_conc_ch4 => ch4_inst%conc_ch4_surface_layer_unsat_col ! Output: [real(r8) (:) ] CH4 conc surface soil layer (mol/m3) + surface_layer_conc_o2 => ch4_inst%conc_o2_surface_layer_unsat_col ! Output: [real(r8) (:) ] O2 conc surface soil layer (mol/m3) + surface_layer_thickness => ch4_inst%surface_layer_thickness_unsat_col ! Output: [real(r8) (:) ] surface soil layer thickness (m) else ! saturated o2_decomp_depth => ch4_inst%o2_decomp_depth_sat_col ! Output: [real(r8) (:,:) ] O2 consumption during decomposition in each soil layer (nlevsoi) (mol/m3/s) o2stress => ch4_inst%o2stress_sat_col ! Output: [real(r8) (:,:) ] Ratio of oxygen available to that demanded by roots, aerobes, & methanotrophs (nlevsoi) @@ -3511,6 +3772,9 @@ subroutine ch4_tran (bounds, & co2_decomp_depth => ch4_inst%co2_decomp_depth_sat_col ! Output: [real(r8) (:,:) ] CO2 production during decomposition in each soil layer (nlevsoi) (mol/m3/s) conc_ch4 => ch4_inst%conc_ch4_sat_col ! Output: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) conc_o2 => ch4_inst%conc_o2_sat_col ! Output: [real(r8) (:,:) ] O2 conc in each soil layer (mol/m3) (nlevsoi) + surface_layer_conc_ch4 => ch4_inst%conc_ch4_surface_layer_sat_col ! Output: [real(r8) (:) ] CH4 conc surface water layer (mol/m3) (nlevsoi) + surface_layer_conc_o2 => ch4_inst%conc_o2_surface_layer_sat_col ! Output: [real(r8) (:) ] O2 conc surface water layer (mol/m3) (nlevsoi) + surface_layer_thickness => ch4_inst%surface_layer_thickness_sat_col ! Output: [real(r8) (:) ] surface water layer thickness (m) endif ! Get land model time step @@ -3527,6 +3791,28 @@ subroutine ch4_tran (bounds, & ! Set shared constant organic_max = CNParamsShareInst%organic_max + ! Set the Henry's Law coefficients + ! k_h_inv uses (4.12) Wania (L atm/mol) + ! k_h_cc uses (4.21) Wania [(mol/m3w) / (mol/m3g)] + do j = -1,nlevsoi + do fc = 1, num_methc + c = filter_methc (fc) + + do s=1,2 + if (j == -1) then + k_h_inv = exp(-c_h_inv(s) * (1._r8 / t_grnd(c) - 1._r8 / kh_tbase) + log (kh_theta(s))) + k_h_cc(c,j,s) = t_grnd(c) / k_h_inv * rgasLatm + elseif (j == 0) then ! use j=1 temperature for surface layer + k_h_inv = exp(-c_h_inv(s) * (1._r8 / t_soisno(c,1) - 1._r8 / kh_tbase) + log (kh_theta(s))) + k_h_cc(c,j,s) = t_soisno(c,1) / k_h_inv * rgasLatm + else + k_h_inv = exp(-c_h_inv(s) * (1._r8 / t_soisno(c,j) - 1._r8 / kh_tbase) + log (kh_theta(s))) + k_h_cc(c,j,s) = t_soisno(c,j) / k_h_inv * rgasLatm + end if + end do + end do + end do + ! Perform competition for oxygen and methane in each soil layer if demands over the course of the timestep ! exceed that available. Assign to each process in proportion to the quantity demanded in the absense of ! the limitation. @@ -3568,7 +3854,7 @@ subroutine ch4_tran (bounds, & ! Reset oxidation ch4_oxid_depth(c,j) = ch4_oxid_depth(c,j) * ch4stress(c,j) o2_oxid_depth(c,j) = o2_oxid_depth(c,j) * ch4stress(c,j) - else ! oxygen limited + else ! oxygen limited if (ch4stress(c,j) < 1._r8) then ! Recalculate methane limitation ch4demand = ch4_aere_depth(c,j) + ch4_ebul_depth(c,j) @@ -3603,32 +3889,152 @@ subroutine ch4_tran (bounds, & enddo enddo + ! Calculate surface layer oxidation including potential ebullition input from soil layers - ! Set the Henry's Law coefficients - do j = 0,nlevsoi - do fc = 1, num_methc - c = filter_methc (fc) + ! Also set values for surface water layer (no oxidation, but calculate thickness and concentrations + ! In the surface layer, unsat oxidation is a sink only, for sat no sinks/sources, just a diffusive resistance - do s=1,2 - if (j == 0) then - k_h_inv = exp(-c_h_inv(s) * (1._r8 / t_grnd(c) - 1._r8 / kh_tbase) + log (kh_theta(s))) - ! (4.12) Wania (L atm/mol) - k_h_cc(c,j,s) = t_grnd(c) / k_h_inv * rgasLatm ! (4.21) Wania [(mol/m3w) / (mol/m3g)] - else - k_h_inv = exp(-c_h_inv(s) * (1._r8 / t_soisno(c,j) - 1._r8 / kh_tbase) + log (kh_theta(s))) - ! (4.12) Wania (L atm/mol) - k_h_cc(c,j,s) = t_soisno(c,j) / k_h_inv * rgasLatm ! (4.21) Wania [(mol/m3w) / (mol/m3g)] - end if - end do - end do - end do + do fc = 1, num_methc + c = filter_methc(fc) + g = col%gridcell(c) + + ! Initialize variables + surface_layer_ch4_oxid(c) = 0._r8 + surface_layer_o2_oxid(c) = 0._r8 + surface_layer_thickness_old(c) = surface_layer_thickness(c) + + ! calculate ch4, o2 concentration in virtual layer + if (sat==0) then + if (jwt(c)==0) then + ! if water table is at the soil surface, include ebullition flux + call surface_layer_oxid(watsat(c,1), sucsat(c,1), bsw(c,1), & + h2osoi_vol(c,1), smp_l(c,1), t_soisno(c,1), & + surface_layer_conc_ch4(c)+ch4_ebul_total(c)/surface_layer_thickness(c)*dtime, surface_layer_conc_o2(c), & + surface_layer_ch4_oxid(c), surface_layer_o2_oxid(c)) + else + call surface_layer_oxid(watsat(c,1), sucsat(c,1), bsw(c,1), & + h2osoi_vol(c,1), smp_l(c,1), t_soisno(c,1), & + surface_layer_conc_ch4(c), surface_layer_conc_o2(c), & + surface_layer_ch4_oxid(c), surface_layer_o2_oxid(c)) + endif + + ! assign surface (soil) layer thickness (min 5 cm) + !surface_layer_thickness(c) = 0.05_r8 ! static value + ! crude estimate of mean unsaturated thickness + surface_layer_thickness(c) = max((3._r8*col%micro_sigma(c)-1e-3_r8*h2osfc(c)),0.05_r8) + endif + + ! assume a surface water (not soil) layer for sat=1 + ! use instead of using resistance later (removed h2osfc resistance) + if (sat==1) then + surface_layer_ch4_oxid(c) = 0._r8 + surface_layer_o2_oxid(c) = 0._r8 + ! assign surface (water) layer thickness (min 1 cm) + surface_layer_thickness(c) = max(1e-3_r8*h2osfc(c),0.01_r8) + endif + ! adjust concentration for changes in surface layer thickness + surface_layer_conc_ch4(c) = surface_layer_conc_ch4(c)*(surface_layer_thickness_old(c)/surface_layer_thickness(c)) + surface_layer_conc_o2(c) = surface_layer_conc_o2(c)*(surface_layer_thickness_old(c)/surface_layer_thickness(c)) - ! Set the source term for each species (no need to do j=0, since epsilon_t and source not used there) + ! calculate stress for unsaturated layer + o2demand = surface_layer_o2_oxid(c) + if (o2demand > 0._r8) then + if ( (surface_layer_conc_o2(c) / dtime) > o2demand )then + surface_layer_o2stress = 1._r8 + else + surface_layer_o2stress = (surface_layer_conc_o2(c) / dtime) / o2demand + end if + else + surface_layer_o2stress= 1._r8 + end if + + ch4demand = surface_layer_ch4_oxid(c) + if (ch4demand > 0._r8) then + ! include ebullition flux if water table is at the soil surface + if (jwt(c)==0) then + surface_layer_ch4stress = min((surface_layer_conc_ch4(c)/dtime + ch4_ebul_total(c)/surface_layer_thickness(c)) / ch4demand, 1._r8) + else + surface_layer_ch4stress = min((surface_layer_conc_ch4(c) / dtime ) / ch4demand, 1._r8) + endif + else + surface_layer_ch4stress = 1._r8 + end if + + ! Resolve methane oxidation + if (surface_layer_o2stress < 1._r8 .or. surface_layer_ch4stress < 1._r8) then + if (surface_layer_ch4stress <= surface_layer_o2stress) then ! methane limited + ! Reset oxidation + surface_layer_ch4_oxid(c) = surface_layer_ch4_oxid(c) * surface_layer_ch4stress + surface_layer_o2_oxid(c) = surface_layer_o2_oxid(c) * surface_layer_ch4stress + else ! oxygen limited + ! Reset oxidation + surface_layer_ch4_oxid(c) = surface_layer_ch4_oxid(c) * surface_layer_o2stress + surface_layer_o2_oxid(c) = surface_layer_o2_oxid(c) * surface_layer_o2stress + end if + end if + + ! Set source values in surface layer and perform error checks + j = 0 + source(c,j,1) = -surface_layer_ch4_oxid(c) ! CH4 [mol/m3-total/s] + source(c,j,2) = -surface_layer_o2_oxid(c) ! O2 [mol/m3/s] + + ! methane check + if (jwt(c)==0) then + ! water table at soil surface (add ebullition flux) + if ((source(c,j,1) + surface_layer_conc_ch4(c)/dtime + ch4_ebul_total(c)/surface_layer_thickness(c)) < -1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') & + 'Methane demands exceed methane available. Error in methane competition (mol/m^3/s), c,j:', & + source(c,j,1) + surface_layer_conc_ch4(c) / dtime & + + ch4_ebul_total(c)/surface_layer_thickness(c), c, j + else if (surface_layer_ch4stress < 1._r8 .and. & + (source(c,j,1) + surface_layer_conc_ch4(c) / dtime & + + ch4_ebul_total(c)/surface_layer_thickness(c)) > 1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') 'Methane limited, yet some left over. Error in methane competition (mol/m^3/s), c,j:', & + source(c,j,1) + surface_layer_conc_ch4(c) / dtime + ch4_ebul_total(c)/surface_layer_thickness(c), c, j + endif + else + ! water table below soil surface + if (source(c,j,1) + surface_layer_conc_ch4(c) / dtime < -1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') 'Methane demands exceed methane available. Error in methane competition (mol/m^3/s), c,j:', & + source(c,j,1) + surface_layer_conc_ch4(c) / dtime, c, j + else if (surface_layer_ch4stress < 1._r8 .and. source(c,j,1) + surface_layer_conc_ch4(c) / dtime > 1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') 'Methane limited, yet some left over. Error in methane competition (mol/m^3/s), c,j:', & + source(c,j,1) + surface_layer_conc_ch4(c) / dtime, c, j + end if + endif + + ! Endrun if any CH4 errors identified + if (trim(err_msg) /= 'unset') then + call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & + msg=err_msg//errMsg(sourcefile, __LINE__)) + endif + + ! oxygen check + if (source(c,j,2) + surface_layer_conc_o2(c) / dtime < -1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') 'Oxygen demands exceed oxygen available. Error in oxygen competition (mol/m^3/s), c,j:', & + source(c,j,2) + surface_layer_conc_o2(c) / dtime, c, j + else if (surface_layer_o2stress < 1._r8 .and. source(c,j,2) + surface_layer_conc_o2(c) / dtime > 1.e-12_r8) then + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') 'Oxygen limited, yet some left over. Error in oxygen competition (mol/m^3/s), c,j:', & + source(c,j,2) + surface_layer_conc_o2(c) / dtime, c, j + end if + + ! Endrun if any O2 errors identified + if (trim(err_msg) /= 'unset') then + call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & + msg=err_msg//errMsg(sourcefile, __LINE__)) + endif + + conc_ch4_bef(c,j) = surface_layer_conc_ch4(c) !For Balance Check + + enddo + + ! Set the source term for each species (no need to do j=-1, since epsilon_t and source not used there) ! Note that because of the semi-implicit diffusion and the 30 min timestep combined with explicit ! sources, occasionally negative concentration will result. In this case it is brought to zero and the ! surface flux is adjusted to conserve. This results in some inaccuracy as compared to a shorter timestep ! or iterative solution. + do j = 1,nlevsoi do fc = 1, num_methc c = filter_methc (fc) @@ -3641,44 +4047,43 @@ subroutine ch4_tran (bounds, & source(c,j,1) = ch4_prod_depth(c,j) - ch4_oxid_depth(c,j) - & ch4_aere_depth(c,j) - ch4_ebul_depth(c,j) ! [mol/m3-total/s] + ! aerenchyma added to surface flux below ! ebul added to soil depth just above WT + if (source(c,j,1) + conc_ch4(c,j) / dtime < -1.e-12_r8)then - write(iulog,*) 'Methane demands exceed methane available. Error in methane competition (mol/m^3/s), c,j:', & + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') & + 'Methane demands exceed methane available. Error in methane competition (mol/m^3/s), c,j:', & source(c,j,1) + conc_ch4(c,j) / dtime, c, j - g = col%gridcell(c) - write(iulog,*)'Latdeg,Londeg=',grc%latdeg(g),grc%londeg(g) - call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & - msg=' ERROR: Methane demands exceed methane available.'& - //errMsg(sourcefile, __LINE__)) else if (ch4stress(c,j) < 1._r8 .and. source(c,j,1) + conc_ch4(c,j) / dtime > 1.e-12_r8) then - write(iulog,*) 'Methane limited, yet some left over. Error in methane competition (mol/m^3/s), c,j:', & + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') & + 'Methane limited, yet some left over. Error in methane competition (mol/m^3/s), c,j:', & source(c,j,1) + conc_ch4(c,j) / dtime, c, j - g = col%gridcell(c) - write(iulog,*)'Latdeg,Londeg=',grc%latdeg(g),grc%londeg(g) - call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & - msg=' ERROR: Methane limited, yet some left over.'//& - errMsg(sourcefile, __LINE__)) end if + ! Endrun if any CH4 errors identified + if (trim(err_msg) /= 'unset') then + call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & + msg=err_msg//errMsg(sourcefile, __LINE__)) + endif + source(c,j,2) = -o2_oxid_depth(c,j) - o2_decomp_depth(c,j) + o2_aere_depth(c,j) ! O2 [mol/m3/s] if (source(c,j,2) + conc_o2(c,j) / dtime < -1.e-12_r8) then - write(iulog,*) 'Oxygen demands exceed oxygen available. Error in oxygen competition (mol/m^3/s), c,j:', & + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') & + 'Oxygen demands exceed oxygen available. Error in oxygen competition (mol/m^3/s), c,j:', & source(c,j,2) + conc_o2(c,j) / dtime, c, j - g = col%gridcell(c) - write(iulog,*)'Latdeg,Londeg=',grc%latdeg(g),grc%londeg(g) - call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & - msg=' ERROR: Oxygen demands exceed oxygen available.'//& - errMsg(sourcefile, __LINE__) ) else if (o2stress(c,j) < 1._r8 .and. source(c,j,2) + conc_o2(c,j) / dtime > 1.e-12_r8) then - write(iulog,*) 'Oxygen limited, yet some left over. Error in oxygen competition (mol/m^3/s), c,j:', & + write(unit=err_msg,fmt='(a,e14.8,i8,i8)') & + 'Oxygen limited, yet some left over. Error in oxygen competition (mol/m^3/s), c,j:', & source(c,j,2) + conc_o2(c,j) / dtime, c, j - g = col%gridcell(c) - write(iulog,*)'Latdeg,Londeg=',grc%latdeg(g),grc%londeg(g) - call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & - msg=' ERROR: Oxygen limited, yet some left over.'//errMsg(sourcefile, __LINE__)) end if + ! Endrun if any O2 errors identified + if (trim(err_msg) /= 'unset') then + call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & + msg=err_msg//errMsg(sourcefile, __LINE__)) + endif + conc_ch4_bef(c,j) = conc_ch4(c,j) !For Balance Check enddo ! fc enddo ! j @@ -3697,19 +4102,45 @@ subroutine ch4_tran (bounds, & c = filter_methc(fc) if (jwt(c) /= 0) then source(c,jwt(c),1) = source(c,jwt(c),1) + ch4_ebul_total(c)/dz(c,jwt(c)) + elseif (jwt(c) == 0 .and. sat==0) then + source(c,jwt(c),1) = source(c,jwt(c),1) + ch4_ebul_total(c)/surface_layer_thickness(c) endif enddo ! fc ! Calculate concentration relative to m^3 of air or water: needed for the diffusion - do j = 0,nlevsoi + do j = -1,nlevsoi do fc = 1, num_methc c = filter_methc (fc) g = col%gridcell(c) - if (j == 0) then + if (j == -1) then conc_ch4_rel(c,j) = c_atm(g,1) conc_o2_rel(c,j) = c_atm(g,2) - else + + else if (j == 0) then + if (ch4frzout) then + liqfrac(c,j) = max(0.05_r8, (h2osoi_liq(c,1)/denh2o+smallnumber)/ & + (h2osoi_liq(c,1)/denh2o+h2osoi_ice(c,1)/denice+smallnumber)) + else + liqfrac(c,j) = 1._r8 + end if + + if (sat == 0) then + ! maximum water content in unsaturated layer is field capacity; use j=1 h2osoi liq/ice values to describe surface layer conditions + h2osoi_vol_min(c,j) = min(watsat(c,1)*(sucsat(c,1)/psi_fc)**(1._r8/bsw(c,1)), h2osoi_vol(c,1)) + + do s=1,2 + epsilon_t(c,j,s) = watsat(c,1)-(1._r8-k_h_cc(c,j,s))*h2osoi_vol_min(c,j)*liqfrac(c,j) + end do + else ! sat=1; surface water layer (no soil media) + do s=1,2 + epsilon_t(c,j,s) = liqfrac(c,j) + end do + endif + conc_ch4_rel(c,j) = surface_layer_conc_ch4(c)/epsilon_t(c,j,1) + conc_o2_rel(c,j) = surface_layer_conc_o2(c)/epsilon_t(c,j,2) + + else ! resolved soil layers (1:nlevsoi) h2osoi_vol_min(c,j) = min(watsat(c,j), h2osoi_vol(c,j)) if (ch4frzout) then liqfrac(c,j) = max(0.05_r8, (h2osoi_liq(c,j)/denh2o+smallnumber)/ & @@ -3719,7 +4150,7 @@ subroutine ch4_tran (bounds, & end if if (j <= jwt(c)) then ! Above the WT do s=1,2 - epsilon_t(c,j,s) = watsat(c,j)- (1._r8-k_h_cc(c,j,s))*h2osoi_vol_min(c,j)*liqfrac(c,j) + epsilon_t(c,j,s) = watsat(c,j)-(1._r8-k_h_cc(c,j,s))*h2osoi_vol_min(c,j)*liqfrac(c,j) end do ! Partition between the liquid and gas phases. The gas phase will drive the diffusion. else ! Below the WT @@ -3733,11 +4164,9 @@ subroutine ch4_tran (bounds, & end do end do - ! Loop over species do s = 1, 2 ! 1=CH4; 2=O2; 3=CO2 - ! Adjust the grnd_ch4_cond to keep it positive, and add the snow resistance & pond resistance do j = -nlevsno + 1,0 do fc = 1, num_methc @@ -3792,123 +4221,195 @@ subroutine ch4_tran (bounds, & pondres = pondz / ponddiff end if - ! Now add new h2osfc form - if (.not. lake .and. sat == 1 .and. frac_h2osfc(c) > 0._r8) then - if (t_h2osfc(c) >= tfrz) then - t_soisno_c = t_h2osfc(c) - tfrz - ponddiff = (d_con_w(s,1) + d_con_w(s,2)*t_soisno_c + d_con_w(s,3)*t_soisno_c**2) * 1.e-9_r8 & - * scale_factor_liqdiff - pondz = h2osfc(c) / 1000._r8 / frac_h2osfc(c) ! Assume all h2osfc corresponds to sat area - ! mm / mm/m - pondres = pondres + pondz / ponddiff - else if (h2osfc(c)/frac_h2osfc(c) > capthick) then - ! assume surface ice is impermeable - pondres = 1/smallnumber - end if - end if - spec_grnd_cond(c,s) = 1._r8/(1._r8/grnd_ch4_cond(c) + snowres(c) + pondres) - end if - + end if end do ! fc end do ! j ! Determine gas diffusion and fraction of open pore (f_a) - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) g = col%gridcell(c) - t_soisno_c = t_soisno(c,j) - tfrz + if (j == 0) then ! surface layer, use j=1 values for soil properties + t_soisno_c = t_soisno(c,1) - tfrz + if (sat == 0) then + f_a = 1._r8 - h2osoi_vol_min(c,j) / watsat(c,1) + ! Provisionally calculate diffusivity as linear combination of the Millington-Quirk + ! expression in Wania (for peat) & Moldrup (for mineral soil) + eps = watsat(c,1)-h2osoi_vol_min(c,j) ! Air-filled fraction of total soil volume + if (organic_max > 0._r8) then + om_frac = min(params_inst%om_frac_sf*cellorg(c,1)/organic_max, 1._r8) + ! Use first power, not square as in iniTimeConst + else + om_frac = 1._r8 + end if + diffus(c,j) = (d_con_g(s,1) + d_con_g(s,2)*t_soisno_c) * 1.e-4_r8 * & + (om_frac * f_a**(10._r8/3._r8) / watsat(c,1)**2._r8 + & + (1._r8-om_frac) * (eps*eps) * f_a**(3._r8 / bsw(c,1)) ) & + * scale_factor_gasdiff - if (j <= jwt(c)) then ! Above the WT - f_a = 1._r8 - h2osoi_vol_min(c,j) / watsat(c,j) - ! Provisionally calculate diffusivity as linear combination of the Millington-Quirk - ! expression in Wania (for peat) & Moldrup (for mineral soil) - eps = watsat(c,j)-h2osoi_vol_min(c,j) ! Air-filled fraction of total soil volume - if (organic_max > 0._r8) then - om_frac = min(params_inst%om_frac_sf*cellorg(c,j)/organic_max, 1._r8) - ! Use first power, not square as in iniTimeConst - else - om_frac = 1._r8 - end if - diffus (c,j) = (d_con_g(s,1) + d_con_g(s,2)*t_soisno_c) * 1.e-4_r8 * & - (om_frac * f_a**(10._r8/3._r8) / watsat(c,j)**2._r8 + & - (1._r8-om_frac) * eps**2._r8 * f_a**(3._r8 / bsw(c,j)) ) & - * scale_factor_gasdiff - else ! Below the WT use saturated diffusivity and only water in epsilon_t - ! Note the following is not currently corrected for the effect on diffusivity of excess ice in soil under - ! lakes (which is currently experimental only). - eps = watsat(c,j) ! Water-filled fraction of total soil volume - diffus (c,j) = eps**satpow * (d_con_w(s,1) + d_con_w(s,2)*t_soisno_c + d_con_w(s,3)*t_soisno_c**2) * 1.e-9_r8 & - * scale_factor_liqdiff - if (t_soisno(c,j)<=tfrz) then - diffus(c,j) = diffus(c,j)*(h2osoi_liq(c,j)/denh2o+smallnumber)/ & - (h2osoi_liq(c,j)/denh2o+h2osoi_ice(c,j)/denice+smallnumber) - end if - endif ! Above/below the WT + endif + if (sat == 1) then + eps = watsat(c,1) ! Water-filled fraction of total soil volume + diffus (c,j) = eps**satpow * (d_con_w(s,1) + d_con_w(s,2)*t_soisno_c + d_con_w(s,3)*t_soisno_c**2) * 1.e-9_r8 & + * scale_factor_liqdiff + if (t_soisno(c,1)<=tfrz) then + diffus(c,j) = diffus(c,j)*liqfrac(c,j) + end if + + endif + else ! resolved soil layers + t_soisno_c = t_soisno(c,j) - tfrz + + if (j <= jwt(c)) then ! Above the WT + f_a = 1._r8 - h2osoi_vol_min(c,j) / watsat(c,j) + ! Provisionally calculate diffusivity as linear combination of the Millington-Quirk + ! expression in Wania (for peat) & Moldrup (for mineral soil) + eps = watsat(c,j)-h2osoi_vol_min(c,j) ! Air-filled fraction of total soil volume + if (organic_max > 0._r8) then + om_frac = min(params_inst%om_frac_sf*cellorg(c,j)/organic_max, 1._r8) + ! Use first power, not square as in iniTimeConst + else + om_frac = 1._r8 + end if + diffus (c,j) = (d_con_g(s,1) + d_con_g(s,2)*t_soisno_c) * 1.e-4_r8 * & + (om_frac * f_a**(10._r8/3._r8) / watsat(c,j)**2._r8 + & + (1._r8-om_frac) * (eps*eps) * f_a**(3._r8 / bsw(c,j)) ) & + * scale_factor_gasdiff + + else ! Below the WT use saturated diffusivity and only water in epsilon_t + ! Note the following is not currently corrected for the effect on diffusivity of excess ice in soil under + ! lakes (which is currently experimental only). + eps = watsat(c,j) ! Water-filled fraction of total soil volume + diffus (c,j) = eps**satpow * (d_con_w(s,1) + d_con_w(s,2)*t_soisno_c + d_con_w(s,3)*t_soisno_c**2) * 1.e-9_r8 & + * scale_factor_liqdiff + if (t_soisno(c,j)<=tfrz) then + diffus(c,j) = diffus(c,j)*(h2osoi_liq(c,j)/denh2o+smallnumber)/ & + (h2osoi_liq(c,j)/denh2o+h2osoi_ice(c,j)/denice+smallnumber) + end if + endif ! Above/below the WT + endif ! j=0 diffus(c,j) = max(diffus(c,j), smallnumber) ! Prevent overflow enddo ! fp enddo ! j - do j = 1,nlevsoi + ! Set up coefficients for tridiagonal solver + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) - ! Set up coefficients for tridiagonal solver. - if (j == 1 .and. j /= jwt(c) .and. j /= jwt(c)+1) then - dm1_zm1(c,j) = 1._r8/(1._r8/spec_grnd_cond(c,s)+dz(c,j)/(diffus(c,j)*2._r8)) - ! replace Diffusivity / Delta_z by conductance (grnd_ch4_cond) for top layer - dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - else if (j == 1 .and. j == jwt(c)) then - dm1_zm1(c,j) = 1._r8/(1._r8/spec_grnd_cond(c,s)+dz(c,j)/(diffus(c,j)*2._r8)) - ! layer resistance mult. by k_h_cc for dp1_zp1 term - dp1_zp1(c,j) = 2._r8/(dz(c,j)*k_h_cc(c,j,s)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - else if (j == 1) then ! water table at surface: multiply ground resistance by k_h_cc - dm1_zm1(c,j) = 1._r8/(k_h_cc(c,j-1,s)/spec_grnd_cond(c,s)+dz(c,j)/(diffus(c,j)*2._r8)) - ! air concentration will be mult. by k_h_cc below - dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - else if (j <= nlevsoi-1 .and. j /= jwt(c) .and. j /= jwt(c)+1) then - dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j-1)/diffus(c,j-1)) - dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - else if (j <= nlevsoi-1 .and. j == jwt(c)) then ! layer resistance mult. by k_h_cc for dp1_zp1 term - dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j-1)/diffus(c,j-1)) - dp1_zp1(c,j) = 2._r8/(dz(c,j)*k_h_cc(c,j,s)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - ! Concentration in layer will be mult. by k_h_cc below - else if (j <= nlevsoi-1) then ! j==jwt+1: layer above resistance mult. by k_h_cc for dm1_zm1 term - dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j-1)*k_h_cc(c,j-1,s)/diffus(c,j-1)) - ! Concentration in layer above will be mult. by k_h_cc below - dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) - else if (j /= jwt(c)+1) then ! j ==nlevsoi - dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j-1)/diffus(c,j-1)) - else ! jwt == nlevsoi-1: layer above resistance mult. by k_h_cc for dm1_zm1 term - dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j-1)*k_h_cc(c,j-1,s)/diffus(c,j-1)) - end if + ! set dz(j-1) + if (j==1) then + dzm1 = surface_layer_thickness(c) + elseif (j > 1) then + dzm1 = dz(c,j-1) + endif + + if (sat==0) then ! non-inundated conditions + + if (j == 0) then + ! for sat=0, there is an unsaturated surface soil layer, + ! even if all soil layers are saturated + + ! replace Diffusivity by conductance (grnd_ch4_cond) in dm1_zm1 + dm1_zm1(c,j) = 1._r8/(1._r8/spec_grnd_cond(c,s)+surface_layer_thickness(c)/(diffus(c,j)*2._r8)) + if (jwt(c) > 0) then + dp1_zp1(c,j) = 2._r8/(surface_layer_thickness(c)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + else ! jwt=0 (top soil layer saturated) + ! layer resistance mult. by k_h_cc for dp1_zp1 term + dp1_zp1(c,j) = 2._r8/(surface_layer_thickness(c)*k_h_cc(c,j,s)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + endif + else if (j <= nlevsoi-1 .and. j /= jwt(c) .and. j /= jwt(c)+1) then + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1/diffus(c,j-1)) + dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + else if (j <= nlevsoi-1 .and. j == jwt(c)) then ! layer resistance mult. by k_h_cc for dp1_zp1 term + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1/diffus(c,j-1)) + dp1_zp1(c,j) = 2._r8/(dz(c,j)*k_h_cc(c,j,s)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + ! Concentration in layer will be mult. by k_h_cc below + else if (j <= nlevsoi-1) then ! j==jwt+1: layer above resistance mult. by k_h_cc for dm1_zm1 term + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1*k_h_cc(c,j-1,s)/diffus(c,j-1)) + ! Concentration in layer above will be mult. by k_h_cc below + dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + else if (j /= jwt(c)+1) then ! j ==nlevsoi + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1/diffus(c,j-1)) + else ! jwt == nlevsoi-1: layer above resistance mult. by k_h_cc for dm1_zm1 term + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1*k_h_cc(c,j-1,s)/diffus(c,j-1)) + end if + + else ! sat=1, jwt=0 by definition + + if (j == 0) then + ! soil layers saturated and overlain by a layer of water + + ! replace Diffusivity by conductance (grnd_ch4_cond) in dm1_zm1, and multiply ground resistance by k_h_cc (water table at surface) + dm1_zm1(c,j) = 1._r8/(k_h_cc(c,j-1,s)/spec_grnd_cond(c,s)+surface_layer_thickness(c)/(diffus(c,j)*2._r8)) + dp1_zp1(c,j) = 2._r8/(surface_layer_thickness(c)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + else if (j <= nlevsoi-1) then + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1/diffus(c,j-1)) + dp1_zp1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dz(c,j+1)/diffus(c,j+1)) + else if (j == nlevsoi) then ! j=nlevsoi + dm1_zm1(c,j) = 2._r8/(dz(c,j)/diffus(c,j)+dzm1/diffus(c,j-1)) + end if + + endif ! sat; enddo ! fp; patch end do ! j; nlevsoi - + ! Perform a second loop for the tridiagonal coefficients since need dp1_zp1 and dm1_z1 at each depth - do j = 0,nlevsoi + do j = -1,nlevsoi do fc = 1, num_methc c = filter_methc (fc) g = col%gridcell(c) conc_ch4_rel_old(c,j) = conc_ch4_rel(c,j) - if (j > 0) dzj = dz(c,j) - if (j == 0) then ! top layer (atmosphere) doesn't change regardless of where WT is + ! assign layer thickness (ignore j = -1) + if (j == 0) dzj = surface_layer_thickness(c) + if (j > 0) dzj = dz(c,j) + + if (j == -1) then ! top layer (atmosphere) doesn't change regardless of where WT is at(c,j) = 0._r8 bt(c,j) = 1._r8 ct(c,j) = 0._r8 - rt(c,j) = c_atm(g,s) ! 0th level stays at constant atmospheric conc + rt(c,j) = c_atm(g,s) ! top level stays at constant atmospheric conc + + elseif (j == 0) then + if (sat == 0) then ! non-inundated conditions + ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) + if (jwt(c) == j) then + bt(c,j) = epsilon_t(c,j,s) / dtime_ch4 + 0.5_r8 / dzj * (dp1_zp1(c,j)*k_h_cc(c,j,s) + dm1_zm1(c,j)) + else + bt(c,j) = epsilon_t(c,j,s) / dtime_ch4 + 0.5_r8 / dzj * (dp1_zp1(c,j) + dm1_zm1(c,j)) + endif + ct(c,j) = -0.5_r8 / dzj * dp1_zp1(c,j) + endif + if (sat == 1) then ! inundated conditions + ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term + !at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) * k_h_cc(c,j-1,s) + at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) * k_h_cc(c,j,s) + bt(c,j) = epsilon_t(c,j,s) / dtime_ch4 + 0.5_r8 / dzj * (dp1_zp1(c,j) + dm1_zm1(c,j)) + ct(c,j) = -0.5_r8 / dzj * dp1_zp1(c,j) + endif + elseif (j < nlevsoi .and. j == jwt(c)) then ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) + ! no if-block needed; jwt always 0 for sat=1 bt(c,j) = epsilon_t(c,j,s) / dtime_ch4 + 0.5_r8 / dzj * (dp1_zp1(c,j)*k_h_cc(c,j,s) + dm1_zm1(c,j)) ct(c,j) = -0.5_r8 / dzj * dp1_zp1(c,j) elseif (j < nlevsoi .and. j == jwt(c)+1) then ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term - at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) * k_h_cc(c,j-1,s) + ! at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) * k_h_cc(c,j-1,s) + + if (sat==0) then + at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) * k_h_cc(c,j-1,s) + else ! j=1, jwt=0 -> j=0 also saturated so no k_h applied + at(c,j) = -0.5_r8 / dzj * dm1_zm1(c,j) + endif + bt(c,j) = epsilon_t(c,j,s) / dtime_ch4 + 0.5_r8 / dzj * (dp1_zp1(c,j) + dm1_zm1(c,j)) ct(c,j) = -0.5_r8 / dzj * dp1_zp1(c,j) elseif (j < nlevsoi) then @@ -3930,13 +4431,13 @@ subroutine ch4_tran (bounds, & do fc = 1, num_methc c = filter_methc (fc) - jtop(c) = 0 + jtop(c) = -1 end do if (s == 1) then ! CH4 - + ! Set rt, since it depends on conc - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) @@ -3946,18 +4447,56 @@ subroutine ch4_tran (bounds, & ! source_old could be removed later epsilon_t_old(c,j,s) = epsilon_t(c,j,s) ! epsilon_t acts like source also - dzj = dz(c,j) - if (j < nlevsoi .and. j == jwt(c)) then ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + + ! assign layer thickness + if (j == 0) then + dzj = surface_layer_thickness(c) + else + dzj = dz(c,j) + endif + + if (j == 0) then + if (sat == 0) then ! non-inundated conditions + if (jwt(c) > 0) then + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & + dm1_zm1(c,j) * (conc_ch4_rel(c,j) -conc_ch4_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + else ! jwt=0 (top soil layer saturated) + ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)*k_h_cc(c,j,s)) - & + dm1_zm1(c,j) * (conc_ch4_rel(c,j)-conc_ch4_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif + endif + if (sat == 1) then ! inundated conditions + ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & + dm1_zm1(c,j) * (conc_ch4_rel(c,j)-conc_ch4_rel(c,j-1)*k_h_cc(c,j-1,s))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif + + elseif (j < nlevsoi .and. j == jwt(c)) then ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)*k_h_cc(c,j,s)) - & dm1_zm1(c,j) * (conc_ch4_rel(c,j) -conc_ch4_rel(c,j-1))) + & 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) elseif (j < nlevsoi .and. j == jwt(c)+1) then ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term - rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & - 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & - dm1_zm1(c,j) * (conc_ch4_rel(c,j) -conc_ch4_rel(c,j-1)*k_h_cc(c,j-1,s))) + & - 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + if (sat==0) then + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & + dm1_zm1(c,j) * (conc_ch4_rel(c,j) -conc_ch4_rel(c,j-1)*k_h_cc(c,j-1,s))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + else !j=1/jwt=0 -> no k_h applied + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & + dm1_zm1(c,j) * (conc_ch4_rel(c,j) -conc_ch4_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif + elseif (j < nlevsoi) then rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_ch4_rel(c,j) + & 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_ch4_rel(c,j+1)-conc_ch4_rel(c,j)) - & @@ -3979,41 +4518,58 @@ subroutine ch4_tran (bounds, & enddo ! fc; column enddo ! j; nlevsoi - call Tridiagonal(bounds, 0, nlevsoi, & + call Tridiagonal(bounds, -1, nlevsoi, & jtop(bounds%begc:bounds%endc), & num_methc, filter_methc, & at(bounds%begc:bounds%endc, :), & bt(bounds%begc:bounds%endc, :), & ct(bounds%begc:bounds%endc, :), & rt(bounds%begc:bounds%endc, :), & - conc_ch4_rel(bounds%begc:bounds%endc, 0:nlevsoi)) + conc_ch4_rel(bounds%begc:bounds%endc, :)) ! Calculate net ch4 flux to the atmosphere from the surface (+ to atm) do fc = 1, num_methc c = filter_methc (fc) g = col%gridcell(c) if (jwt(c) /= 0) then ! WT not at the surface - ch4_surf_diff(c) = dm1_zm1(c,1) * ( (conc_ch4_rel(c,1)+conc_ch4_rel_old(c,1))/2._r8 & + ch4_surf_diff(c) = dm1_zm1(c,0) & + * ( (conc_ch4_rel(c,0)+conc_ch4_rel_old(c,0))/2._r8 & - c_atm(g,s)) ! [mol/m2/s] ch4_surf_ebul(c) = 0._r8 ! all the ebullition has already come out in the soil column (added to source) + ! Try adding directly to atm. to prevent destabilization of diffusion !ch4_surf_ebul(c) = ch4_ebul_total(c) ! [mol/m2/s] else ! WT at the surface; i.e., jwt(c)==0 - ch4_surf_diff(c) = dm1_zm1(c,1) * ( (conc_ch4_rel(c,1)+conc_ch4_rel_old(c,1))/2._r8 & - - c_atm(g,s)*k_h_cc(c,0,s)) ! [mol/m2/s] - ! atmospheric concentration gets mult. by k_h_cc as above - ch4_surf_ebul(c) = ch4_ebul_total(c) ! [mol/m2/s] + if (sat==0) then + ch4_surf_diff(c) = dm1_zm1(c,0) & + * ( (conc_ch4_rel(c,0)+conc_ch4_rel_old(c,0))/2._r8 & + - c_atm(g,s)) ! [mol/m2/s] + ch4_surf_ebul(c) = 0._r8 ! all the ebullition has already come out in the soil column (added to source) + + endif + if (sat==1) then + ch4_surf_diff(c) = dm1_zm1(c,0) & + * ( (conc_ch4_rel(c,0)+conc_ch4_rel_old(c,0))/2._r8 & + - c_atm(g,s)*k_h_cc(c,-1,s)) ! [mol/m2/s] + + ! atmospheric concentration gets mult. by k_h_cc as above + ch4_surf_ebul(c) = ch4_ebul_total(c) ! [mol/m2/s] + endif endif enddo ! Ensure that concentrations stay above 0 ! This should be done after the flux, so that the flux calculation is consistent. - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) if (conc_ch4_rel(c,j) < 0._r8) then - deficit = - conc_ch4_rel(c,j)*epsilon_t(c,j,1)*dz(c,j) ! Mol/m^2 added + if (j == 0) then + deficit = - conc_ch4_rel(c,j)*epsilon_t(c,j,1)*surface_layer_thickness(c) ! Mol/m^2 added + else + deficit = - conc_ch4_rel(c,j)*epsilon_t(c,j,1)*dz(c,j) ! Mol/m^2 added + endif if (deficit > 1.e-3_r8 * scale_factor_gasdiff) then if (deficit > 1.e-2_r8) then write(iulog,*)'Note: sink > source in ch4_tran, sources are changing '// & @@ -4038,7 +4594,7 @@ subroutine ch4_tran (bounds, & elseif (s == 2) then ! O2 ! Set rt, since it depends on conc - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) @@ -4047,18 +4603,56 @@ subroutine ch4_tran (bounds, & ! source_old could be removed later epsilon_t_old(c,j,s) = epsilon_t(c,j,s) ! epsilon_t acts like source also - dzj = dz(c,j) - if (j < nlevsoi .and. j == jwt(c)) then ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + + ! assign layer thickness + if (j == 0) then + dzj = surface_layer_thickness(c) + else + dzj = dz(c,j) + endif + + if (j == 0) then + if (sat == 0) then ! non-inundated conditions + if (jwt(c) > 0) then + ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & + dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + else ! jwt=0 (top soil layer saturated) + ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)*k_h_cc(c,j,s)) - & + dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif + endif + if (sat == 1) then ! inundated conditions + ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & + dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1)*k_h_cc(c,j-1,s))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif + + elseif (j < nlevsoi .and. j == jwt(c)) then ! concentration inside needs to be mult. by k_h_cc for dp1_zp1 term rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)*k_h_cc(c,j,s)) - & dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1))) + & 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) elseif (j < nlevsoi .and. j == jwt(c)+1) then ! concentration above needs to be mult. by k_h_cc for dm1_zm1 term - rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & - 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & - dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1)*k_h_cc(c,j-1,s))) + & - 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + if (sat==0) then + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & + dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1)*k_h_cc(c,j-1,s))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + else !j=1/jwt=0 -> no k_h applied + rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & + 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & + dm1_zm1(c,j) * (conc_o2_rel(c,j) -conc_o2_rel(c,j-1))) + & + 0.5_r8 * (source(c,j,s) + source_old(c,j,s)) + endif elseif (j < nlevsoi) then rt(c,j) = epsilon_t_old(c,j,s) / dtime_ch4 * conc_o2_rel(c,j) + & 0.5_r8 / dzj * (dp1_zp1(c,j) * (conc_o2_rel(c,j+1)-conc_o2_rel(c,j)) - & @@ -4080,16 +4674,16 @@ subroutine ch4_tran (bounds, & enddo ! fc; column enddo ! j; nlevsoi - call Tridiagonal(bounds, 0, nlevsoi, jtop(bounds%begc:bounds%endc), & + call Tridiagonal(bounds, -1, nlevsoi, jtop(bounds%begc:bounds%endc), & num_methc, filter_methc, & at(bounds%begc:bounds%endc, :), & bt(bounds%begc:bounds%endc, :), & ct(bounds%begc:bounds%endc, :), & rt(bounds%begc:bounds%endc, :), & - conc_o2_rel(bounds%begc:bounds%endc,0:nlevsoi)) + conc_o2_rel(bounds%begc:bounds%endc,:)) ! Ensure that concentrations stay above 0 - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) g = col%gridcell(c) @@ -4106,25 +4700,34 @@ subroutine ch4_tran (bounds, & enddo ! species ! Update absolute concentrations per unit volume - do j = 1,nlevsoi ! No need to update the atm. level concentrations + do j = 0,nlevsoi ! No need to update the atm. level concentrations do fc = 1, num_methc c = filter_methc (fc) - - conc_ch4(c,j) = conc_ch4_rel(c,j)*epsilon_t(c,j,1) - conc_o2(c,j) = conc_o2_rel(c,j) *epsilon_t(c,j,2) + if (j==0) then + surface_layer_conc_ch4(c) = conc_ch4_rel(c,j)*epsilon_t(c,j,1) + surface_layer_conc_o2(c) = conc_o2_rel(c,j)*epsilon_t(c,j,2) + else + conc_ch4(c,j) = conc_ch4_rel(c,j)*epsilon_t(c,j,1) + conc_o2(c,j) = conc_o2_rel(c,j) *epsilon_t(c,j,2) + endif end do end do ! Do Balance Check and absorb small ! discrepancy into surface flux. - do j = 1,nlevsoi + do j = 0,nlevsoi do fc = 1, num_methc c = filter_methc (fc) - if (j == 1) errch4(c) = 0._r8 - errch4(c) = errch4(c) + (conc_ch4(c,j) - conc_ch4_bef(c,j))*dz(c,j) - errch4(c) = errch4(c) - ch4_prod_depth(c,j)*dz(c,j)*dtime - errch4(c) = errch4(c) + ch4_oxid_depth(c,j)*dz(c,j)*dtime + if (j == 0) then + errch4(c) = 0._r8 + errch4(c) = errch4(c) + (surface_layer_conc_ch4(c) - conc_ch4_bef(c,j))*surface_layer_thickness(c) + errch4(c) = errch4(c) + surface_layer_ch4_oxid(c)*surface_layer_thickness(c)*dtime + else + errch4(c) = errch4(c) + (conc_ch4(c,j) - conc_ch4_bef(c,j))*dz(c,j) + errch4(c) = errch4(c) - ch4_prod_depth(c,j)*dz(c,j)*dtime + errch4(c) = errch4(c) + ch4_oxid_depth(c,j)*dz(c,j)*dtime + endif end do end do @@ -4133,14 +4736,19 @@ subroutine ch4_tran (bounds, & ! For history make sure that grnd_ch4_cond includes snow, for methane diffusivity grnd_ch4_cond(c) = spec_grnd_cond(c,1) - + errch4(c) = errch4(c) + (ch4_surf_aere(c) + ch4_surf_ebul(c) + ch4_surf_diff(c))*dtime if (abs(errch4(c)) < 1.e-8_r8) then ch4_surf_diff(c) = ch4_surf_diff(c) - errch4(c)/dtime + + ! add surface layer oxidation to top soil layer for ch4 conservation checking + ch4_oxid_depth(c,1) = ch4_oxid_depth(c,1) + surface_layer_ch4_oxid(c)*surface_layer_thickness(c)/dz(c,1) ! normalize by layer 1 thickness + else ! errch4 > 1e-8 mol / m^2 / timestep write(iulog,*)'CH4 Conservation Error in CH4Mod during diffusion, nstep, c, errch4 (mol /m^2.timestep)', & nstep,c,errch4(c) + g = col%gridcell(c) write(iulog,*)'Latdeg,Londeg=',grc%latdeg(g),grc%londeg(g) call endrun(subgrid_index=c, subgrid_level=subgrid_level_column, & @@ -4212,6 +4820,7 @@ subroutine get_jwt (bounds, num_methc, filter_methc, jwt, & if (jwt(c) == perch .and. h2osoi_vol(c,1) > f_sat * watsat(c,1)) then ! missed that the top layer is saturated jwt(c) = 0 endif + end do end associate @@ -4362,7 +4971,11 @@ subroutine ch4_totcolch4(bounds, num_nolakec, filter_nolakec, num_lakec, filter_ dz => col%dz , & ! Input: [real(r8) (:,:) ] layer thickness (m) (-nlevsno+1:nlevsoi) finundated => ch4_inst%finundated_col , & ! Input: [real(r8) (:) ] fractional inundated area in soil column (excluding dedicated wetland columns) conc_ch4_sat => ch4_inst%conc_ch4_sat_col , & ! Input: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) - conc_ch4_unsat => ch4_inst%conc_ch4_unsat_col & ! Input: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) + conc_ch4_unsat => ch4_inst%conc_ch4_unsat_col , & ! Input: [real(r8) (:,:) ] CH4 conc in each soil layer (mol/m3) (nlevsoi) + surface_layer_conc_ch4_sat => ch4_inst%conc_ch4_surface_layer_sat_col , & ! Input: [real(r8) (:) ] CH4 conc surface water layer (mol/m3) + surface_layer_conc_ch4_unsat => ch4_inst%conc_ch4_surface_layer_unsat_col , & ! Input: [real(r8) (:) ] CH4 conc surface soil layer (mol/m3) + surface_layer_thickness_sat => ch4_inst%surface_layer_thickness_sat_col , & ! Input: [real(r8) (:) ] surface water layer thickness (m) + surface_layer_thickness_unsat => ch4_inst%surface_layer_thickness_unsat_col & ! Input: [real(r8) (:) ] surface soil layer thickness (m) ) do fc = 1, num_nolakec @@ -4375,6 +4988,15 @@ subroutine ch4_totcolch4(bounds, num_nolakec, filter_nolakec, num_lakec, filter_ totcolch4(c) = 0._r8 end do + ! Add surface layer ch4 and convert mol CH4 to g C + do fc = 1, num_nolakec + c = filter_nolakec(fc) + totcolch4(c) = totcolch4(c) & + + (finundated(c)*surface_layer_conc_ch4_sat(c)*surface_layer_thickness_sat(c) & + + (1._r8-finundated(c))*surface_layer_conc_ch4_unsat(c)*surface_layer_thickness_unsat(c)) & + *catomw + end do + do j = 1, nlevsoi do fc = 1, num_nolakec c = filter_nolakec(fc) @@ -4383,11 +5005,13 @@ subroutine ch4_totcolch4(bounds, num_nolakec, filter_nolakec, num_lakec, filter_ ! urban landunits. if (lun%itype(l) .eq. isturb_tbd .or. lun%itype(l) .eq. isturb_hd .or. lun%itype(l) .eq. isturb_md) then totcolch4(c) = totcolch4(c) + & - (finundated(c)*conc_ch4_sat(c,j) + (1._r8-finundated(c))*conc_ch4_unsat(c,j)) * & + (finundated(c)*conc_ch4_sat(c,j) & + + (1._r8-finundated(c))*conc_ch4_unsat(c,j)) * & dzsoi_decomp(j)*catomw else totcolch4(c) = totcolch4(c) + & - (finundated(c)*conc_ch4_sat(c,j) + (1._r8-finundated(c))*conc_ch4_unsat(c,j)) * & + (finundated(c)*conc_ch4_sat(c,j) & + + (1._r8-finundated(c))*conc_ch4_unsat(c,j)) * & dz(c,j)*catomw end if ! mol CH4 --> g C diff --git a/src/biogeophys/HillslopeHydrologyMod.F90 b/src/biogeophys/HillslopeHydrologyMod.F90 index 8fccd762f0..1c0b107a28 100644 --- a/src/biogeophys/HillslopeHydrologyMod.F90 +++ b/src/biogeophys/HillslopeHydrologyMod.F90 @@ -53,6 +53,8 @@ module HillslopeHydrologyMod integer, private, parameter :: soil_profile_set_lowland_upland = 2 integer, private, parameter :: soil_profile_linear = 3 + real(r8), private, parameter :: floodplain_slope = 1e-3_r8 + !----------------------------------------------------------------------- contains @@ -962,17 +964,32 @@ subroutine HillslopeStreamOutflow(bounds, & integer :: c, l, g, i, j integer :: nstep - real(r8) :: dtime ! land model time step (sec) + real(r8) :: dtime ! land model time step (sec) + real(r8) :: cross_sectional_area_channel ! cross sectional area of stream channel (m2) real(r8) :: cross_sectional_area ! cross sectional area of stream water (m2) real(r8) :: stream_depth ! depth of stream water (m) real(r8) :: hydraulic_radius ! cross sectional area divided by wetted perimeter (m) real(r8) :: flow_velocity ! flow velocity (m/s) real(r8) :: overbank_area ! area of water above bankfull (m2) + real(r8) :: overbank_depth ! depth of water above bankfull (m) + real(r8) :: dynamic_viscosity + real(r8) :: bankfull_flow_velocity + real(r8) :: kinematic_viscosity + real(r8) :: length_scale + real(r8) :: reynolds real(r8), parameter :: manning_roughness = 0.03_r8 ! manning roughness real(r8), parameter :: manning_exponent = 0.667_r8 ! manning exponent integer, parameter :: overbank_method = 1 ! method to treat overbank stream storage; 1 = increase dynamic slope, 2 = increase flow area cross section, 3 = remove instantaneously logical :: active_stream + real(r8), parameter :: Acoef = 0.02939 !mPa·s + real(r8), parameter :: Bcoef = 507.88 !K + real(r8), parameter :: Ccoef = 149.30 !K + real(r8), parameter :: Tavg = 293.15 !20C + 273.15 + real(r8), parameter :: rhow = 1e3 + real(r8), parameter :: reynolds_thresh = 1000 + real(r8), parameter :: stream_conductivity = 1e-3_r8 + real(r8), parameter :: water_surface_slope = 5e-5_r8 character(len=*), parameter :: subname = 'HillslopeStreamOutflow' !----------------------------------------------------------------------- @@ -998,25 +1015,71 @@ subroutine HillslopeStreamOutflow(bounds, & if (lun%active(l) .and. active_stream) then ! Streamflow calculated from Manning equation if (streamflow_method == streamflow_manning) then +cross_sectional_area_channel = lun%stream_channel_width(l)*lun%stream_channel_depth(l) + cross_sectional_area = stream_water_volume(l) & /lun%stream_channel_length(l) - stream_depth = cross_sectional_area & - /lun%stream_channel_width(l) - hydraulic_radius = cross_sectional_area & - /(lun%stream_channel_width(l) + 2*stream_depth) + + ! calculate hydraulic radius + if(cross_sectional_area <= cross_sectional_area_channel) then + stream_depth = cross_sectional_area & + /lun%stream_channel_width(l) + overbank_depth = 0._r8 + hydraulic_radius = cross_sectional_area & + /(lun%stream_channel_width(l) + 2*stream_depth) + else ! overbank conditions exist + stream_depth = lun%stream_channel_depth(l) + overbank_area = cross_sectional_area & + -cross_sectional_area_channel + ! use small positive slope for floodplain + overbank_depth = sqrt(floodplain_slope*overbank_area) + + hydraulic_radius = cross_sectional_area & + /(lun%stream_channel_width(l) & + + 2*stream_depth & + +2*overbank_depth/floodplain_slope) ! sin ~ tan for small values + endif if (hydraulic_radius <= 0._r8) then volumetric_streamflow(l) = 0._r8 else - flow_velocity = (hydraulic_radius)**manning_exponent & + ! check reynolds number using bankfull channel properties (static) + bankfull_flow_velocity = (lun%stream_channel_depth(l))**manning_exponent & * sqrt(lun%stream_channel_slope(l)) & / manning_roughness + + ! check reynolds number dynamically + flow_velocity = (hydraulic_radius)**manning_exponent & + * sqrt(lun%stream_channel_slope(l)) & + / manning_roughness + + ! Vogel-Fulcher-Tammann equation + dynamic_viscosity = Acoef * exp((Bcoef/(Tavg-Ccoef))) + kinematic_viscosity = dynamic_viscosity/rhow + + ! let length scale equal bankfull depth + !length_scale = lun%stream_channel_depth(l) + !reynolds = bankfull_flow_velocity*length_scale/kinematic_viscosity + ! let length scale equal dynamic depth + length_scale = stream_depth + overbank_depth + reynolds = flow_velocity*length_scale/kinematic_viscosity + + if( reynolds < reynolds_thresh) then + ! use exponent = 2 instead of 0.66, and specified water surface slope + flow_velocity = (hydraulic_radius)**2 & + * water_surface_slope & + / (manning_roughness) + else + flow_velocity = (hydraulic_radius)**manning_exponent & + * sqrt(lun%stream_channel_slope(l)) & + / manning_roughness + endif ! overbank flow - if (stream_depth > lun%stream_channel_depth(l)) then + if (overbank_depth > 0._r8) then if (overbank_method == 1) then - ! try increasing dynamic slope - volumetric_streamflow(l) = cross_sectional_area * flow_velocity & - *(stream_depth/lun%stream_channel_depth(l)) + ! flow velocity already accounts for overbank conditions + + volumetric_streamflow(l) = cross_sectional_area * flow_velocity else if (overbank_method == 2) then ! try increasing flow area cross section overbank_area = (stream_depth -lun%stream_channel_depth(l)) * 30._r8 * lun%stream_channel_width(l) @@ -1080,6 +1143,7 @@ subroutine HillslopeUpdateStreamWater(bounds, waterstatebulk_inst, & real(r8) :: qflx_drain_vol ! volumetric saturated drainage (m3/s) real(r8) :: dtime ! land model time step (sec) logical :: active_stream + real(r8) :: overbank_area, overbank_depth character(len=*), parameter :: subname = 'HillslopeUpdateStreamWater' @@ -1138,6 +1202,14 @@ subroutine HillslopeUpdateStreamWater(bounds, waterstatebulk_inst, & /lun%stream_channel_length(l) & /lun%stream_channel_width(l) + ! recalculate for floodplain approximation + if (stream_water_depth(l)>lun%stream_channel_depth(l)) then + overbank_area = (stream_water_volume(l) & + /lun%stream_channel_length(l)) & + -lun%stream_channel_width(l)*lun%stream_channel_depth(l) + overbank_depth = sqrt(floodplain_slope*overbank_area) + stream_water_depth(l) = lun%stream_channel_depth(l) + overbank_depth + endif end if enddo diff --git a/src/biogeophys/HydrologyNoDrainageMod.F90 b/src/biogeophys/HydrologyNoDrainageMod.F90 index 8f294c652e..c8cea33c9e 100644 --- a/src/biogeophys/HydrologyNoDrainageMod.F90 +++ b/src/biogeophys/HydrologyNoDrainageMod.F90 @@ -322,7 +322,7 @@ subroutine HydrologyNoDrainage(bounds, & call UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & infiltration_excess_runoff_inst, & - energyflux_inst, soilhydrology_inst, & + energyflux_inst, soilhydrology_inst, soilstate_inst, & b_waterflux_inst, b_waterstate_inst, b_waterdiagnostic_inst) call Infiltration(bounds, num_hydrologyc, filter_hydrologyc, & diff --git a/src/biogeophys/SaturatedExcessRunoffMod.F90 b/src/biogeophys/SaturatedExcessRunoffMod.F90 index 9956a7dfb8..5349c34063 100644 --- a/src/biogeophys/SaturatedExcessRunoffMod.F90 +++ b/src/biogeophys/SaturatedExcessRunoffMod.F90 @@ -267,14 +267,17 @@ subroutine SaturatedExcessRunoff (this, bounds, num_hydrologyc, filter_hydrology endif ! ------------------------------------------------------------------------ - ! Set fsat to zero for upland hillslope columns + ! Set fsat to zero for hillslope columns + ! + ! hillslope explicitly represents lateral water movement and + ! convergence, so the implicit topmodel approach is unnecessary + ! ! ------------------------------------------------------------------------ if (hillslope_fsat_equals_zero) then do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if(col%is_hillslope_column(c) .and. col%active(c)) then - ! Set fsat to zero for upland columns - if (col%cold(c) /= ispval) fsat(c) = 0._r8 + fsat(c) = 0._r8 endif end do endif diff --git a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 index 5f38830891..7e30ee91a6 100644 --- a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 +++ b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 @@ -77,6 +77,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst use column_varcon , only : icol_shadewall, icol_road_perv, icol_road_imperv, icol_roof, icol_sunwall use fileutils , only : getfil use ncdio_pio , only : file_desc_t, ncd_io, ncd_pio_openfile, ncd_pio_closefile + use SurfaceWaterMod , only : pc_hillslope ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -99,6 +100,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst real(r8), pointer :: sandcol (:,:) ! column level sand fraction for calculating VIC parameters real(r8), pointer :: claycol (:,:) ! column level clay fraction for calculating VIC parameters real(r8), pointer :: om_fraccol (:,:) ! column level organic matter fraction for calculating VIC parameters + real(r8) :: pc ! column-level pc parameter !----------------------------------------------------------------------- ! Initialize VIC variables @@ -215,12 +217,19 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst ! determine h2osfc threshold ("fill & spill" concept) ! set to zero for no h2osfc (w/frac_infclust =large) - + + if(col%is_hillslope_column(c)) then + pc = max(pc_hillslope(col%hill_slope(c)), 0.01_r8) + else + ! use standard value + pc = params_inst%pc + endif + soilhydrology_inst%h2osfc_thresh_col(c) = 0._r8 if (micro_sigma(c) > 1.e-6_r8 .and. (soilhydrology_inst%h2osfcflag /= 0)) then d = 0.0_r8 do p = 1,4 - fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - params_inst%pc + fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - pc dfdd = exp(-d**2/(2.0_r8*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0_r8*shr_const_pi)) d = d - fd/dfdd enddo diff --git a/src/biogeophys/SurfaceWaterMod.F90 b/src/biogeophys/SurfaceWaterMod.F90 index 562c64cc18..89e1768a05 100644 --- a/src/biogeophys/SurfaceWaterMod.F90 +++ b/src/biogeophys/SurfaceWaterMod.F90 @@ -32,6 +32,7 @@ module SurfaceWaterMod ! !PUBLIC MEMBER FUNCTIONS: public :: UpdateFracH2oSfc ! Determine fraction of land surfaces which are submerged public :: UpdateH2osfc ! Calculate fluxes out of h2osfc and update the h2osfc state + public :: pc_hillslope ! Calculate pc parameter for hillslope columns public :: readParams ! !PRIVATE MEMBER FUNCTIONS: @@ -71,6 +72,21 @@ subroutine readParams( ncid ) end subroutine readParams + !----------------------------------------------------------------------- + real(r8) function pc_hillslope( hill_slope ) + !------------------------------------------------------------------------------ + ! + ! Calculate the pc parameter using hillslope geoparameter as input + ! + !------------------------------------------------------------------------------ + real(r8), intent(in) :: hill_slope ! Hillslope column slope + real(r8), parameter :: pc_a = 0.6_r8 + real(r8), parameter :: pc_b = -20._r8 + + pc_hillslope = pc_a*exp(pc_b*hill_slope) + + end function pc_hillslope + !----------------------------------------------------------------------- subroutine UpdateFracH2oSfc(bounds, num_soilc, filter_soilc, & water_inst) @@ -217,7 +233,8 @@ subroutine BulkDiag_FracH2oSfc(bounds, num_soilc, filter_soilc, & real(r8):: d,fd,dfdd ! temporary variable for frac_h2o iteration real(r8):: sigma ! microtopography pdf sigma in mm real(r8):: min_h2osfc - !----------------------------------------------------------------------- + real(r8):: f_overlap + !----------------------------------------------------------------------- SHR_ASSERT_FL((ubound(micro_sigma, 1) == bounds%endc), sourcefile, __LINE__) SHR_ASSERT_FL((ubound(h2osno_total, 1) == bounds%endc), sourcefile, __LINE__) @@ -276,7 +293,11 @@ subroutine BulkDiag_FracH2oSfc(bounds, num_soilc, filter_soilc, & if (frac_sno(c) > (1._r8 - frac_h2osfc(c)) .and. h2osno_total(c) > 0) then if (frac_h2osfc(c) > 0.01_r8) then - frac_h2osfc(c) = max(1.0_r8 - frac_sno(c),0.01_r8) + f_overlap = frac_sno(c) + frac_h2osfc(c) - 1._r8 + ! split overlap between snow and surface water rather than reducing only surface water + if (f_overlap > 0._r8) then + frac_h2osfc(c) = frac_h2osfc(c) - 0.5_r8*f_overlap + endif frac_sno(c) = 1.0_r8 - frac_h2osfc(c) else frac_sno(c) = 1.0_r8 - frac_h2osfc(c) @@ -344,17 +365,20 @@ end subroutine UpdateState_TooSmallH2osfcToSoil !----------------------------------------------------------------------- subroutine UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & infiltration_excess_runoff_inst, & - energyflux_inst, soilhydrology_inst, & + energyflux_inst, soilhydrology_inst, soilstate_inst, & waterfluxbulk_inst, waterstatebulk_inst, waterdiagnosticbulk_inst) ! ! !DESCRIPTION: ! Calculate fluxes out of h2osfc and update the h2osfc state ! + use SoilStateType , only : soilstate_type + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(infiltration_excess_runoff_type), intent(in) :: infiltration_excess_runoff_inst + type(soilstate_type) , intent(in) :: soilstate_inst type(energyflux_type) , intent(in) :: energyflux_inst type(soilhydrology_type) , intent(in) :: soilhydrology_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst @@ -378,6 +402,7 @@ subroutine UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & qflx_h2osfc_surf => waterfluxbulk_inst%qflx_h2osfc_surf_col , & ! Output: [real(r8) (:) ] surface water runoff (mm H2O /s) qflx_h2osfc_drain => waterfluxbulk_inst%qflx_h2osfc_drain_col , & ! Output: [real(r8) (:) ] bottom drainage from h2osfc (mm H2O /s) + icefrac => soilhydrology_inst%icefrac_col , & ! Input: [real(r8) (:,:) ] fraction of ice h2osfc_thresh => soilhydrology_inst%h2osfc_thresh_col, & ! Input: [real(r8) (:) ] level at which h2osfc "percolates" h2osfcflag => soilhydrology_inst%h2osfcflag & ! Input: integer ) @@ -412,6 +437,8 @@ subroutine UpdateH2osfc(bounds, num_hydrologyc, filter_hydrologyc, & h2osfc = h2osfc_partial(bounds%begc:bounds%endc), & frac_h2osfc = frac_h2osfc(bounds%begc:bounds%endc), & qinmax = qinmax(bounds%begc:bounds%endc), & + icefrac = icefrac(bounds%begc:bounds%endc,1), & + soilstate_inst = soilstate_inst, & qflx_h2osfc_drain = qflx_h2osfc_drain(bounds%begc:bounds%endc)) ! Update h2osfc based on fluxes @@ -456,6 +483,7 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & real(r8) :: dtime ! land model time step (sec) real(r8) :: frac_infclust ! fraction of submerged area that is connected real(r8) :: k_wet ! linear reservoir coefficient for h2osfc + real(r8) :: pc ! column-level pc parameter real(r8),parameter :: min_hill_slope = 1e-3_r8! minimum value of hillslope for outflow character(len=*), parameter :: subname = 'QflxH2osfcSurf' @@ -473,10 +501,17 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & c = filter_hydrologyc(fc) if (h2osfcflag==1) then - if (frac_h2osfc_nosnow(c) <= params_inst%pc) then + + if(col%is_hillslope_column(c)) then + pc = max(pc_hillslope(col%hill_slope(c)), 0.01_r8) + else + pc = params_inst%pc + endif + + if (frac_h2osfc_nosnow(c) <= pc) then frac_infclust=0.0_r8 else - frac_infclust=(frac_h2osfc_nosnow(c)-params_inst%pc)**params_inst%mu + frac_infclust=(frac_h2osfc_nosnow(c)-pc)**params_inst%mu endif endif @@ -485,8 +520,7 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & ! spatially variable k_wet k_wet=1.0e-4_r8 * sin((rpi/180._r8) * topo_slope(c)) if (col%is_hillslope_column(c)) then - ! require a minimum value to ensure non-zero outflow - k_wet = 1e-4_r8 * max(col%hill_slope(c),min_hill_slope) + k_wet = max(col%hill_slope(c),0._r8)/dtime endif qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c)) @@ -506,8 +540,8 @@ end subroutine QflxH2osfcSurf !----------------------------------------------------------------------- subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & - h2osfcflag, h2osfc, frac_h2osfc, qinmax, & - qflx_h2osfc_drain) + h2osfcflag, h2osfc, frac_h2osfc, qinmax, icefrac, & + soilstate_inst, qflx_h2osfc_drain) ! ! !DESCRIPTION: ! Compute qflx_h2osfc_drain @@ -515,6 +549,11 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & ! Note that, if h2osfc is negative, then qflx_h2osfc_drain will be negative - acting ! to exactly restore h2osfc to 0. ! + ! !USES: + use ColumnType , only: col + use SoilHydrologyMod, only: params_inst + use SoilStateType , only: soilstate_type + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter @@ -523,11 +562,19 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & real(r8) , intent(in) :: h2osfc( bounds%begc: ) ! surface water (mm) real(r8) , intent(in) :: frac_h2osfc( bounds%begc: ) ! fraction of ground covered by surface water (0 to 1) real(r8) , intent(in) :: qinmax( bounds%begc: ) ! maximum infiltration rate (mm H2O /s) + real(r8) , intent(in) :: icefrac( bounds%begc: ) ! soil ice fraction in top soil layer () + type(soilstate_type), intent(in) :: soilstate_inst real(r8) , intent(inout) :: qflx_h2osfc_drain( bounds%begc: ) ! bottom drainage from h2osfc (mm H2O /s) ! ! !LOCAL VARIABLES: integer :: fc, c real(r8) :: dtime ! land model time step (sec) + real(r8) :: ice_imped ! ice impedance () + real(r8) :: den ! temporary variable + real(r8) :: hk1 ! surface hydraulic conductivity (mm/s) + real(r8) :: smp1 ! minimum soil suction (mm) + + real(r8),parameter :: m_to_mm = 1.e3_r8 !convert meters to mm character(len=*), parameter :: subname = 'QflxH2osfcDrain' !----------------------------------------------------------------------- @@ -537,21 +584,43 @@ subroutine QflxH2osfcDrain(bounds, num_hydrologyc, filter_hydrologyc, & SHR_ASSERT_ALL_FL((ubound(qinmax) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(qflx_h2osfc_drain) == (/bounds%endc/)), sourcefile, __LINE__) - dtime = get_step_size_real() + associate(& + z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) + sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) + smp => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm] + hk => soilstate_inst%hk_l_col & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s) + ) ! end associate statement - do fc = 1, num_hydrologyc - c = filter_hydrologyc(fc) + dtime = get_step_size_real() - if (h2osfc(c) < 0.0) then - qflx_h2osfc_drain(c) = h2osfc(c)/dtime - else - qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax(c),h2osfc(c)/dtime) - if(h2osfcflag==0) then - ! ensure no h2osfc - qflx_h2osfc_drain(c)= max(0._r8,h2osfc(c)/dtime) - end if - end if - end do + do fc = 1, num_hydrologyc + c = filter_hydrologyc(fc) + + if (h2osfc(c) < 0.0) then + qflx_h2osfc_drain(c) = h2osfc(c)/dtime + else + + ! calculate surface soil layer hydraulic properties + ice_imped = 10._r8**(-params_inst%e_ice*icefrac(c)) + den = z(c,1)* m_to_mm + hk1 = frac_h2osfc(c)*ice_imped*hk(c,1) + smp1 = -sucsat(c,1) + + ! compute the flux into the top soil layer + qflx_h2osfc_drain(c) = -hk1*(smp(c,1) - smp1)/den + hk1 + + ! bound drainage values + qflx_h2osfc_drain(c)=min(frac_h2osfc(c)*qinmax(c),qflx_h2osfc_drain(c)) + qflx_h2osfc_drain(c)=min(h2osfc(c)/dtime,qflx_h2osfc_drain(c)) + + if(h2osfcflag==0) then + ! ensure no h2osfc + qflx_h2osfc_drain(c)= max(0._r8,h2osfc(c)/dtime) + end if + end if + end do + + end associate end subroutine QflxH2osfcDrain diff --git a/src/cpl/share_esmf/ch4FInundatedStreamType.F90 b/src/cpl/share_esmf/ch4FInundatedStreamType.F90 index 901f54b31f..0c0ec3c257 100644 --- a/src/cpl/share_esmf/ch4FInundatedStreamType.F90 +++ b/src/cpl/share_esmf/ch4FInundatedStreamType.F90 @@ -202,13 +202,15 @@ logical function UseStreams(this) ! Return true if ! ! !USES: + use ch4varcon , only : finundation_mtd_h2osfc ! ! !ARGUMENTS: implicit none class(ch4finundatedstream_type) :: this ! ! !LOCAL VARIABLES: - if ( trim(control%stream_fldFileName_ch4finundated) == '' )then + if ( trim(control%stream_fldFileName_ch4finundated) == '' & + .or. finundation_mtd == finundation_mtd_h2osfc )then UseStreams = .false. else UseStreams = .true. diff --git a/src/dyn_subgrid/dynColumnStateUpdaterMod.F90 b/src/dyn_subgrid/dynColumnStateUpdaterMod.F90 index b5373247ff..7dc11ff4b6 100644 --- a/src/dyn_subgrid/dynColumnStateUpdaterMod.F90 +++ b/src/dyn_subgrid/dynColumnStateUpdaterMod.F90 @@ -50,20 +50,40 @@ module dynColumnStateUpdaterMod ! update_column_state_fill_using_fixed_values. This method is appropriate when ! all special landunits contribute the same fixed value. ! - ! All update methods accept optional fractional areas. If provided, these give the - ! fraction of the column over which the state variable applies. For example, this can be - ! used if each column is split into an inundated and uninundated fraction; if the state - ! variable gives mass per unit area of the inundated portion of the column, then - ! fractional_area should specify the inundated fraction of each column. (If no - ! fractional areas are provided, then the state variable is assumed to apply to the - ! entire column, which is the typical case.) If fractional areas are provided, both - ! 'old' (before area updates) and 'new' (after area updates) versions must be provided. - ! The determination of the new fractional areas should be done via a call to - ! update_column_state_no_special_handling: the code for working with fractional areas - ! has no capability to do any special handling, and there could be inconsistencies if - ! the new fractional areas were determined using a method that has some special - ! handling. Thus, fractional_area must be valid for all columns, even if the state - ! variable itself is filled using one of the special-handling methods here. + ! All update methods accept optional volume multipliers. If provided, these specify any + ! extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that vary in + ! space or time for the given state variable (in addition to col%wtgcell, which is + ! already accounted for). This is needed to convert variables that have units of per-m3 + ! or per-m2 to the actual masses that are needed for conservation. (The elements of + ! volume or area that are constant in space and time - such as col%dz - do not need to + ! be included here, although it doesn't hurt to include them.) Here are some examples of + ! the use of these volume multipliers: + ! + ! - An additional area multiplier. e.g., if each column is split into an inundated + ! and uninundated fraction: For a state variable that gives mass per unit volume of + ! the inundated portion of the column (and which uses level thicknesses that are + ! constant in space and time), volume_multiplier should specify the inundated + ! fraction of each column. + ! + ! - A level thickness: If the given state variable applies to a vertical level that + ! has a variable thickness, volume_multiplier should specify the thickness of this + ! level for each column. + ! + ! - If there is more than one multiplier that applies to the given variable, these + ! should be combined multiplicatively. e.g., for a state variable that uses + ! variable level thicknesses and applies just over the inundated portion of the + ! column, volume_multiplier should be the product of the inundated fraction of each + ! column and the thickness of this level for each column. + ! + ! If volume multipliers are provided, both 'old' (before area updates) and 'new' (after + ! area updates) versions must be provided. Elements of the volume multipliers don't + ! necessarily need to be updated with column area changes, but if they are updated, this + ! update should be done via a call to update_column_state_no_special_handling: the code + ! for working with volume multipliers has no capability to do any special handling, and + ! there could be inconsistencies if the new volume multipliers were determined using a + ! method that has some special handling. Thus, volume_multiplier must be valid for all + ! columns, even if the state variable itself is filled using one of the special-handling + ! methods here. ! ! For methods other than update_column_state_no_special_handling, an additional inout ! argument (non_conserved_mass_grc) accumulates the non-conserved mass due to shrinking @@ -77,9 +97,9 @@ module dynColumnStateUpdaterMod ! All methods accept an optional output argument (adjustment) which gives the apparent ! state adjustment in each column - that is, the per-column-area value after area ! changes minus the per-column-area value before area changes. This adjustment will be 0 - ! for shrinking columns, and can be either positive or negative for growing - ! columns. If there are fractional areas, then this adjustment is defined as (val_new * - ! fractional_area_new - val_old * fractional_area_old). + ! for shrinking columns, and can be either positive or negative for growing columns. If + ! there are volume multipliers, then this adjustment is defined as (val_new * + ! volume_multiplier_new - val_old * volume_multiplier_old). ! ! NOTE(wjs, 2017-02-24) The implementation involves a few levels of calls to get to the ! real work routine (update_column_state). This design made sense (in terms of removing @@ -142,9 +162,9 @@ module dynColumnStateUpdaterMod ! Private routines ! intermediate routine between the public routines and the real work routine - ! (update_column_state); this routine determines the fractional_areas to use in + ! (update_column_state); this routine determines the volume_multipliers to use in ! update_column_state - procedure, private :: update_column_state_with_optional_fractions + procedure, private :: update_column_state_with_optional_multipliers ! do the work of updating a column state procedure, private :: update_column_state @@ -270,7 +290,7 @@ end subroutine set_new_weights !----------------------------------------------------------------------- subroutine update_column_state_no_special_handling(this, bounds, clump_index, & - var, fractional_area_old, fractional_area_new, adjustment) + var, volume_multiplier_old, volume_multiplier_new, adjustment) ! ! !DESCRIPTION: ! Adjust the values of a column-level state variable due to changes in subgrid @@ -290,14 +310,15 @@ subroutine update_column_state_no_special_handling(this, bounds, clump_index, & real(r8), intent(inout) :: var( bounds%begc: ) ! column-level variable - ! Fraction of each column over which the state variable applies. See module-level - ! documentation for details. You must provide both old & new fractional areas, or - ! neither: it is invalid to provide just one. Fractional areas should be valid for all - ! columns, and fractional_area_new should have been computed based on a call to - ! update_column_state_no_special_handling: code that works with these fractional areas - ! is not able to do any special handling. - real(r8), optional, intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), optional, intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable. See module-level documentation + ! for details. You must provide both old & new volume multipliers, or neither: it is + ! invalid to provide just one. Volume multipliers should be valid for all columns, and + ! any elements of volume_multiplier_new that change due to area changes should have + ! been computed based on a call to update_column_state_no_special_handling: code that + ! works with these volume multipliers is not able to do any special handling. + real(r8), optional, intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), optional, intent(in) :: volume_multiplier_new( bounds%begc: ) ! Apparent state adjustment in each column real(r8), optional, intent(out) :: adjustment( bounds%begc: ) @@ -330,17 +351,17 @@ subroutine update_column_state_no_special_handling(this, bounds, clump_index, & non_conserved_mass(bounds%begg:bounds%endg) = 0._r8 ! explicit bounds not needed on any of these arguments - and specifying explicit - ! bounds defeats some later bounds checking (for fractional_area_old and - ! fractional_area_new) - call this%update_column_state_with_optional_fractions(& + ! bounds defeats some later bounds checking (for volume_multiplier_old and + ! volume_multiplier_new) + call this%update_column_state_with_optional_multipliers(& bounds = bounds, & vals_input = vals_input, & vals_input_valid = vals_input_valid, & has_prognostic_state = has_prognostic_state, & var = var, & non_conserved_mass = non_conserved_mass, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) ! Since there is no special handling in this routine, the non_conserved_mass variable @@ -355,7 +376,7 @@ end subroutine update_column_state_no_special_handling !----------------------------------------------------------------------- subroutine update_column_state_fill_special_using_natveg(this, bounds, clump_index, & - var, non_conserved_mass_grc, fractional_area_old, fractional_area_new, adjustment) + var, non_conserved_mass_grc, volume_multiplier_old, volume_multiplier_new, adjustment) ! ! !DESCRIPTION: ! Adjust the values of a column-level state variable due to changes in subgrid @@ -392,14 +413,15 @@ subroutine update_column_state_fill_special_using_natveg(this, bounds, clump_ind ! out appropriately. real(r8) , intent(inout) :: non_conserved_mass_grc( bounds%begg: ) - ! Fraction of each column over which the state variable applies. See module-level - ! documentation for details. You must provide both old & new fractional areas, or - ! neither: it is invalid to provide just one. Fractional areas should be valid for all - ! columns, and fractional_area_new should have been computed based on a call to - ! update_column_state_no_special_handling: code that works with these fractional areas - ! is not able to do any special handling. - real(r8), optional, intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), optional, intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable. See module-level documentation + ! for details. You must provide both old & new volume multipliers, or neither: it is + ! invalid to provide just one. Volume multipliers should be valid for all columns, and + ! any elements of volume_multiplier_new that change due to area changes should have + ! been computed based on a call to update_column_state_no_special_handling: code that + ! works with these volume multipliers is not able to do any special handling. + real(r8), optional, intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), optional, intent(in) :: volume_multiplier_new( bounds%begc: ) ! Apparent state adjustment in each column real(r8), optional, intent(out) :: adjustment( bounds%begc: ) @@ -447,17 +469,17 @@ subroutine update_column_state_fill_special_using_natveg(this, bounds, clump_ind end do ! explicit bounds not needed on any of these arguments - and specifying explicit - ! bounds defeats some later bounds checking (for fractional_area_old and - ! fractional_area_new) - call this%update_column_state_with_optional_fractions(& + ! bounds defeats some later bounds checking (for volume_multiplier_old and + ! volume_multiplier_new) + call this%update_column_state_with_optional_multipliers(& bounds = bounds, & vals_input = vals_input, & vals_input_valid = vals_input_valid, & has_prognostic_state = has_prognostic_state, & var = var, & non_conserved_mass = non_conserved_mass_grc, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) end if @@ -468,7 +490,7 @@ end subroutine update_column_state_fill_special_using_natveg subroutine update_column_state_fill_using_fixed_values(this, bounds, clump_index, & var, & landunit_values, non_conserved_mass_grc, & - fractional_area_old, fractional_area_new, & + volume_multiplier_old, volume_multiplier_new, & adjustment) ! ! !DESCRIPTION: @@ -521,14 +543,15 @@ subroutine update_column_state_fill_using_fixed_values(this, bounds, clump_index ! appropriately. real(r8) , intent(inout) :: non_conserved_mass_grc( bounds%begg: ) - ! Fraction of each column over which the state variable applies. See module-level - ! documentation for details. You must provide both old & new fractional areas, or - ! neither: it is invalid to provide just one. Fractional areas should be valid for all - ! columns, and fractional_area_new should have been computed based on a call to - ! update_column_state_no_special_handling: code that works with these fractional areas - ! is not able to do any special handling. - real(r8), optional, intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), optional, intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable. See module-level documentation + ! for details. You must provide both old & new volume multipliers, or neither: it is + ! invalid to provide just one. Volume multipliers should be valid for all columns, and + ! any elements of volume_multiplier_new that change due to area changes should have + ! been computed based on a call to update_column_state_no_special_handling: code that + ! works with these volume multipliers is not able to do any special handling. + real(r8), optional, intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), optional, intent(in) :: volume_multiplier_new( bounds%begc: ) ! Apparent state adjustment in each column real(r8), optional, intent(out) :: adjustment( bounds%begc: ) @@ -576,17 +599,17 @@ subroutine update_column_state_fill_using_fixed_values(this, bounds, clump_index end do ! explicit bounds not needed on any of these arguments - and specifying explicit - ! bounds defeats some later bounds checking (for fractional_area_old and - ! fractional_area_new) - call this%update_column_state_with_optional_fractions( & + ! bounds defeats some later bounds checking (for volume_multiplier_old and + ! volume_multiplier_new) + call this%update_column_state_with_optional_multipliers( & bounds = bounds, & vals_input = vals_input, & vals_input_valid = vals_input_valid, & has_prognostic_state = has_prognostic_state, & var = var, & non_conserved_mass = non_conserved_mass_grc, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) end if @@ -597,7 +620,7 @@ end subroutine update_column_state_fill_using_fixed_values subroutine update_column_state_fill_special_using_fixed_value(this, bounds, clump_index, & var, & special_value, non_conserved_mass_grc, & - fractional_area_old, fractional_area_new, & + volume_multiplier_old, volume_multiplier_new, & adjustment) ! ! !DESCRIPTION: @@ -632,14 +655,15 @@ subroutine update_column_state_fill_special_using_fixed_value(this, bounds, clum ! out appropriately. real(r8) , intent(inout) :: non_conserved_mass_grc( bounds%begg: ) - ! Fraction of each column over which the state variable applies. See module-level - ! documentation for details. You must provide both old & new fractional areas, or - ! neither: it is invalid to provide just one. Fractional areas should be valid for all - ! columns, and fractional_area_new should have been computed based on a call to - ! update_column_state_no_special_handling: code that works with these fractional areas - ! is not able to do any special handling. - real(r8), optional, intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), optional, intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable. See module-level documentation + ! for details. You must provide both old & new volume multipliers, or neither: it is + ! invalid to provide just one. Volume multipliers should be valid for all columns, and + ! any elements of volume_multiplier_new that change due to area changes should have + ! been computed based on a call to update_column_state_no_special_handling: code that + ! works with these volume multipliers is not able to do any special handling. + real(r8), optional, intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), optional, intent(in) :: volume_multiplier_new( bounds%begc: ) ! Apparent state adjustment in each column real(r8), optional, intent(out) :: adjustment( bounds%begc: ) @@ -665,8 +689,8 @@ subroutine update_column_state_fill_special_using_fixed_value(this, bounds, clum var = var, & landunit_values = landunit_values, & non_conserved_mass_grc = non_conserved_mass_grc, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) end subroutine update_column_state_fill_special_using_fixed_value @@ -678,15 +702,15 @@ end subroutine update_column_state_fill_special_using_fixed_value ! ======================================================================== !----------------------------------------------------------------------- - subroutine update_column_state_with_optional_fractions(this, bounds, & + subroutine update_column_state_with_optional_multipliers(this, bounds, & vals_input, vals_input_valid, has_prognostic_state, & var, non_conserved_mass, & - fractional_area_old, fractional_area_new, & + volume_multiplier_old, volume_multiplier_new, & adjustment) ! ! !DESCRIPTION: ! Intermediate routine between the public routines and the real work routine - ! (update_column_state). This routine determines the fractional areas to use in the + ! (update_column_state). This routine determines the volume multipliers to use in the ! call to update_column_state, and then does the call to update_column_state. ! ! !USES: @@ -717,45 +741,46 @@ subroutine update_column_state_with_optional_fractions(this, bounds, & ! the grid cell, negative denotes mass gained by the grid cell. real(r8), intent(inout) :: non_conserved_mass( bounds%begg: ) - ! Fraction of each column over which the state variable applies. See module-level - ! documentation for details. You must provide both old & new fractional areas, or - ! neither: it is invalid to provide just one. Fractional areas should be valid for all - ! columns, and fractional_area_new should have been computed based on a call to - ! update_column_state_no_special_handling: code that works with these fractional areas - ! is not able to do any special handling. - real(r8), optional, intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), optional, intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable. See module-level documentation + ! for details. You must provide both old & new volume multipliers, or neither: it is + ! invalid to provide just one. Volume multipliers should be valid for all columns, and + ! any elements of volume_multiplier_new that change due to area changes should have + ! been computed based on a call to update_column_state_no_special_handling: code that + ! works with these volume multipliers is not able to do any special handling. + real(r8), optional, intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), optional, intent(in) :: volume_multiplier_new( bounds%begc: ) ! Apparent state adjustment in each column real(r8), optional, intent(inout) :: adjustment( bounds%begc: ) ! ! !LOCAL VARIABLES: - real(r8) :: my_fractional_area_old(bounds%begc:bounds%endc) - real(r8) :: my_fractional_area_new(bounds%begc:bounds%endc) + real(r8) :: my_volume_multiplier_old(bounds%begc:bounds%endc) + real(r8) :: my_volume_multiplier_new(bounds%begc:bounds%endc) - character(len=*), parameter :: subname = 'update_column_state_with_optional_fractions' + character(len=*), parameter :: subname = 'update_column_state_with_optional_multipliers' !----------------------------------------------------------------------- - if (present(fractional_area_old) .and. .not. present(fractional_area_new)) then - call endrun(subname//' ERROR: If fractional_area_old is provided, then fractional_area_new must be provided, too') + if (present(volume_multiplier_old) .and. .not. present(volume_multiplier_new)) then + call endrun(subname//' ERROR: If volume_multiplier_old is provided, then volume_multiplier_new must be provided, too') end if - if (present(fractional_area_new) .and. .not. present(fractional_area_old)) then - call endrun(subname//' ERROR: If fractional_area_new is provided, then fractional_area_old must be provided, too') + if (present(volume_multiplier_new) .and. .not. present(volume_multiplier_old)) then + call endrun(subname//' ERROR: If volume_multiplier_new is provided, then volume_multiplier_old must be provided, too') end if - if (present(fractional_area_old)) then - SHR_ASSERT_ALL_FL((ubound(fractional_area_old) == (/bounds%endc/)), sourcefile, __LINE__) - my_fractional_area_old(bounds%begc:bounds%endc) = fractional_area_old(bounds%begc:bounds%endc) + if (present(volume_multiplier_old)) then + SHR_ASSERT_ALL_FL((ubound(volume_multiplier_old) == (/bounds%endc/)), sourcefile, __LINE__) + my_volume_multiplier_old(bounds%begc:bounds%endc) = volume_multiplier_old(bounds%begc:bounds%endc) else - my_fractional_area_old(bounds%begc:bounds%endc) = 1._r8 + my_volume_multiplier_old(bounds%begc:bounds%endc) = 1._r8 end if - if (present(fractional_area_new)) then - SHR_ASSERT_ALL_FL((ubound(fractional_area_new) == (/bounds%endc/)), sourcefile, __LINE__) - my_fractional_area_new(bounds%begc:bounds%endc) = fractional_area_new(bounds%begc:bounds%endc) + if (present(volume_multiplier_new)) then + SHR_ASSERT_ALL_FL((ubound(volume_multiplier_new) == (/bounds%endc/)), sourcefile, __LINE__) + my_volume_multiplier_new(bounds%begc:bounds%endc) = volume_multiplier_new(bounds%begc:bounds%endc) else - my_fractional_area_new(bounds%begc:bounds%endc) = 1._r8 + my_volume_multiplier_new(bounds%begc:bounds%endc) = 1._r8 end if call this%update_column_state(& @@ -763,19 +788,19 @@ subroutine update_column_state_with_optional_fractions(this, bounds, & vals_input = vals_input(bounds%begc:bounds%endc), & vals_input_valid = vals_input_valid(bounds%begc:bounds%endc), & has_prognostic_state = has_prognostic_state(bounds%begc:bounds%endc), & - fractional_area_old = my_fractional_area_old(bounds%begc:bounds%endc), & - fractional_area_new = my_fractional_area_new(bounds%begc:bounds%endc), & + volume_multiplier_old = my_volume_multiplier_old(bounds%begc:bounds%endc), & + volume_multiplier_new = my_volume_multiplier_new(bounds%begc:bounds%endc), & var = var(bounds%begc:bounds%endc), & non_conserved_mass = non_conserved_mass(bounds%begg:bounds%endg), & adjustment = adjustment) - end subroutine update_column_state_with_optional_fractions + end subroutine update_column_state_with_optional_multipliers !----------------------------------------------------------------------- subroutine update_column_state(this, bounds, & vals_input, vals_input_valid, has_prognostic_state, & - fractional_area_old, fractional_area_new, & + volume_multiplier_old, volume_multiplier_new, & var, non_conserved_mass, adjustment) ! ! !DESCRIPTION: @@ -800,10 +825,11 @@ subroutine update_column_state(this, bounds, & ! determines whether it can accept mass of this variable) logical, intent(in) :: has_prognostic_state( bounds%begc: ) - ! Fraction of each column over which the state variable applies, for both the old and - ! new subgrid weights - real(r8), intent(in) :: fractional_area_old( bounds%begc: ) - real(r8), intent(in) :: fractional_area_new( bounds%begc: ) + ! Any extra pieces of volume (for 3-d quantities) or area (for 2-d quantities) that + ! vary in space or time for the given state variable, for both the old and new subgrid + ! weights + real(r8), intent(in) :: volume_multiplier_old( bounds%begc: ) + real(r8), intent(in) :: volume_multiplier_new( bounds%begc: ) ! column-level variable of interest, updated in-place real(r8), intent(inout) :: var( bounds%begc: ) @@ -859,8 +885,8 @@ subroutine update_column_state(this, bounds, & SHR_ASSERT_ALL_FL((ubound(var) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(vals_input) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(has_prognostic_state) == (/bounds%endc/)), sourcefile, __LINE__) - SHR_ASSERT_ALL_FL((ubound(fractional_area_old) == (/bounds%endc/)), sourcefile, __LINE__) - SHR_ASSERT_ALL_FL((ubound(fractional_area_new) == (/bounds%endc/)), sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(volume_multiplier_old) == (/bounds%endc/)), sourcefile, __LINE__) + SHR_ASSERT_ALL_FL((ubound(volume_multiplier_new) == (/bounds%endc/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(non_conserved_mass) == (/bounds%endg/)), sourcefile, __LINE__) if (present(adjustment)) then SHR_ASSERT_ALL_FL((ubound(adjustment) == (/bounds%endc/)), sourcefile, __LINE__) @@ -895,7 +921,7 @@ subroutine update_column_state(this, bounds, & end if area_lost = -1._r8 * this%area_gained_col(c) total_area_lost_grc(g) = total_area_lost_grc(g) + area_lost - area_weighted_loss = area_lost * vals_input(c) * fractional_area_old(c) + area_weighted_loss = area_lost * vals_input(c) * volume_multiplier_old(c) total_loss_grc(g) = total_loss_grc(g) + area_weighted_loss if (.not. has_prognostic_state(c)) then @@ -926,19 +952,19 @@ subroutine update_column_state(this, bounds, & if (has_prognostic_state(c)) then val_old = var(c) - ! Need to make sure fractional_area_new /= 0 to avoid divide-by-zero. Note - ! that fractional_area_new == 0 can only happen if both - ! fractional_area_old(c) == 0 and the fractional_areas of the shrinking + ! Need to make sure volume_multiplier_new /= 0 to avoid divide-by-zero. Note + ! that volume_multiplier_new == 0 can only happen if both + ! volume_multiplier_old(c) == 0 and the volume_multipliers of the shrinking ! columns were all 0 - in which case the value of var is irrelevant for ! conservation purposes. - if (fractional_area_new(c) /= 0._r8) then - var(c) = (this%cwtgcell_old(c) * var(c) * fractional_area_old(c) + mass_gained) / & - (this%cwtgcell_new(c) * fractional_area_new(c)) + if (volume_multiplier_new(c) /= 0._r8) then + var(c) = (this%cwtgcell_old(c) * var(c) * volume_multiplier_old(c) + mass_gained) / & + (this%cwtgcell_new(c) * volume_multiplier_new(c)) end if if (present(adjustment)) then - adjustment(c) = var(c) * fractional_area_new(c) - & - val_old * fractional_area_old(c) + adjustment(c) = var(c) * volume_multiplier_new(c) - & + val_old * volume_multiplier_old(c) end if else non_conserved_mass(g) = non_conserved_mass(g) + mass_gained diff --git a/src/dyn_subgrid/test/dynColumnStateUpdater_test/test_column_state_updater.pf b/src/dyn_subgrid/test/dynColumnStateUpdater_test/test_column_state_updater.pf index 0758b84ed2..be067878e7 100644 --- a/src/dyn_subgrid/test/dynColumnStateUpdater_test/test_column_state_updater.pf +++ b/src/dyn_subgrid/test/dynColumnStateUpdater_test/test_column_state_updater.pf @@ -835,19 +835,19 @@ contains ! ------------------------------------------------------------------------ - ! Tests with fractional_areas + ! Tests with volume_multipliers ! ------------------------------------------------------------------------ @Test - subroutine withFractionalAreas_oneShrinksOneGrows(this) - ! This tests the updating of a variable with fractional area terms, with one column + subroutine withVolumeMultipliers_oneShrinksOneGrows(this) + ! This tests the updating of a variable with volume multiplier terms, with one column ! shrinking, one growing, and two staying the same in area class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: cwtlunit_old(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] real(r8), parameter :: cwtlunit_new(4) = [0.1_r8, 0.3_r8, 0.2_r8, 0.4_r8] - real(r8), parameter :: fractional_area_old(4) = [0.0_r8, 0.1_r8, 0.7_r8, 0.5_r8] - real(r8) :: myvar(4), myvar_orig(4), expected(4), fractional_area_new(4) + real(r8), parameter :: volume_multiplier_old(4) = [0.0_r8, 0.1_r8, 0.7_r8, 0.5_r8] + real(r8) :: myvar(4), myvar_orig(4), expected(4), volume_multiplier_new(4) real(r8) :: myvar_orig_weighted(4), myvar_weighted(4) ! Setup @@ -860,37 +860,37 @@ contains myvar = [2._r8, 3._r8, 7._r8, 8._r8] myvar_orig = myvar expected = myvar - expected(2) = (this%cwtgcell_old(2) * fractional_area_old(2) * myvar(2) + & - this%cwtgcell_new_minus_old(2) * fractional_area_old(3) * myvar(3)) / & - (this%cwtgcell_old(2)*fractional_area_old(2) + & - this%cwtgcell_new_minus_old(2) * fractional_area_old(3)) + expected(2) = (this%cwtgcell_old(2) * volume_multiplier_old(2) * myvar(2) + & + this%cwtgcell_new_minus_old(2) * volume_multiplier_old(3) * myvar(3)) / & + (this%cwtgcell_old(2)*volume_multiplier_old(2) + & + this%cwtgcell_new_minus_old(2) * volume_multiplier_old(3)) - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=myvar, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new) + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new) ! Verify @assertEqual(expected, myvar, tolerance=tol) - myvar_orig_weighted = myvar_orig*fractional_area_old - myvar_weighted = myvar*fractional_area_new - call this%assertConservation(myvar_orig_weighted, myvar_weighted, 'withFractionalAreas_oneShrinksOneGrows') - end subroutine withFractionalAreas_oneShrinksOneGrows + myvar_orig_weighted = myvar_orig*volume_multiplier_old + myvar_weighted = myvar*volume_multiplier_new + call this%assertConservation(myvar_orig_weighted, myvar_weighted, 'withVolumeMultipliers_oneShrinksOneGrows') + end subroutine withVolumeMultipliers_oneShrinksOneGrows @Test - subroutine withFractionalAreas_fraction0_oneShrinksOneGrows(this) - ! This tests the updating of a variable with fractional area terms, with one column + subroutine withVolumeMultipliers_multiplier0_oneShrinksOneGrows(this) + ! This tests the updating of a variable with volume multiplier terms, with one column ! shrinking, one growing, and two staying the same in area. Here the growing column - ! has a new fractional area of 0. + ! has a new volume multiplier of 0. class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: cwtlunit_old(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] real(r8), parameter :: cwtlunit_new(4) = [0.1_r8, 0.3_r8, 0.2_r8, 0.4_r8] - real(r8), parameter :: fractional_area_old(4) = [0.0_r8, 0.0_r8, 0.0_r8, 0.5_r8] - real(r8) :: myvar(4), myvar_orig(4), expected(4), fractional_area_new(4) + real(r8), parameter :: volume_multiplier_old(4) = [0.0_r8, 0.0_r8, 0.0_r8, 0.5_r8] + real(r8) :: myvar(4), myvar_orig(4), expected(4), volume_multiplier_new(4) real(r8) :: myvar_orig_weighted(4), myvar_weighted(4) ! Setup @@ -902,35 +902,35 @@ contains myvar = [2._r8, 3._r8, 7._r8, 8._r8] myvar_orig = myvar - ! If the growing column has a new fractional area of 0, then we expect no change in + ! If the growing column has a new volume multiplier of 0, then we expect no change in ! my_var. expected = myvar - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=myvar, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new) + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new) ! Verify @assertEqual(expected, myvar, tolerance=tol) - end subroutine withFractionalAreas_fraction0_oneShrinksOneGrows + end subroutine withVolumeMultipliers_multiplier0_oneShrinksOneGrows @Test - subroutine withFractionalAreas_twoShrinkTwoGrow(this) - ! This tests the updating of a variable with fractional area terms, with one column + subroutine withVolumeMultipliers_twoShrinkTwoGrow(this) + ! This tests the updating of a variable with volume multiplier terms, with one column ! shrinking to 0 area, one shrinking to non-zero area, one growing from 0 area, and ! one growing from non-zero area class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: cwtlunit_old(4) = [0.1_r8, 0.2_r8, 0._r8, 0.7_r8] real(r8), parameter :: cwtlunit_new(4) = [0._r8 , 0.15_r8, 0.08_r8, 0.77_r8] - real(r8), parameter :: fractional_area_old(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] + real(r8), parameter :: volume_multiplier_old(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] real(r8), parameter :: new_minus_old(4) = cwtlunit_new - cwtlunit_old real(r8) :: mass_gain_per_unit_area - real(r8) :: myvar(4), myvar_orig(4), expected(4), fractional_area_new(4) + real(r8) :: myvar(4), myvar_orig(4), expected(4), volume_multiplier_new(4) real(r8) :: myvar_orig_weighted(4), myvar_weighted(4) ! Setup @@ -943,40 +943,40 @@ contains myvar = [2._r8, 3._r8, 4._r8, 5._r8] myvar_orig = myvar - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=myvar, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new) + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new) ! Verify expected(1:2) = myvar_orig(1:2) mass_gain_per_unit_area = ( & - new_minus_old(1) * fractional_area_old(1) * myvar_orig(1) + & - new_minus_old(2) * fractional_area_old(2) * myvar_orig(2)) / & + new_minus_old(1) * volume_multiplier_old(1) * myvar_orig(1) + & + new_minus_old(2) * volume_multiplier_old(2) * myvar_orig(2)) / & (new_minus_old(1) + new_minus_old(2)) - expected(3) = mass_gain_per_unit_area / fractional_area_new(3) - expected(4) = (cwtlunit_old(4) * fractional_area_old(4) * myvar_orig(4) + & + expected(3) = mass_gain_per_unit_area / volume_multiplier_new(3) + expected(4) = (cwtlunit_old(4) * volume_multiplier_old(4) * myvar_orig(4) + & new_minus_old(4) * mass_gain_per_unit_area) / & - (cwtlunit_new(4) * fractional_area_new(4)) + (cwtlunit_new(4) * volume_multiplier_new(4)) @assertEqual(expected, myvar, tolerance=tol) - myvar_orig_weighted = myvar_orig*fractional_area_old - myvar_weighted = myvar*fractional_area_new - call this%assertConservation(myvar_orig_weighted, myvar_weighted, 'withFractionalAreas_twoShrinkTwoGrow') - end subroutine withFractionalAreas_twoShrinkTwoGrow + myvar_orig_weighted = myvar_orig*volume_multiplier_old + myvar_weighted = myvar*volume_multiplier_new + call this%assertConservation(myvar_orig_weighted, myvar_weighted, 'withVolumeMultipliers_twoShrinkTwoGrow') + end subroutine withVolumeMultipliers_twoShrinkTwoGrow @Test - subroutine nonConservedMass_withFractionalAreas_specialShrinks(this) - ! This tests the computation of non_conserved_mass with fractional area terms, when + subroutine nonConservedMass_withVolumeMultipliers_specialShrinks(this) + ! This tests the computation of non_conserved_mass with volume multiplier terms, when ! a special column shrinks class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: input_from_special = 11._r8 - real(r8), parameter :: fractional_area_old(3) = [0.1_r8, 0.2_r8, 0.3_r8] - real(r8) :: fractional_area_new(3) + real(r8), parameter :: volume_multiplier_old(3) = [0.1_r8, 0.2_r8, 0.3_r8] + real(r8) :: volume_multiplier_new(3) real(r8) :: myvar(3), myvar_orig(3) real(r8) :: expected_non_conserved_mass real(r8) :: non_conserved_mass_grc(1) = 0._r8 @@ -988,40 +988,40 @@ contains myvar = [2._r8, 3._r8, 4._r8] myvar_orig = myvar - expected_non_conserved_mass = this%cwtgcell_new_minus_old(3) * fractional_area_old(3) & + expected_non_conserved_mass = this%cwtgcell_new_minus_old(3) * volume_multiplier_old(3) & * input_from_special - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_fill_special_using_fixed_value(bounds, clump_index=1, var=myvar, & special_value=input_from_special, non_conserved_mass_grc=non_conserved_mass_grc, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new) + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new) ! Verify @assertEqual(expected_non_conserved_mass, non_conserved_mass_grc(1), tolerance=tol) ! Note that the specially-treated column is excluded from the conservation check: the ! value in that column is ignored for conservation. - myvar_orig_weighted = myvar_orig * fractional_area_old - myvar_weighted = myvar * fractional_area_new + myvar_orig_weighted = myvar_orig * volume_multiplier_old + myvar_weighted = myvar * volume_multiplier_new call this%assertConservation(myvar_orig_weighted, myvar_weighted, & - 'nonConservedMass_withFractionalAreas_specialShrinks', & + 'nonConservedMass_withVolumeMultipliers_specialShrinks', & begc = 1, endc = 2, non_conserved_mass = non_conserved_mass_grc(1)) - end subroutine nonConservedMass_withFractionalAreas_specialShrinks + end subroutine nonConservedMass_withVolumeMultipliers_specialShrinks @Test - subroutine nonConservedMass_withFractionalAreas_specialGrows(this) - ! This tests the computation of non_conserved_mass with fractional area terms, when + subroutine nonConservedMass_withVolumeMultipliers_specialGrows(this) + ! This tests the computation of non_conserved_mass with volume multiplier terms, when ! a special column grows class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: input_from_special = 11._r8 ! irrelevant - real(r8), parameter :: fractional_area_old(3) = [0.1_r8, 0.2_r8, 0.3_r8] - real(r8) :: fractional_area_new(3) + real(r8), parameter :: volume_multiplier_old(3) = [0.1_r8, 0.2_r8, 0.3_r8] + real(r8) :: volume_multiplier_new(3) real(r8) :: myvar(3), myvar_orig(3) real(r8) :: expected_non_conserved_mass real(r8) :: non_conserved_mass_grc(1) = 0._r8 @@ -1034,30 +1034,30 @@ contains myvar = [2._r8, 3._r8, 4._r8] myvar_orig = myvar expected_non_conserved_mass = -1._r8 * & - (myvar(1) * this%cwtgcell_new_minus_old(1) * fractional_area_old(1) + & - myvar(2) * this%cwtgcell_new_minus_old(2) * fractional_area_old(2)) + (myvar(1) * this%cwtgcell_new_minus_old(1) * volume_multiplier_old(1) + & + myvar(2) * this%cwtgcell_new_minus_old(2) * volume_multiplier_old(2)) - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_fill_special_using_fixed_value(bounds, clump_index=1, var=myvar, & special_value=input_from_special, non_conserved_mass_grc=non_conserved_mass_grc, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new) + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new) ! Verify @assertEqual(expected_non_conserved_mass, non_conserved_mass_grc(1), tolerance=tol) ! Note that the specially-treated column is excluded from the conservation check: the ! value in that column is ignored for conservation. - myvar_orig_weighted = myvar_orig * fractional_area_old - myvar_weighted = myvar * fractional_area_new + myvar_orig_weighted = myvar_orig * volume_multiplier_old + myvar_weighted = myvar * volume_multiplier_new call this%assertConservation(myvar_orig_weighted, myvar_weighted, & - 'nonConservedMass_withFractionalAreas_specialGrows', & + 'nonConservedMass_withVolumeMultipliers_specialGrows', & begc = 1, endc = 2, non_conserved_mass = non_conserved_mass_grc(1)) - end subroutine nonConservedMass_withFractionalAreas_specialGrows + end subroutine nonConservedMass_withVolumeMultipliers_specialGrows ! ------------------------------------------------------------------------ ! Tests of the optional adjustment output @@ -1133,13 +1133,13 @@ contains end subroutine adjustment_withSpecialGrowing_isUnchanged @Test - subroutine adjustment_withFractionalAreas(this) + subroutine adjustment_withVolumeMultipliers(this) class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater real(r8), parameter :: cwtlunit_old(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] real(r8), parameter :: cwtlunit_new(4) = [0.1_r8, 0.3_r8, 0.2_r8, 0.4_r8] - real(r8), parameter :: fractional_area_old(4) = [0.0_r8, 0.1_r8, 0.7_r8, 0.5_r8] - real(r8) :: myvar(4), myvar_orig(4), fractional_area_new(4) + real(r8), parameter :: volume_multiplier_old(4) = [0.0_r8, 0.1_r8, 0.7_r8, 0.5_r8] + real(r8) :: myvar(4), myvar_orig(4), volume_multiplier_new(4) real(r8) :: myvar_orig_weighted(4), myvar_weighted(4) real(r8) :: adjustment(4) @@ -1154,33 +1154,33 @@ contains myvar_orig = myvar adjustment(:) = 0._r8 - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=myvar, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) ! Verify - myvar_orig_weighted = myvar_orig*fractional_area_old - myvar_weighted = myvar*fractional_area_new + myvar_orig_weighted = myvar_orig*volume_multiplier_old + myvar_weighted = myvar*volume_multiplier_new @assertEqual(myvar_orig_weighted, myvar_weighted - adjustment) - end subroutine adjustment_withFractionalAreas + end subroutine adjustment_withVolumeMultipliers @Test - subroutine adjustment_withFractionalAreas_newFrac0(this) - ! If the new fractional area is 0, then the adjustment should be equal to - ! -(myvar_orig * fractional_area_old) + subroutine adjustment_withVolumeMultipliers_newMultiplier0(this) + ! If the new volume multiplier is 0, then the adjustment should be equal to + ! -(myvar_orig * volume_multiplier_old) class(TestColumnStateUpdater), intent(inout) :: this type(column_state_updater_type) :: cs_updater - ! col #2 grows from 0, taking over some of col #3; col #2's new fractional area - ! should be 0 because col #3 had a fractional area of 0 + ! col #2 grows from 0, taking over some of col #3; col #2's new volume multiplier + ! should be 0 because col #3 had a volume multiplier of 0 real(r8), parameter :: cwtlunit_old(4) = [0.1_r8, 0.0_r8, 0.5_r8, 0.4_r8] real(r8), parameter :: cwtlunit_new(4) = [0.1_r8, 0.2_r8, 0.3_r8, 0.4_r8] - real(r8), parameter :: fractional_area_old(4) = [0.0_r8, 0.1_r8, 0.0_r8, 0.5_r8] - real(r8) :: myvar(4), myvar_orig(4), fractional_area_new(4) + real(r8), parameter :: volume_multiplier_old(4) = [0.0_r8, 0.1_r8, 0.0_r8, 0.5_r8] + real(r8) :: myvar(4), myvar_orig(4), volume_multiplier_new(4) real(r8) :: myvar_orig_weighted(4), myvar_weighted(4) real(r8) :: adjustment(4) real(r8) :: expected_adjustment @@ -1196,18 +1196,18 @@ contains myvar_orig = myvar adjustment(:) = 0._r8 - fractional_area_new(:) = fractional_area_old(:) - call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=fractional_area_new) + volume_multiplier_new(:) = volume_multiplier_old(:) + call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=volume_multiplier_new) ! Exercise call cs_updater%update_column_state_no_special_handling(bounds, clump_index=1, var=myvar, & - fractional_area_old = fractional_area_old, & - fractional_area_new = fractional_area_new, & + volume_multiplier_old = volume_multiplier_old, & + volume_multiplier_new = volume_multiplier_new, & adjustment = adjustment) ! Verify - expected_adjustment = -1._r8 * (myvar_orig(2) * fractional_area_old(2)) + expected_adjustment = -1._r8 * (myvar_orig(2) * volume_multiplier_old(2)) @assertEqual(expected_adjustment, adjustment(2)) - end subroutine adjustment_withFractionalAreas_newFrac0 + end subroutine adjustment_withVolumeMultipliers_newMultiplier0 end module test_column_state_updater