diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index e9b882a14a..255f30f40a 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -7,6 +7,7 @@ FAIL SMS.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN FAIL SMS_D.f10_f10_musgs.I2000Clm50BgcCrop.hobart_pgi.clm-crop RUN FAIL ERS_D_Ln9_P480x3.f19_g16.I2000Clm50SpGs.cheyenne_intel.clm-waccmx_offline COMPARE_base_rest + FAIL PFS_Ld20.f09_g17.I2000Clm50BgcCrop.cheyenne_intel GENERATE FAIL ERS_Ld60.f45_f45_mg37.I2000Clm45Fates.hobart_nag.clm-Fates COMPARE_base_rest diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index fa14f1eb4f..8847e1f0bc 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1980,6 +1980,16 @@ + + + + + + + + + + diff --git a/doc/.ChangeLog_template b/doc/.ChangeLog_template index bb8bf7f789..93fb12d18f 100644 --- a/doc/.ChangeLog_template +++ b/doc/.ChangeLog_template @@ -46,7 +46,7 @@ Changes made to namelist defaults (e.g., changed parameter values): Changes to the datasets (e.g., parameter, surface or initial files): -Substantial timing or memory changes: +Substantial timing or memory changes: [For timing changes, can check PFS test(s) in the test suite] Notes of particular relevance for developers: (including Code reviews and testing) --------------------------------------------- diff --git a/doc/ChangeLog b/doc/ChangeLog index aa447b7b21..78aa2e3f05 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,147 @@ =============================================================== +Tag name: ctsm1.0.dev018 +Originator(s): sacks (Bill Sacks) +Date: Thu Nov 29 16:03:50 MST 2018 +One-line Summary: Water tracer updates for initial things in driver loop + +Purpose of changes +------------------ + +Update water tracers for initial stuff done in driver loop. This includes +atm2lnd forcings (non-downscaled and downscaled), balance check initialization, +and dyn subgrid updates. + +Broadly speaking, the changes here are: + +(1) Reworked WaterType to make it easier / more robust for other code to loop + over tracers or bulk+tracers + +(2) The most interesting changes are probably the code to update the atm2lnd + water tracers (in Wateratm2lndType.F90 and WaterTracerUtils.F90) + +(3) In various other places, do some infrastructurey stuff (initializing water + balance, doing dyn subgrid stuff) for tracers as well as bulk + +(4) Supporting unit tests and unit test infrastructure + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): +Resolves ESCOMP/ctsm#487 +Resolves ESCOMP/ctsm#488 +Resolves ESCOMP/ctsm#489 + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm5_0 + +[ ] clm4_5 + +[ ] clm4_0 + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): none + +Changes to CTSM's user interface (e.g., new/renamed XML or namelist variables): none + +Changes made to namelist defaults (e.g., changed parameter values): none + +Changes to the datasets (e.g., parameter, surface or initial files): none + +Substantial timing or memory changes: none + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- +NOTE: Be sure to review the steps in ../CTSMMasterChecklist as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): +- We have chosen not to set all water tracers as soon as possible, but instead + to do these tracer settings later in the driver loop. This keeps the driver + loop cleaner, but means that you cannot arbitrarily sprinkle calls to + TracerConsistencyCheck throughout the driver. Specifically for this tag: the + non-downscaled, gridcell-level atm2lnd water tracers are not updated until + after the call to downscale_forcings, so tracer consistency checks before that + point would fail. + +Changes to tests or testing: +- Added a PFS test + +Code reviewed by: Portions of the design (and possibly code) have been reviewed +by Mat Rothstein, David Noone and Mariana Vertenstein + + +CTSM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - not run + + tools-tests (test/tools): + + cheyenne - not run + + PTCLM testing (tools/shared/PTCLM/test): + + cheyenne - not run + + regular tests (aux_clm): + + cheyenne ---- ok + hobart ------ ok + + ok means tests pass, some answers change as expected + +CTSM tag used for the baseline comparisons: ctsm1.0.dev017 + + +Answer changes +-------------- + +Changes answers relative to baseline: YES + + If a tag changes answers relative to baseline comparison the + following should be filled in (otherwise remove this section): + + Summarize any changes to answers, i.e., + - what code configurations: many + - what platforms/compilers: all + - nature of change (roundoff; larger than roundoff/same climate; new climate): + roundoff-level changes in sensible heat flux from precip conversion due to + refactoring this calculation; everything else bit-for-bit + + If bitwise differences were observed, how did you show they were no worse + than roundoff? via summarize_cprnc_diffs to see differences in the test suite + + If this tag changes climate describe the run(s) done to evaluate the new + climate (put details of the simulations in the experiment database) + - casename: N/A + + URL for LMWG diagnostics output used to validate new climate: N/A + + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none + +Pull Requests that document the changes (include PR ids): +https://github.com/ESCOMP/ctsm/pull/572 + +=============================================================== +=============================================================== Tag name: ctsm1.0.dev017 Originator(s): slevis (Samuel Levis, Slevis Consulting LLC,303-665-1310) Date: Wed Nov 28 14:27:50 MST 2018 diff --git a/doc/ChangeSum b/doc/ChangeSum index 2e0a8eea8c..ebdd943719 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm1.0.dev018 sacks 11/29/2018 Water tracer updates for initial things in driver loop ctsm1.0.dev017 slevis 11/28/2018 Merge the collapse2gencrop branch ctsm1.0.dev016 sacks 11/01/2018 Update cime, fix FATES DEBUG token, add script to easily run system tests ctsm1.0.dev015 sacks 10/28/2018 CMIP6 compset modifiers, usermods for typical output, and other output enhancements diff --git a/doc/design/water_tracers.rst b/doc/design/water_tracers.rst index b265e72011..b52d653393 100644 --- a/doc/design/water_tracers.rst +++ b/doc/design/water_tracers.rst @@ -113,6 +113,48 @@ of understanding the suite of water flux variables, and reducing the number of i that need to get passed all around (e.g., we won't need to pass irrigation_inst to as many places). +======================== + Loops over all tracers +======================== + +Initially, I was hoping that we could keep loops over tracers in WaterType, for the +following reasons: + +1. To keep this complexity out of other modules + +2. To make it easier to change the details of how the bulk and tracer instances are + stored, if we ever need to: By keeping as many of the loops as possible in WaterType, + we reduce the number of places that would need to be changed + +However, it was starting to get too awkward to require all loops over tracers to happen in +WaterType (or some other centralized location): I had originally imagined that we wouldn't +need too many loops over tracers, but it turns out that we need loops over tracers in a +lot of places. Requiring all of these loops over tracers to be in WaterType both (a) +bloats that module, and (b) adds extra indirection (which makes it harder to understand +the code, because you're bouncing back and forth between more modules, and has possible +performance implications as we break routines into tiny pieces for this purpose). + +So we allow loops over tracers (or bulk plus tracers) anywhere in the code. See comments +at the top of WaterType.F90 for example code showing how to write these loops. + +Note that the bulk instances (``waterfluxbulk_inst``, etc.) can be obtained in two ways: + +1. Using ``water_inst%water*bulk_inst`` + +2. As one of the indices in ``water_inst%bulk_and_tracers(:)%water*_inst`` + +Method (2) is just meant to be used when you are doing the same operation on bulk water +and all water tracers. Reasons why it is better, or necessary, to use method (1) when you +are really just working with bulk water are: + +- This makes it more explicit and clear that you are working with bulk water + +- When passing an argument to a subroutine whose dummy argument is of type + ``water*bulk_type``, method (2) only works if you surround the call with a ``select + type`` statement, which is awkward, to say the least. (Subroutines that expect bulk + water should generally declare their dummy arguments to be of type ``water*bulk_type``, + as discussed in `Object-oriented design for water tracer types`_.) + ============================================== Infrastructure for looping through variables ============================================== diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index d7b7fcd0ae..f214251a81 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -17,13 +17,15 @@ module BalanceCheckMod use atm2lndType , only : atm2lnd_type use EnergyFluxType , only : energyflux_type use SolarAbsorbedType , only : solarabs_type - use SoilHydrologyType , only : soilhydrology_type - use WaterStateBulkType , only : waterstatebulk_type - use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type - use Wateratm2lndBulkType , only : wateratm2lndbulk_type - use WaterBalanceType , only : waterbalance_type - use SurfaceAlbedoType , only : surfalb_type - use WaterFluxBulkType , only : waterfluxbulk_type + use SoilHydrologyType , only : soilhydrology_type + use SurfaceAlbedoType , only : surfalb_type + use WaterStateType , only : waterstate_type + use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use WaterDiagnosticType, only : waterdiagnostic_type + use Wateratm2lndType , only : wateratm2lnd_type + use WaterBalanceType , only : waterbalance_type + use WaterFluxType , only : waterflux_type + use WaterType , only : water_type use TotalWaterAndHeatMod, only : ComputeWaterMassNonLake, ComputeWaterMassLake use GridcellType , only : grc use LandunitType , only : lun @@ -51,6 +53,10 @@ module BalanceCheckMod real(r8), private, parameter :: skip_size = 3600.0_r8 ! Time steps to skip the balance check at startup (sec) integer, private :: skip_steps = -999 ! Number of time steps to skip the balance check for + ! + ! !PRIVATE MEMBER FUNCTIONS: + private :: BeginWaterBalanceSingle ! Initialize water balance check for bulk or a single tracer + character(len=*), parameter, private :: sourcefile = & __FILE__ !----------------------------------------------------------------------- @@ -117,10 +123,47 @@ end function GetBalanceCheckSkipSteps !----------------------------------------------------------------------- subroutine BeginWaterBalance(bounds, & num_nolakec, filter_nolakec, num_lakec, filter_lakec, & - soilhydrology_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst) + water_inst, soilhydrology_inst) + ! + ! !DESCRIPTION: + ! Initialize column-level water balance at beginning of time step, for bulk water and + ! each water tracer + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_nolakec ! number of column non-lake points in column filter + integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points + integer , intent(in) :: num_lakec ! number of column lake points in column filter + integer , intent(in) :: filter_lakec(:) ! column filter for lake points + type(water_type) , intent(inout) :: water_inst + type(soilhydrology_type) , intent(in) :: soilhydrology_inst + ! + ! !LOCAL VARIABLES: + integer :: i + + character(len=*), parameter :: subname = 'BeginWaterBalance' + !----------------------------------------------------------------------- + + do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end + call BeginWaterBalanceSingle(bounds, & + num_nolakec, filter_nolakec, & + num_lakec, filter_lakec, & + soilhydrology_inst, & + water_inst%bulk_and_tracers(i)%waterstate_inst, & + water_inst%bulk_and_tracers(i)%waterdiagnostic_inst, & + water_inst%bulk_and_tracers(i)%waterbalance_inst) + end do + + end subroutine BeginWaterBalance + + !----------------------------------------------------------------------- + subroutine BeginWaterBalanceSingle(bounds, & + num_nolakec, filter_nolakec, num_lakec, filter_lakec, & + soilhydrology_inst, waterstate_inst, waterdiagnostic_inst, waterbalance_inst) ! ! !DESCRIPTION: - ! Initialize column-level water balance at beginning of time step + ! Initialize column-level water balance at beginning of time step, for bulk or a + ! single tracer ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds @@ -128,20 +171,22 @@ subroutine BeginWaterBalance(bounds, & integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points - type(soilhydrology_type) , intent(inout) :: soilhydrology_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst + type(soilhydrology_type) , intent(in) :: soilhydrology_inst + class(waterstate_type) , intent(inout) :: waterstate_inst + class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst + class(waterbalance_type) , intent(inout) :: waterbalance_inst ! ! !LOCAL VARIABLES: integer :: c, j, fc ! indices !----------------------------------------------------------------------- - associate( & - zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) - zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) - wa => waterstatebulk_inst%wa_col , & ! Output: [real(r8) (:) ] water in the unconfined aquifer (mm) - begwb => waterbalancebulk_inst%begwb_col & ! Output: [real(r8) (:) ] water mass begining of the time step + associate( & + zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) + zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) + h2osno => waterstate_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) + wa => waterstate_inst%wa_col , & ! Output: [real(r8) (:) ] water in the unconfined aquifer (mm) + begwb => waterbalance_inst%begwb_col , & ! Output: [real(r8) (:) ] water mass begining of the time step + h2osno_old => waterbalance_inst%h2osno_old_col & ! Output: [real(r8) (:) ] snow water (mm H2O) at previous time step ) do fc = 1, num_nolakec @@ -154,19 +199,23 @@ subroutine BeginWaterBalance(bounds, & end do call ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, begwb(bounds%begc:bounds%endc)) + waterstate_inst, waterdiagnostic_inst, begwb(bounds%begc:bounds%endc)) call ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstatebulk_inst, begwb(bounds%begc:bounds%endc)) + waterstate_inst, begwb(bounds%begc:bounds%endc)) + + do c = bounds%begc, bounds%endc + h2osno_old(c) = h2osno(c) + end do end associate - end subroutine BeginWaterBalance + end subroutine BeginWaterBalanceSingle !----------------------------------------------------------------------- subroutine BalanceCheck( bounds, & - atm2lnd_inst, solarabs_inst, waterfluxbulk_inst, waterstatebulk_inst, & - waterdiagnosticbulk_inst, waterbalancebulk_inst, wateratm2lndbulk_inst, & + atm2lnd_inst, solarabs_inst, waterflux_inst, waterstate_inst, & + waterdiagnosticbulk_inst, waterbalance_inst, wateratm2lnd_inst, & surfalb_inst, energyflux_inst, canopystate_inst) ! ! !DESCRIPTION: @@ -194,12 +243,12 @@ subroutine BalanceCheck( bounds, & type(bounds_type) , intent(in) :: bounds type(atm2lnd_type) , intent(in) :: atm2lnd_inst type(solarabs_type) , intent(in) :: solarabs_inst - type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst - type(wateratm2lndbulk_type) , intent(inout) :: wateratm2lndbulk_inst - type(surfalb_type) , intent(inout) :: surfalb_inst + class(waterflux_type) , intent(in) :: waterflux_inst + class(waterstate_type), intent(in) :: waterstate_inst + type(waterdiagnosticbulk_type), intent(in) :: waterdiagnosticbulk_inst + class(waterbalance_type), intent(inout) :: waterbalance_inst + class(wateratm2lnd_type) , intent(in) :: wateratm2lnd_inst + type(surfalb_type) , intent(in) :: surfalb_inst type(energyflux_type) , intent(inout) :: energyflux_inst type(canopystate_type), intent(inout) :: canopystate_inst ! @@ -215,59 +264,50 @@ subroutine BalanceCheck( bounds, & !----------------------------------------------------------------------- associate( & - volr => wateratm2lndbulk_inst%volr_grc , & ! Input: [real(r8) (:) ] river water storage (m3) forc_solad => atm2lnd_inst%forc_solad_grc , & ! Input: [real(r8) (:,:) ] direct beam radiation (vis=forc_sols , nir=forc_soll ) forc_solai => atm2lnd_inst%forc_solai_grc , & ! Input: [real(r8) (:,:) ] diffuse radiation (vis=forc_solsd, nir=forc_solld) - forc_rain => wateratm2lndbulk_inst%forc_rain_downscaled_col , & ! Input: [real(r8) (:) ] rain rate [mm/s] - forc_snow => wateratm2lndbulk_inst%forc_snow_downscaled_col , & ! Input: [real(r8) (:) ] snow rate [mm/s] + forc_rain => wateratm2lnd_inst%forc_rain_downscaled_col , & ! Input: [real(r8) (:) ] rain rate [mm/s] + forc_snow => wateratm2lnd_inst%forc_snow_downscaled_col , & ! Input: [real(r8) (:) ] snow rate [mm/s] forc_lwrad => atm2lnd_inst%forc_lwrad_downscaled_col , & ! Input: [real(r8) (:) ] downward infrared (longwave) radiation (W/m**2) - h2osno => waterstatebulk_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) - h2osno_old => waterbalancebulk_inst%h2osno_old_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) at previous time step + h2osno => waterstate_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) + h2osno_old => waterbalance_inst%h2osno_old_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) at previous time step frac_sno_eff => waterdiagnosticbulk_inst%frac_sno_eff_col , & ! Input: [real(r8) (:) ] effective snow fraction frac_sno => waterdiagnosticbulk_inst%frac_sno_col , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) snow_depth => waterdiagnosticbulk_inst%snow_depth_col , & ! Input: [real(r8) (:) ] snow height (m) - begwb => waterbalancebulk_inst%begwb_col , & ! Input: [real(r8) (:) ] water mass begining of the time step - errh2o => waterbalancebulk_inst%errh2o_col , & ! Output: [real(r8) (:) ] water conservation error (mm H2O) - errh2osno => waterbalancebulk_inst%errh2osno_col , & ! Output: [real(r8) (:) ] error in h2osno (kg m-2) - endwb => waterbalancebulk_inst%endwb_col , & ! Output: [real(r8) (:) ] water mass end of the time step - total_plant_stored_h2o_col => waterdiagnosticbulk_inst%total_plant_stored_h2o_col, & ! Input: [real(r8) (:) ] water mass in plant tissues (kg m-2) - qflx_rootsoi_col => waterfluxbulk_inst%qflx_rootsoi_col , & ! Input [real(r8) (:) ] water loss in soil layers to root uptake (mm H2O/s) - ! (ie transpiration demand, often = transpiration) - qflx_rain_grnd_col => waterfluxbulk_inst%qflx_rain_grnd_col , & ! Input: [real(r8) (:) ] rain on ground after interception (mm H2O/s) [+] - qflx_snow_grnd_col => waterfluxbulk_inst%qflx_snow_grnd_col , & ! Input: [real(r8) (:) ] snow on ground after interception (mm H2O/s) [+] - qflx_evap_soi => waterfluxbulk_inst%qflx_evap_soi_col , & ! Input: [real(r8) (:) ] soil evaporation (mm H2O/s) (+ = to atm) - qflx_snwcp_liq => waterfluxbulk_inst%qflx_snwcp_liq_col , & ! Input: [real(r8) (:) ] excess liquid h2o due to snow capping (outgoing) (mm H2O /s) [+]` - qflx_snwcp_ice => waterfluxbulk_inst%qflx_snwcp_ice_col , & ! Input: [real(r8) (:) ] excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]` - qflx_snwcp_discarded_liq => waterfluxbulk_inst%qflx_snwcp_discarded_liq_col, & ! Input: [real(r8) (:) ] excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]` - qflx_snwcp_discarded_ice => waterfluxbulk_inst%qflx_snwcp_discarded_ice_col, & ! Input: [real(r8) (:) ] excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]` - qflx_evap_tot => waterfluxbulk_inst%qflx_evap_tot_col , & ! Input: [real(r8) (:) ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg - qflx_dew_snow => waterfluxbulk_inst%qflx_dew_snow_col , & ! Input: [real(r8) (:) ] surface dew added to snow pack (mm H2O /s) [+] - qflx_sub_snow => waterfluxbulk_inst%qflx_sub_snow_col , & ! Input: [real(r8) (:) ] sublimation rate from snow pack (mm H2O /s) [+] - qflx_evap_grnd => waterfluxbulk_inst%qflx_evap_grnd_col , & ! Input: [real(r8) (:) ] ground surface evaporation rate (mm H2O/s) [+] - qflx_dew_grnd => waterfluxbulk_inst%qflx_dew_grnd_col , & ! Input: [real(r8) (:) ] ground surface dew formation (mm H2O /s) [+] - qflx_prec_grnd => waterfluxbulk_inst%qflx_prec_grnd_col , & ! Input: [real(r8) (:) ] water onto ground including canopy runoff [kg/(m2 s)] - qflx_snow_h2osfc => waterfluxbulk_inst%qflx_snow_h2osfc_col , & ! Input: [real(r8) (:) ] snow falling on surface water (mm/s) - qflx_h2osfc_to_ice => waterfluxbulk_inst%qflx_h2osfc_to_ice_col , & ! Input: [real(r8) (:) ] conversion of h2osfc to ice - qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s) - qflx_floodc => waterfluxbulk_inst%qflx_floodc_col , & ! Input: [real(r8) (:) ] total runoff due to flooding - qflx_snow_drain => waterfluxbulk_inst%qflx_snow_drain_col , & ! Input: [real(r8) (:) ] drainage from snow pack - qflx_surf => waterfluxbulk_inst%qflx_surf_col , & ! Input: [real(r8) (:) ] surface runoff (mm H2O /s) - qflx_qrgwl => waterfluxbulk_inst%qflx_qrgwl_col , & ! Input: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes - qflx_drain => waterfluxbulk_inst%qflx_drain_col , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s) - qflx_runoff => waterfluxbulk_inst%qflx_runoff_col , & ! Input: [real(r8) (:) ] total runoff (mm H2O /s) - qflx_ice_runoff_snwcp => waterfluxbulk_inst%qflx_ice_runoff_snwcp_col, & ! Input: [real(r8) (:) ] solid runoff from snow capping (mm H2O /s) - qflx_ice_runoff_xs => waterfluxbulk_inst%qflx_ice_runoff_xs_col , & ! Input: [real(r8) (:) ] solid runoff from excess ice in soil (mm H2O /s) - qflx_top_soil => waterfluxbulk_inst%qflx_top_soil_col , & ! Input: [real(r8) (:) ] net water input into soil from top (mm/s) - qflx_sl_top_soil => waterfluxbulk_inst%qflx_sl_top_soil_col , & ! Input: [real(r8) (:) ] liquid water + ice from layer above soil to top soil layer or sent to qflx_qrgwl (mm H2O/s) - qflx_liq_dynbal => waterfluxbulk_inst%qflx_liq_dynbal_grc , & ! Input: [real(r8) (:) ] liq runoff due to dynamic land cover change (mm H2O /s) - qflx_ice_dynbal => waterfluxbulk_inst%qflx_ice_dynbal_grc , & ! Input: [real(r8) (:) ] ice runoff due to dynamic land cover change (mm H2O /s) - snow_sources => waterfluxbulk_inst%snow_sources_col , & ! Output: [real(r8) (:) ] snow sources (mm H2O /s) - snow_sinks => waterfluxbulk_inst%snow_sinks_col , & ! Output: [real(r8) (:) ] snow sinks (mm H2O /s) - - qflx_irrig => waterfluxbulk_inst%qflx_irrig_col , & ! Input: [real(r8) (:) ] irrigation flux (mm H2O /s) - - qflx_glcice_dyn_water_flux => waterfluxbulk_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system) + begwb => waterbalance_inst%begwb_col , & ! Input: [real(r8) (:) ] water mass begining of the time step + errh2o => waterbalance_inst%errh2o_col , & ! Output: [real(r8) (:) ] water conservation error (mm H2O) + errh2osno => waterbalance_inst%errh2osno_col , & ! Output: [real(r8) (:) ] error in h2osno (kg m-2) + endwb => waterbalance_inst%endwb_col , & ! Output: [real(r8) (:) ] water mass end of the time step + snow_sources => waterbalance_inst%snow_sources_col , & ! Output: [real(r8) (:) ] snow sources (mm H2O /s) + snow_sinks => waterbalance_inst%snow_sinks_col , & ! Output: [real(r8) (:) ] snow sinks (mm H2O /s) + qflx_rain_grnd_col => waterflux_inst%qflx_rain_grnd_col , & ! Input: [real(r8) (:) ] rain on ground after interception (mm H2O/s) [+] + qflx_snow_grnd_col => waterflux_inst%qflx_snow_grnd_col , & ! Input: [real(r8) (:) ] snow on ground after interception (mm H2O/s) [+] + qflx_snwcp_liq => waterflux_inst%qflx_snwcp_liq_col , & ! Input: [real(r8) (:) ] excess liquid h2o due to snow capping (outgoing) (mm H2O /s) [+]` + qflx_snwcp_ice => waterflux_inst%qflx_snwcp_ice_col , & ! Input: [real(r8) (:) ] excess solid h2o due to snow capping (outgoing) (mm H2O /s) [+]` + qflx_snwcp_discarded_liq => waterflux_inst%qflx_snwcp_discarded_liq_col, & ! Input: [real(r8) (:) ] excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]` + qflx_snwcp_discarded_ice => waterflux_inst%qflx_snwcp_discarded_ice_col, & ! Input: [real(r8) (:) ] excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) [+]` + qflx_evap_tot => waterflux_inst%qflx_evap_tot_col , & ! Input: [real(r8) (:) ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg + qflx_dew_snow => waterflux_inst%qflx_dew_snow_col , & ! Input: [real(r8) (:) ] surface dew added to snow pack (mm H2O /s) [+] + qflx_sub_snow => waterflux_inst%qflx_sub_snow_col , & ! Input: [real(r8) (:) ] sublimation rate from snow pack (mm H2O /s) [+] + qflx_evap_grnd => waterflux_inst%qflx_evap_grnd_col , & ! Input: [real(r8) (:) ] ground surface evaporation rate (mm H2O/s) [+] + qflx_dew_grnd => waterflux_inst%qflx_dew_grnd_col , & ! Input: [real(r8) (:) ] ground surface dew formation (mm H2O /s) [+] + qflx_prec_grnd => waterflux_inst%qflx_prec_grnd_col , & ! Input: [real(r8) (:) ] water onto ground including canopy runoff [kg/(m2 s)] + qflx_snow_h2osfc => waterflux_inst%qflx_snow_h2osfc_col , & ! Input: [real(r8) (:) ] snow falling on surface water (mm/s) + qflx_h2osfc_to_ice => waterflux_inst%qflx_h2osfc_to_ice_col , & ! Input: [real(r8) (:) ] conversion of h2osfc to ice + qflx_drain_perched => waterflux_inst%qflx_drain_perched_col , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s) + qflx_floodc => waterflux_inst%qflx_floodc_col , & ! Input: [real(r8) (:) ] total runoff due to flooding + qflx_snow_drain => waterflux_inst%qflx_snow_drain_col , & ! Input: [real(r8) (:) ] drainage from snow pack + qflx_surf => waterflux_inst%qflx_surf_col , & ! Input: [real(r8) (:) ] surface runoff (mm H2O /s) + qflx_qrgwl => waterflux_inst%qflx_qrgwl_col , & ! Input: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes + qflx_drain => waterflux_inst%qflx_drain_col , & ! Input: [real(r8) (:) ] sub-surface runoff (mm H2O /s) + qflx_ice_runoff_snwcp => waterflux_inst%qflx_ice_runoff_snwcp_col, & ! Input: [real(r8) (:) ] solid runoff from snow capping (mm H2O /s) + qflx_ice_runoff_xs => waterflux_inst%qflx_ice_runoff_xs_col , & ! Input: [real(r8) (:) ] solid runoff from excess ice in soil (mm H2O /s) + qflx_sl_top_soil => waterflux_inst%qflx_sl_top_soil_col , & ! Input: [real(r8) (:) ] liquid water + ice from layer above soil to top soil layer or sent to qflx_qrgwl (mm H2O/s) + + qflx_irrig => waterflux_inst%qflx_irrig_col , & ! Input: [real(r8) (:) ] irrigation flux (mm H2O /s) + + qflx_glcice_dyn_water_flux => waterflux_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system) eflx_lwrad_out => energyflux_inst%eflx_lwrad_out_patch , & ! Input: [real(r8) (:) ] emitted infrared (longwave) radiation (W/m**2) eflx_lwrad_net => energyflux_inst%eflx_lwrad_net_patch , & ! Input: [real(r8) (:) ] net infrared (longwave) rad (W/m**2) [+ = to atm] @@ -296,8 +336,8 @@ subroutine BalanceCheck( bounds, & fabd => surfalb_inst%fabd_patch , & ! Input: [real(r8) (:,:)] flux absorbed by canopy per unit direct flux fabi => surfalb_inst%fabi_patch , & ! Input: [real(r8) (:,:)] flux absorbed by canopy per unit indirect flux - albd => surfalb_inst%albd_patch , & ! Output: [real(r8) (:,:)] surface albedo (direct) - albi => surfalb_inst%albi_patch , & ! Output: [real(r8) (:,:)] surface albedo (diffuse) + albd => surfalb_inst%albd_patch , & ! Input: [real(r8) (:,:)] surface albedo (direct) + albi => surfalb_inst%albi_patch , & ! Input: [real(r8) (:,:)] surface albedo (diffuse) ftdd => surfalb_inst%ftdd_patch , & ! Input: [real(r8) (:,:)] down direct flux below canopy per unit direct flux ftid => surfalb_inst%ftid_patch , & ! Input: [real(r8) (:,:)] down diffuse flux below canopy per unit direct flux ftii => surfalb_inst%ftii_patch , & ! Input: [real(r8) (:,:)] down diffuse flux below canopy per unit diffuse flux @@ -395,8 +435,6 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'qflx_ice_runoff_xs = ',qflx_ice_runoff_xs(indexc)*dtime write(iulog,*)'qflx_snwcp_discarded_ice = ',qflx_snwcp_discarded_ice(indexc)*dtime write(iulog,*)'qflx_snwcp_discarded_liq = ',qflx_snwcp_discarded_liq(indexc)*dtime - write(iulog,*)'qflx_rootsoi_col(1:nlevsoil) = ',qflx_rootsoi_col(indexc,:)*dtime - write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc) write(iulog,*)'deltawb = ',endwb(indexc)-begwb(indexc) write(iulog,*)'deltawb/dtime = ',(endwb(indexc)-begwb(indexc))/dtime write(iulog,*)'deltaflux = ',forc_rain_col(indexc)+forc_snow_col(indexc) - (qflx_evap_tot(indexc) + & @@ -412,7 +450,6 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'errh2o = ',errh2o(indexc) write(iulog,*)'forc_rain = ',forc_rain_col(indexc)*dtime write(iulog,*)'forc_snow = ',forc_snow_col(indexc)*dtime - write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc) write(iulog,*)'endwb = ',endwb(indexc) write(iulog,*)'begwb = ',begwb(indexc) @@ -428,7 +465,6 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'qflx_glcice_dyn_water_flux = ', qflx_glcice_dyn_water_flux(indexc)*dtime write(iulog,*)'qflx_snwcp_discarded_ice = ',qflx_snwcp_discarded_ice(indexc)*dtime write(iulog,*)'qflx_snwcp_discarded_liq = ',qflx_snwcp_discarded_liq(indexc)*dtime - write(iulog,*)'qflx_rootsoi_col(1:nlevsoil) = ',qflx_rootsoi_col(indexc,:)*dtime write(iulog,*)'clm model is stopping' call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__)) end if @@ -520,7 +556,6 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'qflx_sub_snow = ',qflx_sub_snow(indexc)*dtime write(iulog,*)'qflx_snow_drain = ',qflx_snow_drain(indexc)*dtime write(iulog,*)'qflx_evap_grnd = ',qflx_evap_grnd(indexc)*dtime - write(iulog,*)'qflx_top_soil = ',qflx_top_soil(indexc)*dtime write(iulog,*)'qflx_dew_snow = ',qflx_dew_snow(indexc)*dtime write(iulog,*)'qflx_dew_grnd = ',qflx_dew_grnd(indexc)*dtime write(iulog,*)'qflx_snwcp_ice = ',qflx_snwcp_ice(indexc)*dtime diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index 03da5c3da3..afc976a223 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -15,8 +15,10 @@ module TotalWaterAndHeatMod use LandunitType , only : lun use subgridAveMod , only : p2c use SoilHydrologyType , only : soilhydrology_type - use WaterStateBulkType , only : waterstatebulk_type - use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type + use WaterStateBulkType , only : waterstatebulk_type + use WaterStateType , only : waterstate_type + use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type + use WaterDiagnosticType, only : waterdiagnostic_type use UrbanParamsType , only : urbanparams_type use SoilStateType , only : soilstate_type use TemperatureType , only : temperature_type @@ -84,17 +86,20 @@ module TotalWaterAndHeatMod !----------------------------------------------------------------------- subroutine ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, water_mass) + waterstate_inst, waterdiagnostic_inst, water_mass) ! ! !DESCRIPTION: ! Compute total water mass for all non-lake columns ! + ! This can also be used to compute the mass of a specific water tracer (rather than + ! bulk water). + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_nolakec ! number of column non-lake points in column filter integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points - type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst + class(waterstate_type) , intent(in) :: waterstate_inst + class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2) ! ! !LOCAL VARIABLES: @@ -111,8 +116,8 @@ subroutine ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, & bounds = bounds, & num_nolakec = num_nolakec, & filter_nolakec = filter_nolakec, & - waterstatebulk_inst = waterstatebulk_inst, & - waterdiagnosticbulk_inst = waterdiagnosticbulk_inst, & + waterstate_inst = waterstate_inst, & + waterdiagnostic_inst = waterdiagnostic_inst, & liquid_mass = liquid_mass(bounds%begc:bounds%endc), & ice_mass = ice_mass(bounds%begc:bounds%endc)) @@ -125,16 +130,19 @@ end subroutine ComputeWaterMassNonLake !----------------------------------------------------------------------- subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & - waterstatebulk_inst, water_mass) + waterstate_inst, water_mass) ! ! !DESCRIPTION: ! Compute total water mass for all lake columns ! + ! This can also be used to compute the mass of a specific water tracer (rather than + ! bulk water). + ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points - type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst + class(waterstate_type) , intent(in) :: waterstate_inst real(r8) , intent(inout) :: water_mass( bounds%begc: ) ! computed water mass (kg m-2) ! ! !LOCAL VARIABLES: @@ -151,7 +159,7 @@ subroutine ComputeWaterMassLake(bounds, num_lakec, filter_lakec, & bounds = bounds, & num_lakec = num_lakec, & filter_lakec = filter_lakec, & - waterstatebulk_inst = waterstatebulk_inst, & + waterstate_inst = waterstate_inst, & liquid_mass = liquid_mass(bounds%begc:bounds%endc), & ice_mass = ice_mass(bounds%begc:bounds%endc)) @@ -165,11 +173,14 @@ end subroutine ComputeWaterMassLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, liquid_mass, ice_mass) + waterstate_inst, waterdiagnostic_inst, liquid_mass, ice_mass) ! ! !DESCRIPTION: ! Compute total water mass for all non-lake columns, separated into liquid and ice ! + ! This can also be used to compute the mass of a specific water tracer (rather than + ! bulk water). + ! ! Note: Changes to this routine should generally be accompanied by similar changes ! to ComputeHeatNonLake ! @@ -177,8 +188,8 @@ subroutine ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_nolakec ! number of column non-lake points in column filter integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points - type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst + class(waterstate_type) , intent(in) :: waterstate_inst + class(waterdiagnostic_type), intent(in) :: waterdiagnostic_inst real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2) real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! computed ice mass (kg m-2) ! @@ -198,15 +209,14 @@ subroutine ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & associate( & snl => col%snl , & ! Input: [integer (:) ] negative number of snow layers - h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) - h2osno => waterstatebulk_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) - h2ocan_patch => waterstatebulk_inst%h2ocan_patch , & ! Input: [real(r8) (:) ] canopy water (mm H2O) - snocan_patch => waterstatebulk_inst%snocan_patch , & ! Input: [real(r8) (:) ] canopy snow water (mm H2O) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) - h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) - total_plant_stored_h2o => waterdiagnosticbulk_inst%total_plant_stored_h2o_col, & - ! Input: [real(r8) (:,:) ] plant internal stored water (mm H2O) - wa => waterstatebulk_inst%wa_col & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm) + h2osfc => waterstate_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) + h2osno => waterstate_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) + h2ocan_patch => waterstate_inst%h2ocan_patch , & ! Input: [real(r8) (:) ] canopy water (mm H2O) + snocan_patch => waterstate_inst%snocan_patch , & ! Input: [real(r8) (:) ] canopy snow water (mm H2O) + h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) + h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + total_plant_stored_h2o => waterdiagnostic_inst%total_plant_stored_h2o_col, & ! Input: [real(r8) (:,:) ] plant internal stored water (mm H2O) + wa => waterstate_inst%wa_col & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm) ) do fc = 1, num_nolakec @@ -226,7 +236,7 @@ subroutine ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & do fc = 1, num_nolakec c = filter_nolakec(fc) - ! waterstatebulk_inst%snocan_patch and waterstatebulk_inst%liqcan_patch are only set if + ! waterstate_inst%snocan_patch and waterstate_inst%liqcan_patch are only set if ! we're using snow-on-veg; otherwise they are 0. However, we can rely on ! h2ocan_patch being set in all cases, so we can always determine the liquid mass ! as (h2ocan - snocan). @@ -303,11 +313,14 @@ end subroutine ComputeLiqIceMassNonLake !----------------------------------------------------------------------- subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & - waterstatebulk_inst, liquid_mass, ice_mass) + waterstate_inst, liquid_mass, ice_mass) ! ! !DESCRIPTION: ! Compute total water mass for all lake columns, separated into liquid and ice ! + ! This can also be used to compute the mass of a specific water tracer (rather than + ! bulk water). + ! ! Note: Changes to this routine should generally be accompanied by similar changes ! to ComputeHeatLake ! @@ -315,7 +328,7 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_lakec ! number of column lake points in column filter integer , intent(in) :: filter_lakec(:) ! column filter for lake points - type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst + class(waterstate_type), intent(in) :: waterstate_inst real(r8) , intent(inout) :: liquid_mass( bounds%begc: ) ! computed liquid water mass (kg m-2) real(r8) , intent(inout) :: ice_mass( bounds%begc: ) ! computed ice mass (kg m-2) ! @@ -331,9 +344,9 @@ subroutine ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & associate( & snl => col%snl , & ! Input: [integer (:) ] negative number of snow layers - h2osno => waterstatebulk_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) - h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) - h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + h2osno => waterstate_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) + h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) + h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) ) do fc = 1, num_lakec diff --git a/src/biogeophys/WaterBalanceType.F90 b/src/biogeophys/WaterBalanceType.F90 index 7cf0ca5a06..e62ddc4eb9 100644 --- a/src/biogeophys/WaterBalanceType.F90 +++ b/src/biogeophys/WaterBalanceType.F90 @@ -33,6 +33,8 @@ module WaterBalanceType real(r8), pointer :: ice1_grc (:) ! grc initial gridcell total h2o ice content real(r8), pointer :: ice2_grc (:) ! grc post land cover change total ice content + real(r8), pointer :: snow_sources_col (:) ! col snow sources (mm H2O/s) + real(r8), pointer :: snow_sinks_col (:) ! col snow sinks (mm H2O/s) ! Balance Checks @@ -105,6 +107,14 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%ice2_grc, name = 'ice2_grc', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_GRIDCELL) + + call AllocateVar1d(var = this%snow_sources_col, name = 'snow_sources_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%snow_sinks_col, name = 'snow_sinks_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%begwb_col, name = 'begwb_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) @@ -149,6 +159,28 @@ subroutine InitHistory(this, bounds) begc = bounds%begc; endc= bounds%endc begg = bounds%begg; endg= bounds%endg + ! As defined here, snow_sources - snow_sinks will equal the change in h2osno at any + ! given time step but only if there is at least one snow layer (for all landunits + ! except lakes). Also note that monthly average files of snow_sources and snow_sinks + ! sinks must be weighted by number of days in the month to diagnose, for example, an + ! annual value of the change in h2osno. + + this%snow_sources_col(begc:endc) = spval + call hist_addfld1d ( & + fname=this%info%fname('SNOW_SOURCES'), & + units='mm/s', & + avgflag='A', & + long_name=this%info%lname('snow sources (liquid water)'), & + ptr_col=this%snow_sources_col, c2l_scale_type='urbanf') + + this%snow_sinks_col(begc:endc) = spval + call hist_addfld1d ( & + fname=this%info%fname('SNOW_SINKS'), & + units='mm/s', & + avgflag='A', & + long_name=this%info%lname('snow sinks (liquid water)'), & + ptr_col=this%snow_sinks_col, c2l_scale_type='urbanf') + this%liq1_grc(begg:endg) = spval call hist_addfld1d ( & fname=this%info%fname('LIQUID_CONTENT1'), & diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 1fea963eda..55bcc6984c 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -45,12 +45,6 @@ module WaterDiagnosticBulkType real(r8), pointer :: snounload_patch (:) ! Canopy snow unloading (mm H2O) real(r8), pointer :: swe_old_col (:,:) ! col initial snow water - real(r8), pointer :: total_plant_stored_h2o_col(:) ! col water that is bound in plants, including roots, sapwood, leaves, etc - ! in most cases, the vegetation scheme does not have a dynamic - ! water storage in plants, and thus 0.0 is a suitable for the trivial case. - ! When FATES is coupled in with plant hydraulics turned on, this storage - ! term is set to non-zero. (kg/m2 H2O) - real(r8), pointer :: snw_rds_col (:,:) ! col snow grain radius (col,lyr) [m^-6, microns] real(r8), pointer :: snw_rds_top_col (:) ! col snow grain radius (top layer) [m^-6, microns] real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] @@ -155,8 +149,6 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%snounload_patch (begp:endp)) ; this%snounload_patch (:) = nan allocate(this%swe_old_col (begc:endc,-nlevsno+1:0)) ; this%swe_old_col (:,:) = nan - allocate(this%total_plant_stored_h2o_col(begc:endc)) ; this%total_plant_stored_h2o_col(:) = nan - allocate(this%snw_rds_col (begc:endc,-nlevsno+1:0)) ; this%snw_rds_col (:,:) = nan allocate(this%snw_rds_top_col (begc:endc)) ; this%snw_rds_top_col (:) = nan allocate(this%h2osno_top_col (begc:endc)) ; this%h2osno_top_col (:) = nan @@ -175,7 +167,7 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%frac_h2osfc_col (begc:endc)) ; this%frac_h2osfc_col (:) = nan allocate(this%frac_h2osfc_nosnow_col (begc:endc)) ; this%frac_h2osfc_nosnow_col (:) = nan allocate(this%wf_col (begc:endc)) ; this%wf_col (:) = nan - allocate(this%wf2_col (begc:endc)) ; + allocate(this%wf2_col (begc:endc)) ; this%wf2_col (:) = nan allocate(this%fwet_patch (begp:endp)) ; this%fwet_patch (:) = nan allocate(this%fcansno_patch (begp:endp)) ; this%fcansno_patch (:) = nan allocate(this%fdry_patch (begp:endp)) ; this%fdry_patch (:) = nan @@ -502,11 +494,6 @@ subroutine InitBulkCold(this, bounds, & end do - ! Water Stored in plants is almost always a static entity, with the exception - ! of when FATES-hydraulics is used. As such, this is trivially set to 0.0 (rgk 03-2017) - this%total_plant_stored_h2o_col(bounds%begc:bounds%endc) = 0.0_r8 - - associate(snl => col%snl) this%snounload_patch(bounds%begp:bounds%endp) = 0._r8 diff --git a/src/biogeophys/WaterDiagnosticType.F90 b/src/biogeophys/WaterDiagnosticType.F90 index 57c8dceca6..c6a01d0445 100644 --- a/src/biogeophys/WaterDiagnosticType.F90 +++ b/src/biogeophys/WaterDiagnosticType.F90 @@ -33,6 +33,12 @@ module WaterDiagnosticType real(r8), pointer :: snowice_col (:) ! col average snow ice lens real(r8), pointer :: snowliq_col (:) ! col average snow liquid water + real(r8), pointer :: total_plant_stored_h2o_col(:) ! col water that is bound in plants, including roots, sapwood, leaves, etc + ! in most cases, the vegetation scheme does not have a dynamic + ! water storage in plants, and thus 0.0 is a suitable for the trivial case. + ! When FATES is coupled in with plant hydraulics turned on, this storage + ! term is set to non-zero. (kg/m2 H2O) + real(r8), pointer :: h2osoi_liqice_10cm_col (:) ! col liquid water + ice lens in top 10cm of soil (kg/m2) real(r8), pointer :: tws_grc (:) ! grc total water storage (mm H2O) real(r8), pointer :: q_ref2m_patch (:) ! patch 2 m height surface specific humidity (kg/kg) @@ -100,6 +106,9 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%snowliq_col, name = 'snowliq_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%total_plant_stored_h2o_col, name = 'total_plant_stored_h2o_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) call AllocateVar1d(var = this%h2osoi_liqice_10cm_col, name = 'h2osoi_liqice_10cm_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) @@ -263,7 +272,9 @@ subroutine InitCold(this, bounds) character(len=256) :: locfn !----------------------------------------------------------------------- - + ! Water Stored in plants is almost always a static entity, with the exception + ! of when FATES-hydraulics is used. As such, this is trivially set to 0.0 (rgk 03-2017) + this%total_plant_stored_h2o_col(bounds%begc:bounds%endc) = 0.0_r8 do l = bounds%begl, bounds%endl if (lun%urbpoi(l)) then diff --git a/src/biogeophys/WaterFluxBulkType.F90 b/src/biogeophys/WaterFluxBulkType.F90 index 156c40491e..d0d31ade39 100644 --- a/src/biogeophys/WaterFluxBulkType.F90 +++ b/src/biogeophys/WaterFluxBulkType.F90 @@ -14,7 +14,6 @@ module WaterFluxBulkType use clm_varcon , only : spval use decompMod , only : bounds_type use CNSharedParamsMod , only : use_fun - use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_gridcell use WaterFluxType , only : waterflux_type use WaterInfoBaseType, only : water_info_base_type use WaterTracerContainerType, only : water_tracer_container_type @@ -29,12 +28,6 @@ module WaterFluxBulkType real(r8), pointer :: qflx_phs_neg_col (:) ! col sum of negative hydraulic redistribution fluxes (mm H2O/s) [+] - ! In the snow capping parametrization excess mass above h2osno_max is removed. A breakdown of mass into liquid - real(r8), pointer :: qflx_snwcp_discarded_liq_col(:) ! col excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) - real(r8), pointer :: qflx_snwcp_discarded_ice_col(:) ! col excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) - real(r8), pointer :: qflx_glcice_dyn_water_flux_col(:) ! col water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system); valid for all columns - - real(r8), pointer :: qflx_snowindunload_patch (:) ! patch canopy snow wind unloading (mm H2O /s) real(r8), pointer :: qflx_snowindunload_col (:) ! col canopy snow wind unloading (mm H2O /s) real(r8), pointer :: qflx_snotempunload_patch (:) ! patch canopy snow temp unloading (mm H2O /s) @@ -59,22 +52,9 @@ module WaterFluxBulkType real(r8), pointer :: qflx_h2osfc_drain_col (:) ! col bottom drainage from h2osfc (mm/s) real(r8), pointer :: qflx_top_soil_to_h2osfc_col(:) ! col portion of qflx_top_soil going to h2osfc, minus evaporation (mm/s) real(r8), pointer :: qflx_in_h2osfc_col(:) ! col total surface input to h2osfc - real(r8), pointer :: qflx_h2osfc_to_ice_col (:) ! col conversion of h2osfc to ice - real(r8), pointer :: qflx_snow_h2osfc_col (:) ! col snow falling on surface water real(r8), pointer :: qflx_deficit_col (:) ! col water deficit to keep non-negative liquid water content (mm H2O) real(r8), pointer :: qflx_snomelt_lyr_col (:,:) ! col snow melt in each layer (mm H2O /s) - real(r8), pointer :: qflx_snow_drain_col (:) ! col drainage from snow pack - real(r8), pointer :: qflx_ice_runoff_snwcp_col(:) ! col solid runoff from snow capping (mm H2O /s) - real(r8), pointer :: qflx_ice_runoff_xs_col (:) ! col solid runoff from excess ice in soil (mm H2O /s) real(r8), pointer :: qflx_drain_vr_col (:,:) ! col liquid water losted as drainage (m /time step) - real(r8), pointer :: snow_sources_col (:) ! col snow sources (mm H2O/s) - real(r8), pointer :: snow_sinks_col (:) ! col snow sinks (mm H2O/s) - - - ! Objects that help convert once-per-year dynamic land cover changes into fluxes - ! that are dribbled throughout the year - type(annual_flux_dribbler_type) :: qflx_liq_dynbal_dribbler - type(annual_flux_dribbler_type) :: qflx_ice_dynbal_dribbler ! ET accumulation real(r8), pointer :: AnnEt (:) ! Annual average ET flux mmH20/s @@ -137,11 +117,7 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%qflx_snowindunload_col (begp:endp)) ; this%qflx_snowindunload_col (:) = nan allocate(this%qflx_snotempunload_patch (begp:endp)) ; this%qflx_snotempunload_patch (:) = nan allocate(this%qflx_snotempunload_col (begp:endp)) ; this%qflx_snotempunload_col (:) = nan - - allocate(this%qflx_snwcp_discarded_liq_col(begc:endc)) ; this%qflx_snwcp_discarded_liq_col(:) = nan - allocate(this%qflx_snwcp_discarded_ice_col(begc:endc)) ; this%qflx_snwcp_discarded_ice_col(:) = nan - allocate(this%qflx_glcice_dyn_water_flux_col(begc:endc)) ; this%qflx_glcice_dyn_water_flux_col (:) = nan allocate(this%qflx_phs_neg_col (begc:endc)) ; this%qflx_phs_neg_col (:) = nan allocate( this%qflx_ev_snow_patch (begp:endp)) ; this%qflx_ev_snow_patch (:) = nan @@ -162,30 +138,13 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%qflx_h2osfc_drain_col (begc:endc)) ; this%qflx_h2osfc_drain_col (:) = nan allocate(this%qflx_top_soil_to_h2osfc_col(begc:endc)) ; this%qflx_top_soil_to_h2osfc_col(:) = nan allocate(this%qflx_in_h2osfc_col (begc:endc)) ; this%qflx_in_h2osfc_col(:) = nan - allocate(this%qflx_h2osfc_to_ice_col (begc:endc)) ; this%qflx_h2osfc_to_ice_col (:) = nan allocate(this%qflx_infl_excess_surf_col(begc:endc)) ; this%qflx_infl_excess_surf_col(:) = nan allocate(this%qflx_h2osfc_surf_col (begc:endc)) ; this%qflx_h2osfc_surf_col (:) = nan - allocate(this%qflx_snow_h2osfc_col (begc:endc)) ; this%qflx_snow_h2osfc_col (:) = nan allocate(this%qflx_snomelt_lyr_col (begc:endc,-nlevsno+1:0)) ; this%qflx_snomelt_lyr_col (:,:) = nan - allocate(this%qflx_snow_drain_col (begc:endc)) ; this%qflx_snow_drain_col (:) = nan allocate(this%qflx_deficit_col (begc:endc)) ; this%qflx_deficit_col (:) = nan - allocate(this%qflx_ice_runoff_snwcp_col(begc:endc)) ; this%qflx_ice_runoff_snwcp_col(:) = nan - allocate(this%qflx_ice_runoff_xs_col (begc:endc)) ; this%qflx_ice_runoff_xs_col (:) = nan - allocate(this%snow_sources_col (begc:endc)) ; this%snow_sources_col (:) = nan - allocate(this%snow_sinks_col (begc:endc)) ; this%snow_sinks_col (:) = nan allocate(this%AnnET (begc:endc)) ; this%AnnET (:) = nan - this%qflx_liq_dynbal_dribbler = annual_flux_dribbler_gridcell( & - bounds = bounds, & - name = 'qflx_liq_dynbal', & - units = 'mm H2O') - - this%qflx_ice_dynbal_dribbler = annual_flux_dribbler_gridcell( & - bounds = bounds, & - name = 'qflx_ice_dynbal', & - units = 'mm H2O') - end subroutine InitBulkAllocate !------------------------------------------------------------------------ @@ -210,21 +169,6 @@ subroutine InitBulkHistory(this, bounds) begc = bounds%begc; endc= bounds%endc begg = bounds%begg; endg= bounds%endg - this%qflx_snow_drain_col(begc:endc) = spval - call hist_addfld1d ( & - fname=this%info%fname('QFLX_SNOW_DRAIN'), & - units='mm/s', & - avgflag='A', & - long_name=this%info%lname('drainage from snow pack'), & - ptr_col=this%qflx_snow_drain_col, c2l_scale_type='urbanf') - - call hist_addfld1d ( & - fname=this%info%fname('QFLX_SNOW_DRAIN_ICE'), & - units='mm/s', & - avgflag='A', & - long_name=this%info%lname('drainage from snow pack melt (ice landunits only)'), & - ptr_col=this%qflx_snow_drain_col, c2l_scale_type='urbanf', l2g_scale_type='ice') - this%qflx_snomelt_lyr_col(begc:endc,-nlevsno+1:0) = spval data2dptr => this%qflx_snomelt_lyr_col(begc:endc,-nlevsno+1:0) call hist_addfld2d ( & @@ -243,14 +187,6 @@ subroutine InitBulkHistory(this, bounds) l2g_scale_type='ice', default='inactive') - this%qflx_h2osfc_to_ice_col(begc:endc) = spval - call hist_addfld1d ( & - fname=this%info%fname('QH2OSFC_TO_ICE'), & - units='mm/s', & - avgflag='A', & - long_name=this%info%lname('surface water converted to ice'), & - ptr_col=this%qflx_h2osfc_to_ice_col, default='inactive') - call hist_addfld2d ( & fname=this%info%fname('QROOTSINK'), & units='mm/s', type2d='levsoi', & @@ -292,28 +228,6 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('net negative hydraulic redistribution flux'), & ptr_col=this%qflx_phs_neg_col, default='inactive') - ! As defined here, snow_sources - snow_sinks will equal the change in h2osno at any - ! given time step but only if there is at least one snow layer (for all landunits - ! except lakes). Also note that monthly average files of snow_sources and snow_sinks - ! sinks must be weighted by number of days in the month to diagnose, for example, an - ! annual value of the change in h2osno. - - this%snow_sources_col(begc:endc) = spval - call hist_addfld1d ( & - fname=this%info%fname('SNOW_SOURCES'), & - units='mm/s', & - avgflag='A', & - long_name=this%info%lname('snow sources (liquid water)'), & - ptr_col=this%snow_sources_col, c2l_scale_type='urbanf') - - this%snow_sinks_col(begc:endc) = spval - call hist_addfld1d ( & - fname=this%info%fname('SNOW_SINKS'), & - units='mm/s', & - avgflag='A', & - long_name=this%info%lname('snow sinks (liquid water)'), & - ptr_col=this%snow_sinks_col, c2l_scale_type='urbanf') - this%AnnET(begc:endc) = spval call hist_addfld1d ( & fname=this%info%fname('AnnET'), & @@ -340,18 +254,6 @@ subroutine InitBulkCold(this, bounds) this%qflx_phs_neg_col(bounds%begc:bounds%endc) = 0.0_r8 this%qflx_h2osfc_surf_col(bounds%begc:bounds%endc) = 0._r8 - this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 - - ! This variable only gets set in the hydrology filter; need to initialize it to 0 for - ! the sake of columns outside this filter - this%qflx_ice_runoff_xs_col(bounds%begc:bounds%endc) = 0._r8 - - ! Initialize qflx_glcice_dyn_water_flux_col to 0 for all columns because we want this - ! flux to remain 0 for columns where is is never set, including non-glacier columns. - ! - ! Other qflx_glcice fluxes intentionally remain unset (spval) outside the do_smb - ! filter, so that they are flagged as missing value outside that filter. - this%qflx_glcice_dyn_water_flux_col(bounds%begc:bounds%endc) = 0._r8 end subroutine InitBulkCold @@ -479,26 +381,9 @@ subroutine RestartBulk(this, bounds, ncid, flag) character(len=*) , intent(in) :: flag ! 'read' or 'write' ! ! !LOCAL VARIABLES: - logical :: readvar ! determine if variable is on initial file !----------------------------------------------------------------------- call this%restart ( bounds, ncid, flag=flag ) - ! needed for SNICAR - - call restartvar(ncid=ncid, flag=flag, & - varname=this%info%fname('qflx_snow_drain')//':'//this%info%fname('qflx_snow_melt'), & - xtype=ncd_double, & - dim1name='column', & - long_name=this%info%lname('drainage from snow column'), & - units='mm/s', & - interpinic_flag='interp', readvar=readvar, data=this%qflx_snow_drain_col) - if (flag == 'read' .and. .not. readvar) then - ! initial run, not restart: initialize qflx_snow_drain to zero - this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 - endif - - call this%qflx_liq_dynbal_dribbler%Restart(bounds, ncid, flag) - call this%qflx_ice_dynbal_dribbler%Restart(bounds, ncid, flag) end subroutine RestartBulk end module WaterFluxBulkType diff --git a/src/biogeophys/WaterFluxType.F90 b/src/biogeophys/WaterFluxType.F90 index c1fe47354a..a477eea874 100644 --- a/src/biogeophys/WaterFluxType.F90 +++ b/src/biogeophys/WaterFluxType.F90 @@ -14,6 +14,7 @@ module WaterFluxType use decompMod , only : BOUNDS_SUBGRID_PATCH, BOUNDS_SUBGRID_COLUMN, BOUNDS_SUBGRID_GRIDCELL use LandunitType , only : lun use ColumnType , only : col + use AnnualFluxDribbler, only : annual_flux_dribbler_type, annual_flux_dribbler_gridcell use WaterInfoBaseType, only : water_info_base_type use WaterTracerContainerType, only : water_tracer_container_type use WaterTracerUtils, only : AllocateVar1d, AllocateVar2d @@ -51,10 +52,12 @@ module WaterFluxType ! and solid fluxes is done, these are represented by qflx_snwcp_liq_col and qflx_snwcp_ice_col. real(r8), pointer :: qflx_snwcp_liq_col (:) ! col excess liquid h2o due to snow capping (outgoing) (mm H2O /s) real(r8), pointer :: qflx_snwcp_ice_col (:) ! col excess solid h2o due to snow capping (outgoing) (mm H2O /s) + real(r8), pointer :: qflx_snwcp_discarded_liq_col(:) ! col excess liquid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) + real(r8), pointer :: qflx_snwcp_discarded_ice_col(:) ! col excess solid h2o due to snow capping, which we simply discard in order to reset the snow pack (mm H2O /s) real(r8), pointer :: qflx_glcice_col(:) ! col net flux of new glacial ice (growth - melt) (mm H2O/s), passed to GLC; only valid inside the do_smb_c filter real(r8), pointer :: qflx_glcice_frz_col (:) ! col ice growth (positive definite) (mm H2O/s); only valid inside the do_smb_c filter real(r8), pointer :: qflx_glcice_melt_col(:) ! col ice melt (positive definite) (mm H2O/s); only valid inside the do_smb_c filter - + real(r8), pointer :: qflx_glcice_dyn_water_flux_col(:) ! col water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system); valid for all columns real(r8), pointer :: qflx_tran_veg_patch (:) ! patch vegetation transpiration (mm H2O/s) (+ = to atm) real(r8), pointer :: qflx_tran_veg_col (:) ! col vegetation transpiration (mm H2O/s) (+ = to atm) @@ -80,6 +83,12 @@ module WaterFluxType real(r8), pointer :: qflx_rsub_sat_col (:) ! col soil saturation excess [mm/s] real(r8), pointer :: qflx_snofrz_lyr_col (:,:) ! col snow freezing rate (positive definite) (col,lyr) [kg m-2 s-1] real(r8), pointer :: qflx_snofrz_col (:) ! col column-integrated snow freezing rate (positive definite) (col) [kg m-2 s-1] + real(r8), pointer :: qflx_snow_drain_col (:) ! col drainage from snow pack + real(r8), pointer :: qflx_ice_runoff_snwcp_col(:) ! col solid runoff from snow capping (mm H2O /s) + real(r8), pointer :: qflx_ice_runoff_xs_col (:) ! col solid runoff from excess ice in soil (mm H2O /s) + + real(r8), pointer :: qflx_h2osfc_to_ice_col (:) ! col conversion of h2osfc to ice + real(r8), pointer :: qflx_snow_h2osfc_col (:) ! col snow falling on surface water ! Dynamic land cover change real(r8), pointer :: qflx_liq_dynbal_grc (:) ! grc liq dynamic land cover change conversion runoff flux @@ -88,6 +97,10 @@ module WaterFluxType real(r8), pointer :: qflx_irrig_patch (:) ! patch irrigation flux (mm H2O/s) [+] real(r8), pointer :: qflx_irrig_col (:) ! col irrigation flux (mm H2O/s) [+] + ! Objects that help convert once-per-year dynamic land cover changes into fluxes + ! that are dribbled throughout the year + type(annual_flux_dribbler_type) :: qflx_liq_dynbal_dribbler + type(annual_flux_dribbler_type) :: qflx_ice_dynbal_dribbler contains @@ -184,6 +197,12 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%qflx_snwcp_ice_col, name = 'qflx_snwcp_ice_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_snwcp_discarded_liq_col, name = 'qflx_snwcp_discarded_liq_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_snwcp_discarded_ice_col, name = 'qflx_snwcp_discarded_ice_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) call AllocateVar1d(var = this%qflx_glcice_col, name = 'qflx_glcice_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) @@ -193,6 +212,9 @@ subroutine InitAllocate(this, bounds, tracer_vars) call AllocateVar1d(var = this%qflx_glcice_melt_col, name = 'qflx_glcice_melt_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_glcice_dyn_water_flux_col, name = 'qflx_glcice_dyn_water_flux_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) call AllocateVar1d(var = this%qflx_tran_veg_col, name = 'qflx_tran_veg_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) @@ -258,6 +280,15 @@ subroutine InitAllocate(this, bounds, tracer_vars) container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN, & dim2beg = -nlevsno+1, dim2end = 0) + call AllocateVar1d(var = this%qflx_snow_drain_col, name = 'qflx_snow_drain_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_ice_runoff_snwcp_col, name = 'qflx_ice_runoff_snwcp_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_ice_runoff_xs_col, name = 'qflx_ice_runoff_xs_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) call AllocateVar1d(var = this%qflx_qrgwl_col, name = 'qflx_qrgwl_col', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) @@ -280,6 +311,13 @@ subroutine InitAllocate(this, bounds, tracer_vars) container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_h2osfc_to_ice_col, name = 'qflx_h2osfc_to_ice_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_snow_h2osfc_col, name = 'qflx_snow_h2osfc_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%qflx_liq_dynbal_grc, name = 'qflx_liq_dynbal_grc', & container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_GRIDCELL) @@ -294,6 +332,16 @@ subroutine InitAllocate(this, bounds, tracer_vars) container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + this%qflx_liq_dynbal_dribbler = annual_flux_dribbler_gridcell( & + bounds = bounds, & + name = this%info%fname('qflx_liq_dynbal'), & + units = 'mm H2O') + + this%qflx_ice_dynbal_dribbler = annual_flux_dribbler_gridcell( & + bounds = bounds, & + name = this%info%fname('qflx_ice_dynbal'), & + units = 'mm H2O') + end subroutine InitAllocate !------------------------------------------------------------------------ @@ -461,6 +509,21 @@ subroutine InitHistory(this, bounds) ptr_col=data2dptr, c2l_scale_type='urbanf',no_snow_behavior=no_snow_normal, & l2g_scale_type='ice', default='inactive') + this%qflx_snow_drain_col(begc:endc) = spval + call hist_addfld1d ( & + fname=this%info%fname('QFLX_SNOW_DRAIN'), & + units='mm/s', & + avgflag='A', & + long_name=this%info%lname('drainage from snow pack'), & + ptr_col=this%qflx_snow_drain_col, c2l_scale_type='urbanf') + + call hist_addfld1d ( & + fname=this%info%fname('QFLX_SNOW_DRAIN_ICE'), & + units='mm/s', & + avgflag='A', & + long_name=this%info%lname('drainage from snow pack melt (ice landunits only)'), & + ptr_col=this%qflx_snow_drain_col, c2l_scale_type='urbanf', l2g_scale_type='ice') + this%qflx_prec_intr_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('QINTR'), & @@ -636,6 +699,14 @@ subroutine InitHistory(this, bounds) long_name=this%info%lname('saturation excess drainage'), & ptr_col=this%qflx_rsub_sat_col, c2l_scale_type='urbanf') + this%qflx_h2osfc_to_ice_col(begc:endc) = spval + call hist_addfld1d ( & + fname=this%info%fname('QH2OSFC_TO_ICE'), & + units='mm/s', & + avgflag='A', & + long_name=this%info%lname('surface water converted to ice'), & + ptr_col=this%qflx_h2osfc_to_ice_col, default='inactive') + this%qflx_irrig_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('QIRRIG'), & @@ -670,6 +741,18 @@ subroutine InitCold(this, bounds) this%qflx_evap_grnd_col(bounds%begc:bounds%endc) = 0.0_r8 this%qflx_dew_grnd_col (bounds%begc:bounds%endc) = 0.0_r8 this%qflx_dew_snow_col (bounds%begc:bounds%endc) = 0.0_r8 + this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 + + ! This variable only gets set in the hydrology filter; need to initialize it to 0 for + ! the sake of columns outside this filter + this%qflx_ice_runoff_xs_col(bounds%begc:bounds%endc) = 0._r8 + + ! Initialize qflx_glcice_dyn_water_flux_col to 0 for all columns because we want this + ! flux to remain 0 for columns where is is never set, including non-glacier columns. + ! + ! Other qflx_glcice fluxes intentionally remain unset (spval) outside the do_smb + ! filter, so that they are flagged as missing value outside that filter. + this%qflx_glcice_dyn_water_flux_col(bounds%begc:bounds%endc) = 0._r8 ! needed for CNNLeaching do c = bounds%begc, bounds%endc @@ -690,7 +773,7 @@ subroutine Restart(this, bounds, ncid, flag) use restUtilMod ! ! !ARGUMENTS: - class(waterflux_type), intent(in):: this + class(waterflux_type), intent(inout) :: this type(bounds_type), intent(in) :: bounds type(file_desc_t), intent(inout) :: ncid ! netcdf id character(len=*) , intent(in) :: flag ! 'read' or 'write' @@ -712,6 +795,21 @@ subroutine Restart(this, bounds, ncid, flag) this%qflx_snofrz_lyr_col(bounds%begc:bounds%endc,-nlevsno+1:0) = 0._r8 endif + call restartvar(ncid=ncid, flag=flag, & + varname=this%info%fname('qflx_snow_drain')//':'//this%info%fname('qflx_snow_melt'), & + xtype=ncd_double, & + dim1name='column', & + long_name=this%info%lname('drainage from snow column'), & + units='mm/s', & + interpinic_flag='interp', readvar=readvar, data=this%qflx_snow_drain_col) + if (flag == 'read' .and. .not. readvar) then + ! initial run, not restart: initialize qflx_snow_drain to zero + this%qflx_snow_drain_col(bounds%begc:bounds%endc) = 0._r8 + endif + + call this%qflx_liq_dynbal_dribbler%Restart(bounds, ncid, flag) + call this%qflx_ice_dynbal_dribbler%Restart(bounds, ncid, flag) + end subroutine Restart end module WaterFluxType diff --git a/src/biogeophys/WaterInfoBaseType.F90 b/src/biogeophys/WaterInfoBaseType.F90 index 7ac575ee71..f0fa1bc420 100644 --- a/src/biogeophys/WaterInfoBaseType.F90 +++ b/src/biogeophys/WaterInfoBaseType.F90 @@ -25,6 +25,9 @@ module WaterInfoBaseType ! Get a history/restart long name for this tracer (or bulk) procedure(lname_interface), public, deferred :: lname + + procedure(is_communicated_with_coupler_interface), public, deferred :: is_communicated_with_coupler + procedure :: get_ratio procedure :: set_metadata @@ -52,6 +55,16 @@ pure function lname_interface(this, basename) result(lname) class(water_info_base_type) , intent(in) :: this character(len=*) , intent(in) :: basename end function lname_interface + + pure function is_communicated_with_coupler_interface(this) result(coupled) + ! Returns true if this tracer is received from and sent to the coupler. Returns + ! false if this tracer is just used internally in CTSM, and is set to some fixed + ! ratio times the bulk water. + import :: water_info_base_type + logical :: coupled + class(water_info_base_type), intent(in) :: this + end function is_communicated_with_coupler_interface + end interface contains diff --git a/src/biogeophys/WaterInfoBulkType.F90 b/src/biogeophys/WaterInfoBulkType.F90 index 5ced69112b..ac942234a4 100644 --- a/src/biogeophys/WaterInfoBulkType.F90 +++ b/src/biogeophys/WaterInfoBulkType.F90 @@ -20,6 +20,7 @@ module WaterInfoBulkType contains procedure, public :: fname ! Get a history/restart field name procedure, public :: lname ! Get a history/restart long name + procedure, public :: is_communicated_with_coupler end type water_info_bulk_type interface water_info_bulk_type @@ -83,4 +84,12 @@ pure function lname(this, basename) end function lname + !----------------------------------------------------------------------- + pure function is_communicated_with_coupler(this) result(coupled) + logical :: coupled + class(water_info_bulk_type), intent(in) :: this + + coupled = .true. + end function is_communicated_with_coupler + end module WaterInfoBulkType diff --git a/src/biogeophys/WaterInfoIsotopeType.F90 b/src/biogeophys/WaterInfoIsotopeType.F90 index cdf42e3cee..1f387918a6 100644 --- a/src/biogeophys/WaterInfoIsotopeType.F90 +++ b/src/biogeophys/WaterInfoIsotopeType.F90 @@ -38,17 +38,21 @@ module WaterInfoIsotopeType contains - function constructor(tracer_name,ratio) result(this) + function constructor(tracer_name, ratio, communicated_with_coupler) result(this) ! Create a water_info_isotope_type object ! ! Eventually, this will either (a) accept various arguments specifying information ! about this isotope (molecular weight, diffusivity ratio, etc.), or (b) look up this ! information from a lookup table defined here or elsewhere, based on the tracer_name. type(water_info_isotope_type) :: this ! function result - character(len=*), intent(in) :: tracer_name + character(len=*), intent(in) :: tracer_name real(r8), intent(in) :: ratio + logical , intent(in) :: communicated_with_coupler ! see documentation in WaterInfoTracerType.F90 - this%water_info_tracer_type = water_info_tracer_type(tracer_name,ratio) + this%water_info_tracer_type = water_info_tracer_type( & + tracer_name = tracer_name, & + ratio = ratio, & + communicated_with_coupler = communicated_with_coupler) end function constructor end module WaterInfoIsotopeType diff --git a/src/biogeophys/WaterInfoTracerType.F90 b/src/biogeophys/WaterInfoTracerType.F90 index e5d8dbd93e..da03f41223 100644 --- a/src/biogeophys/WaterInfoTracerType.F90 +++ b/src/biogeophys/WaterInfoTracerType.F90 @@ -18,9 +18,15 @@ module WaterInfoTracerType type, extends(water_info_base_type), public :: water_info_tracer_type private character(len=:), allocatable :: tracer_name + + ! If true, this tracer is received from and sent to the coupler. If false, this + ! tracer is just used internally in CTSM, and is set to some fixed ratio times the + ! bulk water. + logical :: communicated_with_coupler contains procedure, public :: fname ! Get a history/restart field name for this tracer procedure, public :: lname ! Get a history/restart long name for this tracer + procedure, public :: is_communicated_with_coupler end type water_info_tracer_type interface water_info_tracer_type @@ -32,13 +38,19 @@ module WaterInfoTracerType contains - function constructor(tracer_name,ratio) result(this) + function constructor(tracer_name, ratio, communicated_with_coupler) result(this) ! Create a water_info_tracer_type object type(water_info_tracer_type) :: this ! function result character(len=*), intent(in) :: tracer_name real(r8), intent(in) :: ratio + ! If true, this tracer is received from and sent to the coupler. If false, this tracer + ! is just used internally in CTSM, and is set to some fixed ratio times the bulk + ! water. + logical, intent(in) :: communicated_with_coupler + this%tracer_name = trim(tracer_name) + this%communicated_with_coupler = communicated_with_coupler call this%set_metadata(ratio = ratio) end function constructor @@ -86,4 +98,26 @@ pure function lname(this, basename) end function lname + !----------------------------------------------------------------------- + pure function is_communicated_with_coupler(this) result(coupled) + ! + ! !DESCRIPTION: + ! Returns true if this tracer is received from and sent to the coupler. Returns false + ! if this tracer is just used internally in CTSM, and is set to some fixed ratio + ! times the bulk water. + ! + ! !ARGUMENTS: + logical :: coupled ! function result + class(water_info_tracer_type), intent(in) :: this + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'is_communicated_with_coupler' + !----------------------------------------------------------------------- + + coupled = this%communicated_with_coupler + + end function is_communicated_with_coupler + + end module WaterInfoTracerType diff --git a/src/biogeophys/WaterStateType.F90 b/src/biogeophys/WaterStateType.F90 index ab62bf141f..ffdcac7cdd 100644 --- a/src/biogeophys/WaterStateType.F90 +++ b/src/biogeophys/WaterStateType.F90 @@ -10,6 +10,7 @@ module WaterStateType ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun use decompMod , only : bounds_type use decompMod , only : BOUNDS_SUBGRID_PATCH, BOUNDS_SUBGRID_COLUMN use clm_varctl , only : iulog, use_bedrock @@ -43,11 +44,11 @@ module WaterStateType contains - procedure :: Init - procedure :: Restart - procedure, private :: InitAllocate - procedure, private :: InitHistory - procedure, private :: InitCold + procedure, public :: Init + procedure, public :: Restart + procedure, private :: InitAllocate + procedure, private :: InitHistory + procedure, private :: InitCold end type waterstate_type @@ -389,6 +390,10 @@ subroutine InitCold(this, bounds, & do j = 1, nlevs this%h2osoi_vol_col(c,j) = 1.0_r8 end do + else + write(iulog,*) 'water_state_type InitCold: unhandled landunit type ', lun%itype(l) + call endrun(msg = 'unhandled landunit type', & + additional_msg = errMsg(sourcefile, __LINE__)) endif do j = 1, nlevs this%h2osoi_vol_col(c,j) = min(this%h2osoi_vol_col(c,j), watsat_col(c,j)) diff --git a/src/biogeophys/WaterTracerUtils.F90 b/src/biogeophys/WaterTracerUtils.F90 index 079ab450da..6719000fd8 100644 --- a/src/biogeophys/WaterTracerUtils.F90 +++ b/src/biogeophys/WaterTracerUtils.F90 @@ -22,6 +22,8 @@ module WaterTracerUtils ! !PUBLIC MEMBER FUNCTIONS: public :: AllocateVar1d public :: AllocateVar2d + public :: CalcTracerFromBulk + public :: CalcTracerFromBulkFixedRatio public :: CompareBulkToTracer ! !PRIVATE MEMBER DATA: @@ -120,6 +122,113 @@ subroutine AllocateVar2d(var, name, container, bounds, subgrid_level, dim2beg, d end subroutine AllocateVar2d + !----------------------------------------------------------------------- + subroutine CalcTracerFromBulk(lb, num_pts, filter_pts, & + bulk_source, bulk_val, tracer_source, tracer_val) + ! + ! !DESCRIPTION: + ! Calculate a tracer variable from a corresponding bulk variable when the ratio of + ! the tracer to the bulk should be based on the ratio in some 'source' variable. + ! + ! Typically this is used to calculate a tracer flux given a bulk flux. In this case, + ! the source variable should be the source of the flux - since the tracer ratio of + ! the flux will be the same as the tracer ratio in the source state variable. + ! + ! tracer_val will be updated within the given filter, and will remain at its original + ! values elsewhere + ! + ! !ARGUMENTS: + integer , intent(in) :: lb ! lower bound for arrays + integer , intent(in) :: num_pts ! number of points in the filter + integer , intent(in) :: filter_pts(:) ! filter in which tracer_val should be updated + real(r8) , intent(in) :: bulk_source(lb:) ! values of the source for this variable, for bulk + real(r8) , intent(in) :: bulk_val(lb:) ! values of the variable of interest, for bulk + real(r8) , intent(in) :: tracer_source(lb:) ! values of the source for this variable, for the tracer + real(r8) , intent(inout) :: tracer_val(lb:) ! output values of the variable of interest, for the tracer + ! + ! !LOCAL VARIABLES: + integer :: num + integer :: fn, n + + character(len=*), parameter :: subname = 'CalcTracerFromBulk' + !----------------------------------------------------------------------- + + num = size(bulk_val) + SHR_ASSERT((size(bulk_source) == num), errMsg(sourcefile, __LINE__)) + SHR_ASSERT((size(tracer_source) == num), errMsg(sourcefile, __LINE__)) + SHR_ASSERT((size(tracer_val) == num), errMsg(sourcefile, __LINE__)) + + do fn = 1, num_pts + n = filter_pts(fn) + + if (bulk_source(n) /= 0._r8) then + ! Standard case + tracer_val(n) = bulk_val(n) * (tracer_source(n) / bulk_source(n)) + else if (bulk_val(n) == 0._r8 .and. tracer_source(n) == 0._r8) then + ! This is acceptable: bulk_source, bulk_val and tracer_source are all 0 + tracer_val(n) = 0._r8 + else if (bulk_val(n) /= 0._r8) then + write(iulog,*) subname//' ERROR: Non-zero bulk val despite zero bulk source:' + write(iulog,*) 'bulk_val = ', bulk_val(n) + write(iulog,*) 'at n = ', n + write(iulog,*) 'This would lead to an indeterminate tracer val.' + call endrun(msg=subname//': Non-zero bulk val despite zero bulk source', & + additional_msg=errMsg(sourcefile, __LINE__)) + else if (tracer_source(n) /= 0._r8) then + ! NOTE(wjs, 2018-09-28) To avoid this error, we might need code elsewhere that + ! sets tracer state variables to 0 if the corresponding bulk state variable is 0 + ! and the tracer state is originally near 0 (within roundoff) - in order deal + ! with roundoff issues arising during state updates. (There's a bit of + ! discussion of this point in https://github.com/ESCOMP/ctsm/issues/487.) + write(iulog,*) subname//' ERROR: Non-zero tracer source despite zero bulk source:' + write(iulog,*) 'tracer_source = ', tracer_source(n) + write(iulog,*) 'at n = ', n + write(iulog,*) 'This would lead to an indeterminate tracer val.' + call endrun(msg=subname//': Non-zero tracer source despite zero bulk source', & + additional_msg=errMsg(sourcefile, __LINE__)) + else + write(iulog,*) subname//' ERROR: unhandled condition; we should never get here.' + write(iulog,*) 'This indicates a programming error in this subroutine.' + write(iulog,*) 'bulk_val = ', bulk_val(n) + write(iulog,*) 'bulk_source = ', bulk_source(n) + write(iulog,*) 'tracer_source = ', tracer_source(n) + write(iulog,*) 'at n = ', n + call endrun(msg=subname//': unhandled condition; we should never get here', & + additional_msg=errMsg(sourcefile, __LINE__)) + end if + end do + + end subroutine CalcTracerFromBulk + + + !----------------------------------------------------------------------- + subroutine CalcTracerFromBulkFixedRatio(bulk, ratio, tracer) + ! + ! !DESCRIPTION: + ! Calculate a tracer variable from a corresponding bulk variable when the tracer + ! should be a fixed ratio times the bulk + ! + ! !ARGUMENTS: + real(r8), intent(in) :: bulk(:) + real(r8), intent(in) :: ratio ! ratio of tracer to bulk + real(r8), intent(inout) :: tracer(:) + ! + ! !LOCAL VARIABLES: + integer :: num + integer :: i + + character(len=*), parameter :: subname = 'CalcTracerFromBulkFixedRatio' + !----------------------------------------------------------------------- + + num = size(bulk) + SHR_ASSERT((size(tracer) == num), errMsg(sourcefile, __LINE__)) + do i = 1, num + tracer(i) = bulk(i) * ratio + end do + + end subroutine CalcTracerFromBulkFixedRatio + + !----------------------------------------------------------------------- subroutine CompareBulkToTracer(bounds_beg, bounds_end, & bulk, tracer, caller_location, name) diff --git a/src/biogeophys/WaterType.F90 b/src/biogeophys/WaterType.F90 index 20d20cb883..fcaa8edbcf 100644 --- a/src/biogeophys/WaterType.F90 +++ b/src/biogeophys/WaterType.F90 @@ -5,46 +5,79 @@ module WaterType ! Container for derived types relating to water, both for bulk water and for isotopes ! and other tracers. ! - ! To loop through all tracers, use code like this: - ! do i = 1, water_inst%num_tracers - ! call some_subroutine(..., water_inst%waterflux_tracer_inst(i), ...) + ! Variables pertaining to bulk water can be accessed in two ways: + ! + ! (1) Using water_inst%water*bulk_inst + ! + ! (2) As one of the indices in water_inst%bulk_and_tracers(:)%water*_inst + ! + ! Method (1) is greatly preferable when you are just operating on bulk water. Method + ! (2) is just meant to be used when you are doing the same operation on bulk water + ! and all water tracers. + ! + ! To loop through bulk and all tracers, use code like this: + ! do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end + ! associate( & + ! waterflux_inst => water_inst%bulk_and_tracers(i)%waterflux_inst, & + ! [and other associations, as necessary]) + ! [Do calculations involving waterflux_inst, etc.] + ! end associate ! end do ! - ! To loop through all isotopes, use code like this: + ! To loop through all tracers (not bulk), use code like this: + ! do i = water_inst%tracers_beg, water_inst%tracers_end + ! associate( & + ! waterflux_inst => water_inst%bulk_and_tracers(i)%waterflux_inst, & + ! [and other associations, as necessary]) + ! [Do calculations involving waterflux_inst, etc.] + ! end associate + ! end do + ! + ! To loop through all isotopes (not bulk or other water tracers), use code like this: ! type(water_info_isotope_type), pointer :: iso_info ! - ! do i = 1, water_inst%num_tracers + ! do i = water_inst%tracers_beg, water_inst%tracers_end ! if (water_inst%IsIsotope(i)) then ! call water_inst%GetIsotopeInfo(i, iso_info) - ! call some_subroutine(..., iso_info, water_inst%waterflux_tracer_inst(i), ...) + ! associate( & + ! waterflux_inst => water_inst%bulk_and_tracers(i)%waterflux_inst, & + ! [and other associations, as necessary]) + ! [Do calculations involving iso_info, waterflux_inst, etc.] + ! end associate ! end if ! end do ! + ! The associate statements given above aren't crucial. If the block of code refers to + ! multiple instances (waterstate, waterflux, etc.), but only refers to each one once or + ! twice, it can be best to just have: + ! associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) + ! ! !USES: #include "shr_assert.h" - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_log_mod , only : errMsg => shr_log_errMsg - use abortutils , only : endrun - use decompMod , only : bounds_type - use clm_varctl , only : iulog - use clm_varpar , only : nlevsno - use ncdio_pio , only : file_desc_t - use WaterFluxBulkType , only : waterfluxbulk_type - use WaterFluxType , only : waterflux_type - use WaterStateBulkType , only : waterstatebulk_type - use WaterStateType , only : waterstate_type - use WaterDiagnosticType , only : waterdiagnostic_type - use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type - use WaterBalanceType , only : waterbalance_type - use WaterInfoBaseType , only : water_info_base_type - use WaterInfoBulkType , only : water_info_bulk_type - use WaterInfoIsotopeType , only : water_info_isotope_type - use Waterlnd2atmType , only : waterlnd2atm_type - use Waterlnd2atmBulkType , only : waterlnd2atmbulk_type - use Wateratm2lndType , only : wateratm2lnd_type - use Wateratm2lndBulkType , only : wateratm2lndbulk_type - use WaterTracerContainerType, only : water_tracer_container_type - use WaterTracerUtils, only : CompareBulkToTracer + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg + use abortutils , only : endrun + use decompMod , only : bounds_type + use clm_varctl , only : iulog + use clm_varpar , only : nlevsno + use ncdio_pio , only : file_desc_t + use WaterFluxBulkType , only : waterfluxbulk_type + use WaterFluxType , only : waterflux_type + use WaterStateBulkType , only : waterstatebulk_type + use WaterStateType , only : waterstate_type + use WaterDiagnosticType , only : waterdiagnostic_type + use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type + use WaterBalanceType , only : waterbalance_type + use WaterInfoBaseType , only : water_info_base_type + use WaterInfoBulkType , only : water_info_bulk_type + use WaterInfoTracerType , only : water_info_tracer_type + use WaterInfoIsotopeType , only : water_info_isotope_type + use Waterlnd2atmType , only : waterlnd2atm_type + use Waterlnd2atmBulkType , only : waterlnd2atmbulk_type + use Wateratm2lndType , only : wateratm2lnd_type + use Wateratm2lndBulkType , only : wateratm2lndbulk_type + use WaterTracerContainerType , only : water_tracer_container_type + use WaterTracerUtils , only : CompareBulkToTracer implicit none private @@ -52,14 +85,30 @@ module WaterType ! ! !PRIVATE TYPES: - ! This type is a container for objects of class water_info_base_type, to facilitate - ! having an array of polymorphic entities. - type, private :: water_info_container_type + ! This type holds instances needed for bulk water or for a single tracer + type, private :: bulk_or_tracer_type private - ! 'info' needs to be a pointer so other pointers can point to it (since a derived - ! type component cannot have the target attribute) - class(water_info_base_type), pointer :: info - end type water_info_container_type + + ! ------------------------------------------------------------------------ + ! Public data members + ! ------------------------------------------------------------------------ + + class(waterflux_type) , pointer, public :: waterflux_inst + class(waterstate_type) , pointer, public :: waterstate_inst + class(waterdiagnostic_type) , pointer, public :: waterdiagnostic_inst + class(waterbalance_type) , pointer, public :: waterbalance_inst + class(waterlnd2atm_type) , pointer, public :: waterlnd2atm_inst + class(wateratm2lnd_type) , pointer, public :: wateratm2lnd_inst + + ! ------------------------------------------------------------------------ + ! Private data members + ! ------------------------------------------------------------------------ + + logical :: is_isotope = .false. + class(water_info_base_type) , pointer :: info + type(water_tracer_container_type) :: vars + + end type bulk_or_tracer_type ! ! !PUBLIC TYPES: @@ -82,35 +131,27 @@ module WaterType ! Public data members ! ------------------------------------------------------------------------ - integer, public :: num_tracers + ! indices into the bulk_and_tracers array + integer, public :: bulk_and_tracers_beg ! first index when looping over bulk & tracers + integer, public :: bulk_and_tracers_end ! last index when looping over bulk & tracers + integer, public :: tracers_beg ! first index when looping over just tracers + integer, public :: tracers_end ! last index when looping over just tracers + integer, public :: i_bulk ! index of bulk in arrays that contain both bulk and tracers - type(waterfluxbulk_type), public :: waterfluxbulk_inst - type(waterstatebulk_type), public :: waterstatebulk_inst - type(waterdiagnosticbulk_type), public :: waterdiagnosticbulk_inst - type(waterbalance_type), public :: waterbalancebulk_inst - type(waterlnd2atmbulk_type), public :: waterlnd2atmbulk_inst - type(wateratm2lndbulk_type), public :: wateratm2lndbulk_inst + type(waterfluxbulk_type), pointer, public :: waterfluxbulk_inst + type(waterstatebulk_type), pointer, public :: waterstatebulk_inst + type(waterdiagnosticbulk_type), pointer, public :: waterdiagnosticbulk_inst + type(waterbalance_type), pointer, public :: waterbalancebulk_inst + type(waterlnd2atmbulk_type), pointer, public :: waterlnd2atmbulk_inst + type(wateratm2lndbulk_type), pointer, public :: wateratm2lndbulk_inst - type(waterflux_type), allocatable, public :: waterflux_tracer_inst(:) - type(waterstate_type), allocatable, public :: waterstate_tracer_inst(:) - type(waterdiagnostic_type), allocatable, public :: waterdiagnostic_tracer_inst(:) - type(waterbalance_type), allocatable, public :: waterbalance_tracer_inst(:) - type(waterlnd2atm_type), allocatable, public :: waterlnd2atm_tracer_inst(:) - type(wateratm2lnd_type), allocatable, public :: wateratm2lnd_tracer_inst(:) + type(bulk_or_tracer_type), allocatable, public :: bulk_and_tracers(:) ! ------------------------------------------------------------------------ ! Private data members ! ------------------------------------------------------------------------ type(water_params_type) :: params - - ! bulk_info needs to be a pointer so other pointers can point to it (since a derived - ! type component cannot have the target attribute) - type(water_info_bulk_type), pointer :: bulk_info - type(water_tracer_container_type) :: bulk_vars ! water tracer variables for bulk water (note that this only includes variables that are also included for water tracers) - logical, allocatable :: is_isotope(:) - type(water_info_container_type), allocatable :: tracer_info(:) - type(water_tracer_container_type), allocatable :: tracer_vars(:) integer :: bulk_tracer_index ! index of the tracer that replicates bulk water (-1 if it doesn't exist) contains @@ -131,6 +172,8 @@ module WaterType procedure, private :: DoInit procedure, private :: ReadNamelist procedure, private :: SetupTracerInfo + procedure, private :: AllocateBulk + procedure, private :: AllocateTracer end type water_type interface water_params_type @@ -253,84 +296,84 @@ subroutine DoInit(this, bounds, & SHR_ASSERT_ALL((ubound(watsat_col, 1) == endc), errMsg(sourcefile, __LINE__)) SHR_ASSERT_ALL((ubound(t_soisno_col, 1) == endc), errMsg(sourcefile, __LINE__)) - allocate(this%bulk_info, source = water_info_bulk_type()) - call this%bulk_vars%init() + call this%SetupTracerInfo() + + call this%AllocateBulk() + + associate( & + bulk_info => this%bulk_and_tracers(this%i_bulk)%info, & + bulk_vars => this%bulk_and_tracers(this%i_bulk)%vars & + ) + + call bulk_vars%init() call this%waterstatebulk_inst%InitBulk(bounds, & - this%bulk_info, & - this%bulk_vars, & + bulk_info, & + bulk_vars, & h2osno_input_col = h2osno_col(begc:endc), & watsat_col = watsat_col(begc:endc, 1:), & t_soisno_col = t_soisno_col(begc:endc, -nlevsno+1:) ) call this%waterdiagnosticbulk_inst%InitBulk(bounds, & - this%bulk_info, & - this%bulk_vars, & + bulk_info, & + bulk_vars, & snow_depth_input_col = snow_depth_col(begc:endc), & waterstatebulk_inst = this%waterstatebulk_inst ) call this%waterbalancebulk_inst%Init(bounds, & - this%bulk_info, & - this%bulk_vars) + bulk_info, & + bulk_vars) call this%waterfluxbulk_inst%InitBulk(bounds, & - this%bulk_info, & - this%bulk_vars) + bulk_info, & + bulk_vars) call this%waterlnd2atmbulk_inst%InitBulk(bounds, & - this%bulk_info, & - this%bulk_vars) + bulk_info, & + bulk_vars) call this%wateratm2lndbulk_inst%InitBulk(bounds, & - this%bulk_info, & - this%bulk_vars) + bulk_info, & + bulk_vars) - call this%bulk_vars%complete_setup() + call bulk_vars%complete_setup() - call this%SetupTracerInfo() + end associate - if (this%num_tracers > 0) then - allocate(this%tracer_vars(this%num_tracers)) - allocate(this%waterflux_tracer_inst(this%num_tracers)) - allocate(this%waterstate_tracer_inst(this%num_tracers)) - allocate(this%waterdiagnostic_tracer_inst(this%num_tracers)) - allocate(this%waterbalance_tracer_inst(this%num_tracers)) - allocate(this%waterlnd2atm_tracer_inst(this%num_tracers)) - allocate(this%wateratm2lnd_tracer_inst(this%num_tracers)) - end if + do i = this%tracers_beg, this%tracers_end - do i = 1, this%num_tracers + call this%AllocateTracer(i) - call this%tracer_vars(i)%init() + call this%bulk_and_tracers(i)%vars%init() - call this%waterstate_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i), & + call this%bulk_and_tracers(i)%waterstate_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars, & h2osno_input_col = h2osno_col(begc:endc), & watsat_col = watsat_col(begc:endc, 1:), & t_soisno_col = t_soisno_col(begc:endc, -nlevsno+1:) ) - call this%waterdiagnostic_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i)) + call this%bulk_and_tracers(i)%waterdiagnostic_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars) - call this%waterbalance_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i)) + call this%bulk_and_tracers(i)%waterbalance_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars) - call this%waterflux_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i)) + call this%bulk_and_tracers(i)%waterflux_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars) - call this%waterlnd2atm_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i)) + call this%bulk_and_tracers(i)%waterlnd2atm_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars) - call this%wateratm2lnd_tracer_inst(i)%Init(bounds, & - this%tracer_info(i)%info, & - this%tracer_vars(i)) + call this%bulk_and_tracers(i)%wateratm2lnd_inst%Init(bounds, & + this%bulk_and_tracers(i)%info, & + this%bulk_and_tracers(i)%vars) - call this%tracer_vars(i)%complete_setup() + call this%bulk_and_tracers(i)%vars%complete_setup() end do @@ -442,38 +485,115 @@ subroutine SetupTracerInfo(this) num_tracers = num_tracers + 2 end if - this%num_tracers = num_tracers - if (this%num_tracers > 0) then - allocate(this%tracer_info(this%num_tracers)) - allocate(this%is_isotope(this%num_tracers)) - end if + this%bulk_and_tracers_beg = 0 + this%tracers_beg = 1 + this%bulk_and_tracers_end = num_tracers + this%tracers_end = num_tracers + this%i_bulk = 0 + + allocate(this%bulk_and_tracers(this%bulk_and_tracers_beg:this%bulk_and_tracers_end)) + + allocate(this%bulk_and_tracers(this%i_bulk)%info, source = water_info_bulk_type()) tracer_num = 1 if (enable_bulk_tracer) then - allocate(this%tracer_info(tracer_num)%info, source = water_info_isotope_type('H2OTR',1._r8)) - this%is_isotope(tracer_num) = .true. + allocate(this%bulk_and_tracers(tracer_num)%info, source = water_info_isotope_type( & + tracer_name = 'H2OTR', & + ratio = 1._r8, & + communicated_with_coupler = .false.)) + this%bulk_and_tracers(tracer_num)%is_isotope = .true. this%bulk_tracer_index = tracer_num tracer_num = tracer_num + 1 end if if (this%params%enable_isotopes) then - allocate(this%tracer_info(tracer_num)%info, source = water_info_isotope_type('HDO',0.9_r8)) - this%is_isotope(tracer_num) = .true. + allocate(this%bulk_and_tracers(tracer_num)%info, source = water_info_isotope_type( & + tracer_name = 'HDO', & + ratio = 0.9_r8, & + communicated_with_coupler = .false.)) + this%bulk_and_tracers(tracer_num)%is_isotope = .true. tracer_num = tracer_num + 1 - allocate(this%tracer_info(tracer_num)%info, source = water_info_isotope_type('H218O',0.5_r8)) - this%is_isotope(tracer_num) = .true. + allocate(this%bulk_and_tracers(tracer_num)%info, source = water_info_isotope_type( & + tracer_name = 'H218O', & + ratio = 0.5_r8, & + communicated_with_coupler = .false.)) + this%bulk_and_tracers(tracer_num)%is_isotope = .true. tracer_num = tracer_num + 1 end if - if (tracer_num - 1 /= this%num_tracers) then + if (tracer_num - 1 /= num_tracers) then write(iulog,*) subname//' ERROR: tracer_num discrepancy' - write(iulog,*) 'num_tracers = ', this%num_tracers + write(iulog,*) 'num_tracers = ', num_tracers write(iulog,*) 'but added ', tracer_num - 1, ' tracers' call endrun(msg='tracer_num discrepancy '//errMsg(sourcefile, __LINE__)) end if end subroutine SetupTracerInfo + !----------------------------------------------------------------------- + subroutine AllocateBulk(this) + ! + ! !DESCRIPTION: + ! Allocate each of the bulk objects + ! + ! !ARGUMENTS: + class(water_type), intent(inout) :: this + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'AllocateBulk' + !----------------------------------------------------------------------- + + associate( & + i_bulk => this%i_bulk & + ) + + allocate(this%waterfluxbulk_inst) + this%bulk_and_tracers(i_bulk)%waterflux_inst => this%waterfluxbulk_inst + + allocate(this%waterstatebulk_inst) + this%bulk_and_tracers(i_bulk)%waterstate_inst => this%waterstatebulk_inst + + allocate(this%waterdiagnosticbulk_inst) + this%bulk_and_tracers(i_bulk)%waterdiagnostic_inst => this%waterdiagnosticbulk_inst + + allocate(this%waterbalancebulk_inst) + this%bulk_and_tracers(i_bulk)%waterbalance_inst => this%waterbalancebulk_inst + + allocate(this%waterlnd2atmbulk_inst) + this%bulk_and_tracers(i_bulk)%waterlnd2atm_inst => this%waterlnd2atmbulk_inst + + allocate(this%wateratm2lndbulk_inst) + this%bulk_and_tracers(i_bulk)%wateratm2lnd_inst => this%wateratm2lndbulk_inst + + end associate + + end subroutine AllocateBulk + + !----------------------------------------------------------------------- + subroutine AllocateTracer(this, i) + ! + ! !DESCRIPTION: + ! Allocate each of the tracer objects for tracer i + ! + ! !ARGUMENTS: + class(water_type), intent(inout) :: this + integer, intent(in) :: i ! tracer number + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'AllocateTracer' + !----------------------------------------------------------------------- + + allocate(waterflux_type :: this%bulk_and_tracers(i)%waterflux_inst) + allocate(waterstate_type :: this%bulk_and_tracers(i)%waterstate_inst) + allocate(waterdiagnostic_type :: this%bulk_and_tracers(i)%waterdiagnostic_inst) + allocate(waterbalance_type :: this%bulk_and_tracers(i)%waterbalance_inst) + allocate(waterlnd2atm_type :: this%bulk_and_tracers(i)%waterlnd2atm_inst) + allocate(wateratm2lnd_type :: this%bulk_and_tracers(i)%wateratm2lnd_inst) + + end subroutine AllocateTracer + !----------------------------------------------------------------------- subroutine InitAccBuffer(this, bounds) ! @@ -568,16 +688,16 @@ subroutine Restart(this, bounds, ncid, flag, & call this%wateratm2lndbulk_inst%restartBulk (bounds, ncid, flag=flag) - do i = 1, this%num_tracers + do i = this%tracers_beg, this%tracers_end - call this%waterflux_tracer_inst(i)%Restart(bounds, ncid, flag=flag) + call this%bulk_and_tracers(i)%waterflux_inst%Restart(bounds, ncid, flag=flag) - call this%waterstate_tracer_inst(i)%Restart(bounds, ncid, flag=flag, & + call this%bulk_and_tracers(i)%waterstate_inst%Restart(bounds, ncid, flag=flag, & watsat_col=watsat_col(bounds%begc:bounds%endc,:)) - call this%waterdiagnostic_tracer_inst(i)%Restart(bounds, ncid, flag=flag) + call this%bulk_and_tracers(i)%waterdiagnostic_inst%Restart(bounds, ncid, flag=flag) - call this%wateratm2lnd_tracer_inst(i)%Restart(bounds, ncid, flag=flag) + call this%bulk_and_tracers(i)%wateratm2lnd_inst%Restart(bounds, ncid, flag=flag) end do @@ -589,9 +709,7 @@ function IsIsotope(this, i) ! !DESCRIPTION: ! Returns true if tracer i is an isotope ! - ! i must be <= this%num_tracers - ! - ! !USES: + ! i must be >= this%tracers_beg and <= this%tracers_end ! ! !ARGUMENTS: logical :: IsIsotope ! function result @@ -603,9 +721,10 @@ function IsIsotope(this, i) character(len=*), parameter :: subname = 'IsIsotope' !----------------------------------------------------------------------- - SHR_ASSERT(i <= this%num_tracers, errMsg(sourcefile, __LINE__)) + SHR_ASSERT(i >= this%tracers_beg, errMsg(sourcefile, __LINE__)) + SHR_ASSERT(i <= this%tracers_end, errMsg(sourcefile, __LINE__)) - IsIsotope = this%is_isotope(i) + IsIsotope = this%bulk_and_tracers(i)%is_isotope end function IsIsotope @@ -617,7 +736,8 @@ subroutine GetIsotopeInfo(this, i, isotope_info) ! ! This provides a mechanism for passing the isotope info to subroutines that need it. ! - ! i must be <= this%num_tracers, and this%IsIsotope(i) must be true + ! i must be >= this%tracers_beg and <= this%tracers_end, and this%IsIsotope(i) must be + ! true ! ! Assumes that the 'isotope_info' pointer is not currently allocated. (Otherwise this ! will result in a memory leak. It is okay for the isotope_info pointer to be @@ -634,9 +754,10 @@ subroutine GetIsotopeInfo(this, i, isotope_info) character(len=*), parameter :: subname = 'GetIsotopeInfo' !----------------------------------------------------------------------- - SHR_ASSERT(i <= this%num_tracers, errMsg(sourcefile, __LINE__)) + SHR_ASSERT(i >= this%tracers_beg, errMsg(sourcefile, __LINE__)) + SHR_ASSERT(i <= this%tracers_end, errMsg(sourcefile, __LINE__)) - select type(info => this%tracer_info(i)%info) + select type(info => this%bulk_and_tracers(i)%info) type is(water_info_isotope_type) isotope_info => info class default @@ -719,17 +840,22 @@ subroutine TracerConsistencyCheck(this, bounds, caller_location) call endrun(msg=errMsg(sourcefile, __LINE__)) end if - num_vars = this%tracer_vars(i)%get_num_vars() - SHR_ASSERT(num_vars == this%bulk_vars%get_num_vars(), errMsg(sourcefile, __LINE__)) + associate( & + tracer_vars => this%bulk_and_tracers(i)%vars, & + bulk_vars => this%bulk_and_tracers(this%i_bulk)%vars & + ) + + num_vars = tracer_vars%get_num_vars() + SHR_ASSERT(num_vars == bulk_vars%get_num_vars(), errMsg(sourcefile, __LINE__)) do var_num = 1, num_vars - name = this%tracer_vars(i)%get_description(var_num) - SHR_ASSERT(name == this%bulk_vars%get_description(var_num), errMsg(sourcefile, __LINE__)) + name = tracer_vars%get_description(var_num) + SHR_ASSERT(name == bulk_vars%get_description(var_num), errMsg(sourcefile, __LINE__)) - call this%tracer_vars(i)%get_bounds(var_num, bounds, begi, endi) + call tracer_vars%get_bounds(var_num, bounds, begi, endi) - call this%bulk_vars%get_data(var_num, bulk) - call this%tracer_vars(i)%get_data(var_num, tracer) + call bulk_vars%get_data(var_num, bulk) + call tracer_vars%get_data(var_num, tracer) call CompareBulkToTracer(begi, endi, & bulk = bulk(begi:endi), & @@ -739,6 +865,8 @@ subroutine TracerConsistencyCheck(this, bounds, caller_location) end do + end associate + end subroutine TracerConsistencyCheck end module WaterType diff --git a/src/biogeophys/Wateratm2lndBulkType.F90 b/src/biogeophys/Wateratm2lndBulkType.F90 index 758e65e454..573c7d4193 100644 --- a/src/biogeophys/Wateratm2lndBulkType.F90 +++ b/src/biogeophys/Wateratm2lndBulkType.F90 @@ -474,7 +474,6 @@ subroutine Clean(this) ! atm->lnd deallocate(this%forc_rh_grc) - ! atm->lnd not downscaled deallocate(this%forc_q_not_downscaled_grc) deallocate(this%forc_rain_not_downscaled_grc) diff --git a/src/biogeophys/Wateratm2lndType.F90 b/src/biogeophys/Wateratm2lndType.F90 index ad202c0d2f..14d78684bb 100644 --- a/src/biogeophys/Wateratm2lndType.F90 +++ b/src/biogeophys/Wateratm2lndType.F90 @@ -8,13 +8,15 @@ module Wateratm2lndType ! !USES: #include "shr_assert.h" use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use decompMod , only : BOUNDS_SUBGRID_COLUMN, BOUNDS_SUBGRID_GRIDCELL use clm_varctl , only : iulog use clm_varcon , only : spval + use ColumnType , only : col use WaterInfoBaseType, only : water_info_base_type use WaterTracerContainerType, only : water_tracer_container_type - use WaterTracerUtils, only : AllocateVar1d + use WaterTracerUtils, only : AllocateVar1d, CalcTracerFromBulk, CalcTracerFromBulkFixedRatio ! implicit none save @@ -33,10 +35,17 @@ module Wateratm2lndType real(r8), pointer :: forc_rain_downscaled_col (:) ! downscaled atm rain rate [mm/s] real(r8), pointer :: forc_snow_downscaled_col (:) ! downscaled atm snow rate [mm/s] + real(r8), pointer :: rain_to_snow_conversion_col (:) ! amount of rain converted to snow via precipitation repartitioning (mm/s) + real(r8), pointer :: snow_to_rain_conversion_col (:) ! amount of snow converted to rain via precipitation repartitioning (mm/s) + contains procedure, public :: Init procedure, public :: Restart + procedure, public :: IsCommunicatedWithCoupler + procedure, public :: SetNondownscaledTracers + procedure, public :: SetDownscaledTracers + procedure, public :: Clean procedure, private :: InitAllocate procedure, private :: InitHistory procedure, private :: InitCold @@ -114,6 +123,12 @@ subroutine InitAllocate(this, bounds, tracer_vars) container = tracer_vars, & bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN, & ival=ival) + call AllocateVar1d(var = this%rain_to_snow_conversion_col, name = 'rain_to_snow_conversion_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) + call AllocateVar1d(var = this%snow_to_rain_conversion_col, name = 'snow_to_rain_conversion_col', & + container = tracer_vars, & + bounds = bounds, subgrid_level = BOUNDS_SUBGRID_COLUMN) end subroutine InitAllocate @@ -240,4 +255,228 @@ subroutine Restart(this, bounds, ncid, flag) end subroutine Restart + !----------------------------------------------------------------------- + pure function IsCommunicatedWithCoupler(this) result(coupled) + ! + ! !DESCRIPTION: + ! Returns true if this tracer is received from the coupler. Returns false if this + ! tracer is just used internally in CTSM, and should be set to some fixed ratio times + ! the bulk water. + ! + ! !ARGUMENTS: + logical :: coupled ! function result + class(wateratm2lnd_type), intent(in) :: this + !----------------------------------------------------------------------- + + coupled = this%info%is_communicated_with_coupler() + + end function IsCommunicatedWithCoupler + + + !----------------------------------------------------------------------- + subroutine SetNondownscaledTracers(this, bounds, bulk) + ! + ! !DESCRIPTION: + ! Set tracer values for the non-downscaled atm2lnd water quantities from the bulk quantities + ! + ! This should only be called for tracers that are not communicated with the coupler + ! (i.e., for which this%IsCommunicatedWithCoupler() is false). Note that the tracer + ! values are set to a fixed ratio times the bulk (because we don't have any other + ! information to go on for these fields). + ! + ! !ARGUMENTS: + class(wateratm2lnd_type), intent(inout) :: this + type(bounds_type), intent(in) :: bounds + class(wateratm2lnd_type), intent(in) :: bulk + ! + ! !LOCAL VARIABLES: + real(r8) :: ratio + + character(len=*), parameter :: subname = 'SetNondownscaledTracers' + !----------------------------------------------------------------------- + + associate( & + begg => bounds%begg, & + endg => bounds%endg & + ) + + ratio = this%info%get_ratio() + + call CalcTracerFromBulkFixedRatio( & + bulk = bulk%forc_q_not_downscaled_grc(begg:endg), & + ratio = ratio, & + tracer = this%forc_q_not_downscaled_grc(begg:endg)) + + call CalcTracerFromBulkFixedRatio( & + bulk = bulk%forc_rain_not_downscaled_grc(begg:endg), & + ratio = ratio, & + tracer = this%forc_rain_not_downscaled_grc(begg:endg)) + + call CalcTracerFromBulkFixedRatio( & + bulk = bulk%forc_snow_not_downscaled_grc(begg:endg), & + ratio = ratio, & + tracer = this%forc_snow_not_downscaled_grc(begg:endg)) + + call CalcTracerFromBulkFixedRatio( & + bulk = bulk%forc_flood_grc(begg:endg), & + ratio = ratio, & + tracer = this%forc_flood_grc(begg:endg)) + + end associate + + end subroutine SetNondownscaledTracers + + !----------------------------------------------------------------------- + subroutine SetDownscaledTracers(this, bounds, num_allc, filter_allc, & + bulk) + ! + ! !DESCRIPTION: + ! Set tracer values for the downscaled atm2lnd water quantities from the bulk quantities + ! + ! !ARGUMENTS: + class(wateratm2lnd_type) , intent(inout) :: this + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_allc ! number of column points in filter_allc + integer , intent(in) :: filter_allc(:) ! column filter for all points + class(wateratm2lnd_type) , intent(in) :: bulk + ! + ! !LOCAL VARIABLES: + integer :: fc, c, g + + character(len=*), parameter :: subname = 'SetDownscaledTracers' + !----------------------------------------------------------------------- + + associate( & + begg => bounds%begg, & + endg => bounds%endg, & + begc => bounds%begc, & + endc => bounds%endc & + ) + + call SetOneDownscaledTracer( & + bulk_not_downscaled = bulk%forc_q_not_downscaled_grc(begg:endg), & + tracer_not_downscaled = this%forc_q_not_downscaled_grc(begg:endg), & + bulk_downscaled = bulk%forc_q_downscaled_col(begc:endc), & + tracer_downscaled = this%forc_q_downscaled_col(begc:endc)) + + call SetOneDownscaledTracer( & + bulk_not_downscaled = bulk%forc_rain_not_downscaled_grc(begg:endg), & + tracer_not_downscaled = this%forc_rain_not_downscaled_grc(begg:endg), & + bulk_downscaled = bulk%rain_to_snow_conversion_col(begc:endc), & + tracer_downscaled = this%rain_to_snow_conversion_col(begc:endc)) + + call SetOneDownscaledTracer( & + bulk_not_downscaled = bulk%forc_snow_not_downscaled_grc(begg:endg), & + tracer_not_downscaled = this%forc_snow_not_downscaled_grc(begg:endg), & + bulk_downscaled = bulk%snow_to_rain_conversion_col(begc:endc), & + tracer_downscaled = this%snow_to_rain_conversion_col(begc:endc)) + + do fc = 1, num_allc + c = filter_allc(fc) + g = col%gridcell(c) + this%forc_rain_downscaled_col(c) = AdjustPrecipTracer( & + not_downscaled = this%forc_rain_not_downscaled_grc(g), & + addition = this%snow_to_rain_conversion_col(c), & + subtraction = this%rain_to_snow_conversion_col(c)) + this%forc_snow_downscaled_col(c) = AdjustPrecipTracer( & + not_downscaled = this%forc_snow_not_downscaled_grc(g), & + addition = this%rain_to_snow_conversion_col(c), & + subtraction = this%snow_to_rain_conversion_col(c)) + end do + + end associate + + contains + subroutine SetOneDownscaledTracer(bulk_not_downscaled, tracer_not_downscaled, & + bulk_downscaled, tracer_downscaled) + real(r8), intent(in) :: bulk_not_downscaled( bounds%begg: ) + real(r8), intent(in) :: tracer_not_downscaled( bounds%begg: ) + real(r8), intent(in) :: bulk_downscaled( bounds%begc: ) + real(r8), intent(inout) :: tracer_downscaled( bounds%begc: ) + + integer :: fc, c, g + real(r8) :: bulk_not_downscaled_col(bounds%begc:bounds%endc) + real(r8) :: tracer_not_downscaled_col(bounds%begc:bounds%endc) + + SHR_ASSERT_ALL((ubound(bulk_not_downscaled) == [bounds%endg]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(tracer_not_downscaled) == [bounds%endg]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(bulk_downscaled) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + SHR_ASSERT_ALL((ubound(tracer_downscaled) == [bounds%endc]), errMsg(sourcefile, __LINE__)) + + associate( & + begc => bounds%begc, & + endc => bounds%endc & + ) + + do fc = 1, num_allc + c = filter_allc(fc) + g = col%gridcell(c) + ! Note that this copying of bulk_not_downscaled to bulk_not_downscaled_col will + ! be repeated for every tracer. At some point we might want to optimize this so + ! that it's just done once and shared for all tracers (probably by doing this + ! copy outside of the loop over tracers that calls SetDownscaledTracers). + bulk_not_downscaled_col(c) = bulk_not_downscaled(g) + tracer_not_downscaled_col(c) = tracer_not_downscaled(g) + end do + + call CalcTracerFromBulk( & + lb = begc, & + num_pts = num_allc, & + filter_pts = filter_allc, & + bulk_source = bulk_not_downscaled_col(begc:endc), & + bulk_val = bulk_downscaled(begc:endc), & + tracer_source = tracer_not_downscaled_col(begc:endc), & + tracer_val = tracer_downscaled(begc:endc)) + + end associate + + end subroutine SetOneDownscaledTracer + + pure function AdjustPrecipTracer(not_downscaled, addition, subtraction) result(downscaled) + real(r8) :: downscaled + real(r8), intent(in) :: not_downscaled + real(r8), intent(in) :: addition + real(r8), intent(in) :: subtraction + + real(r8), parameter :: tolerance = 1.e-13_r8 + + downscaled = not_downscaled + addition - subtraction + if (not_downscaled /= 0._r8) then + if (abs(downscaled / not_downscaled) < tolerance) then + ! Roundoff correction: If it seems that the downscaled quantity is supposed + ! to go to exactly 0, then make sure it is indeed exactly 0 rather than + ! roundoff-level different from 0. + downscaled = 0._r8 + end if + end if + end function AdjustPrecipTracer + + end subroutine SetDownscaledTracers + + !----------------------------------------------------------------------- + subroutine Clean(this) + ! + ! !DESCRIPTION: + ! Deallocate memory associated with this instance + ! + ! !ARGUMENTS: + class(wateratm2lnd_type), intent(inout) :: this + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'Clean' + !----------------------------------------------------------------------- + + deallocate(this%forc_q_not_downscaled_grc) + deallocate(this%forc_rain_not_downscaled_grc) + deallocate(this%forc_snow_not_downscaled_grc) + deallocate(this%forc_q_downscaled_col) + deallocate(this%forc_flood_grc) + deallocate(this%forc_rain_downscaled_col) + deallocate(this%forc_snow_downscaled_col) + deallocate(this%rain_to_snow_conversion_col) + deallocate(this%snow_to_rain_conversion_col) + + end subroutine Clean + end module Wateratm2lndType diff --git a/src/biogeophys/test/CMakeLists.txt b/src/biogeophys/test/CMakeLists.txt index 7a5df8008d..ac43260605 100644 --- a/src/biogeophys/test/CMakeLists.txt +++ b/src/biogeophys/test/CMakeLists.txt @@ -4,6 +4,7 @@ add_subdirectory(HumanStress_test) add_subdirectory(SnowHydrology_test) add_subdirectory(Balance_test) add_subdirectory(TotalWaterAndHeat_test) +add_subdirectory(Wateratm2lnd_test) add_subdirectory(WaterTracerContainerType_test) add_subdirectory(WaterTracerUtils_test) -add_subdirectory(WaterType_test) \ No newline at end of file +add_subdirectory(WaterType_test) diff --git a/src/biogeophys/test/WaterTracerUtils_test/CMakeLists.txt b/src/biogeophys/test/WaterTracerUtils_test/CMakeLists.txt index 586e787afa..ed16d6200f 100644 --- a/src/biogeophys/test/WaterTracerUtils_test/CMakeLists.txt +++ b/src/biogeophys/test/WaterTracerUtils_test/CMakeLists.txt @@ -1,4 +1,6 @@ set (pfunit_sources + test_calc_tracer_from_bulk_fixed_ratio.pf + test_calc_tracer_from_bulk.pf test_compare_bulk_to_tracer.pf) # extra sources used for this test, which are not .pf files diff --git a/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk.pf b/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk.pf new file mode 100644 index 0000000000..1c5218705c --- /dev/null +++ b/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk.pf @@ -0,0 +1,174 @@ +module test_calc_tracer_from_bulk + + ! Tests of WaterTracerUtils: CalcTracerFromBulk + + use pfunit_mod + use WaterTracerUtils + use shr_kind_mod , only : r8 => shr_kind_r8 + use unittestFilterBuilderMod, only : filter_from_range + use unittestUtils, only : endrun_msg + + implicit none + + @TestCase + type, extends(TestCase) :: TestCalc + contains + procedure :: setUp + procedure :: tearDown + procedure :: doCall3pt + end type TestCalc + + ! Most / all of the tests in this module are 3-point tests, where point 2 is the point + ! of interest. Points 1 and 3 are there to make sure that whatever happens to point 2 + ! doesn't affect other points (e.g., to catch mistakes like whole-array + ! assignment). These are the inputs and expected outputs for those other points. + real(r8), parameter :: bulk_val_other = 10._r8 + real(r8), parameter :: bulk_source_other = 1000._r8 + real(r8), parameter :: tracer_source_other = 500._r8 + real(r8), parameter :: tracer_val_other = 5._r8 + + real(r8), parameter :: tol = 1.e-13_r8 + +contains + + subroutine setUp(this) + class(TestCalc), intent(inout) :: this + end subroutine setUp + + subroutine tearDown(this) + class(TestCalc), intent(inout) :: this + end subroutine tearDown + + subroutine doCall3pt(this, bulk_source, bulk_val, tracer_source, tracer_vals, & + num_pts, filter_pts) + ! Wraps a call to CalcTracerFromBulk with 3 points. Point 2 is the point of interest, + ! and its inputs and outputs are given by the arguments to this routine. Points 1 and + ! 3 have hard-coded inputs. Tests should do assertions on all 3 points, to ensure + ! that (for example) point 1 isn't accidentally overwritten when setting point 2. + class(TestCalc), intent(inout) :: this + real(r8), intent(in) :: bulk_source ! bulk source value in point 2 + real(r8), intent(in) :: bulk_val ! bulk value in point 2 + real(r8), intent(in) :: tracer_source ! tracer source value in point 2 + real(r8), intent(out) :: tracer_vals(17:) ! tracer values in all 3 points + + ! If provided, num_pts and filter_pts give the filter to run over. If not provided, + ! this routine will build a filter that runs over all 3 points. + integer, intent(in), optional :: num_pts + integer, intent(in), optional :: filter_pts(:) + + integer :: l_num_pts + integer, allocatable :: l_filter_pts(:) + + ! Use an arbitrary lower bound that isn't 1, in order to make sure the routine + ! properly handles non-1 lower bounds. Note that this needs to agree with the lower + ! bound of tracer_vals. + real(r8) :: bulk_sources(17:19) + real(r8) :: bulk_vals(17:19) + real(r8) :: tracer_sources(17:19) + + @assertEqual(3, size(tracer_vals)) + if (present(num_pts) .and. present(filter_pts)) then + l_num_pts = num_pts + l_filter_pts = filter_pts + else + @assertFalse(present(num_pts), message = "Must provide both num_pts and filter_pts or neither") + @assertFalse(present(filter_pts), message = "Must provide both num_pts and filter_pts or neither") + call filter_from_range(1, 3, l_num_pts, l_filter_pts) + end if + + bulk_sources(:) = [bulk_source_other, bulk_source, bulk_source_other] + bulk_vals(:) = [bulk_val_other, bulk_val, bulk_val_other] + tracer_sources(:) = [tracer_source_other, tracer_source, tracer_source_other] + call CalcTracerFromBulk( & + lb = 1, & + num_pts = l_num_pts, & + filter_pts = l_filter_pts, & + bulk_source = bulk_sources, & + bulk_val = bulk_vals, & + tracer_source = tracer_sources, & + tracer_val = tracer_vals) + + end subroutine doCall3pt + + @Test + subroutine basic(this) + class(TestCalc), intent(inout) :: this + real(r8) :: tracer_vals(3) + + call this%doCall3pt( & + bulk_source = 300._r8, & + bulk_val = 30._r8, & + tracer_source = 60._r8, & + tracer_vals = tracer_vals) + + @assertEqual([tracer_val_other, 6._r8, tracer_val_other], tracer_vals) + end subroutine basic + + @Test + subroutine outsideFilter_valsUnchanged(this) + class(TestCalc), intent(inout) :: this + real(r8) :: tracer_vals(3) + integer :: num_pts + integer, allocatable :: filter_pts(:) + + tracer_vals(:) = [-1._r8, -2._r8, -3._r8] + call filter_from_range(1, 1, num_pts, filter_pts) + + call this%doCall3pt( & + bulk_source = 300._r8, & + bulk_val = 30._r8, & + tracer_source = 60._r8, & + tracer_vals = tracer_vals, & + num_pts = num_pts, & + filter_pts = filter_pts) + + @assertEqual([tracer_val_other, -2._r8, -3._r8], tracer_vals) + end subroutine outsideFilter_valsUnchanged + + @Test + subroutine bulkSource0_tracerSource0_bulkVal0_yieldsTracerVal0(this) + class(TestCalc), intent(inout) :: this + real(r8) :: tracer_vals(3) + + call this%doCall3pt( & + bulk_source = 0._r8, & + bulk_val = 0._r8, & + tracer_source = 0._r8, & + tracer_vals = tracer_vals) + + @assertEqual([tracer_val_other, 0._r8, tracer_val_other], tracer_vals) + end subroutine bulkSource0_tracerSource0_bulkVal0_yieldsTracerVal0 + + @Test + subroutine bulkSource0_tracerSource0_bulkValNon0_aborts(this) + class(TestCalc), intent(inout) :: this + real(r8) :: tracer_vals(3) + character(len=:), allocatable :: expected_msg + + call this%doCall3pt( & + bulk_source = 0._r8, & + bulk_val = 1._r8, & + tracer_source = 0._r8, & + tracer_vals = tracer_vals) + + expected_msg = endrun_msg('CalcTracerFromBulk: Non-zero bulk val despite zero bulk source') + @assertExceptionRaised(expected_msg) + end subroutine bulkSource0_tracerSource0_bulkValNon0_aborts + + @Test + subroutine bulkSource0_tracerSourceNon0_bulkVal0_aborts(this) + class(TestCalc), intent(inout) :: this + real(r8) :: tracer_vals(3) + character(len=:), allocatable :: expected_msg + + call this%doCall3pt( & + bulk_source = 0._r8, & + bulk_val = 0._r8, & + tracer_source = 1._r8, & + tracer_vals = tracer_vals) + + expected_msg = endrun_msg('CalcTracerFromBulk: Non-zero tracer source despite zero bulk source') + @assertExceptionRaised(expected_msg) + end subroutine bulkSource0_tracerSourceNon0_bulkVal0_aborts + +end module test_calc_tracer_from_bulk diff --git a/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk_fixed_ratio.pf b/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk_fixed_ratio.pf new file mode 100644 index 0000000000..5bde146ac3 --- /dev/null +++ b/src/biogeophys/test/WaterTracerUtils_test/test_calc_tracer_from_bulk_fixed_ratio.pf @@ -0,0 +1,43 @@ +module test_calc_tracer_from_bulk_fixed_ratio + + ! Tests of WaterTracerUtils: CalcTracerFromBulkFixedRatio + + use pfunit_mod + use WaterTracerUtils, only : CalcTracerFromBulkFixedRatio + use shr_kind_mod , only : r8 => shr_kind_r8 + + implicit none + + @TestCase + type, extends(TestCase) :: TestCalcFixedRatio + contains + procedure :: setUp + procedure :: tearDown + end type TestCalcFixedRatio + + real(r8), parameter :: tol = 1.e-13_r8 + +contains + + subroutine setUp(this) + class(TestCalcFixedRatio), intent(inout) :: this + end subroutine setUp + + subroutine tearDown(this) + class(TestCalcFixedRatio), intent(inout) :: this + end subroutine tearDown + + @Test + subroutine basic(this) + class(TestCalcFixedRatio), intent(inout) :: this + real(r8) :: tracer(3) + + call CalcTracerFromBulkFixedRatio( & + bulk = [2._r8, 6._r8, 4._r8], & + ratio = 0.5_r8, & + tracer = tracer) + + @assertEqual([1._r8, 3._r8, 2._r8], tracer, tolerance=tol) + end subroutine basic + +end module test_calc_tracer_from_bulk_fixed_ratio diff --git a/src/biogeophys/test/WaterType_test/test_water_type.pf b/src/biogeophys/test/WaterType_test/test_water_type.pf index 7adbda4278..dded573efe 100644 --- a/src/biogeophys/test/WaterType_test/test_water_type.pf +++ b/src/biogeophys/test/WaterType_test/test_water_type.pf @@ -7,16 +7,15 @@ module test_water_type use shr_kind_mod , only : r8 => shr_kind_r8 use unittestSubgridMod, only : bounds, unittest_subgrid_teardown use unittestSimpleSubgridSetupsMod, only : setup_single_veg_patch - use unittestArrayMod, only : col_array use unittestUtils, only : endrun_msg - use clm_varpar, only : nlevsoi, nlevgrnd, nlevsno - use ColumnType, only : col + use unittestWaterTypeFactory, only : unittest_water_type_factory_type implicit none @TestCase type, extends(TestCase) :: TestWaterType type(water_type) :: water_inst + type(unittest_water_type_factory_type) :: factory contains procedure :: setUp procedure :: tearDown @@ -30,19 +29,17 @@ contains subroutine setUp(this) class(TestWaterType), intent(inout) :: this - ! Arbitrarily set nlevs (needed for allocating multi-level variables) - nlevsoi = 3 - nlevgrnd = nlevsoi + 1 - nlevsno = 3 + call this%factory%init() + call this%factory%setup_before_subgrid( & + my_nlevsoi = 3, & + nlevgrnd_additional = 1, & + my_nlevsno = 3) end subroutine setUp subroutine tearDown(this) class(TestWaterType), intent(inout) :: this - ! Ideally this would call water_inst%Clean. But we don't yet have a clean method in - ! water_type or most of the types it contains. So for now we have a small memory leak - ! in each test. - + call this%factory%teardown(this%water_inst) call unittest_subgrid_teardown() end subroutine tearDown @@ -51,29 +48,12 @@ contains ! this%water_inst class(TestWaterType), intent(inout) :: this - type(water_params_type) :: params - real(r8), allocatable :: watsat_col(:,:) - real(r8), allocatable :: t_soisno_col(:,:) - call setup_single_veg_patch(pft_type=1) - - col%snl(:) = 0 - col%dz(:,:) = 1._r8 - - params = water_params_type( & + call this%factory%setup_after_subgrid(snl = 0, dz = 1._r8) + call this%factory%create_water_type(this%water_inst, & enable_consistency_checks = .true., & enable_isotopes = .false.) - allocate(watsat_col(bounds%begc:bounds%endc, nlevgrnd)) - watsat_col(:,:) = 0._r8 - allocate(t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) - t_soisno_col(:,:) = 275._r8 - call this%water_inst%InitForTesting(bounds, params, & - h2osno_col = col_array(0._r8), & - snow_depth_col = col_array(0._r8), & - watsat_col = watsat_col, & - t_soisno_col = t_soisno_col) - end subroutine init_water_inst_1pt ! ------------------------------------------------------------------------ @@ -129,8 +109,8 @@ contains bulk_tracer = this%water_inst%GetBulkTracerIndex() @assertGreaterThan(bulk_tracer, 0) - this%water_inst%waterstate_tracer_inst(bulk_tracer)%h2osfc_col(bounds%begc) = & - this%water_inst%waterstate_tracer_inst(bulk_tracer)%h2osfc_col(bounds%begc) + 1._r8 + this%water_inst%bulk_and_tracers(bulk_tracer)%waterstate_inst%h2osfc_col(bounds%begc) = & + this%water_inst%bulk_and_tracers(bulk_tracer)%waterstate_inst%h2osfc_col(bounds%begc) + 1._r8 @assertTrue(this%water_inst%DoConsistencyCheck()) @@ -161,7 +141,7 @@ contains ! We set both the bulk and tracer values explicitly, because we can't rely on these ! being set to reasonable values in all layers in initialization this%water_inst%waterstatebulk_inst%h2osoi_liq_col(bounds%begc, lev) = 1._r8 - this%water_inst%waterstate_tracer_inst(bulk_tracer)%h2osoi_liq_col(bounds%begc, lev) = 2._r8 + this%water_inst%bulk_and_tracers(bulk_tracer)%waterstate_inst%h2osoi_liq_col(bounds%begc, lev) = 2._r8 @assertTrue(this%water_inst%DoConsistencyCheck()) diff --git a/src/biogeophys/test/Wateratm2lnd_test/CMakeLists.txt b/src/biogeophys/test/Wateratm2lnd_test/CMakeLists.txt new file mode 100644 index 0000000000..35e12aea3e --- /dev/null +++ b/src/biogeophys/test/Wateratm2lnd_test/CMakeLists.txt @@ -0,0 +1,12 @@ +set (pfunit_sources + test_set_tracers.pf) + +# extra sources used for this test, which are not .pf files +# (currently none) +set (extra_sources + ) + +create_pFUnit_test(water_atm2lnd test_water_atm2lnd_exe + "${pfunit_sources}" "${extra_sources}") + +target_link_libraries(test_water_atm2lnd_exe clm csm_share esmf_wrf_timemgr) \ No newline at end of file diff --git a/src/biogeophys/test/Wateratm2lnd_test/test_set_tracers.pf b/src/biogeophys/test/Wateratm2lnd_test/test_set_tracers.pf new file mode 100644 index 0000000000..54d7491b71 --- /dev/null +++ b/src/biogeophys/test/Wateratm2lnd_test/test_set_tracers.pf @@ -0,0 +1,190 @@ +module test_set_tracers + + ! Tests of Wateratm2lndType: routines that set tracers + + use pfunit_mod + use Wateratm2lndType + use WaterInfoBulkType, only : water_info_bulk_type + use WaterInfoTracerType, only : water_info_tracer_type + use WaterTracerContainerType, only : water_tracer_container_type + use shr_kind_mod , only : r8 => shr_kind_r8 + use unittestSimpleSubgridSetupsMod, only : setup_single_veg_patch + use unittestSubgridMod, only : bounds, unittest_subgrid_teardown + use unittestFilterBuilderMod, only : filter_from_range + + implicit none + + @TestCase + type, extends(TestCase) :: TestSetTracers + type(wateratm2lnd_type) :: wateratm2lnd_bulk + type(wateratm2lnd_type) :: wateratm2lnd_tracer + integer :: num_allc + integer, allocatable :: filter_allc(:) + contains + procedure :: setUp + procedure :: tearDown + procedure :: setInputs1Col + end type TestSetTracers + + real(r8), parameter :: tol = 1.e-13_r8 + +contains + + ! ======================================================================== + ! Helper routines + ! ======================================================================== + + subroutine setUp(this) + class(TestSetTracers), intent(inout) :: this + type(water_tracer_container_type) :: my_vars + + call my_vars%init() + + call setup_single_veg_patch(pft_type = 1) + + call this%wateratm2lnd_bulk%Init( & + bounds = bounds, & + info = water_info_bulk_type(), & + tracer_vars = my_vars) + + call this%wateratm2lnd_tracer%Init( & + bounds = bounds, & + + info = water_info_tracer_type(& + tracer_name = 'foo', & + ratio = 0.456_r8, & + communicated_with_coupler = .false.), & + + tracer_vars = my_vars) + + call filter_from_range(bounds%begc, bounds%endc, this%num_allc, this%filter_allc) + end subroutine setUp + + subroutine tearDown(this) + class(TestSetTracers), intent(inout) :: this + + call this%wateratm2lnd_bulk%Clean() + call this%wateratm2lnd_tracer%Clean() + call unittest_subgrid_teardown() + end subroutine tearDown + + subroutine setInputs1Col(this, & + bulk_rain_not_downscaled, bulk_snow_not_downscaled, & + tracer_rain_not_downscaled, tracer_snow_not_downscaled, & + bulk_rain_to_snow, bulk_snow_to_rain) + class(TestSetTracers), intent(inout) :: this + real(r8), intent(in) :: bulk_rain_not_downscaled + real(r8), intent(in) :: bulk_snow_not_downscaled + real(r8), intent(in) :: tracer_rain_not_downscaled + real(r8), intent(in) :: tracer_snow_not_downscaled + real(r8), intent(in) :: bulk_rain_to_snow + real(r8), intent(in) :: bulk_snow_to_rain + + this%wateratm2lnd_bulk%forc_rain_not_downscaled_grc(bounds%begg) = bulk_rain_not_downscaled + this%wateratm2lnd_bulk%forc_snow_not_downscaled_grc(bounds%begg) = bulk_snow_not_downscaled + this%wateratm2lnd_tracer%forc_rain_not_downscaled_grc(bounds%begg) = tracer_rain_not_downscaled + this%wateratm2lnd_tracer%forc_snow_not_downscaled_grc(bounds%begg) = tracer_snow_not_downscaled + this%wateratm2lnd_bulk%rain_to_snow_conversion_col(bounds%begc) = bulk_rain_to_snow + this%wateratm2lnd_bulk%snow_to_rain_conversion_col(bounds%begc) = bulk_snow_to_rain + + ! These aren't used for the unit tests here, but need to be set to something reasonable + this%wateratm2lnd_bulk%forc_q_not_downscaled_grc(bounds%begg) = 0._r8 + this%wateratm2lnd_tracer%forc_q_not_downscaled_grc(bounds%begg) = 0._r8 + this%wateratm2lnd_bulk%forc_q_downscaled_col(bounds%begc) = 0._r8 + + end subroutine setInputs1Col + + ! ======================================================================== + ! Begin tests + ! ======================================================================== + + @Test + subroutine downscaled_precip_noConversion(this) + ! Test downscaled rain and snow when there is no rain-snow conversion + class(TestSetTracers), intent(inout) :: this + + call this%setInputs1Col( & + bulk_rain_not_downscaled = 10._r8, & + bulk_snow_not_downscaled = 90._r8, & + tracer_rain_not_downscaled = 10._r8, & + tracer_snow_not_downscaled = 9._r8, & + bulk_rain_to_snow = 0._r8, & + bulk_snow_to_rain = 0._r8) + + call this%wateratm2lnd_tracer%SetDownscaledTracers(bounds, & + this%num_allc, this%filter_allc, & + this%wateratm2lnd_bulk) + + @assertEqual(10._r8, this%wateratm2lnd_tracer%forc_rain_downscaled_col(bounds%begc)) + @assertEqual(9._r8, this%wateratm2lnd_tracer%forc_snow_downscaled_col(bounds%begc)) + end subroutine downscaled_precip_noConversion + + @Test + subroutine downscaled_precip_snowToRain(this) + ! Test downscaled rain and snow when there is snow-to-rain conversion + class(TestSetTracers), intent(inout) :: this + + call this%setInputs1Col( & + bulk_rain_not_downscaled = 10._r8, & + bulk_snow_not_downscaled = 90._r8, & + tracer_rain_not_downscaled = 10._r8, & + tracer_snow_not_downscaled = 9._r8, & + bulk_rain_to_snow = 0._r8, & + bulk_snow_to_rain = 40._r8) + + call this%wateratm2lnd_tracer%SetDownscaledTracers(bounds, & + this%num_allc, this%filter_allc, & + this%wateratm2lnd_bulk) + + ! A conversion of 40 in the bulk results with 4 for the tracer, since + ! tracer_snow/bulk_snow = 1/10 + @assertEqual(14._r8, this%wateratm2lnd_tracer%forc_rain_downscaled_col(bounds%begc), tolerance=tol) + @assertEqual(5._r8, this%wateratm2lnd_tracer%forc_snow_downscaled_col(bounds%begc), tolerance=tol) + end subroutine downscaled_precip_snowToRain + + @Test + subroutine downscaled_precip_rainToSnow(this) + ! Test downscaled rain and snow when there is rain-to-snow conversion + class(TestSetTracers), intent(inout) :: this + + call this%setInputs1Col( & + bulk_rain_not_downscaled = 90._r8, & + bulk_snow_not_downscaled = 10._r8, & + tracer_rain_not_downscaled = 9._r8, & + tracer_snow_not_downscaled = 10._r8, & + bulk_rain_to_snow = 40._r8, & + bulk_snow_to_rain = 0._r8) + + call this%wateratm2lnd_tracer%SetDownscaledTracers(bounds, & + this%num_allc, this%filter_allc, & + this%wateratm2lnd_bulk) + + ! A conversion of 40 in the bulk results with 4 for the tracer, since + ! tracer_snow/bulk_snow = 1/10 + @assertEqual(5._r8, this%wateratm2lnd_tracer%forc_rain_downscaled_col(bounds%begc), tolerance=tol) + @assertEqual(14._r8, this%wateratm2lnd_tracer%forc_snow_downscaled_col(bounds%begc), tolerance=tol) + end subroutine downscaled_precip_rainToSnow + + @Test + subroutine downscaled_precip_rainToZero(this) + ! Test downscaled rain and snow when all rain is converted to snow + class(TestSetTracers), intent(inout) :: this + + call this%setInputs1Col( & + bulk_rain_not_downscaled = 90._r8, & + bulk_snow_not_downscaled = 10._r8, & + tracer_rain_not_downscaled = 9._r8, & + tracer_snow_not_downscaled = 10._r8, & + bulk_rain_to_snow = 90._r8 - 1.e-14_r8, & + bulk_snow_to_rain = 0._r8) + + call this%wateratm2lnd_tracer%SetDownscaledTracers(bounds, & + this%num_allc, this%filter_allc, & + this%wateratm2lnd_bulk) + + ! tracer rain should end up EXACTLY 0 even if there are small roundoff errors in rain_to_snow + @assertEqual(0._r8, this%wateratm2lnd_tracer%forc_rain_downscaled_col(bounds%begc)) + @assertEqual(19._r8, this%wateratm2lnd_tracer%forc_snow_downscaled_col(bounds%begc), tolerance=tol) + end subroutine downscaled_precip_rainToZero + +end module test_set_tracers diff --git a/src/dyn_subgrid/dynConsBiogeophysMod.F90 b/src/dyn_subgrid/dynConsBiogeophysMod.F90 index 377a51acc3..a143728722 100644 --- a/src/dyn_subgrid/dynConsBiogeophysMod.F90 +++ b/src/dyn_subgrid/dynConsBiogeophysMod.F90 @@ -9,37 +9,41 @@ module dynConsBiogeophysMod ! cover. ! ! !USES: - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_log_mod , only : errMsg => shr_log_errMsg - use decompMod , only : bounds_type - use UrbanParamsType , only : urbanparams_type - use EnergyFluxType , only : energyflux_type - use SoilHydrologyType , only : soilhydrology_type - use SoilStateType , only : soilstate_type - use TemperatureType , only : temperature_type - use WaterFluxBulkType , only : waterfluxbulk_type - use WaterStateBulkType , only : waterstatebulk_type - use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type - use WaterBalanceType , only : waterbalance_type - use TotalWaterAndHeatMod, only : ComputeLiqIceMassNonLake, ComputeLiqIceMassLake - use TotalWaterAndHeatMod, only : ComputeHeatNonLake, ComputeHeatLake - use TotalWaterAndHeatMod, only : AdjustDeltaHeatForDeltaLiq - use TotalWaterAndHeatMod, only : heat_base_temp - use subgridAveMod , only : p2c, c2g - use LandunitType , only : lun - use ColumnType , only : col - use PatchType , only : patch - use clm_varcon , only : tfrz, cpliq, hfus - use dynSubgridControlMod, only : get_for_testing_zero_dynbal_fluxes + use shr_kind_mod , only : r8 => shr_kind_r8 + use shr_log_mod , only : errMsg => shr_log_errMsg + use decompMod , only : bounds_type + use UrbanParamsType , only : urbanparams_type + use EnergyFluxType , only : energyflux_type + use SoilHydrologyType , only : soilhydrology_type + use SoilStateType , only : soilstate_type + use TemperatureType , only : temperature_type + use WaterFluxType , only : waterflux_type + use WaterStateBulkType , only : waterstatebulk_type + use WaterStateType , only : waterstate_type + use WaterDiagnosticType , only : waterdiagnostic_type + use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type + use WaterBalanceType , only : waterbalance_type + use WaterType , only : water_type + use TotalWaterAndHeatMod , only : ComputeLiqIceMassNonLake, ComputeLiqIceMassLake + use TotalWaterAndHeatMod , only : ComputeHeatNonLake, ComputeHeatLake + use TotalWaterAndHeatMod , only : AdjustDeltaHeatForDeltaLiq + use TotalWaterAndHeatMod , only : heat_base_temp + use subgridAveMod , only : p2c, c2g + use LandunitType , only : lun + use ColumnType , only : col + use PatchType , only : patch + use clm_varcon , only : tfrz, cpliq, hfus + use dynSubgridControlMod , only : get_for_testing_zero_dynbal_fluxes ! ! !PUBLIC MEMBER FUNCTIONS: implicit none private ! - public :: dyn_hwcontent_init ! compute grid-level heat and water content, before land cover change - public :: dyn_hwcontent_final ! compute grid-level heat and water content, after land cover change; also compute dynbal fluxes + public :: dyn_hwcontent_init ! compute grid-level heat and water content, before land cover change + public :: dyn_hwcontent_final ! compute grid-level heat and water content and dynbal fluxes after land cover change ! ! !PRIVATE MEMBER FUNCTIONS + private :: dyn_water_content_final ! compute grid-level water content and dynbal fluxes after landcover change, for a single water tracer or bulk water private :: dyn_water_content ! compute gridcell total liquid and ice water contents private :: dyn_heat_content ! compute gridcell total heat contents @@ -54,12 +58,10 @@ subroutine dyn_hwcontent_init(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst, & - waterfluxbulk_inst, temperature_inst, energyflux_inst) + water_inst, temperature_inst) ! ! !DESCRIPTION: - ! Initialize variables used for dyn_hwcontent, and compute grid cell-level heat - ! and water content before land cover change + ! Compute grid cell-level heat and water content before land cover change ! ! Should be called BEFORE any subgrid weight updates this time step ! @@ -71,32 +73,41 @@ subroutine dyn_hwcontent_init(bounds, & integer , intent(in) :: filter_lakec(:) type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst - type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst + type(water_type) , intent(inout) :: water_inst type(temperature_type) , intent(inout) :: temperature_inst - type(energyflux_type) , intent(inout) :: energyflux_inst ! ! !LOCAL VARIABLES: - integer :: g ! grid cell index + integer :: i !------------------------------------------------------------------------------- - - call dyn_water_content(bounds, & - num_nolakec, filter_nolakec, & - num_lakec, filter_lakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, & - liquid_mass = waterbalancebulk_inst%liq1_grc(bounds%begg:bounds%endg), & - ice_mass = waterbalancebulk_inst%ice1_grc(bounds%begg:bounds%endg)) + + associate( & + begg => bounds%begg, & + endg => bounds%endg & + ) + + do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end + associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) + call dyn_water_content(bounds, & + num_nolakec, filter_nolakec, & + num_lakec, filter_lakec, & + bulk_or_tracer%waterstate_inst, & + bulk_or_tracer%waterdiagnostic_inst, & + liquid_mass = bulk_or_tracer%waterbalance_inst%liq1_grc(begg:endg), & + ice_mass = bulk_or_tracer%waterbalance_inst%ice1_grc(begg:endg)) + end associate + end do call dyn_heat_content( bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, & - heat_grc = temperature_inst%heat1_grc(bounds%begg:bounds%endg), & - liquid_water_temp_grc = temperature_inst%liquid_water_temp1_grc(bounds%begg:bounds%endg)) + temperature_inst, & + water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & + heat_grc = temperature_inst%heat1_grc(begg:endg), & + liquid_water_temp_grc = temperature_inst%liquid_water_temp1_grc(begg:endg)) + + end associate end subroutine dyn_hwcontent_init @@ -105,12 +116,11 @@ subroutine dyn_hwcontent_final(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst, & - waterfluxbulk_inst, temperature_inst, energyflux_inst) + water_inst, & + temperature_inst, energyflux_inst) ! ! !DESCRIPTION: - ! Compute grid cell-level heat and water content after land cover change, and compute - ! the dynbal fluxes + ! Compute grid cell-level heat and water content and dynbal fluxes after land cover change ! ! Should be called AFTER all subgrid weight updates this time step ! @@ -122,97 +132,163 @@ subroutine dyn_hwcontent_final(bounds, & integer , intent(in) :: filter_lakec(:) type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst - type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst + type(water_type) , intent(inout) :: water_inst type(temperature_type) , intent(inout) :: temperature_inst type(energyflux_type) , intent(inout) :: energyflux_inst ! ! !LOCAL VARIABLES: - integer :: begg, endg + integer :: i integer :: g ! grid cell index - real(r8) :: delta_liq(bounds%begg:bounds%endg) ! change in gridcell h2o liq content - real(r8) :: delta_ice(bounds%begg:bounds%endg) ! change in gridcell h2o ice content + real(r8) :: this_delta_liq(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk or one tracer + real(r8) :: delta_liq_bulk(bounds%begg:bounds%endg) ! change in gridcell h2o liq content for bulk water real(r8) :: delta_heat(bounds%begg:bounds%endg) ! change in gridcell heat content !--------------------------------------------------------------------------- - begg = bounds%begg - endg = bounds%endg - - call dyn_water_content(bounds, & - num_nolakec, filter_nolakec, & - num_lakec, filter_lakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, & - liquid_mass = waterbalancebulk_inst%liq2_grc(bounds%begg:bounds%endg), & - ice_mass = waterbalancebulk_inst%ice2_grc(bounds%begg:bounds%endg)) + associate( & + begg => bounds%begg, & + endg => bounds%endg) + + do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end + associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) + call dyn_water_content_final(bounds, & + num_nolakec, filter_nolakec, & + num_lakec, filter_lakec, & + bulk_or_tracer%waterstate_inst, & + bulk_or_tracer%waterdiagnostic_inst, & + bulk_or_tracer%waterbalance_inst, & + bulk_or_tracer%waterflux_inst, & + delta_liq = this_delta_liq(begg:endg)) + if (i == water_inst%i_bulk) then + delta_liq_bulk(begg:endg) = this_delta_liq(begg:endg) + end if + end associate + end do call dyn_heat_content( bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & urbanparams_inst, soilstate_inst, & - temperature_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, & - heat_grc = temperature_inst%heat2_grc(bounds%begg:bounds%endg), & - liquid_water_temp_grc = temperature_inst%liquid_water_temp2_grc(bounds%begg:bounds%endg)) + temperature_inst, & + water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & + heat_grc = temperature_inst%heat2_grc(begg:endg), & + liquid_water_temp_grc = temperature_inst%liquid_water_temp2_grc(begg:endg)) if (get_for_testing_zero_dynbal_fluxes()) then do g = begg, endg - delta_liq(g) = 0._r8 - delta_ice(g) = 0._r8 delta_heat(g) = 0._r8 end do else do g = begg, endg - delta_liq(g) = waterbalancebulk_inst%liq2_grc(g) - waterbalancebulk_inst%liq1_grc(g) - delta_ice(g) = waterbalancebulk_inst%ice2_grc(g) - waterbalancebulk_inst%ice1_grc(g) delta_heat(g) = temperature_inst%heat2_grc(g) - temperature_inst%heat1_grc(g) end do end if call AdjustDeltaHeatForDeltaLiq( & bounds, & - delta_liq = delta_liq(bounds%begg:bounds%endg), & - liquid_water_temp1 = temperature_inst%liquid_water_temp1_grc(bounds%begg:bounds%endg), & - liquid_water_temp2 = temperature_inst%liquid_water_temp2_grc(bounds%begg:bounds%endg), & - delta_heat = delta_heat(bounds%begg:bounds%endg)) - - call waterfluxbulk_inst%qflx_liq_dynbal_dribbler%set_curr_delta(bounds, & - delta_liq(begg:endg)) - call waterfluxbulk_inst%qflx_liq_dynbal_dribbler%get_curr_flux(bounds, & - waterfluxbulk_inst%qflx_liq_dynbal_grc(begg:endg)) - - call waterfluxbulk_inst%qflx_ice_dynbal_dribbler%set_curr_delta(bounds, & - delta_ice(begg:endg)) - call waterfluxbulk_inst%qflx_ice_dynbal_dribbler%get_curr_flux(bounds, & - waterfluxbulk_inst%qflx_ice_dynbal_grc(begg:endg)) + delta_liq = delta_liq_bulk(begg:endg), & + liquid_water_temp1 = temperature_inst%liquid_water_temp1_grc(begg:endg), & + liquid_water_temp2 = temperature_inst%liquid_water_temp2_grc(begg:endg), & + delta_heat = delta_heat(begg:endg)) call energyflux_inst%eflx_dynbal_dribbler%set_curr_delta(bounds, & delta_heat(begg:endg)) call energyflux_inst%eflx_dynbal_dribbler%get_curr_flux(bounds, & energyflux_inst%eflx_dynbal_grc(begg:endg)) + end associate + end subroutine dyn_hwcontent_final + !----------------------------------------------------------------------- + subroutine dyn_water_content_final(bounds, & + num_nolakec, filter_nolakec, & + num_lakec, filter_lakec, & + waterstate_inst, waterdiagnostic_inst, & + waterbalance_inst, waterflux_inst, & + delta_liq) + ! + ! !DESCRIPTION: + ! Compute grid cell-level water content and dynbal fluxes after landcover change, for + ! a single water tracer or bulk water + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_nolakec + integer , intent(in) :: filter_nolakec(:) + integer , intent(in) :: num_lakec + integer , intent(in) :: filter_lakec(:) + class(waterstate_type) , intent(in) :: waterstate_inst + class(waterdiagnostic_type) , intent(in) :: waterdiagnostic_inst + class(waterbalance_type) , intent(inout) :: waterbalance_inst + class(waterflux_type) , intent(inout) :: waterflux_inst + real(r8) , intent(out) :: delta_liq(bounds%begg:) ! change in gridcell h2o liq content + ! + ! !LOCAL VARIABLES: + integer :: g + real(r8) :: delta_ice(bounds%begg:bounds%endg) ! change in gridcell h2o ice content + + character(len=*), parameter :: subname = 'dyn_water_content_final' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL((ubound(delta_liq) == [bounds%endg]), errMsg(sourcefile, __LINE__)) + + associate( & + begg => bounds%begg, & + endg => bounds%endg) + + call dyn_water_content(bounds, & + num_nolakec, filter_nolakec, & + num_lakec, filter_lakec, & + waterstate_inst, waterdiagnostic_inst, & + liquid_mass = waterbalance_inst%liq2_grc(bounds%begg:bounds%endg), & + ice_mass = waterbalance_inst%ice2_grc(bounds%begg:bounds%endg)) + + if (get_for_testing_zero_dynbal_fluxes()) then + do g = begg, endg + delta_liq(g) = 0._r8 + delta_ice(g) = 0._r8 + end do + else + do g = begg, endg + delta_liq(g) = waterbalance_inst%liq2_grc(g) - waterbalance_inst%liq1_grc(g) + delta_ice(g) = waterbalance_inst%ice2_grc(g) - waterbalance_inst%ice1_grc(g) + end do + end if + + call waterflux_inst%qflx_liq_dynbal_dribbler%set_curr_delta(bounds, & + delta_liq(begg:endg)) + call waterflux_inst%qflx_liq_dynbal_dribbler%get_curr_flux(bounds, & + waterflux_inst%qflx_liq_dynbal_grc(begg:endg)) + + call waterflux_inst%qflx_ice_dynbal_dribbler%set_curr_delta(bounds, & + delta_ice(begg:endg)) + call waterflux_inst%qflx_ice_dynbal_dribbler%get_curr_flux(bounds, & + waterflux_inst%qflx_ice_dynbal_grc(begg:endg)) + + end associate + + end subroutine dyn_water_content_final + !----------------------------------------------------------------------- subroutine dyn_water_content(bounds, & num_nolakec, filter_nolakec, & num_lakec, filter_lakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, & + waterstate_inst, waterdiagnostic_inst, & liquid_mass, ice_mass) ! ! !DESCRIPTION: ! Compute gridcell total liquid and ice water contents ! ! !ARGUMENTS: - type(bounds_type) , intent(in) :: bounds - integer , intent(in) :: num_nolakec - integer , intent(in) :: filter_nolakec(:) - integer , intent(in) :: num_lakec - integer , intent(in) :: filter_lakec(:) - type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst - real(r8) , intent(out) :: liquid_mass( bounds%begg: ) ! kg m-2 - real(r8) , intent(out) :: ice_mass( bounds%begg: ) ! kg m-2 + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_nolakec + integer , intent(in) :: filter_nolakec(:) + integer , intent(in) :: num_lakec + integer , intent(in) :: filter_lakec(:) + class(waterstate_type) , intent(in) :: waterstate_inst + class(waterdiagnostic_type) , intent(in) :: waterdiagnostic_inst + real(r8) , intent(out) :: liquid_mass( bounds%begg: ) ! kg m-2 + real(r8) , intent(out) :: ice_mass( bounds%begg: ) ! kg m-2 ! ! !LOCAL VARIABLES: real(r8) :: liquid_mass_col(bounds%begc:bounds%endc) ! kg m-2 @@ -225,12 +301,12 @@ subroutine dyn_water_content(bounds, & SHR_ASSERT_ALL((ubound(ice_mass) == (/bounds%endg/)), errMsg(sourcefile, __LINE__)) call ComputeLiqIceMassNonLake(bounds, num_nolakec, filter_nolakec, & - waterstatebulk_inst, waterdiagnosticbulk_inst, & + waterstate_inst, waterdiagnostic_inst, & liquid_mass_col(bounds%begc:bounds%endc), & ice_mass_col(bounds%begc:bounds%endc)) call ComputeLiqIceMassLake(bounds, num_lakec, filter_lakec, & - waterstatebulk_inst, & + waterstate_inst, & liquid_mass_col(bounds%begc:bounds%endc), & ice_mass_col(bounds%begc:bounds%endc)) diff --git a/src/dyn_subgrid/dynInitColumnsMod.F90 b/src/dyn_subgrid/dynInitColumnsMod.F90 index 62c06c0165..90f6be421f 100644 --- a/src/dyn_subgrid/dynInitColumnsMod.F90 +++ b/src/dyn_subgrid/dynInitColumnsMod.F90 @@ -14,7 +14,7 @@ module dynInitColumnsMod use clm_varctl , only : iulog use clm_varcon , only : namec use TemperatureType , only : temperature_type - use WaterStateBulkType , only : waterstatebulk_type + use WaterType , only : water_type use SoilHydrologyType , only : soilhydrology_type use GridcellType , only : grc use LandunitType , only : lun @@ -46,7 +46,7 @@ module dynInitColumnsMod !----------------------------------------------------------------------- subroutine initialize_new_columns(bounds, cactive_prior, & - temperature_inst, waterstatebulk_inst) + temperature_inst, water_inst) ! ! !DESCRIPTION: ! Do initialization for all columns that are newly-active in this time step @@ -57,8 +57,8 @@ subroutine initialize_new_columns(bounds, cactive_prior, & ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds ! bounds logical , intent(in) :: cactive_prior( bounds%begc: ) ! column-level active flags from prior time step - type(temperature_type) , intent(inout) :: temperature_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst + type(temperature_type) , intent(inout) :: temperature_inst + type(water_type) , intent(inout) :: water_inst ! ! !LOCAL VARIABLES: integer :: c ! column index @@ -75,7 +75,7 @@ subroutine initialize_new_columns(bounds, cactive_prior, & c_template = initial_template_col_dispatcher(bounds, c, cactive_prior(bounds%begc:bounds%endc)) if (c_template /= TEMPLATE_NONE_FOUND) then call copy_state(c, c_template, & - temperature_inst, waterstatebulk_inst) + temperature_inst, water_inst) else write(iulog,*) subname// ' WARNING: No template column found to initialize newly-active column' write(iulog,*) '-- keeping the state that was already in memory, possibly from arbitrary initialization' @@ -211,7 +211,7 @@ end function initial_template_col_crop !----------------------------------------------------------------------- subroutine copy_state(c_new, c_template, & - temperature_inst, waterstatebulk_inst) + temperature_inst, water_inst) ! ! !DESCRIPTION: ! Copy a subset of state variables from a template column (c_template) to a newly- @@ -223,9 +223,10 @@ subroutine copy_state(c_new, c_template, & integer, intent(in) :: c_new ! index of newly-active column integer, intent(in) :: c_template ! index of column to use as a template type(temperature_type) , intent(inout) :: temperature_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst + type(water_type) , intent(inout) :: water_inst ! ! !LOCAL VARIABLES: + integer :: i character(len=*), parameter :: subname = 'copy_state' !----------------------------------------------------------------------- @@ -242,20 +243,27 @@ subroutine copy_state(c_new, c_template, & temperature_inst%t_soisno_col(c_new,1:) = temperature_inst%t_soisno_col(c_template,1:) - ! TODO(wjs, 2016-08-31) If we had more general uses of this initial template col - ! infrastructure (copying state between very different landunits), then we might need - ! to handle bedrock layers - e.g., zeroing out any water that would be added to a - ! bedrock layer(?). But for now we just use this initial template col infrastructure - ! for nat veg -> crop, for which the bedrock will be the same, so we're not dealing - ! with that complexity for now. - waterstatebulk_inst%h2osoi_liq_col(c_new,1:) = waterstatebulk_inst%h2osoi_liq_col(c_template,1:) - waterstatebulk_inst%h2osoi_ice_col(c_new,1:) = waterstatebulk_inst%h2osoi_ice_col(c_template,1:) - waterstatebulk_inst%h2osoi_vol_col(c_new,1:) = waterstatebulk_inst%h2osoi_vol_col(c_template,1:) + do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end - waterstatebulk_inst%wa_col(c_new) = waterstatebulk_inst%wa_col(c_template) + associate( & + waterstate_inst => water_inst%bulk_and_tracers(i)%waterstate_inst) - end subroutine copy_state + ! TODO(wjs, 2016-08-31) If we had more general uses of this initial template col + ! infrastructure (copying state between very different landunits), then we might need + ! to handle bedrock layers - e.g., zeroing out any water that would be added to a + ! bedrock layer(?). But for now we just use this initial template col infrastructure + ! for nat veg -> crop, for which the bedrock will be the same, so we're not dealing + ! with that complexity for now. + waterstate_inst%h2osoi_liq_col(c_new,1:) = waterstate_inst%h2osoi_liq_col(c_template,1:) + waterstate_inst%h2osoi_ice_col(c_new,1:) = waterstate_inst%h2osoi_ice_col(c_template,1:) + waterstate_inst%h2osoi_vol_col(c_new,1:) = waterstate_inst%h2osoi_vol_col(c_template,1:) + + waterstate_inst%wa_col(c_new) = waterstate_inst%wa_col(c_template) + end associate + end do + + end subroutine copy_state end module dynInitColumnsMod diff --git a/src/dyn_subgrid/dynSubgridDriverMod.F90 b/src/dyn_subgrid/dynSubgridDriverMod.F90 index ed0ca7a026..73c3f9e4b0 100644 --- a/src/dyn_subgrid/dynSubgridDriverMod.F90 +++ b/src/dyn_subgrid/dynSubgridDriverMod.F90 @@ -36,10 +36,7 @@ module dynSubgridDriverMod use PhotosynthesisMod , only : photosyns_type use SoilHydrologyType , only : soilhydrology_type use SoilStateType , only : soilstate_type - use WaterFluxBulkType , only : waterfluxbulk_type - use WaterStateBulkType , only : waterstatebulk_type - use WaterDiagnosticBulkType , only : waterdiagnosticbulk_type - use WaterBalanceType , only : waterbalance_type + use WaterType , only : water_type use TemperatureType , only : temperature_type use CropType , only : crop_type use glc2lndMod , only : glc2lnd_type @@ -155,9 +152,8 @@ end subroutine dynSubgrid_init !----------------------------------------------------------------------- subroutine dynSubgrid_driver(bounds_proc, & - urbanparams_inst, soilstate_inst, & - waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst, & - waterfluxbulk_inst, temperature_inst, energyflux_inst, & + urbanparams_inst, soilstate_inst, water_inst, & + temperature_inst, energyflux_inst, & canopystate_inst, photosyns_inst, crop_inst, glc2lnd_inst, bgc_vegetation_inst, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, & c13_soilbiogeochem_carbonstate_inst, c14_soilbiogeochem_carbonstate_inst, & @@ -184,10 +180,7 @@ subroutine dynSubgrid_driver(bounds_proc, type(bounds_type) , intent(in) :: bounds_proc ! processor-level bounds type(urbanparams_type) , intent(in) :: urbanparams_inst type(soilstate_type) , intent(in) :: soilstate_inst - type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst - type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst + type(water_type) , intent(inout) :: water_inst type(temperature_type) , intent(inout) :: temperature_inst type(energyflux_type) , intent(inout) :: energyflux_inst type(canopystate_type) , intent(inout) :: canopystate_inst @@ -228,8 +221,7 @@ subroutine dynSubgrid_driver(bounds_proc, filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_lakec, filter(nc)%lakec, & urbanparams_inst, soilstate_inst, & - waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst, & - waterfluxbulk_inst, temperature_inst, energyflux_inst) + water_inst, temperature_inst) call prior_weights%set_prior_weights(bounds_clump) call patch_state_updater%set_old_weights(bounds_clump) @@ -293,14 +285,14 @@ subroutine dynSubgrid_driver(bounds_proc, call initialize_new_columns(bounds_clump, & prior_weights%cactive(bounds_clump%begc:bounds_clump%endc), & - temperature_inst, waterstatebulk_inst) + temperature_inst, water_inst) call dyn_hwcontent_final(bounds_clump, & filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_lakec, filter(nc)%lakec, & urbanparams_inst, soilstate_inst, & - waterstatebulk_inst, waterdiagnosticbulk_inst, waterbalancebulk_inst, & - waterfluxbulk_inst, temperature_inst, energyflux_inst) + water_inst, & + temperature_inst, energyflux_inst) if (use_cn) then call bgc_vegetation_inst%DynamicAreaConservation(bounds_clump, nc, & diff --git a/src/dyn_subgrid/test/dynInitColumns_test/CMakeLists.txt b/src/dyn_subgrid/test/dynInitColumns_test/CMakeLists.txt index 59e6d13bb3..15647e0517 100644 --- a/src/dyn_subgrid/test/dynInitColumns_test/CMakeLists.txt +++ b/src/dyn_subgrid/test/dynInitColumns_test/CMakeLists.txt @@ -1,4 +1,4 @@ create_pFUnit_test(dynInitColumns test_dynInitColumns_exe "test_init_columns.pf" "") -target_link_libraries(test_dynInitColumns_exe clm csm_share) +target_link_libraries(test_dynInitColumns_exe clm csm_share esmf_wrf_timemgr) diff --git a/src/dyn_subgrid/test/dynInitColumns_test/test_init_columns.pf b/src/dyn_subgrid/test/dynInitColumns_test/test_init_columns.pf index 8ebf6f4849..829d9e8260 100644 --- a/src/dyn_subgrid/test/dynInitColumns_test/test_init_columns.pf +++ b/src/dyn_subgrid/test/dynInitColumns_test/test_init_columns.pf @@ -7,37 +7,54 @@ module test_init_columns use dynInitColumnsMod use ColumnType , only : col use LandunitType , only : lun - use GridcellType , only : grc + use GridcellType , only : grc + use landunit_varcon , only : istwet use decompMod , only : bounds_type use clm_varcon , only : ispval use clm_varpar , only : nlevsno, nlevgrnd use shr_kind_mod , only : r8 => shr_kind_r8 use TemperatureType , only : temperature_type - use WaterstateBulkType , only : waterstatebulk_type + use WaterType , only : water_type + use unittestWaterTypeFactory, only : unittest_water_type_factory_type + use WaterstateType , only : waterstate_type use dynColumnTemplateMod, only : TEMPLATE_NONE_FOUND implicit none - save - logical, allocatable, private :: cactive_prior(:) - integer :: c_new ! column index of the new column to initialize in some tests - integer :: l1 ! index of the landunit with landunit type 1 - integer :: l2 ! index of the landunit with landunit type 2 - - ! TODO(wjs, 2016-09-01) This test should be changed to use a class, with these being - ! instances of the class. Then the setup and cleanup routines here should be turned - ! into setup & teardown methods of the class. - type(temperature_type) :: temperature_vars - type(waterstatebulk_type) :: waterstatebulk_vars + @TestCase + type, extends(TestCase) :: TestInit + logical, allocatable :: cactive_prior(:) + integer :: c_new ! column index of the new column to initialize in some tests + integer :: l1 ! index of the landunit with landunit type 1 + integer :: l2 ! index of the landunit with landunit type 2 + type(temperature_type) :: temperature_inst + type(water_type) :: water_inst + type(unittest_water_type_factory_type) :: water_factory + contains + procedure :: setUp + procedure :: tearDown + procedure :: init_waterstate + end type TestInit contains - subroutine setup() - ! Set up variables needed for tests: various subgrid type variables, along with - ! bounds and cactive_prior. + subroutine setUp(this) + ! Set up variables needed for tests: various subgrid type variables, along with + ! bounds and this%cactive_prior. + ! + ! There is guaranteed to be at least one water tracer. ! - ! col%active and cactive_prior need to be set by specific tests + ! col%active and this%cactive_prior need to be set by specific tests + class(TestInit), intent(inout) :: this + integer :: c, lev + integer :: tracer_num + + call this%water_factory%init() + call this%water_factory%setup_before_subgrid( & + my_nlevsoi = 6, & + nlevgrnd_additional = 4, & + my_nlevsno = 3) ! Set up subgrid structure ! The weights (of both landunits and columns) and column types in the following are @@ -47,100 +64,123 @@ contains call unittest_add_gridcell() - call unittest_add_landunit(my_gi=gi, ltype=3, wtgcell=0.25_r8) + ! The first landunit is neither natural veg nor crop + call unittest_add_landunit(my_gi=gi, ltype=istwet, wtgcell=0.25_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.5_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.5_r8) call unittest_add_landunit(my_gi=gi, ltype=1, wtgcell=0.5_r8) - l1 = li + this%l1 = li call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) ! This column (the second column on the landunit with ltype=1) will be the target for ! some tests of initialization of a new column - c_new = ci + this%c_new = ci call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) call unittest_add_landunit(my_gi=gi, ltype=2, wtgcell=0.25_r8) - l2 = li + this%l2 = li call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.25_r8) call unittest_add_column(my_li=li, ctype=1, wtlunit=0.5_r8) call unittest_subgrid_setup_end() + call this%water_factory%setup_after_subgrid(snl = 0, dz = 1._r8) + call this%water_factory%create_water_type(this%water_inst, & + enable_isotopes = .true.) + col%active(begc:endc) = .false. - allocate(cactive_prior(bounds%begc:bounds%endc), source=.false.) - - nlevgrnd=10 - allocate(temperature_vars%t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) - allocate(waterstatebulk_vars%h2osoi_liq_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) - allocate(waterstatebulk_vars%h2osoi_ice_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) - allocate(waterstatebulk_vars%h2osoi_vol_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) - allocate(waterstatebulk_vars%wa_col(bounds%begc:bounds%endc)) + allocate(this%cactive_prior(bounds%begc:bounds%endc), source=.false.) + + allocate(this%temperature_inst%t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) do lev = -nlevsno+1, nlevgrnd do c = bounds%begc, bounds%endc - temperature_vars%t_soisno_col(c, lev) = c*1000 + lev - - ! Also need to initialize some waterstate_type variables, but we don't have any - ! assertions on them in this test - waterstatebulk_vars%h2osoi_liq_col(c, lev) = 0._r8 - waterstatebulk_vars%h2osoi_ice_col(c, lev) = 0._r8 - waterstatebulk_vars%h2osoi_vol_col(c, lev) = 0._r8 + this%temperature_inst%t_soisno_col(c, lev) = c*1000 + lev end do end do - ! Also need to initialize some other variables for which we don't have any assertions - do c = bounds%begc, bounds%endc - waterstatebulk_vars%wa_col(c) = 0._r8 + call this%init_waterstate(this%water_inst%waterstatebulk_inst, & + tracer_num = 0) + @assertGreaterThan(this%water_inst%tracers_end, this%water_inst%tracers_beg) + do tracer_num = this%water_inst%tracers_beg, this%water_inst%tracers_end + call this%init_waterstate(this%water_inst%bulk_and_tracers(tracer_num)%waterstate_inst, & + tracer_num = tracer_num) end do - end subroutine setup - subroutine cleanup() + end subroutine setUp + + subroutine tearDown(this) ! clean up stuff set up in setup() + class(TestInit), intent(inout) :: this + call this%water_factory%teardown(this%water_inst) + deallocate(this%temperature_inst%t_soisno_col) call unittest_subgrid_teardown() - deallocate(cactive_prior) - deallocate(temperature_vars%t_soisno_col) - deallocate(waterstatebulk_vars%h2osoi_liq_col) - deallocate(waterstatebulk_vars%h2osoi_ice_col) - deallocate(waterstatebulk_vars%h2osoi_vol_col) - deallocate(waterstatebulk_vars%wa_col) + end subroutine tearDown + + subroutine init_waterstate(this, waterstate_inst, tracer_num) + ! Initialize variables that need to be set for these tests in one waterstate instance + class(TestInit), intent(in) :: this + class(waterstate_type), intent(inout) :: waterstate_inst + integer, intent(in) :: tracer_num ! used to give each tracer unique values - end subroutine cleanup + integer :: lev, c + + do lev = -nlevsno+1, nlevgrnd + do c = bounds%begc, bounds%endc + ! 137 is a nice arbitrary number, larger than any column index (so we won't + ! have the same value in any column) + waterstate_inst%h2osoi_liq_col(c, lev) = & + c*(137._r8 + 37._r8*tracer_num) + lev + + ! Also need to initialize some other waterstate_type variables, but we don't + ! have any assertions on them in this test, so just initialize them to 0 + waterstate_inst%h2osoi_ice_col(c, lev) = 0._r8 + end do + end do + + ! Also need to initialize some other variables for which we don't have any assertions + do lev = 1, nlevgrnd + do c = bounds%begc, bounds%endc + waterstate_inst%h2osoi_vol_col(c, lev) = 0._r8 + end do + end do + do c = bounds%begc, bounds%endc + waterstate_inst%wa_col(c) = 0._r8 + end do + end subroutine init_waterstate ! ------------------------------------------------------------------------ ! Tests of initial_template_col_crop ! ------------------------------------------------------------------------ @Test - subroutine test_crop_active_in_soil_and_crop() + subroutine test_crop_active_in_soil_and_crop(this) ! there are active columns both on the soil & crop landunits; should take the soil one - call setup() - cactive_prior(:) = .true. - @assertEqual(lun%coli(l1), initial_template_col_crop(bounds, c_new, cactive_prior)) - call cleanup() + class(TestInit), intent(inout) :: this + this%cactive_prior(:) = .true. + @assertEqual(lun%coli(this%l1), initial_template_col_crop(bounds, this%c_new, this%cactive_prior)) end subroutine test_crop_active_in_soil_and_crop @Test - subroutine test_crop_no_soil() + subroutine test_crop_no_soil(this) ! no soil landunit, should take a crop column - call setup() - cactive_prior(:) = .true. + class(TestInit), intent(inout) :: this + this%cactive_prior(:) = .true. grc%landunit_indices(1,gi) = ispval - @assertEqual(lun%coli(l2), initial_template_col_crop(bounds, c_new, cactive_prior)) - call cleanup() + @assertEqual(lun%coli(this%l2), initial_template_col_crop(bounds, this%c_new, this%cactive_prior)) end subroutine test_crop_no_soil @Test - subroutine test_crop_no_soil_or_crop() + subroutine test_crop_no_soil_or_crop(this) ! no soil or crop landunits, should give TEMPLATE_NONE_FOUND - call setup() - cactive_prior(:) = .true. + class(TestInit), intent(inout) :: this + this%cactive_prior(:) = .true. grc%landunit_indices(1:2,gi) = ispval - @assertEqual(TEMPLATE_NONE_FOUND, initial_template_col_crop(bounds, c_new, cactive_prior)) - call cleanup() + @assertEqual(TEMPLATE_NONE_FOUND, initial_template_col_crop(bounds, this%c_new, this%cactive_prior)) end subroutine test_crop_no_soil_or_crop ! ------------------------------------------------------------------------ @@ -151,67 +191,81 @@ contains ! ------------------------------------------------------------------------ @Test - subroutine test_initialize_new_columns_none() + subroutine test_initialize_new_columns_none(this) ! Nothing to initialize + class(TestInit), intent(inout) :: this real(r8), allocatable :: t_soisno_expected(:,:) - call setup() ! col%active and cactive_prior are a mix of true/true, false/false and false/true, so ! there's nothing to initialize col%active(:) = .true. - cactive_prior(:) = .true. - col%active(lun%coli(l2)+1) = .false. - cactive_prior(lun%coli(l2)+1) = .false. - col%active(lun%coli(l2)+2) = .false. - t_soisno_expected = temperature_vars%t_soisno_col - call initialize_new_columns(bounds, cactive_prior, & - temperature_vars, waterstatebulk_vars) - @assertEqual(t_soisno_expected, temperature_vars%t_soisno_col) - call cleanup() + this%cactive_prior(:) = .true. + col%active(lun%coli(this%l2)+1) = .false. + this%cactive_prior(lun%coli(this%l2)+1) = .false. + col%active(lun%coli(this%l2)+2) = .false. + t_soisno_expected = this%temperature_inst%t_soisno_col + call initialize_new_columns(bounds, this%cactive_prior, & + this%temperature_inst, this%water_inst) + @assertEqual(t_soisno_expected, this%temperature_inst%t_soisno_col) end subroutine test_initialize_new_columns_none @Test - subroutine test_initialize_new_columns_TEMPLATE_NONE_FOUND() + subroutine test_initialize_new_columns_TEMPLATE_NONE_FOUND(this) ! Something to initialize, but template_col results in TEMPLATE_NONE_FOUND: state should remain ! the same as before + class(TestInit), intent(inout) :: this real(r8), allocatable :: t_soisno_expected(:,:) - call setup() col%active(:) = .false. - col%active(lun%coli(l2)+1) = .true. + col%active(lun%coli(this%l2)+1) = .true. ! all cactive_prior points were false, so there's nothing to use as a template: - cactive_prior(:) = .false. - t_soisno_expected = temperature_vars%t_soisno_col - call initialize_new_columns(bounds, cactive_prior, & - temperature_vars, waterstatebulk_vars) - @assertEqual(t_soisno_expected, temperature_vars%t_soisno_col) - call cleanup() + this%cactive_prior(:) = .false. + t_soisno_expected = this%temperature_inst%t_soisno_col + call initialize_new_columns(bounds, this%cactive_prior, & + this%temperature_inst, this%water_inst) + @assertEqual(t_soisno_expected, this%temperature_inst%t_soisno_col) end subroutine test_initialize_new_columns_TEMPLATE_NONE_FOUND @Test - subroutine test_initialize_new_columns_copy_state() + subroutine test_initialize_new_columns_copy_state(this) ! Something to initialize, which results in a state copy + class(TestInit), intent(inout) :: this real(r8), allocatable :: t_soisno_expected(:,:) + real(r8), allocatable :: h2osoi_liq_expected(:,:) + real(r8), allocatable :: h2osoi_liq_tracer_expected(:,:) integer :: source_col, dest_col - call setup() - col%active(:) = .false. - dest_col = lun%coli(l2) + 1 + dest_col = lun%coli(this%l2) + 1 col%active(dest_col) = .true. - cactive_prior(:) = .false. - source_col = lun%coli(l1) + 1 - cactive_prior(source_col) = .true. + this%cactive_prior(:) = .false. + source_col = lun%coli(this%l1) + 1 + this%cactive_prior(source_col) = .true. - t_soisno_expected = temperature_vars%t_soisno_col + ! We assume that, if the first tracer is handled correctly, then all of them are + associate( & + tracer_h2osoi_liq_col => this%water_inst%bulk_and_tracers(this%water_inst%tracers_beg)%waterstate_inst%h2osoi_liq_col) + + t_soisno_expected = this%temperature_inst%t_soisno_col + h2osoi_liq_expected = this%water_inst%waterstatebulk_inst%h2osoi_liq_col + h2osoi_liq_tracer_expected = tracer_h2osoi_liq_col ! In the following, note that only the belowground portion (starting with level 1) is ! copied: - t_soisno_expected(dest_col,1:) = temperature_vars%t_soisno_col(source_col,1:) + t_soisno_expected(dest_col,1:) = & + this%temperature_inst%t_soisno_col(source_col,1:) + h2osoi_liq_expected(dest_col,1:) = & + this%water_inst%waterstatebulk_inst%h2osoi_liq_col(source_col,1:) + h2osoi_liq_tracer_expected(dest_col,1:) = & + tracer_h2osoi_liq_col(source_col,1:) + + call initialize_new_columns(bounds, this%cactive_prior, & + this%temperature_inst, this%water_inst) + + @assertEqual(t_soisno_expected, this%temperature_inst%t_soisno_col) + @assertEqual(h2osoi_liq_expected, this%water_inst%waterstatebulk_inst%h2osoi_liq_col) + @assertEqual(h2osoi_liq_tracer_expected, tracer_h2osoi_liq_col) - call initialize_new_columns(bounds, cactive_prior, & - temperature_vars, waterstatebulk_vars) + end associate - @assertEqual(t_soisno_expected, temperature_vars%t_soisno_col) - call cleanup() end subroutine test_initialize_new_columns_copy_state end module test_init_columns diff --git a/src/main/atm2lndMod.F90 b/src/main/atm2lndMod.F90 index 8a13c87686..7a5bd00303 100644 --- a/src/main/atm2lndMod.F90 +++ b/src/main/atm2lndMod.F90 @@ -22,6 +22,7 @@ module atm2lndMod use LandunitType , only : lun use ColumnType , only : col use landunit_varcon, only : istice_mec + use WaterType , only : water_type use Wateratm2lndBulkType, only : wateratm2lndbulk_type ! ! !PUBLIC TYPES: @@ -30,19 +31,20 @@ module atm2lndMod save ! ! !PUBLIC MEMBER FUNCTIONS: - public :: downscale_forcings ! Downscale atm forcing fields from gridcell to column + public :: set_atm2lnd_water_tracers ! Set tracer values for the atm2lnd water quantities + public :: downscale_forcings ! Downscale atm forcing fields from gridcell to column - ! The following routines are public for the sake of unit testing; they should not be + ! The following routine is public for the sake of unit testing; it should not be ! called by production code outside this module - public :: partition_precip ! Partition precipitation into rain/snow - public :: sens_heat_from_precip_conversion ! Compute sensible heat flux needed to compensate for rain-snow conversion + public :: partition_precip ! Partition precipitation into rain/snow ! ! !PRIVATE MEMBER FUNCTIONS: - private :: rhos ! calculate atmospheric density - private :: repartition_rain_snow_one_col ! Re-partition precipitation for a single column - private :: downscale_longwave ! Downscale longwave radiation from gridcell to column - private :: build_normalization ! Compute normalization factors so that downscaled fields are conservative - private :: check_downscale_consistency ! Check consistency of downscaling + private :: rhos ! calculate atmospheric density + private :: repartition_rain_snow_one_col ! Re-partition precipitation for a single column + private :: sens_heat_from_precip_conversion ! Compute sensible heat flux needed to compensate for rain-snow conversion + private :: downscale_longwave ! Downscale longwave radiation from gridcell to column + private :: build_normalization ! Compute normalization factors so that downscaled fields are conservative + private :: check_downscale_consistency ! Check consistency of downscaling character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -50,6 +52,43 @@ module atm2lndMod contains + !----------------------------------------------------------------------- + subroutine set_atm2lnd_water_tracers(bounds, num_allc, filter_allc, water_inst) + ! + ! !DESCRIPTION: + ! Set tracer values for the atm2lnd water quantities + ! + ! Should be called after downscale_forcings + ! + ! !ARGUMENTS: + type(bounds_type) , intent(in) :: bounds + integer , intent(in) :: num_allc ! number of column points in filter_allc + integer , intent(in) :: filter_allc(:) ! column filter for all points + type(water_type) , intent(in) :: water_inst + ! + ! !LOCAL VARIABLES: + integer :: i + + character(len=*), parameter :: subname = 'set_atm2lnd_water_tracers' + !----------------------------------------------------------------------- + + do i = water_inst%tracers_beg, water_inst%tracers_end + associate( & + wateratm2lnd_inst => water_inst%bulk_and_tracers(i)%wateratm2lnd_inst) + + if (.not. wateratm2lnd_inst%IsCommunicatedWithCoupler()) then + call wateratm2lnd_inst%SetNondownscaledTracers( & + bounds, water_inst%wateratm2lndbulk_inst) + end if + + call wateratm2lnd_inst%SetDownscaledTracers( & + bounds, num_allc, filter_allc, water_inst%wateratm2lndbulk_inst) + + end associate + end do + + end subroutine set_atm2lnd_water_tracers + !----------------------------------------------------------------------- subroutine downscale_forcings(bounds, & topo_inst, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh_precip_conversion) @@ -263,8 +302,8 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh ! ! !LOCAL VARIABLES: integer :: c,l,g ! indices - real(r8) :: rain_old ! rain before conversion - real(r8) :: snow_old ! snow before conversion + real(r8) :: rain_orig ! rain before conversion + real(r8) :: snow_orig ! snow before conversion real(r8) :: all_snow_t ! temperature at which all precip falls as snow (K) real(r8) :: frac_rain_slope ! slope of the frac_rain vs. temperature relationship @@ -279,9 +318,11 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh forc_snow_g => wateratm2lndbulk_inst%forc_snow_not_downscaled_grc , & ! Input: [real(r8) (:)] snow rate [mm/s] ! Column-level downscaled fields: - forc_t_c => atm2lnd_inst%forc_t_downscaled_col , & ! Input: [real(r8) (:)] atmospheric temperature (Kelvin) - forc_rain_c => wateratm2lndbulk_inst%forc_rain_downscaled_col , & ! Output: [real(r8) (:)] rain rate [mm/s] - forc_snow_c => wateratm2lndbulk_inst%forc_snow_downscaled_col & ! Output: [real(r8) (:)] snow rate [mm/s] + forc_t_c => atm2lnd_inst%forc_t_downscaled_col , & ! Input: [real(r8) (:)] atmospheric temperature (Kelvin) + forc_rain_c => wateratm2lndbulk_inst%forc_rain_downscaled_col , & ! Output: [real(r8) (:)] rain rate [mm/s] + forc_snow_c => wateratm2lndbulk_inst%forc_snow_downscaled_col , & ! Output: [real(r8) (:)] snow rate [mm/s] + rain_to_snow_conversion_c => wateratm2lndbulk_inst%rain_to_snow_conversion_col , & ! Output: [real(r8) (:)] amount of rain converted to snow via precipitation repartition [mm/s] + snow_to_rain_conversion_c => wateratm2lndbulk_inst%snow_to_rain_conversion_col & ! Output: [real(r8) (:)] amount of snow converted to rain via precipitation repartition [mm/s] ) ! Initialize column forcing @@ -290,6 +331,8 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh g = col%gridcell(c) forc_rain_c(c) = forc_rain_g(g) forc_snow_c(c) = forc_snow_g(g) + rain_to_snow_conversion_c(c) = 0._r8 + snow_to_rain_conversion_c(c) = 0._r8 eflx_sh_precip_conversion(c) = 0._r8 end if end do @@ -299,8 +342,8 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh do c = bounds%begc, bounds%endc if (col%active(c)) then l = col%landunit(c) - rain_old = forc_rain_c(c) - snow_old = forc_snow_c(c) + rain_orig = forc_rain_c(c) + snow_orig = forc_snow_c(c) if (lun%itype(l) == istice_mec) then all_snow_t = atm2lnd_inst%params%precip_repartition_glc_all_snow_t frac_rain_slope = atm2lnd_inst%params%precip_repartition_glc_frac_rain_slope @@ -314,11 +357,15 @@ subroutine partition_precip(bounds, atm2lnd_inst, wateratm2lndbulk_inst, eflx_sh frac_rain_slope = frac_rain_slope, & rain = forc_rain_c(c), & snow = forc_snow_c(c)) + if (forc_rain_c(c) > rain_orig) then + snow_to_rain_conversion_c(c) = forc_rain_c(c) - rain_orig + end if + if (forc_snow_c(c) > snow_orig) then + rain_to_snow_conversion_c(c) = forc_snow_c(c) - snow_orig + end if call sens_heat_from_precip_conversion(& - rain_old = rain_old, & - snow_old = snow_old, & - rain_new = forc_rain_c(c), & - snow_new = forc_snow_c(c), & + rain_to_snow = rain_to_snow_conversion_c(c), & + snow_to_rain = snow_to_rain_conversion_c(c), & sens_heat_flux = eflx_sh_precip_conversion(c)) end if end do @@ -363,40 +410,28 @@ subroutine repartition_rain_snow_one_col(temperature, all_snow_t, frac_rain_slop end subroutine repartition_rain_snow_one_col !----------------------------------------------------------------------- - subroutine sens_heat_from_precip_conversion(rain_old, snow_old, rain_new, snow_new, & - sens_heat_flux) + subroutine sens_heat_from_precip_conversion(rain_to_snow, snow_to_rain, sens_heat_flux) ! ! !DESCRIPTION: - ! Given old and new rain and snow amounts, compute the sensible heat flux needed to - ! compensate for the rain-snow conversion. + ! Given conversion fluxes from rain to snow and snow to rain, compute the sensible + ! heat flux needed to compensate for the rain-snow conversion. ! ! !USES: ! ! !ARGUMENTS: - real(r8), intent(in) :: rain_old ! [mm/s] - real(r8), intent(in) :: snow_old ! [(mm water equivalent)/s] - real(r8), intent(in) :: rain_new ! [mm/s] - real(r8), intent(in) :: snow_new ! [(mm water equivalent)/s] + real(r8), intent(in) :: rain_to_snow ! amount of rain converted to snow [mm/s] + real(r8), intent(in) :: snow_to_rain ! amount of snow converted to rain [mm/s] real(r8), intent(out) :: sens_heat_flux ! [W/m^2] ! ! !LOCAL VARIABLES: - real(r8) :: total_old - real(r8) :: total_new - real(r8) :: rain_to_snow ! net conversion of rain to snow - real(r8), parameter :: mm_to_m = 1.e-3_r8 ! multiply by this to convert from mm to m real(r8), parameter :: tol = 1.e-13_r8 ! relative tolerance for error checks character(len=*), parameter :: subname = 'sens_heat_from_precip_conversion' !----------------------------------------------------------------------- - total_old = rain_old + snow_old - total_new = rain_new + snow_new - SHR_ASSERT(abs(total_new - total_old) <= (tol * total_old), subname//' ERROR: mismatch between old and new totals') - ! rain to snow releases energy, so results in a positive heat flux to atm - rain_to_snow = snow_new - snow_old - sens_heat_flux = rain_to_snow * mm_to_m * denh2o * hfus + sens_heat_flux = (rain_to_snow - snow_to_rain) * mm_to_m * denh2o * hfus end subroutine sens_heat_from_precip_conversion diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index ea2035cb76..d597ab40f0 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -61,7 +61,7 @@ module clm_driver ! use filterMod , only : setFilters ! - use atm2lndMod , only : downscale_forcings + use atm2lndMod , only : downscale_forcings, set_atm2lnd_water_tracers use lnd2atmMod , only : lnd2atm use lnd2glcMod , only : lnd2glc_type ! @@ -304,9 +304,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call t_startf('dyn_subgrid') call dynSubgrid_driver(bounds_proc, & - urbanparams_inst, soilstate_inst, & - water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & - water_inst%waterbalancebulk_inst, water_inst%waterfluxbulk_inst, & + urbanparams_inst, soilstate_inst, water_inst, & temperature_inst, energyflux_inst, & canopystate_inst, photosyns_inst, crop_inst, glc2lnd_inst, bgc_vegetation_inst, & soilbiogeochem_state_inst, soilbiogeochem_carbonstate_inst, & @@ -338,8 +336,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call BeginWaterBalance(bounds_clump, & filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_lakec, filter(nc)%lakec, & - soilhydrology_inst, water_inst%waterstatebulk_inst, & - water_inst%waterdiagnosticbulk_inst, water_inst%waterbalancebulk_inst) + water_inst, soilhydrology_inst) call t_stopf('begwbal') @@ -406,8 +403,8 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro filter(nc)%num_nolakep, filter(nc)%nolakep, & filter(nc)%num_soilp , filter(nc)%soilp, & canopystate_inst, water_inst%waterstatebulk_inst, & - water_inst%waterbalancebulk_inst, water_inst%waterdiagnosticbulk_inst, & - water_inst%waterfluxbulk_inst, energyflux_inst) + water_inst%waterdiagnosticbulk_inst, & + energyflux_inst) call topo_inst%UpdateTopo(bounds_clump, & filter(nc)%num_icemecc, filter(nc)%icemecc, & @@ -418,6 +415,20 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro topo_inst, atm2lnd_inst, water_inst%wateratm2lndbulk_inst, & eflx_sh_precip_conversion = energyflux_inst%eflx_sh_precip_conversion_col(bounds_clump%begc:bounds_clump%endc)) + call set_atm2lnd_water_tracers(bounds_clump, & + filter(nc)%num_allc, filter(nc)%allc, & + water_inst) + + if (water_inst%DoConsistencyCheck()) then + ! BUG(wjs, 2018-09-05, ESCOMP/ctsm#498) Eventually do tracer consistency checks + ! every time step + if (get_nstep() == 0) then + call t_startf("tracer_consistency_check") + call water_inst%TracerConsistencyCheck(bounds_clump, 'after downscale_forcings') + call t_stopf("tracer_consistency_check") + end if + end if + ! Update filters that depend on variables set in clm_drv_init call setExposedvegpFilter(bounds_clump, & @@ -1223,7 +1234,7 @@ subroutine clm_drv_init(bounds, & num_nolakep, filter_nolakep, & num_soilp , filter_soilp, & canopystate_inst, waterstatebulk_inst, & - waterbalancebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, energyflux_inst) + waterdiagnosticbulk_inst, energyflux_inst) ! ! !DESCRIPTION: ! Initialization of clm driver variables needed from previous timestep @@ -1249,9 +1260,7 @@ subroutine clm_drv_init(bounds, & integer , intent(in) :: filter_soilp(:) ! patch filter for soil points type(canopystate_type), intent(inout) :: canopystate_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst - type(waterbalance_type) , intent(inout) :: waterbalancebulk_inst type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst - type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst type(energyflux_type) , intent(inout) :: energyflux_inst ! ! !LOCAL VARIABLES: @@ -1262,10 +1271,8 @@ subroutine clm_drv_init(bounds, & associate( & snl => col%snl , & ! Input: [integer (:) ] number of snow layers - h2osno => waterstatebulk_inst%h2osno_col , & ! Input: [real(r8) (:) ] snow water (mm H2O) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) - h2osno_old => waterbalancebulk_inst%h2osno_old_col , & ! Output: [real(r8) (:) ] snow water (mm H2O) at previous time step frac_iceold => waterdiagnosticbulk_inst%frac_iceold_col , & ! Output: [real(r8) (:,:) ] fraction of ice relative to the tot water elai => canopystate_inst%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow @@ -1286,9 +1293,6 @@ subroutine clm_drv_init(bounds, & end do do c = bounds%begc,bounds%endc - ! Save snow mass at previous time step - h2osno_old(c) = h2osno(c) - ! Reset flux from beneath soil/ice column eflx_bot(c) = 0._r8 end do diff --git a/src/main/lnd2atmMod.F90 b/src/main/lnd2atmMod.F90 index 71d3909515..dd9c080a76 100644 --- a/src/main/lnd2atmMod.F90 +++ b/src/main/lnd2atmMod.F90 @@ -70,7 +70,7 @@ subroutine lnd2atm_minimal(bounds, & ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds - type(water_type) , intent(in) :: water_inst + type(water_type) , intent(inout) :: water_inst type(surfalb_type) , intent(in) :: surfalb_inst type(energyflux_type) , intent(in) :: energyflux_inst type(lnd2atm_type) , intent(inout) :: lnd2atm_inst @@ -84,27 +84,20 @@ subroutine lnd2atm_minimal(bounds, & real(r8), parameter :: convertgC2kgCO2 = 1.0e-3_r8 * (amCO2/amC) !------------------------------------------------------------------------ - call c2g(bounds, & - water_inst%waterstatebulk_inst%h2osno_col (bounds%begc:bounds%endc), & - water_inst%waterlnd2atmbulk_inst%h2osno_grc (bounds%begg:bounds%endg), & - c2l_scale_type= 'urbanf', l2g_scale_type='unity') - - do g = bounds%begg,bounds%endg - water_inst%waterlnd2atmbulk_inst%h2osno_grc(g) = water_inst%waterlnd2atmbulk_inst%h2osno_grc(g)/1000._r8 - end do + do i = water_inst%bulk_and_tracers_beg, water_inst%bulk_and_tracers_end + associate(bulk_or_tracer => water_inst%bulk_and_tracers(i)) + call c2g(bounds, & + bulk_or_tracer%waterstate_inst%h2osno_col(bounds%begc:bounds%endc), & + bulk_or_tracer%waterlnd2atm_inst%h2osno_grc(bounds%begg:bounds%endg), & + c2l_scale_type= 'urbanf', l2g_scale_type='unity') - do i = 1, water_inst%num_tracers - call c2g(bounds, & - water_inst%waterstate_tracer_inst(i)%h2osno_col (bounds%begc:bounds%endc), & - water_inst%waterlnd2atm_tracer_inst(i)%h2osno_grc (bounds%begg:bounds%endg), & - c2l_scale_type= 'urbanf', l2g_scale_type='unity') - - do g = bounds%begg,bounds%endg - water_inst%waterlnd2atm_tracer_inst(i)%h2osno_grc(g) = water_inst%waterlnd2atm_tracer_inst(i)%h2osno_grc(g)/1000._r8 - end do + do g = bounds%begg,bounds%endg + bulk_or_tracer%waterlnd2atm_inst%h2osno_grc(g) = & + bulk_or_tracer%waterlnd2atm_inst%h2osno_grc(g)/1000._r8 + end do + end associate end do - call c2g(bounds, nlevgrnd, & water_inst%waterstatebulk_inst%h2osoi_vol_col (bounds%begc:bounds%endc, :), & water_inst%waterlnd2atmbulk_inst%h2osoi_vol_grc (bounds%begg:bounds%endg, :), & diff --git a/src/main/test/atm2lnd_test/CMakeLists.txt b/src/main/test/atm2lnd_test/CMakeLists.txt index 018b875a27..e42192b45b 100644 --- a/src/main/test/atm2lnd_test/CMakeLists.txt +++ b/src/main/test/atm2lnd_test/CMakeLists.txt @@ -1,7 +1,6 @@ set(pfunit_sources test_downscale_forcings.pf - test_partition_precip.pf - test_sens_heat_from_precip_conversion.pf) + test_partition_precip.pf) create_pFUnit_test(atm2lnd test_atm2lnd_exe "${pfunit_sources}" "") diff --git a/src/main/test/atm2lnd_test/test_partition_precip.pf b/src/main/test/atm2lnd_test/test_partition_precip.pf index 750b7bc2f8..878bdac4f6 100644 --- a/src/main/test/atm2lnd_test/test_partition_precip.pf +++ b/src/main/test/atm2lnd_test/test_partition_precip.pf @@ -102,9 +102,13 @@ contains associate(& rain_col => this%wateratm2lndbulk_inst%forc_rain_downscaled_col, & - snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) + snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col, & + rain_to_snow_conversion_col => this%wateratm2lndbulk_inst%rain_to_snow_conversion_col, & + snow_to_rain_conversion_col => this%wateratm2lndbulk_inst%snow_to_rain_conversion_col) @assertEqual(0._r8, rain_col(begc), tolerance=tol) @assertEqual(3._r8, snow_col(begc), tolerance=tol) + @assertEqual(1._r8, rain_to_snow_conversion_col(begc), tolerance=tol) + @assertEqual(0._r8, snow_to_rain_conversion_col(begc)) end associate end subroutine lowTemp_resultsInCorrectPartitioning @@ -119,28 +123,54 @@ contains associate(& rain_col => this%wateratm2lndbulk_inst%forc_rain_downscaled_col, & - snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) + snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col, & + rain_to_snow_conversion_col => this%wateratm2lndbulk_inst%rain_to_snow_conversion_col, & + snow_to_rain_conversion_col => this%wateratm2lndbulk_inst%snow_to_rain_conversion_col) @assertEqual(3._r8, rain_col(begc), tolerance=tol) @assertEqual(0._r8, snow_col(begc), tolerance=tol) + @assertEqual(2._r8, snow_to_rain_conversion_col(begc), tolerance=tol) + @assertEqual(0._r8, rain_to_snow_conversion_col(begc)) end associate end subroutine highTemp_resultsInCorrectPartitioning @Test - subroutine intermediateTemp_resultsInCorrectPartitioning(this) + subroutine intermediateTemp_resultsInCorrectPartitioningAndHeatFlux(this) + ! Unlike other tests, where we check for either correct partitioning or correct heat + ! flux in a given test: This test checks both for convenience (to avoid duplication + ! between two tests which would require us to update both of them if the ramp for + ! rain-snow conversion changed). class(TestPartitionPrecip), intent(inout) :: this + real(r8), parameter :: rain_orig = 1._r8 + real(r8), parameter :: snow_orig = 2._r8 + real(r8) :: tot_precip + real(r8) :: expected_rain + real(r8) :: expected_snow + real(r8) :: expected_heat_flux call setup_single_veg_patch(pft_type=1) - call this%set_inputs(rain=[1._r8], snow=[2._r8], temperature=[SHR_CONST_TKFRZ + 1.5_r8]) + call this%set_inputs(rain=[rain_orig], snow=[snow_orig], temperature=[SHR_CONST_TKFRZ + 1.5_r8]) call partition_precip(bounds, this%atm2lnd_inst, this%wateratm2lndbulk_inst, this%sh_from_conversion) associate(& rain_col => this%wateratm2lndbulk_inst%forc_rain_downscaled_col, & - snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) - @assertEqual(3._r8 * 0.75_r8, rain_col(begc), tolerance=tol) - @assertEqual(3._r8 * 0.25_r8, snow_col(begc), tolerance=tol) + snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col, & + rain_to_snow_conversion_col => this%wateratm2lndbulk_inst%rain_to_snow_conversion_col, & + snow_to_rain_conversion_col => this%wateratm2lndbulk_inst%snow_to_rain_conversion_col) + + tot_precip = rain_orig + snow_orig + expected_rain = tot_precip * 0.75_r8 + expected_snow = tot_precip * 0.25_r8 + @assertEqual(expected_rain, rain_col(begc), tolerance=tol) + @assertEqual(expected_snow, snow_col(begc), tolerance=tol) + @assertEqual((expected_rain - rain_orig), snow_to_rain_conversion_col(begc), tolerance=tol) + @assertEqual(0._r8, rain_to_snow_conversion_col(begc)) + + ! Snow to rain extracts energy, so results in a negative heat flux to atm + expected_heat_flux = (rain_orig - expected_rain) * mm_to_m * denh2o * hfus + @assertEqual([expected_heat_flux], this%sh_from_conversion, tolerance=tol) end associate - end subroutine intermediateTemp_resultsInCorrectPartitioning + end subroutine intermediateTemp_resultsInCorrectPartitioningAndHeatFlux @Test subroutine intermediateTemp_glacier_resultsInCorrectPartitioning(this) @@ -162,6 +192,9 @@ contains snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) @assertEqual(3._r8 * 0.75_r8, rain_col(begc), tolerance=tol) @assertEqual(3._r8 * 0.25_r8, snow_col(begc), tolerance=tol) + ! Don't bother checking rain_to_snow_conversion_col and snow_to_rain_conversion_col + ! here: we don't gain anything from this beyond what's already checked in + ! intermediateTemp_resultsInCorrectPartitioningAndHeatFlux end associate end subroutine intermediateTemp_glacier_resultsInCorrectPartitioning @@ -221,9 +254,13 @@ contains associate(& rain_col => this%wateratm2lndbulk_inst%forc_rain_downscaled_col, & - snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) + snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col, & + rain_to_snow_conversion_col => this%wateratm2lndbulk_inst%rain_to_snow_conversion_col, & + snow_to_rain_conversion_col => this%wateratm2lndbulk_inst%snow_to_rain_conversion_col) @assertEqual([1._r8], rain_col) @assertEqual([2._r8], snow_col) + @assertEqual([0._r8], rain_to_snow_conversion_col) + @assertEqual([0._r8], snow_to_rain_conversion_col) @assertEqual([0._r8], this%sh_from_conversion) end associate end subroutine repartitionFlagFalse_resultsInNoChange @@ -236,6 +273,8 @@ contains real(r8), parameter :: temp(2) = [290._r8, 250._r8] real(r8), parameter :: rain_expected(2) = [3._r8, 0._r8] real(r8), parameter :: snow_expected(2) = [0._r8, 7._r8] + real(r8), parameter :: rain_to_snow_expected(2) = [0._r8, 3._r8] + real(r8), parameter :: snow_to_rain_expected(2) = [2._r8, 0._r8] call setup_ncells_single_veg_patch(ncells=2, pft_type=1) call this%set_inputs(rain=rain, snow=snow, temperature=temp) @@ -244,9 +283,13 @@ contains associate(& rain_col => this%wateratm2lndbulk_inst%forc_rain_downscaled_col, & - snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col) + snow_col => this%wateratm2lndbulk_inst%forc_snow_downscaled_col, & + rain_to_snow_conversion_col => this%wateratm2lndbulk_inst%rain_to_snow_conversion_col, & + snow_to_rain_conversion_col => this%wateratm2lndbulk_inst%snow_to_rain_conversion_col) @assertEqual(rain_expected, rain_col, tolerance=tol) @assertEqual(snow_expected, snow_col, tolerance=tol) + @assertEqual(rain_to_snow_expected, rain_to_snow_conversion_col, tolerance=tol) + @assertEqual(snow_to_rain_expected, snow_to_rain_conversion_col, tolerance=tol) end associate end subroutine multiPoint_resultsInCorrectPartitioning diff --git a/src/main/test/atm2lnd_test/test_sens_heat_from_precip_conversion.pf b/src/main/test/atm2lnd_test/test_sens_heat_from_precip_conversion.pf deleted file mode 100644 index 44ca36b9ec..0000000000 --- a/src/main/test/atm2lnd_test/test_sens_heat_from_precip_conversion.pf +++ /dev/null @@ -1,43 +0,0 @@ -module test_sens_heat_from_precip_conversion - - ! Tests of atm2lndMod: sens_heat_from_precip_conversion - ! This module just tests edge cases that would be difficult to test from the - ! multi-point wrapper. - - use pfunit_mod - use atm2lndMod - use shr_kind_mod, only : r8 => shr_kind_r8 - use clm_varcon, only : hfus ! latent heat of fusion for ice [J/kg] - use clm_varcon, only : denh2o ! density of liquid water [kg/m3] - - implicit none - - real(r8), parameter :: tol = 1.e-13_r8 - real(r8), parameter :: mm_to_m = 1.e-3_r8 ! multiply by this to convert from mm to m - -contains - - @Test - subroutine partialConversion_resultsInCorrectHeatFlux() - real(r8), parameter :: rain_old = 2._r8 ! [mm] - real(r8), parameter :: snow_old = 5._r8 ! [mm] - real(r8), parameter :: rain_new = 6._r8 ! [mm] - real(r8), parameter :: snow_new = 1._r8 ! [mm] - real(r8) :: sens_heat_flux ! [W/m2 to atm] - real(r8) :: expected - - call sens_heat_from_precip_conversion( & - rain_old = rain_old, & - snow_old = snow_old, & - rain_new = rain_new, & - snow_new = snow_new, & - sens_heat_flux = sens_heat_flux) - - ! Snow to rain extracts energy, so results in a negative heat flux to atm - expected = -4._r8 * mm_to_m * denh2o * hfus - @assertEqual(expected, sens_heat_flux, tolerance=tol) - - end subroutine partialConversion_resultsInCorrectHeatFlux - -end module test_sens_heat_from_precip_conversion - diff --git a/src/unit_test_shr/CMakeLists.txt b/src/unit_test_shr/CMakeLists.txt index 8e3494f8c7..d0172f5931 100644 --- a/src/unit_test_shr/CMakeLists.txt +++ b/src/unit_test_shr/CMakeLists.txt @@ -9,12 +9,13 @@ sourcelist_to_parent(clm_genf90_sources) list(APPEND clm_sources "${clm_genf90_sources}") list(APPEND clm_sources - unittestTimeManagerMod.F90 - unittestSubgridMod.F90 - unittestSimpleSubgridSetupsMod.F90 unittestFilterBuilderMod.F90 unittestGlcMec.F90 + unittestSimpleSubgridSetupsMod.F90 + unittestSubgridMod.F90 + unittestTimeManagerMod.F90 unittestUtils.F90 + unittestWaterTypeFactory.F90 ) sourcelist_to_parent(clm_sources) diff --git a/src/unit_test_shr/unittestWaterTypeFactory.F90 b/src/unit_test_shr/unittestWaterTypeFactory.F90 new file mode 100644 index 0000000000..4060c21dbc --- /dev/null +++ b/src/unit_test_shr/unittestWaterTypeFactory.F90 @@ -0,0 +1,142 @@ +module unittestWaterTypeFactory + + ! This module contains a class and associated methods that assist unit tests that need + ! to create a water_type instance + ! + ! This is object-oriented in case we want to remember initial values for the sake of + ! doing a clean teardown, or have any other memory between calls (currently we don't + ! make use of this, but we could in the future). + ! + ! Typical usage is: + ! + ! - Include an instance of this class in the class used for unit testing + ! + ! - Before calling any other methods: call init() + ! + ! - Before doing any subgrid setup: call setup_before_subgrid() + ! + ! - After doing any subgrid setup: call setup_after_subgrid() + ! + ! - Then get the water_type inst with: call create_water_type() + ! + ! - In the unit test tearDown method, before unittest_subgrid_teardown: call teardown() + + use shr_kind_mod , only : r8 => shr_kind_r8 + use clm_varpar, only : nlevsoi, nlevgrnd, nlevsno + use ColumnType, only : col + use WaterType, only : water_type, water_params_type + use unittestArrayMod, only : col_array + use unittestSubgridMod, only : bounds + + implicit none + private + + type, public :: unittest_water_type_factory_type + contains + procedure, public :: init + procedure, public :: setup_before_subgrid + procedure, public :: setup_after_subgrid + procedure, public :: create_water_type + procedure, public :: teardown + end type unittest_water_type_factory_type + +contains + + subroutine init(this) + ! Initialize this instance (like a constructor, but Fortran constructors can be + ! problematic / buggy, so we do it via an init subroutine instead) + class(unittest_water_type_factory_type), intent(inout) :: this + + ! For now, nothing to do: just a placeholder + end subroutine init + + subroutine setup_before_subgrid(this, my_nlevsoi, nlevgrnd_additional, my_nlevsno) + ! Do initializations that need to happen before setting up the subgrid structure. + ! + ! If you have already initialized any of these variables yourself, then simply pass + ! in the current values to make this a routine a no-op. + class(unittest_water_type_factory_type), intent(inout) :: this + integer, intent(in) :: my_nlevsoi + integer, intent(in) :: nlevgrnd_additional ! nlevgrnd = nlevsoi + nlevgrnd_additional + integer, intent(in) :: my_nlevsno + + nlevsoi = my_nlevsoi + nlevgrnd = nlevsoi + nlevgrnd_additional + nlevsno = my_nlevsno + end subroutine setup_before_subgrid + + subroutine setup_after_subgrid(this, snl, dz) + ! Do initializations that need to happen after setting up the subgrid structure + class(unittest_water_type_factory_type), intent(inout) :: this + + ! For now, set all snl and dz values to the single input + integer, intent(in) :: snl + real(r8), intent(in) :: dz + + col%snl(:) = snl + col%dz(:,:) = dz + end subroutine setup_after_subgrid + + subroutine create_water_type(this, water_inst, & + enable_consistency_checks, enable_isotopes) + ! Initialize water_inst + ! + ! Assumes that setup_before_subgrid and setup_after_subgrid have been called and that + ! subgrid setup has been done. + ! + ! Arbitrary values are used for some of the inputs to the water_inst init routine + ! + ! If enable_consistency_checks or enable_isotopes are missing, they are assumed to be + ! false. + class(unittest_water_type_factory_type), intent(in) :: this + type(water_type), intent(inout) :: water_inst + logical, intent(in), optional :: enable_consistency_checks + logical, intent(in), optional :: enable_isotopes + + logical :: l_enable_consistency_checks + logical :: l_enable_isotopes + type(water_params_type) :: params + real(r8), allocatable :: watsat_col(:,:) + real(r8), allocatable :: t_soisno_col(:,:) + + if (present(enable_consistency_checks)) then + l_enable_consistency_checks = enable_consistency_checks + else + l_enable_consistency_checks = .false. + end if + + if (present(enable_isotopes)) then + l_enable_isotopes = enable_isotopes + else + l_enable_isotopes = .false. + end if + + params = water_params_type( & + enable_consistency_checks = l_enable_consistency_checks, & + enable_isotopes = l_enable_isotopes) + + allocate(watsat_col(bounds%begc:bounds%endc, nlevgrnd)) + watsat_col(:,:) = 0._r8 + allocate(t_soisno_col(bounds%begc:bounds%endc, -nlevsno+1:nlevgrnd)) + t_soisno_col(:,:) = 275._r8 + + call water_inst%InitForTesting(bounds, params, & + h2osno_col = col_array(0._r8), & + snow_depth_col = col_array(0._r8), & + watsat_col = watsat_col, & + t_soisno_col = t_soisno_col) + end subroutine create_water_type + + subroutine teardown(this, water_inst) + ! Should be called from the unittest tearDown method, before unittest_subgrid_teardown + class(unittest_water_type_factory_type), intent(inout) :: this + type(water_type), intent(in) :: water_inst + + ! For now, nothing to do: just a placeholder. + ! + ! Ideally this would call water_inst%Clean. But we don't yet have a clean method in + ! water_type or most of the types it contains. So for now we have a small memory leak + ! in each test. + end subroutine teardown + +end module unittestWaterTypeFactory