Civil Engineering Reference
In-Depth Information
CALL seep4(kp,coord,kay); CALL fsparv(kv,kp,num,kdiag)
END DO elements_2; kvh=kv
!-----------------------specify boundary values---------------------------
loads=zero; kv(kdiag(node_up))=kv(kdiag(node_up))+penalty
loads(node_up)=kv(kdiag(node_up))*hup
kv(kdiag(node_down))=kv(kdiag(node_down))+penalty
loads(node_down)=kv(kdiag(node_down))*hdown
kv(kdiag(node_seep))=kv(kdiag(node_seep))+penalty
DO i=1,fixed_seep
loads(node_seep(i))=kv(kdiag(node_seep(i)))*
&
(hdown+(surf(1)-hdown)*(fixed_seep+1-i)/(fixed_seep+1))
END DO
!-----------------------equation solution---------------------------------
CALL sparin(kv,kdiag); CALL spabac(kv,loads,kdiag)
surf(1:nxe)=loads(1:nxe)
!-----------------------check convergence---------------------------------
CALL checon(loads,oldpot,tol,converged)
IF(converged.OR.iters==limit)EXIT
END DO its; CALL linmul_sky(kvh,loads,disps,kdiag)
WRITE(11,'(/A)')" Node Total Head Flow rate"
DO k=1,nn; WRITE(11,'(I5,2E12.4)')k,loads(k),disps(k); END DO
disps(0)=zero; WRITE(11,'(/A)')"
Inflow
Outflow"
WRITE(11,'(5X,2E12.4)')
&
SUM(disps,MASK=disps<zero),SUM(disps,MASK=disps>zero)
WRITE(11,'(/A,I3,A)')" Converged in",iters," iterations"
CALL mesh(g_coord,g_num,12)
READ(10,*)nci; CALL contour(loads,g_coord,g_num,nci,13)
STOP
END PROGRAM p73
New scalar integers:
fixed down
number of nodes on downstream side
number of nodes on seepage surface
fixed seep
number of nodes on upstream side
fixed up
counts free-surface iterations
iters
iteration ceiling
limit
New scalar reals:
hdown
fixed total head on downstream side
hup
fixed total head on upstream side
d180
set to 180
.
0
initial height
initial height of free surface to start process
tol
convergence tolerance
Scalar logical:
converged
set to
.TRUE.
if mesh has converged
New dynamic integer arrays:
node down
nodes fixed on downstream side
nodes fixed on downstream seepage surface
node seep
nodes fixed on upstream side
node up
