Skip to content

Commit

Permalink
Merge commit eCLM-ParFlow coupling (ESCOMP#17) '81383687999c70d529a43…
Browse files Browse the repository at this point in the history 
…9ca00749912639cdc43' into eclm
  • Loading branch information
@AGonzalezNicolas
AGonzalezNicolas committed Jun 27, 2024
2 parents 279151b + 8138368 commit 0b9c215
Show file tree
Hide file tree
Showing 16 changed files with 376 additions and 127 deletions.
62 changes: 41 additions & 21 deletions src/clm5/biogeophys/BalanceCheckMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -179,11 +179,12 @@ subroutine BalanceCheck( bounds, &
! error = abs(precipitation - change of water storage - evaporation - runoff)
!
! !USES:
use clm_varcon , only : spval
use clm_time_manager , only : get_step_size, get_nstep
use clm_time_manager , only : get_nstep_since_startup_or_lastDA_restart_or_pause
use CanopyStateType , only : canopystate_type
use SurfaceAlbedoType , only : surfalb_type
use clm_varcon , only : spval
use clm_time_manager , only : get_step_size, get_nstep
use clm_time_manager , only : get_nstep_since_startup_or_lastDA_restart_or_pause
use CanopyStateType , only : canopystate_type
use SurfaceAlbedoType , only : surfalb_type
use SoilWaterMovementMod , only : use_parflow_soilwater
use subgridAveMod
!
! !ARGUMENTS:
Expand Down Expand Up @@ -227,6 +228,8 @@ subroutine BalanceCheck( bounds, &
errh2osno => waterstate_inst%errh2osno_col , & ! Output: [real(r8) (:) ] error in h2osno (kg m-2)
endwb => waterstate_inst%endwb_col , & ! Output: [real(r8) (:) ] water mass end of the time step
total_plant_stored_h2o_col => waterstate_inst%total_plant_stored_h2o_col, & ! Input: [real(r8) (:) ] water mass in plant tissues (kg m-2)
pfl_psi => waterstate_inst%pfl_psi_col , & ! Input: [real(r8) (:,:) ] COUP_OAS_PFL

qflx_rootsoi_col => waterflux_inst%qflx_rootsoi_col , & ! Input [real(r8) (:) ] water loss in soil layers to root uptake (mm H2O/s)
! (ie transpiration demand, often = transpiration)
qflx_rain_grnd_col => waterflux_inst%qflx_rain_grnd_col , & ! Input: [real(r8) (:) ] rain on ground after interception (mm H2O/s) [+]
Expand Down Expand Up @@ -260,7 +263,6 @@ subroutine BalanceCheck( bounds, &
qflx_ice_dynbal => waterflux_inst%qflx_ice_dynbal_grc , & ! Input: [real(r8) (:) ] ice runoff due to dynamic land cover change (mm H2O /s)
snow_sources => waterflux_inst%snow_sources_col , & ! Output: [real(r8) (:) ] snow sources (mm H2O /s)
snow_sinks => waterflux_inst%snow_sinks_col , & ! Output: [real(r8) (:) ] snow sinks (mm H2O /s)

qflx_irrig => irrigation_inst%qflx_irrig_col , & ! Input: [real(r8) (:) ] irrigation flux (mm H2O /s)

qflx_glcice_dyn_water_flux => glacier_smb_inst%qflx_glcice_dyn_water_flux_col, & ! Input: [real(r8) (:)] water flux needed for balance check due to glc_dyn_runoff_routing (mm H2O/s) (positive means addition of water to the system)
Expand Down Expand Up @@ -346,7 +348,6 @@ subroutine BalanceCheck( bounds, &
- qflx_ice_runoff_xs(c) &
- qflx_snwcp_discarded_liq(c) &
- qflx_snwcp_discarded_ice(c)) * dtime

else

errh2o(c) = 0.0_r8
Expand All @@ -356,11 +357,12 @@ subroutine BalanceCheck( bounds, &
end do

found = .false.

do c = bounds%begc, bounds%endc
if (abs(errh2o(c)) > 1.e-9_r8) then
found = .true.
indexc = c
end if
if (abs(errh2o(c)) > 1.e-9_r8) then
found = .true.
indexc = c
end if
end do

if ( found ) then
Expand Down Expand Up @@ -393,19 +395,26 @@ subroutine BalanceCheck( bounds, &
write(iulog,*)'qflx_ice_runoff_xs = ',qflx_ice_runoff_xs(indexc)*dtime
write(iulog,*)'qflx_snwcp_discarded_ice = ',qflx_snwcp_discarded_ice(indexc)*dtime
write(iulog,*)'qflx_snwcp_discarded_liq = ',qflx_snwcp_discarded_liq(indexc)*dtime
write(iulog,*)'qflx_floodc = ',qflx_floodc(indexc)*dtime
write(iulog,*)'qflx_drain_perched = ',qflx_drain_perched(indexc)*dtime
write(iulog,*)'qflx_glcice_dyn_water_flux = ',qflx_glcice_dyn_water_flux(indexc)*dtime
write(iulog,*)'qflx_rootsoi_col(1:nlevsoil) = ',qflx_rootsoi_col(indexc,:)*dtime
write(iulog,*)'total_plant_stored_h2o_col = ',total_plant_stored_h2o_col(indexc)
write(iulog,*)'deltawb = ',endwb(indexc)-begwb(indexc)
write(iulog,*)'deltawb/dtime = ',(endwb(indexc)-begwb(indexc))/dtime
write(iulog,*)'deltaflux = ',forc_rain_col(indexc)+forc_snow_col(indexc) - (qflx_evap_tot(indexc) + &
qflx_surf(indexc)+qflx_h2osfc_surf(indexc)+qflx_drain(indexc))

write(iulog,*)'clm model is stopping'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))

if (use_parflow_soilwater()) then
! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
write(iulog,*)'Ignoring water balance error...'
else
write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
end if
else if (abs(errh2o(indexc)) > 1.e-5_r8 .and. (DAnstep > skip_steps) ) then

write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'

write(iulog,*)'nstep = ',nstep
write(iulog,*)'errh2o = ',errh2o(indexc)
write(iulog,*)'forc_rain = ',forc_rain_col(indexc)*dtime
Expand All @@ -427,9 +436,17 @@ subroutine BalanceCheck( bounds, &
write(iulog,*)'qflx_glcice_dyn_water_flux = ', qflx_glcice_dyn_water_flux(indexc)*dtime
write(iulog,*)'qflx_snwcp_discarded_ice = ',qflx_snwcp_discarded_ice(indexc)*dtime
write(iulog,*)'qflx_snwcp_discarded_liq = ',qflx_snwcp_discarded_liq(indexc)*dtime
write(iulog,*)'qflx_floodc = ',qflx_floodc(indexc)*dtime
write(iulog,*)'qflx_drain_perched = ',qflx_drain_perched(indexc)*dtime
write(iulog,*)'qflx_glcice_dyn_water_flux = ',qflx_glcice_dyn_water_flux(indexc)*dtime
write(iulog,*)'qflx_rootsoi_col(1:nlevsoil) = ',qflx_rootsoi_col(indexc,:)*dtime
write(iulog,*)'clm model is stopping'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
if (use_parflow_soilwater()) then
! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
write(iulog,*)'Ignoring water balance error...'
else
write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
end if
end if
end if

Expand Down Expand Up @@ -676,7 +693,7 @@ subroutine BalanceCheck( bounds, &
end if

! Soil energy balance check

! COUP_OAS_PFL
found = .false.
do c = bounds%begc,bounds%endc
if (col%active(c)) then
Expand All @@ -690,9 +707,12 @@ subroutine BalanceCheck( bounds, &
write(iulog,*)'WARNING: BalanceCheck: soil balance error (W/m2)'
write(iulog,*)'nstep = ',nstep
write(iulog,*)'errsoi_col = ',errsoi_col(indexc)
if (abs(errsoi_col(indexc)) > 1.e-4_r8 .and. (DAnstep > skip_steps) ) then
write(iulog,*)'clm model is stopping'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
if (use_parflow_soilwater()) then
! TODO: Balance errors must be fixed for fully coupled model (ICON-eCLM-ParFlow)
write(iulog,*)'Ignoring soil balance error...'
else
write(iulog,*)'clm model is stopping - error is greater than 1e-5 (mm)'
call endrun(decomp_index=indexc, clmlevel=namec, msg=errmsg(sourcefile, __LINE__))
end if
end if

Expand Down
16 changes: 12 additions & 4 deletions src/clm5/biogeophys/CanopyTemperatureMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,11 +152,12 @@ subroutine CanopyTemperature(bounds, &
qg => waterstate_inst%qg_col , & ! Output: [real(r8) (:) ] ground specific humidity [kg/kg]
qg_h2osfc => waterstate_inst%qg_h2osfc_col , & ! Output: [real(r8) (:) ] specific humidity at h2osfc surface [kg/kg]
dqgdT => waterstate_inst%dqgdT_col , & ! Output: [real(r8) (:) ] d(qg)/dT
pfl_psi => waterstate_inst%pfl_psi_col , & ! Input: [real(r8) (:,:) ] COUP_OAS_PFL

qflx_evap_tot => waterflux_inst%qflx_evap_tot_patch , & ! Output: [real(r8) (:) ] qflx_evap_soi + qflx_evap_can + qflx_tran_veg
qflx_evap_veg => waterflux_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm)
qflx_tran_veg => waterflux_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm)

htvp => energyflux_inst%htvp_col , & ! Output: [real(r8) (:) ] latent heat of vapor of water (or sublimation) [j/kg]
cgrnd => energyflux_inst%cgrnd_patch , & ! Output: [real(r8) (:) ] deriv. of soil energy flux wrt to soil temp [w/m2/k]
cgrnds => energyflux_inst%cgrnds_patch , & ! Output: [real(r8) (:) ] deriv. of soil sensible heat flux wrt soil temp [w/m2/k]
Expand Down Expand Up @@ -197,7 +198,7 @@ subroutine CanopyTemperature(bounds, &
rootr_road_perv => soilstate_inst%rootr_road_perv_col , & ! Input: [real(r8) (:,:) ] effective fraction of roots in each soil layer for urban pervious road
soilalpha => soilstate_inst%soilalpha_col , & ! Output: [real(r8) (:) ] factor that reduces ground saturated specific humidity (-)
soilalpha_u => soilstate_inst%soilalpha_u_col , & ! Output: [real(r8) (:) ] Urban factor that reduces ground saturated specific humidity (-)

t_h2osfc => temperature_inst%t_h2osfc_col , & ! Input: [real(r8) (:) ] surface water temperature
t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin)
beta => temperature_inst%beta_col , & ! Output: [real(r8) (:) ] coefficient of convective velocity [-]
Expand Down Expand Up @@ -256,8 +257,15 @@ subroutine CanopyTemperature(bounds, &
wx = (h2osoi_liq(c,1)/denh2o+h2osoi_ice(c,1)/denice)/dz(c,1)
fac = min(1._r8, wx/watsat(c,1))
fac = max( fac, 0.01_r8 )
psit = -sucsat(c,1) * fac ** (-bsw(c,1))
psit = max(smpmin(c), psit)
! clm3.5/bld/usr.src/Biogeophysics1Mod.F90
! if COUP_OAS_PFL
if (pfl_psi(c,1)>= 0.0_r8) psit = 0._r8
if (pfl_psi(c,1) < 0.0_r8) psit = pfl_psi(c,1)
!else
!psit = -sucsat(c,1) * fac ** (-bsw(c,1))
!psit = max(smpmin(c), psit)
!end if

! modify qred to account for h2osfc
hr = exp(psit/roverg/t_soisno(c,1))
qred = (1._r8 - frac_sno(c) - frac_h2osfc(c))*hr &
Expand Down
57 changes: 24 additions & 33 deletions src/clm5/biogeophys/HydrologyDrainageMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,8 @@ subroutine HydrologyDrainage(bounds, &
use clm_varctl , only : use_vichydro
use clm_varpar , only : nlevgrnd, nlevurb, nlevsoi
use clm_time_manager , only : get_step_size, get_nstep
use SoilHydrologyMod , only : CLMVICMap, Drainage, PerchedLateralFlow, LateralFlowPowerLaw
use SoilWaterMovementMod , only : use_aquifer_layer
use SoilHydrologyMod , only : CLMVICMap, Drainage, PerchedLateralFlow, LateralFlowPowerLaw, ParFlowDrainage
use SoilWaterMovementMod , only : use_aquifer_layer, use_parflow_soilwater
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
Expand Down Expand Up @@ -107,8 +107,6 @@ subroutine HydrologyDrainage(bounds, &
qflx_drain => waterflux_inst%qflx_drain_col , & ! sub-surface runoff (mm H2O /s)
qflx_surf => waterflux_inst%qflx_surf_col , & ! surface runoff (mm H2O /s)
qflx_infl => waterflux_inst%qflx_infl_col , & ! infiltration (mm H2O /s)
qflx_rootsoi => waterflux_inst%qflx_rootsoi_col , & ! root and soil water exchange [mm H2O/s] [+ into root]
qflx_parflow => waterflux_inst%qflx_parflow_col , & ! source/sink flux passed to ParFlow for each soil layer
qflx_qrgwl => waterflux_inst%qflx_qrgwl_col , & ! qflx_surf at glaciers, wetlands, lakes
qflx_runoff => waterflux_inst%qflx_runoff_col , & ! total runoff
! (qflx_drain+qflx_surf+qflx_qrgwl) (mm H2O /s)
Expand All @@ -128,24 +126,29 @@ subroutine HydrologyDrainage(bounds, &
soilhydrology_inst, waterstate_inst)
endif

if (use_aquifer_layer()) then
call Drainage(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
temperature_inst, soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)
if (use_parflow_soilwater()) then
! clm3.5/bld/usr.src/SoilHydrologyMod.F90
! ignore drainage calculations in eCLM and instead pass these fluxes to ParFlow
call ParFlowDrainage(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc, waterflux_inst)
else

call PerchedLateralFlow(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)


call LateralFlowPowerLaw(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)

if (use_aquifer_layer()) then
call Drainage(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
temperature_inst, soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)
else
call PerchedLateralFlow(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)

call LateralFlowPowerLaw(bounds, num_hydrologyc, filter_hydrologyc, &
num_urbanc, filter_urbanc,&
soilhydrology_inst, soilstate_inst, &
waterstate_inst, waterflux_inst)

endif
endif

do j = 1, nlevgrnd
Expand Down Expand Up @@ -219,18 +222,6 @@ subroutine HydrologyDrainage(bounds, &
end if

end do

! Calculate here the source/sink term for ParFlow
do j = 1, nlevsoi
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (j == 1) then
qflx_parflow(c,j) = qflx_infl(c) - qflx_rootsoi(c,j)
else
qflx_parflow(c,j) = -qflx_rootsoi(c,j)
end if
end do
end do
end associate

end subroutine HydrologyDrainage
Expand Down
30 changes: 19 additions & 11 deletions src/clm5/biogeophys/HydrologyNoDrainageMod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ subroutine HydrologyNoDrainage(bounds, &
use SnowHydrologyMod , only : SnowWater, BuildSnowFilter
use SoilHydrologyMod , only : CLMVICMap, SurfaceRunoff, Infiltration, WaterTable, PerchedWaterTable
use SoilHydrologyMod , only : ThetaBasedWaterTable, RenewCondensation
use SoilWaterMovementMod , only : SoilWater
use SoilWaterMovementMod , only : SoilWater, use_parflow_soilwater
use SoilWaterRetentionCurveMod, only : soil_water_retention_curve_type
use SoilWaterMovementMod , only : use_aquifer_layer
use SoilWaterPlantSinkMod , only : Compute_EffecRootFrac_And_VertTranSink
Expand Down Expand Up @@ -156,13 +156,15 @@ subroutine HydrologyNoDrainage(bounds, &
h2osno_top => waterstate_inst%h2osno_top_col , & ! Output: [real(r8) (:) ] mass of snow in top layer (col) [kg]
wf => waterstate_inst%wf_col , & ! Output: [real(r8) (:) ] soil water as frac. of whc for top 0.05 m
wf2 => waterstate_inst%wf2_col , & ! Output: [real(r8) (:) ] soil water as frac. of whc for top 0.17 m
pfl_psi => waterstate_inst%pfl_psi_col , & ! Input: [real(r8) (:,:) ] COUP_OAS_PFL

watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity)
sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm)
bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b"
smp_l => soilstate_inst%smp_l_col , & ! Input: [real(r8) (:,:) ] soil matrix potential [mm]
smpmin => soilstate_inst%smpmin_col , & ! Input: [real(r8) (:) ] restriction for min of soil potential (mm)
soilpsi => soilstate_inst%soilpsi_col & ! Output: [real(r8) (:,:) ] soil water potential in each soil layer (MPa)

)

! Determine step size
Expand Down Expand Up @@ -198,9 +200,12 @@ subroutine HydrologyNoDrainage(bounds, &

if ( use_fates ) call clm_fates%ComputeRootSoilFlux(bounds, num_hydrologyc, filter_hydrologyc, soilstate_inst, waterflux_inst)

! -- clm3.5/src/biogeophys/SoilHydrologyMod.F90
! if not COUP_OAS_PFL
call SoilWater(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filter_urbanc, &
soilhydrology_inst, soilstate_inst, waterflux_inst, waterstate_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
soilhydrology_inst, soilstate_inst, waterflux_inst, waterstate_inst, temperature_inst, &
canopystate_inst, energyflux_inst, soil_water_retention_curve)
! end if

if (use_vichydro) then
! mapping soilmoist from CLM to VIC layers for runoff calculations
Expand Down Expand Up @@ -471,17 +476,20 @@ subroutine HydrologyNoDrainage(bounds, &
! ZMS: Note, this form, which seems to be the same as used in SoilWater, DOES NOT distinguish between
! ice and water volume, in contrast to the soilpsi calculation above. It won't be used in ch4Mod if
! t_soisno <= tfrz, though.
do j = 1, nlevgrnd
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)
if (.not. use_parflow_soilwater()) then
do j = 1, nlevgrnd
do fc = 1, num_hydrologyc
c = filter_hydrologyc(fc)

s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
s_node = min(1.0_r8, s_node)
s_node = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8)
s_node = min(1.0_r8, s_node)

smp_l(c,j) = -sucsat(c,j)*s_node**(-bsw(c,j))
smp_l(c,j) = max(smpmin(c), smp_l(c,j))
smp_l(c,j) = -sucsat(c,j)*s_node**(-bsw(c,j))
smp_l(c,j) = max(smpmin(c), smp_l(c,j))
end do
end do
end do
endif


! if (use_cn) then
! Available soil water up to a depth of 0.05 m.
Expand Down
Loading

0 comments on commit 0b9c215

Please sign in to comment.