c
      program geoadj_2d
c
c
c   Compute geostrophic adjustment by simple, numerical methods.    
c   This is public domain for all educational purposes.
c   Jim Price, WHOI, 508-289-2436
c   jprice@whoi.edu
c
c   All subroutines are included here.
c
c   This program may be used to compute the response of an 
c   ocean basin or channel to a prescribed mass field. 
c   The ocean model is layered in the vertical, and uses a 
c   staggered grid (a C grid) in the horizontal. It uses 
c   an energy-conserving finite differencing scheme. 
c   Time stepping is a leap-frog with an occasional
c   trapezoidal correction in order to avoid time-splitting.
c
c   The variables are laid out on a C grid as follows:
c
c      h u | h u | h u
c      v s | v s | v s
c      ---------------
c      h u | h u | h u 
c      v s | v s | v s 
c      ---------------
c      h u | h u | h u
c      v s | v s | v s 
c  
c   h, u, v have their obvious meanings; pressure is defined at the
c   same point as is h, and s is the stream function. Wind stress is
c   defined at the corresponding u,v points. A tracer, t, is defined at
c   the h points, as is the Coriolis parameter, f. The x,y coordinates
c   are presumed to be centered on the h points.
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   The top most row (h and u only) and the left most column (h and v only) 
c   are junk (not meanignful) if the BCs are for solid boundaries. 
c                                                                             
c   u and v are the east and north transports in their respective layers. To
c   get the speeds (m/sec) you must divide by the appropriately averaged
c   layer thickness.  
c
c   The variables are stored in the u,v,h,t, arrays that have dimensions
c   u(nx,ny,nx,3), where the last index refers to time level or otherwise
c   as follows:
c     u(...1) is the present time 
c     u(...2) is the past time
c     u(...3) is the tendency for this variable (i.e., du/dt).
c
c
c.................................. Data Allocation ..........................
c
      parameter (nx=251, ny=251, nz=1, nxp=nx+1, nyp=ny+1,
     1    nf=1, nf2=2*nf)
c
      dimension u(nx,ny,nz,3),v(nx,ny,nz,3),
     1  h(nx,ny,nz,3),t(nx,ny,nz,3)
c     
      dimension h0(nz),delr0(nz),gprime(nz),x(nx),y(ny),
     1 p(nx,ny,nz),wa(nz),
     2 f(ny),taux(nx,ny),tauy(nx,ny),work7(nx,ny),work8(nx,ny),
     3 kpa(ny),jpa(nx),kma(ny),jma(nx),work5(nx,ny),
     4 work6(nx,ny),work4(nx,ny),ctu(nx,ny,nz),ctv(nx,ny,nz),
     5 kmd(ny),kpd(ny),jmd(nx),jpd(nx),eta(nx,ny,nz)
c     
      dimension ui(nx,ny,nz,3),vi(nx,ny,nz,3),uuhath(nx,ny,nz),
     1 vvhath(nx,ny,nz),uvhats(nx,ny,nz),ud24(nx,ny,nz),
     2 vd24(nx,ny,nz),td24(nx,ny,nz),work24(nx,ny,nz)
c
c  The arrays below are required solely for the rigid lid.
c
      dimension ua(nx,ny),va(nx,ny),del2st(nx,ny),st(nx,ny),
     1 pr(nxp,nyp),d2pr(nxp,nyp),prm(nx,ny),
     2 jpap(nxp),jmap(nxp),kmap(nyp),kpap(nyp)     
c
c  The arrays below store float data.
c
      dimension xf(nf),yf(nf),xyfk(nf2)
c      
      integer frslip, cyclic, rgrav, rlid 
      character*24 ifile 
c
      ncx = nx/2 + 1 
      nxm = nx - 1
      ncy = ny/2 + 1
      nym = ny - 1
      if(nym.eq.0) nym = 1
      nzm = nz - 1
      if(nz.eq.1) nzm = 1
c
      g = 9.8
      r0 = 1025.      
      tpi = 2.*4.*atan(1.)
      time = 0. 
      wworks = 0.
      fworks = 0. 
      timd = 0.
      timm = 0.
      timt = 0.
c 
c .............................. Open Files ....................                       
c      
      open (unit=1,file='ga_eta.dat',status='unknown')
      open (unit=2,file='ga_u.dat',status='unknown')
      open (unit=3,file='ga_v.dat',status='unknown')
	open (unit=4,file='ga_t.dat',status='unknown') 
      open (unit=32,file='ga_time.dat',status='unknown')	    
      open (unit=35,file='ga_float.dat',status='unknown')     
      open (unit=30,file='ga_what.dat',status='unknown')
c 
c ............................. Initialize Boundary Conditions ...............
c
c  set a flag that determines the boundary conditions. 
c
      frslip = 1
      cyclic = 0 
c
      ibc = 3    ! for solid walls all the way around
	ibc = 2    ! for open boundaries on the north and south, and solid on east 
c                    and west; used for an equatorial beta-plane
	ibc = 1    !   for radiation on all sides   (the default)
c
c      write (6,9922) frslip, cyclic
c 9922 format
c     1 (1x,/,'  define the BCs: frslip, cyclic; = '3i5)
c      read (5,*) frslip, cyclic
c  
c ............................ Initialize the Pressure Model .................
c
c  The pressure model is set by the following flags:
c   
c  rgrav = 0 for a free surface model,
c  rgrav = 1 for a reduced gravity model (the default),
c  rlid = 1 if you want to apply a rigid lid pressure correction.
c
c  rnonl can be set to zero to remove nonlinear terms from the momentum and 
c  tracer equations (be careful with this !).
c     
      rgrav = 0
      rlid = 0
      rnonl = 1.
c     
c 9135 continue
c      write (6,9923) rgrav, rlid, rnonl 
c 9923 format (1x,/, 
c     1 '  define pressure: rgrav, rlid, rnonl; = '2i5,f5.1)
c      read (5,*) rgrav, rlid, rnonl 
c
c      if(nz.eq.1.and.rlid.eq.1) then
c      write (6,9134) 
c 9134 format (1x,///,
c     1 '   It does not make sense to have nz and rlid = 1. Try again.')
c      go to 9135
c      end if      
c 
c ....................... Initialize the Density Profile ..................
c
c  This next section sets the initial density 'profile' by defining the mean 
c  layer thicknesses and density differences as follows: 
c
c  htot is the total water column depth,
c  rtcl is the density difference across the thermocline,
c  delradd is a density that is added to the upper layer (to make
c   a free surface model with a realistic air-sea density difference),
c  h0 are the nominal (initial) layer depths,
c  delr0 are the fixed density differences across the layer interfaces
c   (which layer depends upon the kind of pressure model so pay attention).
c 
      htot = 4000.
      hsl = 500.
      if(nz.eq.1) htot = hsl
      rtcl = 2.0 
      delradd = 0.
c                 
      write (6,9924) htot
 9924 format (1x,/, 
     1 ' Enter layer thickness[m] = ', f8.0)
      read (5,*) hsl   
      if(nz.eq.1) htot = hsl
      write (6,9925) rtcl
 9925 format (1x,/,
     1 ' Enter delta rho [kg/m**3]; = ', f8.2)
      read (5,*) rtcl
c            
c  Set the density profile for a reduced gravity model. In this case the density
c  differences are associated with the bottom of the layer.
c      
      if(rgrav.eq.1.and.nz.gt.1) then                   
      hs = 0.
      do 780 m=1,nzm
      h0(m) = hsl
      hs = hs + h0(m)
      delr0(m) = rtcl/float(nzm)
  780 continue     
      h0(nz) = htot - hs
      delr0(nz) = 0.
      end if
c
      if(nz.eq.1) then
      h0(1) = hsl
      delr0(1) = rtcl
      end if
c
c  Check whether this is a free surface model. In that case, the density
c  differences are associated with the top of the layer and delradd is added to
c  the topmost layer (the sea surface).
c      
      if(rgrav.eq.0) then 
      delr0(1) = rtcl/(float(nz)) + delradd
      if(nz.gt.1) then
      do 781 m=2,nz
      delr0(m) = rtcl/float(nz)
  781 continue 
      end if
      end if  
c 
c  Compute the reduced gravity for each layer (used in pressure calculations).
c
      do 785 m=1,nz
      gprime(m) = g*delr0(m)/r0
c
c      write (6,8822) m,delr0(m),gprime(m),h0(m)
c 8822 format (1x,'   m, delr, gp, h0', i4,3e12.4)
c             
  785 continue          
c
c ................... Initialize the Time and Space Intervals .................
c 
c  Set the time step, dt,  and the horizontal resolution, dx, dy.
c
 1177 continue
c 
      dt = 600.
      dx = 15.e3 
      dy = 15.e3    
c
      dxk = dx/1000.
      dyk = dy/1000.
      write (6,442) dt, dxk, dyk
  442 format (1x,/,
     1 ' Enter resolution: dt[s], dx and dy[km]; = ',f8.0,2f7.1)
      read (5,*) dt, dxk, dyk 
c 
      twodt = 2.*dt
      halfdt = 0.5*dt      
c
      dx = dxk*1000.
      dy = dyk*1000. 
      sdx = 1./dx
      sdy = 1./dy
c
c  Estimate the CFL limit on dt from the long wave phase speed only.
c
      cph = sqrt(g*rtcl*h0(1)/r0) + 0.1
      if(rgrav.eq.0) cph = sqrt(g*(delradd + rtcl)*htot/r0) + 0.01 
