c=======================================================================
c
SUBROUTINE inter2
c
c=======================================================================
c
c calculation of interception and drainage of rainfall and snow
c incorporating effects of patchy snow cover and temperature
c adjustments.
c
c----------------------------------------------------------------------
c
c (1) non-uniform precipitation
c convective ppn. is described by area-intensity
c relationship :-
c
c f(x) = a*exp(-b*x)+c
c
c throughfall, interception and infiltration
c excess are functional on this relationship
c and proportion of large-scale ppn.
c reference: sato et al.(1989b), appendix.
c
c (2) reorganisation of snowmelt and rain-freeze procedures.
c subroutine adjust
c
c (3) additional calculation for partial snow-cover case.
c subroutine patchs
c
c (4) reorganisation of overland flow.
c reference: SA-89B, appendix.
c
c (5) modified calaculation of soil heat capacity and
c conductivity.
c
c=======================================================================
c
c subroutines in this block : snow1
c ------------------------- adjust
c patchs
c snow1
c
c++++++++++++++++++++++++++++++output+++++++++++++++++++++++++++++++++++
c
c roff runoff (mm)
c tc canopy temperature (K)
c tg ground surface temperature (K)
c www(1) ground wetness of surface layer
c capac(2) canopy/ground liquid interception store (m)
c snoww(2) canopy/ground snow interception store (m)
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
INCLUDE 'COMSIBC.H'
c
real pcoefs(2,2)
c data pcoefs(1,1)/ 20. /, pcoefs(1,2)/ .206e-8 /,
c & pcoefs(2,1)/ 0.0001 /, pcoefs(2,2)/ 0.9999 /, bp /20. /
c! data pcoefs(1,1)/ 25. /, pcoefs(1,2)/ 1.38879E-11 /,
c! & pcoefs(2,1)/ 25. /, pcoefs(2,2)/ 1.38879E-11 /, bp /25. /
c! tatsch
data pcoefs(1,1)/ 0.0001 /, pcoefs(1,2)/ 0.9999 /,
& pcoefs(2,1)/ 0.0001 /, pcoefs(2,2)/ 0.9999 /, bp /20. /
c
c
c-----------------------------------------------------------------------
c
CALL snow1
c
c-----------------------------------------------------------------------
c
c prec ( pi-x ) : equation (c.3), SA-89B
c
c-----------------------------------------------------------------------
c
pcoefs(1,1) = app
pcoefs(2,1) = app
pcoefs(1,2) = cpp
pcoefs(2,2) = cpp
bp = bpp
ap = pcoefs(2,1)
cp = pcoefs(2,2)
totalp = ppc + ppl
IF( snoww(1) .gt. 0. .or. snoww(2) .gt. 0. .or. tm .lt. tf )
& ppc = 0.
ppl = totalp - ppc
IF(totalp.ge.1.e-8) THEN
ap = ppc/totalp * pcoefs(1,1) + ppl/totalp * pcoefs(2,1)
cp = ppc/totalp * pcoefs(1,2) + ppl/totalp * pcoefs(2,2)
END IF
c
roff = 0.
roff1 = 0.
roff2 = 0.
roffGA = 0. ! tatsch
roff3 = 0.
roffq3g = 0.
roffq3g2 = 0.
roff4 = 0.
gwsoil = 0.
thru = 0.
fpi = 0.
finfil = 0.
c
c----------------------------------------------------------------------
c heat capacity of the soil, as used in force-restore heat flux
c description. dependence of csoil on porosity and wetness is
c based on CS-81.
c----------------------------------------------------------------------
c
slamda = ( 1.5*(1.-poros) + 1.3*www(1)*poros ) /
& ( 0.75 + 0.65*poros - 0.4*www(1)*poros ) * 0.4186
shcap = ( 0.5*(1.-poros) + www(1)*poros ) * 4.186 * 1.e6
csoil = sqrt( slamda * shcap * 86400./pie ) / 2.
c
c----------------------------------------------------------------------
c input precipitation is given in mm, converted to m to give p0.
c----------------------------------------------------------------------
c
p0 = totalp * 0.001
c
DO iveg = 1, 2
c
realc = 2. - iveg
realg = iveg - 1.
c
xsc = amax1(0., capac(iveg) - satcap(iveg) )
capac(iveg) = capac(iveg) - xsc
xss = amax1(0., snoww(iveg) - satcap(iveg) ) * realc
snoww(iveg) = snoww(iveg) - xss
if (iveg.eq.1) then !Because freezing, snoww(1)>satcap, put to snoww(2)
snoww(2) =snoww(2) + xss
xss = 0.0
endif
p0 = p0 + xsc + xss !tang@2005/12/21
c
capacp = capac(iveg)
snowwp = snoww(iveg)
c
spechc = amin1( 0.05, ( capac(iveg) + snoww(iveg) ) ) * cw
& + realc * zlt * clai + realg * csoil
ts = tc * realc + tg * realg
c
c----------------------------------------------------------------------
c proportional saturated area (xs) and leaf drainage(tex)
c
c tex ( d-c ) : equation (c.8), SA-89B
c xs ( x-s ) : equation (c.7), SA-89B
c tex ( d-c ) : equation (c.8), SA-89B
c
c-----------------------------------------------------------------------
c
chiv = chil ! Leaf area distribution factor
IF ( abs(chiv) .le. 0.01 ) chiv = 0.01
aa = 0.5 - 0.633 * chiv - 0.33 * chiv * chiv
bb = 0.877 * ( 1. - 2. * aa )
exrain = aa + bb
c
zload = capac(iveg) + snoww(iveg)
fpi = ( 1.-exp( - exrain*zlt/vcover ) )*vcover*realc + realg
tti = p0 * ( 1.-fpi )
xs = 1.
IF ( p0 .ge. 1.e-9 ) THEN
arg = ( satcap(iveg)-zload )/( p0*fpi*ap ) -cp/ap
IF ( arg .ge. 1.e-9 ) THEN
xs = -1./bp * alog( arg )
xs = amin1( xs, 1. )
xs = amax1( xs, 0. )
END IF
END IF
Pxi = amax1(0.0,satcap(iveg) - zload)
tex = p0*fpi * ( ap/bp*( 1.- exp( -bp*xs )) + cp*xs ) -
& ( Pxi ) * xs
tex = amax1( tex, 0. )
c
c----------------------------------------------------------------------
c total throughfall (thru) and store augmentation
c----------------------------------------------------------------------
c
IF ( iveg .eq. 1 ) THEN
c
thru = tti + tex
pinf = p0 - thru
IF( tm .gt. tf ) capac(iveg) = capac(iveg) + pinf
IF( tm .le. tf ) snoww(iveg) = snoww(iveg) + pinf
c
CALL adjust ( tc, spechc, capacp, snowwp, iveg )
c
p0 = thru
c
ELSE IF ( tg .gt. tf .and. snoww(2) .gt. 0. ) THEN
CALL patchs ( p0 )
ELSE
thru = tti + tex
IF ( tg .le. tf .or. tm .le. tf ) thru = 0.
pinf = p0 - thru
IF( tm .gt. tf ) capac(iveg) = capac(iveg) + pinf
IF( tm .le. tf ) snoww(iveg) = snoww(iveg) + pinf
IF( tm .gt. tf ) THEN
c www(1) = www(1) + thru/ ( poros*zdepth(1) )
finfil = finfil + thru
END IF
c
CALL adjust ( tg, spechc, capacp, snowwp, iveg )
c
END IF
END DO
c
c----------------------------------------------------------------------
c
c instantaneous overland flow contribution ( roff )
c roff( r-i ) : equation (c.13), SA-89B
c finfil : through fall and snow melt etc which reach the soil top surface
c-----------------------------------------------------------------------
xw1 = amax1(0., www(1) - 1.0 )
www(1) = www(1) - xw1
finfil = finfil + xw1 * ( poros*zdepth(1) )
c----------------------------------------------------------------------
c if froze, reduce infiltration
c-----------------------------------------------------------------------
tsnow = amin1 ( tf-0.01, tg )
tgs = tsnow*areas + tg*(1.-areas)
props = ( tgs-(tf-1.) ) / 1.
props = amax1( 0.0, amin1( 1.0, props ) )
if (idirr.ne.1) then
ap = app
bp = bpp
cp = cpp
c! equdep = satco * dtt * props
equdep = satcoz(1) * dtt * props
xs = 1.
IF( finfil .ge. 1.e-9 ) THEN
arg = equdep / ( finfil * ap ) -cp/ap
IF( arg .ge. 1.e-9 ) THEN
xs = -1./bp * alog( arg )
xs = amin1( xs, 1. )
xs = amax1( xs, 0. )
END IF
END IF
roffo = finfil * ( ap/bp * ( 1.-exp( -bp*xs )) + cp*xs )
& - equdep*xs
roffo = amax1 ( roffo, 0. )
roff = roff + roffo
roff1= roff1+ roffo
finfil = finfil - roffo
endif
c if (finfil.gt.0.) then
c xw1 = amin1(finfil, (1.0-www(1) )* ( poros*zdepth(1) ) )
c www(1) = www(1) + xw1 / ( poros*zdepth(1) )
c finfil = finfil -xw1
c endif
CALL Green_AMPT(finfil)
c thru_re= finfil
roff = roff + finfil
roff4= roff4+ finfil
finfil=0.
RETURN
END