Skip to content

Commit 886964b

Browse files
chuckaustinchuckaustin
committed
New routine to calculate surface runoff in polygonal tundra and some small fixes elsewhere
1 parent edb5419 commit 886964b

File tree

3 files changed

+78
-44
lines changed

3 files changed

+78
-44
lines changed

‎components/elm/src/biogeophys/CanopyHydrologyMod.F90

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -795,18 +795,18 @@ subroutine FracH2OSfc(bounds, num_h2osfc, filter_h2osfc, &
795795
!
796796
! !LOCAL VARIABLES:
797797
integer :: c,f,l,k ! indices
798-
real(r8):: d,fd,dfdd ! temporary variable for frac_h2oscs iteration
799-
real(r8):: sigma ! microtopography pdf sigma in mm
800-
real(r8):: swc ! surface water content in m
798+
real(r8):: d,fd,dfdd ! temporary variable for frac_h2oscs iteration
799+
real(r8):: sigma ! microtopography pdf sigma in mm
800+
real(r8):: swc ! surface water content in m
801801
real(r8):: min_h2osfc
802802
!-----------------------------------------------------------------------
803803

804804
associate( &
805805
micro_sigma => col_pp%micro_sigma , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m)
806806

807807
iwp_microrel => col_pp%iwp_microrel , & ! Input: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
808-
iwp_exclvol => col_pp%iwp_exclvol , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m)
809-
iwp_ddep => col_pp%iwp_ddep , & ! Input: [real(r8) (:) ] microtopography pdf sigma (m)
808+
iwp_exclvol => col_pp%iwp_exclvol , & ! Input: [real(r8) (:) ] ice wedge polygon excluded volume (m)
809+
iwp_ddep => col_pp%iwp_ddep , & ! Input: [real(r8) (:) ] ice wedge polygon depression depth (m)
810810

811811
h2osno => col_ws%h2osno , & ! Input: [real(r8) (:) ] snow water (mm H2O)
812812

‎components/elm/src/biogeophys/SoilHydrologyMod.F90

Lines changed: 70 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -201,13 +201,18 @@ subroutine SurfaceRunoff (bounds, num_hydrologyc, filter_hydrologyc, &
201201

202202
do fc = 1, num_hydrologyc
203203
c = filter_hydrologyc(fc)
204-
205-
! assume qinmax large relative to qflx_top_soil in control
206-
if (origflag == 1) then
207-
qflx_surf(c) = fcov(c) * qflx_top_soil(c)
204+
l = col_pp%landunit(c)
205+
! no qflx_surf in polygonal ground
206+
if (lun_pp%ispolygon(l)) then
207+
qflx_surf(c) = 0
208208
else
209-
! only send fast runoff directly to streams
210-
qflx_surf(c) = fsat(c) * qflx_top_soil(c)
209+
! assume qinmax large relative to qflx_top_soil in control
210+
if (origflag == 1) then
211+
qflx_surf(c) = fcov(c) * qflx_top_soil(c)
212+
else
213+
! only send fast runoff directly to streams
214+
qflx_surf(c) = fsat(c) * qflx_top_soil(c)
215+
endif
211216
endif
212217
end do
213218

@@ -267,7 +272,7 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
267272
use elm_varcon , only : denh2o, denice, roverg, wimp, mu, tfrz
268273
use elm_varcon , only : pondmx, watmin
269274
use column_varcon , only : icol_roof, icol_road_imperv, icol_sunwall, icol_shadewall, icol_road_perv
270-
use landunit_varcon , only : istsoil, istcrop
275+
use landunit_varcon , only : istsoil, istcrop, ilowcenpoly
271276
use elm_time_manager , only : get_step_size, get_nstep
272277
use atm2lndType , only : atm2lnd_type ! land river two way coupling
273278
use lnd2atmType , only : lnd2atm_type
@@ -313,6 +318,9 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
313318
real(r8) :: d
314319
real(r8) :: h2osoi_vol
315320
real(r8) :: basis ! temporary, variable soil moisture holding capacity
321+
real(r8) :: vdep ! temporary, ice wedge polygon volumetric depression depth (m)
322+
real(r8) :: phi_eff ! temporary, polygonal ground effective subsidence (maxes out at 0.4) (m)
323+
real(r8) :: swc ! temporary, polygonal surface water content in m
316324
! in top VIC layers for runoff calculation
317325
real(r8) :: rsurf_vic ! temp VIC surface runoff
318326
real(r8) :: top_moist(bounds%begc:bounds%endc) ! temporary, soil moisture in top VIC layers
@@ -328,7 +336,11 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
328336
dz => col_pp%dz , & ! Input: [real(r8) (:,:) ] layer depth (m)
329337
nlev2bed => col_pp%nlevbed , & ! Input: [integer (:) ] number of layers to bedrock
330338
cgridcell => col_pp%gridcell , & ! Input: [integer (:) ] column's gridcell
331-
wtgcell => col_pp%wtgcell , & ! Input: [real(r8) (:) ] weight (relative to gridcell)
339+
wtgcell => col_pp%wtgcell , & ! Input: [real(r8) (:) ] weight (relative to gridcell)
340+
iwp_microrel => col_pp%iwp_microrel , & ! Input: [real(r8) (:) ] ice wedge polygon microtopographic relief (m)
341+
iwp_exclvol => col_pp%iwp_exclvol , & ! Input: [real(r8) (:) ] ice wedge polygon excluded volume (m)
342+
iwp_ddep => col_pp%iwp_ddep , & ! Input: [real(r8) (:) ] ice wedge polygon depression depth (m)
343+
meangradz => col_pp%meangradz , & ! Input: [real(r8) (:) ] mean topographic gradient at the column level (unitless)
332344

333345
t_soisno => col_es%t_soisno , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin)
334346

