c
      program peomod
c      
c
c   The main program of a layered, primitive equation ocean model. 
c   The necessary subrotuines are in basinsubs.f.
c
c   This model may be used to compute the response of an 
c   ocean basin or a channel to a prescribed wind field. It may
c   also be used to study geostrophic adjustment of an 
c   initial density anomaly.  This model is intended to be more an
c   educational tool than a first line research tool.
c 
c   The 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 by a leap-frog with an occasional
c   trapezoidal correction in order to supress time-splitting.
c   Hydrostatic pressure is computed by a reduced gravity or
c   a free surface pressure model. A rigid lid 'pressure' can 
c   also be applied.
c
c   The code is intended to be completely portable Fortran 77.
c
c
c   Written by Jim Price beginning in March, 1994. Changes were 
c   ongoing at least as late as December, 2005. 
c
c   NOTE: This model code is developmental, and has not been fully
c   tested.  
c
c   The variables are:
c   
c   (u, v) are the (east, north) components of the transport,
c   h is the layer thickness,
c   p is the hydrostatic pressure,
c   s is a streamfunction for the transport (or the net force if rigid lid), 
c   t is a passive tracer,
c   (taux, tauy) are (east, north) wind stress components,
c   (x, y) are east and north coordinates.
c
c   These 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   Pressure is defined at the h points, as is the tracer, t, and
c   and the Coriolis parameter, f. The x,y coordinates
c   are also centered on the h points.  Wind stress is
c   defined at the corresponding u,v points.  This grid scheme is 
c   a crucial aspect of the numerical method to which you have to pay
c   rigorous attention if you change any aspect of the finite differencing
c   or the boundary conditions!  For example, u is the east component of 
c   transport; to compute the velocity, ui, you must divide by the
c   appropriately averaged (horizontally) layer thickness as follows:
c   ui(j,k,m) = u(j,k,m)/(0.5*(h(j,k,m) + h(j+1,k,m))).  This is highly
c   tedious, but necessary to achieve accurate solutions. 
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   The variables are stored in four-dimensional arrays, i.e.,  
c   u(nx,ny,nz,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   As of December 1995, the following problems remained:
c
c   1) The most troubling (problematic) aspects of the model seem to be
c   almost anything to do with the rigid lid calculations.  The rigid lid
c   version works much better when the sum of the layer thicknesses are
c   required to equal the water column depth. This may be 
c   a warning of a small error in the rigid lid code. In any case, the long-
c   term behavior can be sensitive to the details of how the rigid lid 
c   pressure is solved (for times in excess of 1500 days), and e.g., the
c   rigid lid calculation works much better when double precision
c   arithmetic is used throughout (a must, actually). 
c
c   2) Part of the whining above probably comes from the 
c   rather simple elliptic solver included here.  These errors are made
c   noticably worse as the size of the model domain (number of points
c   in the x and y directions) is increased beyond about 40.    
c 
c   3) The various boundary condition options (i.e., basin, channel or fully
c   symmetric) and pressure models (reduced gravity, free surface and
c   rigid lid) have not all been checked in every possible configuration.
c
c
c............................. Data Allocation ..........................
c
      parameter (nx=51, ny=51, nz=2, nxp=nx+1, nyp=ny+1, 
     1    nf=4, nf2=2*nz*nf)
c
      dimension h(nx,ny,nz,3),u(nx,ny,nz,3),v(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),f(ny),taux(nx,ny),tauy(nx,ny),
     3 kpa(ny),jpa(nx),kma(ny),jma(nx),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     
      dimension work1(nx,ny),work2(nx,ny),work3(nx,ny),work4(nx,ny),
     1 work5(nx,ny),work6(nx,ny),work7(nx,ny),work8(nx,ny)         
c
c  The arrays below are required solely for the rigid lid calculation.
c
      dimension ua(nx,ny),va(nx,ny),del2st(nx,ny),st(nx,ny),
     1 prp(nxp,nyp),d2pr(nxp,nyp),prm(nx,ny),prx(nxp,nyp),pry(nxp,nyp),
     2 jpap(nxp),jmap(nxp),kmap(nyp),kpap(nyp),uapr(nx,ny),vapr(nx,ny),
     3 pr(nx,ny),dprdx(nx,ny),dprdy(nx,ny)     
c
c  The arrays below are needed for float tracking and float data storage.
c
      dimension xf(nf,nz),yf(nf,nz),xyfk(nf2),
     1 ufl(nx,ny,nz),vfl(nx,ny,nz)
c
c  These arrays are used to store the time averages and other diagnsotics. 
c
      dimension savg(nx,ny),pavg(nx,ny,nz),pvavg(nx,ny,nz),
     1 eavg(nx,ny,nz),avprof(ny,nz),pvprof(ny,nz),rvprof(ny,nz),
     2 uavg(nx,ny,nz),vavg(nx,ny,nz)
c      
      integer sidebc, basin, rgrav, rlid 
      character*60 ifile, alf 
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. 
      bworks = 0.
      timd = 0.
      timm = 0.
      timt = 0.
      iblow = 0
c      
c  rnonl can be set to zero to remove nonlinear terms from the momentum and 
c   tracer equations (but you had better be careful with this !).
      rnonl = 1.0
c      
c 
c ........................ Restart, and Open Files .........................                       
c
c  Check to see if you want to begin from a restart file.
c     
      irest = 0
      iapp = 0
      write (6,7000)
 7000 format (1x,//,
     1 '  Should we restart from a previous run?',/,
     2 '   restart?                    (0 = no (default), 1 = yes)',/,
     3 '   append data to old files?   (0 = no (default), 1 = yes)',/,
     4 '   (enter both data or default with ,,)')  
      read (5,*) irest, iapp 
c 
c  The files opened below are used to save results for later Matlab 
c  plotting program. Add data on to some of these these files if iapp 
c  is set on (= 1). 
c
      if(iapp.eq.1) then
      open (unit=32,file='looktim.dat',access='append',status='unknown')
      open (unit=1,file='basin1.dat',access='append',status='unknown')
      open (unit=2,file='basin2.dat',access='append',status='unknown')
      open (unit=3,file='basin3.dat',access='append',status='unknown')
      open (unit=4,file='basin4.dat',access='append',status='unknown')
      open (unit=15,file='basin5.dat',access='append',status='unknown')
      open (unit=16,file='basin6.dat',access='append',status='unknown')      
c       
      open (unit=10,file='basin10.dat',access='append',status='unknown')
      open (unit=40,file='basin40.dat',access='append',status='unknown')      
c      
      open (unit=31,file='movtim.dat',access='append',status='unknown')
      open (unit=11,file='basin11.dat',access='append',status='unknown')
c      
      open (unit=35,file='bfloat.dat',access='append',status='unknown')
      open (unit=36,file='pvprof.dat',access='append',status='unknown')      
c      
      else 
      open (unit=32,file='looktim.dat') 
      open (unit=1,file='basin1.dat')
      open (unit=2,file='basin2.dat')
      open (unit=3,file='basin3.dat')
      open (unit=4,file='basin4.dat')
      open (unit=15,file='basin5.dat')
      open (unit=16,file='basin6.dat')       
c      
      open (unit=10,file='basin10.dat')
      open (unit=40,file='basin40.dat')      
c      
      open (unit=31,file='movtim.dat')
      open (unit=11,file='basin11.dat')
c      
      open (unit=35,file='bfloat.dat')
      open (unit=36,file='pvprof.dat')      
c      
      end if 
c
      open (unit=12,file='basin12.dat',status='unknown')
      open (unit=30,file='what.dat')
c 
c  Unit 39 is connected to a file called stopfile.dat that has a single
c  number that can be set to 1 to gracefully stop a model run under DOS. 
c 
c      open (unit=39,file='stopfile.dat',status='unknown',
c     1        form='formatted')
c      istop = 0
c      write (39,3939) istop
c 3939 format (i1)
c      close (39)
c      rewind(39)            
c
c  Open the restart file here, if need be.
c      
      if(irest.eq.1) then   
      write (6,7052)
      read (5,7001) ifile
 7001 format (a24)
      open(unit=20,file=ifile,form='unformatted',status='old')