c      
      dtcfl = (dx/2.)/cph
      if(dtcfl.lt.dt) then
      write (6,3378) dt, dtcfl
 3378 format (1x,/,'  This dt exceeds the CFL limit:', 2f8.1,/)
      go to 1177
      end if      
c      
c ....................... Initialize Spatially Dependent Fields ................
c       
c  Set the latitude and beta (on or off by way of the flag ibeta). The choice
c  rlat = 0. should be OK, and note that beta could be made negative by
c  setting ibeta = -1
c 
      rlat = 20.
      ibeta = 1 
	Rearth = 6370.e3
c
      write (6,4467) rlat, ibeta 
 4467 format (1x,/,
     c  ' Enter rotation: rlat[deg], ibeta (0, 1, or 2) ', f6.0, i3)
      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 = 2.*7.292e-5*cos(rlat*tpi/360.)/Rearth 
	beta2 = -2.*7.292e-5*sin(rlat*tpi/360.)/(Rearth*Rearth)
	if(ibeta.ne.2) beta2 = 0. 
c 
c  Define x and y, the east and north coordinates in meters, and the variable
c  Coriolis parameter, f.
c                     
      do 10 j=1,nx
      x(j) = float(j-ncx)*dx
c      
      do 10 k=1,ny
      y(k) = float(k-ncy)*dy
      if(ibeta.eq.0) f(k) = f0
	if(ibeta.eq.1) f(k) = f0 + y(k)*beta 
	if(ibeta.eq.2) f(k) = f0 + y(k)*beta  + y(k)*y(k)*beta2/2.
   10 continue
c
c  Estimate the nondispersive baroclinic long and Rossby wave phase speeds.
c     
      if(abs(f0).gt.1.e-9) then
      rcph = -(cph**2)*beta/f0**2
      rcphkd = rcph*8.64e4/1000.
      rdef = (cph/f0)/1000. 
      write (6,3045) cph, rcph, rcphkd, rdef 
 3045 format (1x,' C, Ro wave, Rokmperday and Rd are ',3f8.3,f8.1)
      end if
c       
c  Define the field of wind stress. The form can be a single gyre, as in 
c  Stommel 1948 in which case use cos(y), or a double gyre. In either case,
c  the amplitude and sign are set by tau, which should be positive for an 
c  anticyclonic curl in the single gyre case. tau is immediately divided by
c  density to give the friction velocity. tramp is the time over which the 
c  wind stress is brought up to full strength.
c 
      tau = 0.0
      tramp = 0.01
c
c      write(6,4001) tau, tramp
c 4001 format (1x,/,
c     1 '  define wind stress: tau[Pa], tramp[days]; = ',2f7.2)
c      read (5,*) tau, tramp
c      tau = tau/r0 
c                  
      do 1080 j=1,nx
      do 1080 k=1,ny
c
      tauy(j,k) = 0.      
      yphas = 3.1415*(y(k)-y(1))/(y(ny) - y(1))
      taux(j,k) = -tau*cos(yphas)     
c
c  The form sin(y) below gives two gyres.     
c      taux(j,k) = 2.*tau*(sin(yphas) - 0.5) 
c
 1080 continue 
c 

c  Initialize the layer thicknesses, and as an option include an initial 
c  thickness anomaly (eddy) in the h field if this is a reduced gravity model.
c
      delh = 20.          !  m
      hscale = 200.       ! radius,  km
c
      write (6,4007) delh, hscale
 4007 format (1x,/,' Enter h anomaly: delh[m], hscale[km]; = ',2f7.0)
      read (5,*) delh, hscale                             
c
      do 1081 j=1,nx
      do 1081 k=1,ny
c      
      rr = sqrt( (x(j)/1000.)**2 + (y(k)/1000.)**2 )
      r9 = rr - hscale
	rtaper = 25.
      eddyh = 
     c  delh*0.5*(1. - tanh(r9/rtaper)) + 
     c        1.0*(1. - tanh(r9/(100.*rtaper)))  !  to make a definite max
c
      if(rr.gt.2*hscale) eddyh = 0.
c            
      hs = 0.
      do 1050 m=1,nzm
      h(j,k,m,1) = h0(m) + eddyh
      hs = hs + h(j,k,m,1)
      eta(j,k,m) = eddyh
 1050 continue 
c
c  Now set the lower layer thickness to make the total thickness equal htot.      
c             
      if(nz.gt.1) h(j,k,nz,1) = h(j,k,nz,1) - (hs - htot)
c
c  Define an initial tracer field (note that t is tracer times h). t1 is
c  the horizontal scale of the tracer patch.
c      
      tl = 100.
      rr = (x(j)/1000. + 500.)**2  
      t(j,k,1,1) = h(j,k,1,1)
      if(sqrt(rr).gt.(2.*tl)) t(j,k,1,1) = 0.  
c  
 1081 continue
c
c  Specify the initial position of the floats. 
c     
      if(nf.gt.1) then
      nx3 = nx/3 
      nx32 = 2*nx3
      ny3 = ny/3  
      ny32 = 2*ny3
      xf(1) = x(nx3)
      xf(2) = x(nx32)
      xf(3) = x(nx3)
      xf(4) = x(nx32)
      yf(1) = y(ny3)
      yf(2) = y(ny3)
      yf(3) = y(ny32)
      yf(4) = y(ny32)   
      end if
c  
c .......................  Some Miscellaneous Initialization  ..............
c
c  Define the subgridscale diffusivity, and whether it is to be applied by 
c  a Laplacian or a biharmonic operator. Also enter the e-folding time for
c  a linear bottom drag. 
c   
      diff = 3.0e2
      ibhm = 0
      efold = 9999.
c      
c      write (6,1188) diff, ibhm, efold
c 1188 format (1x,/,
c     1 '  Dissipation: diff[m*m/s], biharm, efold ',e7.1,i3, e9.2)
c      read (5,*) diff, ibhm, efold
      cbd = 1./(efold*8.64e4)
c
c Initialize the old fields.
c
      do 16 j=1,nx
      do 16 k=1,ny
      do 16 m=1,nz
      u(j,k,m,1) = 0.
      v(j,k,m,1) = 0.
      u(j,k,m,2) = u(j,k,m,1)
      v(j,k,m,2) = v(j,k,m,1)
      t(j,k,m,2) = t(j,k,m,1)
      h(j,k,m,2) = h(j,k,m,1)
   16 continue  
c
c  Come here to begin integrating, whether you have done a restart or not.
c 
 7100 continue
c 
c  wa are weights used to apply the wind stress to the upper layer only.
c          
      do 886 m=1,nz
      wa(m) = 0. 
  886 continue    
      wa(1) = 1.0
c
c  Set up some arrays used to specify indices when horizontal differences
c  are taken. These are set up differently for solid or cyclic boundaries.
c  kma, jpa, etc are for the advection operations;
c  kmd, jpd, etc are for the Laplacian.
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 In the case of a channel, u(1,k) = u(nx,k), and so:    
c
      if(cyclic.eq.1) 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 For a completely symmetric domain, v(j,1) = v(j,ny), and so:
c
      if(cyclic.eq.2) then
      kpa(ny) = 1 
      kma(1) = ny 
      kpd(ny) = 1
      kmd(1) = ny
      end if
c
      if(rnonl.ne.1.0) write (6,8234) rnonl
 8234 format (1x,///,
     1  '   did you know that rnonl is not equal 1.0, =', f5.2)      
c 
c ........................... Define Some Control Variables ...............
c 
c  Define some times used to control the integration and the ouput
c  to disk files.
c
c  runto is the time in days at which to stop the integration,
c  wrtd is the time interval (days) between writing to the Matlab disk files, 
c  wrtm is the time between writing to a Matlab move file, and, 
c  wrtt is time between writing to the screen. 
c
      time99 = 0.
      runto = 160.
      wrtd = 0.5
      wrtm = 1000.
      wrtt = 0.25
c           
      write (6,777)  
  777 format (1x,/,' Enter the run and output times: runto (160.),',
     1 ' wrtdisk (0.5)')
      read (5,*) runto, wrtd 
c
c  ncycle is the number of time steps taken in this integration.
c
      ncycle = 0
	tim99 = 0
c
c .............................. Start a Full Time Step ..........................         
c
   90 continue  
c                   
c  Reset the flag leap to zero.
c
      leap = 0
      ncycle = ncycle + 1
c                                            
      inow = 0
      if(ncycle.eq.1.and.irest.eq.0) inow = 1      
c
      time = time + dt
      dtime = time/(24.*3600.) 
c
c  Compute the time-dependent wind stress amplitude. The amplitude can be 
c  ramped up to avoid inertial oscillations, but this doesn't seem to be
c  important for most cases.
c
      wtfac = dtime/(tramp + 0.05)
      if(wtfac.gt.1.) wtfac = 1.
c                                                                      
c ...................... Diagnostics and Data Storage ........................
c
c  Compute KE and PE as a diagnsotic of the energy budget.
c                                            
      call kepe(h,ui,vi,eta,h0,delr0,nx,ny,nz,
     1            rgrav,cyclic,rke1,rketc,rpe,ah1,ah2)
      enet = rke1 + rketc + rpe 
c
c  Compute some tracer statistics.
c    
      do 9921 j=1,nx
      do 9921 k=1,ny
 9921 work4(j,k) = t(j,k,1,1)/h(j,k,1,1)
