c     Main loop goes here
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine runSim()

      implicit none

c     User-modifiable parameters come first
c     Regular variables come later

      integer nX,nY,nMax
      real*8 eps, tStep, 
     *     xMin, xMax,
     *     yMin, yMax

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     User-modifiable parameters are eitehr here
c     or in the file iniVals.f

      integer bufSize           !number of past timesteps to buffer
c     The buffer of past non-linear terms for
c     the adams bashforth method works like a
c     "ring buffer"  Essentially the current
c     head of the buffer is tracked by an
c     integer, "bufIdx", which is taken modulo
c     bufSize.  Previous entries are bufIdx-1 mod bufsize, etc.
c     Obviously this works better with an index that starts at zero

      parameter (nX=400)
      parameter (nY=100)
      parameter (nMax=400)
      parameter (bufSize=2)     ! Size of ring buffer
      parameter (xMin=0.d0)
      parameter (xMax=4.d0)
      parameter (yMin=0.d0)
      parameter (yMax=1.d0)
      parameter (eps=1.d-20)

c     xStep approx 1/100 = 0.01
c     yStep approx the same
c     xStep*xStep approx 0.0001
c     Set tStep to a quarter of that (2.5d-5)
      parameter (tStep = 2.5d-5) ! restore to 2.5d-5
c     Thats a really small timestep!

      integer maxIter, printFreq
c      parameter(maxIter=1000000) ! Thats a lotta iters!
      parameter(maxIter=9000)   ! For testing
      parameter(printFreq=5)

c     End of user-modifiable stuff
ccccccccccccccccccccccccccccccccccccc

      real*8 u, uv, uu, w, psi
      real*8 dwdtNL, dwdtLin
      real*8 x,y

      dimension dwdtNL(0:nMax+1, 0:nMax+1, 
     *     0:bufSize-1)         ! non-linear part of deriv
      dimension dwdtLin(0:nMax+1, 0:nMax+1, 
     *     0:bufSize-1)         ! linear part of deriv
      dimension x(0:nX+1), y(0:nY+1)

      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1, 0:nMax+1)
      dimension u(2,0:nMax+1,0:nMax+1)

      integer it, bufIdx, nPrintCnt
c     bufIdx is the index into past deriv buffers


      real*8 currTime
      real*8 vMax,vMin

      call step(x, xMin, xMax, nx, .TRUE.)
      call step(y, yMin, yMax, ny, .FALSE.)

      currTime = 0.d0

      nPrintCnt = 0

c     Set initial conditions
c     
      call initializeFlow(
     *     u,  nX, nY, 
     *     nMax, x, 
     *     y)

      call printUVresult(x,y,u,nx,ny,nmax,
     *     currTime,nPrintCnt,
     *     .FALSE., .FALSE.)

      nPrintCnt = nPrintCnt+1


      call wFromU(w, u, nX, nY, 
     *     nMax, x(1)-x(0), y(1)-y(0), eps)


      call safePsiFromW(psi, 
     *     w, nX, nY, nMax,
     *     x(1)-x(0), y(1)-y(0), eps)

      call printPsi(x,y,psi,nx,ny,nmax,
     *     currTime,0,
     *     .FALSE., .FALSE.)


      call uFromPsi(u, psi, nX, nY, 
     *     nMax, x(1)-x(0), y(1)-y(0), eps)


      call printUVresult(x,y,u,nx,ny,nmax,
     *     currTime,nPrintCnt,
     *     .FALSE., .FALSE.)

      call printPsi(x,y,psi,nx,ny,nmax,
     *     currTime,nPrintCnt,
     *     .FALSE., .FALSE.)

      nPrintCnt = nPrintCnt+1
      

      do it=0,bufSize
         bufIdx = mod(it,bufSize)

         write (*,*) "Starting an RK step iter ",it

         call DoRkStep(x,y,psi,w,
     *        dwdtLin, dwdtNL, tStep,
     *        nX,nY,nMax,
     *        bufIdx, bufSize, eps)


         currTime = currTime+tStep
      enddo


            call safePsiFromW(psi, 
     *           w, nX, nY, nMax,
     *           x(1)-x(0), y(1)-y(0), eps)


            call uFromPsi(u, psi, nX, nY, 
     *           nMax, x(1)-x(0), y(1)-y(0), eps)


            call printUVresult(x,y,u,nx,ny,nmax,
     *           currTime,nPrintCnt,
     *           .FALSE., .FALSE.)
            call printPsi(x,y,psi,nx,ny,nmax,
     *           currTime,nPrintCnt,
     *           .FALSE., .FALSE.)
            
            nPrintCnt = nPrintCnt+1