c      
      read (20) nx4,ny4,nz4,nf4
c
c  Verify that the dimensions are consistent with the restart file data set.
c
      if(nx4.ne.nx.or.ny4.ne.ny.or.nz4.ne.nz.or.nf4.ne.nf) then
      write (6,7010) nx4,ny4,nz4,nf4,nx,ny,nz,nf
 7010 format (1x,///,
     1'  Sorry bud, the dimensions in startup do not match.', 6i5)
      go to 9999
      end if           
c      
      read (20) time9,timd,timm,timt,
     1 sidebc,basin,rgrav,rlid,ibhm,cdb,
     2 rnonl,dt,dx,dy,diff,hbot,h0,delr0,rtcl,x,y,f,f0,beta,tau,gprime,
     3 wworks,fworks,bworks,xf,yf,taux,tauy,u,v,h,t,pr,ui,vi 
c
c  Write out some information on this case just as a reminder.
c     
      write (6,8811) time9, sidebc, basin, rgrav, rlid
 8811 format (1x,/,' time, sidebc, basin, rgrav, rlid',
     1 /, 1x, f12.2, 6i5,/)
      write (6,8812) dt, dx, dy, diff, hbot
 8812 format (1x,/,' dt, dx, dy, diff, hbot',/,1x,5f9.2,/)
c
      twodt = 2.*dt
      halfdt = 0.5*dt     
      sdx = 1./dx
      sdy = 1./dy
c
      go to 7100
      end if           
c 
c ......... Initialize Basin Geometry and Boundary Conditions ..............
c
c  The geometry and boundary conditions are specified via the flags:
c  
c  The basin geometry is set by: 
c  basin = 0 for solid walls all the way around (a closed basin, the default),
c  basin = 1 for open BCs on east and west sides and walls on the 
c   north and south sides (this makes a channel),
c  basin = 2 for open BCs all the way around (probably makes sense only if
c   beta is zero). This makes an open ocean.
c
c  The BC for tangential velocity on the sidewalls is set by:
c  sidebc = 0 for no-slip BCs on the sidewalls, or, 
c  sidebc = 1 for free-slip BCs on the sidewalls (the default).
c
c  The BC for open boundaries (if there are any) is set by:
c  openbc = 0 for wrap-around BCs, or,
c  openbc = 1 for radiation BCs on all open sides.

c  For example, the default is a closed basin with free-slip BCs, and thus the 
c  default is basin = 0 and sidebc = 1. To have no-slip boundary conditions 
c  in a channel you would need to set basin = 1 and sidebc = 0.  To have a
c  completely open domain with radiation BCs (say for a geostrophic adjustment 
c  problem) set basin = 2 and openbc =1. 
c
      basin = 0 
      sidebc = 1
      openbc = 0 
c
      write (6,9929) 
 9929 format
     1 (1x,/,'  Basin configuration: ',/,
     2 '   closed basin       = 0 (default), or,',/,
     3 '   east-west channel  = 1, or,',/,
     4 '   open all sides     = 2',/,
     5 '   (enter one, or default with ,)')
      read (5,*) basin            
c
      if(basin.ne.2) then
      write (6,9922) 
 9922 format
     1 (1x,/,'  Sidewall boundary conditions: '/,
     2 '   no-slip            = 0 ',/,
     3 '   free-slip          = 1 (default)',/,
     4 '   (enter one, or default with ,)')
      read (5,*) sidebc
      end if
c
      if(basin.ge.1) then
      write (6,9928) 
 9928 format (1x,/,'  Open boundary boundary conditions:',/,
     1 '  wrap-around         = 0 (default), or,',/,
     2 '  radiation BC        = 1',/,
     3 '  (enter one, or default with ,)')
      read (5,*) openbc
      end if      
c  
c ....................... Specify the Pressure Model .....................
c
c  The pressure model is set by the following flags:
c   
c  rgrav = 0 for a free surface pressure model, or,
c  rgrav = 1 for a reduced gravity pressure model (the default), and,
c  rlid = 0 (or 1) if you do (or do not) want to apply a rigid lid pressure
c   correction.  This is allowed only with a reduced gravity model.
c
c  Most of the time you will probably use a reduced gravity model.  The free
c  surface model is included as well mainly as a possible future alternative
c  to the damned rigid lid approximation.
c
c
 9135 continue 
c     
      rgrav = 1
      rlid = 2
c     
      write (6,9923)  
 9923 format (1x,/, 
     1 '  Pressure model: ',/,
     2 '   free surface       = 0',/,
     3 '   reduced gravity    = 1 (default)',/,
     4 '   (enter one, or default with ,)' )
      read (5,*) rgrav 
      if(rgrav.eq.1) then
      write (6,9972) 
 9972 format (1x,/,
     1 '  Apply a rigid lid ?',/,
     2 '   no                           = 0',/,
     3 '   yes                          = 1',/,
     4 '   yes, w/ corrected thickness  = 2 (default)',/,
     5 '   (enter one, or default with ,)' ) 
      read (5,*) rlid
      end if 
c
c 
c ....................... Initialize the Density Profile .....................
c
c  This next section sets the initial density 'profile' by defining the initial 
c  layer thicknesses and density differences as follows: 
c
c  hbot is the depth of the water column (m),
c  htcl is the thickness of the surface and thermocline layer(s) (m),
c  rtcl is the density difference across the thermocline (kg/m**3),
c  delradd is a density that may be added to the upper layer (i.e., to make
c   a free surface model with a realistic air-sea density difference you
c   should set delradd = 1025.).
c
c  From these input data the model then computes what it needs:
c
c  h0 are the initial layer thicknesses (m),
c  delr0 are the fixed density differences across the layer interfaces.
c 
c  The details of all of this are a little involved and depend upon the
c  kind of pressure model.  Be sure to check the resulting h0 and delr0.
c 
      hbot = 4000.
      htcl = 500.
      rtcl = 4.0 
      delradd = 0.
c                 
      write (6,9924) htcl, hbot 
 9924 format (1x,/, 
     1 '  Layer thicknesses and bottom depth: ',/,
     2 '   surface and thermocline layer    (default = ',f6.0,' m)',/,
     3 '   bottom depth                     (default = ',f6.0,' m)',/,
     4 '   (enter both data, or default with  ,,)') 
      read (5,*) htcl, hbot
c        
      write (6,9925) rtcl, delradd
 9925 format (1x,/,
     1 '  Density differences: ',/,
     2 '   thermocline density diff.  (default = ',f6.1,' kg/m**3)',/,
     3 '   surface density diff.      (default = ',f6.1,' kg/m**3)',/,
     4 '   (enter both data, or default with  ,,)') 
      read (5,*) rtcl, delradd
c            
c  Set the density profile for a reduced gravity model. In this case the
c  density differences are associated with the lower interface of a layer.
c  The layer thicknesses are all set to htcl/nz, and the abyssal layer below 
c  (which is not explicitly represented) is presumed motionless.
c          
      if(rgrav.eq.1.and.nz.eq.1) then
      h0(1) = htcl
      delr0(1) = rtcl
      end if
c
      if(rgrav.eq.1.and.nz.gt.1) then
      hs = 0.
      do 780 m=1,nz
      h0(m) = htcl/float(nz)
      hs = hs + h0(m)
      delr0(m) = rtcl/float(nz)
  780 continue     
      end if
c
c  Set the density profile for a reduced gravity model that will also have
c  a rigid lid condition applied.  In this case, the sum of the layers has
c  to sum up to hbot, and there is no density difference across the lowest
c  layer. 
c
      if(rlid.ge.1) then
c
      if(nz.eq.1) then
      h0(1) = htcl
      delr0(1) = 0.
      end if
c
      if(nz.gt.1) then
      nzm = nz - 1
      hs = 0.
      do 781 m=1,nzm
      delr0(m) = rtcl/float(nzm)
      h0(m) = htcl/float(nzm)
      hs = hs + h0(m)
  781 continue
      delr0(nz) = 0.
      h0(nz) = hbot - hs
      end if
