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