c     pause "Done RK Steps"

      do it=bufSize+1,maxIter
         bufIdx = mod(it,bufSize)

c     SIAB stands for
c     Semi Implicit Adams Bashforth
c
c     I don't know if its a standard acronym
c     but DoSemiImplicitAdamsBashforthStep
c     is an awfully long function name

         call DoSIABStep(x,y,psi,w,
     *        dwdtLin, dwdtNL, tStep,
     *        nX,nY,nMax,
     *        bufIdx, bufSize, eps)


         call getUvMinMax(u,nX,nY,nMax,
     *        vMin, vMax,2)

         write (*,*) "on iteration ", it,
     *        " Min/max are ", vMin, vMax

         if(mod(it,printFreq).EQ.0) then

            call safePsiFromW(psi, 
     *           w, nX, nY, nMax,
     *           x(1)-x(0), y(1)-y(0), eps)


            call uFromPsi(u, psi, nX, nY, 
     *           nMax, x(1)-x(0), y(1)-y(0), eps)


            call printUVresult(x,y,u,nx,ny,nmax,
     *           currTime,nPrintCnt,
     *           .FALSE., .FALSE.)
            call printPsi(x,y,psi,nx,ny,nmax,
     *           currTime,nPrintCnt,
     *           .FALSE., .FALSE.)
            
            nPrintCnt = nPrintCnt+1

         endif

         currTime = currTime+tStep
      enddo
      

      return
      end


ccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Steps the system once using RK4
c     Used to jump-start semi-implicit Adams Bashforth
c
c     stands for DoRkStep as in
c     Do Runge Kutta step
c     not as in "DORKstep"
c     
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine DoRkStep(x,y,psi,w,
     *     linDrvs, nlDrvs, dt,
     *     nX,nY,nMax,
     *     nBufIdx, nBufSize, eps)
      implicit none
c     linDrvs means "LINear components of DeRiVativeS"
c     nlDrvs means "NonLinear components of DeRiVativeS"
      
      real*8 x,y,psi,w
      real*8 linDrvs, nlDrvs
      real*8 dt
      real*8 tmpArr
      real*8 eps
      
      integer nX, nY, nMax,
     *     nBufIdx, nBufSize
      
      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1,0:nMax+1)
      dimension linDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
      dimension nlDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
      dimension x(0:nX+1), y(0:nY+1)

      dimension tmpArr(0:nMax+1,0:nMax+1)
      
      real*8 deltaX, deltaY 

      deltaX = x(1)-x(0)
      deltaY = y(1)-y(0)


      call psiFromW(psi, 
     *     w, tmpArr, nX, nY, nMax,
     *     deltaX, deltaY, eps)


      call GetLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMax,
     *     linDrvs(0,0,nBufIdx))

      call GetNonLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMax,
     *     nlDrvs(0,0,nBufIdx))

c     For now we're just going to Euler-step these ones
c     Hopefully Euler-stepping for the first couple
c     steps is Okay for jump-starting Adams Bashforth

      call makeLinComb(
     *     linDrvs(0,0,nBufIdx),
     *     nlDrvs(0,0,nBufIdx),tmpArr,
     *     nX,nY,nMax,nMax,nMax,
     *     dt,dt, .FALSE.)

      call makeLinComb(
     *     w,
     *     tmpArr,w,
     *     nX,nY,nMax,nMax,nMax,
     *     1.d0,1.d0, .FALSE.)


      return
      end


ccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Steps the system using
c     semi-implicit Adams Bashforth
c     
cccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine DoSIABStep(x,y,psi,w,
     *     linDrvs, nlDrvs, dt,
     *     nX,nY,nMax,
     *     nBufIdx, nBufSize, eps)
      implicit none