@@ -514,15 +526,36 @@ subroutine Infiltration(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, f
514526
endif
515527
endif
516528

517-
! limit runoff to value of storage above S(pc)
518-
if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
519-
! spatially variable k_wet
520-
k_wet=1.0_r8 * sin((rpi/180.) * col_pp%topo_slope(c))
521-
qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c))
522-
523-
qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime)
529+
if (lun_pp%ispolygon(col_pp%landunit(c))) then
530+
vdep = (2_r8*iwp_exclvol(c) - iwp_microrel(c)) * (iwp_ddep(c)/iwp_microrel(c))**3_r8 &
531+
+ (2_r8*iwp_microrel(c) - 3_r8*iwp_exclvol(c)) * (iwp_ddep(c)/iwp_microrel(c))**2_r8
532+
phi_eff = min(subsidence, 0.4) !fix this variable when available to pull from alt calculations
533+
swc = h2osfc(c)/1000_r8 ! convert to m
534+
535+
if (swc >= vdep) then
536+
if (lun_pp%polygontype(col_pp%landunit(c)) == ilowcenpoly) then
537+
k_wet = (2890_r8*phi_eff**4 - 1171.1_r8*phi_eff**3 + 144.94_r8*phi_eff**2 + 1.682_r8*phi_eff + 2.028) &
538+
* (710.3_r8*meangradz(c)**2 - 28.736_r8*meangradz(c) + 12.74_r8)
539+
else
540+
k_wet = 24.925_r8 * (710.3_r8*meangradz(c)**2 - 28.736_r8*meangradz(c) + 12.74_r8)
541+
endif
542+
qflx_h2osfc_surf(c) = k_wet * (swc - vdep)
543+
qflx_h2osfc_surf(c) = min(qflx_h2osfc_surf(c), (swc - vdep)*1000_r8/dtime)
544+
else
545+
qflx_h2osfc_surf(c) = 0._r8
546+
endif
547+
524548
else
525-
qflx_h2osfc_surf(c)= 0._r8
549+
! limit runoff to value of storage above S(pc)
550+
if(h2osfc(c) >= h2osfc_thresh(c) .and. h2osfcflag/=0) then
551+
! spatially variable k_wet
552+
k_wet=1.0_r8 * sin((rpi/180.) * col_pp%topo_slope(c))
553+
qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c))
554+
555+
qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime)
556+
else
557+
qflx_h2osfc_surf(c)= 0._r8
558+
endif
526559
endif
527560