c
      call tracer(work4,nx,ny,ta,tsa,dtsa)      
c
      if(timt.gt.wrtt.or.inow.eq.1) then
      nxw = nx/6 + 1 
      nyw = 2*(ny/3) + 1
      u9 = u(nxw,nyw,1,1)
      v9 = v(nxw,nyw,1,1)
      rn9 = eta(nxw,nyw,1)
	if(ncycle.eq.1) then
	write (6, 2254)
 2254 format (1x, 
     c  '  these data are: time (days), eta, u, v, total energy')
	end if
      write (6,653) dtime, rn9, u9, v9, enet  
  653 format (1x,f7.2,5e12.3)
c 
c 
      if(ncycle.ne.1) timt = timt - wrtt
      end if
      timt = timt + dt/8.64e4
c 
c  Write data to disk files for later analysis and plotting.
c  (the first part makes a time series that could get pretty long)
c
c  Write some basic data defining the case to the file ga2d30.dat.
c
      if(ncycle.eq.1) then
      write (30,2020) nx,ny,nz,nf,dx,dy,x(1),x(nx),y(1),y(ny),
     1  f0,beta,tau,rtcl,htot
 2020 format(4i5,12e20.3)
      end if
c
c
      itimep = 1
      nq = mod(ncycle,3)
      if(itimep.eq.1.and.nq.eq.1) then
c
c  write out some time series data from a few points
c
      nxw = int(nx/2) 
      nyw = int(ny/2 + ny/10)
      u9 = ui(nxw,nyw,1,1)
      v9 = vi(nxw,nyw,1,1)
      h9 = h(nxw,nyw,1,1) - h0(1)
      t9 = eta(nxw,nyw,1)
c       
      nyw = int(ny/2 + 2*ny/10) 
      u10 = ui(nxw,nyw,1,1)    
      v10 = vi(nxw,nyw,1,1)
      t10 = eta(nxw,nyw,1)
c
      nyw = int(ny/2 + 3*ny/10) 
      u11 = ui(nxw,nyw,1,1)    
      v11 = vi(nxw,nyw,1,1)
      t11 = eta(nxw,nyw,1)
c
      nyw = int(ny/2 + 4*ny/10) 
      u12 = ui(nxw,nyw,1,1)    
      v12 = vi(nxw,nyw,1,1)
      t12 = eta(nxw,nyw,1)
c
c      write (10,2234) dtime, u9, v9, t9, t9, rke1, rketc, rpe, enet,
c     c  u10, v10, t10, u11, v11, t11, u12, v12, t12  
 2234 format (1x, 25e15.5)
c
      mm = 0
      do 2235 m=1,nz 
      u9 = ui(nxw,nyw,m,1)
      v9 = vi(nxw,nyw,m,1)
      h9 = h(nxw,nyw,m,1) - h0(m)
      t9 = eta(nxw,nyw,m)
      mm = mm + 1
      work7(mm,1) = u9
      mm = mm + 1
      work7(mm,1) = v9
      mm = mm + 1
      work7(mm,1) = t9
 2235 continue
c      write (40,2234) dtime, (work7(n,1),n=1,mm) 

      end if
c   
c  Write out some data to disk files for general plotting.
c  
      if (dtime.gt.5.0) then ! reduce the output frequency after 5 days
	wrtd = 5.0
	end if 
c	 
c  Write out some matrices that are later read in for Matlab plotting.
c
      time99 = time99 + dt/8.64e4
      if (time99.gt.wrtd.or.ncycle.eq.1) then
      time99 = time99 - wrtd
c 
      write (32,2121) dtime
 2121 format (f15.3)            
c
      do 1206 k=1,ny
      write (1,1298) (eta(j,k,1),j=1,nx)
      write (2,1298) (ui(j,k,1,1),j=1,nx)
      write (3,1298) (vi(j,k,1,1),j=1,nx)       
      write (4,1298) (t(j,k,1,1),j=1,nx) 
c	write (9,1298) (prm(j,k,1,1),j=1,nx)
 1206 continue
 1298 format (900e18.5)
      end if

c 
c      if(timd.gt.wrtd.or.inow.eq.1) then
c           
c      
c  Load the relative vorticity into work7 and work5, then solve for
c  the  streamfunction. The Laplacian is computed over the entire domain,
c  however, it is not used on the solid boundaries where s = 0. 
c 
      do 1301 j=1,nx 
      jp = jpa(j)
      do 1301 k=1,ny
      km = kma(k) 
      work5(j,k) = sdx*(v(jp,k,1,1) - v(j,k,1,1)) - 
     1  sdy*(u(j,k,1,1) - u(j,km,1,1))
      work7(j,k) = sdx*(v(jp,k,nz,1) - v(j,k,nz,1)) - 
     1  sdy*(u(j,k,nz,1) - u(j,km,nz,1))
c
c  To check on a geo adjustment problem output u and v.
c       
      work6(j,k) = ui(j,k,1,1)
      work8(j,k) = vi(j,k,1,1)                     
c
 1301 continue 
c
c      call sorpsi(1,work5,cyclic,dx,dy,nx,ny,jma,jpa,kma,kpa,work6) 
c      call sorpsi(1,work7,cyclic,dx,dy,nx,ny,jma,jpa,kma,kpa,work8)
c
c             
c....................... Come Here to Do a Trapezoidal Correction ............
c
c  Come here to start the second pass of a leap-frog tapezoidal correction
c  to avoid adding to the time or to ncycle.
c
  901 continue 
c
c ............................ Hydrostatic Pressure ..........................
c
c  Compute eta and hydrostatic pressure, just how depending upon the kind 
c  of pressure model.
c 
c  The following is for a reduced gravity model, in which case
c  eta is defined to be at the lower boundary of a layer. 
c     
      if(rgrav.eq.1) then
      do 370 j=1,nx 
      do 370 k=1,ny 
c
c  Compute eta for a reduced gravity model. Start from the surface, and
c  integrate downwards. 
c       
      etas = 0.
      do 371 m=1,nz
      etas = etas - (h(j,k,m,1) - h0(m)) 
      eta(j,k,m) = etas 
  371 continue
c
c  Compute the pressure for a reduced gravity model. Note that this 
c  starts from the bottom, where pressure is assumed to vanish in
c  a deep abyssal layer, and integrates upward.
c 
      ps = 0.
      do 372 m=1,nz
      mm = nz + 1 - m
      ps = ps - gprime(mm)*eta(j,k,mm)
      p(j,k,mm) = ps 
c
c      if(j.eq.ncx.and.k.eq.ncy) then
c      if(mod(ncycle,5).eq.1) then
c      write (6,3546) ncycle,j,k,m,mm,h(j,k,m,1),eta(j,k,m),p(j,k,m)
c 3546 format (1x, ' h,eta,p', 5i4,3e12.5)
c      end if
c      end if
c      
  372 continue
c       
  370 continue
      end if      
c
c  The free surface versions of eta and pressure are computed below.
c  Note that the direction of itegration and the identification of
c  eta with the layers is different from that above. 
c                                  
      if(rgrav.eq.0) then
      do 379 j=1,nx 
      do 379 k=1,ny 
c      
      etas = 0.
      do 377 m=1,nz
      mm = nz + 1 - m
      etas = etas + (h(j,k,mm,1) - h0(mm)) 
      eta(j,k,mm) = etas 
  377 continue
c
      ps = 0.
      do 378 m=1,nz
      ps = ps + gprime(m)*eta(j,k,m)
      p(j,k,m) = ps 
  378 continue 
c
  379 continue
      end if               
c
c  The end of the pressure calculations.
c
c ........................... Evaluate the Budget Equations ..................
c
c  Compute some fields that will need to be differentiated to 
c  estimate horizontal advection terms in the budget equations.
c 
c  Compute the intensive velocities, ui, vi from the mass transports.
c
      do 5680 j=1,nx
      jp = jpa(j)
      do 5680 k=1,ny
      km = kma(k)
      do 5680 m=1,nz
      do 5680 n=1,2
c     
      ui(j,k,m,n) = 2.*u(j,k,m,n)/(h(jp,k,m,n) + h(j,k,m,n))
      vi(j,k,m,n) = 2.*v(j,k,m,n)/(h(j,k,m,n) + h(j,km,m,n))
c      
 5680 continue
c
      do 100 j=1,nx
      jp = jpa(j)
      jm = jma(j)
c      
      do 100 k=1,ny
      kp = kpa(k)
      km = kma(k) 
c
      do 100 m=1,nz
c
c  Compute the tracer transports at the velocity points. 
c
      ctu(j,k,m) = 0.5*(t(j,k,m,1) + t(jp,k,m,1))*ui(j,k,m,1)
      ctv(j,k,m) = 0.5*(t(j,k,m,1) + t(j,km,m,1))*vi(j,k,m,1)
c
c  Compute the several different velocity and transport fields needed
c  for momentum advection.
c
      uuhath(j,k,m) = (h(j,k,m,1)*(ui(j,k,m,1) + ui(jm,k,m,1))**2)/4. 
      vvhath(j,k,m) = (h(j,k,m,1)*(vi(j,k,m,1) + vi(j,kp,m,1))**2)/4. 
      hats = (h(j,k,m,1) + h(jp,k,m,1) + h(j,km,m,1) + h(jp,km,m,1))/4.
      uvhats(j,k,m) = hats*
     1     (ui(j,k,m,1) + ui(j,km,m,1))*(vi(j,k,m,1) + vi(jp,k,m,1))/4.
