brettm
/
3DRWALK
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
3DRWALK/source/3dwalking.for

773 lines
25 KiB
Plaintext
Raw Permalink Normal View History Unescape Escape

First Commit
5 years ago
c---------------------------------------------------------------walkem
subroutine walkem(cord,nop,itype,ebox,vel1,vel2,vel,dudt,ao,
* nxyz,egrp,egps,nblk,nelb,nbb,nee,poll1,poll2,
* poll,PDOT,massp,pollwin,partE,partN,partZ,
* partI,partM,partA,partV,opbd,no_open,NUMPARTS,
* t,t_p,v_time,diffs)
c---------------------------------------------------------------------c
c purpose: c
c To perform random walking. c
c---------------------------------------------------------------------c
include '3duparms.cb'
include '3dpolls.cb'
include '3dgeom.cb'
c poll(1,i) = Vol of discharge at window i
c poll(2,i) = Conc of discharge at window i
c poll(3,i) = Depth of discharge at window i
c poll(4,i) = Radius of discharge at window i
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
common /etype/inode(4)
common /randm/isd, isd1
LOGICAL dead, walkjustone
INTEGER isd, isd1, p, pr, m, w, i, j, key, it, nen, num
INTEGER oldparts, newparts, deadparts, NUMPARTS, nxyz, pil
REAL r, ran3, tZ, tR, q, c, pel, pnl, pzl, pa, pv
REAL pe,pn,pz,xi,eta,zeta,shap
INTEGER t,t_p,v_time,no_open
INTEGER*2 nop(20,1)
INTEGER partI(1),nblk(2,1),nbb(1),nee(1),nelb(1),itype(1)
REAL partE(1),partN(1),partZ(1),partM(1),partA(1),vel(3,1)
REAL cord(3,1),egrp(6,1),egps(6,1),poll(4,1),pollwin(4,1)
REAL ebox(6,1),vel1(3,1),vel2(3,1),poll1(4,1),poll2(4,1)
REAL partV(1),opbd(4,1),dudt(3,1),PDOT(1),massp(1)
REAL TOFDAY1,DAYOFY1
REAL*8 ao(1)
REAL diffs(2,1)
oldparts=NUMPARTS
newparts=0
deadparts=0
num_steps=outputTS/dt
if(.not.hsteadystate.or.(.not.psteadystate)) then
num_steps=min(num_steps,newcondTS/dt)
endif
ibd=0
do i=1,num_steps
c
c----------calculate dieoff
c
if(dieoffc) then
TOFDAY=TOFDAY+dt/3600.0
if(TOFDAY.ge.24.0) then
TOFDAY=TOFDAY-24.0
DAYOFY=DAYOFY+1.0
if(DAYOFY.gt.366.0) DAYOFY=DAYOFY-366.0
endif
DELT=float(dt)
TOFDAY1=TOFDAY
DAYOFY1=DAYOFY
call DIEOFFS(TOFDAY1,DAYOFY1,T90_S)
cdrc
call DIEOFFS(0.0,DAYOFY1,T90_D)
c print *,'T90_S=',T90_S
crdc
c print *,'T90_D=',T90_D
T90_S=3600.0*T90_S
crdc
T90_D=3600.0*T90_D
endif
c
c----------Temporal interpolation of velocities
c
if(.not.hsteadystate.or.(.not.psteadystate)) then
call interp_vels(vel1,vel2,vel,dudt,poll1,
* poll2,poll,v_time)
endif
c
c----------Old particles
c
m=1
do while (m.le.NUMPARTS)
c
c---------------The particle is "active" and in the system.
c
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dt,diffs,ibd)
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partM(m)=pm
IF (dead.and.(.not.tracking)) THEN
c
c-----------------------particle is dead or out of system
c
NUMPARTS=NUMPARTS-1
do j=m,NUMPARTS
partI(j)=partI(j+1)
partE(j)=partE(j+1)
partN(j)=partN(j+1)
partZ(j)=partZ(j+1)
partM(j)=partM(j+1)
partA(j)=partA(j+1)
partV(j)=partV(j+1)
enddo
deadparts=deadparts+1
m=m-1
END IF
m=m+1
enddo
c
c----------New Particles.
c
NUMPARTS0=NUMPARTS
do ns=1,int(dt/dtnew)
c
c------------Temporal interpolation of velocities
c
if(ns.ne.1) then
if(.not.hsteadystate.or.(.not.psteadystate)) then
call interp_vels(vel1,vel2,vel,dudt,poll1,
* poll2,poll,v_time)
endif
endif
ccc num_s=num_steps-ns+1
c
c------------Old particles related to time increment dtnew
c
m=1+NUMPARTS0
do while (m.le.NUMPARTS)
c
c---------------The particle is "active" and in the system.
c
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd)
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partM(m)=pm
IF (dead.and.(.not.tracking)) THEN
c
c-----------------------particle is dead or out of system
c
NUMPARTS=NUMPARTS-1
do j=m,NUMPARTS
partI(j)=partI(j+1)
partE(j)=partE(j+1)
partN(j)=partN(j+1)
partZ(j)=partZ(j+1)
partM(j)=partM(j+1)
partA(j)=partA(j+1)
partV(j)=partV(j+1)
enddo
deadparts=deadparts+1
m=m-1
END IF
m=m+1
enddo
do w=1,num_pollwins
c
c-------------Establish the number of particles to be released.
c
q=poll(1,w)
c=poll(2,w)
PDOT(w)=PDOT(w)+(c*q*dtnew)/massp(w)
P=INT(PDOT(w))
PDOT(w)=PDOT(w)-P
c p=100000
c if(i.gt.1.or.ns.gt.1) p=0
do pr=1,P
NUMPARTS=NUMPARTS+1
m=NUMPARTS
r=ran3(isd)
pe=pollwin(1,w)+r*(pollwin(3,w)-pollwin(1,w))
pn=pollwin(2,w)+r*(pollwin(4,w)-pollwin(2,w))
partA(m)=0.
pm=massp(w)
tZ=-poll(3,w)
tR=poll(4,w)
r=2.0*(ran3(isd)-0.5)
pz=tZ+r*tR
pz_keep=pz
c
c---------------Check if the particle is in the water column.
c
IF (pz.gt.0.0) pz=0.0
c i.e. if it is above the water surface then put it on the surface
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,
* nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key)
num_keep1=num
if(key.eq.1) then
c i.e. the pe,pn,pz coordinate was inside of element num
partE(m)=xi
partN(m)=eta
partZ(m)=zeta
partI(m)=num
else
c try and see if the plan form (pe,pn) is inside an element
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,
* nelb,nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta,
* zeta,num,key)
num_keep2=num
if(key.eq.1) then
c the (pe, pn) point was in element num, but we knew that (pe,pn,pz) wasn't so the
c old pz value must have been below the bed.
pz=0.0
call xn3(it,nen,xi,eta,zeta)
do j=1,nen
pz=pz+shap(j)*ao(nop(j,num))
enddo
pz=0.99*pz
c set the pz value to be at 99% of the depth and all is well.
else
c the (pe,pn) pair was not found to be in an element.
c This may happen if the (pe, pn) pair falls directly onto a element boundary.
write(*,'(a)') ' *** Window for discharging '//
* 'plumes out of the modelling domain (1)***'
write(*,'(i5,3f15.5)') 0,pe,pn,pz
stop
endif
num=-1
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,
* nelb,nbb,nee,nxyz,nen,it,pe,pn,pz,xi,eta,
* zeta,num,key)
if(key.eq.1) then
partE(m)=xi
partN(m)=eta
partZ(m)=zeta
partI(m)=num
else
write(*,'(a)') ' *** Window for discharging '//
* 'plumes out of the modelling domain (2)***'
write(*,'(i5,3f15.5)') 1,pe,pn,pz
write(*,*) 'orig pz, num_keep1, num_keep2',pz_keep
* ,num_keep1, num_keep2
stop
endif
endif
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
call rise_fall_vel(pv)
c following line for plume only
if(plumeonly) pv=0.0
c following few lines for floatables only
if(floatableonly) then
if(pv.le.0.0) then
NUMPARTS = NUMPARTS - 1
goto 70
endif
endif
c------------------------------------------
if(settles.and.pv.ge.0.0) then
NUMPARTS = NUMPARTS - 1
goto 70
endif
dead=walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,
* ebox,nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,
* nee,nxyz,pe,pn,pz,no_open,opbd,dtnew,diffs,ibd)
newparts=newparts+1
partE(m)=pel
partN(m)=pnl
partZ(m)=pzl
partI(m)=pil
partA(m)=pa
partV(m)=pv
partm(m)=pm
IF (dead.and.(.not.tracking)) THEN
NUMPARTS=NUMPARTS-1
deadparts=deadparts+1
END IF
70 enddo
enddo
v_time=v_time+dtnew
enddo
c---------------------
c
c----------output tracking path
c
if(tracking.and.t+dt.eq.outputTS) then
ctmp inquire(35,NEXTREC=ntrack)
ctmp write(35,rec=ntrack) NUMPARTS,t_p
write(35,*) NUMPARTS,t_p
ntrack=ntrack+1
do m=1,NUMPARTS
pel=partE(m)
pnl=partN(m)
pzl=partZ(m)
pil=partI(m)
pa=partA(m)
pv=partV(m)
pm=partM(m)
it=itype(pil)
nen=inode(it)
call xn3(it,nen,pel,pnl,pzl)
pe=0.0
pn=0.0
pz=0.0
do j=1,nen
pe=pe+shap(j)*cord(1,nop(j,pil))
pn=pn+shap(j)*cord(2,nop(j,pil))
pz=pz+shap(j)*cord(3,nop(j,pil))
enddo
call calc_ab(aE,bN,pe,pn)
ctmp write(35,rec=ntrack) aE*resXY,bN*resXY
write(35,*) pe,pn,pz,pm
ntrack=ntrack+1
enddo
endif
t=t+dt
t_p=t_p+dt
enddo
print *,'ibd=',ibd
if(.not.hsteadystate.or.(.not.psteadystate)) then
WRITE(*,200) min(outputTS,newcondTS), oldparts,
* newparts,deadparts
else
WRITE(*,200) outputTS, oldparts,
* newparts,deadparts
endif
200 FORMAT('During the last',I5,' seconds:- [ Old:',I8,
1 ' ; New:',I7,' ; Dead:',I7,' ]')
RETURN
END
c-----------------------------------------------------------walkjustone
FUNCTION walkjustone(pel,pnl,pzl,pil,pa,pv,pm,cord,itype,ebox,
* nop,vel,dudt,ao,egrp,egps,nblk,nelb,nbb,nee,nxyz,pe,
* pn,pz,no_open,opbd,dtt,diffs,ibd)
c---------------------------------------------------------------------c
c purpose: c
c To walk one step. c
c---------------------------------------------------------------------c
include '3duparms.cb'
common /shape/shap(20),shpx(20),shpy(20),shpz(20)
common /etype/inode(4)
common /randm/isd, isd1
LOGICAL walkjustone
REAL cord(3,1),vel(3,1),egrp(6,1),egps(6,1),ebox(6,1)
REAL opbd(4,1),dudt(3,1)
REAL pe,pn,pz,pa,pel,pnl,pzl,xi,eta,zeta
REAL peold,pnold,pzold,peold1,pnold1,pzold1
REAL*8 ao(1)
INTEGER*2 nop(20,1)
INTEGER nblk(2,1),nelb(1),nbb(1),nee(1),itype(1)
INTEGER num,pil,key,nen,nxyz,it,no_open
REAL ran3, normdis
LOGICAL boundaries
INTEGER isd, isd1
REAL len,ang
REAL vE,vN,vZ
REAL pv
real*8 a(3,3),b(3),ai(3,3),du(3,3),a1,a2,a3,dudt0,dvdt0,dwdt0
real*8 detj,ass
integer dtt
real diffs(2,1)
walkjustone=.FALSE.
ccc do i=1,num_steps
pa=pa+dtt
if(pa.gt.tot_dead_age.or.pm/mass_pp.le.0.001) then
c if(pa.gt.tot_dead_age) then
walkjustone=.TRUE.
return
endif
c
c----------gradually die-off
c
if(dieoff.and.pv.eq.0.0) then
c print*,tdpth
if(abs(pz).le.tdpth) then
c print*,'surface dieoff=',10.0**(-dtt/t90_s)
c print*,pm
pm=pm*10.0**(-dtt/t90_s)
c print*,pm
else
c print*,'deep dieoff=',10.0**(-dtt/t90_d)
c print*,pm
pm=pm*10.0**(-dtt/t90_d)
c print*,pm
endif
endif
c
c----------settleable particles on bed -- return
c
if(pv.eq.-9999.0) return
c
c----------interpolate velocity from FE mesh
c
vE=0.0
vN=0.0
vZ=0.0
it=itype(pil)
nen=inode(it)
ccc call xn3(it,nen,pel,pnl,pzl)
do j=1,nen
vE=vE+shap(j)*vel(1,nop(j,pil))
vN=vN+shap(j)*vel(2,nop(j,pil))
vZ=vZ+shap(j)*vel(3,nop(j,pil))
c print *,'pil,vel=',pil,j,vel(3,nop(j,pil))
enddo
c print *,'vE,vN,vZ1=',vE,vN,vZ
c
c----------predict the trajectory of the particle
c
if(.not.euler) then
call xn3x(it,nen,pel,pnl,pzl)
c
c----------calculate the Jacobian matrix
c
do j=1,3
a(j,1)=0.0
a(j,2)=0.0
a(j,3)=0.0
du(1,j)=0.0
du(2,j)=0.0
du(3,j)=0.0
do k=1,nen
nele=nop(k,pil)
a(j,1)=a(j,1)+shpx(k)*cord(j,nele)
a(j,2)=a(j,2)+shpy(k)*cord(j,nele)
a(j,3)=a(j,3)+shpz(k)*cord(j,nele)
du(1,j)=du(1,j)+shpx(k)*vel(j,nele)
du(2,j)=du(2,j)+shpy(k)*vel(j,nele)
du(3,j)=du(3,j)+shpz(k)*vel(j,nele)
enddo
enddo
a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
c if(detj.lt.1.0e-10) then
c print *,'detj,pil,nen,pel,pnl,pzl=',detj,pil,nen,pel,pnl,pzl
c stop
c endif
c
c----------calculate the inverse Jacobi matrix
c
ai(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2))/detj
ai(1,2)=(a(3,1)*a(2,3)-a(2,1)*a(3,3))/detj
ai(1,3)=(a(2,1)*a(3,2)-a(3,1)*a(2,2))/detj
ai(2,1)=(a(3,2)*a(1,3)-a(1,2)*a(3,3))/detj
ai(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1))/detj
ai(2,3)=(a(3,1)*a(1,2)-a(1,1)*a(3,2))/detj
ai(3,1)=(a(1,2)*a(2,3)-a(2,2)*a(1,3))/detj
ai(3,2)=(a(2,1)*a(1,3)-a(1,1)*a(2,3))/detj
ai(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1))/detj
do m=1,3
do j=1,3
a(m,j)=0.0
do k=1,3
a(m,j)=a(m,j)+dtt*ai(j,k)*du(k,m)
enddo
c a(m,j)=dtt*a(m,j)
enddo
c print *,'a=',(a(m,j),j=1,3)
enddo
do k=1,3
a(k,k)=a(k,k)-2.0
enddo
dudt0=0.0
dvdt0=0.0
dwdt0=0.0
do k=1,nen
dudt0=dudt0+shap(k)*dudt(1,nop(k,pil))
dvdt0=dvdt0+shap(k)*dudt(2,nop(k,pil))
dwdt0=dwdt0+shap(k)*dudt(3,nop(k,pil))
enddo
b(1)=-2.0*vE-dudt0*dtt
b(2)=-2.0*vN-dvdt0*dtt
b(3)=-2.0*vZ-dwdt0*dtt
c do i=1,3
c ass=max(dabs(a(i,1)),dabs(a(i,2)),dabs(a(i,3)))
c do j=1,3
c a(i,j)=a(i,j)/ass
c enddo
c b(i)=b(i)/ass
c 25 enddo
a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
c detj=detj*dtt
vE=(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(3,1)*a(2,3)-a(3,3)*a(2,1)
a2=a(1,1)*a(3,3)-a(1,3)*a(3,1)
a3=a(2,1)*a(1,3)-a(2,3)*a(1,1)
vN=(a1*b(1)+a2*b(2)+a3*b(3))/detj
a1=a(2,1)*a(3,2)-a(2,2)*a(3,1)
a2=a(3,1)*a(1,2)-a(3,2)*a(1,1)
a3=a(1,1)*a(2,2)-a(1,2)*a(2,1)
vZ=(a1*b(1)+a2*b(2)+a3*b(3))/detj
c print *,'vE,vN,vZ2=',vE,vN,vZ,dtt
c stop
endif
c-------------------------------
peold=pe
pnold=pn
pzold=pz
len=sqrt(6*diffs(1,pil)*dtt)*ran3(isd)
ang=6.2831853*ran3(isd)
angz=6.2831853*ran3(isd)
c if(vdiff.ne.0.0) then
c angz=6.2831853*ran3(isd)
c else
c angz=0.0
c endif
peold1 = pe + vE*dtt
pnold1 = pn + vN*dtt
pe = peold1 + len*cos(ang)*cos(angz)
pn = pnold1 + len*sin(ang)*cos(angz)
if(pv.eq.0.0) then
pzold1 = pz + vZ*dtt
pz = pzold1 + ran3(isd)*sqrt(6*diffs(2,pil)*dtt)*sin(angz)
else
pzold1 = pz + (vZ+pv)*dtt
pz = pzold1
endif
IF(boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd)) THEN
walkjustone=.TRUE.
RETURN
ENDIF
c -->
IF (pz.gt.0.0) pz=0.0
num=pil
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb,
* nee,nxyz,nen,it,pe,pn,pz,xi,eta,zeta,num,key)
if(key.eq.1) then
pel=xi
pnl=eta
pzl=zeta
pil=num
else
c
c-------------neglect the random displacements
c
num=pil
if(pnold1.gt.0.0) pnold1=0.0
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,nbb,
* nee,nxyz,nen,it,peold1,pnold1,pzold1,xi,eta,zeta,num,key)
if(key.eq.1) then
pe=peold1
pn=pnold1
pz=pzold1
pel=xi
pnl=eta
pzl=zeta
pil=num
c
c-------------neutral bouyant and floatable particales
c
else if(pv.ge.0.0) then
pe=peold
pn=pnold
pz=pzold
ibd=ibd+1
c
c-------------settleable particles
c
else
num=pil
call trinvs3(nop,cord,itype,ebox,egrp,nblk,egps,nelb,
* nbb,nee,nxyz,nen,it,pe,pn,0.0,xi,eta,zeta,num,key)
if(key.eq.1) then
if(settles) then
no_settles=no_settles+1
write(30,rec=no_settles) pe,pn
walkjustone=.true.
else
pz=0.0
call xn3(it,nen,xi,eta,zeta)
do j=1,nen
pz=pz+shap(j)*ao(nop(j,num))
enddo
pv=-9999.0
endif
else
pe=peold
pn=pnold
pz=pzold
endif
endif
endif
ccc enddo
RETURN
END
c------------------------------------------------------------boundaries
FUNCTION boundaries(pe,pn,pz,peold,pnold,pzold,no_open,opbd)
c---------------------------------------------------------------------c
c purpose: c
c To check whether the particle within the caculation region. c
c---------------------------------------------------------------------c
include '3duparms.cb'
include '3dgeom.cb'
REAL pe, pn, pz, peold, pnold, pzold
REAL a, b
INTEGER no_open
LOGICAL boundaries
real opbd(4,1),s,r
real*8 deta
c This is very rudimentary. It must be modified later.
c A particle can leave the domain in two ways:-
c i) Leave the (x,y) range.
c ii) Enter a part of the (x,y) range that is "land".
c Any particle leaving by mode (i) are considered to be "gone"
c Particles leaving by (ii) are returned to their old position.
c
c To particle is considered to entered a land area if:-
c a) It is beneath the bathymetry.
c b) It is in a grid cell that doesn't have all four
c corners within the ocean area.
c
c The function returns true only if the particle is "gone".
boundaries=.FALSE.
c CALL calc_ab(a, b, pe, pn)
c Check if particle has left the array region. Case i
c IF ((a.le.0.0).or.(a.ge.float(numX)).or.
c 1 (b.le.0.0).or.(b.ge.float(numY))) THEN
c boundaries=.TRUE.
c RETURN
c ENDIF
c-------check if particle has left from specified open boundaries
if(openbd) then
do i=1,no_open
deta=(pnold-pn)*(opbd(3,i)-opbd(1,i))
* -(peold-pe)*(opbd(4,i)-opbd(2,i))
if(dabs(deta).lt.1.0e-8) goto 10
r=((pnold-opbd(2,i))*(opbd(3,i)-opbd(1,i))
* -(peold-opbd(1,i))*(opbd(4,i)-opbd(2,i)))/deta
if(r.lt.-1.0e-5.or.r.gt.1.0+1.0e-5) goto 10
s=((pnold-pn)*(peold-opbd(1,i))
* -(peold-pe)*(pnold-opbd(2,i)))/deta
if(s.le.1.0+1.0e-5.and.s.ge.-1.0e-5) then
boundaries=.TRUE.
return
endif
10 enddo
endif
c Check for particles below the bathymetry. (Case (iia))
RETURN
END
c-----------------------------------------------------------interp_vels
subroutine interp_vels(vel1,vel2,vel,dudt,poll1,poll2,poll,
* v_time)
c---------------------------------------------------------------------c
c purpose: c
c To perform temporal interpolation of velocities. c
c---------------------------------------------------------------------c
INCLUDE '3duparms.cb'
INCLUDE '3dpolls.cb'
INCLUDE '3dmesh.cb'
dimension vel1(3,1),vel2(3,1),vel(3,1),poll1(4,1),poll2(4,1),
* poll(4,1),dudt(3,1)
REAL r
INTEGER i,j,v_time
c
r=float(v_time)/float(newcondTS)
do i=1,np
do j=1,3
vel(j,i)=vel1(j,i)+r*(vel2(j,i)-vel1(j,i))
enddo
enddo
do i=1,num_pollwins
do j=1,4
poll(j,i)=poll1(j,i)+r*(poll2(j,i)-poll1(j,i))
enddo
enddo
c
c-------calculate dudt
c
if(.not.euler) then
do i=1,np
do j=1,3
dudt(j,i)=(vel2(j,i)-vel1(j,i))/float(newcondTS)
enddo
enddo
endif
return
end