c
      end if      
c
c  Set the density profile for a free surface model. In this 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
      h0(1) = hbot/(float(nz))
      if(nz.gt.1) then
      hs = 0.
      do 782 m=1,nzm
      delr0(m) = rtcl/float(nz)
      h0(m) = htcl/float(nzm)
      hs = hs + h0(m)
  782 continue 
      delr0(1) = delr0(1) + delradd
      delr0(nz) = rtcl/float(nz)
      h0(nz) = hbot - hs
      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
 785  continue          
c
c  Print out the density 'profile'.
c
      write (6,8821) 
 8821 format (1x, 
     1 ' These parameters gave the following layered profile:',/,
     2 '   layer no.     density diff.     thickness ')
      do 7854 m=1,nz
      write (6,8822) m, delr0(m), h0(m)
 8822 format (1x,5x, i4, 3x, f10.2, 8x, f7.0)             
 7854 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 = 25.e3 
      dy = 25.e3    
c
      dxk = dx/1000.
      dyk = dy/1000.
      write (6,442) dt, dxk, dyk
  442 format (1x,/,
     1 '  Time and space resolution:',/,
     2 '   time step                     (default = ',f7.0,' sec)',/,
     3 '   east-west space difference     (default = ',f7.0,' km)',/,
     4 '   north-south space difference   (default = ',f7.0,' km)',/,
     5 '   (enter all three data, or default with  ,,,)')          
      read (5,*) dt, dxk, dyk 
      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  (the present form of this works for two layers only)
c
      cph = sqrt(g*rtcl*(h0(1)*h0(nz)/hbot)/r0) + 0.1
      if(rgrav.eq.0) cph = sqrt(g*(delradd + rtcl)*hbot/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     
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 = 30.
      ibeta = 1 
c
      write (6,4467) rlat 
 4467 format (1x,/,
     1 '  Earths rotation:',/,
     2 '   latitude                     (default = ',f8.1,' deg)',/,
     3 '   beta   (0 = off, 1 = on (default), -1 (change sign) )',/,
     4 '   (enter both data or default with  ,,)')
      read (5,*) rlat, ibeta
c
      f0 = 2.*7.292e-5*sin(rlat*tpi/360.)
      beta = float(ibeta)*2.*7.292e-5*cos(rlat*tpi/360.)/6370.e3 
c 
c  Define east and north coordinates,  x and y, in meters, and the variable
c  Coriolis parameter, f.  Note that the origin of x, y is centered on 
c  model domain. 
c                     
      do 10 j=1,nx
      x(j) = float(j-ncx)*dx
   10 continue     
c      
      do 107 k=1,ny
      y(k) = float(k-ncy)*dy
      f(k) = f0 + y(k)*beta  
  107 continue
c
c  Estimate some long baroclinic 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.
c      
      bsize = sqrt(1./(x(nx)-x(1))**2 + 1./(y(ny)-y(1))**2 )
      btrm1 = ( (2.*3.1415)**2*bsize/beta )/8.64e4 
c
c  compute and print some higher Ro modes as well
c
      do 2259 m1 = 1,3
      do 2259 m2 = 1,3
      rmn = float(m1)
      rnn = float(m2)
      bsize = sqrt(1./
     c   ( rmn**2*(x(nx)-x(1))**2) + 1./(rnn**2*(y(ny)-y(1))**2) )
      btrm = ( (2.*3.1415)**2*bsize/beta )/8.64e4 
      write (6,2256) m1, m2, btrm
 2256 format (1x,' m1, m2, btrm1m2', 2i5, f10.3)
 2259 continue
c
      write (6,3045) f0, beta, cph, rdef, rcph, btrm1 
 3045 format (1x,//,1x,
     1' Given this latitude and density profile, the:       ',//,
     2'  middle f and beta are               ',    2e12.3,/,
     3'  long (nondispersive) gravity wave speed (m/sec) is ',f9.2,/,
     4'  radius of deformation (km) is                      ',f9.2,/,
     5'  long (nondispersive) baroclinic Ro wave speed is   ',f9.3,/,
     6'  period of the first barotropic Ro mode (days) is   ',f9.2)
      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.  The nominal wind stress
c  pattern is multiplied by tautim to give a seasonally varying wind;
c  tautim = 0 gives no seasonal variation; tautim = 1. gives 100% variation
c  with no steady part.
c 
      tau = 0.1
      igyres = 1
      tautim = 0.000
      wtfac = 0.      
c
      write(6,4001) tau
 4001 format (1x,/,
     1 '  Wind stress: ',/,
     2 '   stress magnitude                (default = ',f6.1,' Pa)',/,
     3 '   seasonal variation  (0. = none (default), 1.0 = 100%)',/,
     4 '   number of gyres         (0 = none, 1 (default), or 2)',/,
     5 '   (enter all three data, or default with ,,,)') 
      read (5,*) tau, tautim, igyres
      tau = tau/r0 
c                  
      do 1080 j=1,nx
      do 1080 k=1,ny
c
      yphas = 3.1415*(y(k)-y(1))/(y(ny) - y(1))
c
c  No gyres, i.e., a spatially uniform wind:
c
      if(igyres.eq.0) then
      tauy(j,k) = 0.
      taux(j,k) = tau
      end if 
c
c  One gyre, a la the subtropical gyre, say:
c
      if(igyres.eq.1) then
      tauy(j,k) = 0.      
      taux(j,k) = -tau*cos(yphas)
      end if     
c
c  Two gyres, with the boundary current squirting out into the middle:
c     
      if(igyres.eq.2) then
      tauy(j,k) = 0.
      taux(j,k) = tau*(sin(yphas) - 0.5)
      end if 
c            
 1080 continue 
c 
c  Zero out the wind stress on the west edge (left most column) where u
c  and v are set to zero normal derivative (junk). This is done solely to
c  refine the energy budget. (Depends upon the BCs, of course.)
c   
      if(basin.eq.2) go to 8323
      do 8320 k=1,ny
 8320 tauy(1,k) = 0.
      if(basin.eq.1) go to 8323
      do 8321 j=1,nx
 8321 taux(j,ny) = 0.                   
 8323 continue     
c
c
c ........  Initialize the Layer Thicknesses and Tracer Field  .............
c
c  Initialize the layer thicknesses. 
c
      do 8301 j=1,nx
      do 8301 k=1,ny
      do 8301 m=1,nz
      h(j,k,m,1) = h0(m)
 8301 continue
c
c  As an option, you can add on an initial thickness anomaly (a front or
c  an eddy) in the h field.
c
      delh = 0.
      hscale = 300.
c
      write (6,4007) delh, hscale
 4007 format (1x,/,
     1 '  Thickness anomaly of an initial eddy:',/,
     2 '   amplitude                      (default = ',f7.0,' m)',/,
     3 '   horizontal scale              (default = ',f7.0,' km)',/,
     4 '   (enter both data or default with ,,)')
      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 )