c
  100 continue 
c
c  Compute the diffusion terms, using either a Laplacian or biharmonic 
c  (if ibhm=1) form.
c
      call d24(ui,nx,ny,nz,dx,dy,diff,ibhm,work24,jmd,jpd,kmd,kpd,ud24)
      call d24(vi,nx,ny,nz,dx,dy,diff,ibhm,work24,jmd,jpd,kmd,kpd,vd24)
      call d24(t,nx,ny,nz,dx,dy,diff,ibhm,work24,jmd,jpd,kmd,kpd,td24)  
c
      wwork = 0.
      fwork = 0.
c
      do 200 j=1,nx
      jp = jpa(j)
      jm = jma(j)
c
      do 200 k=1,ny
      kp = kpa(k)
      km = kma(k)
c
      ua(j,k) = 0.
      va(j,k) = 0.
c      
c 
c  Apply wind stress and bottom drag to the upper and lower layers.
c      
      wx = wtfac*taux(j,k)  
      wy = wtfac*tauy(j,k)
      u(j,k,1,3) = u(j,k,1,3) + wx
      v(j,k,1,3) = v(j,k,1,3) + wy
c
      bdx = cbd*ui(j,k,nz,1)*(h(j,k,nz,1) + h(jp,k,nz,1))/2.
      bdy = cbd*vi(j,k,nz,1)*(h(j,k,nz,1) + h(j,km,nz,1))/2.
      u(j,k,nz,3) = u(j,k,nz,3) - bdx
      v(j,k,nz,3) = v(j,k,nz,3) - bdy 
c                 
            
c
c  Evaluate the east-west momentum budget. Note that when evaluating 
c  the Coriolis terms, you have to recall that f is defined at the h points.
c
      fp = 0.5*(f(k) + f(kp))
      fm = 0.5*(f(k) + f(km)) 
c            
      do 200 m=1,nz
c      
      hxa = (h(j,k,m,1) + h(jp,k,m,1))/2. 
      hya = (h(j,k,m,1) + h(j,km,m,1))/2. 
      udiff = ud24(j,k,m)*hxa 
      vdiff = vd24(j,k,m)*hya
c      
      avf = 0.25*(  fm*(v(j,k,m,1) + v(jp,k,m,1)) + 
     1  fp*(v(j,kp,m,1) + v(jp,kp,m,1)) )
      uac = avf + udiff - sdx*(p(jp,k,m) - p(j,k,m))*hxa
      uadv = sdx*(uuhath(jp,k,m) - uuhath(j,k,m)) + 
     1  sdy*(uvhats(j,kp,m) - uvhats(j,k,m)) 
c      
      u(j,k,m,3) = u(j,k,m,3) + uac - rnonl*uadv 
      ua(j,k) = ua(j,k) + uac - rnonl*uadv
c
c  North-south momentum budget.
c
      auf = 0.25*( f(k)*(u(j,k,m,1) + u(jm,k,m,1)) + 
     1  f(km)*(u(jm,km,m,1) + u(j,km,m,1)) )
      vac = vdiff - ( auf + sdy*(p(j,k,m) - p(j,km,m))*hya )
      vadv = sdx*(uvhats(j,k,m) - uvhats(jm,k,m)) + 
     1 sdy*(vvhath(j,k,m) - vvhath(j,km,m))                  
c     
      v(j,k,m,3) = v(j,k,m,3) + vac - rnonl*vadv 
      va(j,k) = va(j,k) + vac - rnonl*vadv  
c
c  A tracer budget.
c
      t(j,k,m,3) = t(j,k,m,3) + 0.25*td24(j,k,m) - 
     1 ( sdx*(ctu(j,k,m) - ctu(jm,k,m)) 
     2       + sdy*(ctv(j,kp,m) - ctv(j,k,m)) )   ! note the  reduced diffusivity          
c
c  The continuity equation.
c 
      h(j,k,m,3) = h(j,k,m,3)  - 
     1 ( sdx*(u(j,k,m,1) - u(jm,k,m,1)) + 
     2    sdy*(v(j,kp,m,1) - v(j,k,m,1)) ) 
c 
c  Compute the work by the windstress on the surface current.
c
      wwork = wwork + r0*wa(m)*( wx*ui(j,k,m,1) + wy*vi(j,k,m,1) )
c 
c  And compute the frictional work on both layers (ignores bottom drag)
c
      fwork = fwork +  r0*( udiff*ui(j,k,m,1) +  vdiff*vi(j,k,m,1) )    
  200 continue 
c
c .......................... Apply a Rigid Lid Condition ...................
c 
      if(rlid.ne.1) go to 6000
c
c  The next section applies the rigid lid constraint in one of two ways. The
c  first and evidently the better way is to require that the net force have 
c  zero divergence; the second and conventional way is to require that a
c  barotropic potential vorticity  conservation holds.
c      
c  just a reminder ....
c
c      h u | h u | h u | h
c      v s | v s | v s | v
c      --------------- . 
c      h u | h u | h u | h
c      v s | v s | v s | v
c      --------------- . 
c      h u | h u | h u | h
c      v s | v s | v s | v
c      _______________ .
c      h u | h u | h u | h
c
c  set ivort = 1 to do the curl version.
c
      ivort = 1
      if(ivort.eq.1) go to 6090
c      
c  Compute the rigid lid pressure by requiring that the divergence of the net
c  force, ua + dprdx (for x component), must vanish if the column height is to
c  remain constant for all times.
c      
c  Apply some solid wall boundary conditions on ua, va (the 'force') before 
c  computing the divergence if this is not a channel or symmetric domain
c     
      if(cyclic.ne.2) then       
      do 6999 j=1,nx
      va(j,1) = 0.
      va(j,ny) = 0.
 6999 continue 
      end if
      if(cyclic.lt.1) then
      do 6998 k=1,ny
      ua(1,k) = 0.
      ua(nx,k) = 0.
 6998 continue 
      end if
c
c  Compute the divergence of the force. This will be the source term in a 
c  Poisson equation for the rigid lid pressure, pr. pr is solved for on a 
c  grid that is one row and one column wider than the regular grid so that h
c  is on the outside all the way around the regular grid (see above). The
c  extra row is along the bottom, and the extra column is on the right side.
c  This also allows sensible boundary conditions (zero normal derivative) to 
c  be placed upon pr (apparently). 
c     
      j1 = 2
      if(cyclic.ge.1) j1 = 1
      do 7700 j=j1,nx
      do 7700 k=2,ny
      d2pr(j,k) = sdx*(ua(j,k-1) - ua(jma(j),k-1)) + 
     1   sdy*(va(j,kpa(k-1)) - va(j,k-1))
 7700 continue                             
c
c  Use a successive relaxation method to solve for the rigid lid pressure.
c  This subr assumes that the boundary condition is zero normal derivative
c  on all sides (i.e., this works for an enclosed basin only).
c       
cccccc      call
cccccc     1  sorpr(ncycle,d2pr,cyclic,dx,dy,nxp,nyp,jmap,jpap,kmap,kpap,pr) 
c
c  prm is the surface displacement that corresponds to the rigid lid pressure.
c      
      do 7760 j=1,nx
      do 7760 k=1,ny
      prm(j,k) = pr(j,k+1)/(htot*g*(rtcl/r0))
 7760 continue
c
c  Now add the rigid lid pressure gradient onto the force so that the net force
c  will have zero divergence. 
c
      twoh = 2.*htot
c              
      do 7750 j=2,nxm
      jp = jpa(j)
      do 7750 k=2,nym 
      km = kma(k)
c
c  In the lines below the k index is shifted +1 to account for the 
c  extra row in pr.
c      
      dprdx = sdx*(pr(jp,k+1) - pr(j,k+1))
      dprdy = sdy*(pr(j,k+1) - pr(j,k))
c                     
      do 7750 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) -
     1   dprdx*(h(j,k,m,1) + h(jp,k,m,1))/twoh 
      v(j,k,m,3) = v(j,k,m,3) - 
     1   dprdy*(h(j,k,m,1) + h(j,km,m,1))/twoh
c
 7750 continue 
c
c  One row and one column were not specified by the bcs and were missed 
c  in the above; do the u component along the first row.
c     
      k = 1
      do 7751 j=2,nxm
      dprdx = sdx*(pr(j+1,k+1) - pr(j,k+1))
      do 7751 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) - dprdx*(h(j,k,m,1) + h(j+1,k,m,1))/twoh
 7751 continue
c
c  Now do the v component in the last column. 
c     
      j = nx
      do 7752 k=2,nym
      dprdy = sdy*(pr(j,k+1) - pr(j,k))
      do 7752 m=1,nz
      v(j,k,m,3) = v(j,k,m,3) - dprdy*(h(j,k,m,1) + h(j,k-1,m,1))/twoh
 7752 continue      
c 
c 
c  That takes care of the pressue-inferred version of the rigid lid 
c  skip over the vorticity version below.
c
      go to 6000                      
c
 6090 continue
c  
c 
c  Come here to apply a rigid lid constraint by requiring that the barotropic
c  vorticity follow a constant depth conservation law.
c
c  Compute the curl of the force saved in ua,va.
c
      do 6060 j=1,nx
      jp = jpa(j)
      do 6060 k=1,ny
      km = kma(k)
      del2st(j,k) = sdx*(va(jp,k) - va(j,k)) - sdy*(ua(j,k) - ua(j,km)) 
