www.pudn.com > cylindar.zip > cylindar.f90, change:2013-11-19,size:21061b
MODULE cellConst
INTEGER, PARAMETER:: fluid = 0, wall = 1, inlet = 10, outlet = 11, incylinder=100
END MODULE cellConst
!========================================================
! Lattice constants for the D2Q9 lattice
!========================================================
MODULE D2Q9Const
! D2Q9 Weights
REAL(8),PARAMETER:: t(0:8) = (/4.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0,1.0d0/9.0d0&
&,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0,1.0d0/36.0d0/)
! D2Q9 Directions
INTEGER:: v(0:8,0:1)
! = (/(/0,1,0,-1,0,1,-1,-1,1/),(/0,0,1,0,-1,1,1,-1,-1/)/)
INTEGER, PARAMETER:: opposite(0:8) = (/0,3,4,1,2,7,8,5,6/)
END MODULE D2Q9Const
!========================================================
! Constants for simulation setup
!========================================================
MODULE simParam
!INTEGER, PARAMETER:: xdim = 4500
!INTEGER, PARAMETER:: ydim = 800
!INTEGER, PARAMETER:: obstX = xdim/5
!INTEGER, PARAMETER:: obstY = ydim/2
!INTEGER, PARAMETER:: obstR = ydim/10+1
INTEGER, PARAMETER:: tmax = 100000
REAL(8), PARAMETER:: umax = 0.2d0
REAL(8), PARAMETER:: Re1 = 500.0d0
INTEGER, PARAMETER:: xdim = 2200
INTEGER, PARAMETER:: ydim = 320
INTEGER, PARAMETER:: obstX = 150
INTEGER, PARAMETER:: obstY = 155
INTEGER, PARAMETER:: obstR = 50
END MODULE simParam
!========================================================
! The main program, implementing a flow past a cylinder
!========================================================
PROGRAM unsteady
USE simParam, ONLY: xdim, ydim, tmax
USE D2Q9Const, ONLY: v
implicit none
REAL(8):: omega, time1, time2, timeTot
REAL(8), dimension(:,:,:), allocatable:: f, fEq, u, meq
REAL(8), dimension(:,:), allocatable:: rho, uSqr, streamline, vorticity
INTEGER, dimension(:,:), allocatable:: image
INTEGER:: tStep=0
allocate(f(ydim,xdim,0:8))
allocate(fEq(ydim,xdim,0:8))
allocate(meq(ydim,xdim,0:8))
allocate(u(ydim,xdim,0:1))
allocate(uSqr(ydim,xdim))
allocate(rho(ydim,xdim))
allocate(image(ydim,xdim))
allocate(streamline(ydim,xdim))
allocate(vorticity(ydim,xdim))
CALL constructImage(image)
CALL computeOmega(omega)
CALL writeInput(omega)
CALL initMacro(rho,u,uSqr,image)
CALL computeFeq(fEq,rho,u,uSqr)
f = fEq
CALL computevorticity(u, vorticity, image)
!CALL computestreamline(u,streamline,image)
CALL writeOutput(u,tstep,vorticity)
timeTot = 0.0d0
do tStep = 1, tmax
CALL CPU_TIME(time1)
CALL inletOutlet(f,rho,u,image)
CALL boundaries(f,image)
CALL computeMacros(f,rho,u,uSqr)
! CALL computeFeq(fEq,rho,u,uSqr)
! CALL collide(f,fEq,omega,image)
CALL computemeq(meq,rho,u,uSqr)
CALL coll_mrt(f,meq,omega,image)
CALL streaming(f,v,image)
! CALL writeImage(image)
IF (MOD(tstep, 1000)==0) THEN
!CALL computestreamline(u,streamline,image)
CALL computevorticity(u, vorticity, image)
CALL writeOutput(u,tstep,vorticity)
END IF
CALL CPU_TIME(time2)
timeTot = timeTot + (time2-time1)
end do
!CALL computestreamline(u,streamline,image)
CALL computevorticity(u, vorticity, image)
CALL writeOutput(u,tmax,vorticity)
write(*,*) dble(tmax) * (dble(ydim * xdim)) / timeTot ,'cells per second'
deallocate(f)
deallocate(fEq)
deallocate(meq)
deallocate(u)
deallocate(uSqr)
deallocate(rho)
deallocate(image)
deallocate(streamline)
deallocate(vorticity)
END PROGRAM unsteady
!========================================================
! Compute the relaxation PARAMETER from the Reynolds number
!========================================================
SUBROUTINE computeOmega(omega)
USE simParam, ONLY: Re1,umax,obstR
implicit none
REAL(8), INTENT(INOUT):: omega
REAL(8):: nu
nu = umax * 2.0d0 * dble(obstR) / Re1
omega = 1.0d0 / (3.0d0*nu+0.5d0)
END SUBROUTINE computeOmega
!========================================================
! Construct an array the defines the flow geometry
!========================================================
SUBROUTINE constructImage(image)
USE cellConst
USE simParam, ONLY: xdim, ydim, obstX, obstY, obstR
USE D2Q9Const, ONLY: v
implicit none
INTEGER, INTENT(INOUT):: image(ydim,xdim)
INTEGER:: x,y
v(0:8,0) = (/0,1,0,-1,0,1,-1,-1,1/)
v(0:8,1) = (/0,0,1,0,-1,1,1,-1,-1/)
image = fluid
image(:,1) = inlet
image(:,xdim) = outlet
image(1,:) = wall
image(ydim,:) = wall
do x=1,xdim
do y=1,ydim
if(((x-obstX)**2+(y-obstY)**2) <= (obstR**2)) image(y,x)=incylinder
enddo
enddo
END SUBROUTINE constructImage
!========================================================
! Initialize the simulation to Poiseuille profile at
! an equilibrium distribution
!========================================================
SUBROUTINE initMacro(rho,u,uSqr,image)
USE simParam, ONLY: xdim, ydim, umax
USE cellConst, ONLY: incylinder,wall
implicit none
REAL(8), INTENT(INOUT):: rho(ydim,xdim), u(ydim,xdim,0:1), uSqr(ydim,xdim)
!REAL(8):: uProf
!REAL(8):: randomnumber
INTEGER, INTENT(IN):: image(ydim,xdim)
INTEGER:: x,y
DO x=1,xdim
do y=1,ydim
IF (image(y,x)/=incylinder .and. image(y,x)/=wall) THEN
If (MOD(x, 10)==0 .and. MOD(y, 10)==0) THEN
!CALL RANDOM_SEED()
!CALL RANDOM_NUMBER(randomnumber)
u(y,x,0)=umax*0.5
ELSE
u(y,x,0)=umax
END IF
!u(y,x,0)=uProf(y)
!u(y,x,0)=umax
u(y,x,1)=0.0d0
END IF
enddo
END DO
rho =1.0d0
uSqr=u(:,:,0)*u(:,:,0)+u(:,:,1)*u(:,:,1)
END SUBROUTINE initMacro
!========================================================
! Compute equilibrium distribution
!========================================================
SUBROUTINE computeFeq(fEq,rho,u,uSqr)
USE D2Q9COnst, ONLY: t, v
USE simParam, ONLY: xdim, ydim
implicit none
REAL(8), INTENT(IN):: rho(ydim,xdim), uSqr(ydim,xdim), u(ydim,xdim,0:1)
REAL(8), INTENT(INOUT):: fEq(ydim,xdim,0:8)
INTEGER:: i, x, y
REAL(8):: uxy
do i = 0, 8
do x = 1, xdim
do y = 1, ydim
uxy=u(y,x,0)*v(i,0)+u(y,x,1)*v(i,1)
fEq(y,x,i)=t(i)*rho(y,x)*(1.0d0+3.0d0*uxy+4.5d0*uxy*uxy-1.5d0*uSqr(y,x))
enddo
enddo
enddo
END SUBROUTINE computeFeq
!========================================================
! Compute density and velocity from distribution functions
!========================================================
SUBROUTINE computeMacros(f,rho,u,uSqr)
USE simParam, ONLY: xdim, ydim
implicit none
REAL(8), INTENT(IN):: f(ydim,xdim,0:8)
REAL(8), INTENT(INOUT):: u(ydim,xdim,0:1), rho(ydim, xdim), uSqr(ydim, xdim)
INTEGER:: x,y
do x = 1, xdim
do y = 1, ydim
rho(y,x)=f(y,x,0)+f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4)+f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8)
u(y,x,0)=(f(y,x,1)-f(y,x,3)+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8))/rho(y,x)
u(y,x,1)=(f(y,x,2)-f(y,x,4)+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8))/rho(y,x)
uSqr(y,x)=u(y,x,0) * u(y,x,0) + u(y,x,1) * u(y,x,1)
enddo
enddo
END SUBROUTINE computeMacros
!========================================================
! Implement Bounce-back on upper/lower boundaries
!========================================================
SUBROUTINE boundaries(f,image)
USE D2Q9Const, ONLY: opposite
USE cellConst, ONLY: wall
USE simParam, ONLY: xdim, ydim
implicit none
INTEGER, INTENT(IN):: image(ydim,xdim)
REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8)
REAL(8):: fTmp(0:8)
INTEGER:: i, x, y
do x = 1, xdim
do y = 1, ydim
if (image(y,x) == wall) then
do i = 0, 8
fTmp(i) = f(y,x,opposite(i))
end do
do i = 0, 8
f(y,x,i) = fTmp(i)
end do
end if
end do
end do
END SUBROUTINE boundaries
!========================================================
! Use Zou/He boundary condition to implement Dirichlet
! boundaries on inlet/outlet
!========================================================
SUBROUTINE inletOutlet(f,rho,u,image)
USE cellConst, ONLY: inlet, outlet
USE simParam
implicit none
REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8),u(ydim,xdim,0:1),rho(ydim,xdim)
INTEGER, INTENT(IN):: image(ydim,xdim)
REAL(8):: uProf
INTEGER:: x,y
do x=1,xdim
do y=1,ydim
if(image(y,x) == inlet) then
u(y,x,0) = uProf(y)
u(y,x,1) = 0.0d0
CALL inletZou(f(y,x,:),u(y,x,:),rho(y,x))
elseif(image(y,x) == outlet) then
u(y,x,0) = uProf(y)
u(y,x,1) = 0.0d0
CALL outletZou(f(y,x,:),u(y,x,:),rho(y,x))
endif
enddo
enddo
CONTAINS
!========================================================
! Zou/He boundary on inlet
!========================================================
SUBROUTINE inletZou(f,u,rho)
implicit none
REAL(8), INTENT(INOUT):: f(0:8),rho
REAL(8), INTENT(IN):: u(0:1)
REAL(8):: fInt, fInt2
fInt = f(0) + f(2) + f(4)
fInt2 = f(3) + f(6) + f(7)
rho = (fInt + 2.0d0 * fInt2) / (1.0d0 - u(0))
CALL zouWestWall(f,rho,u)
END SUBROUTINE inletZou
SUBROUTINE zouWestWall(f,rho,u)
implicit none
REAL(8), INTENT(INOUT):: f(0:8)
REAL(8), INTENT(IN):: rho, u(0:1)
REAL(8):: fDiff, rhoUx, rhoUy
fDiff = 0.5d0 * (f(2) - f(4))
rhoUx = rho * u(0) / 6.0d0
rhoUy = 0.5d0 * rho * u(1)
f(1) = f(3) + 4.0d0 * rhoUx
f(5) = f(7) - fDiff + rhoUx + rhoUy
f(8) = f(6) + fDiff + rhoUx - rhoUy
END SUBROUTINE zouWestWall
!========================================================
! Zou/He boundary on outlet
!========================================================
SUBROUTINE outletZou(f,u,rho)
implicit none
REAL(8), INTENT(INOUT):: f(0:8),rho,u(0:1)
REAL(8):: fInt, fInt2
fInt = f(0) + f(2) + f(4)
fInt2 = f(1) + f(8) + f(5)
rho = (fInt + 2.0d0 * fInt2) / (1.0d0 + u(0))
CALL zouEastWall(f,rho,u)
END SUBROUTINE outletZou
SUBROUTINE zouEastWall(f,rho,u)
implicit none
REAL(8), INTENT(INOUT):: f(0:8)
REAL(8), INTENT(IN):: rho, u(0:1)
REAL(8):: fDiff, rhoUx, rhoUy
fDiff = 0.5d0 * (f(2) - f(4))
rhoUx = rho * u(0) / 6.0d0
rhoUy = 0.5d0 * rho * u(1)
f(3) = f(1) - 4.0d0 * rhoUx
f(7) = f(5) + fDiff - rhoUx - rhoUy
f(6) = f(8) - fDiff - rhoUx + rhoUy
END SUBROUTINE zouEastWall
END SUBROUTINE inletOutlet
!========================================================
! Computation of Poiseuille profile for the inlet/outlet
!========================================================
FUNCTION uProf(y)
USE simParam, ONLY: ydim, umax
implicit none
INTEGER, INTENT(IN):: y
REAL(8):: radius, uProf
radius = dble(ydim-1) * 0.5d0
uProf = -umax * ((abs(1 - dble(y-1) / radius))**2 - 1.0d0)
END FUNCTION uProf
!========================================================
! Streaming step: the population functions are shifted
! one site along their corresponding lattice direction
! (no temporary memory is needed)
!========================================================
SUBROUTINE streaming(f,v,image)
USE simParam
USE cellConst, ONLY: incylinder
implicit none
INTEGER::y,x,yd,xd,k
REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8)
INTEGER, INTENT(IN):: v(0:8,0:1)
REAL(8):: f_temp(ydim,xdim,0:8)
INTEGER, INTENT(IN):: image(ydim,xdim)
! -------------------------------------
f_temp = f
DO k = 0,8
DO x = 1,xdim
DO y = 1,ydim
yd = y - v(k,1)
xd = x - v(k,0)
IF (yd>=1 .and. yd<=ydim .and. xd>=1 .and. xd<=xdim .and. image(y,x) /= incylinder) THEN
f(y,x,k)=f_temp(yd,xd,k)
END IF
END DO
END DO
END DO
END SUBROUTINE streaming
!========================================================
! LBGK collision step
!========================================================
SUBROUTINE collide(f,fEq,omega,image)
USE simParam, ONLY: xdim, ydim
USE cellConst, ONLY: incylinder
implicit none
INTEGER, INTENT(IN):: image(ydim,xdim)
REAL(8), INTENT(IN):: fEq(ydim,xdim,0:8), omega
REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8)
INTEGER:: x,y,i
do i=0,8
do x=1,xdim
do y=1,ydim
if(image(y,x) /= incylinder) f(y,x,i)=(1.0d0-omega)*f(y,x,i)+omega*feq(y,x,i)
enddo
enddo
enddo
END SUBROUTINE collide
!========================================================
! Write the components of the velocity to a text file,
! with indices (x,y)
!========================================================
SUBROUTINE writeOutput(u,tStep,vorticity)
USE simParam, ONLY: xdim, ydim, Re1, umax
implicit none
INTEGER, INTENT(IN):: tStep
REAL(8), INTENT(IN):: u(ydim,xdim,0:1)
REAL(8), INTENT(IN):: vorticity(ydim,xdim)
INTEGER:: x,y
character (LEN=100):: fileName
101 format(2i10,3f20.10)
write(fileName, '(I7,A8,F7.1,A8,F7.3,A8,I7,A8,I7)') tStep, " RE = ", Re1,&
" umax = ", umax, " xdim = ", xdim, " ydim = ", ydim
fileName = adjustl(fileName)
open(14,file='outputUx_'//trim(fileName)//'.plt')
write(14,*)'TITLE=LBfor'
write(14,*)'VARIABLES="x","y","u","v","w"'
write(14,*)'ZONE F=POINT, I=',ydim,', J=',xdim
do x=1,xdim
do y=1,ydim
write(14,101) x,y,u(y,x,0),u(y,x,1),vorticity(y,x)
enddo
enddo
close(14)
END SUBROUTINE writeOutput
!========================================================
! Write the flow geometry to a file
!========================================================
SUBROUTINE writeImage(image)
USE simParam, ONLY: xdim, ydim
implicit none
INTEGER, INTENT(IN):: image(ydim,xdim)
INTEGER:: x,y
open(13,file='outputImage.dat')
do x=1, xdim
do y=1, ydim
write(13,102) image(y,x)
end do
end do
102 format(3i10)
close(15)
END SUBROUTINE writeImage
!========================================================
! Print out simulation PARAMETERs to screen
!========================================================
SUBROUTINE writeInput(omega)
USE simParam
implicit none
REAL(8), INTENT(IN):: omega
write(*,*) 'xdim = ', xdim
write(*,*) 'ydim = ', ydim
write(*,*) 'Obstacle X = ', obstX
write(*,*) 'Obstacle Y = ', obstY
write(*,*) 'Obstacle Radius = ', obstR
write(*,*) 'tmax = ', tmax
write(*,*) 'umax = ', umax
write(*,*) 'Re = ', Re1
write(*,*) 'omega = ', omega
END SUBROUTINE writeInput
SUBROUTINE computestreamline(u, streamline,image)
USE simParam, ONLY: xdim, ydim, obstY
USE cellConst, ONLY: incylinder
implicit none
REAL(8), INTENT(IN):: u(ydim,xdim,0:1)
REAL(8), INTENT(INOUT):: streamline(ydim,xdim)
INTEGER, INTENT(IN):: image(ydim,xdim)
INTEGER:: x,y
streamline(1,1)=0.0d0
DO x = 2, xdim
streamline(1,x) = streamline(1, x-1) - u(1, x-1, 1) * 1.0
END DO
DO x = 1, xdim
DO y = 2, obstY
If(image(y-1,x) /= incylinder) THEN
streamline(y,x) = streamline(y-1,x) + u(y-1, x, 0) * 1.0
END IF
END DO
END DO
streamline(ydim,1)=0.0d0
DO x = 2, xdim
streamline(ydim,x) = streamline(ydim, x-1) - u(ydim, x-1, 1) * 1.0
END DO
DO x = 1, xdim
DO y = ydim-1, obstY,-1
If(image(y,x) /= incylinder) THEN
streamline(y,x) = streamline(y+1,x) + u(y+1, x, 0) * 1.0
END IF
END DO
END DO
END SUBROUTINE computestreamline
SUBROUTINE computemeq(meq,rho,u,uSqr)
USE simParam, ONLY: xdim, ydim
implicit none
REAL(8), INTENT(IN):: rho(ydim,xdim), uSqr(ydim,xdim), u(ydim,xdim,0:1)
REAL(8), INTENT(INOUT):: meq(ydim,xdim,0:8)
INTEGER:: i, x, y
do i = 0, 8
do x = 1, xdim
do y = 1, ydim
SELECT CASE(i)
CASE(0)
meq(y,x,i)=rho(y,x)
CASE(1)
meq(y,x,i)=rho(y,x)*(-2.0d0+3.0d0*usqr(y,x))
CASE(2)
meq(y,x,i)=rho(y,x)*(1.0d0-3.0d0*usqr(y,x))
CASE(3)
meq(y,x,i)=rho(y,x)*u(y,x,0)
CASE(4)
meq(y,x,i)=-rho(y,x)*u(y,x,0)
CASE(5)
meq(y,x,i)=rho(y,x)*u(y,x,1)
CASE(6)
meq(y,x,i)=-rho(y,x)*u(y,x,1)
CASE(7)
meq(y,x,i)=rho(y,x)*(u(y,x,0)*u(y,x,0)-u(y,x,1)*u(y,x,1))
CASE(8)
meq(y,x,i)=rho(y,x)*u(y,x,0)*u(y,x,1)
CASE DEFAULT
meq(y,x,i)=0.0d0
END SELECT
END DO
enddo
enddo
END SUBROUTINE computemeq
SUBROUTINE coll_mrt(f,meq,omega,image)
USE simParam, ONLY: xdim, ydim
USE cellConst, ONLY: incylinder
implicit none
INTEGER:: x,y,i
INTEGER, INTENT(IN):: image(ydim,xdim)
REAL(8), INTENT(IN):: meq(ydim,xdim,0:8), omega
REAL(8), INTENT(INOUT):: f(ydim,xdim,0:8)
REAL(8):: m(0:8),s(0:8)
REAL(8):: d(0:8)=(/9.0d0,36.0d0,36.0d0,6.0d0,12.0d0,6.0d0,12.0d0,4.0d0,4.0d0/)
s(7:8)=omega
s(0)=0.0d0
s(1)=1.6d0
s(2)=1.8d0
s(3:5:2)=0.0d0
s(4:6:2)=8*(2-s(7))/(8-s(7))
do x=1,xdim
do y=1,ydim
If(image(y,x) /= incylinder) THEN
m(0)=f(y,x,0)+f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4)&
+f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8)
m(1)=-4.0d0*f(y,x,0)-f(y,x,1)-f(y,x,2)-f(y,x,3)-f(y,x,4)&
+2.0d0*(f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8))
m(2)=4.0d0*f(y,x,0)-2.0d0*(f(y,x,1)+f(y,x,2)+f(y,x,3)+f(y,x,4))&
+f(y,x,5)+f(y,x,6)+f(y,x,7)+f(y,x,8)
m(3)=f(y,x,1)-f(y,x,3)+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8)
m(4)=-2.0d0*(f(y,x,1)-f(y,x,3))+f(y,x,5)-f(y,x,6)-f(y,x,7)+f(y,x,8)
m(5)=f(y,x,2)-f(y,x,4)+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8)
m(6)=-2.0d0*(f(y,x,2)-f(y,x,4))+f(y,x,5)+f(y,x,6)-f(y,x,7)-f(y,x,8)
m(7)=f(y,x,1)-f(y,x,2)+f(y,x,3)-f(y,x,4)
m(8)=f(y,x,5)-f(y,x,6)+f(y,x,7)-f(y,x,8)
DO i=0,8
m(i)=m(i)-s(i)*(m(i)-meq(y,x,i))
m(i)=m(i)/d(i)
END DO
f(y,x,0)=m(0)-4*(m(1)-m(2))
f(y,x,1)=m(0)-m(1)-2*(m(2)+m(4))+m(3)+m(7)
f(y,x,2)=m(0)-m(1)-2*(m(2)+m(6))+m(5)-m(7)
f(y,x,3)=m(0)-m(1)-2*(m(2)-m(4))-m(3)+m(7)
f(y,x,4)=m(0)-m(1)-2*(m(2)-m(6))-m(5)-m(7)
f(y,x,5)=m(0)+m(1)+m(1)+m(2)+m(3)+m(4)+m(5)+m(6)+m(8)
f(y,x,6)=m(0)+m(1)+m(1)+m(2)-m(3)-m(4)+m(5)+m(6)-m(8)
f(y,x,7)=m(0)+m(1)+m(1)+m(2)-m(3)-m(4)-m(5)-m(6)+m(8)
f(y,x,8)=m(0)+m(1)+m(1)+m(2)+m(3)+m(4)-m(5)-m(6)-m(8)
END IF
enddo
enddo
END SUBROUTINE coll_mrt
SUBROUTINE computevorticity(u, vorticity, image)
USE simParam, ONLY: xdim, ydim, obstY, obstX, obstR
USE cellConst, ONLY: wall
implicit none
REAL(8), INTENT(IN):: u(ydim,xdim,0:1)
REAL(8), INTENT(INOUT):: vorticity(ydim,xdim)
INTEGER, INTENT(IN):: image(ydim,xdim)
INTEGER:: x,y
DO x = 2, xdim
DO y = 2, obstY-obstR-1
vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0- (u(y, x, 0) - u(y-1, x, 0)) / 1.0
END DO
END DO
DO x = 2, xdim
DO y = obstY-obstR-1, obstY+obstR+1
If(image(y,x) /= wall) THEN
vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0 - (u(y, x, 0) - u(y-1, x, 0)) / 1.0
END IF
END DO
END DO
DO x = 2, xdim
DO y = obstY+obstR+1, ydim-1
vorticity(y,x) = (u(y, x, 1) - u(y, x-1, 1)) / 1.0 - (u(y, x, 0) - u(y-1, x, 0)) / 1.0
END DO
END DO
END SUBROUTINE computevorticity