c
c  The form below gives a very sharp cutoff of the eddy amplitude; a Gaussian
c  form is also possible (but c'd out for now).
c
      eddyh = 0.
      if(rr.lt.hscale) eddyh = delh
c      
c      eddyh = delh*exp(-rr**/hscale**2)
c            
      do 1050 m=1,nz
      h(j,k,m,1) = h(j,k,m,1) + eddyh
 1050 continue 
c
c  If this is a rigid lid model, then reset the lower layer thickness
c  to make the sum of the layers equal to hbot.      
c             
      if(nz.gt.1.and.rlid.ge.1) then
      hs = 0.
      do 1058 m=1,nzm
      hs = hs + h(j,k,m,1)
 1058 continue
      h(j,k,nz,1) = hbot - hs
      end if
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 = 150.
      rr = (x(j)/1000.)**2  
      t(j,k,1,1) = h(j,k,1,1)*exp(-rr/tl**2)
      if(sqrt(rr).gt.(2.*tl)) t(j,k,1,1) = 0.  
c  
 1081 continue
c
c  
c ......................... Diffusion and Bottom Drag ....................
c
c  Define the subgridscale diffusivity, and whether it is to be applied by 
c  a Laplacian or a biharmonic operator.  Also define the bottom drag
c  coefficient; should be suitable for a gesotrophic velocity.
c   
      diff = 3.0e2
      ibhm = 0
      cdb = 0.
c
      write (6,1188) diff, cdb
 1188 format (1x,/,
     1 '  Horizontal diffusion and bottom drag:',/,
     2 '   eddy diffusion coeff.     (default =',e10.2,' m**2/sec)',/,
     3 '   biharmonic friction?         (0 = no (default), 1 = yes)',/,
     4 '   bottom drag coeffic.              (default = ',e10.2, ')',/,
     5 '   (enter all three data or default with ,,,)' )
      read (5,*) diff, ibhm, cdb
c
c
c.........................  Miscellaneous Initialization ......................
c
c
c  Specify the initial position of the floats.
c  This scheme lays them out in a circle of radius 1/4 the domain size.
c  
      if(nf.ge.1) then
      radf = (y(ny) - y(ncy))/2.
      angf = 0.
      do 9305 i=1,nf
      angf = float(i-1)*360./float(nf)
      do 9305 m=1,nz
      xf(i,m) = radf*cos(angf*3.1415/180.)
      yf(i,m) = radf*sin(angf*3.1415/180.)
 9305 continue   
      end if
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  Set up some arrays used to specify indices when horizontal differences
c  are taken. These are set up differently for solid and open 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(basin.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 open domain, v(j,1) = v(j,ny), and so:
c
      if(basin.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
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
      runto = 21.
      wrtd = 10.
      wrtm = 25.
      wrtt = 1.
c     
      if(irest.eq.1) then
      write (6,7020) time9
 7020 format (1x,//' Time up to now is  ',f9.2)
      time = time9*8.64e4
      end if       
      write (6,777) runto, wrtd, wrtm, wrtt  
  777 format (1x,/,
     1 '  Set the time to run and frequency of plotting :',/,
     2 '   time to run to               (default = ',f5.1,' days)',/,
     3 '   time between data saves      (default = ',f5.1,'  " )',/,
     4 '   time between movie frames    (default = ',f5.1,'  " )',/,
     5 '   time between terminal output (default = ',f5.1,'  " )',/,
     6 '   (enter all four or default with ,,,,)') 
      read (5,*) runto, wrtd, wrtm, wrtt
c
c  Define when and how frequently to sum for time averages:
c  tavg1 is the day on which to start the sum,
c  tavg2 is the day on which to end summing,
c  avgint is the time interval between summing.
c 
      tavg1 = 0.
      tavg2 = 0.
      avgint = 1.0
      write (6,776) 
  776 format (1x,/,
     1 '  Start and stop times for time-averaging (days):',/,
     2 '   starting day         (default = 0)',/,
     3 '   ending day           (default = 0)',/,
     4 '   (enter both data or default with ,,)')
      read (5,*) tavg1, tavg2
c
      navg = 0
      do 4507 j=1,nx
      do 4507 k=1,ny
      savg(j,k) = 0.
      do 4507 m=1,nz
      pavg(j,k,m) = 0.
      pvavg(j,k,m) = 0.
      eavg(j,k,m) = 0.
      uavg(j,k,m) = 0.
      vavg(j,k,m) = 0.
 4507 continue        
c
c  nstep is the number of time steps actually taken in this integration.
c
      nstep = 0
c
c ......................... Start a Full Time Step ........................         
c
   90 continue  
c                   
c  Reset the flag leap to zero.
c
      leap = 0
      nstep = nstep + 1
c                                            
      inow = 0
c
      if(nstep.eq.1.and.irest.eq.0) inow = 1      
c
      time = time + dt
      dtime = time/(24.*3600.) 
c
c                                                                      
c ...................... Diagnostics and Data Storage ........................
c
c  Write some basic data defining this integration to the file what.dat.
c
      if(nstep.eq.1) then
      nff = nf*nz
      write (30,2020) nx,ny,nz,nff,dx,dy,x(1),x(nx),y(1),y(ny),
     1  f0,beta,tau,rtcl,hbot
 2020 format(4i5,12e15.5)
      close (30)
      end if
c
c  Write to a disk file for time series plots, if wrtt has passed.
c
      if(timt.gt.wrtt.or.inow.eq.1) then
c
c  Compute KE and PE to check on the energy budget.
c                                            
      call kepe(h,ui,vi,eta,h0,delr0,nx,ny,nz,
     1            rgrav,basin,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 work5(j,k) = t(j,k,1,1)/h(j,k,1,1)
      call tracer(work5,nx,ny,ta,tsa,dtsa)      
c
c  Write out a time series of data from a single point.
c      
      nxw = 0.8*nx      !  this point is in a quiet region
      nyw = 0.3*ny
c
      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)
      write (10,2234) dtime, u9, v9, h9, t9, rke1, rketc,
     1 rpe, enet, wworks, fworks, bworks, ah1, ah2,
     2   ta, tsa, dtsa
 2234 format (1x, 25e15.5)
c      call flush(10)
c
c 
c
c  Compute and store the transport at three interior positions.  This is
c  useful for diagnosing the Sverdrup transport, for example.
c
c  First compute the stress curl and expected Sverdrup transport.
c
      if(nstep.eq.1) then
      nyw = int(ny*1/4)
      nx2 = nx/2
      curlt = r0*(taux(nx2,nyw+1) - taux(nx2,nyw-1))/
     1                                   (y(nyw+1) - y(nyw-1))
      svtran = curlt/(beta*r0)
      write(6,2229) curlt, beta, svtran
 2229 format (1x,'  Curlt, Beta and Svtran are ..', 3e12.3)
      end if 
c      
      mm = 0
      do 2235 m=1,3
      nyw = int(ny*1/4)
      nxw = nx - int(nx*m/4)
      u9 = 0.
      v9 = 0.
      t9 = 0.
      do 2236 kk = 1,nz 
      u9 = u9 + u(nxw,nyw,kk,1)
      v9 = v9 + v(nxw,nyw,kk,1)
      t9 = t9 + eta(nxw,nyw,kk)
 2236 continue
      mm = mm + 1
      work7(mm,1) = u9
      mm = mm + 1
      work7(mm,1) = v9
      mm = mm + 1
      work7(mm,1) = t9
 2235 continue
c 
c  Now compute the average over the eastern half of the basin.
c
      nyw = int(ny*1/2)
      nxw = nx - int(nx*2/4)
      u9 = 0.
      v9 = 0.
      t9 = 0.
      nx2s = 0
      do 2245 nx2 = nxw,nx
      nx2s = nx2s + 1
      do 2245 kk = 1,nz 
      u9 = u9 + u(nx2,nyw,kk,1)
      v9 = v9 + v(nx2,nyw,kk,1)
      t9 = t9 + eta(nx2,nyw,kk)
 2245 continue
      u9 = u9/float(nx2s)
      v9 = v9/float(nx2s)
      t9 = t9/float(nx2s)
      mm = mm + 1
      work7(mm,1) = u9
      mm = mm + 1
      work7(mm,1) = v9
      mm = mm + 1
      work7(mm,1) = t9
c 
c  Now compute the average over the western one fifth of the basin.
c
      nyw = int(ny*1/2)
      nxw = int(nx/5) 
      u9 = 0.
      v9 = 0.
      t9 = 0.
      nx2s = 0
      do 2247 nx2 = 1,nxw
      nx2s = nx2s + 1
      do 2247 kk = 1,nz 
      u9 = u9 + u(nx2,nyw,kk,1)
      v9 = v9 + v(nx2,nyw,kk,1)
      t9 = t9 + eta(nx2,nyw,kk)
 2247 continue
      u9 = u9/float(nx2s)
      v9 = v9/float(nx2s)
      t9 = t9/float(nx2s)
      mm = mm + 1
      work7(mm,1) = u9
      mm = mm + 1
      work7(mm,1) = v9
      mm = mm + 1
      work7(mm,1) = t9
c
c  Now compute the overall curl, too (use the old wtfac).
c
      curlt = wtfac*(taux(10,ny) - taux(10,1))/(y(ny) - y(1))       
c       
      write (40,2234) dtime, (work7(n,1),n=1,mm), curlt, wtfac 
c      call flush(40)
c
c  Write to the terminal occasionally to show that something is happening.
c
      nxw = nx/6 + 1 
      nyw = 2*(ny/3) + 1
      u9 = ui(nxw,nyw,1,1)
      v9 = vi(nxw,nyw,1,1)
      rn9 = eta(nxw,nyw,1)
      epse = -(wworks + fworks + bworks - enet)/(wworks + 0.001)
      write (6,653) 
     1   dtime, nits, rn9, u9, v9, rke1, rketc, rpe,
     2  enet, wworks, fworks, bworks, epse  
  653 format (1x,f7.2,i7,3f7.2,2x,3f7.2,2x,4f7.2,f8.3)
c  
c
c  At this point, check to see if the stopfile flag has been set on (=1).
c  (This serves as a way to stop this program gracefully when run on
c   a DOS machine (provided it can share files)).
c
c      open (unit=39,file='stopfile.dat',status='old',
c     1        form='formatted')
c      read (39,3939) istop
c      close (39)
c      
c      if(istop.eq.1) go to 9999
c
      if(nstep.ne.1) timt = timt - wrtt
      end if
      timt = timt + dt/8.64e4    
c 
c  Check to see if it is time to sum for time averages.
c
      if(timavg.gt.avgint) then
      if(dtime.gt.tavg1.and.dtime.lt.tavg2) then
c
c  Load the relative vorticity of the depth-integarted transport into work5, 
c  then solve for the streamfunction. The Laplacian is computed over
c  the entire domain, however, it is not used on the solid boundaries
c  where s = 0. 
c 
      do 7307 j=1,nx
      do 7307 k=1,ny
 7307 work5(j,k) = 0.
      do 7301 j=1,nx 
      jp = jpa(j)
      do 7301 k=1,ny
      km = kma(k)
      do 7302 m=1,nz 
      work5(j,k) = work5(j,k) + ( sdx*(v(jp,k,m,1) - v(j,k,m,1)) - 
     1  sdy*(u(j,k,m,1) - u(j,km,m,1)) )
 7302 continue 
 7301 continue 
c
      call sorpsi(1,work5,basin,dx,dy,nx,ny,jma,jpa,kma,kpa,work6,nnn) 
c      
      do 5001 j=1,nx
      jp = jpa(j)
      do 5001 k=1,ny
      km = kma(k)
      savg(j,k) = savg(j,k) + work6(j,k)
      do 5001 m=1,nz
      pavg(j,k,m) = pavg(j,k,m) + (p(j,k,m) + pr(j,k))/g
      eavg(j,k,m) = eavg(j,k,m) + eta(j,k,m)
      uavg(j,k,m) = uavg(j,k,m) + ui(j,k,m,1)
      vavg(j,k,m) = vavg(j,k,m) + vi(j,k,m,1)      
      pvavg(j,k,m) = pvavg(j,k,m) +  
     1  ( sdx*(vi(jp,k,m,1) - vi(j,k,m,1)) - 
     2         sdy*(ui(j,k,m,1) - ui(j,km,m,1)) )
 5001 continue
      navg = navg + 1
      end if 
c      
      timavg = timt - avgint
      end if
      timavg = timavg + dt/8.64e4 
c  
c  End of summing for the time averages.           
c      
c
c  Write out some data to disk files for general plotting.     
c 
      if(timd.gt.wrtd.or.inow.eq.1) then
c
 9998 continue
c 
      write (32,2121) dtime
 2121 format (f10.3)
c      call flush(32)            
c      
c  Load the relative vorticity of the depth-integrated transport into work5, 
c  then solve for the streamfunction. The Laplacian is computed over
c  the entire domain but 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) = 0.
      do 1302 m=1,nz 
      work5(j,k) = work5(j,k) + ( sdx*(v(jp,k,m,1) - v(j,k,m,1)) - 
     1  sdy*(u(j,k,m,1) - u(j,km,m,1)) )
 1302 continue 
 1301 continue 
c
      call sorpsi(1,work5,basin,dx,dy,nx,ny,jma,jpa,kma,kpa,work6,nnn) 
c 
c  Load pressures (divided by g) into work1, work2
c
      if(rlid.ge.1) then
      do 1308 j=1,nx
      do 1308 k=1,ny
      work1(j,k) = (pr(j,k) + p(j,k,1))/g
      work2(j,k) = (pr(j,k) + p(j,k,2))/g
 1308 continue
      else
      do 1318 j=1,nx
      do 1318 k=1,ny
      work1(j,k) =  p(j,k,1)/g
      work2(j,k) =  p(j,k,2)/g
 1318 continue
      end if      
c  
c  Write out some matrices that can be used for Matlab plotting.
c 
      do 1206 k=1,ny
c            
      write (1,1298) (eta(j,k,1),j=1,nx)
      write (2,1298) (work1(j,k),j=1,nx) 
      write (3,1298) (work5(j,k),j=1,nx) 
      write (4,1298) (work2(j,k),j=1,nx)
c
 1206 continue
 1298 format (120e15.5)
c      call flush(1)
c      call flush(2)
c      call flush(3)
c      call flush(4)
c
c  Compute the potential vorticity and absolute vorticity along a meridional
c  section and output to a disk file for plotting.
c
      do 9080 k=1,ny
      do 9080 m=1,nz
      rv = sdx*(vi(ncx+1,k,m,1) - vi(ncx,k,m,1)) - 
     1         sdy*(ui(ncx,k,m,1) - ui(ncx,kma(k),m,1))
      pvprof(k,m) = (rv + f(k))/h(ncx,k,m,1)
      avprof(k,m) = rv + f(k)
      rvprof(k,m) = rv
 9080 continue
      do 9081 k=1,ny 
      write (36,1298) (pvprof(k,m),m=1,nz), (avprof(k,m),m=1,nz), 
     1    (rvprof(k,m),m=1,nz)
 9081 continue
c      call flush(36)
c       
c      write (6,3680) 
c 3680 format (1x,'   Wrote data to disk for general plotting')
c 
      if(iblow.eq.1) go to 9999    
      if(inow.ne.1) timd = timd - wrtd 
      end if
      timd = timd + dt/8.64e4 
c
c  Write out some matrices for a matlab movie.
c     
      if(timm.gt.wrtm.or.inow.eq.1) then
c
c  Load pressures (divided by g) into work1, work2
c
      do 1346 j=1,nx
      do 1346 k=1,ny
ccccc      work1(j,k) = eta(j,k,1)
      work2(j,k) = eta(j,k,2)
 1346 continue

      write (31,2121) dtime
c      call flush(31)
c    
      do 9300 k=1,ny
      write (11,1298) (work1(j,k),j=1,nx)
 9300 continue 
c      call flush(11)
c
      if(nstep.ne.1) timm = timm - wrtm
      end if
      timm = timm + dt/8.64e4     
c
c      alf = 'at end of diag' 
c      call blowup(nstep,alf,nx,ny,nz,h0,h,ui,vi,iblow)
c             
c.................. Come Here to Begin a Trapezoidal Correction ...............
c
c  Come here to start the second pass of a leap-frog tapezoidal correction.
c  This avoids adding to the time or to nstep.
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 eta 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 
  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 layer boundaries is different from that used 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 have 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
      if(diff.gt.1.e-4) then
      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)
      end if  
c
      wwork = 0.
      fwork = 0.
      bwork = 0.
c
c  Compute the time-dependent wind stress amplitude.   Note that tautim is
c  a factor that determines the size of the seasonal cycle of the wind stress.
c
      wtfac = (1. - tautim) + tautim*sin(dtime*2.*3.14156/365.)
c
      do 2001 j=1,nx
      jp = jpa(j)
      jm = jma(j)
c
      do 2001 k=1,ny
      kp = kpa(k)
      km = kma(k)
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
      ua(j,k) = ua(j,k) + wx
      va(j,k) = va(j,k) + wy      
c
      bdx = cdb*ui(j,k,nz,1)*abs(ui(j,k,nz,1))
      bdy = cdb*vi(j,k,nz,1)*abs(vi(j,k,nz,1))
      u(j,k,nz,3) = u(j,k,nz,3) - bdx
      v(j,k,nz,3) = v(j,k,nz,3) - bdy 
c 
      ua(j,k) = ua(j,k) - bdx
      va(j,k) = va(j,k) - bdx
c      
      do 2000 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. 
c      
      udiff = ud24(j,k,m)*hxa 
      vdiff = vd24(j,k,m)*hya
c 
c  East-west momentum budget.
c 
c  Note that you must account for the C grid when evaluating the 
c  Coriolis term in the east-west momentum budget, i.e., 
c  that f is defined at the h points.
c
      fp = 0.5*(f(k) + f(kp))
      fm = 0.5*(f(k) + f(km))     
      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)) )
