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