subMuddx.f90 3.13 KB
``````
subroutine computeRHS(stage, rhs)
! Set the ringht and side field
implicit none
integer i,j
integer stage, nx,ny
real dlx,dly,Pi , x, y
real rhs(nx,ny)

do i=1,nx
do j=1,ny
x = float(i-1)*dx
y = float(j-1)*dy

if (stage == 1)
rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * sin(2*Pi*y)
else if (stage == 2)
rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * sin(2*Pi*y)
else if (stage ==3)
rhs(i,j)    =  -8* Pi * cos(2*Pi * x) * cos(Pi*y)
end if
end do
end do
end

subroutine setBndy(stage,psi)
! set boundary condition for a specific stage
implicit none
integer i,j
integer stage, nx, ny
real dlx,dly,Pi , x, y
real psi(nx,ny)

! ---------------- zeta ----------------
! set boundary condition flags
if (stage == 1) then
! set specified boundaries in phi
x = 0
do j=1,ny
y = float(j-1)*dy
psi(1,j) = sin(2*Pi*y)
end do

x = Lx
do j=1,ny
y = float(j-1)*dy
psi(nx,j) = sin(2*Pi*y)
end do

y = 0
do i=1,nx
x = float(i-1)*dx
psi(i,1) = 0.0
end do

y = Ly
do i=1,nx
x = float(i-1)*dx
psi(i,ny) = 0.0
end do

! ---------------- phi ----------------
! set boundary condition flags
else if (stage == 2) then
! set specified boundaries in phi
x = xa
do j=1,ny
y = float(j-1)*dy
psi(1,j) = sin(2*Pi*y)
end do

x = xb
do j=1,ny
y = float(j-1)*dy
psi(nx,j) = sin(2*Pi*y)
end do

y = yc
do i=1,nx
x = float(i-1)*dx
psi(i,1) = 0.0
end do

y = yd
do i=1,nx
x = float(i-1)*dx
psi(i,ny) = 0.0
end do

! --------------------- v ------------------------
! set boundary condition flags
else if (stage ==3) then
! set specified boundaries in phi
x = xa
do j=1,ny
y = yc+float(j-1)*dly
psi(1,j) = cos(Pi*y)
end do

x = xb
do j=1,ny
y = yc+float(j-1)*dly
psi(nx,j) = cos(Pi*y)
end do

y = yc
do i=1,nx
x = xa+float(i-1)*dlx
psi(i,1) = cos(2*pi*x)
end do

y = yd
do i=1,nx
x = xa+float(i-1)*dlx
psi(i,ny) =-cos(2*pi*x)
end do
!------------ else error ---------------
else
write(*,*) 'Wrong stage for setBndy'
call exit(2)
end if
end

subroutine bndc(kbdy,xory,alfa,gbdy)
! input mixed derivative b.c. to mud2
implicit none
integer kbdy
real xory,alfa,gbdy,x,y,pe,px,py,pxx,pyy
real xa,xb,yc,yd,tolmax,relmax
common/ftmud2/xa,xb,yc,yd,tolmax,relmax
real tau, nu

tau  = 0.1
nu = 1E-6

if (kbdy.eq.3) then  ! y=yc tboundary
x = xory
y = yc
alfa = 0
gbdy = 0
return
end if
if (kbdy.eq.4) then  ! y=yd boundary
y = yd
x = xory
alfa = 0
gbdy = tau/nu
return
end if
end

subroutine cof(x,y,cxx,cyy,cx,cy,ce)
! input pde coefficients at any grid point (x,y) in the solution region

implicit none
real x,y,cxx,cyy,cx,cy,ce
cxx = 1
cyy = 1
cx = 0.
cy = 0.
ce = 0.
return
end

!END MODULE
``````