subroutine swrun(uw,qw,vw,dw,u1,q1,v1,n,h & ,ux1,vx1,qx1,ux2,vx2,qx2,t) real*8 u1(n),q1(n),v1(n),vw(n),uw(n),h(n),qw(n),dw(n) real*8 u1p(3000),q1p(3000) cA The above array should be at least the same size as the u1 array real*8 oldhl(500),oldl(500),oldhr(500),oldr(500) real*8 dwl(500),dwr(500),qwl(500),qwr(500),uwl(500),uwr(500) real*8 lam1,lam2,qww,dww,uuu,d1,d2,t1,t,time,fi,fr real*8 ux1,vx1,qx1,ux2,vx2,qx2,PI,RE,GA,groff,grnd real*8 dwp1,qwp1,h0,h01,h02,hnew,hm1 integer nbeg(500),nend(500) logical wet1,wetn,wet2,wetn1 common /run/time,nstep,npics,groff,grnd,fi,fr common /const/ PI,RE,GA c.....fr is n**2, where n is Manning coefficient CC no narr nests cA the code has been modified to keep old characteristic information in cA different variables, corrected errors in bottom friction and modified cA velocities at the wet/dry interface (Arun - 1/30/06) do i=1,500 nbeg(i)=0 nend(i)=0 end do CC check to see if for this row the start (1) is wet if (qw(1).GT.grnd)then wet1=.TRUE. else wet1=.FALSE. end if CC check to see if for this row the end (n) is wet if (qw(n).GT.grnd)then wetn=.TRUE. else wetn=.FALSE. end if CC initialize the counter ks (segment counter) qww dww and fix the last metric CC to the next to last this does change h(n) ks=1 qww=qw(1) dww=dw(1) h(n)=h(n-1) c***** Checking shoreline ******************** do 111 i=1,n CC set the i plus one value for look ahead if (i.NE.n) then dwp1=dw(i+1) qwp1=qw(i+1) else qwp1=qw(n) dwp1=dw(n) end if if(qw(i).GT.grnd)then if(nbeg(ks).EQ.0) nbeg(ks)=i qww=qw(i) dww=dw(i) u1(i)=uw(i)+2.*sqrt(GA*qw(i)) q1(i)=uw(i)-2.*sqrt(GA*qw(i)) q1p(i)=q1(i) u1p(i)=u1(i) v1(i)=vw(i) else if((qww.LE.grnd).AND. 1 (qwp1.LE.grnd))then c******* Ground ****************************** c qw(i)=0. c uw(i)=0. c vw(i)=0. qww=qw(i) dww=dw(i) goto111 else if(qww.LE.grnd)then c******* Shore on the left ******************* nbeg(ks)=i oldhl(ks)=h(i) oldl(ks)=dw(i) if(i.NE.1) then hm1=h(i-1) else hm1=h(1) end if if( dwp1-dw(i) .GE. qwp1 )then hnew=oldhl(ks)*qwp1/(dwp1-dw(i)) dw(nbeg(ks))=dwp1-qwp1 h(nbeg(ks))=hnew h01=t*(abs(uw(i+1))+sqrt(GA*qwp1)) h02=t*GA*abs(dwp1-dw(i))/(sqrt(GA*qwp1)) C h0 = max(h01,h02) h0 = h01 if (h0.GT.hnew) then qw(i+1)=0. uw(i+1)=0. vw(i+1)=0. dw(nbeg(ks))=oldl(ks) h(nbeg(ks))=oldhl(ks) nbeg(ks)=0 end if else if(dwp1-dww .GT. qwp1) then hnew=oldhl(ks)+ 1 hm1*(qwp1-(dwp1-dw(i)))/(dw(i)-dww) dw(i)=dwp1-qwp1 h(i)=hnew c ... Case three is handled below c else c hnew=oldhl(ks)+h(i) c dw(i)=dww c h(i)=hnew c$$$ else ! fix ET c$$$ dw(i)=0 end if qww=0. dww=oldl(ks) if(uw(i+1).LT.0.)then uw(i)=uw(i+1) vw(i)=vw(i+1) q1p(i)=uw(i) u1p(i)=uw(i) else uw(i)=0. q1p(i)=0. u1p(i)=0. vw(i)=0. end if qw(i)=0. q1(i)=q1p(i) u1(i)=u1p(i) v1(i)=vw(i) else if(qwp1.LE.grnd) then c******* Shore on the right ****************** nend(ks)=i oldhr(ks)=h(i-1) oldr(ks)=dw(i) if(dww-dw(i) .GE. qww)then hnew=oldhr(ks)*qww/(dww-dw(i)) h01=t*(abs(uw(i-1))+sqrt(GA*qww)) h02=t*GA*abs(dww-dw(i))/(sqrt(GA*qww)) C h0 = max(h01,h02) h0 = h01 if (hnew.LT.h0) then nend(ks)=i-1 hnew=hnew+h(i-2) oldr(ks)=dw(i-1) if(nbeg(ks).NE.i-2)then oldhr(ks)=h(i-2) else oldhr(ks)=oldhl(ks) end if end if dw(nend(ks))=dww-qww h(nend(ks)-1)=hnew else if(dww-dwp1 .GT. qww) then hnew=oldhr(ks)+ 1 h(i)*(qww-(dww-dw(i)))/(dw(i)-dwp1) dw(nend(ks))=dww-qww h(nend(ks)-1)=hnew c ... Case three is handled below c else c hnew=oldhr(ks)+h(i) c dw(nend(ks))=dwp1 c h(nend(ks)-1)=hnew c$$$ else ! ET c$$$ dw(i)=0 end if qww=0. dww=oldr(ks) if(uw(nend(ks)-1).GT.0.)then uw(nend(ks))=uw(nend(ks)-1) vw(nend(ks))=vw(nend(ks)-1) q1p(nend(ks))=uw(nend(ks)) u1p(nend(ks))=uw(nend(ks)) else uw(nend(ks))=0. vw(nend(ks))=0. q1p(nend(ks))=0. u1p(nend(ks))=0. end if C set eta at the shore to minus dw qw(nend(ks))=0. v1(nend(ks))=vw(nend(ks)) q1(nend(ks))=q1p(nend(ks)) u1(nend(ks))=u1p(nend(ks)) ks=ks+1 else c******* Isolated dry point ********************** nend(ks)=i oldr(ks)=dw(i) oldhr(ks)=h(i-1) dwr(ks)=dw(i-1) qwr(ks)=-2.*sqrt(GA*qww) uwr(ks)=2.*sqrt(GA*qww) nbeg(ks+1)=i oldl(ks+1)=dw(i) oldhl(ks+1)=h(i) dwl(ks+1)=dw(i+1) qwl(ks+1)=-2.*sqrt(GA*qw(i+1)) uwl(ks+1)=2.*sqrt(GA*qw(i+1)) vw(i)=0. qww=qw(i) dww=dw(i) ks=ks+1 end if 111 continue if(nbeg(ks).NE.0) then nend(ks)=n else ks=ks-1 end if if (qw(2).LE.grnd) then wet2 = .false. else wet2 = .true. endif if (qw(n-1).LE.grnd) then wetn1 = .false. else wetn1 = .true. endif c*********Sea Boundary ***************************** if((dw(1)+qx1).GT.grnd.AND.wet1.AND.wet2)then c ****** deep water on left boundary q1(1)=q1p(1)-(u1p(1)+u1p(2)+3.*(q1p(1)+q1p(2)))/8.* & (q1p(2)-q1p(1))*t/h(1)+GA*t/h(1)*(dw(2)-dw(1)) u1(1)=ux1+2.*sqrt(GA*(dw(1)+qx1)) if(uw(1).LT.0.) then v1(1)=vw(1)-t/h(1)*(vw(2)-vw(1))*uw(1) else v1(1)=vx1 end if else if (wet1.AND.wet2) then c ****** wet shore on left boundary q1(1)=q1p(1)-(u1p(1)+u1p(2)+3.*(q1p(1)+q1p(2)))/8.* & (q1p(2)-q1p(1))*t/h(1)+GA*t/h(1)*(dw(2)-dw(1)) u1(1)=-q1(1) if(uw(1).GT.0.) then v1(1)=vw(1)-t/h(1)*(vw(2)-vw(1))*uw(1) else v1(1)=0. end if else if ((dw(1)+qx1).GT.grnd) then u1(1)=ux1+2.*sqrt(GA*(dw(1)+qx1)) q1(1)=ux1-2.*sqrt(GA*(dw(1)+qx1)) v1(1)=vx1 qw(1)=dw(1)+qx1 end if if((dw(n)+qx2).GT.grnd.AND.wetn.AND.wetn1)then c ****** deep water on right if(uw(n).GT.0) then v1(n)=vw(n)-t/h(n-1)*(vw(n)-vw(n-1))*uw(n) else v1(n)=vx2 end if q1(n)=ux2-2.*sqrt(GA*(dw(n)+qx2)) u1(n)=u1p(n)-(3.*(u1p(n)+u1p(n-1))+q1p(n-1)+q1p(n))/8.* & (u1p(n)-u1p(n-1))*t/h(n-1)+GA*t/h(n-1)*(dw(n)-dw(n-1)) else if (wetn.AND.wetn1) then c ****** wet shore on right if(uw(n).GT.0) then v1(n)=vw(n)-t/h(n-1)*(vw(n)-vw(n-1))*uw(n) else v1(n)=0. end if u1(n)=u1p(n)-(3.*(u1p(n)+u1p(n-1))+q1p(n-1)+q1p(n))/8.* & (u1p(n)-u1p(n-1))*t/h(n-1)+GA*t/h(n-1)*(dw(n)-dw(n-1)) q1(n)=-u1(n) else if ((dw(n)+qx2).GT.grnd) then u1(n)=ux2+2.*sqrt(GA*(dw(n)+qx2)) q1(n)=ux2-2.*sqrt(GA*(dw(n)+qx2)) v1(n)=vx2 qw(n)=dw(n)+qx2 end if c******* Space integration in 'wet' points ********** do 100 j=1,ks CC write(*,*) ks,j,nbeg(j),nend(j) CC Aurelio's catch if (j.gt.1) then if (nbeg(j).EQ.nend(j-1)) then dw(nbeg(j))=dwl(j) q1p(nbeg(j))=qwl(j) u1p(nbeg(j))=uwl(j) end if end if if (j.lt.ks) then if (nend(j).EQ.nbeg(j+1)) then dw(nend(j))=dwr(j) q1p(nend(j))=qwr(j) u1p(nend(j))=uwr(j) end if end if do 10 i=nbeg(j)+1,nend(j)-1 C this check enhances the shock wave generation at low C delta t so skipping it C if ((nbeg(j)+1).GE.(nend(j)-1)) then C u1(i) = u1p(i) C q1(i) = q1p(i) C v1(i) = vw(i) C else lam1=3.*u1p(i)+q1p(i) lam2=u1p(i)+3.*q1p(i) d1=GA*t*((dw(i+1)-dw(i))/h(i)-(dw(i)-dw(i-1))/h(i-1)) d2=GA*(dw(i+1)-dw(i-1)) & -fr*uw(i)*(h(i-1)+h(i))*sqrt(uw(i)*uw(i)+vw(i)*vw(i)) & *GA/(qw(i))**1.3333333333333D0 t1=t/(h(i-1)+h(i)) u1(i)=u1p(i)-t1*((3.*u1p(i-1)+q1p(i-1)+3.*u1p(i+1) & +q1p(i+1))/8.*(u1p(i+1)-u1p(i-1))-d2+d1* & lam1/4.-t/32.*lam1* & ((3.*u1p(i+1)+q1p(i+1)+lam1)*(u1p(i+1)-u1p(i))/h(i)- & (3.*u1p(i-1)+q1p(i-1)+lam1)*(u1p(i)-u1p(i-1))/h(i-1))) q1(i)=q1p(i)-t1*((u1p(i-1)+3.*q1p(i-1)+u1p(i+1)+ & 3.*q1p(i+1))/8.*(q1p(i+1)-q1p(i-1))-d2+d1* & lam2/4.-t/32.*lam2* & ((u1p(i+1)+3.*q1p(i+1)+lam2)*(q1p(i+1)-q1p(i))/h(i)- & (lam2+u1p(i-1)+3.*q1p(i-1))*(q1p(i)-q1p(i-1))/h(i-1))) v1(i)=vw(i)-t1*(uw(i)*(vw(i+1)- & vw(i-1))-uw(i)/2.*t*((uw(i+1)+uw(i))* & (vw(i+1)-vw(i))/h(i)-(uw(i)+uw(i-1))* & (vw(i)-vw(i-1))/h(i-1))) C endif 10 continue do 110 i=nbeg(j)+1,nend(j)-1 qw(i)=(u1(i)-q1(i))**2/(16.*GA) if (qw(i).GE.grnd) then uw(i)=(u1(i)+q1(i))*.5 vw(i)=v1(i) c froude number limit for all wet points c if (abs(uw(i)).GE.1.0*sqrt(GA*qw(i))) then c uw(i) = sign(1.0*sqrt(GA*qw(i)),uw(i)) c endif c if (abs(vw(i)).GE.1.0*sqrt(GA*qw(i))) then c vw(i) = sign(1.0*sqrt(GA*qw(i)),vw(i)) c endif else qw(i)=grnd uw(i)=0. vw(i)=0. end if 110 continue c***** Right shoreline ************************************** if((nend(j).EQ.n).AND.wetn) goto99 if(nend(j).NE.n) then if(dw(nend(j)-1)-oldr(j) .GE. qw(nend(j)-1))then qw(nend(j))=qw(nend(j)-1)-dw(nend(j)-1)+oldr(j) else if(dw(nend(j)-1)-dw(nend(j)+1).GE.qw(nend(j)-1))then qw(nend(j))=qw(nend(j)-1)-dw(nend(j)-1)+oldr(j) else qw(nend(j))=oldr(j)-dw(nend(j)+1)+ 1 h(nend(j))/(h(nend(j))+oldhr(j))* 2 (qw(nend(j)-1)-dw(nend(j)-1)+dw(nend(j)+1)) if (qw(nend(j)).LT.0.) 1 qw(nend(j))=qw(nend(j)-1)*uw(nend(j)-1)*t/ 2 (uw(nend(j)-1)*t+oldhr(j)) end if else if( dw(nend(j)-1)-oldr(j) .GE. qw(nend(j)-1)/2. )then qw(nend(j))=qw(nend(j)-1)-dw(nend(j)-1)+oldr(j) else qw(nend(j))=qw(nend(j)-1)*.5 end if end if c old method uw a chacteristic C if (uw(nend(j)-1).GT.0)then if ((uw(nend(j)-1).GT.0).AND.(qw(nend(j)).GT.grnd))then vw(nend(j))=vw(nend(j)-1) uw(nend(j))=uw(nend(j)-1) C froude near the beach if (abs(uw(nend(j))).GE.1.0*sqrt(GA*qw(nend(j)))) then uw(nend(j)) = sign(1.0*sqrt(GA*qw(nend(j))),uw(nend(j))) endif if (abs(vw(nend(j))).GE.1.0*sqrt(GA*qw(nend(j)))) then vw(nend(j)) = sign(1.0*sqrt(GA*qw(nend(j))),vw(nend(j))) endif else vw(nend(j))=0. uw(nend(j))=0. end if C if (qw(nend(j)).GT.grnd.AND.qw(nend(j)-1).GT.grnd)then c vw(nend(j))=vw(nend(j)-1)*qw(nend(j)-1)/qw(nend(j)) C vw(nend(j))=vw(nend(j)-1) C uw(nend(j))=uw(nend(j)-1)*qw(nend(j)-1)/qw(nend(j)) C if (abs(uw(nend(j))).GE.1.0*sqrt(GA*qw(nend(j)))) then C uw(nend(j)) = sign(1.0*sqrt(GA*qw(nend(j))),uw(nend(j))) C endif C if (abs(vw(nend(j))).GE.1.0*sqrt(GA*qw(nend(j)))) then C vw(nend(j)) = sign(1.0*sqrt(GA*qw(nend(j))),vw(nend(j))) C endif C else C vw(nend(j))=0. C uw(nend(j))=0. C end if dw(nend(j))=oldr(j) h(nend(j)-1)=oldhr(j) if (nend(j).EQ.nbeg(j+1)) qwr(j)=qw(nend(j)) c***** Left shoreline *************************************** 99 if((nbeg(j).EQ.1).AND.wet1) goto100 if (nbeg(j).NE.1) then if(dw(nbeg(j)+1)-oldl(j) .GE. qw(nbeg(j)+1))then qw(nbeg(j))=qw(nbeg(j)+1)-dw(nbeg(j)+1)+oldl(j) else if(dw(nbeg(j)+1)-dw(nbeg(j)-1).GE.qw(nbeg(j)+1))then qw(nbeg(j))=qw(nbeg(j)+1)-dw(nbeg(j)+1)+oldl(j) else qw(nbeg(j))=oldl(j)-dw(nbeg(j)-1)+ 1 h(nbeg(j)-1)/(h(nbeg(j)-1)+oldhl(j))* 2 (qw(nbeg(j)+1)-dw(nbeg(j)+1)+dw(nbeg(j)-1)) if (qw(nbeg(j)).LT.0.) 1 qw(nbeg(j))=qw(nbeg(j)+1)*uw(nbeg(j)+1)*t/ 2 (uw(nbeg(j)+1)*t+oldhl(j)) end if else if( dw(nbeg(j)+1)-oldl(j) .GE. qw(nbeg(j)+1)/2. )then qw(nbeg(j))=qw(nbeg(j)+1)-dw(nbeg(j)+1)+oldl(j) else qw(nbeg(j))=qw(nbeg(j)+1)*.5 end if end if c old method uw a characteristic C if (uw(nbeg(j)+1).LT.0)then if ((uw(nbeg(j)+1).LT.0).AND.(qw(nbeg(j)).GT.grnd))then vw(nbeg(j))=vw(nbeg(j)+1) uw(nbeg(j))=uw(nbeg(j)+1) C froude near the beach if (abs(uw(nbeg(j))).GE.1.0*sqrt(GA*qw(nbeg(j)))) then uw(nbeg(j)) = sign(1.0*sqrt(GA*qw(nbeg(j))),uw(nbeg(j))) endif if (abs(vw(nbeg(j))).GE.1.0*sqrt(GA*qw(nbeg(j)))) then vw(nbeg(j)) = sign(1.0*sqrt(GA*qw(nbeg(j))),vw(nbeg(j))) endif else vw(nbeg(j))=0. uw(nbeg(j))=0. end if C if (qw(nbeg(j)).GT.grnd.AND.qw(nbeg(j)+1).GT.grnd)then c vw(nbeg(j))=vw(nbeg(j)+1)*qw(nbeg(j)+1)/qw(nbeg(j)) C vw(nbeg(j))=vw(nbeg(j)+1) C uw(nbeg(j))=uw(nbeg(j)+1)*qw(nbeg(j)+1)/qw(nbeg(j)) C if (abs(uw(nbeg(j))).GE.1.0*sqrt(GA*qw(nbeg(j)))) then C uw(nbeg(j)) = sign(1.0*sqrt(GA*qw(nbeg(j))),uw(nbeg(j))) C endif C if (abs(vw(nbeg(j))).GE.1.0*sqrt(GA*qw(nbeg(j)))) then C vw(nbeg(j)) = sign(1.0*sqrt(GA*qw(nbeg(j))),vw(nbeg(j))) C endif C else C vw(nbeg(j))=0. C uw(nbeg(j))=0. C end if dw(nbeg(j))=oldl(j) h(nbeg(j))=oldhl(j) qwl(j)=qw(nbeg(j)) if (j.GT.1) then if (nend(j-1).EQ.nbeg(j)) then qw(nbeg(j))=max(qwl(j),qwr(j-1)) uw(nbeg(j))=0. vw(nbeg(j))=0. end if end if 100 continue if (wet1) then qw(1)=(u1(1)-q1(1))**2/(16.*GA) uw(1)=(u1(1)+q1(1))*.5 vw(1)=v1(1) end if if (wetn) then qw(n)=(u1(n)-q1(n))**2/(16.*GA) uw(n)=(u1(n)+q1(n))*.5 vw(n)=v1(n) end if 999 return end