c
 6060 continue
c
c  Solve a Poisson equation for the time rate change of the barotropic
c  streamfunction, st. This subroutine assumes that st is zero on the
c  boundaries.
c
c CCCCCC      call sor(del2st,dx,dy,st,nx,ny,nits)
c
c  Require that the PE momentum equation be consistent with this
c  barotropic streamfunction.
c 
      do 6061 j=2,nx
      jm = jma(j)
      jp = jpa(j) 
      do 6061 k=1,nym 
      km = kma(k)
      kp = kpa(k)
      dsdx = sdx*(st(j,k) - st(jm,k))                                          
      dsdy = sdy*(st(j,kp) - st(j,k))
      dprdx = (ua(j,k) + dsdy)/htot
      dprdy = (va(j,k) - dsdx)/htot
      do 6061 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) - dprdx*(h(j,k,m,1) + h(jp,k,m,1))/2.
      v(j,k,m,3) = v(j,k,m,3) - dprdy*(h(j,k,m,1) + h(j,km,m,1))/2.
 6061 continue
c 
c  End of the rigid lid calculation.
c
 6000 continue
c
c  If you want to insure a rigid lid then require that the layer
c  tendencies should sum to zero (set ijk = 1 to skip this).  
c     
      ijk = 0
      if(ijk.eq.1.or.nz.eq.1) go to 4682
      do 4680 j=1,nx
      do 4680 k=1,ny 
      hs = 0.
      do 4681 m=1,nzm
 4681 hs = hs + h(j,k,m,3)
      h(j,k,nz,3) = -hs 
 4680 continue               
 4682 continue
c    
c  The end of the rigid lid (tg).       
c 
      ibc = 1    !  this for radiation bc all the way around
c	
c  Apply boundary conditions.
c
c     
      cfact = 1.05
      cphr = cfact*cph    
      do 3500 k=1,ny 
      do 3500 m=1,nz
c
c  Apply a radiation BC to the east and west sides, assuming that cph is 
c  the relevent speed. 
c
      u(1,k,m,3) = cphr*( u(2,k,m,1) - u(1,k,m,1) )/dx
      u(nx,k,m,3) = -cphr*( u(nx,k,m,1) - u(nxm,k,m,1) )/dx
      v(2,k,m,3) = cphr*( v(3,k,m,1) - v(2,k,m,1) )/dx
      v(nx,k,m,3) = -cphr*( v(nx,k,m,1) - v(nxm,k,m,1) )/dx      
      h(2,k,m,3) = cphr*( h(3,k,m,1) - h(2,k,m,1) )/dx
      h(nx,k,m,3) = -cphr*( h(nx,k,m,1) - h(nxm,k,m,1) )/dx
c
c  for solid wall BCs on the east and west for an equatorial box (sort of) 
c 
      ibc = 2
      if(ibc.eq.2) then
	u(1,k,m,3) = 0.
	u(nx,k,m,3) = 0.
	v(2,k,m,3) = v(3,k,m,3)
	v(nx,k,m,3) = v(nxm,k,m,3)
	h(2,k,m,3) = h(3,k,m,3)
	h(nx,k,m,3) = h(nxm,k,m,3)
	end if
c
 3500 continue
c
      do 3501 j=1,nx
      do 3501 m=1,nz 
c
c  apply a radiation BC to the north and south sides too.
c
      u(j,1,m,3) = cph*( u(j,2,m,1) - u(j,1,m,1) )/dy
      u(j,nym,m,3) = -cph*( u(j,nym,m,1) - u(j,nym-1,m,1) )/dy
      v(j,1,m,3) = cph*( v(j,2,m,1) - v(j,1,m,1) )/dy
      v(j,ny,m,3) = -cph*( v(j,ny,m,1) - v(j,nym,m,1) )/dy
      h(j,1,m,3) = cph*( h(j,2,m,1) - h(j,1,m,1) )/dy
      h(j,nym,m,3) = -cph*( h(j,nym,m,1) - h(j,nym-1,m,1) )/dy 
c
c  or solid walls on the north and south sides....
c
      if(ibc.eq.3) then
      u(j,1,m,3) = u(j,2,m,3)
	u(j,nym,m,3) = u(j,nym-1,m,3)
	v(j,1,m,3) = 0.
	v(j,ny,m,3) = 0.
	h(j,1,m,3) = h(j,2,m,3)
	h(j,nym,m,3) = h(j,nym-1,m,3)
	end if
c
 3501 continue


c 
c
c ................................ Time Step ................................
c
c  Step ahead, noting whether this step is a:
c   straight leap-frog (leap = 0), or a 
c   leap-frog to be followed by a trapeozoidal correction (leap = 1), or
c   the trapezoidal correction itself (leap = 2).
c
      if(mod(ncycle,15).le.2) leap = leap + 1 
c     
      if(leap.eq.0) then
c      
      do 300 j=1,nx
      do 300 k=1,ny
      do 300 m=1,nz 
c      
      uold = u(j,k,m,2) 
      u(j,k,m,2) = u(j,k,m,1)
      u(j,k,m,1) = uold + twodt*u(j,k,m,3)
      u(j,k,m,3) = 0.
c      
      vold = v(j,k,m,2) 
      v(j,k,m,2) = v(j,k,m,1)
      v(j,k,m,1) = vold + twodt*v(j,k,m,3)
      v(j,k,m,3) = 0.      
c      
      hold = h(j,k,m,2) 
      h(j,k,m,2) = h(j,k,m,1)
      h(j,k,m,1) = hold + twodt*h(j,k,m,3)
      h(j,k,m,3) = 0.
c
      told = t(j,k,m,2) 
      t(j,k,m,2) = t(j,k,m,1)
      t(j,k,m,1) = told + twodt*t(j,k,m,3)
      t(j,k,m,3) = 0.
  300 continue
c 
c  Sum and normalize the wind and frictional work. 
c     
      rnxy = float(nx*ny)
      if(cyclic.eq.0) rnxy = float((nx-1)*(ny-1))
      wworks = wworks + dt*wwork/(rnxy*r0)
      fworks = fworks + dt*fwork/(rnxy*r0)
c      
      go to 303
      end if
c
      if(leap.eq.1) then 
c      
      do 301 j=1,nx
      do 301 k=1,ny
      do 301 m=1,nz 
c  
      u(j,k,m,2) = u(j,k,m,1)
      u(j,k,m,1) = u(j,k,m,2) + twodt*u(j,k,m,3)
c  
      v(j,k,m,2) = v(j,k,m,1)
      v(j,k,m,1) = v(j,k,m,2) + twodt*v(j,k,m,3)
c  
      h(j,k,m,2) = h(j,k,m,1)
      h(j,k,m,1) = h(j,k,m,2) + twodt*h(j,k,m,3)
c 
      t(j,k,m,2) = t(j,k,m,1)
      t(j,k,m,1) = t(j,k,m,2) + twodt*t(j,k,m,3)
  301 continue 
      go to 303
      end if
c
      if(leap.eq.2) then
c      
      do 302 j=1,nx
      do 302 k=1,ny
      do 302 m=1,nz 
c      
      uold = u(j,k,m,2) 
      u(j,k,m,1) = uold + halfdt*u(j,k,m,3)
      u(j,k,m,3) = 0.
c      
      vold = v(j,k,m,2) 
      v(j,k,m,1) = vold + halfdt*v(j,k,m,3)
      v(j,k,m,3) = 0.
c      
      hold = h(j,k,m,2) 
      h(j,k,m,1) = hold + halfdt*h(j,k,m,3)
      h(j,k,m,3) = 0.
c
      told = t(j,k,m,2) 
      t(j,k,m,1) = told + halfdt*t(j,k,m,3)
      t(j,k,m,3) = 0.
  302 continue 
c
c  Sum and normalize the wind and frictional work. 
c     
      rnxy = float(nx*ny)
      if(cyclic.eq.0) rnxy = float((nx-1)*(ny-1))
      wworks = wworks + dt*wwork/(rnxy*r0) 
      fworks = fworks + dt*fwork/(rnxy*r0)
c      
      end if  
c
  303 continue
c
c............................... Boundary Conditions .........................        
c        
c  Apply free slip or no slip BCs to all solid boundaries.
c
      if(cyclic.eq.2) go to 3100  
c
c  Apply bcs to the east and west sides of a fully enclosed basin.
c             
      if(cyclic.eq.1) go to 3003
c      
      do 3000 k=1,ny 
      do 3000 m=1,nz
c            
c      u(1,k,m,1) = 0. 
c      u(nx,k,m,1) = 0. 
c      
c  Apply no slip BCs to v if required.
c
c      if(frslip.eq.0) then 
c      v(2,k,m,1) = 0.
c      v(nx,k,m,1) = 0. 
c     end if 
c      
c  set the junk column 1 h and v to zero normal derivative      
c
      h(1,k,m,1) = h(2,k,m,1)
      t(1,k,m,1) = t(2,k,m,1)
      v(1,k,m,1) = v(2,k,m,1)
 3000 continue
c
 3003 continue
c
c  South and north sides, done for channel or basin.
c
      do 3001 j=1,nx
      do 3001 m=1,nz 