c 
      ulin = avf + udiff
     1  - 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) + ulin - rnonl*uadv 
      ua(j,k) = ua(j,k) + ulin - 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)) )
c     
      vlin =  vdiff - ( auf +
     1  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) + vlin - rnonl*vadv 
      va(j,k) = va(j,k) + vlin - rnonl*vadv   
c
c  A tracer budget.
c
      t(j,k,m,3) = t(j,k,m,3) + td24(j,k,m) - 
     1 ( sdx*(ctu(j,k,m) - ctu(jm,k,m)) 
     2       + sdy*(ctv(j,kp,m) - ctv(j,k,m)) )           
c
c  Evaluate 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  Sum the viscous frictional work on all the layers.
c
      fwork = fwork + r0*( udiff*ui(j,k,m,1) +  vdiff*vi(j,k,m,1) ) 
c 
 2000 continue
c 
c  Sum the work by bottom drag on the lower layer.
c
      bwork = bwork + r0*( bdx*ui(j,k,nz,1) + bdy*vi(j,k,nz,1) ) 
c  
c  Sum the work by the windstress on the surface current.
c
      wwork = wwork + r0*( wx*ui(j,k,1,1) + wy*vi(j,k,1,1) )
c 
c   
 2001 continue 
c
c ................... Apply a Rigid Lid Pressure Condition ....................
c 
      if(rlid.le.0) go to 6000