c     linDrvs means "LINear components of DeRiVativeS"
c     nlDrvs means "NonLinear components of DeRiVativeS"
      
      real*8 x,y,psi,w
      real*8 linDrvs, nlDrvs
      real*8 dt
      real*8 abCoeff
      real*8 tmpArr
      real*8 diag, diagx, diagy
      real*8 eps, deltax, deltay
      real*8 getViscosity
      
      dimension abCoeff(0:2)    !Adams Bashforth coeffs
      
      integer nX, nY, nMax,
     *     nBufIdx, nBufSize,
     *     nBTmp                ! nBTmp is used in place
c     of nBufIdx in places where we want to subtract
c     and mod
      
      dimension psi(0:nMax+1,0:nMax+1)
      dimension w(0:nMax+1,0:nMax+1)
      dimension linDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
      dimension nlDrvs(0:nMax+1,0:nMax+1,0:nBufSize-1)
      dimension tmpArr(0:nMax+1, 0:nMax+1)
      dimension x(0:nX+1), y(0:nY+1)

      deltaX = x(1)-x(0)
      deltaY = y(1)-y(0)
      call psiFromW(psi, 
     *     w, tmpArr, nX, nY, nMax,
     *     deltaX, deltaY, eps)

      nBTmp = nBufIdx+nBufSize

      abCoeff(0) = 23.d0/12.d0
      abCoeff(1) = -16.d0/12.d0
      abCoeff(2) = 5.d0/12.d0
      
      call GetLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMax,
     *     linDrvs(0,0,nBufIdx))

      call GetNonLinTerm(
     *     psi, w, x,y,
     *     nX, nY, nMax, 
     *     nMax,
     *     nlDrvs(0,0,nBufIdx))

      call mul_and_copy(w, 
     *     tmpArr, 1.d0, nX, nY,
     *     nMax, nMax,
     *     .TRUE.)
c     Now the boundaries on tmpArr are correct
c     also tmpArr =w
      
      call makeLinComb(
     *     tmpArr,
     *     nlDrvs(0,0,mod(nBTmp,nBufSize)),
     *     tmpArr,
     *     nX,nY,nMax,nMax,nMax,
     *     1.d0,
     *     dt*abCoeff(0), .FALSE.)

      call makeLinComb(
     *     tmpArr,
     *     nlDrvs(0,0,mod(nBTmp-1,nBufSize)),
     *     tmpArr,
     *     nX,nY,nMax,nMax,nMax,
     *     1.d0,
     *     dt*abCoeff(1), .FALSE.)

      call makeLinComb(
     *     tmpArr,
     *     nlDrvs(0,0,mod(nBTmp-2,nBufSize)),
     *     tmpArr,
     *     nX,nY,nMax,nMax,nMax,
     *     1.d0,
     *     dt*abCoeff(2), .FALSE.)

c     Now we have the non-linear term,
c     it is time for the semi-implicit part

c     First add half deltaT times
c     the linear term in the derivative

      call makeLinComb(
     *     tmpArr,
     *     linDrvs(0,0,mod(nBTmp,nBufSize)),
     *     tmpArr,
     *     nX,nY,nMax,nMax,nMax,
     *     1.d0,
     *     dt*5.d-1, .FALSE.)

c     Now we have to use the conjugate-gradient
c     solver to find the next w

      diagx = -dt*5.d-1/(deltaX*deltaX)
      diagy = -dt*5.d-1/(deltaY*deltaY)

      diagx = diagx*getViscosity()
      diagy = diagy*getViscosity()

      diag = 1.d0-2.d0*diagx-2.d0*diagy

      write (*,*) "Calling conjgrad "
      write (*,*) "dt is ", dt

      call conjgrad_2d_lat(diag,diagx,diagy,
     *     w, tmpArr,
     *     nX, nY, nMax, eps, .FALSE.)

c     subroutine conjgrad_2d_lat(diag,xdiag,ydiag,x0,b,
c     *     nX, nY, nMax, eps, bYBoundaries)

      return 
      end