528561
! cutoff lower limit
@@ -696,7 +729,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
696729
real(r8) :: q_perch
697730
real(r8) :: q_perch_max
698731
real(r8) :: dflag=0._r8
699-
real(r8) :: qcharge_temp
732+
real(r8) :: qcharge_temp
700733
!-----------------------------------------------------------------------
701734

702735
associate( &
@@ -783,23 +816,23 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
783816
! Water table changes due to qcharge
784817
do fc = 1, num_hydrologyc
785818
c = filter_hydrologyc(fc)
786-
nlevbed = nlev2bed(c)
819+
nlevbed = nlev2bed(c)
787820

788821
!scs: use analytical expression for aquifer specific yield
789822
rous = watsat(c,nlevbed) &
790823
* ( 1. - (1.+1.e3*zwt(c)/sucsat(c,nlevbed))**(-1./bsw(c,nlevbed)))
791824
rous=max(rous,0.02_r8)
792825

793826
!-- water table is below the soil column --------------------------------------
794-
g = col_pp%gridcell(c)
827+
g = col_pp%gridcell(c)
795828
l = col_pp%landunit(c)
796829
qcharge_temp = qcharge(c)
797830

798831
wa(c) = wa(c) - qflx_grnd_irrig_col(c) * dtime
799832
zwt(c) = zwt(c) + (qflx_grnd_irrig_col(c) * dtime)/1000._r8/rous
800833

801834
if(jwt(c) == nlevbed) then
802-
if (.not. (zengdecker_2009_with_var_soil_thick)) then
835+
if (.not. (zengdecker_2009_with_var_soil_thick)) then
803836
wa(c) = wa(c) + qcharge(c) * dtime
804837
zwt(c) = zwt(c) - (qcharge(c) * dtime)/1000._r8/rous
805838
end if
@@ -865,7 +898,7 @@ subroutine WaterTable(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, fil
865898
! perched water table code
866899
do fc = 1, num_hydrologyc
867900
c = filter_hydrologyc(fc)
868-
nlevbed = nlev2bed(c)
901+
nlevbed = nlev2bed(c)
869902

870903
! define frost table as first frozen layer with unfrozen layer above it
871904
if(t_soisno(c,1) > tfrz) then
@@ -1132,7 +1165,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
11321165

11331166
do fc = 1, num_hydrologyc
11341167
c = filter_hydrologyc(fc)
1135-
nlevbed = nlev2bed(c)
1168+
nlevbed = nlev2bed(c)
11361169
jwt(c) = nlevbed
11371170
! allow jwt to equal zero when zwt is in top layer
11381171
do j = 1,nlevbed
@@ -1142,7 +1175,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
11421175
else
11431176
jwt(c) = j-1
11441177
exit
1145-
end if
1178+
end if
11461179
end if
11471180
enddo
11481181
end do
@@ -1153,7 +1186,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
11531186
! perched water table code
11541187
do fc = 1, num_hydrologyc
11551188
c = filter_hydrologyc(fc)
1156-
nlevbed = nlev2bed(c)
1189+
nlevbed = nlev2bed(c)
11571190

11581191
! specify maximum drainage rate
11591192
q_perch_max = 1.e-5_r8 * sin(col_pp%topo_slope(c) * (rpi/180._r8))
@@ -1357,11 +1390,11 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
13571390
! make sure baseflow isn't negative
13581391
rsub_top(c) = max(0._r8, rsub_top(c))
13591392
else
1360-
if (jwt(c) == nlevbed .and. zengdecker_2009_with_var_soil_thick) then
1393+
if (jwt(c) == nlevbed .and. zengdecker_2009_with_var_soil_thick) then
13611394
rsub_top(c) = 0._r8
13621395
else
13631396
rsub_top(c) = imped * rsub_top_max* exp(-fff(c)*zwt(c))
1364-
end if
1397+
end if
13651398
end if
13661399

13671400
if (use_vsfm) rsub_top(c) = 0._r8
@@ -1373,10 +1406,10 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
13731406

13741407
!-- water table is below the soil column --------------------------------------
13751408
if(jwt(c) == nlevbed) then
1376-
if (zengdecker_2009_with_var_soil_thick) then
1377-
if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
1378-
zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
1379-
end if
1409+
if (zengdecker_2009_with_var_soil_thick) then
1410+
if (-1._r8 * smp_l(c,nlevbed) < 0.5_r8 * dzmm(c,nlevbed)) then
1411+
zwt(c) = z(c,nlevbed) - (smp_l(c,nlevbed) / 1000._r8)
1412+
end if
13801413
rsub_top(c) = imped * rsub_top_max * exp(-fff(c) * zwt(c))
13811414
rsub_top_tot = - rsub_top(c) * dtime
13821415
s_y = watsat(c,nlevbed) &
@@ -1391,8 +1424,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
13911424
else
13921425
zwt(c) = zi(c,nlevbed)
13931426
end if
1394-
if (rsub_top_tot < 0.) then
1395-
rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
1427+
if (rsub_top_tot < 0.) then
1428+
rsub_top(c) = rsub_top(c) + rsub_top_tot / dtime
13961429
rsub_top_tot = 0.
13971430
end if
13981431
else
@@ -1488,7 +1521,7 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
14881521

14891522
do fc = 1, num_hydrologyc
14901523
c = filter_hydrologyc(fc)
1491-
nlevbed = nlev2bed(c)
1524+
nlevbed = nlev2bed(c)
14921525
do j = nlevbed,2,-1
14931526
xsi(c) = max(h2osoi_liq(c,j)-eff_porosity(c,j)*dzmm(c,j),0._r8)
14941527
if (use_vsfm) then
@@ -1536,8 +1569,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
15361569

15371570
do fc = 1, num_hydrologyc
15381571
c = filter_hydrologyc(fc)
1539-
nlevbed = nlev2bed(c)
1540-
do j = 1, nlevbed-1
1572+
nlevbed = nlev2bed(c)
1573+
do j = 1, nlevbed-1
15411574
if (h2osoi_liq(c,j) < watmin) then
15421575
xs(c) = watmin - h2osoi_liq(c,j)
15431576
! deepen water table if water is passed from below zwt layer
@@ -1555,8 +1588,8 @@ subroutine Drainage(bounds, num_hydrologyc, filter_hydrologyc, num_urbanc, filte
15551588
! Get water for bottom layer from layers above if possible
15561589
do fc = 1, num_hydrologyc
15571590
c = filter_hydrologyc(fc)
1558-
nlevbed = nlev2bed(c)
1559-
j = nlevbed
1591+
nlevbed = nlev2bed(c)
1592+
j = nlevbed
15601593
if (h2osoi_liq(c,j) < watmin) then
15611594
xs(c) = watmin-h2osoi_liq(c,j)
15621595
searchforwater: do i = nlevbed-1, 1, -1

‎components/elm/src/main/landunit_varcon.F90

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,8 +38,8 @@ module landunit_varcon
3838

3939
! land unit polygonal ground types
4040
integer, parameter, public :: ilowcenpoly = 1 ! low-centered polygons
41-
integer, parameter, public :: iflatcenpoly = 2 ! flat-centered polygons
42-
integer, parameter, public :: ihighcenpoly = 3 ! high-centered polygons
41+
integer, parameter, public :: iflatcenpoly = 2 ! flat-centered polygons
42+
integer, parameter, public :: ihighcenpoly = 3 ! high-centered polygons
4343

4444
integer, parameter, public :: max_poly = 3 !maximum value that lun_pp%polygontype can have
4545
!(i.e., largest value in the above list)
@@ -82,6 +82,7 @@ subroutine landunit_varcon_init()
8282
!-----------------------------------------------------------------------
8383

8484
call set_landunit_names()
85+
call set_polygon_names()
8586

8687
end subroutine landunit_varcon_init
8788

0 commit comments

Comments
 (0)