c
c  The next section applies a rigid lid condition in one of two ways. The
c  first method attempted here is to require that the depth-averaged force have 
c  zero divergence (the so-called 'pressure version'); the second and the more 
c  conventional way is to require that the depth-averaged force be computed from
c  a streamfunction, in which case it also has zero divergence. As implemented
c  here, the latter method seems to be the better one, (it works better),
c  though in principal the former should be prefferred since the correction 
c  made is smaller. 
c 
c  Apply some solid wall boundary conditions on ua, va (the 'force') before 
c  computing the curl or divergence, if this is not a channel or symmetric
c  domain. (Boundary conditions are not so clear for this.  Oh well.)
c     
      if(basin.ne.2) then       
      do 6999 j=1,nx
      va(j,1) = 0.
      va(j,ny) = 0.
 6999 continue 
      end if
      if(basin.lt.1) then
      do 6998 k=1,ny
      ua(1,k) = 0.
      ua(nx,k) = 0.
 6998 continue 
      end if
c
      ivort = 1
      if(ivort.eq.1) go to 3311
c 
c  The so-called pressure version follows next.
c    
c  Just a reminder of the grid. The rightmost column and the lowest row are
c  added on here in order to allow application of suitable boundary conditions.
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    
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
c  Compute the divergence of the net force. This will be the source term in a 
c  Poisson equation for the rigid lid pressure, pr, that 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 used in the elliptic solver for pr (apparently). 
c  
c  Try setting dx, dy to 1 for this purpose.  Also, divide by hbot.
c
      dxp = 1.
      sdxp = 1./dxp
      dyp = dy/dx
      sdyp = 1./dyp
c 
      prsc = 1.