c      v(j,1,m,1) = 0.
c      v(j,ny,m,1) = 0. 
c      
c  Apply no slip BCs to u if required.      
c
c      if(frslip.eq.0) then
c      u(j,1,m,1) = 0.
c      u(j,nym,m,1) = 0.
c      end if
c      
c  Set the junk row ny u and h to zero normal derivative. 
c
      h(j,ny,m,1) = h(j,nym,m,1)
      u(j,ny,m,1) = u(j,nym,m,1)
      t(j,ny,m,1) = t(j,nym,m,1) 
 3001 continue
c 
 3100 continue
c 
c  End of the boundary conditions.
c
c  Go back up and re-do calculations for the trapezoidal correction if this
c  step has leap = 1.
c
      if(leap.eq.1) go to 901  
c
c................... Housekeeping Before the End of Time Step ...............
c
c  Step ahead the floats.
c                        
      if(nf.gt.0.and.ncycle.gt.1) then
c      
      lev = 1
      do 7300 i=1,nf
      xx = xf(i)
      yy = yf(i)
      call floatxy(xx,yy,ui,vi,x,y,nx,ny,nz,lev,dt,xxn,yyn,uf,vf)
      xf(i) = xxn
      yf(i) = yyn 
 7300 continue      
c
c  Store the float data on unit 35 if sufficient time has passed.
c
      if(mod(ncycle,50).eq.2) then
      do 7303 i=1,nf
      xyfk(i) = xf(i)
      xyfk(i+nf) = yf(i)
 7303 continue     
      write (35,7302) (xyfk(i),i=1,2*nf)
 7302 format (20f12.3)
      end if
c
      end if
c           
c                  
c  Check to see if time is up, or if you need to do the starting timestep.
c
      if(dtime.gt.runto) then 
      write (6,2255) 
 2255 format (1x,' reached end of time ')
      go to 9999
      end if
      if(ncycle.eq.1.and.irest.ne.1) go to 310 
c
c  Check for blowups, occasionally.
c
      ublow = 8.
c     
      if(mod(ncycle,10).eq.1) then
      umax = 0.
      etamax = 0.     
      do 800 j=1,nx
      do 800 k=1,ny
      do 800 m=1,nz
      eta8 = h(j,k,m,1) - h0(m)
      if(abs(eta8).gt.etamax) then
      jemax = j
      kemax = k
      memax = m
      etamax = eta8
      end if
      v8 = sqrt(ui(j,k,m,1)**2 + vi(j,k,m,1)**2)
      if(v8.gt.umax) then
      jvmax = j
      kvmax = k
      mvmax = m
      umax = v8
      end if
c
      if(umax.gt.ublow) then
      write (6,803) jvmax, kvmax, mvmax, ui(j,k,m,1), vi(j,k,m,1),
     1    h(j,k,m,1) 
  803 format (1x,' j, k, m, ui, vi, h are ', 3i5, 3f10.3)
      go to 9999
      end if
c
      etablow = 0.8*h0(m)
c      
      if(etamax.gt.etablow) then
      write (6,804) jemax, kemax, memax,
     1    h(j,k,m,1), u(j,k,m,1), v(j,k,m,1)
  804 format (1x,'etablow;  j, k, m, h, u, v', 3i6, 4f10.3)
      go to 9999
      end if
c
  800 continue
      end if 
c
c  End of the blowup check.
c
c
c  Evaluate the  voritcity and store in the tracer matrix 

      do 5207 j=2,nxm
	do 5207 k=2,nym
c	t(j,k,1,1) = ( f(k) + (vi(j+1, k,1,1) - vi(j-1,k,1,1))/(2*dx) - 
c     x  (ui(j,k+1, 1,1) - ui(j,k-1,1,1))/(2*dy) )/h(j,k,1,1)

	t(j,k,1,1) = (vi(j+1, k,1,1) - vi(j-1,k,1,1))/(2*dx) - 
     x  (ui(j,k+1, 1,1) - ui(j,k-1,1,1))/(2*dy) 

 5207 continue
c
c 
c................................ End of a Time Step .........................
c
c  Go back up to start another time step.
c
      go to 90
c
c  Come here only once at the start of the integration. This gives an old and a
c  present time step needed to begin leap-frogging.
c 
  310 continue
c    
      if(ncycle.eq.1) then
      do 400 j=1,nx
      do 400 k=1,ny
      do 400 m=1,nz  
      u(j,k,m,2) = u(j,k,m,1)
      v(j,k,m,2) = v(j,k,m,1)
      h(j,k,m,2) = h(j,k,m,1)
      t(j,k,m,2) = t(j,k,m,1)
      u(j,k,m,1) = u(j,k,m,1) + dt*u(j,k,m,3)
      v(j,k,m,1) = v(j,k,m,1) + dt*v(j,k,m,3)
      t(j,k,m,1) = t(j,k,m,1) + dt*t(j,k,m,3)
      h(j,k,m,1) = h(j,k,m,1) + dt*h(j,k,m,3)
  400 continue            
      end if
c           
c  This time step is complete; go back and do another one.
c
      go to 90
c
c................................. Close Up ..................................
c
 9999 continue
c                    
      irest = 0
c      write (6,7050)
c 7050 format (1x,' do you want to write a restart file ? (1=yes)')
c      read (5,*) irest
      if(irest.eq.1) then
      write (6,7052)
 7052 format (1x,' enter the name of restart file  ')     
      read (5,7051) ifile
 7051 format (a24)
      open(unit=20,file=ifile,form='unformatted',status='unknown')
c     
      rewind (20) 
      write (20) nx,ny,nz,nf     
      write (20) dtime,timd,timm,timt,
     1 frslip,cyclic,rgrav,rlid,ibhm,
     2 rnonl,dt,dx,dy,diff,htot,h0,delr0,rtcl,x,y,f,f0,beta,tau,gprime,
     3 wworks,fworks,xf,yf,taux,tauy,u,v,h,t,pr,ui,vi 
      end if
c
      stop 
      end              !   end of the main program 
c  
c  From here on down are subroutines needed to complete this monster.
c  Some of these are not used, of course.
c     
       subroutine kepe(h,u,v,eta,h0,delr0,nx,ny,nz,
     1 rgrav,cyclic,rke1,rketc,rpe,ah1,ah2)
c      
c  Compute kinetic and potential energies, and the mean layer thicknesses.
c
      dimension h(nx,ny,nz,3), u(nx,ny,nz,3), v(nx,ny,nz,3),
     1   h0(nz), delr0(nz), eta(nx,ny,nz)
c
      integer  rgrav, cyclic
      save
c                                                                        
      r0 = 1025.
      g = 9.8
c
      ah1s = 0.
      ah2s = 0.
      rke = 0.
      rke1 = 0.
      rpe = 0.
c
      rxy = 0.
      rke = 0.
      rpe = 0.
c 
c  Determine the portion of the domain to integrate over; if this is
c  a basin, then the left column and top row are not included, since they
c  acquire their values solely from the bcs. 
c                                  
      j1 = 1
      j2 = nx
      if(cyclic.eq.0) j1 = 2
      k1 = 1
      k2 = ny
      if(cyclic.le.1) k2 = ny - 1
c      
      do 659 j=j1,j2
      do 659 k=k1,k2
c
      rxy = rxy + 1.
      etas = 0. 
c 
c  If this is a free surface model then:
c     
      if(rgrav.eq.0) then
      do 640 m=nz,1,-1
      etas = etas + (h(j,k,m,1) - h0(m))
      eta(j,k,m) = etas
  640 continue 
      end if
c
c  If this is a reduced gravity model then:
c
      if(rgrav.eq.1) then
      do 642 m=1,nz-1
      etas = etas - (h(j,k,m,1) - h0(m))
      eta(j,k,m) = etas
  642 continue
      end if              
c
      do 650 m=1,nz
      rpe = rpe + delr0(m)*(eta(j,k,m)**2)
      rke = rke + h(j,k,m,1)*(u(j,k,m,1)**2 + v(j,k,m,1)**2)
  650 continue
c
      rke1 = rke1 + h(j,k,1,1)*(u(j,k,1,1)**2 + v(j,k,1,1)**2)
      ah1s = ah1s + h(j,k,1,1)
      ah2s = ah2s + h(j,k,nz,1)
c
  659 continue 
c
      rketc = rke - rke1
c
      rke = 0.5*rke/rxy
      rpe = 0.5*g*rpe/(rxy*r0)
      rke1 = 0.5*rke1/rxy
      rketc = 0.5*rketc/rxy
      ah1 = ah1s/rxy            
      ah2 = ah2s/rxy
c
      return
      end
c  
      subroutine tracer(t,nx,ny,ta,tsa,dtsa)
c
c  Compute some statistics on the tracer, t.
c
      dimension t(nx,ny)
c
      rn = 0.
      tas = 0.
      tsas = 0.
      dtsas = 0.
c
      do 1 j=2,nx
      do 1 k=1,ny-1
      rn = rn + 1.0
      tas = tas + t(j,k)
      tsas = tsas + t(j,k)**2
      dtsas = dtsas + sqrt( (t(j,k) - t(j,k+1))**2 + 
     1  (t(j,k) - t(j-1,k))**2 )
    1 continue
c
      ta = tas/rn
      tsa = tsas/rn
      dtsa = dtsas/rn
c
c  Normalize by the average t.
c
      tsa = tsa/(ta**2)
      dtsa = dtsa/(ta**2)      
c
      return
      end
c     
      subroutine sor(f,dx,dy,u,nx,ny,nits)
