c
program bvort
c
c Integrate the barotropic vorticity equation on a
c beta plane. This is meant to be simple and easily modified, and
c not highly efficient. Spatial differencing is second
c order; time stepping is by a leap-frog, and occasional trapezoidal
c correction in order to supress time splitting.
c
c The basic case is very much like Veronis' classic 1966 model study
c (Deep Sea Research, 1966, p 17-55) in that the domain is a square basin,
c and the circulation of a single layer is driven by a wind stress. The
c model can also be configured as a channel, or as a completely cyclic
c domain, however, neither of these cases has been checked carefully
c as of August 94.
c
c Output from the model goes primarily to disk files that are later read in
c and plotted by Matlab programs (btime.m, blook.m). A very small amount
c of data goes to the terminal during the integration.
c
c Written by Jim Price beginning in March 1994. Last modified in March 99.
c
c NOTICE: This model code is only partially tested and debugged and very
c likely has errors embedded within it. You are responsible for any and all
c results you may get.
c
c For compilation on an HP C180 (note use of double precision)
c
c f77 +O2 +DA2.0N +DS2.0A +autodblpad bvort.f -o bvort.x
c
c Rows and columns and the x and y directions are defined as follows:
c
c (1,ny) .... (nx,ny)
c . .
c . .
c . . . . . . . .
c . .
c . .
c (1,1) .... (nx,1)
c
c
c The variables stored in arrays are as follows:
c
c s is the mass transport streamfunction, the thing you are after!,
c sold and snew are the old and new time levels of s,
c st is the time rate of change of s ,
c u and v are the east and north velocty components, m/sec
c del2s is the Laplacian of s (the vorticity), del2st is its time rate change,
c kpa, kma etc. are defined below,
c taux and tauy are the east and north wind stress/density components,
c x and y are the east and north coordinates, m
c bc, del2p, etc are used to compute the surface pressure (defined below),
c bv, curlt etc. are used to hold the terms of the vorticity equation,
c sd24 is used to hold the diffusion term (not used at present).
c
c Units are dimensional, MKS.
c
parameter (nx=51, ny=51)
c
dimension s(nx,ny), sold(nx,ny), snew(nx,ny), st(nx,ny)
dimension u(nx,ny), v(nx,ny)
dimension del2s(nx,ny), del2st(nx,ny)
dimension kpa(nx), kma(nx), jma(ny), jpa(ny),
1 kpd(ny), kmd(ny), jpd(nx), jmd(nx)
dimension taux(nx,ny), tauy(nx,ny)
dimension x(nx), y(ny), f(ny)
dimension bc(nx,ny), del2p(nx,ny), rlp(nx,ny)
dimension bv(nx,ny), curlt(nx,ny), vdecay(nx,ny), rv(nx,ny),
1 rjsv(nx,ny), vort(nx,ny), spd(nx,ny)
c
integer cyclic
c
c These files are used to save results that are later read in by Matlab
c plotting programs. The names have ceased to have any meaning, and could be
c changed here and in the Matlab programs.
c
open(unit=1,file='basin1.dat',status='unknown')
open(unit=2,file='basin2.dat',status='unknown')
open(unit=3,file='basin3.dat',status='unknown')
open(unit=4,file='basin4.dat',status='unknown')
open(unit=14,file='basin5.dat',status='unknown')
open(unit=16,file='basin6.dat',status='unknown')
open(unit=11,file='what.dat',status='unknown')
open(unit=12,file='looktim.dat',status='unknown')
open(unit=10,file='basin10.dat',status='unknown')
c
ncx = nx/2 + 1
nxm = nx - 1
ncy = ny/2 + 1
nym = ny - 1
if(nym.eq.0) nym = 1
c
c To specify the model run there are a number of parameters than have to
c be input (or defaulted). The default case is a basin with nominal wind
c forcing and linear dissipation.
c
write (6,1900)
1900 format
1 (1x,/,' Are you a novice user of this program? (1=yes, def=0)')
nov = 0
read (5,*) nov
c
if(nov.eq.1) then
write (6,1920)
1920 format (1x,//,
1 1x,' This model integrates a barotropic vorticity ',/,
1 1x,' equation over a square domain. The default case is a',/,
4 1x,' wind-driven basin much like the classic Veronis (1966,',/,
5 1x,' DSR, pp 17) study (and which is the best single ',/,
6 1x,' reference for the model dynamics). Other model ',/,
6 1x,' configurations are also possible, including an',/,
6 1x,' east-west channel, or a completely symmetric domain. In',/,
7 1x,' the latter case, the streamfunction field can be ',/,
8 1x,' initializd with a blob.',/)
write (6,1921)
1921 format (
9 1x,' There is a very small list of output to the terminal.',/,
1 1x,' Most of the data goes to disk files for later plotting',/,
2 1x,' by Matlab programs (see the listing for details).',///,
2 1x,' The program will now ask for data that you can default',/,
3 1x,' to values in parentheses by entering commas.',/)
c
write (6,1901)
1901 format (1x,/,' RUNTIM defines the number of days to run,',//,
11x,' IWRTD is the interval (days) between writing output data ',/,
21x,' to disk files for later plotting',/)
end if
c
c runtim is the number of days to run, wrtt is the time interval between
c writing to the screen and wrtd is the interval for writing to the disk
c files.
c
runtim = 100.
wrtt = 0.5
wrtd = 10.
c
write (6,777)
777 format (1x,/,' Enter RUNTIM (def=100.), IWRTD (def=10.)')
read (5,*) runtim, wrtd
c
c The boundary conditions are specified via the following flags:
c
c noslip = 0 for free slip BCs on the sidewalls of a basin (the default)
c " = 1 for no slip BCs on the sidewalls of a basin (note: this
c option has not been fully implemented as of Aug 94)
c cyclic = 0 gives a fully enclosed domain (a square basin)
c " = 1 for cyclic BCs on east and west sides, and slip or noslip on the
c north and south sides according to the flags above
c (this makes an east-west channel)
c cyclic = 2 for cyclic BCs all the way around (probably makes sense only if
c beta effect is zero)
c Note: only the basin case (cyclic = 0) has been fully tested.
c
c if(nov.eq.1) then
c write (6,1902)
c 1902 format (1x,//,
c c ' Boundary conditions are specified by 2 flags:',//,
c 1 1x,' NOSLIP = 0 gives free slip boundary conditions on walls',/,
c 2 1x,' " = 1 gives no slip boundary conditions on walls,',//,
c 3 1x,' CYCLIC = 0 gives a basin (walls all around),',/,
c 4 1x,' " = 1 gives an east-west channel,',/,
c 5 1x,' " = 2 gives a completely symmetric (open) domain',/)
c end if
c
noslip = 0
cyclic = 0
c
c write (6,9922)
c 9922 format
c 1 (1x,' Enter bc flags: NOSLIP (def=0), CYCLIC (def=0)')
c read (5,*) noslip, cyclic
c
c Set up some arrays used to specify indices in difference equations.
c kpa, kma etc. are used in the advection terms; kpd, kmd are used in the
c diffusion terms.
c
c These are set up differently for solid or cyclic boundaries.
c
do 81 j=1,nx
jpa(j) = j + 1
jma(j) = j - 1
jpd(j) = j + 1
jmd(j) = j - 1
81 continue
jma(1) = 1
jpa(nx) = nx
jmd(1) = 2
jpd(nx) = nx -1
c
c If this is a channel or completely cyclic domain then
c
if(cyclic.eq.1.or.cyclic.eq.2) then
jma(1) = nx
jpa(nx) = 1
jmd(1) = nx
jpd(nx) = 1
end if
c
do 80 k=1,ny
kpa(k) = k + 1
kma(k) = k - 1
kpd(k) = k + 1
kmd(k) = k - 1
80 continue
kpa(ny) = ny
kma(1) = 1
kpd(ny) = ny - 1
kmd(1) = 2
c
c If this is a cyclic domain then
c
if(cyclic.eq.2) then
kpa(ny) = 1
kma(1) = ny
kpd(ny) = 1
kmd(1) = ny
end if
c
c Define the first and last points over which the interior calculation
c will be done.
c
jf = 2
jl = nx - 1
if(cyclic.eq.1.or.cyclic.eq.2) then
jf = 1
jl = nx
end if
c
kf = 2
kl = ny - 1
if(cyclic.eq.2) then
kf = 1
kl = ny
end if
c
c htot is the total water depth, and tau is the wind stress amplitude.
c rnonl can be set to zero to remove nonlinear terms (be careful!)
c
if(nov.eq.1) then
write (6,1903)
1903 format (1x,//,1x,' HTOT is the water column depth (m),',//,
1 1x,' TAU is the wind stress amplitude (Pa),',//,
2 1x,' RNONL multiplies the nonlinear term in the vort. eqn.',//)
end if
c
htot = 500.
tau = 0.1
rnonl = 1.
c
write (6,436)
436 format (1x,
1 ' Enter HTOT (def=500.), TAU (def=0.1), RNONL (def=1.0)')
read (5,*) htot, tau, rnonl
c
c Make tau into a friction velocity squared so that density does not appear.
c
r0 = 1025.
tau = tau/r0 ! divide tau by nominal seawater density
c
g = 9.8
c
c Set the horizontal and temporal resolution.
c
if(nov.eq.1) then
write (6,1904)
1904 format
1 (1x,//,1x,' DX is the horizontal grid resolution (km),',//,
2 1x,' DT is the time step (sec)',//)
end if
c
dx = 20.
dy = 20.
dt = 2.e4
c
write (6,778)
778 format (1x,' Enter DX (def=20.), DT (def=2.e4.)')
read (5,*) dx, dt
c
dx = dx*1000.
dy = dx
sdx = 1./(2.*dx)
sdy = 1./(2.*dy)
c
c
c Set the horizontal diffusivity and decay time scale. The latter is used
c to calcualte a linear drag coefficient, as in Stommel 1948 and Veronis 1966.
c To set this linear decay process to zero, set decayt < 0.
c
if(nov.eq.1) then
write (6,1906)
1906 format
1 (1x,//,1x,' DIFFUS is the diffusion coefficient (m*m/s),',//,
2 1x,' DECAYT is the decay time in a linear drag form (days)',//)
end if
c
diff = 0.
decayt = 20.
difmax = 0.25*(dx*dy)/dt
write (6,442)
442 format (1x,' Enter DECAYT (def=20. 1/days) ')
read (5,*) decayt
c
decay = 0.
if(decayt.gt.0.1) decay = 1./(decayt*8.64e4)
c
time = 0.
c
pi = 4.*atan(1.)
tpi = 2.*pi
di = tpi/float(nx - 1)
dj = tpi/float(ny - 1)
c
rlat = 30.
ibeta = 1
c
c Set the latitude and beta (0, 1, or -1).
c
if(nov.eq.1) then
write (6,1907)
1907 format (1x,//,1x,' RLAT is the latitude (deg),',//,
1 1x,' IBETA multiplies beta (=0 to give beta=0)',//)
end if
c
write (6,4467)
4467 format (1x,' Enter RLAT (def=30.), IBETA (def=1) ')
read (5,*) rlat, ibeta
c
ieq = 0
if(abs(rlat).lt.0.01) ieq = 1
f0 = 2.*7.292e-5*sin(rlat*tpi/360.)
beta = float(ibeta)*2.*7.292e-5*cos(rlat*tpi/360.)/6370.e3
c
c Compute the x and y coordinates (km).
c
do 9 j=1,nx
x(j) = float(j-ncx)*dx/1000.
9 continue
do 10 k=1,ny
y(k) = float(k-ncy)*dy/1000.
f(k) = f0 + y(k)*1000.*beta
10 continue
c
c Set the wind stress fields here. Use sin(y) to have downward pumping
c over the entire basin, with a maximum stress on the north and south
c boundaries (as in Stommel, 1948). Use cos(y) to have a max wind stress
c at y = 0. To add a 'doldrums' to the sin(y) form, set dolamp to about 0.5.
c
dolamp = 0.0
doltau = dolamp*tau
ndol = ny - int(0.8*ny)
do 8010 j=1,nx
do 8010 k=1,ny
taux(j,k) = tau*sin(y(k)*3.1415/(2.*abs(y(1))))
taux(j,k) = taux(j,k) +
c doltau*exp( -(abs((y(ndol)-y(k)))/(100.))**2 )
tauy(j,k) = 0.
c
c try inserting regions of zero stress curl by setting squez = 2, say
c
squez = 1.0
ysq = squez*y(k)*3.1415/(2.*abs(y(1)))
taux(j,k) = tau*sin(ysq)
if(ysq.gt.(3.14159/2.)) taux(j,k) = tau
if(ysq.lt.(-3.14159/2.)) taux(j,k) = -tau
if(j.eq.ncx) write (6,9932) y(k), ysq, taux(j,k)
9932 format (1x,' y, ysq, taux', 5e12.4)
c
8010 continue
c
c Write the wind stress to a file.
c
do 8012 k=1,ny
write (16,8013) y(k), taux(ncx,k), tauy(ncx,k)
8012 continue
8013 format (10e15.5)
c
c
c If this is a completely cyclic domain then initialize s (it may
c not make sense to have tau set nonzero).
c
c if(cyclic.eq.2) then
c
if(nov.eq.1) then
write (6,1908)
1908 format
1 (1x,//,' SAMP is the amplitude of an initial s blob (m*m/s),',//,
2 1x,' SSIZE is the horizontal size of the blob (km)',//)
end if
c
samp = 0.
ssize = 400.
write (6,3400)
3400 format
c (1x,' Enter SAMP (def= 0. (try 2.e7)), SSIZE (def = 400.)')
read (5,*) samp, ssize
c
do 3401 j=1,nx
do 3401 k=1,ny
rs = (x(j))**2 + (y(k))**2
s(j,k) = samp*exp(-rs/ssize**2)
3401 continue
c end if
c
c Initialize the old s fields.
c
do 16 j=1,nx
do 16 k=1,ny
sold(j,k) = s(j,k)
16 continue
c
c Evaluate the important non-d parameters; sno is the streamfunction
c amplitude (dimensional), ekno is the Ekman number, rno is the Rossby
c number. This scaling holds in the interior, and hence will lead to
c a huge underestimate of the Rossby and Ekman numbers in the western
c boundary layer.
c
rl = float(nx-1)*dx
ekno = decay/f0
blth = 3.1415*(decay/beta)/1000.
rno = tau/(htot*(beta**2)*(rl**3))
sno = tau/(htot*beta)
if(cyclic.eq.2) sno = samp
write (6,2244) blth, ekno, rno
2244 format
1 (1x,///,' The blthick (km), Ekman and Rossby nos are:',
2 /,6x,3f12.4 //)
c
c
c Write some data defining this run to file what.dat for use by later analysis.
c
nnn = 1
xnnx = 0.
write (11,1110) nx, ny, nnn, dt, dx, dy, x(1), y(1), htot,
c f0, beta, tau, decay, sno, ekno, rno
1110 format (3i8, 20e15.5)
c
c ncycle is the number of time steps taken.
c
ncycle = 0
c
c Begin time stepping here.
c
90 continue
c
leap = 0
ncycle = ncycle + 1
c
time = time + dt
dtime = time/(24.*3600.)
c
c Compute the average kinetic energy and enstrophy.
c
ske = 0.
rke = 0.
se = 0.
do 2000 j=jf,jl
do 2000 k=kf,kl
ske = ske + htot*(u(j,k)**2 + v(j,k)**2)/2.
se = se + del2s(j,k)**2
rke = rke + 1.
2000 continue
ske = ske/rke
se = se/rke
c
c Check to see if it's time to write to the screen.
c
if(timt.gt.wrtt.or.ncycle.eq.1) then
if(ncycle.eq.1) write (6,1990)
1990 format (1x,
1' The data listed below are:',/,
2 1x,' time (days), kinetic energy, enstrophy')
write (6,653) dtime, ske, se
653 format (1x,f10.3,2e12.3)
timt = timt - wrtt
end if
timt = timt + dt/8.64e4
c
c Write data to disk for later analysis and plotting.
c
c Write out some data as a time series.
c
if(mod(ncycle,5).eq.1) then
c
c Evaluate the basin-averaged vorticity balance
c
vorta = 0.
curlta = 0.
rsjva = 0.
bva = 0.
vdeca = 0.
npa = 0.
c
do 248 j=2,nxm
do 248 k=2,nym
npa = npa + 1
vorta = vorta + vort(j,k)
curlta = curlta + curlt(j,k)
rjsva = rjsva + rjsv(j,k)
bva = bva + bv(j,k)
vdeca = vdeca + vdecay(j,k)
248 continue
c
do 249 k=2,nym
bva = bva + 0.5*(bv(2,k) + bv(nxm,k))
vdeca = vdeca + 0.5*(vdecay(2,k) + vdecay(nxm,k))
249 continue
c
rnpa = float(npa)
rnpap = rnpa + float(nym-1)
vorta = vorta/rnpa
bva = bva/rnpap
vdeca = vdeca/rnpap
curlta = curlta/rnpa
rjsva = rsjva/rnpa
c
c A time series of many variables
c
nbx = 3
ndl = ndol + 3
c
write (10,2234) dtime, ske, se, vort(ncx,ncy), bv(ncx,ncy),
1 curlt(ncx,ncy),vdecay(ncx,ncy),rjsv(ncx,ncy),
2 vort(nbx,ncy), bv(nbx,ncy),
3 curlt(nbx,ncy),vdecay(nbx,ncy),rjsv(nbx,ncy),
4 vorta, bva, curlta, vdeca, rjsva,
5 vort(ncx,ndl), bv(ncx,ndl),
6 curlt(ncx,ndl),vdecay(ncx,ndl),rjsv(ncx,ndl)
2234 format (1x, 100e10.3)
c
c The streamfunction profile across the central latitude
c
write (14,2214) (s(j,ncy), j=1,nx)
2214 format(120e12.4)
end if
c
c Write to disk for contour plots of some fields. This data is used by a
c Matlab plotting program called basin.m
c
if(timd.gt.wrtd.or.ncycle.eq.1) then
m = 1
c
do 700 j=jf,jl
jp = jpa(j)
jm = jma(j)
do 700 k=kf,kl
kp = kpa(k)
km = kma(k)
c
u(j,k) = -sdy*(s(j,kp) - s(j,km))/htot
v(j,k) = sdx*(s(jp,k) - s(jm,k))/htot
sxx = (1./dx**2)*(s(jp,k) + s(jm,k) - 2.*s(j,k))
syy = (1./dy**2)*(s(j,kp) + s(j,km) - 2.*s(j,k))
del2p(j,k) = f(k)*(sxx + syy)
rv(j,k) = sxx + syy
spd(j,k) = sqrt(u(j,k)**2 + v(j,k)**2)
bc(j,k) = 0.
c
700 continue
c
c Estimate the surface pressure that is consistent with this s field (if
c it were steady). This is a long and complex thing of marginal use as a
c diagnostic only (not needed for the integration).
c
c Evaluate the normal derivative on the east and west sides.
c
do 800 k=1,ny
dsdx1 = (s(2,k) - s(1,k))/dx
dsdx2 = (s(nx,k) - s(nxm,k))/dy
bc(1,k) = f(k)*dsdx1 + taux(1,k)/htot
bc(nx,k) = f(k)*dsdx2 + taux(nx,k)/htot
800 continue
c
c South and north sides.
c
do 801 j=1,nx
bc(j,1) = f(1)*(s(j,2) - s(j,1))/dx + tauy(j,1)/htot
bc(j,ny) = f(ny)*(s(j,ny) - s(j,nym))/dy + tauy(j,ny)/htot
801 continue
c
ibc = 1
idmean = 0
c
call sorbc
1 (del2p,bc,ibc,cyclic,idmean,dx,dy,nx,ny,jma,jpa,kma,kpa,rlp)
c
do 1206 k=1,ny
write (1,1298) (rv(j,k),j=1,nx)
write (2,1298) (v(j,k),j=1,nx)
write (3,1298) (s(j,k),j=1,nx)
write (4,1298) (rlp(j,k),j=1,nx)
1206 continue
1298 format (200e12.4)
timd = timd - wrtd
c
c Write the present time (days) to file looktim.dat
c
write (12,1212) dtime
1212 format (f10.3)
end if
c
timd = timd + dt/8.64e4
c
c Come here to start the second pass of a leap-frog tapezoidal correction
c (this avoids adding to time or ncycle).
c
901 continue
c
c Compute the velocity and the vorticity.
c
do 101 j=jf,jl
jp = jpa(j)
jm = jma(j)
c
do 101 k=kf,kl
kp = kpa(k)
km = kma(k)
c
u(j,k) = -sdy*(s(j,kp) - s(j,km))/htot
v(j,k) = sdx*(s(jp,k) - s(jm,k))/htot
c
sxx = s(jp,k) + s(jm,k) - 2.*s(j,k)
syy = s(j,kp) + s(j,km) - 2.*s(j,k)
del2s(j,k) = sxx/dx**2 + syy/dy**2
101 continue
c
c OK, now you are ready to evaluate the vorticity eqn.
c
c Compute the terms in the vorticity tendency equation.
c
do 201 j=jf,jl
jp = jpa(j)
jm = jma(j)
c
do 201 k=kf,kl
kp = kpa(k)
km = kma(k)
c
c the beta term
c
bv(j,k) = -beta*v(j,k)*htot
c
c the curl of the wind stress
c
curlt(j,k) = sdx*(tauy(jp,k) - tauy(jm,k)) -
1 sdy*(taux(j,kp) - taux(j,km))
c
c the linear decay form of dissipation as if by bottom friction
c
vdecay(j,k) = -decay*del2s(j,k)
c
c the diffusion form of dissipation
c
c vdiff(j,k) =
c 1 (diff/dx**2)*(del2s(jmd(j),k) + del2s(jpd(j),k) - 2.*del2s(j,k))
c 2+ (diff/dy**2)*(del2s(j,kmd(k)) + del2s(j,kpd(k)) - 2.*del2s(j,k))
c
c the Jacobian, or advection term, written in flux form
c
rjsv(j,k) =
1 - ( sdx*(u(jp,k)*del2s(jp,k) - u(jm,k)*del2s(jm,k)) +
2 sdy*(v(j,kp)*del2s(j,kp) - v(j,km)*del2s(j,km)) )
c
c sum the terms
c
vort(j,k) =
1 curlt(j,k) + vdecay(j,k)
2 + bv(j,k) + rnonl*rjsv(j,k)
c + vdiff(j,k) ! diffusion term c'd out c
c
c add to the tendency (this is needed for the trapezoidal correction)
c
del2st(j,k) = del2st(j,k) + vort(j,k)
c
201 continue
c
c Call the elliptic solver sorbc to calculate st, the time rate change of
c the streamfunction.
c
c Set the matrix bc to zero for this since s is specififed as part of the bc.
c
do 930 j=1,nx
do 930 k=1,ny
bc(j,k) = 0.
930 continue
c
c ibc sets the kind of boundary condition applied on solid boundaries;
c ibc = 0 sets s to zero (free slip boundary conditions)
c ibc = 1 sets both s and the normal derivative of s to zero (for no slip
c boundary conditions). the flag cyclic tells the subroutine which boundaries
c are solid.
c
ibc = noslip
idmean = 0
call sorbc
1 (del2st,bc,ibc,cyclic,idmean,dx,dy,nx,ny,jma,jpa,kma,kpa,st)
c
c Step ahead, noting whether this is a straight leap-frog (leap = 0), a
c leap-frog to be followed by a trapezoidal correction (leap = 1), or
c the trapezoidal correction itself (leap = 2). A trapezoidal
c correction is done for ntrap succesive steps every nncyc steps in order
c to supress time splitting (try ntrap < 0 to see time splitting).
c
nncyc = 15
ntrap = 3
c
if(mod(ncycle,nncyc).lt.ntrap) leap = leap + 1
c
tstep = 2.*dt
if(leap.eq.2) tstep = 0.5*dt
c
do 300 j=1,nx
do 300 k=1,ny
snew(j,k) = sold(j,k) + tstep*st(j,k)
300 continue
c
if(ncycle.eq.1) go to 310
c
c Shuffle time levels around, as required for the kind of time step.
c
c The following is for a straight leap-frog, or the trapezoidal corrector.
c
if(leap.eq.0) then
do 3300 j=1,nx
do 3300 k=1,ny
sold(j,k) = s(j,k)
s(j,k) = snew(j,k)
del2st(j,k) = 0.
3300 continue
end if
c
c This is for a leap-frog that will be followed by a trapezoidal correction.
c
if(leap.eq.1) then
do 3301 j=1,nx
do 3301 k=1,ny
sold(j,k) = s(j,k)
s(j,k) = snew(j,k)
3301 continue
c
c The first stage of a leap-trap time step is done; go back up and re-do
c the calculations for the trapezoidal correction.
c
go to 901
end if
c
if(leap.eq.2) then
do 3500 j=1,nx
do 3500 k=1,ny
s(j,k) = snew(j,k)
del2st(j,k) = 0.
3500 continue
end if
c
c
c The time step is finished. Check the time, and then loop back
c up to start another time step.
c
if(dtime.gt.runtim) go to 9999
c
go to 90
c
310 continue
c
c Come here once at the start of the integration to do a forward time step
c to get started.
c
do 400 j=1,nx
do 400 k=1,ny
sold(j,k) = s(j,k)
s(j,k) = snew(j,k)
400 continue
c
go to 90
c
9999 continue
c
stop
end
c
subroutine
1 sorbc(f,bc,ibc,icyclic,idmean,dx,dy,nx,ny,jm,jp,km,kp,u)
c
c This subroutine solves the Poisson equation, DelSq(u) = f, by simultaneous
c overrelaxation. Taken from Numerical Recipes, modified heavily, and
c documented by Jim Price.
c
c f is input, and is the source term,
c bc is input and is the normal derivative used in evaluating the
c boundary values (input as a matrix for convenience though only the
c boundaries are used),
c ibc = 0 for u = specified bcs, in which case bc is not used
c = 1 for du/dy = bc(y=0,x), etc., bcs.
c icyclic = 0 for solid boundaries,
c = 1 for a channel (cyclic in x direction only),
c = 2 for a fully cyclic domain.
c idmean = 1 to demean the field f. This makes sense if the bcs are
c specified by normal derivatives all the way around.
c dx and dy are input and are the grid intervals.
c ny and nx are the array dimensions. Note that the Poisson equation is
c solved over the interior points only. The boundaries are dealt with
c separately.
c The integer arrays jm, jp, km, kp are input to specify indices needed to
c evaluate the Laplacian (needed to specify bcs; i.e, symmetric or not).
c u is a first guess at the answer on input, and is returned as the
c answer on output. If boundary conditions are u specified on the boundaries,
c then u has those bcs upon input.
c
c
dimension f(nx,ny),bc(nx,ny),u(nx,ny),jm(nx),jp(nx),
1 km(ny),kp(ny)
real omega
c
c Compute rjac, an overrelaxation parameter used to aid convergence.
c
rn = sqrt(float(nx**2) + float(ny**2))
rjac = cos(3.1415/rn)
c
c These maxits and epsilon values are rather arbitrary.
c
maxits = 1000
eps = 1.e-5
anormf = 0.
c
a = 1./(dx**2)
b = 1./(dy**2)
c = -2.*(dx**2 + dy**2)/((dx**2)*(dy**2))
c
omega = 1.0
c
jf = 2
jl = nx - 1
if(icyclic.eq.1.or.icyclic.eq.2) then
jf = 1
jl = nx
end if
c
kfm = 1
kl = ny - 1
if(icyclic.eq.2) then
kfm = 0
kl = ny
end if
kf = kfm + 1
c
do 12 j=jf,jl
do 12 k=kf,kl
anormf = anormf + abs(f(j,k))
12 continue
c
do 16 n=1,maxits
anorm = 0.
jsw = 1
do 15 ipass=1,2
lsw = jsw
c
do 14 j=jf,jl
do 13 k=lsw+kfm,kl,2
c
resid = b*(u(jp(j),k) + u(jm(j),k)) + a*(u(j,kp(k)) + u(j,km(k)))
1 + c*u(j,k) - f(j,k)
anorm = anorm + abs(resid)
u(j,k) = u(j,k) - omega*resid/c
c
13 continue
lsw = 3 - lsw
14 continue
jsw = 3 - jsw
c
if(n.eq.1.and.ipass.eq.1) then
omega = 1.0/(1.0 - 0.5*rjac**2)
anorm1 = anorm
else
omega = 1.0/(1.0 - 0.25*rjac**2*omega)
endif
c
c Apply normal derivative bcs, if ibc = 1. Otherwise, assume
c that u has already been set along the boundaries.
c
if(ibc.eq.1) then
c
c March along the north and south sides, noting whether the bc looks forward
c or backward to the boundary. do this only if not a symmetric domain.
c
if(icyclic.eq.2) go to 299
do 20 j=1,nx
u(j,1) = u(j,2) - bc(j,1)*dy
u(j,ny) = u(j,ny-1) + bc(j,ny)*dy
20 continue
299 continue
c
c March along the east and west sides if this is not a channel.
c
if(icyclic.eq.1.or.icyclic.eq.2) go to 219
do 21 k=2,ny-1
u(1,k) = u(2,k) - bc(1,k)*dx
u(nx,k) = u(nx-1,k) + bc(nx,k)*dx
21 continue
219 continue
c
c
end if
c
c End of the boundary conditions.
c
15 continue
c
c Check for convergence.
c
epsnorm = eps*anormf
if(anorm.lt.epsnorm) go to 99
c
16 continue
c
write (6,1) n, rjac, omega, anorm1, anorm, epsnorm
1 format (1x,///, ' maxits passed in sorbc ', i6, 5e9.3)
c
99 continue
c
c Demean u if the bcs are zero normal derivative all the way around, or if
c the domain is fully symmetric.
c
if(idmean.eq.1) then
do 50 j=1,nx
do 50 k=1,ny
um = um + u(j,k)
50 continue
um = um/float(nx*ny)
do 51 j=1,nx
do 51 k=1,ny
u(j,k) = u(j,k) - um
51 continue
end if
c
return
end