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.
773 lines
25 KiB
Fortran
773 lines
25 KiB
Fortran
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
|
|
|