c
c  This subroutine solves the Poisson equation, DelSq(u) = f, by simultaneous 
c  overrelaxation. Taken from Numerical Recipes, and cleaned up. all errors
c  are due to J. Price.
c
c  f is input, and is the source term, dx and dy are the grid intervals
c  u is a first guess at the answer on input, and is returned as the answer.
c      
      dimension f(nx,ny),u(nx,ny)
      real omega
c 
      rn = sqrt(float(nx**2) + float(ny**2))
      rjac = cos(3.1415/rn)      
      maxits = 600
      eps = 3.0e-5
      anormf = 0.
c
      a = 1./(dx**2)
      b = 1./(dy**2)      
      c = -2.*(dx**2 + dy**2)/((dx**2)*(dy**2)) 
c
      do 12 j=2,nx-1
      do 11 k=2,ny-1         
      anormf = anormf + abs(f(j,k))
   11 continue
   12 continue
c  
c  Check to see if the field f might be essentially zero.
c    
      if(anormf.lt.1.e-12) return
c   
      omega = 1.0
      do 16 n=1,maxits
      anorm = 0.
      jsw = 1
c      
      do 15 ipass=1,2
c 
      k1 = ipass + 1     
      do 14 j=2,nx-1
c                   
      do 13 k=k1,ny-1,2
      resid = b*(u(j+1,k) + u(j-1,k)) + a*(u(j,k+1) + u(j,k-1)) + 
     1  c*u(j,k) - f(j,k)
      anorm = anorm + abs(resid)
      u(j,k) = u(j,k) - omega*resid/c
c      
   13 continue
   14 continue
c      
      if(n.eq.1.and.ipass.eq.1) then
      omega = 1.0/(1.0 - 0.5*rjac**2)
      else
      omega = 1.0/(1.0 - 0.25*rjac**2*omega)
      endif
c      
   15 continue
c     
      nits = n
      if(anorm.lt.(eps*anormf)) return
c          
   16 continue 
c   
       write (6,1) anorm, anormf 
    1 format (1x, '  woops, maxits exceeded within sor ', 2e10.3)
c    
      end 
c  
      subroutine 
     1 d24(a,nx,ny,nz,dx,dy,diff,ibhm,work,jmd,jpd,kmd,kpd,ad24)
c
c  Compute the Laplacian or the biharmonic (if ibhm = 1) of the field a.
c
      dimension a(nx,ny,nz,3),work(nx,ny,nz),ad24(nx,ny,nz),
     1 jmd(nx),jpd(nx),kmd(ny),kpd(ny)
c     
      dx2 = diff/dx**2
      dy2 = diff/dy**2
c 
c  Compute the Laplacian either way.
c     
      do 1 j=1,nx
      jm = jmd(j)
      jp = jpd(j)
      do 1 k=1,ny
      km = kmd(k)
      kp = kpd(k)
      do 1 m=1,nz
c
      ad24(j,k,m) = 
     1   dx2*(a(jm,k,m,2) + a(jp,k,m,2) - 2.*a(j,k,m,2)) + 
     2      dy2*(a(j,km,m,2) + a(j,kp,m,2) - 2.*a(j,k,m,2))
    1 continue
c
      if(ibhm.eq.0) return
c
c  Come here to compute the biharmonic form.
c      
      do 2 j=2,nx-1 
      jm = jmd(j)
      jp = jpd(j)
      do 2 k=2,ny-1
      km = kmd(k)
      kp = kpd(k)
      do 2 m=1,nz
c
      work(j,k,m) = 
     1  0.25*( (ad24(jm,k,m) + ad24(jp,k,m) - 2.*ad24(j,k,m)) + 
     2     (ad24(j,km,m) + ad24(j,kp,m) - 2.*ad24(j,k,m)) )
    2 continue 
c      
      do 3 j=2,nx-1
      do 3 k=2,ny-1
      do 3 m=1,nz
      ad24(j,k,m) = work(j,k,m)
    3 continue  
c
      return
      end
c 
c
      subroutine floatxy(xi,yi,u,v,x,y,nx,ny,nz,lev,dt,xn,yn,uf,vf)
c
      dimension u(nx,ny,nz,3),v(nx,ny,nz,3),x(nx),y(ny)
c
c  This subroutine steps forward a float from an old position xi,yi to
c  a new position xy,yn based upon the current components at the
c  previous time level, and the present time level. The 
c  time step is dt. The basic problem here is to do a 2-D interpolation.
c  The velocity u,v that are input here must be centered on the h points
c  (same as the x,y coordinates).  These are not the C-grid velocities.
c
c      h u | h u | h u
c      v s | v s | v s
c      ---------------
c      h u | h u | h u 
c      v s | v s | v s 
c      ---------------
c      h u | h u | h u
c      v s | v s | v s  
c
c  This assumes that x,y are monotonic and evenly spaced, and that x,t,u 
c  are given in a homogenous system of units.
c
c  Written by Jim Price, August 1994.
c     
      call uvxy(u,v,x,y,nx,ny,nz,lev,xi,yi,uf,vf)
c     
c            
      xn = xi + dt*uf
      yn = yi + dt*vf
c
      return 
      end
c
      subroutine uvxy(u,v,x,y,nx,ny,nz,lev,xi,yi,uf,vf)
c 
c  This subroutine does a 2-D linear interpolation.
c      
      dimension u(nx,ny,nz),v(nx,ny,nz),x(nx),y(ny)
c      
      dx = x(2) - x(1)
      dy = y(2) - y(1) 
c 
c  Figure out where this float is in the grid.
c     
      jn1 = int((xi - x(1))/dx) + 1
      jn2 = jn1 + 1
      kn1 = int((yi - y(1))/dy) + 1
      kn2 = kn1 + 1 
c
      dudx1 = (u(jn1,kn1,lev) - u(jn2,kn1,lev))/(x(jn1) - x(jn2)) 
      dudx2 = (u(jn1,kn2,lev) - u(jn2,kn2,lev))/(x(jn1) - x(jn2))
      dist1 = abs(y(kn1) - yi)/abs(dy)
      dist2 = abs(y(kn2) - yi)/abs(dy)
      dudx = dist2*dudx1 + dist1*dudx2
c           
      dudy1 = (u(jn1,kn1,lev) - u(jn1,kn2,lev))/(y(kn1) - y(kn2)) 
      dudy2 = (u(jn2,kn1,lev) - u(jn2,kn2,lev))/(y(kn1) - y(kn2))
      dist1 = abs(x(jn1) - xi)/abs(dx)
      dist2 = abs(x(jn2) - xi)/abs(dx)
      dudy = dist2*dudy1 + dist1*dudy2
c 
c
      dvdx1 = (v(jn1,kn1,lev) - v(jn2,kn1,lev))/(x(jn1) - x(jn2))
      dvdx2 = (v(jn1,kn2,lev) - v(jn2,kn2,lev))/(x(jn1) - x(jn2))
      dist1 = abs(y(kn1) - yi)/abs(dy)
      dist2 = abs(y(kn2) - yi)/abs(dy)
      dvdx = dist2*dvdx1 + dist1*dvdx2
c           
      dvdy1 = (v(jn1,kn1,lev) - v(jn1,kn2,lev))/(y(kn1) - y(kn2)) 
      dvdy2 = (v(jn2,kn1,lev) - v(jn2,kn2,lev))/(y(kn1) - y(kn2)) 
      dist1 = abs(x(jn1) - xi)/abs(dx)
      dist2 = abs(x(jn2) - xi)/abs(dx)
      dvdy = dist2*dvdy1 + dist1*dvdy2
c
      delx = (xi - x(jn1))
      dely = (yi - y(kn1))
c
      uf = u(jn1,kn1,lev) + delx*dudx + dely*dudy
      vf = v(jn1,kn1,lev) + delx*dvdx + dely*dvdy 
c
      return
      end            
c   
c   
c  

c
      subroutine sorpr
     1  (nstep,nits,eps,f,icyclic,dx,dy,nx,ny,jm,jp,km,kp,pr) 
c
c
c  This subroutine solves the Poisson equation, DelSq(pr) = f, by simultaneous 
c  overrelaxation. Taken from Numerical Recipes, modified heavily, and 
c  documented by Jim Price. This subroutine is specifically intended to 
c  compute the rigid lid pressure, and assumes that the bcs are either zero 
c  normal derivative, or symmetric, depending upon the kind of domain as
c  specified by icyclic. 
c
c  f is input, and is the source term. It is assumed to be available over the
c   entire domain, but in the case of solid boundaries, is not used on the 
c   boundaries. 
c  icyclic = 0 for solid boundaries,
c          = 1 for a channel (cyclic in x direction only),
c          = 2 for a fully cyclic domain
c  dx and dy are input and are the grid intervals.
c  ny and nx are the array dimensions. These are one larger than in the main 
c  program. This is required in the case of all solid boundaries so that zero
c  normal derivative bcs can be used. The Poisson equation is solved over
c  the interior points only for a solid boundary, and over all points for a
c  cyclic domain.  
c  The integer arrays jm, jp, km, kp are dimensioned in the main program, and
c  are defined here in order to specify indices needed to
c  evaluate the Laplacian (needed to specify bcs; i.e, symmetric or not). 
c  pr is a first guess at the answer on input, and is returned as the 
c   answer on output. 
c
c      
      dimension f(nx,ny),pr(nx,ny),jm(nx),jp(nx),km(ny),kp(ny)