c   
      j1 = 2
      if(basin.ge.1) j1 = 1
      do 7700 j=j1,nx
      do 7700 k=2,ny
      d2pr(j,k) = prsc*( sdxp*(ua(j,k-1) - ua(jma(j),k-1)) + 
     1   sdyp*(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 
c  Generally, it makes sense to use the old pr as the first guess for the
c  new one. However, it may be a good idea to occasionally set the first
c  guess (pr) back to zero and start from scratch. (Sounds good, but this
c  didn't help).
c
      go to 8891
      jsetpr = 2000
      if(mod(nstep,jsetpr).eq.1) then
      do 4780 j=1,nxp
      do 4780 k=1,nyp
      prp(j,k) = 0.
 4780 continue
      write (6,4781) nstep
 4781 format (1x, '  I set pr to zero; nstep =', i4)
      end if
 8891 continue
c
      eps = 1.e-5
      eps = 0.5e-5
      if(dtime.gt.25.) eps = 0.5e-5            
      call sorpr
     1  (nstep,nits,eps,d2pr,basin,dxp,dyp,nxp,nyp,
     2        jmap,jpap,kmap,kpap,prp) 
c
c  prm is the isopycnal displacement that corresponds to the rigid
c  lid pressure (used to plot the rigid lid pressure).
c    
      do 7760 j=1,nx
      do 7760 k=1,ny
      prm(j,k) = prp(j,k+1)/(g*rtcl*prsc)
 7760 continue
c
c  Now add on the rigid lid pressure gradient so that the depth-averaged
c  force will have zero divergence. 
c
      twoh = 2.*hbot
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 by +1 to account for the 
c  extra (lower) row added onto pr.
c      
      prx(j,k) = sdxp*(prp(jp,k+1) - prp(j,k+1))/prsc
      pry(j,k) = sdyp*(prp(j,k+1) - prp(j,k))/prsc
c                     
      do 7750 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) -
     1   prx(j,k)*(h(j,k,m,1) + h(jp,k,m,1))/twoh 
      v(j,k,m,3) = v(j,k,m,3) - 
     1   pry(j,k)*(h(j,k,m,1) + h(j,km,m,1))/twoh
c
c  Store some things for diagnostic use.
c
      uapr(j,k) = ua(j,k) - prx(j,k)
      vapr(j,k) = va(j,k) - pry(j,k) 
 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
      prx(j,k) = sdxp*(prp(j+1,k+1) - prp(j,k+1))
      do 7751 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) - 
     1  prx(j,k)*(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
      pry(j,k) = sdyp*(prp(j,k+1) - prp(j,k))
      do 7752 m=1,nz
      v(j,k,m,3) = v(j,k,m,3) - 
     1  pry(j,k)*(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
      go to 6000                      
c  
c 
 3311 continue
c
c  In this next version we attempt to apply the rigid lid constraint by
c  computing the net force from a streamfunction-like field. This has
c  an analogy with barotropic vorticity conservation, but the most direct
c  way to think of this is that the net force has to be divergence-free
c  if the net transport is intended to always be divergence-free.  
c
c  Compute the curl of the net force saved in ua,va. This will be inserted
c  into an elliptic solver, and then differentiated to get a net force that
c  should have not diveregence.  
c
      do 6060 j=2,nxm
      jp = jpa(j)
      do 6060 k=2,nym
      km = kma(k)
      del2st(j,k) = sdx*(va(jp,k) - va(j,k)) - sdy*(ua(j,k) - ua(j,km)) 
 6060 continue
c
c  Solve a Poisson equation for a streamfunction of the total force, st. This 
c  subroutine assumes that st is zero on the boundaries.
c
c  To try to avoid instability, set the first guess field, st, to zero 
c  every nto time steps.
c
c      nto = 300
c      if(mod(nstep,nto).eq.-1) then
c      do 7328 j=1,nx
c      do 7328 k=1,ny
c      st(j,k) = 0.
c 7328 continue
c      write (6,7329)
c 7329 format (1x, ' I set st to zero, just like you said')
c      end if
c
c      if(leap.eq.2) then
c      do 1420 j=1,nx
c      do 1420 k=1,ny
c      del2st(j,k) = 0.5*(del2st(j,k))
c 1420 continue
c      end if
c 
      call sorpsi(nstep,del2st,basin,dx,dy,nx,ny,jma,jpa,
     1 kma,kpa,st,nits)
c
c  Require that the momentum balance be consistent with this
c  streamfunction.
c 
c
c  special diagnostic......................................................
c
      do 5400 j=1,nx
      do 5400 k=1,ny
      work1(j,k) = u(j,k,1,3)*htcl
      work2(j,k) = v(j,k,1,3)*htcl
 5400 continue
c      
c      do 6061 j=2,nx   ***********************************************
      do 6061 j=1,nx
      jm = jma(j)
      jp = jpa(j) 
c      do 6061 k=1,nym    ********************************************
      do 6061 k=1,ny
      km = kma(k)
      kp = kpa(k)
      dsdx = sdx*(st(j,k) - st(jm,k))                                          
      dsdy = sdy*(st(j,kp) - st(j,k))
      dprdx(j,k) = (ua(j,k) + dsdy)/hbot
      dprdy(j,k) = (va(j,k) - dsdx)/hbot
      do 6061 m=1,nz
      u(j,k,m,3) = u(j,k,m,3) - dprdx(j,k)*(h(j,k,m,1) + h(jp,k,m,1))/2.
      v(j,k,m,3) = v(j,k,m,3) - dprdy(j,k)*(h(j,k,m,1) + h(j,km,m,1))/2.
 6061 continue
c
c  Apply BCS that are appropriate to a closed basin.
c
      do 9703 m=1,nz
      do 9701 j=1,nx
      u(j,ny,m,3) = u(j,nym,m,3)
      h(j,ny,m,3) = h(j,nym,m,3)
 9701 continue
      do 9702 k=1,ny
      v(1,k,m,3) = v(2,k,m,3)
      h(1,k,m,3) = h(2,k,m,3)      
 9702 continue
 9703 continue
c 
      do 5401 j=1,nx
      do 5401 k=1,ny
      work3(j,k) = u(j,k,1,3)*htcl
      work4(j,k) = v(j,k,1,3)*htcl
 5401 continue
c
c  special diagnostic......................................................
c  
c
c  Compute the rigid lid pressure field from the pressure derivatives  
c  (this is used for diagnostics only).
c 
      pr(2,1) = 0.
      do 6004 j=2,nx
      do 6003 k=2,ny
      pr(j,k) = pr(j,k-1) + dy*dprdy(j,k)
 6003 continue 
      pr(j,1) = pr(j-1,1) + dx*dprdx(j-1,1)
 6004 continue 
      do 6005 k=1,ny
      pr(1,k) = pr(2,k)
 6005 continue 
c      
 6000 continue
c
c  If you really want to insure a rigid lid then require that the layer
c  thickness tendencies must sum to zero (or set ijk = 0 to skip this).  
c     
      ijk = 1
      if(rlid.eq.2.and.ijk.eq.1) then
      do 4680 j=1,nx
      do 4680 k=1,ny 
      ws = 0.
      do 4681 m=1,nz
 4681 ws = ws + h(j,k,m,3)
      h(j,k,nz,3) = h(j,k,nz,3) - ws 
 4680 continue               
      end if
c    
c  The end of all of the rigid lid calculations.       
c
c  Apply radiation boundary conditions, if required.
c 
      if(basin.ge.1.and.openbc.eq.1) then
c
c  Apply a radiation BC to the east and west sides, assuming that cph is 
c  the relevent speed. 
c           
      do 3500 k=1,ny 
      do 3500 m=1,nz
c
      u(1,k,m,3) = cph*( u(2,k,m,1) - u(1,k,m,1) )/dx
      u(nx,k,m,3) = -cph*( u(nx,k,m,1) - u(nxm,k,m,1) )/dx
      v(1,k,m,3) = cph*( v(2,k,m,1) - v(1,k,m,1) )/dx
      v(nx,k,m,3) = -cph*( v(nx,k,m,1) - v(nxm,k,m,1) )/dx      
      h(1,k,m,3) = cph*( h(2,k,m,1) - h(1,k,m,1) )/dx
      h(nx,k,m,3) = -cph*( h(nx,k,m,1) - h(nxm,k,m,1) )/dx
 3500 continue
      end if
c
      if(basin.eq.2.and.openbc.eq.1) then
c
c  Apply a radiation BC to the north and south sides, if required.
c       
      do 3501 j=1,nx
      do 3501 m=1,nz 
c
      u(j,1,m,3) = cph*( u(j,2,m,1) - u(j,1,m,1) )/dy
      u(j,ny,m,3) = -cph*( u(j,ny,m,1) - u(j,nym,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,ny,m,3) = -cph*( h(j,ny,m,1) - h(j,nym,m,1) )/dy 
 3501 continue
      end if
c 
                                    
c ............................. Time Step Ahead .............................
c
c  Step ahead, noting whether this particular 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  The frequency of the trapezoidal correction is set by n7 and n8, i.e., do
c   a correction on n8 consecutive steps every n7 steps.  Reasonable
c   values seem to be n7 = 15, n8 = 2.  To avoid doing any trapezoidal 
c   correction steps, you could set n8 < 0.
c
      n7 = 25
      n8 = 3
c
      if(mod(nstep,n7).ge.(n7-n8)) leap = leap + 1 
c 
c  A straight leap-frog step.
c    
      if(leap.eq.0) then
c      
      do 300 j=1,nx
      do 300 k=1,ny
c
      ua(j,k) = 0.
      va(j,k) = 0.
c      
      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(basin.eq.0) rnxy = float((nx-1)*(ny-1))
      wworks = wworks + dt*wwork/(rnxy*r0)
      fworks = fworks + dt*fwork/(rnxy*r0)
      bworks = bworks - dt*bwork/(rnxy*r0)
c      
      go to 303
      end if
c
c  A leap-frog that will be followed by a trapezoidal correction.
c  (So in this case you do not zero out the accelerations).
c
      if(leap.eq.1) then 
c      
      do 301 j=1,nx
      do 301 k=1,ny            
      do 301 m=1,nz 
c 
      ua(j,k) = 0.
      va(j,k) = 0.
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)
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)
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)
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)
c
  301 continue 
      go to 303
      end if
c
c  A trapezoidal correction.
c
      if(leap.eq.2) then
c      
      do 302 j=1,nx
      do 302 k=1,ny
c
      ua(j,k) = 0.
      va(j,k) = 0.
c      
      do 302 m=1,nz 
c      
      u(j,k,m,1) = u(j,k,m,2) + halfdt*u(j,k,m,3)
      u(j,k,m,3) = 0.
c       
      v(j,k,m,1) = v(j,k,m,2) + halfdt*v(j,k,m,3)
      v(j,k,m,3) = 0.
c      
      h(j,k,m,1) = h(j,k,m,2) + halfdt*h(j,k,m,3)
      h(j,k,m,3) = 0.
c 
      t(j,k,m,1) = t(j,k,m,2) + 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(basin.eq.0) rnxy = float((nx-1)*(ny-1))
      wworks = wworks + dt*wwork/(rnxy*r0) 
      fworks = fworks + dt*fwork/(rnxy*r0)
      bworks = bworks - dt*bwork/(rnxy*r0)      
c 
c           
      end if  
c
  303 continue
c
c      alf = 'after dt, before bcs'
c      call blowup(nstep,alf,nx,ny,nz,h0,h,ui,vi,iblow)

c  
c
c.......................... Apply Boundary Conditions .........................        
c        

c  Apply BCs to all solid boundaries. Recall that the BCs are specified 
c  via the following flags:
c
c  The basin configuration is set by:
c  basin = 0 for solid walls all the way around,
c  basin = 1 for open BCs on east and west sides, (this makes a channel),
c  basin = 2 for open BCs all the way around. 
c  
c  The tangential velocity on side boundaries is set by:
c  sidebc = 0 for no-slip BCs on the sidewalls of a basin, or, 
c  sidebc = 1 for free-slip BCs on the sidewalls of a basin (the default).
c
      if(basin.eq.2) go to 3100
