module SoilHydrologyMod !----------------------------------------------------------------------- ! !DESCRIPTION: ! Calculate soil hydrology ! #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_vichydro use clm_varcon , only : ispval use clm_varcon , only : denh2o, denice, rpi use clm_varcon , only : pondmx_urban use clm_varpar , only : nlevsoi, nlevgrnd, nlayer, nlayert use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use column_varcon , only : icol_road_imperv use landunit_varcon , only : istsoil, istcrop use clm_time_manager , only : get_step_size_real use NumericsMod , only : truncate_small_values use EnergyFluxType , only : energyflux_type use InfiltrationExcessRunoffMod, only : infiltration_excess_runoff_type use SoilHydrologyType , only : soilhydrology_type use SoilStateType , only : soilstate_type use Wateratm2lndBulkType, only : wateratm2lndbulk_type use WaterFluxType , only : waterflux_type use WaterFluxBulkType , only : waterfluxbulk_type use WaterStateType , only : waterstate_type use WaterStateBulkType, only : waterstatebulk_type use WaterDiagnosticBulkType, only : waterdiagnosticbulk_type use TemperatureType , only : temperature_type use LandunitType , only : lun use ColumnType , only : column_type, col use PatchType , only : patch ! ! !PUBLIC TYPES: implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: SoilHydReadNML ! Read in the Soil hydrology namelist public :: SetSoilWaterFractions ! Set diagnostic variables related to the fraction of water and ice in each layer public :: SetFloodc ! Apply gridcell flood water flux to non-lake columns public :: SetQflxInputs ! Set the flux of water into the soil from the top public :: Infiltration ! Calculate total infiltration public :: TotalSurfaceRunoff ! Calculate total surface runoff public :: UpdateUrbanPonding ! Update the state variable representing ponding on urban surfaces public :: WaterTable ! Calculate water table before imposing drainage public :: Drainage ! Calculate subsurface drainage public :: CLMVICMap public :: PerchedWaterTable ! Calculate perched water table public :: PerchedLateralFlow ! Calculate lateral flow from perched saturated zone public :: PerchedLateralFlowHillslope ! Calculate lateral flow from perched saturated zone in hillslope configuration public :: ThetaBasedWaterTable ! Calculate water table from soil moisture state public :: LateralFlowPowerLaw ! Calculate lateral flow based on power law drainage function public :: LateralFlowHillslope ! Calculate lateral in multi-column hillslope configuration public :: RenewCondensation ! Misc. corrections public :: CalcIrrigWithdrawals ! Calculate irrigation withdrawals from groundwater by layer public :: WithdrawGroundwaterIrrigation ! Remove groundwater irrigation from unconfined and confined aquifers public :: readParams type, private :: params_type real(r8) :: aq_sp_yield_min ! Minimum aquifer specific yield (unitless) real(r8) :: n_baseflow ! Drainage power law exponent (unitless) real(r8) :: perched_baseflow_scalar ! Scalar multiplier for perched base flow rate (kg/m2/s) ! real(r8) :: e_ice ! Soil ice impedance factor (unitless) end type params_type type(params_type), private :: params_inst !----------------------------------------------------------------------- real(r8), private :: baseflow_scalar = 1.e-2_r8 real(r8), parameter :: tolerance = 1.e-12_r8 ! tolerance for checking whether sublimation is greater than ice in top soil layer character(len=*), parameter, private :: sourcefile = & __FILE__ contains !----------------------------------------------------------------------- subroutine readParams( ncid ) ! ! !USES: use ncdio_pio, only: file_desc_t use paramUtilMod, only: readNcdioScalar ! ! !ARGUMENTS: implicit none type(file_desc_t),intent(inout) :: ncid ! pio netCDF file id ! ! !LOCAL VARIABLES: character(len=*), parameter :: subname = 'readParams_SoilHydrology' !-------------------------------------------------------------------- ! Minimum aquifer specific yield (unitless) call readNcdioScalar(ncid, 'aq_sp_yield_min', subname, params_inst%aq_sp_yield_min) ! Drainage power law exponent (unitless) call readNcdioScalar(ncid, 'n_baseflow', subname, params_inst%n_baseflow) ! Scalar multiplier for perched base flow rate (kg/m2/s) call readNcdioScalar(ncid, 'perched_baseflow_scalar', subname, params_inst%perched_baseflow_scalar) ! ! Soil ice impedance factor (unitless) ! call readNcdioScalar(ncid, 'e_ice', subname, params_inst%e_ice) end subroutine readParams !----------------------------------------------------------------------- subroutine soilHydReadNML( NLFilename ) ! ! !DESCRIPTION: ! Read the namelist for soil hydrology ! ! !USES: use fileutils , only : getavu, relavu, opnfil use shr_nl_mod , only : shr_nl_find_group_name use spmdMod , only : masterproc, mpicom use shr_mpi_mod , only : shr_mpi_bcast use clm_varctl , only : iulog use shr_log_mod , only : errMsg => shr_log_errMsg ! ! !ARGUMENTS: character(len=*), intent(in) :: NLFilename ! Namelist filename ! ! !LOCAL VARIABLES: integer :: ierr ! error code integer :: unitn ! unit for namelist file character(len=*), parameter :: subname = 'soilHydReadNML' character(len=*), parameter :: nmlname = 'soilhydrology_inparm' !----------------------------------------------------------------------- namelist /soilhydrology_inparm/ baseflow_scalar ! Initialize options to default values, in case they are not specified in ! the namelist if (masterproc) then unitn = getavu() write(iulog,*) 'Read in '//nmlname//' namelist' call opnfil (NLFilename, unitn, 'F') call shr_nl_find_group_name(unitn, nmlname, status=ierr) if (ierr == 0) then read(unitn, nml=soilhydrology_inparm, iostat=ierr) if (ierr /= 0) then call endrun(msg="ERROR reading "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) end if else call endrun(msg="ERROR could NOT find "//nmlname//"namelist"//errmsg(sourcefile, __LINE__)) end if call relavu( unitn ) end if call shr_mpi_bcast (baseflow_scalar, mpicom) if (masterproc) then write(iulog,*) ' ' write(iulog,*) nmlname//' settings:' write(iulog,nml=soilhydrology_inparm) write(iulog,*) ' ' end if end subroutine soilhydReadNML !----------------------------------------------------------------------- subroutine SetSoilWaterFractions(bounds, num_hydrologyc, filter_hydrologyc, & soilhydrology_inst, soilstate_inst, waterstatebulk_inst) ! ! !DESCRIPTION: ! Set diagnostic variables related to the fraction of water and ice in each layer ! ! !USES: use clm_varcon, only : denice ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(inout) :: soilstate_inst type(waterstatebulk_type), intent(in) :: waterstatebulk_inst ! ! !LOCAL VARIABLES: integer :: j, fc, c real(r8) :: vol_ice(bounds%begc:bounds%endc,1:nlevsoi) !partial volume of ice lens in layer real(r8) :: icefrac_orig ! original formulation for icefrac character(len=*), parameter :: subname = 'SetSoilWaterFractions' !----------------------------------------------------------------------- associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Output: [real(r8) (:,:) ] effective porosity = porosity - vol_ice h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice water (kg/m2) origflag => soilhydrology_inst%origflag , & ! Input: logical icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fracice => soilhydrology_inst%fracice_col & ! Output: [real(r8) (:,:) ] fractional impermeability (-) ) do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Porosity of soil, partial volume of ice and liquid, fraction of ice in each layer, ! fractional impermeability vol_ice(c,j) = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) eff_porosity(c,j) = max(0.01_r8,watsat(c,j)-vol_ice(c,j)) icefrac(c,j) = min(1._r8,vol_ice(c,j)/watsat(c,j)) ! fracice is only used in code with origflag == 1. For this calculation, we use ! the version of icefrac that was used in this original hydrology code. if (h2osoi_ice(c,j) == 0._r8) then ! Avoid possible divide by zero (in case h2osoi_liq(c,j) is also 0) icefrac_orig = 0._r8 else icefrac_orig = min(1._r8,h2osoi_ice(c,j)/(h2osoi_ice(c,j)+h2osoi_liq(c,j))) end if fracice(c,j) = max(0._r8,exp(-3._r8*(1._r8-icefrac_orig))- exp(-3._r8))/(1.0_r8-exp(-3._r8)) end do end do end associate end subroutine SetSoilWaterFractions !----------------------------------------------------------------------- subroutine SetFloodc(num_nolakec, filter_nolakec, col, & wateratm2lndbulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Apply gridcell flood water flux to non-lake columns ! ! !ARGUMENTS: 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(column_type) , intent(in) :: col type(wateratm2lndbulk_type) , intent(in) :: wateratm2lndbulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: fc, c, g character(len=*), parameter :: subname = 'SetFloodc' !----------------------------------------------------------------------- associate( & qflx_floodg => wateratm2lndbulk_inst%forc_flood_grc , & ! Input: [real(r8) (:) ] gridcell flux of flood water from RTM qflx_floodc => waterfluxbulk_inst%qflx_floodc_col & ! Output: [real(r8) (:) ] column flux of flood water from RTM ) do fc = 1, num_nolakec c = filter_nolakec(fc) g = col%gridcell(c) if (col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall) then qflx_floodc(c) = qflx_floodg(g) else qflx_floodc(c) = 0._r8 end if end do end associate end subroutine SetFloodc !----------------------------------------------------------------------- subroutine SetQflxInputs(bounds, num_hydrologyc, filter_hydrologyc, & waterfluxbulk_inst, waterdiagnosticbulk_inst) ! ! !DESCRIPTION: ! Set various input fluxes of water ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst type(waterdiagnosticbulk_type) , intent(in) :: waterdiagnosticbulk_inst ! ! !LOCAL VARIABLES: integer :: fc, c real(r8) :: qflx_evap ! evaporation for this column real(r8) :: fsno ! copy of frac_sno character(len=*), parameter :: subname = 'SetQflxInputs' !----------------------------------------------------------------------- associate( & snl => col%snl , & ! Input: [integer (:) ] minus number of snow layers qflx_top_soil => waterfluxbulk_inst%qflx_top_soil_col , & ! Output: [real(r8) (:)] net water input into soil from top (mm/s) qflx_in_soil => waterfluxbulk_inst%qflx_in_soil_col , & ! Output: [real(r8) (:)] surface input to soil (mm/s) qflx_top_soil_to_h2osfc => waterfluxbulk_inst%qflx_top_soil_to_h2osfc_col , & ! Output: [real(r8) (:)] portion of qflx_top_soil going to h2osfc, minus evaporation (mm/s) qflx_rain_plus_snomelt => waterfluxbulk_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:)] rain plus snow melt falling on the soil (mm/s) qflx_snow_h2osfc => waterfluxbulk_inst%qflx_snow_h2osfc_col , & ! Input: [real(r8) (:)] snow falling on surface water (mm/s) qflx_floodc => waterfluxbulk_inst%qflx_floodc_col , & ! Input: [real(r8) (:)] column flux of flood water from RTM qflx_ev_soil => waterfluxbulk_inst%qflx_ev_soil_col , & ! Input: [real(r8) (:)] evaporation flux from soil (W/m**2) [+ to atm] qflx_liqevap_from_top_layer => waterfluxbulk_inst%qflx_liqevap_from_top_layer_col, & ! Input: [real(r8) (:)] rate of liquid water evaporated from top soil or snow layer (mm H2O/s) [+] qflx_ev_h2osfc => waterfluxbulk_inst%qflx_ev_h2osfc_col , & ! Input: [real(r8) (:)] evaporation flux from h2osfc (W/m**2) [+ to atm] qflx_sat_excess_surf => waterfluxbulk_inst%qflx_sat_excess_surf_col , & ! Input: [real(r8) (:)] surface runoff due to saturated surface (mm H2O /s) frac_sno => waterdiagnosticbulk_inst%frac_sno_eff_col , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1) frac_h2osfc => waterdiagnosticbulk_inst%frac_h2osfc_col & ! Input: [real(r8) (:) ] fraction of ground covered by surface water (0 to 1) ) do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_top_soil(c) = qflx_rain_plus_snomelt(c) + qflx_snow_h2osfc(c) + qflx_floodc(c) ! ------------------------------------------------------------------------ ! Partition surface inputs between soil and h2osfc ! ------------------------------------------------------------------------ if (snl(c) >= 0) then fsno=0._r8 ! if no snow layers, sublimation is removed from h2osoi_ice in drainage qflx_evap=qflx_liqevap_from_top_layer(c) else fsno=frac_sno(c) qflx_evap=qflx_ev_soil(c) endif qflx_in_soil(c) = (1._r8 - frac_h2osfc(c)) * (qflx_top_soil(c) - qflx_sat_excess_surf(c)) qflx_top_soil_to_h2osfc(c) = frac_h2osfc(c) * (qflx_top_soil(c) - qflx_sat_excess_surf(c)) ! remove evaporation (snow treated in SnowHydrology) qflx_in_soil(c) = qflx_in_soil(c) - (1.0_r8 - fsno - frac_h2osfc(c))*qflx_evap qflx_top_soil_to_h2osfc(c) = qflx_top_soil_to_h2osfc(c) - frac_h2osfc(c) * qflx_ev_h2osfc(c) end do end associate end subroutine SetQflxInputs !----------------------------------------------------------------------- subroutine RouteInfiltrationExcess(bounds, num_hydrologyc, filter_hydrologyc, & waterfluxbulk_inst, soilhydrology_inst) ! ! !DESCRIPTION: ! Route the infiltration excess runoff flux ! ! This adjusts infiltration, input to h2osfc, and the associated surface runoff flux. ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst type(soilhydrology_type) , intent(in) :: soilhydrology_inst ! ! !LOCAL VARIABLES: integer :: fc, c character(len=*), parameter :: subname = 'RouteInfiltrationExcess' !----------------------------------------------------------------------- associate( & qflx_in_soil_limited => waterfluxbulk_inst%qflx_in_soil_limited_col , & ! Output: [real(r8) (:) ] surface input to soil, limited by max infiltration rate (mm H2O /s) qflx_in_h2osfc => waterfluxbulk_inst%qflx_in_h2osfc_col , & ! Output: [real(r8) (:) ] total surface input to h2osfc qflx_infl_excess_surf => waterfluxbulk_inst%qflx_infl_excess_surf_col , & ! Output: [real(r8) (:) ] surface runoff due to infiltration excess (mm H2O /s) qflx_in_soil => waterfluxbulk_inst%qflx_in_soil_col , & ! Input: [real(r8) (:) ] surface input to soil (mm/s) qflx_top_soil_to_h2osfc => waterfluxbulk_inst%qflx_top_soil_to_h2osfc_col , & ! Input: [real(r8) (:) ] portion of qflx_top_soil going to h2osfc, minus evaporation (mm/s) qflx_infl_excess => waterfluxbulk_inst%qflx_infl_excess_col , & ! Input: [real(r8) (:) ] infiltration excess runoff (mm H2O /s) h2osfcflag => soilhydrology_inst%h2osfcflag & ! Input: integer ) do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (lun%itype(col%landunit(c)) == istsoil .or. lun%itype(col%landunit(c))==istcrop) then qflx_in_soil_limited(c) = qflx_in_soil(c) - qflx_infl_excess(c) if (h2osfcflag /= 0) then qflx_in_h2osfc(c) = qflx_top_soil_to_h2osfc(c) + qflx_infl_excess(c) qflx_infl_excess_surf(c) = 0._r8 else ! No h2osfc pool, so qflx_infl_excess goes directly to surface runoff qflx_in_h2osfc(c) = qflx_top_soil_to_h2osfc(c) qflx_infl_excess_surf(c) = qflx_infl_excess(c) end if else ! non-vegetated landunits (i.e. urban) use original CLM4 code qflx_in_soil_limited(c) = qflx_in_soil(c) qflx_in_h2osfc(c) = 0._r8 qflx_infl_excess_surf(c) = 0._r8 end if end do end associate end subroutine RouteInfiltrationExcess !----------------------------------------------------------------------- subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, & waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate total infiltration ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: fc, c character(len=*), parameter :: subname = 'Infiltration' !----------------------------------------------------------------------- associate( & qflx_infl => waterfluxbulk_inst%qflx_infl_col , & ! Output: [real(r8) (:) ] infiltration (mm H2O /s) qflx_in_soil_limited => waterfluxbulk_inst%qflx_in_soil_limited_col , & ! Input: [real(r8) (:) ] surface input to soil, limited by max infiltration rate (mm H2O /s) qflx_h2osfc_drain => waterfluxbulk_inst%qflx_h2osfc_drain_col & ! Input: [real(r8) (:) ] bottom drainage from h2osfc (mm H2O /s) ) do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_infl(c) = qflx_in_soil_limited(c) + qflx_h2osfc_drain(c) end do end associate end subroutine Infiltration !----------------------------------------------------------------------- subroutine TotalSurfaceRunoff(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc, & waterfluxbulk_inst, soilhydrology_inst, & waterstatebulk_inst) ! ! !DESCRIPTION: ! Calculate total surface runoff ! ! For hydrologically-active columns, this is the sum of some already-computed terms ! ! In addition, computes total surface runoff for non-hydrologically-active urban ! columns ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type), intent(in) :: waterstatebulk_inst ! ! !LOCAL VARIABLES: integer :: fc, c real(r8) :: dtime ! land model time step (sec) character(len=*), parameter :: subname = 'TotalSurfaceRunoff' !----------------------------------------------------------------------- associate( & snl => col%snl , & ! Input: [integer (:) ] minus number of snow layers qflx_surf => waterfluxbulk_inst%qflx_surf_col , & ! Output: [real(r8) (:) ] total surface runoff (mm H2O /s) qflx_infl_excess_surf => waterfluxbulk_inst%qflx_infl_excess_surf_col, & ! Input: [real(r8) (:) ] surface runoff due to infiltration excess (mm H2O /s) qflx_h2osfc_surf => waterfluxbulk_inst%qflx_h2osfc_surf_col, & ! Input: [real(r8) (:) ] surface water runoff (mm H2O /s) qflx_rain_plus_snomelt => waterfluxbulk_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:) ] rain plus snow melt falling on the soil (mm/s) qflx_liqevap_from_top_layer => waterfluxbulk_inst%qflx_liqevap_from_top_layer_col, & ! Input: [real(r8) (:) ] rate of liquid water evaporated from top soil or snow layer (mm H2O/s) [+] qflx_floodc => waterfluxbulk_inst%qflx_floodc_col , & ! Input: [real(r8) (:) ] column flux of flood water from RTM qflx_sat_excess_surf => waterfluxbulk_inst%qflx_sat_excess_surf_col , & ! Input: [real(r8) (:) ] surface runoff due to saturated surface (mm H2O /s) xs_urban => soilhydrology_inst%xs_urban_col , & ! Output: [real(r8) (:) ] excess soil water above urban ponding limit h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) ) dtime = get_step_size_real() ! ------------------------------------------------------------------------ ! Set qflx_surf for hydrologically-active columns ! ------------------------------------------------------------------------ do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Depending on whether h2osfcflag is 0 or 1, one of qflx_infl_excess or ! qflx_h2osfc_surf will always be 0. But it's safe to just add them both. qflx_surf(c) = qflx_sat_excess_surf(c) + qflx_infl_excess_surf(c) + qflx_h2osfc_surf(c) end do ! ------------------------------------------------------------------------ ! Set qflx_surf for non-hydrologically-active urban columns ! ! Determine water in excess of ponding limit for urban roof and impervious road. ! Excess goes to surface runoff. No surface runoff for sunwall and shadewall. ! ------------------------------------------------------------------------ do fc = 1, num_urbanc c = filter_urbanc(fc) if (col%itype(c) == icol_roof .or. col%itype(c) == icol_road_imperv) then ! If there are snow layers then all qflx_rain_plus_snomelt goes to surface runoff if (snl(c) < 0) then qflx_surf(c) = max(0._r8,qflx_rain_plus_snomelt(c)) else ! NOTE(wjs, 2017-07-11) It looks to me like this use of xs_urban and ! h2osoi_liq(c,1) are roughly analogous to h2osfc in hydrologically-active ! columns. Why not use h2osfc for urban columns, too? xs_urban(c) = max(0._r8, & h2osoi_liq(c,1)/dtime + qflx_rain_plus_snomelt(c) - qflx_liqevap_from_top_layer(c) - & pondmx_urban/dtime) qflx_surf(c) = xs_urban(c) end if ! send flood water flux to runoff for all urban columns qflx_surf(c) = qflx_surf(c) + qflx_floodc(c) else if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall) then qflx_surf(c) = 0._r8 end if end do end associate end subroutine TotalSurfaceRunoff !----------------------------------------------------------------------- subroutine UpdateUrbanPonding(bounds, num_urbanc, filter_urbanc, & waterstatebulk_inst, soilhydrology_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Update the state variable representing ponding on urban surfaces ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(soilhydrology_type) , intent(in) :: soilhydrology_inst type(waterfluxbulk_type) , intent(in) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: fc, c real(r8) :: dtime ! land model time step (sec) character(len=*), parameter :: subname = 'UpdateUrbanPonding' !----------------------------------------------------------------------- associate( & snl => col%snl , & ! Input: [integer (:) ] minus number of snow layers h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) xs_urban => soilhydrology_inst%xs_urban_col , & ! Input: [real(r8) (:) ] excess soil water above urban ponding limit qflx_rain_plus_snomelt => waterfluxbulk_inst%qflx_rain_plus_snomelt_col , & ! Input: [real(r8) (:) ] rain plus snow melt falling on the soil (mm/s) qflx_liqevap_from_top_layer => waterfluxbulk_inst%qflx_liqevap_from_top_layer_col & ! Input: [real(r8) (:) ] rate of liquid water evaporated from top soil or snow layer (mm H2O/s) [+] ) dtime = get_step_size_real() do fc = 1, num_urbanc c = filter_urbanc(fc) if (col%itype(c) == icol_roof .or. col%itype(c) == icol_road_imperv) then if (snl(c) >= 0) then if (xs_urban(c) > 0.) then h2osoi_liq(c,1) = pondmx_urban else h2osoi_liq(c,1) = max(0._r8,h2osoi_liq(c,1)+ & (qflx_rain_plus_snomelt(c)-qflx_liqevap_from_top_layer(c))*dtime) end if end if end if end do end associate end subroutine UpdateUrbanPonding !----------------------------------------------------------------------- subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, & soilhydrology_inst, soilstate_inst, temperature_inst, waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate watertable, considering aquifer recharge but no drainage. ! ! !USES: use clm_varcon , only : pondmx, tfrz, watmin,denice,denh2o use column_varcon , only : icol_roof, icol_road_imperv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(in) :: soilstate_inst type(temperature_type) , intent(in) :: temperature_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: xs(bounds%begc:bounds%endc) ! water needed to bring soil moisture to watmin (mm) real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi) ! layer thickness (mm) integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-) real(r8) :: rsub_top(bounds%begc:bounds%endc) ! subsurface runoff - topographic control (mm/s) real(r8) :: xsi(bounds%begc:bounds%endc) ! excess soil water above saturation at layer i (mm) real(r8) :: rous ! aquifer yield (-) real(r8) :: wh ! smpfz(jwt)-z(jwt) (mm) real(r8) :: ws ! summation of pore space of layers below water table (mm) real(r8) :: s_node ! soil wetness (-) real(r8) :: dzsum ! summation of dzmm of layers below water table (mm) real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-) real(r8) :: fracice_rsub(bounds%begc:bounds%endc) ! fractional impermeability of soil layers (-) real(r8) :: ka ! hydraulic conductivity of the aquifer (mm/s) real(r8) :: available_h2osoi_liq ! available soil liquid water in a layer real(r8) :: imped real(r8) :: rsub_top_tot real(r8) :: rsub_top_layer real(r8) :: qcharge_tot real(r8) :: qcharge_layer real(r8) :: theta_unsat real(r8) :: f_unsat real(r8) :: s_y integer :: k,k_frz,k_perch real(r8) :: sat_lev real(r8) :: s1 real(r8) :: s2 real(r8) :: m real(r8) :: b real(r8) :: q_perch real(r8) :: q_perch_max real(r8) :: dflag=0._r8 real(r8) :: qflx_solidevap_from_top_layer_save ! temporary !----------------------------------------------------------------------- associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] qflx_ev_snow => waterfluxbulk_inst%qflx_ev_snow_col , & ! In/Out: [real(r8) (:) ] evaporation flux from snow (mm H2O/s) [+ to atm] bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice zwt => soilhydrology_inst%zwt_col , & ! Output: [real(r8) (:) ] water table depth (m) zwt_perched => soilhydrology_inst%zwt_perched_col , & ! Output: [real(r8) (:) ] perched water table depth (m) frost_table => soilhydrology_inst%frost_table_col , & ! Output: [real(r8) (:) ] frost table depth (m) wa => waterstatebulk_inst%wa_col , & ! Output: [real(r8) (:) ] water in the unconfined aquifer (mm) qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s) origflag => soilhydrology_inst%origflag , & ! Input: logical qflx_drain => waterfluxbulk_inst%qflx_drain_col , & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s) qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col , & ! Output: [real(r8) (:) ] perched wt sub-surface runoff (mm H2O /s) qflx_rsub_sat => waterfluxbulk_inst%qflx_rsub_sat_col & ! Output: [real(r8) (:) ] soil saturation excess [mm h2o/s] ) ! Get time step dtime = get_step_size_real() ! Convert layer thicknesses from m to mm do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 end do end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain(c) = 0._r8 qflx_rsub_sat(c) = 0._r8 qflx_drain_perched(c) = 0._r8 end do ! The layer index of the first unsaturated layer, i.e., the layer right above ! the water table do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) jwt(c) = nlevsoi ! allow jwt to equal zero when zwt is in top layer do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo end do !============================== QCHARGE ========================================= ! Water table changes due to qcharge do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) rous=max(rous, params_inst%aq_sp_yield_min) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then wa(c) = wa(c) + qcharge(c) * dtime zwt(c) = zwt(c) - (qcharge(c) * dtime)/1000._r8/rous else !-- water table within soil layers 1-9 ------------------------------------- ! try to raise water table to account for qcharge qcharge_tot = qcharge(c) * dtime if(qcharge_tot > 0.) then !rising water table do j = jwt(c)+1, 1,-1 ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y, params_inst%aq_sp_yield_min) qcharge_layer=min(qcharge_tot,(s_y*(zwt(c) - zi(c,j-1))*1.e3)) qcharge_layer=max(qcharge_layer,0._r8) if(s_y > 0._r8) zwt(c) = zwt(c) - qcharge_layer/s_y/1000._r8 qcharge_tot = qcharge_tot - qcharge_layer if (qcharge_tot <= 0.) exit enddo else ! deepening water table (negative qcharge) do j = jwt(c)+1, nlevsoi ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y, params_inst%aq_sp_yield_min) qcharge_layer=max(qcharge_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) qcharge_layer=min(qcharge_layer,0._r8) qcharge_tot = qcharge_tot - qcharge_layer if (qcharge_tot >= 0.) then zwt(c) = zwt(c) - qcharge_layer/s_y/1000._r8 exit else zwt(c) = zi(c,j) endif enddo if (qcharge_tot > 0.) zwt(c) = zwt(c) - qcharge_tot/1000._r8/rous endif !-- recompute jwt for following calculations --------------------------------- ! allow jwt to equal zero when zwt is in top layer jwt(c) = nlevsoi do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo endif enddo !== BASEFLOW ================================================== ! perched water table code do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! define frost table as first frozen layer with unfrozen layer above it if(t_soisno(c,1) > tfrz) then k_frz=nlevsoi else k_frz=1 endif do k=2, nlevsoi if (t_soisno(c,k-1) > tfrz .and. t_soisno(c,k) <= tfrz) then k_frz=k exit endif enddo frost_table(c)=z(c,k_frz) ! initialize perched water table to frost table, and qflx_drain_perched(c) to zero zwt_perched(c)=frost_table(c) !=================== water table above frost table ============================= ! if water table is above frost table, do not use topmodel baseflow formulation if (zwt(c) < frost_table(c) .and. t_soisno(c,k_frz) <= tfrz & .and. origflag == 0) then else !=================== water table below frost table ============================= !-- compute possible perched water table *and* groundwater table afterwards ! locate perched water table from bottom up starting at frost table ! sat_lev is an arbitrary saturation level used to determine perched water table sat_lev=0.9 k_perch=1 do k=k_frz,1,-1 h2osoi_vol(c,k) = h2osoi_liq(c,k)/(dz(c,k)*denh2o) & + h2osoi_ice(c,k)/(dz(c,k)*denice) if (h2osoi_vol(c,k)/watsat(c,k) <= sat_lev) then k_perch=k exit endif enddo ! if frost_table = nlevsoi, only compute perched water table if frozen if (t_soisno(c,k_frz) > tfrz) k_perch=k_frz ! if perched water table exists if (k_frz > k_perch) then ! interpolate between k_perch and k_perch+1 to find perched water table height s1 = (h2osoi_liq(c,k_perch)/(dz(c,k_perch)*denh2o) & + h2osoi_ice(c,k_perch)/(dz(c,k_perch)*denice))/watsat(c,k_perch) s2 = (h2osoi_liq(c,k_perch+1)/(dz(c,k_perch+1)*denh2o) & + h2osoi_ice(c,k_perch+1)/(dz(c,k_perch+1)*denice))/watsat(c,k_perch+1) m=(z(c,k_perch+1)-z(c,k_perch))/(s2-s1) b=z(c,k_perch+1)-m*s2 zwt_perched(c)=max(0._r8,m*sat_lev+b) endif !k_frz > k_perch endif end do end associate end subroutine WaterTable !----------------------------------------------------------------------- subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filter_urbanc, & temperature_inst, soilhydrology_inst, soilstate_inst, waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate subsurface drainage ! ! !USES: use clm_varcon , only : pondmx, tfrz, watmin,rpi, secspday, nlvic use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(temperature_type) , intent(in) :: temperature_inst type(soilstate_type) , intent(in) :: soilstate_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'Drainage' ! subroutine name integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: xs(bounds%begc:bounds%endc) ! water needed to bring soil moisture to watmin (mm) real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi) ! layer thickness (mm) integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-) real(r8) :: rsub_top(bounds%begc:bounds%endc) ! subsurface runoff - topographic control (mm/s) real(r8) :: fff(bounds%begc:bounds%endc) ! decay factor (m-1) real(r8) :: xsi(bounds%begc:bounds%endc) ! excess soil water above saturation at layer i (mm) real(r8) :: xsia(bounds%begc:bounds%endc) ! available pore space at layer i (mm) real(r8) :: xs1(bounds%begc:bounds%endc) ! excess soil water above saturation at layer 1 (mm) real(r8) :: smpfz(1:nlevsoi) ! matric potential of layer right above water table (mm) real(r8) :: wtsub ! summation of hk*dzmm for layers below water table (mm**2/s) real(r8) :: rous ! aquifer yield (-) real(r8) :: wh ! smpfz(jwt)-z(jwt) (mm) real(r8) :: wh_zwt ! water head at the water table depth (mm) real(r8) :: ws ! summation of pore space of layers below water table (mm) real(r8) :: s_node ! soil wetness (-) real(r8) :: dzsum ! summation of dzmm of layers below water table (mm) real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-) real(r8) :: fracice_rsub(bounds%begc:bounds%endc) ! fractional impermeability of soil layers (-) real(r8) :: ka ! hydraulic conductivity of the aquifer (mm/s) real(r8) :: dza ! fff*(zwt-z(jwt)) (-) real(r8) :: available_h2osoi_liq ! available soil liquid water in a layer real(r8) :: rsub_top_max real(r8) :: h2osoi_vol real(r8) :: imped real(r8) :: rsub_top_tot real(r8) :: rsub_top_layer real(r8) :: qcharge_tot real(r8) :: qcharge_layer real(r8) :: theta_unsat real(r8) :: f_unsat real(r8) :: s_y integer :: k,k_frz,k_perch real(r8) :: sat_lev real(r8) :: s1 real(r8) :: s2 real(r8) :: m real(r8) :: b real(r8) :: q_perch real(r8) :: q_perch_max real(r8) :: vol_ice real(r8) :: dsmax_tmp(bounds%begc:bounds%endc) ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rsub_tmp ! temporary variable for ARNO subsurface runoff calculation real(r8) :: frac ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rel_moist ! relative moisture, temporary variable real(r8) :: wtsub_vic ! summation of hk*dzmm for layers in the third VIC layer !----------------------------------------------------------------------- associate( & z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) snl => col%snl , & ! Input: [integer (:) ] number of snow layers t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s) depth => soilhydrology_inst%depth_col , & ! Input: [real(r8) (:,:) ] VIC soil depth c_param => soilhydrology_inst%c_param_col , & ! Input: [real(r8) (:) ] baseflow exponent (Qb) Dsmax => soilhydrology_inst%dsmax_col , & ! Input: [real(r8) (:) ] max. velocity of baseflow (mm/day) max_moist => soilhydrology_inst%max_moist_col , & ! Input: [real(r8) (:,:) ] maximum soil moisture (ice + liq) moist => soilhydrology_inst%moist_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) Ds => soilhydrology_inst%ds_col , & ! Input: [real(r8) (:) ] fracton of Dsmax where non-linear baseflow begins Wsvic => soilhydrology_inst%Wsvic_col , & ! Input: [real(r8) (:) ] fraction of maximum soil moisutre where non-liear base flow occurs icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer hkdepth => soilhydrology_inst%hkdepth_col , & ! Input: [real(r8) (:) ] decay factor (m) frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) zwt_perched => soilhydrology_inst%zwt_perched_col , & ! Input: [real(r8) (:) ] perched water table depth (m) aquifer_water_baseline => waterstatebulk_inst%aquifer_water_baseline, & ! Input: [real(r8)] baseline value for water in the unconfined aquifer (wa_col) (mm) wa => waterstatebulk_inst%wa_col , & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm) ice => soilhydrology_inst%ice_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s) origflag => soilhydrology_inst%origflag , & ! Input: logical h2osfcflag => soilhydrology_inst%h2osfcflag , & ! Input: integer qflx_snwcp_liq => waterfluxbulk_inst%qflx_snwcp_liq_col , & ! Output: [real(r8) (:) ] excess liquid h2o due to snow capping (outgoing) (mm H2O /s) [+] qflx_ice_runoff_xs => waterfluxbulk_inst%qflx_ice_runoff_xs_col , & ! Output: [real(r8) (:) ] solid runoff from excess ice in soil (mm H2O /s) [+] qflx_liqdew_to_top_layer => waterfluxbulk_inst%qflx_liqdew_to_top_layer_col , & ! Output: [real(r8) (:) ] rate of liquid water deposited on top soil or snow layer (dew) (mm H2O /s) [+] qflx_soliddew_to_top_layer => waterfluxbulk_inst%qflx_soliddew_to_top_layer_col , & ! Output: [real(r8) (:) ] rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+] qflx_solidevap_from_top_layer => waterfluxbulk_inst%qflx_solidevap_from_top_layer_col, & ! Output: [real(r8) (:) ] rate of ice evaporated from top soil or snow layer (sublimation) (mm H2O /s) [+] qflx_drain => waterfluxbulk_inst%qflx_drain_col , & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s) qflx_qrgwl => waterfluxbulk_inst%qflx_qrgwl_col , & ! Output: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes (mm H2O /s) qflx_rsub_sat => waterfluxbulk_inst%qflx_rsub_sat_col , & ! Output: [real(r8) (:) ] soil saturation excess [mm h2o/s] qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col , & ! Output: [real(r8) (:) ] perched wt sub-surface runoff (mm H2O /s) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) ) ! Get time step dtime = get_step_size_real() ! Convert layer thicknesses from m to mm do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) end do end do ! Initial set do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain(c) = 0._r8 qflx_rsub_sat(c) = 0._r8 rsub_top(c) = 0._r8 fracice_rsub(c) = 0._r8 end do ! The layer index of the first unsaturated layer, i.e., the layer right above ! the water table do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) jwt(c) = nlevsoi ! allow jwt to equal zero when zwt is in top layer do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo end do rous = 0.2_r8 !== BASEFLOW ================================================== ! perched water table code do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! specify maximum drainage rate q_perch_max = params_inst%perched_baseflow_scalar * sin(col%topo_slope(c) * (rpi/180._r8)) ! if layer containing water table is frozen, compute the following: ! frost table, perched water table, and drainage from perched saturated layer ! define frost table as first frozen layer with unfrozen layer above it if(t_soisno(c,1) > tfrz) then k_frz=nlevsoi else k_frz=1 endif do k=2, nlevsoi if (t_soisno(c,k-1) > tfrz .and. t_soisno(c,k) <= tfrz) then k_frz=k exit endif enddo frost_table(c)=z(c,k_frz) ! initialize perched water table to frost table, and qflx_drain_perched(c) to zero zwt_perched(c)=frost_table(c) qflx_drain_perched(c) = 0._r8 !=================== water table above frost table ============================= ! if water table is above frost table, do not use topmodel baseflow formulation if (zwt(c) < frost_table(c) .and. t_soisno(c,k_frz) <= tfrz & .and. origflag == 0) then ! compute drainage from perched saturated region wtsub = 0._r8 q_perch = 0._r8 do k = jwt(c)+1, k_frz imped=10._r8**(-col%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do if (wtsub > 0._r8) q_perch = q_perch/wtsub qflx_drain_perched(c) = q_perch_max * q_perch & *(frost_table(c) - zwt(c)) ! remove drainage from perched saturated layers rsub_top_tot = - qflx_drain_perched(c) * dtime do k = jwt(c)+1, k_frz rsub_top_layer=max(rsub_top_tot,-(h2osoi_liq(c,k)-watmin)) rsub_top_layer=min(rsub_top_layer,0._r8) rsub_top_tot = rsub_top_tot - rsub_top_layer h2osoi_liq(c,k) = h2osoi_liq(c,k) + rsub_top_layer if (rsub_top_tot >= 0.) then zwt(c) = zwt(c) - rsub_top_layer/eff_porosity(c,k)/1000._r8 exit else zwt(c) = zi(c,k) endif enddo ! if rsub_top_tot is greater than available water (above frost table), ! then decrease qflx_drain_perched by residual amount for water balance qflx_drain_perched(c) = qflx_drain_perched(c) + rsub_top_tot/dtime !-- recompute jwt --------------------------------------------------------- ! allow jwt to equal zero when zwt is in top layer jwt(c) = nlevsoi do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo else !=================== water table below frost table ============================= !-- compute possible perched water table *and* groundwater table afterwards ! locate perched water table from bottom up starting at frost table ! sat_lev is an arbitrary saturation level used to determine perched water table sat_lev=0.9 k_perch=1 do k=k_frz,1,-1 h2osoi_vol = h2osoi_liq(c,k)/(dz(c,k)*denh2o) & + h2osoi_ice(c,k)/(dz(c,k)*denice) if (h2osoi_vol/watsat(c,k) <= sat_lev) then k_perch=k exit endif enddo ! if frost_table = nlevsoi, only compute perched water table if frozen if (t_soisno(c,k_frz) > tfrz) k_perch=k_frz ! if perched water table exists if (k_frz > k_perch) then ! interpolate between k_perch and k_perch+1 to find perched water table height s1 = (h2osoi_liq(c,k_perch)/(dz(c,k_perch)*denh2o) & + h2osoi_ice(c,k_perch)/(dz(c,k_perch)*denice))/watsat(c,k_perch) s2 = (h2osoi_liq(c,k_perch+1)/(dz(c,k_perch+1)*denh2o) & + h2osoi_ice(c,k_perch+1)/(dz(c,k_perch+1)*denice))/watsat(c,k_perch+1) m=(z(c,k_perch+1)-z(c,k_perch))/(s2-s1) b=z(c,k_perch+1)-m*s2 zwt_perched(c)=max(0._r8,m*sat_lev+b) ! compute drainage from perched saturated region wtsub = 0._r8 q_perch = 0._r8 do k = k_perch, k_frz imped=10._r8**(-col%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do if (wtsub > 0._r8) q_perch = q_perch/wtsub qflx_drain_perched(c) = q_perch_max * q_perch & *(frost_table(c) - zwt_perched(c)) ! no perched water table drainage if using original formulation if(origflag == 1) qflx_drain_perched(c) = 0._r8 ! remove drainage from perched saturated layers rsub_top_tot = - qflx_drain_perched(c) * dtime do k = k_perch+1, k_frz rsub_top_layer=max(rsub_top_tot,-(h2osoi_liq(c,k)-watmin)) rsub_top_layer=min(rsub_top_layer,0._r8) rsub_top_tot = rsub_top_tot - rsub_top_layer h2osoi_liq(c,k) = h2osoi_liq(c,k) + rsub_top_layer if (rsub_top_tot >= 0.) then zwt_perched(c) = zwt_perched(c) - rsub_top_layer/eff_porosity(c,k)/1000._r8 exit else zwt_perched(c) = zi(c,k) endif enddo ! if rsub_top_tot is greater than available water (above frost table), ! then decrease qflx_drain_perched by residual amount for water balance qflx_drain_perched(c) = qflx_drain_perched(c) + rsub_top_tot/dtime else qflx_drain_perched(c) = 0._r8 endif !k_frz > k_perch !-- Topographic runoff ---------------------------------------------------------------------- fff(c) = 1._r8/ hkdepth(c) dzsum = 0._r8 icefracsum = 0._r8 do j = max(jwt(c),1), nlevsoi dzsum = dzsum + dzmm(c,j) icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j) end do ! add ice impedance factor to baseflow if(origflag == 1) then if (use_vichydro) then call endrun(msg="VICHYDRO is not available for origflag=1"//errmsg(sourcefile, __LINE__)) else fracice_rsub(c) = max(0._r8,exp(-3._r8*(1._r8-(icefracsum/dzsum))) & - exp(-3._r8))/(1.0_r8-exp(-3._r8)) imped=(1._r8 - fracice_rsub(c)) rsub_top_max = 5.5e-3_r8 end if else if (use_vichydro) then imped=10._r8**(-col%e_ice(c)*min(1.0_r8,ice(c,nlayer)/max_moist(c,nlayer))) dsmax_tmp(c) = Dsmax(c) * dtime/ secspday !mm/day->mm/dtime rsub_top_max = dsmax_tmp(c) else imped=10._r8**(-col%e_ice(c)*(icefracsum/dzsum)) rsub_top_max = 10._r8 * sin((rpi/180.) * col%topo_slope(c)) end if endif if (use_vichydro) then ! ARNO model for the bottom soil layer (based on bottom soil layer ! moisture from previous time step ! use watmin instead for resid_moist to be consistent with default hydrology rel_moist = (moist(c,nlayer) - watmin)/(max_moist(c,nlayer)-watmin) frac = (Ds(c) * rsub_top_max )/Wsvic(c) rsub_tmp = (frac * rel_moist)/dtime if(rel_moist > Wsvic(c))then frac = (rel_moist - Wsvic(c))/(1.0_r8 - Wsvic(c)) rsub_tmp = rsub_tmp + (rsub_top_max * (1.0_r8 - Ds(c)/Wsvic(c)) *frac**c_param(c))/dtime end if rsub_top(c) = imped * rsub_tmp ! make sure baseflow isn't negative rsub_top(c) = max(0._r8, rsub_top(c)) else rsub_top(c) = imped * rsub_top_max* exp(-fff(c)*zwt(c)) end if ! use analytical expression for aquifer specific yield rous = watsat(c,nlevsoi) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevsoi))**(-1./bsw(c,nlevsoi))) rous=max(rous, params_inst%aq_sp_yield_min) !-- water table is below the soil column -------------------------------------- if(jwt(c) == nlevsoi) then wa(c) = wa(c) - rsub_top(c) * dtime zwt(c) = zwt(c) + (rsub_top(c) * dtime)/1000._r8/rous h2osoi_liq(c,nlevsoi) = h2osoi_liq(c,nlevsoi) + max(0._r8,(wa(c)-aquifer_water_baseline)) wa(c) = min(wa(c), aquifer_water_baseline) else !-- water table within soil layers 1-9 ------------------------------------- !============================== RSUB_TOP ========================================= !-- Now remove water via rsub_top rsub_top_tot = - rsub_top(c) * dtime !should never be positive... but include for completeness if(rsub_top_tot > 0.) then !rising water table call endrun(msg="RSUB_TOP IS POSITIVE in Drainage!"//errmsg(sourcefile, __LINE__)) else ! deepening water table if (use_vichydro) then wtsub_vic = 0._r8 do j = (nlvic(1)+nlvic(2)+1), nlevsoi wtsub_vic = wtsub_vic + hk_l(c,j)*dzmm(c,j) end do do j = (nlvic(1)+nlvic(2)+1), nlevsoi rsub_top_layer=max(rsub_top_tot, rsub_top_tot*hk_l(c,j)*dzmm(c,j)/wtsub_vic) rsub_top_layer=min(rsub_top_layer,0._r8) h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer rsub_top_tot = rsub_top_tot - rsub_top_layer end do else do j = jwt(c)+1, nlevsoi ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y, params_inst%aq_sp_yield_min) rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) rsub_top_layer=min(rsub_top_layer,0._r8) h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer rsub_top_tot = rsub_top_tot - rsub_top_layer if (rsub_top_tot >= 0.) then zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8 exit else zwt(c) = zi(c,j) endif enddo end if !-- remove residual rsub_top --------------------------------------------- zwt(c) = zwt(c) - rsub_top_tot/1000._r8/rous wa(c) = wa(c) + rsub_top_tot endif !-- recompute jwt --------------------------------------------------------- ! allow jwt to equal zero when zwt is in top layer jwt(c) = nlevsoi do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo end if! end of jwt if construct zwt(c) = max(0.0_r8,zwt(c)) zwt(c) = min(80._r8,zwt(c)) endif end do ! excessive water above saturation added to the above unsaturated layer like a bucket ! if column fully saturated, excess water goes to runoff do j = nlevsoi,2,-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8) h2osoi_liq(c,j) = min(eff_porosity(c,j)*dzmm(c,j), h2osoi_liq(c,j)) h2osoi_liq(c,j-1) = h2osoi_liq(c,j-1) + xsi(c) end do end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! watmin addition to fix water balance errors xs1(c) = max(max(h2osoi_liq(c,1)-watmin,0._r8)- & max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1)-watmin)),0._r8) h2osoi_liq(c,1) = h2osoi_liq(c,1) - xs1(c) if (lun%urbpoi(col%landunit(c))) then qflx_rsub_sat(c) = xs1(c) / dtime else if(h2osfcflag == 1) then ! send this water up to h2osfc rather than sending to drainage h2osfc(c) = h2osfc(c) + xs1(c) qflx_rsub_sat(c) = 0._r8 else ! use original code to send water to drainage (non-h2osfc case) qflx_rsub_sat(c) = xs1(c) / dtime endif endif ! add in ice check xs1(c) = max(max(h2osoi_ice(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1))),0._r8) h2osoi_ice(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1)), h2osoi_ice(c,1)) qflx_ice_runoff_xs(c) = xs1(c) / dtime end do ! Limit h2osoi_liq to be greater than or equal to watmin. ! Get water needed to bring h2osoi_liq equal watmin from lower layer. ! If insufficient water in soil layers, get from aquifer water do j = 1, nlevsoi-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin - h2osoi_liq(c,j) ! deepen water table if water is passed from below zwt layer if(j == jwt(c)) then zwt(c) = zwt(c) + xs(c)/eff_porosity(c,j)/1000._r8 endif else xs(c) = 0._r8 end if h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c) h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c) end do end do ! Get water for bottom layer from layers above if possible j = nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin-h2osoi_liq(c,j) searchforwater: do i = nlevsoi-1, 1, -1 available_h2osoi_liq = max(h2osoi_liq(c,i)-watmin-xs(c),0._r8) if (available_h2osoi_liq >= xs(c)) then h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) h2osoi_liq(c,i) = h2osoi_liq(c,i) - xs(c) xs(c) = 0._r8 exit searchforwater else h2osoi_liq(c,j) = h2osoi_liq(c,j) + available_h2osoi_liq h2osoi_liq(c,i) = h2osoi_liq(c,i) - available_h2osoi_liq xs(c) = xs(c) - available_h2osoi_liq end if end do searchforwater else xs(c) = 0._r8 end if ! Needed in case there is no water to be found h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) ! Instead of removing water from aquifer where it eventually ! shows up as excess drainage to the ocean, take it back out of ! drainage rsub_top(c) = rsub_top(c) - xs(c)/dtime end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Sub-surface runoff and drainage qflx_drain(c) = qflx_rsub_sat(c) + rsub_top(c) ! Set imbalance for snow capping qflx_qrgwl(c) = qflx_snwcp_liq(c) end do ! No drainage for urban columns (except for pervious road as computed above) do fc = 1, num_urbanc c = filter_urbanc(fc) if (col%itype(c) /= icol_road_perv) then qflx_drain(c) = 0._r8 ! This must be done for roofs and impervious road (walls will be zero) qflx_qrgwl(c) = qflx_snwcp_liq(c) end if end do end associate end subroutine Drainage !----------------------------------------------------------------------- subroutine CLMVICMap(bounds, numf, filter, & soilhydrology_inst, waterstatebulk_inst) ! ! !DESCRIPTION: ! Performs the mapping from CLM layers to VIC layers ! CLM hydrologically active soil layers are mapped to three VIC layers ! by assigning the first nlvic(1) layers to VIC layer 1 ! the next nlvic(2) layers to VIC alyer 2 ! and the remaining to VIC layer 3 ! mapping from VIC to CLM layers, M.Huang ! ! !USES: use clm_varcon , only : denh2o, denice, watmin use decompMod , only : bounds_type ! ! !REVISION HISTORY: ! Created by Maoyi Huang ! 11/13/2012, Maoyi Huang: rewrite the mapping modules in CLM4VIC ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: numf ! number of column soil points in column filter integer , intent(in) :: filter(:) ! column filter for soil points type(waterstatebulk_type) , intent(in) :: waterstatebulk_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst ! ! !LOCAL VARIABLES real(r8) :: ice0(1:nlayer) ! last step ice lens (mm) (new) real(r8) :: moist0(1:nlayer) ! last step soil water (mm) (new) integer :: i, j, c, fc ! note: in CLM4 h2osoil_liq unit is kg/m2, in VIC moist is mm ! h2osoi_ice is actually water equivalent ice content. !----------------------------------------------------------------------- associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer thickness (m) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] (nlevgrnd) depth => soilhydrology_inst%depth_col , & ! Input: [real(r8) (:,:) ] layer depth of upper layer (m) porosity => soilhydrology_inst%porosity_col , & ! Input: [real(r8) (:,:) ] soil porisity (1-bulk_density/soil_density) max_moist => soilhydrology_inst%max_moist_col , & ! Input: [real(r8) (:,:) ] max layer moist + ice (mm) vic_clm_fract => soilhydrology_inst%vic_clm_fract_col , & ! Input: [real(r8) (:,:,:) ] fraction of VIC layers in each CLM layer moist => soilhydrology_inst%moist_col , & ! Output: [real(r8) (:,:) ] liquid water (mm) ice => soilhydrology_inst%ice_col , & ! Output: [real(r8) (:,:) ] ice lens (mm) moist_vol => soilhydrology_inst%moist_vol_col , & ! Output: [real(r8) (:,:) ] volumetric soil moisture for VIC soil layers top_moist => soilhydrology_inst%top_moist_col , & ! Output: [real(r8) (:) ] soil moisture in top VIC layers top_max_moist => soilhydrology_inst%top_max_moist_col, & ! Output: [real(r8) (:) ] maximum soil moisture in top VIC layers top_ice => soilhydrology_inst%top_ice_col , & ! Output: [real(r8) (:) ] ice len in top VIC layers top_moist_limited => soilhydrology_inst%top_moist_limited_col & ! Output: [real(r8) (:) ] soil moisture in top layers, limited to no greater than top_max_moist ) ! map CLM to VIC do fc = 1, numf c = filter(fc) do i = 1, nlayer ice0(i) = ice(c,i) moist0(i) = moist(c,i) ice(c,i) = 0._r8 moist(c,i) = 0._r8 do j = 1, nlevsoi ice(c,i) = ice(c,i) + h2osoi_ice(c,j) * vic_clm_fract(c,i,j) moist(c,i) = moist(c,i) + h2osoi_liq(c,j) * vic_clm_fract(c,i,j) end do ice(c,i) = min((moist0(i) + ice0(i)), ice(c,i)) ice(c,i) = max(0._r8, ice(c,i)) moist(c,i) = max(watmin, moist(c,i)) moist(c,i) = min(max_moist(c,i)-ice(c,i), moist(c,i)) moist_vol(c,i) = moist(c,i)/(depth(c,i)*denice) + ice(c,i)/(depth(c,i)*denh2o) moist_vol(c,i) = min(porosity(c,i), moist_vol(c,i)) moist_vol(c,i) = max(0.01_r8, moist_vol(c,i)) end do ! hydrologic inactive layers ice(c, nlayer+1:nlayert) = h2osoi_ice(c, nlevsoi+1:nlevgrnd) moist(c, nlayer+1:nlayert) = h2osoi_liq(c, nlevsoi+1:nlevgrnd) moist_vol(c, nlayer+1:nlayert) = h2osoi_vol(c, nlevsoi+1:nlevgrnd) end do ! Set values related to top VIC layers do fc = 1, numf c = filter(fc) top_moist(c) = 0._r8 top_ice(c) = 0._r8 top_max_moist(c) = 0._r8 end do do j = 1, nlayer - 1 do fc = 1, numf c = filter(fc) top_ice(c) = top_ice(c) + ice(c,j) top_moist(c) = top_moist(c) + moist(c,j) + ice(c,j) top_max_moist(c) = top_max_moist(c) + max_moist(c,j) end do end do do fc = 1, numf c = filter(fc) top_moist_limited(c) = min(top_moist(c), top_max_moist(c)) end do end associate end subroutine CLMVICMap !#3 !----------------------------------------------------------------------- subroutine PerchedWaterTable(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, & temperature_inst, waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate watertable, considering aquifer recharge but no drainage. ! ! !USES: use clm_varcon , only : pondmx, tfrz, watmin,denice,denh2o use column_varcon , only : icol_roof, icol_road_imperv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(in) :: soilstate_inst type(temperature_type) , intent(in) :: temperature_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: c,j,fc,i ! indices real(r8) :: s_y integer :: k,k_frz,k_perch,k_zwt real(r8) :: sat_lev real(r8) :: s1 real(r8) :: s2 real(r8) :: m real(r8) :: b integer :: sat_flag !----------------------------------------------------------------------- associate( & dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) zwt => soilhydrology_inst%zwt_col , & ! Output: [real(r8) (:) ] water table depth (m) zwt_perched => soilhydrology_inst%zwt_perched_col , & ! Output: [real(r8) (:) ] perched water table depth (m) frost_table => soilhydrology_inst%frost_table_col , & ! Output: [real(r8) (:) ] frost table depth (m) origflag => soilhydrology_inst%origflag & ! Input: logical ) ! calculate perched water table location do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! define frost table as first frozen layer with unfrozen layer above it if(t_soisno(c,1) > tfrz) then k_frz=nlevsoi else k_frz=1 endif do k=2, nlevsoi if (t_soisno(c,k-1) > tfrz .and. t_soisno(c,k) <= tfrz) then k_frz=k exit endif enddo frost_table(c)=z(c,k_frz) ! initialize perched water table to frost table, and qflx_drain_perched(c) to zero zwt_perched(c)=frost_table(c) !======= water table above frost table =================== ! if water table is above frost table, do nothing if (zwt(c) < frost_table(c) .and. t_soisno(c,k_frz) <= tfrz & .and. origflag == 0) then else !========== water table below frost table ============ ! locate perched water table from bottom up starting at ! frost table sat_lev is an arbitrary saturation level ! used to determine perched water table sat_lev=0.9 k_perch=1 do k=k_frz,1,-1 h2osoi_vol(c,k) = h2osoi_liq(c,k)/(dz(c,k)*denh2o) & + h2osoi_ice(c,k)/(dz(c,k)*denice) if (h2osoi_vol(c,k)/watsat(c,k) <= sat_lev) then k_perch=k exit endif enddo ! if frost_table = nlevsoi, only compute perched water table if frozen if (t_soisno(c,k_frz) > tfrz) k_perch=k_frz ! if perched water table exists ! interpolate between k_perch and k_perch+1 to find ! perched water table height if (k_frz > k_perch) then s1 = (h2osoi_liq(c,k_perch)/(dz(c,k_perch)*denh2o) & + h2osoi_ice(c,k_perch)/(dz(c,k_perch)*denice))/watsat(c,k_perch) s2 = (h2osoi_liq(c,k_perch+1)/(dz(c,k_perch+1)*denh2o) & + h2osoi_ice(c,k_perch+1)/(dz(c,k_perch+1)*denice))/watsat(c,k_perch+1) m=(z(c,k_perch+1)-z(c,k_perch))/(s2-s1) b=z(c,k_perch+1)-m*s2 zwt_perched(c)=max(0._r8,m*sat_lev+b) endif !k_frz > k_perch endif end do end associate end subroutine PerchedWaterTable !#4 !----------------------------------------------------------------------- subroutine PerchedLateralFlow(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate subsurface drainage from perched saturated zone ! ! !USES: use clm_varcon , only : pondmx, tfrz, watmin,rpi, secspday, nlvic use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilstate_type) , intent(in) :: soilstate_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'PerchedLateralFlow' ! subroutine name integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi) ! layer thickness (mm) real(r8) :: wtsub ! summation of hk*dzmm for layers below water table (mm**2/s) real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-) real(r8) :: fracice_rsub(bounds%begc:bounds%endc) ! fractional impermeability of soil layers (-) real(r8) :: h2osoi_vol real(r8) :: imped real(r8) :: drainage_tot real(r8) :: drainage_layer real(r8) :: s_y integer :: k,k_frz,k_perch real(r8) :: sat_lev real(r8) :: s1, s2, m, b real(r8) :: q_perch real(r8) :: q_perch_max real(r8) :: vol_ice !----------------------------------------------------------------------- associate( & z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) zwt_perched => soilhydrology_inst%zwt_perched_col , & ! Input: [real(r8) (:) ] perched water table depth (m) origflag => soilhydrology_inst%origflag , & ! Input: logical qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col , & ! Output: [real(r8) (:) ] perched wt sub-surface runoff (mm H2O /s) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) ) ! Get time step dtime = get_step_size_real() ! Compute ice fraction in each layer do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) end do end do ! compute drainage from perched saturated region do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain_perched(c) = 0._r8 if ((frost_table(c) > zwt_perched(c)) .and. origflag == 0) then ! specify maximum drainage rate q_perch_max = 1.e-5_r8 * sin(col%topo_slope(c) * (rpi/180._r8)) ! calculate frost table and perched water table locations do k=1, nlevsoi if (frost_table(c) >= zi(c,k-1) .and. frost_table(c) <= zi(c,k)) then k_frz=k exit endif enddo do k=1, nlevsoi if (zwt_perched(c) >= zi(c,k-1) .and. zwt_perched(c) <= zi(c,k)) then k_perch=k exit endif enddo wtsub = 0._r8 q_perch = 0._r8 do k = k_perch, k_frz imped=10._r8**(-col%e_ice(c)*(0.5_r8*(icefrac(c,k)+icefrac(c,min(nlevsoi, k+1))))) q_perch = q_perch + imped*hksat(c,k)*dzmm(c,k) wtsub = wtsub + dzmm(c,k) end do if (wtsub > 0._r8) q_perch = q_perch/wtsub qflx_drain_perched(c) = q_perch_max * q_perch & *(frost_table(c) - zwt_perched(c)) ! no perched water table drainage if using original formulation if(origflag == 1) qflx_drain_perched(c) = 0._r8 ! if perched water table exists if (k_frz > k_perch) then ! remove drainage from perched saturated layers drainage_tot = - qflx_drain_perched(c) * dtime do k = k_perch+1, k_frz drainage_layer=max(drainage_tot,-(h2osoi_liq(c,k)-watmin)) drainage_layer=min(drainage_layer,0._r8) drainage_tot = drainage_tot - drainage_layer h2osoi_liq(c,k) = h2osoi_liq(c,k) + drainage_layer s_y = watsat(c,k) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,k))**(-1./bsw(c,k))) s_y=max(s_y, params_inst%aq_sp_yield_min) if (drainage_tot >= 0.) then zwt_perched(c) = zwt_perched(c) - drainage_layer/s_y/1000._r8 exit else zwt_perched(c) = zi(c,k) endif enddo ! if drainage_tot is greater than available water ! (above frost table), then decrease qflx_drain_perched ! by residual amount for water balance qflx_drain_perched(c) = qflx_drain_perched(c) + drainage_tot/dtime else qflx_drain_perched(c) = 0._r8 endif !k_frz > k_perch endif enddo end associate end subroutine PerchedLateralFlow !#41 !----------------------------------------------------------------------- subroutine PerchedLateralFlowHillslope(bounds, num_hydrologyc, & filter_hydrologyc, soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterfluxbulk_inst, wateratm2lndbulk_inst) ! ! !DESCRIPTION: ! Calculate subsurface drainage from perched saturated zone ! ! !USES: use clm_varcon , only : pondmx, tfrz, watmin,rpi, secspday, nlvic use landunit_varcon , only : istsoil ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilstate_type) , intent(in) :: soilstate_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst type(wateratm2lndbulk_type), intent(in) :: wateratm2lndbulk_inst ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'PerchedLateralFlowHillslope' ! subroutine name integer :: c,j,fc,i,g ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: drainage_tot real(r8) :: drainage_layer real(r8) :: s_y integer :: k integer :: k_frost(bounds%begc:bounds%endc) integer :: k_perch(bounds%begc:bounds%endc) real(r8) :: wtsub real(r8) :: q_perch real(r8) :: q_perch_max character(len=32) :: transmissivity_method = 'layersum' ! character(len=32) :: baseflow_method = 'kinematic' character(len=32) :: baseflow_method = 'darcy' ! character(len=32) :: baseflow_method = 'mixed' real(r8) :: transmis real(r8) :: dgrad real(r8), parameter :: k_anisotropic = 1._r8 integer :: c0, c_src, c_dst, nstep, nbase real(r8) :: qflx_drain_perched_vol(bounds%begc:bounds%endc) ! volumetric lateral subsurface flow through active layer [m3/s] real(r8) :: qflx_drain_perched_out(bounds%begc:bounds%endc) ! lateral subsurface flow through active layer [mm/s] associate( & nbedrock => col%nbedrock , & ! Input: [real(r8) (:,:) ] depth to bedrock (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) zwt_perched => soilhydrology_inst%zwt_perched_col , & ! Input: [real(r8) (:) ] perched water table depth (m) tdepth => wateratm2lndbulk_inst%tdepth_grc , & ! Input: [real(r8) (:) ] depth of water in tributary channels (m) tdepth_bankfull => wateratm2lndbulk_inst%tdepthmax_grc , & ! Input: [real(r8) (:) ] bankfull depth of tributary channels (m) qflx_drain_perched => waterfluxbulk_inst%qflx_drain_perched_col , & ! Output: [real(r8) (:) ] perched wt sub-surface runoff (mm H2O /s) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) ) ! Get time step dtime = get_step_size_real() ! locate frost table and perched water table do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) k_frost(c)=nbedrock(c) k_perch(c)=nbedrock(c) do k=1, nbedrock(c) if (frost_table(c) >= zi(c,k-1) .and. frost_table(c) <= zi(c,k)) then k_frost(c)=k exit endif enddo do k=1, nbedrock(c) if (zwt_perched(c) >= zi(c,k-1) .and. zwt_perched(c) <= zi(c,k)) then k_perch(c)=k exit endif enddo enddo ! compute drainage from perched saturated region do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain_perched(c) = 0._r8 qflx_drain_perched_out(c) = 0._r8 qflx_drain_perched_vol(c) = 0._r8 if (frost_table(c) > zwt_perched(c)) then ! hillslope columns if (lun%itype(col%landunit(c)) == istsoil) then ! calculate head gradient g = col%gridcell(c) ! kinematic wave approximation if (baseflow_method == 'kinematic') then dgrad = col%hill_slope(c) endif ! darcy's law if (baseflow_method == 'darcy') then if (col%cold(c) /= ispval) then dgrad = (col%hill_elev(c)-zwt_perched(c)) & - (col%hill_elev(col%cold(c))-zwt_perched(col%cold(c))) dgrad = dgrad / (col%hill_distance(c) - col%hill_distance(col%cold(c))) else ! flow between channel and lowest column ! bankfull height is defined to be zero dgrad = (col%hill_elev(c)-zwt_perched(c)) & ! ignore overbankfull storage - min(max((tdepth(g) - tdepth_bankfull(g)), & (col%hill_elev(c)-frost_table(c))),0._r8) dgrad = dgrad / (col%hill_distance(c)) ! dgrad cannot be negative when channel is empty if (tdepth(g) <= 0._r8) then dgrad = max(dgrad, 0._r8) endif endif end if ! Determine source and destination columns if (dgrad >= 0._r8) then c_src = c c_dst = col%cold(c) else c_src = col%cold(c) c_dst = c endif ! Calculate transmissivity of source column transmis = 0._r8 if(c_src /= ispval) then ! sum of layer transmissivities if (transmissivity_method == 'layersum') then do k = k_perch(c_src), k_frost(c_src) if(k == k_perch(c_src)) then transmis = transmis + 1.e-3_r8*hksat(c_src,k)*(zi(c_src,k) - zwt_perched(c_src)) else transmis = transmis + 1.e-3_r8*hksat(c_src,k)*dz(c_src,k) endif end do ! adjust by 'anisotropy factor' transmis = k_anisotropic*transmis endif ! constant conductivity based on shallowest saturated layer hk if (transmissivity_method == 'constant') then transmis = k_anisotropic*(1.e-3_r8*hksat(c_src,k_perch(c_src))) & *(zi(c_src,k_frost(c_src)) - zwt_perched(c_src) ) endif ! power law profile based on shallowest saturated layer hk if (transmissivity_method == 'power') then endif else ! transmissivity of losing stream (c_src == ispval) transmis = k_anisotropic*(1.e-3_r8*hksat(c,k_perch(c_dst)))*tdepth(g) endif qflx_drain_perched_vol(c) = transmis*col%hill_width(c)*dgrad qflx_drain_perched_out(c) = 1.e3_r8*(qflx_drain_perched_vol(c)/col%hill_area(c)) else ! non-hillslope columns ! specify maximum drainage rate q_perch_max = 1.e-5_r8 * sin(col%topo_slope(c) * (rpi/180._r8)) wtsub = 0._r8 q_perch = 0._r8 do k = k_perch(c), k_frost(c) q_perch = q_perch + hksat(c,k)*dz(c,k) wtsub = wtsub + dz(c,k) end do if (wtsub > 0._r8) q_perch = q_perch/wtsub qflx_drain_perched_out(c) = q_perch_max * q_perch & *(frost_table(c) - zwt_perched(c)) endif endif enddo ! compute net drainage from perched saturated region do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! drainage-out qflx_drain_perched(c) = qflx_drain_perched(c) + qflx_drain_perched_out(c) if (lun%itype(col%landunit(c)) == istsoil) then ! drainage-in if (col%cold(c) /= ispval) then qflx_drain_perched(col%cold(c)) = & qflx_drain_perched(col%cold(c)) - & 1.e3_r8*(qflx_drain_perched_vol(c))/col%hill_area(col%cold(c)) endif endif enddo ! remove drainage from soil moisture storage do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! remove drainage from perched saturated layers drainage_tot = qflx_drain_perched(c) * dtime do k = k_perch(c), k_frost(c) s_y = watsat(c,k) & * ( 1. - (1.+1.e3*zwt_perched(c)/sucsat(c,k))**(-1./bsw(c,k))) s_y=max(s_y,0.02_r8) if (k== k_perch(c)) then drainage_layer=min(drainage_tot,(s_y*(zi(c,k) - zwt_perched(c))*1.e3)) else drainage_layer=min(drainage_tot,(s_y*(dz(c,k))*1.e3)) endif drainage_layer=max(drainage_layer,0._r8) drainage_tot = drainage_tot - drainage_layer h2osoi_liq(c,k) = h2osoi_liq(c,k) - drainage_layer enddo ! if drainage_tot is greater than available water ! (above frost table), then decrease qflx_drain_perched ! by residual amount for water balance qflx_drain_perched(c) = qflx_drain_perched(c) - drainage_tot/dtime enddo end associate end subroutine PerchedLateralFlowHillslope !#5 !----------------------------------------------------------------------- subroutine ThetaBasedWaterTable(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate watertable, considering aquifer recharge but no drainage. ! ! !USES: use clm_varcon , only : denice,denh2o use column_varcon , only : icol_roof, icol_road_imperv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(in) :: soilstate_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: c,j,fc,i ! indices integer :: k,k_zwt real(r8) :: sat_lev real(r8) :: s1,s2,m,b ! temporary variables used to interpolate theta integer :: sat_flag !----------------------------------------------------------------------- associate( & nbedrock => col%nbedrock , & ! Input: [real(r8) (:,:) ] depth to bedrock (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) zwt => soilhydrology_inst%zwt_col & ! Output: [real(r8) (:) ] water table depth (m) ) ! calculate water table based on soil moisture state ! this is a simple search for 1st layer with soil moisture ! less than specified threshold (sat_lev) do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! initialize to depth of bottom of lowest layer zwt(c)=zi(c,nlevsoi) ! locate water table from bottom up starting at bottom of soil column ! sat_lev is an arbitrary saturation level used to determine water table sat_lev=0.9 k_zwt=nbedrock(c) sat_flag=1 !will remain unchanged if all layers at saturation do k=nbedrock(c),1,-1 h2osoi_vol(c,k) = h2osoi_liq(c,k)/(dz(c,k)*denh2o) & + h2osoi_ice(c,k)/(dz(c,k)*denice) if (h2osoi_vol(c,k)/watsat(c,k) <= sat_lev) then k_zwt=k sat_flag=0 exit endif enddo if (sat_flag == 1) k_zwt=1 ! if soil column above sat_lev, set water table to lower ! interface of first layer if (k_zwt == 1) then zwt(c)=zi(c,1) else if (k_zwt < nbedrock(c)) then ! interpolate between k_zwt and k_zwt+1 to find water table height s1 = (h2osoi_liq(c,k_zwt)/(dz(c,k_zwt)*denh2o) & + h2osoi_ice(c,k_zwt)/(dz(c,k_zwt)*denice))/watsat(c,k_zwt) s2 = (h2osoi_liq(c,k_zwt+1)/(dz(c,k_zwt+1)*denh2o) & + h2osoi_ice(c,k_zwt+1)/(dz(c,k_zwt+1)*denice))/watsat(c,k_zwt+1) m=(z(c,k_zwt+1)-z(c,k_zwt))/(s2-s1) b=z(c,k_zwt+1)-m*s2 zwt(c)=max(0._r8,m*sat_lev+b) else zwt(c)=zi(c,nbedrock(c)) endif end do end associate end subroutine ThetaBasedWaterTable !#6 !----------------------------------------------------------------------- subroutine LateralFlowPowerLaw(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc,soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate subsurface drainage ! ! !USES: use clm_varcon , only : pondmx, watmin,rpi, secspday, nlvic use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv use GridcellType , only : grc ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilstate_type) , intent(in) :: soilstate_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'LateralFlowPowerLaw' ! subroutine name integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: xs(bounds%begc:bounds%endc) ! water needed to bring soil moisture to watmin (mm) real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi) ! layer thickness (mm) integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-) real(r8) :: rsub_top(bounds%begc:bounds%endc) ! subsurface runoff - topographic control (mm/s) real(r8) :: xsi(bounds%begc:bounds%endc) ! excess soil water above saturation at layer i (mm) real(r8) :: xsia(bounds%begc:bounds%endc) ! available pore space at layer i (mm) real(r8) :: xs1(bounds%begc:bounds%endc) ! excess soil water above saturation at layer 1 (mm) real(r8) :: smpfz(1:nlevsoi) ! matric potential of layer right above water table (mm) real(r8) :: wtsub ! summation of hk*dzmm for layers below water table (mm**2/s) real(r8) :: dzsum ! summation of dzmm of layers below water table (mm) real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-) real(r8) :: fracice_rsub(bounds%begc:bounds%endc) ! fractional impermeability of soil layers (-) real(r8) :: available_h2osoi_liq ! available soil liquid water in a layer real(r8) :: h2osoi_vol real(r8) :: imped real(r8) :: rsub_top_tot real(r8) :: rsub_top_layer real(r8) :: theta_unsat real(r8) :: f_unsat real(r8) :: s_y integer :: k real(r8) :: s1 real(r8) :: s2 real(r8) :: m real(r8) :: b real(r8) :: vol_ice real(r8) :: dsmax_tmp(bounds%begc:bounds%endc) ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rsub_tmp ! temporary variable for ARNO subsurface runoff calculation real(r8) :: frac ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rel_moist ! relative moisture, temporary variable real(r8) :: wtsub_vic ! summation of hk*dzmm for layers in the third VIC layer integer :: g !----------------------------------------------------------------------- associate( & nbedrock => col%nbedrock , & ! Input: [real(r8) (:,:) ] depth to bedrock (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) snl => col%snl , & ! Input: [integer (:) ] number of snow layers h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s) depth => soilhydrology_inst%depth_col , & ! Input: [real(r8) (:,:) ] VIC soil depth c_param => soilhydrology_inst%c_param_col , & ! Input: [real(r8) (:) ] baseflow exponent (Qb) Dsmax => soilhydrology_inst%dsmax_col , & ! Input: [real(r8) (:) ] max. velocity of baseflow (mm/day) max_moist => soilhydrology_inst%max_moist_col , & ! Input: [real(r8) (:,:) ] maximum soil moisture (ice + liq) moist => soilhydrology_inst%moist_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) Ds => soilhydrology_inst%ds_col , & ! Input: [real(r8) (:) ] fracton of Dsmax where non-linear baseflow begins Wsvic => soilhydrology_inst%Wsvic_col , & ! Input: [real(r8) (:) ] fraction of maximum soil moisutre where non-liear base flow occurs icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) wa => waterstatebulk_inst%wa_col , & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm) ice => soilhydrology_inst%ice_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s) origflag => soilhydrology_inst%origflag , & ! Input: logical h2osfcflag => soilhydrology_inst%h2osfcflag , & ! Input: integer qflx_snwcp_liq => waterfluxbulk_inst%qflx_snwcp_liq_col , & ! Output: [real(r8) (:) ] excess rainfall due to snow capping (mm H2O /s) [+] qflx_ice_runoff_xs => waterfluxbulk_inst%qflx_ice_runoff_xs_col , & ! Output: [real(r8) (:) ] solid runoff from excess ice in soil (mm H2O /s) [+] qflx_liqdew_to_top_layer => waterfluxbulk_inst%qflx_liqdew_to_top_layer_col , & ! Output: [real(r8) (:) ] rate of liquid water deposited on top soil or snow layer (dew) (mm H2O /s) [+] qflx_soliddew_to_top_layer => waterfluxbulk_inst%qflx_soliddew_to_top_layer_col , & ! Output: [real(r8) (:) ] rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+] qflx_solidevap_from_top_layer => waterfluxbulk_inst%qflx_solidevap_from_top_layer_col, & ! Output: [real(r8) (:) ] rate of ice evaporated from top soil or snow layer (sublimation) (mm H2O /s) [+] qflx_drain => waterfluxbulk_inst%qflx_drain_col , & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s) qflx_qrgwl => waterfluxbulk_inst%qflx_qrgwl_col , & ! Output: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes (mm H2O /s) qflx_rsub_sat => waterfluxbulk_inst%qflx_rsub_sat_col , & ! Output: [real(r8) (:) ] soil saturation excess [mm h2o/s] h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) ) ! Get time step dtime = get_step_size_real() ! Convert layer thicknesses from m to mm do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) end do end do ! Initial set do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain(c) = 0._r8 qflx_rsub_sat(c) = 0._r8 rsub_top(c) = 0._r8 fracice_rsub(c) = 0._r8 end do ! The layer index of the first unsaturated layer, ! i.e., the layer right above the water table do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) jwt(c) = nlevsoi ! allow jwt to equal zero when zwt is in top layer do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo end do !-- Topographic runoff ------------------------- do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzsum = 0._r8 icefracsum = 0._r8 do j = max(jwt(c),1), nlevsoi dzsum = dzsum + dzmm(c,j) icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j) end do imped=10._r8**(-col%e_ice(c)*(icefracsum/dzsum)) !@@ ! baseflow is power law expression relative to bedrock layer if(zwt(c) <= zi(c,nbedrock(c))) then rsub_top(c) = imped * baseflow_scalar * tan(rpi/180._r8*col%topo_slope(c))* & (zi(c,nbedrock(c)) - zwt(c))**(params_inst%n_baseflow) else rsub_top(c) = 0._r8 endif !-- Now remove water via rsub_top rsub_top_tot = - rsub_top(c)* dtime !should never be positive... but include for completeness if(rsub_top_tot > 0.) then !rising water table call endrun(msg="RSUB_TOP IS POSITIVE in Drainage!"//errmsg(sourcefile, __LINE__)) else ! deepening water table do j = jwt(c)+1, nbedrock(c) ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y, params_inst%aq_sp_yield_min) rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) rsub_top_layer=min(rsub_top_layer,0._r8) h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer rsub_top_tot = rsub_top_tot - rsub_top_layer if (rsub_top_tot >= 0.) then zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8 exit else zwt(c) = zi(c,j) endif enddo !-- remove residual rsub_top -------------------------------- ! make sure no extra water removed from soil column rsub_top(c) = rsub_top(c) - rsub_top_tot/dtime endif zwt(c) = max(0.0_r8,zwt(c)) zwt(c) = min(80._r8,zwt(c)) end do ! excessive water above saturation added to the above unsaturated layer like a bucket ! if column fully saturated, excess water goes to runoff do j = nlevsoi,2,-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8) h2osoi_liq(c,j) = min(eff_porosity(c,j)*dzmm(c,j), h2osoi_liq(c,j)) h2osoi_liq(c,j-1) = h2osoi_liq(c,j-1) + xsi(c) end do end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! watmin addition to fix water balance errors xs1(c) = max(max(h2osoi_liq(c,1)-watmin,0._r8)- & max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1)-watmin)),0._r8) h2osoi_liq(c,1) = h2osoi_liq(c,1) - xs1(c) if (lun%urbpoi(col%landunit(c))) then qflx_rsub_sat(c) = xs1(c) / dtime else ! send this water up to h2osfc rather than sending to drainage h2osfc(c) = h2osfc(c) + xs1(c) qflx_rsub_sat(c) = 0._r8 endif ! add in ice check xs1(c) = max(max(h2osoi_ice(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1))),0._r8) h2osoi_ice(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1)), h2osoi_ice(c,1)) qflx_ice_runoff_xs(c) = xs1(c) / dtime end do ! Limit h2osoi_liq to be greater than or equal to watmin. ! Get water needed to bring h2osoi_liq equal watmin from lower layer. ! If insufficient water in soil layers, get from aquifer water do j = 1, nlevsoi-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin - h2osoi_liq(c,j) ! deepen water table if water is passed from below zwt layer if(j == jwt(c)) then zwt(c) = zwt(c) + xs(c)/eff_porosity(c,j)/1000._r8 endif else xs(c) = 0._r8 end if h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c) h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c) end do end do ! Get water for bottom layer from layers above if possible j = nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin-h2osoi_liq(c,j) searchforwater: do i = nlevsoi-1, 1, -1 available_h2osoi_liq = max(h2osoi_liq(c,i)-watmin-xs(c),0._r8) if (available_h2osoi_liq >= xs(c)) then h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) h2osoi_liq(c,i) = h2osoi_liq(c,i) - xs(c) xs(c) = 0._r8 exit searchforwater else h2osoi_liq(c,j) = h2osoi_liq(c,j) + available_h2osoi_liq h2osoi_liq(c,i) = h2osoi_liq(c,i) - available_h2osoi_liq xs(c) = xs(c) - available_h2osoi_liq end if end do searchforwater else xs(c) = 0._r8 end if ! Needed in case there is no water to be found h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) ! Instead of removing water from aquifer where it eventually ! shows up as excess drainage to the ocean, take it back out of ! drainage qflx_rsub_sat(c) = qflx_rsub_sat(c) - xs(c)/dtime end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Sub-surface runoff and drainage qflx_drain(c) = qflx_rsub_sat(c) + rsub_top(c) ! Set imbalance for snow capping qflx_qrgwl(c) = qflx_snwcp_liq(c) end do ! No drainage for urban columns (except for pervious road as computed above) do fc = 1, num_urbanc c = filter_urbanc(fc) if (col%itype(c) /= icol_road_perv) then qflx_drain(c) = 0._r8 ! This must be done for roofs and impervious road (walls will be zero) qflx_qrgwl(c) = qflx_snwcp_liq(c) end if end do end associate end subroutine LateralFlowPowerLaw !#61 !----------------------------------------------------------------------- subroutine LateralFlowHillslope(bounds, & num_hillslope, filter_hillslopec, & num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc,soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterfluxbulk_inst, wateratm2lndbulk_inst) ! ! !DESCRIPTION: ! Calculate subsurface drainage ! ! !USES: use clm_time_manager , only : get_step_size use clm_varpar , only : nlevsoi, nlevgrnd, nlayer, nlayert use clm_varcon , only : pondmx, watmin,rpi, secspday, nlvic use column_varcon , only : icol_roof, icol_road_imperv, icol_road_perv use abortutils , only : endrun use GridcellType , only : grc use landunit_varcon , only : istsoil, istcrop use clm_time_manager , only : get_nstep ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilstate_type) , intent(in) :: soilstate_inst type(wateratm2lndbulk_type) , intent(in) :: wateratm2lndbulk_inst type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(waterstatebulk_type), intent(inout) :: waterstatebulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst integer , intent(in) :: num_hillslope ! number of hillslope soil cols integer , intent(in) :: filter_hillslopec(:) ! column filter for designating all hillslope cols. ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'LateralFlowHillslope' ! subroutine name integer :: c,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: xs(bounds%begc:bounds%endc) ! water needed to bring soil moisture to watmin (mm) real(r8) :: dzmm(bounds%begc:bounds%endc,1:nlevsoi) ! layer thickness (mm) integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-) real(r8) :: rsub_top(bounds%begc:bounds%endc) ! subsurface runoff - topographic control (mm/s) real(r8) :: fff(bounds%begc:bounds%endc) ! decay factor (m-1) real(r8) :: xsi(bounds%begc:bounds%endc) ! excess soil water above saturation at layer i (mm) real(r8) :: xsia(bounds%begc:bounds%endc) ! available pore space at layer i (mm) real(r8) :: xs1(bounds%begc:bounds%endc) ! excess soil water above saturation at layer 1 (mm) real(r8) :: smpfz(1:nlevsoi) ! matric potential of layer right above water table (mm) real(r8) :: wtsub ! summation of hk*dzmm for layers below water table (mm**2/s) real(r8) :: dzsum ! summation of dzmm of layers below water table (mm) real(r8) :: icefracsum ! summation of icefrac*dzmm of layers below water table (-) real(r8) :: ice_imped_col(bounds%begc:bounds%endc) ! column average hydraulic conductivity reduction due to presence of soil ice (-) real(r8) :: ice_imped(bounds%begc:bounds%endc,1:nlevsoi) ! hydraulic conductivity reduction due to presence of soil ice (-) real(r8) :: available_h2osoi_liq ! available soil liquid water in a layer real(r8) :: h2osoi_vol real(r8) :: rsub_top_tot real(r8) :: rsub_top_layer real(r8) :: theta_unsat real(r8) :: f_unsat real(r8) :: s_y integer :: k real(r8) :: s1 real(r8) :: s2 real(r8) :: m real(r8) :: b real(r8) :: vol_ice real(r8) :: dsmax_tmp(bounds%begc:bounds%endc) ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rsub_tmp ! temporary variable for ARNO subsurface runoff calculation real(r8) :: frac ! temporary variable for ARNO subsurface runoff calculation real(r8) :: rel_moist ! relative moisture, temporary variable real(r8) :: wtsub_vic ! summation of hk*dzmm for layers in the third VIC layer integer :: g character(len=32) :: transmissivity_method = 'layersum' ! character(len=32) :: baseflow_method = 'kinematic' character(len=32) :: baseflow_method = 'darcy' ! character(len=32) :: baseflow_method = 'mixed' logical :: no_lateral_flow = .false. real(r8) :: transmis ! transmissivity real(r8) :: dgrad ! hydraulic head gradient real(r8), parameter :: n_baseflow = 1 ! drainage power law exponent real(r8), parameter :: k_anisotropic = 1._r8 ! anisotropy scalar real(r8) :: qflx_latflow_out_vol(bounds%begc:bounds%endc) real(r8) :: qflx_net_latflow(bounds%begc:bounds%endc) real(r8) :: qflx_latflow_avg(bounds%begc:bounds%endc) real(r8) :: larea integer :: c0, c_src, c_dst, nstep, nbase integer :: l !----------------------------------------------------------------------- associate( & nbedrock => col%nbedrock , & ! Input: [real(r8) (:,:) ] depth to bedrock (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) snl => col%snl , & ! Input: [integer (:) ] number of snow layers h2osfc => waterstatebulk_inst%h2osfc_col , & ! Input: [real(r8) (:) ] surface water (mm) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) eff_porosity => soilstate_inst%eff_porosity_col , & ! Input: [real(r8) (:,:) ] effective porosity = porosity - vol_ice hk_l => soilstate_inst%hk_l_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity (mm/s) qflx_latflow_out => waterfluxbulk_inst%qflx_latflow_out_col, & ! Output: [real(r8) (:) ] lateral saturated outflow (mm/s) qflx_latflow_in => waterfluxbulk_inst%qflx_latflow_in_col, & ! Output: [real(r8) (:) ] lateral saturated inflow (mm/s) qdischarge => waterfluxbulk_inst%qdischarge_col , & ! Output: [real(r8) (:) ] discharge from column (m3/s) tdepth => wateratm2lndbulk_inst%tdepth_grc , & ! Input: [real(r8) (:) ] depth of water in tributary channels (m) tdepth_bankfull => wateratm2lndbulk_inst%tdepthmax_grc , & ! Input: [real(r8) (:) ] bankfull depth of tributary channels (m) depth => soilhydrology_inst%depth_col , & ! Input: [real(r8) (:,:) ] VIC soil depth c_param => soilhydrology_inst%c_param_col , & ! Input: [real(r8) (:) ] baseflow exponent (Qb) Dsmax => soilhydrology_inst%dsmax_col , & ! Input: [real(r8) (:) ] max. velocity of baseflow (mm/day) max_moist => soilhydrology_inst%max_moist_col , & ! Input: [real(r8) (:,:) ] maximum soil moisture (ice + liq) moist => soilhydrology_inst%moist_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) Ds => soilhydrology_inst%ds_col , & ! Input: [real(r8) (:) ] fracton of Dsmax where non-linear baseflow begins Wsvic => soilhydrology_inst%Wsvic_col , & ! Input: [real(r8) (:) ] fraction of maximum soil moisutre where non-liear base flow occurs icefrac => soilhydrology_inst%icefrac_col , & ! Output: [real(r8) (:,:) ] fraction of ice in layer hkdepth => soilhydrology_inst%hkdepth_col , & ! Input: [real(r8) (:) ] decay factor (m) frost_table => soilhydrology_inst%frost_table_col , & ! Input: [real(r8) (:) ] frost table depth (m) zwt => soilhydrology_inst%zwt_col , & ! Input: [real(r8) (:) ] water table depth (m) wa => waterstatebulk_inst%wa_col , & ! Input: [real(r8) (:) ] water in the unconfined aquifer (mm) ice => soilhydrology_inst%ice_col , & ! Input: [real(r8) (:,:) ] soil layer moisture (mm) qcharge => soilhydrology_inst%qcharge_col , & ! Input: [real(r8) (:) ] aquifer recharge rate (mm/s) origflag => soilhydrology_inst%origflag , & ! Input: logical h2osfcflag => soilhydrology_inst%h2osfcflag , & ! Input: logical qflx_snwcp_liq => waterfluxbulk_inst%qflx_snwcp_liq_col , & ! Output: [real(r8) (:) ] excess rainfall due to snow capping (mm H2O /s) [+] qflx_ice_runoff_xs => waterfluxbulk_inst%qflx_ice_runoff_xs_col , & ! Output: [real(r8) (:) ] solid runoff from excess ice in soil (mm H2O /s) [+] qflx_drain => waterfluxbulk_inst%qflx_drain_col , & ! Output: [real(r8) (:) ] sub-surface runoff (mm H2O /s) qflx_qrgwl => waterfluxbulk_inst%qflx_qrgwl_col , & ! Output: [real(r8) (:) ] qflx_surf at glaciers, wetlands, lakes (mm H2O /s) qflx_rsub_sat => waterfluxbulk_inst%qflx_rsub_sat_col , & ! Output: [real(r8) (:) ] soil saturation excess [mm h2o/s] h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) ) ! Get time step dtime = get_step_size_real() nstep = get_nstep() ! Convert layer thicknesses from m to mm do j = 1,nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzmm(c,j) = dz(c,j)*1.e3_r8 vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) icefrac(c,j) = min(1._r8,vol_ice/watsat(c,j)) ice_imped(c,j)=10._r8**(-col%e_ice(c)*icefrac(c,j)) end do end do ! Initial set do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) qflx_drain(c) = 0._r8 qflx_rsub_sat(c) = 0._r8 rsub_top(c) = 0._r8 qflx_latflow_in(c) = 0._r8 qflx_latflow_out(c) = 0._r8 qflx_net_latflow(c) = 0._r8 qdischarge(c) = 0._r8 qflx_latflow_out_vol(c) = 0._r8 end do ! The layer index of the first unsaturated layer, ! i.e., the layer right above the water table do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) jwt(c) = nlevsoi ! allow jwt to equal zero when zwt is in top layer do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo end do ! Calculate ice impedance factor (after jwt calculated) do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) dzsum = 0._r8 icefracsum = 0._r8 do j = max(jwt(c),1), nlevsoi dzsum = dzsum + dzmm(c,j) icefracsum = icefracsum + icefrac(c,j) * dzmm(c,j) end do ice_imped_col(c)=10._r8**(-col%e_ice(c)*(icefracsum/dzsum)) enddo do fc = 1, num_hillslope c = filter_hillslopec(fc) g = col%gridcell(c) ! kinematic wave approximation if (baseflow_method == 'kinematic') then dgrad = col%hill_slope(c) endif ! darcy's law if (baseflow_method == 'darcy') then if (col%cold(c) /= ispval) then dgrad = (col%hill_elev(c)-zwt(c)) & - (col%hill_elev(col%cold(c))-zwt(col%cold(c))) dgrad = dgrad / (col%hill_distance(c) - col%hill_distance(col%cold(c))) else ! flow between channel and lowest column ! bankfull height is defined to be zero dgrad = (col%hill_elev(c)-zwt(c)) & ! - (tdepth(g) - tdepth_bankfull(g)) ! ignore overbankfull storage - min((tdepth(g) - tdepth_bankfull(g)),0._r8) dgrad = dgrad / (col%hill_distance(c)) ! dgrad cannot be negative when channel is empty if (tdepth(g) <= 0._r8) then dgrad = max(dgrad, 0._r8) endif ! add vertical drainage for losing streams ! (this could be a separate term from lateral flow...) if (dgrad < 0._r8) then ! dgrad = dgrad - 1._r8 ! adjust lateral gradient w/ k_anisotropic dgrad = dgrad - 1._r8/k_anisotropic endif endif end if ! Calculate transmissivity of source column if (dgrad >= 0._r8) then c_src = c else c_src = col%cold(c) endif transmis = 0._r8 if(c_src /= ispval) then ! transmissivity non-zero only when saturated conditions exist if(zwt(c_src) <= zi(c_src,nbedrock(c_src))) then ! sum of layer transmissivities if (transmissivity_method == 'layersum') then do j = jwt(c_src)+1, nbedrock(c_src) if(j == jwt(c_src)+1) then transmis = transmis + 1.e-3_r8*ice_imped(c_src,j)*hksat(c_src,j)*(zi(c_src,j) - zwt(c_src)) else transmis = transmis + 1.e-3_r8*ice_imped(c_src,j)*hksat(c_src,j)*dz(c_src,j) endif end do endif ! constant conductivity based on shallowest saturated layer hk if (transmissivity_method == 'constant') then transmis = (1.e-3_r8*ice_imped(c_src,jwt(c_src)+1)*hksat(c_src,jwt(c_src)+1)) & *(zi(c_src,nbedrock(c_src)) - zwt(c_src) ) endif ! power law profile based on shallowest saturated layer hk if (transmissivity_method == 'power') then ! transmis = ice_imped(c_src,jwt(c_src)+1)*hksat(c_src,jwt(c_src)+1)*0.001_r8*dzsumall* & ! ((1-1000._r8*zwt(c_src)/dzsumall)**n_baseflow )/n_baseflow ! (m2/s) endif endif else ! transmissivity of losing stream (c_src == ispval) transmis = (1.e-3_r8*ice_imped(c,jwt(c)+1)*hksat(c,jwt(c)+1))*tdepth(g) endif ! adjust transmissivity by 'anisotropy factor' transmis = k_anisotropic*transmis ! the qflx_latflow_out_vol calculations use the ! transmissivity to determine whether saturated flow ! conditions exist, b/c gradients will be nonzero ! even when no saturated layers are present ! qflx_latflow_out_vol(c) = ice_imped(c)*transmis*col%hill_width(c)*dgrad ! include ice impedance in transmissivity qflx_latflow_out_vol(c) = transmis*col%hill_width(c)*dgrad ! qdischarge from lowest column is qflx_latflow_out_vol ! scaled by total area of column in gridcell divided by column area if (col%cold(c) == ispval) then qdischarge(c) = qflx_latflow_out_vol(c) & *(grc%area(g)*1.e6*col%wtgcell(c)/col%hill_area(c)) endif ! convert volumetric flow to equivalent flux qflx_latflow_out(c) = 1.e3_r8*qflx_latflow_out_vol(c)/col%hill_area(c) ! hilltop column has no inflow if (col%colu(c) == ispval) then qflx_latflow_in(c) = 0._r8 endif ! current outflow is inflow to downhill column normalized by downhill area if (col%cold(c) /= ispval) then qflx_latflow_in(col%cold(c)) = qflx_latflow_in(col%cold(c)) + & 1.e3_r8*qflx_latflow_out_vol(c)/col%hill_area(col%cold(c)) endif enddo ! recalculate average flux for no-lateral flow case if(no_lateral_flow) then if (baseflow_method /= 'kinematic') then call endrun(msg="baseflow_method must be kinematic for no_lateral_flow = .true.! "//errmsg(sourcefile, __LINE__)) endif do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) l = col%landunit(c) !need to sum all columns w/ same hillslope id for each column qflx_latflow_avg(c) = 0._r8 larea = 0._r8 do c0 = lun%coli(l), lun%colf(l) if(col%hillslope_ndx(c0) == col%hillslope_ndx(c)) then qflx_latflow_avg(c) = qflx_latflow_avg(c) + qflx_latflow_out_vol(c0) larea = larea + col%hill_area(c0) endif enddo qflx_latflow_avg(c) = 1.e3_r8*qflx_latflow_avg(c)/larea enddo endif !-- Topographic runoff ------------------------- do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! net lateral flow (positive out) qflx_net_latflow(c) = qflx_latflow_out(c) - qflx_latflow_in(c) if(no_lateral_flow) then qflx_net_latflow(c) = qflx_latflow_avg(c) endif fff(c) = 1._r8/ hkdepth(c) !@@ ! baseflow if(zwt(c) <= zi(c,nbedrock(c))) then if (lun%itype(col%landunit(c)) == istsoil) then ! apply net lateral flow here rsub_top(c) = qflx_net_latflow(c) else rsub_top(c) = ice_imped_col(c) * baseflow_scalar & * tan(rpi/180._r8*col%topo_slope(c))* & (zi(c,nbedrock(c)) - zwt(c))**(n_baseflow) endif else rsub_top(c) = 0._r8 endif !-- Now remove water via rsub_top rsub_top_tot = - rsub_top(c) * dtime if(rsub_top_tot > 0.) then !rising water table do j = jwt(c)+1,1,-1 ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y,0.02_r8) rsub_top_layer=min(rsub_top_tot,(s_y*dz(c,j)*1.e3)) rsub_top_layer=max(rsub_top_layer,0._r8) h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer rsub_top_tot = rsub_top_tot - rsub_top_layer if (rsub_top_tot <= 0.) then zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8 exit else zwt(c) = zi(c,j-1) endif enddo !-- remove residual rsub_top -------------------------------- h2osfc(c) = h2osfc(c) + rsub_top_tot else ! deepening water table do j = jwt(c)+1, nbedrock(c) ! analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y,0.02_r8) rsub_top_layer=max(rsub_top_tot,-(s_y*(zi(c,j) - zwt(c))*1.e3)) rsub_top_layer=min(rsub_top_layer,0._r8) h2osoi_liq(c,j) = h2osoi_liq(c,j) + rsub_top_layer rsub_top_tot = rsub_top_tot - rsub_top_layer if (rsub_top_tot >= 0.) then zwt(c) = zwt(c) - rsub_top_layer/s_y/1000._r8 exit else zwt(c) = zi(c,j) endif enddo !-- remove residual rsub_top --------------------------------------------- ! make sure no extra water removed from soil column rsub_top(c) = rsub_top(c) + rsub_top_tot/dtime endif zwt(c) = max(0.0_r8,zwt(c)) zwt(c) = min(80._r8,zwt(c)) end do ! excessive water above saturation added to the above unsaturated layer like a bucket ! if column fully saturated, excess water goes to runoff do j = nlevsoi,2,-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8) h2osoi_liq(c,j) = min(eff_porosity(c,j)*dzmm(c,j), h2osoi_liq(c,j)) h2osoi_liq(c,j-1) = h2osoi_liq(c,j-1) + xsi(c) end do end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! watmin addition to fix water balance errors xs1(c) = max(max(h2osoi_liq(c,1)-watmin,0._r8)- & max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_ice(c,1)-watmin)),0._r8) h2osoi_liq(c,1) = h2osoi_liq(c,1) - xs1(c) if (lun%urbpoi(col%landunit(c))) then qflx_rsub_sat(c) = xs1(c) / dtime else ! send this water up to h2osfc rather than sending to drainage h2osfc(c) = h2osfc(c) + xs1(c) qflx_rsub_sat(c) = 0._r8 endif ! add in ice check xs1(c) = max(max(h2osoi_ice(c,1),0._r8)-max(0._r8,(pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1))),0._r8) h2osoi_ice(c,1) = min(max(0._r8,pondmx+watsat(c,1)*dzmm(c,1)-h2osoi_liq(c,1)), h2osoi_ice(c,1)) qflx_ice_runoff_xs(c) = xs1(c) / dtime end do ! Limit h2osoi_liq to be greater than or equal to watmin. ! Get water needed to bring h2osoi_liq equal watmin from lower layer. ! If insufficient water in soil layers, get from aquifer water do j = 1, nlevsoi-1 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin - h2osoi_liq(c,j) ! deepen water table if water is passed from below zwt layer if(j == jwt(c)) then zwt(c) = zwt(c) + xs(c)/eff_porosity(c,j)/1000._r8 endif else xs(c) = 0._r8 end if h2osoi_liq(c,j ) = h2osoi_liq(c,j ) + xs(c) h2osoi_liq(c,j+1) = h2osoi_liq(c,j+1) - xs(c) end do end do ! Get water for bottom layer from layers above if possible j = nlevsoi do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) if (h2osoi_liq(c,j) < watmin) then xs(c) = watmin-h2osoi_liq(c,j) searchforwater: do i = nlevsoi-1, 1, -1 available_h2osoi_liq = max(h2osoi_liq(c,i)-watmin-xs(c),0._r8) if (available_h2osoi_liq >= xs(c)) then h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) h2osoi_liq(c,i) = h2osoi_liq(c,i) - xs(c) xs(c) = 0._r8 exit searchforwater else h2osoi_liq(c,j) = h2osoi_liq(c,j) + available_h2osoi_liq h2osoi_liq(c,i) = h2osoi_liq(c,i) - available_h2osoi_liq xs(c) = xs(c) - available_h2osoi_liq end if end do searchforwater else xs(c) = 0._r8 end if ! Needed in case there is no water to be found h2osoi_liq(c,j) = h2osoi_liq(c,j) + xs(c) ! Instead of removing water from aquifer where it eventually ! shows up as excess drainage to the ocean, take it back out of ! drainage qflx_rsub_sat(c) = qflx_rsub_sat(c) - xs(c)/dtime end do do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Sub-surface runoff and drainage qflx_drain(c) = qflx_rsub_sat(c) + rsub_top(c) ! Set imbalance for snow capping qflx_qrgwl(c) = qflx_snwcp_liq(c) end do ! No drainage for urban columns (except for pervious road as computed above) do fc = 1, num_urbanc c = filter_urbanc(fc) if (col%itype(c) /= icol_road_perv) then qflx_drain(c) = 0._r8 ! This must be done for roofs and impervious road (walls will be zero) qflx_qrgwl(c) = qflx_snwcp_liq(c) end if end do end associate end subroutine LateralFlowHillslope !#7 !----------------------------------------------------------------------- subroutine RenewCondensation(bounds, num_hydrologyc, filter_hydrologyc, & num_urbanc, filter_urbanc, soilhydrology_inst, soilstate_inst, & waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst) ! ! !DESCRIPTION: ! Calculate watertable, considering aquifer recharge but no drainage. ! ! !USES: use column_varcon , only : icol_roof, icol_road_imperv ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_hydrologyc ! number of column soil points in column filter integer , intent(in) :: num_urbanc ! number of column urban points in column filter integer , intent(in) :: filter_urbanc(:) ! column filter for urban points integer , intent(in) :: filter_hydrologyc(:) ! column filter for soil points type(soilhydrology_type) , intent(inout) :: soilhydrology_inst type(soilstate_type) , intent(in) :: soilstate_inst type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst type(waterfluxbulk_type) , intent(inout) :: waterfluxbulk_inst ! ! !LOCAL VARIABLES: integer :: c ,j,fc,i ! indices real(r8) :: dtime ! land model time step (sec) real(r8) :: qflx_solidevap_from_top_layer_save ! temporary integer :: num_modifiedc ! number of columns in filter_modifiedc integer :: filter_modifiedc(bounds%endc-bounds%begc+1) ! column filter of points modified in this subroutine real(r8) :: h2osoi_ice_before_evap(bounds%begc:bounds%endc) ! h2osoi_ice in layer 1 before applying solidevap !----------------------------------------------------------------------- associate( & snl => col%snl , & ! Input: [integer (:) ] number of snow layers h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Output: [real(r8) (:,:) ] ice lens (kg/m2) frac_h2osfc => waterdiagnosticbulk_inst%frac_h2osfc_col , & ! Input: [real(r8) (:) ] qflx_liqdew_to_top_layer => waterfluxbulk_inst%qflx_liqdew_to_top_layer_col , & ! Input: [real(r8) (:) ] rate of liquid water deposited on top soil or snow layer (dew) (mm H2O /s) [+] qflx_soliddew_to_top_layer => waterfluxbulk_inst%qflx_soliddew_to_top_layer_col, & ! Input: [real(r8) (:) ] rate of solid water deposited on top soil or snow layer (frost) (mm H2O /s) [+] qflx_ev_snow => waterfluxbulk_inst%qflx_ev_snow_col , & ! In/Out: [real(r8) (:) ] evaporation flux from snow (mm H2O/s) [+ to atm] qflx_solidevap_from_top_layer => waterfluxbulk_inst%qflx_solidevap_from_top_layer_col & ! Output: [real(r8) (:) ] rate of ice evaporated from top soil or snow layer (sublimation) (mm H2O /s) [+] ) ! Get time step dtime = get_step_size_real() num_modifiedc = 0 do fc = 1, num_hydrologyc c = filter_hydrologyc(fc) ! Renew the ice and liquid mass due to condensation if (snl(c)+1 >= 1) then num_modifiedc = num_modifiedc + 1 filter_modifiedc(num_modifiedc) = c ! make consistent with how evap_grnd removed in infiltration h2osoi_liq(c,1) = h2osoi_liq(c,1) + (1._r8 - frac_h2osfc(c))*qflx_liqdew_to_top_layer(c) * dtime h2osoi_ice(c,1) = h2osoi_ice(c,1) + (1._r8 - frac_h2osfc(c))*qflx_soliddew_to_top_layer(c) * dtime h2osoi_ice_before_evap(c) = h2osoi_ice(c,1) h2osoi_ice(c,1) = h2osoi_ice(c,1) - (1._r8 - frac_h2osfc(c)) * qflx_solidevap_from_top_layer(c) * dtime end if end do do fc = 1, num_urbanc c = filter_urbanc(fc) ! Renew the ice and liquid mass due to condensation for urban roof and impervious road if (col%itype(c) == icol_roof .or. col%itype(c) == icol_road_imperv) then if (snl(c)+1 >= 1) then num_modifiedc = num_modifiedc + 1 filter_modifiedc(num_modifiedc) = c h2osoi_liq(c,1) = h2osoi_liq(c,1) + qflx_liqdew_to_top_layer(c) * dtime h2osoi_ice(c,1) = h2osoi_ice(c,1) + (qflx_soliddew_to_top_layer(c) * dtime) h2osoi_ice_before_evap(c) = h2osoi_ice(c,1) h2osoi_ice(c,1) = h2osoi_ice(c,1) - (qflx_solidevap_from_top_layer(c) * dtime) end if end if end do call truncate_small_values( & num_f = num_modifiedc, & filter_f = filter_modifiedc, & lb = bounds%begc, & ub = bounds%endc, & data_baseline = h2osoi_ice_before_evap(bounds%begc:bounds%endc), & data = h2osoi_ice(bounds%begc:bounds%endc, 1), & custom_rel_epsilon = tolerance) do fc = 1, num_modifiedc c = filter_modifiedc(fc) if (h2osoi_ice(c,1) < 0._r8) then write(iulog,*) "ERROR: In RenewCondensation, h2osoi_ice has gone significantly negative" write(iulog,*) "c = ", c write(iulog,*) "h2osoi_ice_before_evap = ", h2osoi_ice_before_evap(c) write(iulog,*) "h2osoi_ice(c,1) = ", h2osoi_ice(c,1) write(iulog,*) "qflx_solidevap_from_top_layer*dtime = ", qflx_solidevap_from_top_layer(c)*dtime call endrun("In RenewCondensation, h2osoi_ice has gone significantly negative") end if end do end associate end subroutine RenewCondensation !#8 !----------------------------------------------------------------------- subroutine CalcIrrigWithdrawals(bounds, & num_soilc, filter_soilc, & soilhydrology_inst, soilstate_inst, & qflx_gw_demand, & qflx_gw_uncon_irrig_lyr, & qflx_gw_con_irrig) ! ! !DESCRIPTION: ! Calculate irrigation withdrawals from groundwater, given groundwater irrigation demand ! ! This routine is called when use_groundwater_irrigation = .true. ! ! It is not compatible with use_aquifer_layer ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_soilc ! number of column soil points in column filter integer , intent(in) :: filter_soilc(:) ! column filter for soil points type(soilhydrology_type) , intent(in) :: soilhydrology_inst type(soilstate_type) , intent(in) :: soilstate_inst real(r8) , intent(in) :: qflx_gw_demand( bounds%begc: ) ! groundwater irrigation demand (mm H2O/s) real(r8) , intent(inout) :: qflx_gw_uncon_irrig_lyr( bounds%begc:, 1: ) ! unconfined aquifer groundwater irrigation withdrawal flux, separated by layer (mm H2O/s) real(r8) , intent(inout) :: qflx_gw_con_irrig( bounds%begc: ) ! total confined aquifer groundwater irrigation withdrawal flux (mm H2O/s) ! ! !LOCAL VARIABLES: character(len=32) :: subname = 'CalcIrrigWithdrawals' ! subroutine name integer :: c,j,fc ! indices real(r8) :: dtime ! land model time step (sec) integer :: jwt(bounds%begc:bounds%endc) ! index of the soil layer right above the water table (-) real(r8) :: s_y real(r8) :: irrig_demand_remaining ! mm H2O real(r8) :: available_water_layer ! mm H2O real(r8) :: irrig_layer ! mm H2O !----------------------------------------------------------------------- SHR_ASSERT_ALL_FL((ubound(qflx_gw_demand) == [bounds%endc]), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(qflx_gw_uncon_irrig_lyr) == [bounds%endc, nlevsoi]), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(qflx_gw_con_irrig) == [bounds%endc]), sourcefile, __LINE__) associate( & nbedrock => col%nbedrock , & ! Input: [real(r8) (:,:) ] depth to bedrock (m) z => col%z , & ! Input: [real(r8) (:,:) ] layer depth (m) zi => col%zi , & ! Input: [real(r8) (:,:) ] interface level below a "z" level (m) dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" hksat => soilstate_inst%hksat_col , & ! Input: [real(r8) (:,:) ] hydraulic conductivity at saturation (mm H2O /s) sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) zwt => soilhydrology_inst%zwt_col & ! Input: [real(r8) (:) ] water table depth (m) ) dtime = get_step_size_real() do j = 1, nlevsoi do fc = 1, num_soilc c = filter_soilc(fc) qflx_gw_uncon_irrig_lyr(c,j) = 0._r8 end do end do !-- Remove groundwater from unconfined aquifer ----------- do fc = 1, num_soilc c = filter_soilc(fc) irrig_demand_remaining = qflx_gw_demand(c)*dtime ! should never be negative... but include for completeness if(irrig_demand_remaining < 0.) then call endrun(msg="negative groundwater irrigation demand! "//errmsg(sourcefile, __LINE__)) else jwt(c) = nlevsoi ! allow jwt to equal zero when zwt is in top layer do j = 1,nlevsoi if(zwt(c) <= zi(c,j)) then jwt(c) = j-1 exit end if enddo do j = jwt(c)+1, nbedrock(c) ! use analytical expression for specific yield s_y = watsat(c,j) & * ( 1. - (1.+1.e3*zwt(c)/sucsat(c,j))**(-1./bsw(c,j))) s_y=max(s_y, params_inst%aq_sp_yield_min) if (j==jwt(c)+1) then available_water_layer=max(0._r8,(s_y*(zi(c,j) - zwt(c))*1.e3)) else available_water_layer=max(0._r8,(s_y*(zi(c,j) - zi(c,j-1))*1.e3)) endif irrig_layer = min(irrig_demand_remaining, available_water_layer) qflx_gw_uncon_irrig_lyr(c,j) = irrig_layer / dtime irrig_demand_remaining = irrig_demand_remaining - irrig_layer if (irrig_demand_remaining <= 0.) then exit endif end do if (irrig_demand_remaining > 0._r8) then ! NOTE(wjs, 2018-12-15) Eventually we could limit irrigation from the ! confined aquifer based on water availability there. qflx_gw_con_irrig(c) = irrig_demand_remaining / dtime else qflx_gw_con_irrig(c) = 0._r8 end if end if end do end associate end subroutine CalcIrrigWithdrawals !#9 !----------------------------------------------------------------------- subroutine WithdrawGroundwaterIrrigation(bounds, & num_soilc, filter_soilc, & waterflux_inst, waterstate_inst) ! ! !DESCRIPTION: ! Remove groundwater irrigation from unconfined and confined aquifers ! This routine is called when use_groundwater_irrigation = .true. ! It is not compatible with use_aquifer_layer ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds integer , intent(in) :: num_soilc ! number of column soil points in column filter integer , intent(in) :: filter_soilc(:) ! column filter for soil points class(waterflux_type) , intent(in) :: waterflux_inst class(waterstate_type) , intent(inout) :: waterstate_inst ! ! !LOCAL VARIABLES: integer :: fc, c, j real(r8) :: dtime ! land model time step (sec) character(len=*), parameter :: subname = 'WithdrawGroundwaterIrrigation' !----------------------------------------------------------------------- associate( & qflx_gw_uncon_irrig_lyr => waterflux_inst%qflx_gw_uncon_irrig_lyr_col, & ! Input: [real(r8) (:,:) ] unconfined groundwater irrigation flux, separated by layer (mm H2O/s) qflx_gw_con_irrig => waterflux_inst%qflx_gw_con_irrig_col , & ! Input: [real(r8) (:,:) ] confined groundwater irrigation flux (mm H2O /s) wa => waterstate_inst%wa_col , & ! Output: [real(r8) (:) ] water in the unconfined aquifer (mm) h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Output: [real(r8) (:,:) ] liquid water (kg/m2) ) dtime = get_step_size_real() do j = 1, nlevsoi do fc = 1, num_soilc c = filter_soilc(fc) h2osoi_liq(c,j) = h2osoi_liq(c,j) - qflx_gw_uncon_irrig_lyr(c,j) * dtime end do end do ! zwt is not being updated, as it will be updated after HydrologyNoDrainage do fc = 1, num_soilc c = filter_soilc(fc) wa(c) = wa(c) - qflx_gw_con_irrig(c) * dtime end do end associate end subroutine WithdrawGroundwaterIrrigation !#0 end module SoilHydrologyMod