Program drain_m
c Vaughan Voller
c University of Minnesota
c October 2002
c program for calculating drainage from an edge drain
c drain program
implicit none
Integer ktime, Kit,ntime
Integer I
integer subbase !material index of subbase
Integer Tri(2999,0:3)
Integer Sup(1599,0:8)
Integer Bound(0:150)
Integer props(1599,8)
Integer Ntri
Integer Nodes
Real Ctri(2999,3,4)
Real Coef(1599,0:8)
Real Area(2999)
Real Narea(1599)
Real NodeA(1599,0:8)
Real Cp(1599)
Real Source(1599)
Real Bval(150),Bcoef(150)
Real P(1599)
Real Pstar(1599)
Real F(1599)
Real Fold(1599)
Real M_D(1599) !nodal mostures in granular base
Real maxF
Real X(1599)
Real Y(1599)
Real fsum,fsumold !tot. moisture in gran. base
Real fsumin !initial moisture in
Real flowout !flow into drain
Real mprops(5,10)
Real pinit
Real delt
Call Readit (Nodes,Ntri,Bound,Bval,Bcoef,Sup,Tri,X,Y,mprops,
& ntime,delt,pinit,subbase)
Call TriCoef(Ntri,Ctri,Tri,Area,Narea,X,Y)
Call Nodeprops(Nodes,Ntri,Tri,Area,Narea,props,NodeA)
Call initial(Nodes,F,Fold,M_D,P,Pstar,pinit)
Call Moist(Ntri,Tri,Area,Narea,P,F,mprops)
Call Moist(Ntri,Tri,Area,Narea,P,Fold,mprops)
C initial moisture
fsum=0
Call Dmoist(Ntri,Tri,Area,P,M_D,mprops,subbase)
Do I=1,Nodes
fsum=fsum + M_D(I)
enddo
fsumin=fsum
open(4,file='dtime.dat')
write(4,*) 0,fsum/fsumin,fsum,0,0
print *, 'PROGRAM RUNNING please wait'
print *,
print *,'time(sec) -- degree of Saturation-cm^2 moist gran base
&convergence~e-4'
print *, ktime*delt,fsum/fsumin,fsum,' ',maxF
Do ktime = 1,ntime
Do Kit = 1,50
Call Coeficient(Ntri,delt,Ctri,Sup,Tri,Cp,Source,Area,Narea,
& Bound,Bval,Bcoef,Coef,F,Fold,P,mprops)
Call Solve(Nodes,Sup,Source,Coef,Pstar,Bound,Bcoef,flowout)
Call AdjustF(Nodes,F,P,Pstar,Cp,props,mprops,NodeA,maxF)
Call Clean(Nodes,Sup,Source,Cp,Coef,M_D)
If (maxF.lt.0.00005.and.kit.gt.3) goto 10
enddo
10 continue
c Drain Moisture
fsumold=fsum
fsum=0
Call Dmoist(Ntri,Tri,Area,P,M_D,mprops,subbase)
Do I=1,Nodes
fsum=fsum + M_D(I)
enddo
print *, ktime*delt,fsum/fsumin,fsum,' ',maxF
write(4,*) ktime *delt,fsum/fsumin,fsum,flowout,
& fsumold-fsum-flowout
c massbal write(4,*) ktime * delt,fsum-fsumin,flowout
Call Swap(Nodes,F,Fold)
enddo
11 close(4)
c Trench moisture for writeup
Call Trench(Nodes,X,Y,P,mprops)
print *, "Program Sucesfully Terminated"
pause
stop
end
c--------------------------
Subroutine Readit(Nodes,Ntri,Bound,Bval,Bcoef,
& Sup,Tri,X,Y,mprops,ntime,delt,pinit,subbase)
Integer I,J
Integer Tri(2999,0:3)
Integer Sup(1599,0:8)
Integer Bound(0:150)
Integer Ntri
Integer Nodes
Integer nmats
Integer subbase
Integer ntime
Real Bval(150),Bcoef(150)
Real X(1599)
Real Y(1599)
Real mprops(5,10)
Real pinit
Real delt
open(3,file='drainin.dat')
read(3,*) Nodes,Ntri
do I=1,Nodes
read(3,*) Sup(I,0)
read(3,*) (Sup(I,J),J=1,Sup(I,0))
enddo
do I=1,Ntri
read(3,*) (Tri(I,J),J=0,3)
enddo
do I=1,Nodes
read(3,*) X(i),Y(i)
enddo
read(3,*) Bound(0)
read(3,*) (Bound(J),J=1,Bound(0))
read(3,*) (Bcoef(j),J=1,Bound(0))
read(3,*) (Bval(j),J=1,Bound(0))
read (3,*) nmats
do i=1,nmats
read(3,*) (mprops(i,j),j=1,6)
enddo
read(3,*) subbase
read(3,*) ntime
read(3,*) delt
read(3,*) pinit
return
end
c----------------------------------
Subroutine TriCoef(Ntri,Ctri,Tri,Area,Narea,X,Y)
Integer Tri(2999,0:3)
Integer Ntri
Real Ctri(2999,3,4)
Real Area(2999)
Real Narea(1599)
Real X(1599)
Real Y(1599)
Real XT(3)
Real YT(3)
Integer I,J
Integer Nod(3)
real SFX1
real SFX2
real SFX3
real SFY1
real SFY2
real SFY3
real DELX
real DELY
Do I=1,Ntri
Do J=1,3
Nod(J)=Tri(I, J)
XT(J)=X(Nod(J))
YT(J)=Y(Nod(J))
enddo
Area(I) = (XT(2)*YT(3)-XT(3)*YT(2)-XT(1)*(YT(3)-YT(2))+
& YT(1)*(XT(3)-XT(2))) /2.
C Nodal Areas
DoJ=1,3
Narea(Nod(J)) = Narea(Nod(J))+Area(I)/3.0
enddo
C Derivatives of shape function
SFX1 = (YT(2) - YT(3)) / (2*Area(I))
SFX2 = (YT(3) - YT(1)) / (2*Area(I))
SFX3 = (YT(1) - YT(2)) / (2*Area(I))
SFY1 = (XT(3) - XT(2)) / (2*Area(I))
SFY2 = (XT(1) - XT(3)) / (2*Area(I))
SFY3 = (XT(2) - XT(1)) / (2*Area(I))
C Length of bisector 1
C Mid point of element - mid point of side 1-2
DELX = (2*XT(3) - XT(1) - XT(2))/6.
DELY = (2*YT(3) - YT(1) - YT(2))/6.
Ctri(I, 1, 1) = SFX1 * DELY - SFY1 * DELX
Ctri(I, 1, 2) = SFX2 * DELY - SFY2 * DELX
Ctri(I, 1, 3) = SFX3 * DELY - SFY3 * DELX
Ctri(I, 1, 4) = -DELX
C Length of bisector 2
C Mid point of element - mid point of side 2-3
DELX = (2*XT(1) - XT(2) - XT(3))/6.
DELY = (2*YT(1) - YT(2) - YT(3))/6.
Ctri(I, 2, 1) = SFX1 * DELY - SFY1 * DELX
Ctri(I, 2, 2) = SFX2 * DELY - SFY2 * DELX
Ctri(I, 2, 3) = SFX3 * DELY - SFY3 * DELX
Ctri(I, 2, 4) = -DELX
C Length of bisector 3
C Mid point of element - mid point of side 3-1
DELX = (2*XT(2) - XT(1) - XT(3))/6.
DELY = (2*YT(2) - YT(1) - YT(3))/6.
Ctri(I, 3, 1) = SFX1 * DELY - SFY1 * DELX
Ctri(I, 3, 2) = SFX2 * DELY - SFY2 * DELX
Ctri(I, 3, 3) = SFX3 * DELY - SFY3 * DELX
Ctri(I, 3, 4) = -DELX
enddo
return
end
c------------------------------------------------
Subroutine Nodeprops(Nodes,Ntri,Tri,Area,Narea,props,NodeA)
Integer Tri(2999,0:3)
Integer props(1599,8)
Integer Ntri
Integer Nodes
Real Area(2999)
Real Narea(1599)
Real NodeA(1599,0:8)
Integer I,K,J
logical Tricon
Do I = 1,Nodes
Do K = 1,Ntri
Tricon = .False.
Do J = 1,3
If(I.eq.Tri(K, J)) Tricon = .True.
enddo
If(Tricon) Then
NodeA(I, 0) = NodeA(I, 0) + 1
NodeA(I, NodeA(I, 0)) = (Area(K)/3.) / Narea(I)
props(I, NodeA(I, 0)) = Tri(K, 0)
endif
enddo
enddo
return
end
c----------------------------------------------------
Subroutine initial(Nodes,F,Fold,M_D,P,Pstar,pinit)
Integer I
Integer Nodes
real P(1599)
real Pstar(1599)
real F(1599)
real Fold(1599)
real M_D(1599)
real pinit
Do I = 1,Nodes
P(I) = pinit
Pstar(I) = pinit
F(i)=0.0
Fold(i)=0.0
M_D(i)=0.0
enddo
return
end
c----------------------------------------------------
Subroutine Dmoist(Ntri,Tri,Area,P,M_D,mprops,subbase)
Integer Tri(2999,0:3)
Integer Ntri
Real Area(2999)
Real P(1599)
Real M_D(1599)
Real mprops(5,10)
Integer I,J
Integer Nod(3)
Integer prop
Integer subbase
Do I = 1,Ntri
prop = Tri(I, 0)
if(prop.ne.subbase) then
Do J = 1,3
Nod(J) = Tri(I, J)
M_D(Nod(J))=M_D(Nod(J))+(FNM(P(Nod(J)),prop,mprops)*Area(I)/3.)
enddo
endif
enddo
return
end
c--------------------------------------------
Subroutine Moist(Ntri,Tri,Area,Narea,P,F,mprops)
Integer Tri(2999,0:3)
Integer Ntri
Real Area(2999)
Real Narea(1599)
Real P(1599)
Real F(1599)
Real mprops(5,10)
Integer I,J
Integer Nod(3)
Integer prop
Do I = 1,Ntri
prop = Tri(I, 0)
Do J = 1,3
Nod(J) = Tri(I, J)
F(Nod(J))=F(Nod(J))+(FNM(P(Nod(J)),prop,mprops)*Area(I)/3.)
& /Narea(Nod(J))
enddo
enddo
return
end
c--------------------------------------------
Subroutine Coeficient(Ntri,delt,Ctri,Sup,Tri,Cp,Source,Area,
& Narea,Bound,Bval,Bcoef,Coef,F,Fold,P,mprops)
Integer Tri(2999,0:3)
Integer Sup(1599,0:8)
Integer Bound(0:150)
Integer Ntri
Real Ctri(2999,3,4)
Real Coef(1599,0:8)
Real Area(2999)
Real Narea(1599)
Real Cp(1599)
Real Source(1599)
Real Bval(150)
real Bcoef(150)
Real F(1599)
real P(1599)
Real Fold(1599)
Real mprops(5,10)
Real delt
Integer I, J
Integer Npos
Integer Nod(3)
c Conductivities
real K12
real K23
real K31
real Kave
c Specific Moisture
real C(3)
c property
integer prop
Do I = 1,Ntri
prop = Tri(I, 0)
Do J = 1,3
Nod(J) = Tri(I, J)
C(J) = FNC(P(Nod(J)),prop,mprops)
enddo
K12 = delt * FNK(P(Nod(1)), P(Nod(2)), prop,mprops)
K23 = delt * FNK(P(Nod(2)), P(Nod(3)), prop,mprops)
K31 = delt * FNK(P(Nod(3)), P(Nod(1)), prop,mprops)
c assumming TRAP integration
Kave = (K12 + K23 + K31)/3.
c Contributions to node 1
Do J = 1,Sup(Nod(1), 0)
If(Nod(2).eq.Sup(Nod(1), J)) Npos = J
enddo
Coef(Nod(1),Npos)=Coef(Nod(1),Npos)+K12*(Ctri(I,1,2)-Ctri(I,3,2))
Coef(Nod(1),0)=Coef(Nod(1),0)+K12*(Ctri(I,1,2)-Ctri(I,3,2))
Do J = 1,Sup(Nod(1), 0)
If (Nod(3).eq.Sup(Nod(1), J)) Npos = J
enddo
Coef(Nod(1),Npos)=Coef(Nod(1),Npos)+K31*(Ctri(I,1,3)-Ctri(I,3,3))
Coef(Nod(1),0)=Coef(Nod(1),0)+K31*(Ctri(I,1,3)-Ctri(I,3,3))
Source(Nod(1)) =Source(Nod(1))+Kave*(Ctri(I,1,4)-Ctri(I,3,4))
c Node 2
Do J = 1,Sup(Nod(2), 0)
If (Nod(1).eq.Sup(Nod(2), J)) Npos = J
enddo
Coef(Nod(2),Npos)=Coef(Nod(2),Npos)+K12*(Ctri(I,2,1)-Ctri(I,1,1))
Coef(Nod(2),0)=Coef(Nod(2),0)+K12*(Ctri(I,2,1)-Ctri(I, 1, 1))
Do J = 1,Sup(Nod(2), 0)
If (Nod(3).eq.Sup(Nod(2), J)) Npos = J
enddo
Coef(Nod(2),Npos)=Coef(Nod(2), Npos)+K23*(Ctri(I,2,3)-Ctri(I,1,3))
Coef(Nod(2),0)=Coef(Nod(2),0)+K23*(Ctri(I,2,3)-Ctri(I,1,3))
Source(Nod(2))=Source(Nod(2))+Kave*(Ctri(I,2,4)-Ctri(I,1,4))
c Node 3
Do J = 1,Sup(Nod(3), 0)
If (Nod(1).eq.Sup(Nod(3), J)) Npos = J
enddo
Coef(Nod(3),Npos)=Coef(Nod(3),Npos)+K31*(Ctri(I,3,1)-Ctri(I,2,1))
Coef(Nod(3), 0)=Coef(Nod(3),0)+K31*(Ctri(I,3,1)-Ctri(I,2,1))
Do J = 1,Sup(Nod(3), 0)
If (Nod(2).eq.Sup(Nod(3), J)) Npos = J
enddo
Coef(Nod(3),Npos)=Coef(Nod(3),Npos)+K23*(Ctri(I,3,2)-Ctri(I,2,2))
Coef(Nod(3),0)=Coef(Nod(3),0)+K23*(Ctri(I,3,2)-Ctri(I,2,2))
Source(Nod(3))=Source(Nod(3))+Kave*(Ctri(I,3,4)-Ctri(I,2,4))
C Sources
Do J = 1,3
Coef(Nod(J), 0) = Coef(Nod(J), 0) + C(J) * Area(I)/3.
Source(Nod(J))=Source(Nod(J))+(Area(I)/3.)*(C(J)*P(Nod(J))+
& Fold(Nod(J)) - F(Nod(J)))
Cp(Nod(J))=Cp(Nod(J))+C(J)*(Area(I)/3.)/Narea(Nod(J))
enddo
enddo
C update boundry coeficients
Do I = 1,Bound(0)
Coef(Bound(I), 0) =Coef(Bound(I),0)+Bcoef(I)
Source(Bound(I)) = Source(Bound(I))+Bval(I)
enddo
return
end
c----------------------------------
Subroutine Solve(Nodes,Sup,Source,Coef,P,Bound,Bcoef,flowout)
Integer Sup(1599,0:8)
Integer Nodes
Integer Bound(0:150)
Real Coef(1599,0:8)
Real Source(1599)
Real P(1599)
Real Bcoef(150)
Integer K,I,J
real RHS
real flowout
real flownode
Do K = 1,50
Do I = 1,Nodes
RHS = 0
Do J = 1,Sup(I, 0)
RHS = RHS + Coef(I, J) * P(Sup(I, J))
enddo
P(I) = P(I) + 1.0 * ((RHS + Source(I))/Coef(I, 0) - P(I))
enddo
enddo
c---- calculation of flow into drain ASSUMMES THAT p = 0 in drain
c--- This will need modification to be a full review boundary
flowout=0
do K=1,Bound(0)
flownode=0
if (Bcoef(K).gt.1e19) then !danger recognize by force term
I=Bound(K)
do J=1,Sup(I,0)
flownode=flownode+Coef(I,J)*P(Sup(I,J))
enddo
flownode=flownode+Source(I)
endif
if (flownode.lt.0) Bcoef(K)=0.0
if (Bcoef(K).gt.1e19) flowout=flowout+flownode
enddo
return
end
c-----------------------------------------------------
Subroutine AdjustF(Nodes,F,P,Pstar,Cp,props,mprops,NodeA,maxF)
Integer props(1599,8)
Integer Nodes
Real NodeA(1599,0:8)
Real Cp(1599)
Real P(1599)
real Pstar(1599)
Real F(1599)
Real mprops(5,10)
Real maxF
Real DelF
Integer I,J,Kit
real Pnode(8)
real Pmax,Pmin,Pmid
real FstarM,FstarL
real ssearch
maxF=0.0
Do I = 1,Nodes
DelF= 0.25* Cp(I) * (Pstar(I) - P(I))
if(abs(DelF)/F(i).lt.0.025) then
F(i)=F(i)+delF
else
f(i)=f(i)*(1+(delF/abs(delF))*0.025)
endif
if(abs(Cp(i)*(Pstar(i)-P(i))).gt.maxF)
& maxF=abs(cp(i)*(Pstar(i)-P(i)))
Pmax = -1000
Pmin = 1000
Do J = 1,NodeA(I, 0)
Pnode(J) = Finv(F(I), props(I, J),mprops)
If(Pnode(J).ge.Pmax) Pmax = Pnode(J)
If(Pnode(J).le.Pmin) Pmin = Pnode(J)
enddo
If(Pmax.eq.Pmin) Then
Pmid = Pmin
Else
FstarL = 0
Do J = 1,NodeA(I, 0)
FstarL = FstarL + FNM(Pmin, props(I, J),mprops) * NodeA(I, J)
enddo
ssearch = (Pmax - Pmin)/2.
Do Kit = 1,20
Pmid = Pmin + ssearch
FstarM = 0
Do J = 1,NodeA(I, 0)
FstarM = FstarM + FNM(Pmid, props(I, J),mprops) * NodeA(I, J)
enddo
If((FstarM-F(I))*(FstarL-F(I)).le.0) Then
ssearch = -ssearch/2.
FstarL = FstarM
endif
Pmin = Pmid
enddo
endif
P(I) = P(I) + (Pmid - P(I))
enddo
return
end
c---------------------------------------
Subroutine Clean(Nodes,Sup,Source,Cp,Coef,M_D)
Integer Sup(1599,0:8)
Integer Nodes
Real Coef(1599,0:8)
Real Cp(1599)
Real Source(1599)
Real M_D(1599)
Integer I,J
Do I = 1,Nodes
Do J = 0,Sup(I, 0)
Coef(I, J) = 0
enddo
Source(I) = 0
Cp(I) = 0
M_D(I)=0
enddo
return
end
c--------------------------------
Subroutine Swap(Nodes,F,Fold)
Integer I,Nodes
real F(1599),Fold(1599)
Do I = 1,Nodes
Fold(I) = F(I)
enddo
return
end
c--------------------------------------------------------
Subroutine Trench(Nodes,X,Y,P,mprops)
Integer I
Integer Nodes
Real X(1599)
Real Y(1599)
Real P(1599)
Real mprops(5,10)
open(5,file='trench.dat')
do I=1,Nodes
if (X(I).ge.150.and.X(I).le.200) then
if (Y(I).ge.-15.and.Y(I).le.15) then
write(5,*) X(I),Y(I),FNM(P(I),2,mprops)
endif
endif
enddo
close(5)
return
end
c--------------------------
Real Function FNK(P1,P2,prop,mprops)
integer prop
real mprops(5,10)
real P1,P2
real fsat,fres,ksat,lamda
real F1,F2
F1 = FNM(P1, prop,mprops)
F2 = FNM(P2, prop,mprops)
real K1,K2
lamda = mprops(prop,6)
ksat = mprops(prop,5) ! 0.000098
fsat = mprops(prop,2)
fres = mprops(prop,3)
K1 = (F1 - fres) / (fsat - fres)
If(K1.gt.1) K1 = 1
K1 = ksat * K1**(lamda)
K2 = (F2 - fres) / (fsat - fres)
If(K2.gt.1) K2 = 1
K2 = ksat * K2**(lamda)
FNK = (K1 + K2) / 2.0
return
end
c -----------------------
real Function FNC(P, prop,mprops)
real P
Integer prop
real Cval
real satval
real hd
real lamda
real fres
real fsat
real F
real mprops(5,10)
F = FNM(P,prop,mprops)
lamda = mprops(prop,1)
fres = mprops(prop,3)
fsat = mprops(prop,2)
hd = mprops(prop,4)
satval = (F - fres) / (fsat - fres)
Cval = (fsat - fres) * (-lamda / hd) * satval**(1 + 1 / lamda)
If(satval.gt.1) Cval = 0.0001
FNC = Cval
return
end
c--------------------------
real Function FNM(P,prop,mprops)
integer prop
real P
real fval
real hd
real lamda
real fres
real fsat
real fhd
real mprops(5,10)
lamda = mprops(prop,1)
fres = mprops(prop,3)
fsat = mprops(prop,2)
hd = mprops(prop,4) !
fhd = mprops(prop,2)-mprops(prop,4)*0.0001
If (P.lt.hd) Then
fval = (fsat - fres) * ((P / hd)**(-lamda)) + fres
Else
fval = 0.0001 * P + fhd
endif
FNM = fval
return
end
c----------------------------------
real Function Finv(F, prop,mprops)
integer prop
real F
real pval
real hd
real lamda
real fres
real fsat
real fhd
real mprops(5,10)
lamda = mprops(prop,1)
fres = mprops(prop,3)
fsat = mprops(prop,2)
hd = mprops(prop,4) !
fhd = mprops(prop,2)-mprops(prop,4)*0.0001
If (F.lt.fsat) Then
pval = hd * ((F - fres) / (fsat - fres))**(-1 / lamda)
Else
pval = 10000*(F - fhd)
End If
Finv=pval
return
end
c--------------------------------