c             
      if(basin.eq.1) go to 3003        
c
c  East and west sides, required for a fully enclosed basin only.
c
c  Apply zero normal flow bcs to the east and west sides.
c      
      do 3000 k=1,ny 
      do 3000 m=1,nz
      u(1,k,m,1) = 0. 
      u(nx,k,m,1) = 0. 
c      
c  Apply no slip BCs to v if so required.
c
      if(sidebc.eq.0) then 
      v(2,k,m,1) = 0.
      v(nx,k,m,1) = 0. 
      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)
c      h(1,k,m,3) = h(2,k,m,3)
c      t(1,k,m,3) = t(2,k,m,3)
c      v(1,k,m,3) = v(2,k,m,3)
      
 3000 continue
c
c  Done with the east and west sides.
c
 3003 continue
c
c  South and north sides, required for a channel or a basin.
c
c  Apply zero normal flow.
c
      do 3001 j=1,nx
      do 3001 m=1,nz 
      v(j,1,m,1) = 0.
      v(j,ny,m,1) = 0. 
c      
c  Apply no slip BCs to u if so required.      
c
      if(sidebc.eq.0) then
      u(j,1,m,1) = 0.
      u(j,nym,m,1) = 0.
      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)
c      h(j,ny,m,3) = h(j,nym,m,3)
c      u(j,ny,m,3) = u(j,nym,m,3)
c      t(j,ny,m,3) = t(j,nym,m,3)
c       
 3001 continue
c 
c  Done with the north and south sides.
c 
c  Now do h and eta in the northwest corner of a closed basin.
c
      if(basin.eq.0) then
      do 3004 m=1,nz
      h(1,ny,m,1) = (h(1,nym,m,1) + h(2,ny,m,1))/2.
 3004 continue
      end if
c 
 3100 continue
c 
c  End of the boundary conditions.
c
c      alf = 'just after bcs'
c      call blowup(nstep,alf,nx,ny,nz,h0,h,ui,vi,iblow)

c
c......................... Control Time Stepping .............................
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 a Time Step ...............
c
c  Step ahead the stupid floats.
c                        
      if(nf.gt.1.and.nstep.gt.1) then
      nof = 0
      if(nof.ne.1) then
c 
c  Compute a time-averaged velocity that is centered on the h points and that
c  knows about the boundary conditions. This new velocity field was found 
c  to be necessary to keep the damned floats from going off of the grid 
c  and causing a blowup of the entire model. 
c
      do 7323 m=1,nz
      do 7320 j=1,nx
      jm = jma(j)
      do 7320 k=1,ny
      kp = kpa(k)
      ufl(j,k,m) = 0.25*( ui(j,k,m,1) + ui(jm,k,m,1) + 
     1                  ui(j,k,m,2) + ui(jm,k,m,2) )
      vfl(j,k,m) = 0.25*( vi(j,k,m,1) + vi(j,kp,m,1) +
     1                  vi(j,k,m,2) + vi(j,kp,m,2) )
 7320 continue
c
c  Apply BCs to this velocity, if required.
c
      if(basin.eq.0) then
      do 7322 k=1,ny
 7322 ufl(nx,k,m) = 0.
      if(basin.eq.1) then 
      do 7321 j=1,nx
 7321 vfl(j,1,m) = 0.
      end if 
      end if
 7323 continue
c     
      do 7300 m=1,nz
      do 7300 i=1,nf
      xx = xf(i,m)
      yy = yf(i,m)
      call floatxy(xx,yy,ufl,vfl,x,y,nx,ny,nz,m,dt,xxn,yyn,uf,vf)
      xf(i,m) = xxn
      yf(i,m) = yyn 
 7300 continue      
c
c  Store the float data on unit 35 if sufficient time has passed.
c
      if(mod(nstep,50).eq.2) then
      madd = nf*nz
      madd2 = 2*madd
      iii = 0
      do 7303 m=1,nz
      do 7303 i=1,nf
      iii = iii + 1
      xyfk(iii) = xf(i,m)/1000.
      xyfk(iii+madd) = yf(i,m)/1000.
 7303 continue     
      write (35,7362) (xyfk(k),k=1,madd2)
 7362 format (50f10.2)
c      call flush(35)
c
      end if
      end if
c
      end if
c
c  End of float tracking.
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,///,1x,
     1  ' You reached the end of the integration. ',//)
      go to 9999
      end if
      if(nstep.eq.1.and.irest.ne.1) go to 310 
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(nstep.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 
c  Normalize the time averages, and write them out to unit 12.
c
      if(navg.gt.2) then
      rnavg = float(navg)
      do 8710 j=1,nx
      do 8710 k=1,ny
      savg(j,k) = savg(j,k)/rnavg
      do 8710 m=1,nz      
      pavg(j,k,m) = pavg(j,k,m)/rnavg
      pvavg(j,k,m) = pvavg(j,k,m)/rnavg
      eavg(j,k,m) = eavg(j,k,m)/rnavg
      uavg(j,k,m) = uavg(j,k,m)/rnavg
      vavg(j,k,m) = vavg(j,k,m)/rnavg      
 8710 continue
      do 8711 k=1,ny
 8711 write (12,1298) (savg(j,k),j=1,nx)
 8719 continue
      do 8712 m=1,nz
      do 8712 k=1,ny
 8712 write (12,1298) (pavg(j,k,m),j=1,nx)
      do 8713 m=1,nz
      do 8713 k=1,ny
 8713 write (12,1298) (pvavg(j,k,m),j=1,nx)
      do 8714 m=1,nz
      do 8714 k=1,ny
 8714 write (12,1298) (eavg(j,k,m),j=1,nx)
      do 8715 m=1,nz
      do 8715 k=1,ny
 8715 write (12,1298) (uavg(j,k,m),j=1,nx)
      do 8716 m=1,nz
      do 8716 k=1,ny
 8716 write (12,1298) (vavg(j,k,m),j=1,nx) 
c      call flush(12)      
      end if
c 
c 
 7053 continue                  
      irest = 1
      write (6,7050)
 7050 format (1x,' Do you want to write a restart file? (1=yes)')
      read (5,7056,err=7053) irest
 7056 format(i3)
      if(irest.eq.1) then
      write (6,7052)
 7052 format (1x,' Enter the name of the 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 sidebc,basin,rgrav,rlid,ibhm,cdb,
     2 rnonl,dt,dx,dy,diff,hbot,h0,delr0,rtcl,x,y,f,f0,beta,tau,gprime,
     3 wworks,fworks,bworks,xf,yf,taux,tauy,u,v,h,t,pr,ui,vi 
      end if
c
      stop 
      end
c
      subroutine blowup(nstep,alf,nx,ny,nz,h0,h,ui,vi,iblow)
c      
c  
c  Check for blowups.
c
      dimension h0(nz),h(nx,ny,nz,3),ui(nx,ny,nz,3),vi(nx,ny,nz,3)
      character*60 alf
c
      iblow = 0
c      
      ublow = 0.2
c     
      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
  800 continue 
c
c  
      if(umax.gt.ublow) then
      write (6,803) jvmax, kvmax, mvmax, 
     1 ui(jvmax,kvmax,mvmax,1), vi(jvmax,kvmax,mvmax,1),
     1    h(jvmax,kvmax,mvmax,1) 
  803 format (1x,' ublow; j, k, m, ui, vi, h are ', 3i5, 3f10.3)
      iblow = 1
      end if
c
c  This limit on eta could be too small.
c
      etablow = 0.8*h0(1)
c 
      if(etamax.gt.etablow) then
      write (6,804) jemax, kemax, memax,
     1    h(jemax,kemax,memax,1), ui(jemax,kemax,memax,1), 
     2  vi(jemax,kemax,memax,1)
  804 format (1x,'etablow;  j, k, m, h, ui, vi', 3i6, 4f10.3)
      iblow = 1
      end if
c
      write (6,1) nstep, alf
    1 format (1x, 1x,i6, '  completed blowup check', a60)    
c 
c
c  End of the blowup check.
c 
      return
      end