c
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 = 1.0e-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 
CCCC     anormf = 0.
      anormf = 1.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