c
      save
c
c      if(nstep.gt.1) go to 100      
c
c  Compute rjac, an over-relaxation parameter used to aid convergence.  
c     
      xx = float(nx)
      yy = float(ny)
      rn = sqrt((xx**2) + (yy**2))
      cp = 3.1415925/rn
      rjac = cos(cp)
c 
c  These maxits and epsilon values are rather arbitrary; the choice of
c  the size of eps can be, however, crucial to the time the thing takes to run
c  (eps = 1.e-5 seems OK for the pe model, x0.5 might be better).
c      
      maxits = 500
      anormf = 0.
c
      a = 1./(dx**2)
      b = 1./(dy**2)      
      c = -2.*(dx**2 + dy**2)/((dx**2)*(dy**2)) 
c   
      romega = 1.0 
c
c  Just a reminder of the grid ....
c
c      h u | h u | h u | h
c      v s | v s | v s | v
c      --------------- . 
c      h u | h u | h u | h
c      v s | v s | v s | u
c      --------------- . 
c      h u | h u | h u | h
c      v s | v s | v s | u
c      _______________ .
c      h u | h u | h u | h
c
c  Make up the indices needed in the Laplacian.
c
      do 300 j=2,nx-1
      jp(j) = j + 1
      jm(j) = j - 1
  300 continue 
c
c  j1 and j2 are the first and last points on which the Laplacian 
c  of the pr field is computed. 
c
      j1 = 2
      j2 = nx - 1
c
c  Now, check to see if this is for a channel
c  (note that you do not have to deal with the last column if it is)
c     
      if(icyclic.ge.1) then
      j1 = 1
      j2 = nx - 1
      jm(1) = nx - 1
      jp(1) = 2   
      jp(nx-1) = 1  
      end if 
c
      do 301 k=2,ny-1
      kp(k) = k + 1
      km(k) = k - 1
  301 continue
c                      
c      k0 = 1
      k1 = 2
      k2 = ny - 1 
c
c  Check to see if this a completely symmetric domain
c  (in this case you don't deal with the first row).
c           
      if(icyclic.eq.2) then
c      k0 = 1
      k1 = 2
      k2 = ny 
      km(ny) = ny - 1
      kp(ny) = 2
      km(2) = ny    
      end if
c
      romega = 1.0/(1.0 - 0.5*rjac**2)
c
c      write (60,400) nstep, icyclic, j1, j2, k1, k2
 400  format (1x, ' ncylce,icyclic,j1,j2,k1,k2 are  ',8i7) 
      do 500 j=1,nx 
      k = j
c      write (60,501) j,jm(j),jp(j),k,km(k),kp(k)     
 501  format (1x,' j,jm,jp,k,km,kp', 7i4)
 500  continue     
c
  100 continue      
c     
      anormf = 0.
      do 12 j=2,nx-1
      do 11 k=2,ny-1         
      anormf = anormf + abs(f(j,k))
   11 continue
   12 continue
c  
c  Check to see if the field f might be essentially zero.
c    
      if(anormf.lt.1.e-19) then
c      write (6,30) anormf
   30 format (1x,' anormf seems to be very small, I quit.', e12.4)
      return
      end if
c          
      do 16 n=1,maxits
      anorm = 0.
      jsw = 1      
      do 15 ipass=1,2
      lsw = jsw
c      
      do 14 j=j1,j2
      jj1 = jm(j)
      jj2 = jp(j)
      if(nstep.eq.35) then  
c      write (60,59) j,j1,j2,jj1,jj2
c      write (60,58) lsw,k1,k2
  59  format (1x,' j,j1,j2,jm,jp', 7i4)
  58  format (1x,' lsw,k1,k2,' 7i6)    
      end if
c      
      do 13 k=k1+lsw-1,k2,2
      kk1 = km(k)
      kk2 = kp(k) 
c     
      resid = b*(pr(jp(j),k) + pr(jm(j),k)) +
     1 a*(pr(j,kp(k)) + pr(j,km(k))) +  c*pr(j,k) - f(j,k)
      anorm = anorm + abs(resid)
      pr(j,k) = pr(j,k) - romega*resid/c
c      
   13 continue
      lsw = 3 - lsw
   14 continue
      jsw = 3 - jsw
c      
      romega = 1.0/(1.0 - 0.25*rjac**2*romega)
c
c  Apply zero normal derivative bcs.
c
c  March along the north and south sides, noting whether this is
c  a cyclic domain or not.
c     
      if(icyclic.eq.2) go to 299 
c
c  Apply zero normal derivative bcs to a basin or a channel. 
c      
      do 20 j=1,nx
      pr(j,1) = pr(j,2) 
      pr(j,ny) = pr(j,ny-1)  
   20 continue   
  299 continue
c 
c  March along the east and west sides if this is not a channel.
c           
      if(icyclic.ge.1) go to 219
      do 21 k=2,ny-1
      pr(1,k) = pr(2,k) 
      pr(nx,k) = pr(nx-1,k) 
   21 continue   
  219 continue
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, romega, anorm, epsnorm  
    1 format (1x,'      maxits passed in sorpr ', i6, 5e9.3)
c
   99 continue
c
c  Demean pr if the bcs are zero normal derivative all the way around, or if
c  the domain is fully symmetric. 
c 
      if(icyclic.eq.0.or.icyclic.eq.2) then
      prm = 0.
      do 50 j=1,nx
      do 50 k=1,ny
      prm = prm + pr(j,k)
   50 continue
      prm = prm/float(nx*ny)
      do 51 j=1,nx
      do 51 k=1,ny
      pr(j,k) = pr(j,k) - prm
   51 continue 
      end if             
c 
      nits = n
c          
      return
      end
c   
c
      subroutine 
     1 sorpsi(nstep,f,icyclic,dx,dy,nx,ny,jm,jp,km,kp,psi,nits)
c
c  This subroutine solves the Poisson equation, DelSq(psi) = f, by 
c  the method of simultaneous overrelaxation. The basis for this code
c  was taken from Numerical Recipes and then modified heavily to be used
c  to solve for a streamfunction (here called psi) in the pe model 
c  written by Jim Price. July, 1994.
c   
c  nstep = 1 on first entering this subroutine.
c  f is input, and is the source term, 
c  icyclic = 0 for solid boundaries, and thus psi = 0 all around,
c          = 1 for a channel (cyclic in x direction only),
c          = 2 for a fully cyclic domain.
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  psi is a first guess at the answer on input, and is returned as the 
c   answer on output.  
c
c      
      dimension f(nx,ny),psi(nx,ny),jm(nx),jp(nx),km(ny),kp(ny)
      real omega 
c      
      save
c
c      if(nstep.gt.1) go to 100      
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 = 700
c      
      eps = 0.3e-5 
c           
c      eps = 0.4e-6
c
c  The very small value of eps used above seems to cause basin modes to
c  appear in the larger domains (nx, ny > 50, say).
c 
      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
      j1 = 2
      j2 = nx - 1
      if(icyclic.ge.1) then
      j1 = 1
      j2 = nx  
      end if 
      k0 = 1
      k2 = ny - 1      
      if(icyclic.eq.2) then
      k0 = 0
      k2 = ny   
      end if 
c     
      omega = 1.0/(1.0 - 0.5*rjac**2)
c      
  100 continue      
c
      do 12 j=2,nx-1
      do 11 k=2,ny-1         
      anormf = anormf + abs(f(j,k))
   11 continue
   12 continue
c  
c check to see if the field f might be essentially zero
c    
      if(anormf.lt.1.e-19) then
      write (6,30) anormf
   30 format (1x,' anormf seemed to be very small', e12.4)
      return
      end if
c          
      do 16 n=1,maxits
      anorm = 0.
      jsw = 1      
      do 15 ipass=1,2
      lsw = jsw
c      
      do 14 j=j1,j2
      jj1 = jm(j)
      jj2 = jp(j)  
c      write (60,59) j,j1,j2,jj1,jj2
      do 13 k=lsw+k0,k2,2
      kk1 = km(k)
      kk2 = kp(k) 
c      write (60,58) k,lsw,k0,k2,kk1,kk2
  59  format (1x,' j,j1,j2,jm,jp', 7i4)
  58  format (1x,' k,lsw,k0,k2,km,kp' 7i6)    
c     
      resid = b*(psi(jp(j),k) + psi(jm(j),k)) + 
     1  a*(psi(j,kp(k)) + psi(j,km(k)))
     1 +  c*psi(j,k) - f(j,k)
      anorm = anorm + abs(resid)
      psi(j,k) = psi(j,k) - omega*resid/c
c      
   13 continue
      lsw = 3 - lsw
   14 continue
      jsw = 3 - jsw
c      
   15 continue  
c
      omega = 1.0/(1.0 - 0.25*rjac**2*omega)   
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, anorm, epsnorm  
    1 format (1x,///, '   maxits passed in sorpsi ', i6, 5e9.3)
c
   99 continue
c
c  demean psi if the domain is fully symmetric. 
c 
      if(icyclic.eq.2) then
      do 50 j=1,nx
      do 50 k=1,ny
      psim = psim + psi(j,k)
   50 continue
      um = um/float(nx*ny)
      do 51 j=1,nx
      do 51 k=1,ny
      psi(j,k) = psi(j,k) - psim
   51 continue 
      end if             
c 
      nits = n
      return
      end