subroutine edgevl(u,v,ud,vd,ux,vx,uy,vy,nx,ny) c**************************************************************** c c determine the edge values of velocity by upwinding, c before making the pressure correction. c c (u,v) -> velocity field c (ux,vx) (uy,vy) <- 1), constructed first derivatives c <- 2), edge values c (ud,vd) <- laplacian of (u,v) c temp space: c (up,vp) (uq,vq) <- first derivatives after c 2nd order upwinding used c in transverse directions. c c**************************************************************** integer i,j,nx,ny,nxm,nym,mxym include 'gridsz' double precision u(nx,ny),v(nx,ny), + ud(nx,ny),vd(nx,ny), + ux(0:nx,0:ny),vx(0:nx,0:ny), + uy(0:nx,0:ny),vy(0:nx,0:ny) double precision up(nxm,nym),vp(nxm,nym), + uq(nxm,nym),vq(nxm,nym) double precision uld,vld,h,hinv,rey,ul,ur,vl,vr,ub,ut,vb,vt, + roa,hf,hdt,dt,uby,vby,utop,ubot,grav parameter (hf=0.5D0) common/timest/dt common/meshsz/h,hinv common/reynod/rey common/density/roa(0:nxm+1,0:nym+1) common/gravity/grav common/uvboud/uby(nym),vby(nym),utop,ubot c initialize the local arrays: do 2 j=1,ny do 2 i=1,nx up(i,j) = 0.D0 vp(i,j) = 0.D0 uq(i,j) = 0.D0 vq(i,j) = 0.D0 2 continue c hdt = hf*dt c c approximate the first derivatives and laplacian of velocities: call evalim(u,ux,uy,uby,utop,ubot,nx,ny) call evalim(v,vx,vy,vby,0.D0,0.D0,nx,ny) call evalap(u,uby,utop,ubot,ud,nx,ny) call evalap(v,vby,0.D0,0.D0,vd,nx,ny) c c *************************************************************** c construct the tangential derivatives along cell edges c c 1): along vertical edges, do 10 j=1,ny up(1,j) = ux(1,j) vp(1,j) = vx(1,j) up(nx,j) = ux(nx,j) vp(nx,j) = vx(nx,j) 10 continue do 20 j=1,ny do 20 i=2,nx-1 call upwsec( u(i,j), u(i-1,j), u(i+1,j), + v(i,j), v(i-1,j), v(i+1,j), + ux(i,j), ux(i-1,j), ux(i+1,j), + vx(i,j), vx(i-1,j), vx(i+1,j), + up(i,j), vp(i,j) ) 20 continue c 2): along horizontal edges, do 40 i=1,nx uq(i,1) = uy(i,1) vq(i,1) = vy(i,1) uq(i,ny) = uy(i,ny) vq(i,ny) = vy(i,ny) 40 continue do 60 j=2,ny-1 do 60 i=1,nx call upwsec( v(i,j), v(i,j-1), v(i,j+1), + u(i,j), u(i,j-1), u(i,j+1), + vy(i,j), vy(i,j-1), vy(i,j+1), + uy(i,j), uy(i,j-1), uy(i,j+1), + vq(i,j), uq(i,j) ) 60 continue c *********************************************************** c construct time derivatives c do 70 j=1,ny do 70 i=1,nx uld = ud(i,j)/rey/roa(i,j) vld = vd(i,j)/rey/roa(i,j) uq(i,j) = uld - hinv*( u(i,j)*ux(i,j) + v(i,j)*uq(i,j) ) up(i,j) = uld - hinv*( u(i,j)*up(i,j) + v(i,j)*uy(i,j) ) vq(i,j) = vld - hinv*( u(i,j)*vx(i,j) + v(i,j)*vq(i,j) ) - grav vp(i,j) = vld - hinv*( u(i,j)*vp(i,j) + v(i,j)*vy(i,j) ) - grav 70 continue c *********************************************************** c Taylor extrapolation and upwinding to obtain edge values, c vertical: do 80 j=1,ny do 80 i=1,nx-1 ul = u(i,j) + hdt*uq(i,j) + ux(i,j)*hf ur = u(i+1,j) + hdt*uq(i+1,j) - ux(i+1,j)*hf vl = v(i,j) + hdt*vq(i,j) + vx(i,j)*hf vr = v(i+1,j) + hdt*vq(i+1,j) - vx(i+1,j)*hf call upwind(ul,ur,vl,vr,ux(i,j),vx(i,j)) 80 continue c horizontal: do 85 j=1,ny-1 do 85 i=1,nx ub = u(i,j) + hdt*up(i,j) + uy(i,j)*hf ut = u(i,j+1) + hdt*up(i,j+1) - uy(i,j+1)*hf vb = v(i,j) + hdt*vp(i,j) + vy(i,j)*hf vt = v(i,j+1) + hdt*vp(i,j+1) - vy(i,j+1)*hf call upwind(vb,vt,ub,ut,vy(i,j),uy(i,j)) 85 continue c--- expanding the array (uy,vy) and (ux,vx) by incorporating c boundary values. do 90 i=1,nx uy(i,0) = ubot vy(i,0) = 0.D0 uy(i,ny) = utop vy(i,ny) = 0.D0 90 continue do 95 j=1,ny ux(0,j) = uby(j) vx(0,j) = 0.D0 ux(nx,j) = uby(j) vx(nx,j) = 0.D0 95 continue end subroutine evalap(u,uby,utop,ubot,ud,nx,ny) c********************************************************************* c c compute the Laplacian of u and add its effect to the velocity. c central differencing in the interior and one-sided differencing c on the boundaries c c u -> u at cell centers c uby -> u at the vertical walls c ud <- Laplacian of u. c c********************************************************************* integer nx,ny,i,j double precision u(nx,ny),ud(nx,ny),uby(ny),utop,ubot, + c2,c3,cb,sv,tn,h,hinv parameter (c2=2.D0,c3=-0.2D0,cb=3.2D0,sv=7.D0,tn=10.D0) common/meshsz/h,hinv c .......................................................... c--- interior do 10 j=2,ny-1 do 10 i=2,nx-1 ud(i,j) = u(i+1,j) + u(i-1,j) + + u(i,j+1) + u(i,j-1) - 4.D0*u(i,j) 10 continue c .......................................................... c boundaries do 20 i=2,nx-1 ud(i,1) = u(i+1,1) + u(i-1,1) + cb*ubot + + c2*u(i,2) + c3*u(i,3) - sv*u(i,1) ud(i,ny) = u(i+1,ny) + u(i-1,ny) + cb*utop + + c2*u(i,ny-1) + c3*u(i,ny-2) - sv*u(i,ny) 20 continue do 30 j=2,ny-1 ud(1,j) = u(1,j+1) + u(1,j-1) + cb*uby(j) + + c2*u(2,j) + c3*u(3,j) - sv*u(1,j) ud(nx,j) = u(nx,j+1) + u(nx,j-1) + cb*uby(j) + + c2*u(nx-1,j) + c3*u(nx-2,j) - sv*u(nx,j) 30 continue c corners ud(1,1) = c2*(u(1,2)+u(2,1)) + + cb*(uby(1)+ubot) + + c3*(u(1,3)+u(3,1)) - + tn*u(1,1) ud(nx,1) = c2*(u(nx-1,1)+u(nx,2)) + + cb*(uby(1)+ubot) + + c3*(u(nx-2,1)+u(nx,3)) - + tn*u(nx,1) ud(1,ny) = c2*(u(1,ny-1)+u(2,ny)) + + cb*(uby(ny)+utop) + + c3*(u(1,ny-2)+u(3,ny)) - + tn*u(1,ny) ud(nx,ny) = c2*(u(nx-1,ny)+u(nx,ny-1)) + + cb*(uby(ny)+utop) + + c3*(u(nx-2,ny)+u(nx,ny-2)) - + tn*u(nx,ny) do 40 j=1,ny do 40 i=1,nx ud(i,j) = hinv*hinv*ud(i,j) 40 continue end subroutine evalim(u,ux,uy,uby,utop,ubot,nx,ny) c******************************************************************* c c construct the first derivatives with limiters c c u -> u at cell centers c ux,uy <- constructed first derivatives c c******************************************************************* integer i,j,nx,ny double precision u(nx,ny),ux(0:nx,0:ny),uy(0:nx,0:ny),uby(ny) double precision blimit,slimit,utop,ubot c .......................................................... c determine ux. do 10 j=1,ny ux(1,j) = blimit( uby(j), u(1,j), u(2,j) ) ux(nx,j) =-blimit( uby(j), u(nx,j), u(nx-1,j) ) 10 continue do 20 j=1,ny do 20 i=2,nx-1 ux(i,j) = slimit( u(i+1,j)-u(i,j), u(i,j)-u(i-1,j) ) 20 continue c .......................................................... c determine uy. do 60 i=1,nx uy(i,1) = blimit( ubot, u(i,1), u(i,2) ) uy(i,ny) =-blimit( utop, u(i,ny), u(i,ny-1) ) 60 continue do 70 j=2,ny-1 do 70 i=1,nx uy(i,j) = slimit( u(i,j+1)-u(i,j), u(i,j)-u(i,j-1) ) 70 continue end