C***********************************************************************
C***********************************************************************
      SUBROUTINE SACO_RUN(CASCO_ITER) 
C-----------------------------------------------------------------------
C 
C The main module for running saco within the saco-moody interation
C
C-----------------------------------------------------------------------
C RELEASE 2.0A (May 19, 1998)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Hacked from QUODDY4 at Dartmouth College by:
C   Charles Hannah 
C   Daniel R. Lynch
C   Christopher E. Naimie
C
C QUODDY4 was Developed at Dartmouth College by:
C   Daniel R. Lynch
C   Justin T.C.Ip
C   Christopher E. Naimie
C   Francisco E. Werner
C
C
C
C   *OLD arrays are initialized with dependent variable data at one
C   time step prior to the beginning of the simulation  
C
C   ZETAOLD(I) - elevation 
C   UOLD(I),VOLD(I) - x,y components of vertically averaged
C               horizontal velocity
C
C   *MID arrays are initialized with dependent variable data at the
C   beginning of the simulation  
C
C   ZETAMID(I) - elevation 
C   UMID(I),VMID(I) - x,y components of vertically averaged
C                     horizontal velocity
C   ZMID(I,J) - 3-D nodal grid 
C   UZMID(I,J),VZMID(I,J) - x,y components of  velocity 
C 
C
C +++++++++++++++++++++++++++++++++++++++++
C + PREPARATION FOR TIME STEPPING SECTION +
C +++++++++++++++++++++++++++++++++++++++++
C
C
C MAIN PROGRAM:
C   
C   KDMID - Set to initial Gregorian Day
C   SECMID - Set to starting time on initial Gregorian Day
C   ITER - Set to 0
C
C +++++++++++++++++++++++++
C + TIME STEPPNIG SECTION +
C +++++++++++++++++++++++++
C
C MAIN PROGRAM:
C
C   SECNEW - Incremented by DELT
C   KDNEW - Incremented as necessary (which also affects SECNEW)
C   ITER - Incremented by 1
C
C
C MAIN PROGRAM:
C
C
C VARIOUS SUBROUTINES (AS NECESSARY): All at time level K:
C
C   AK(I) - Linearized bottom stress coefficients
C
C SUBROUTINE ELEVATIONQ4:  calculation centered at time level K
C
C   ZETANEW(I) - free surface elevation at time level K+1
C
C SUBROUTINE VERTICALS5: calculation centered at time level K+1/2
C
C  UZNEW(I,J),VZNEW(I,J) - x,y components of velocity at time level K+1
C
C MAIN PROGRAM:
C
C  *MID=>*OLD
C  *NEW=>*MID
C
C SUBROUTINES VERTAVG: At the end of the current time step:
C
C  UMID(I,J),VMID(I,J) - vertically averaged horizontal velocity
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C ARRAY DECLARATION SECTION
C
C Objectives: Declare all arrays required by the main program.
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'CASCO.DIM'
      include 'mesh.h'
      include 'saco.h'
      include 'startstop.h'
      include 'statarrays.h'
      include 'friction.h'
      include 'constants.h'
      include 'misc.h'

      integer casco_iter
C
C array declaration
C
C moved to saco.h so that moody can see them 
c     REAL UZNEW(NNDIM,NNVDIM), VZNEW(NNDIM,NNVDIM)
c     REAL ZETANEW(NNDIM)
CGH are unew and vnew used ?? 
c     REAL UNEW(NNDIM),VNEW(NNDIM)
C
C local stuff
      REAL TEMP1(NNDIM),TEMP2(NNDIM)
      INTEGER kdnew, kdmid
      REAL  secnew,secmid
      INTEGER IDAYMID, IYEARMID, IMNTHMID 
      INTEGER ITER 
      INTEGER I,J  

      write(*,*)'SACO_RUN: casco_iter ', casco_iter
      WRITE(2,*)'     Sparse factor is',NFTR
C
 
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C PREPARATION FOR TIME STEPPING SECTION
C
C Basic Initialization should already be done by STARTUP
C Initialization of bottom friction, vertical eddy viscosity
C should have been done by get_friction. 
C
C What is left is to 1) initialize basic arrays 
C                    2) call BCS and OUTPUT to initialize them
C                    3) start the run 
C
C CALL BCS2 and OUTPUTS2 to facilitate initialization of local 
C information.  ATMXNEW and ATMYNEW are used as dummy arrays here.  They will
C be zeroed again before they are used for computation herein.
C
C  Initialization 
        do i = 1, nn
           zetaold(i) = 0.0
           zetamid(i) = 0.0
           zetanew(i) = 0.0
           umid(i) = 0.
           vmid(i) = 0.
        enddo
        do j = 1, nnv
           do i = 1, nn
              uzmid(i,j) = 0.0
              vzmid(i,j) = 0.0
              uznew(i,j) = 0.0
              vznew(i,j) = 0.0
           enddo
        enddo

C 
      ITER = 0
      CALL BCS2(KD0,SEC0,ITER,DEG,TEMP1,TEMP2,0)

C Compute the difference between the model and the observations.
C Write the error field to arrays defined in observes.h
C
      CALL SACO_UVERROR(ITER,IN,NNV,UZMID,VZMID)
      CALL SACO_EERROR(ITER,IN,zetamid)
C
      CALL OUTPUTS2(ITER,KD0,SEC0,DELT,MESHNAME,FILINQ,FILICQ,
     & NN,NNV,NE,X,Y,IN,HDOWN, AK, ENZM,  
     & ZETAMID,UMID,VMID,ZMID,UZMID,VZMID,ZETAOLD,UOLD,VOLD,casco_iter)

C
C
C Finally, initialize timing parameters
C
      KDMID=KD0
      SECMID=SEC0
C
      WRITE(2,'(/A)')'PREPARATION FOR TIME STEPPING SECTION SUCCESSFULLY
     & COMPLETED'

CGH      STOP
      WRITE(*,*)' SACO_RUN about to start time stepping '
      write(*,*)zetaold(1),zetamid(1),zmid(1,1)

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C TIME STEPPING LOOP
C
  123 CONTINUE            
C
C-----------------------------------------------------------------------
C SET UP FOR THIS TIME STEP
C
C KDMID,SECMID are the timing parameters for the beginning of the
C current time step => time level K
C Set the timing parameters for the time at the end of the current
C time step         => time level K+1
C
      KDNEW=KDMID
      SECNEW=SECMID+DELT
      ITER=ITER+1
      IF(SECNEW.GE.86400.)CALL UP_DATE(KDNEW,SECNEW)
C
C
C Zero appropriate arrays,
C
      DO J=1,NNV
         DO I=1,NN
            UZNEW(I,J)=0.0
            VZNEW(I,J)=0.0
         ENDDO
      ENDDO
C
C Evaluate the following explicitly for the current time step (i.e. at
C time level K):
C
C
C-----------------------------------------------------------------------
C SOLVE THE WAVE EQUATION
C
C Objectives:Determine elevation at time level     K+1   (=>ZETANEW)
C            using known info at time levels       K-1   (=>ZETAOLD,*)
C            and                                   K     (=>ZETAMID,*)
C            with user supplied BC's at time level K+1   (=>*)
C            vertical grid is constant at ZMID
C
      CALL ELEVATIONS4(NN,NHBW,NNV,DELT,TAU0,
     &   KDNEW,SECNEW,ITER,DEG,NR,NVN,THETA,HDOWN,
     &   ZETAOLD,UOLD,VOLD,ZETAMID,UMID,VMID,ZMID,UZMID,VZMID,
     &   LRVS,NORM,SUMX,SUMY,DSBDRY,SH,Emat,
     &   Amat,Bmat,Cxmat,Cymat,Dxmat,Dymat,IQ,JQ,
     &   ZETANEW)

c     CALL ELEVATIONS2(NN,NE,NHBW,NNV,DELT,G,TAU0,COR,
c    & KDNEW,SECNEW,ITER,DEG,NR,NVN,THETA,IN,DX,DY,AR,HDOWN,
c    & ZETAOLD,UOLD,VOLD,ZETAMID,UMID,VMID,ZMID,UZMID,VZMID,
c    & AK,LRVS,NORM,SUMX,SUMY,DSBDRY,SV,SH,EMAT,ZETANEW)

C
C Terminate execution if any negative depths exist
C            
      DO I=1,NN
         IF(ZETANEW(I)+HDOWN(I).LE.0.0)THEN
            WRITE(2,'(a)') 'NEGATIVE DEPTH ENCOUNTERED -- ',
     &                     'EXECUTION TERMINATED'
            WRITE(2,'(A,I5,A,F12.4,A,F12.4)')
     &       'NEG DEPTH AT NODE:',I,' X=',X(I),' Y=',Y(I)
	    WRITE(2,*)'ITER, KDNEW, SECNEW ', ITER, KDNEW, SECNEW  
	    WRITE(2,*)' ZETA, HDOWN = ', ZETANEW(I),HDOWN(I)

            CALL DMY(IDAYMID,IMNTHMID,IYEARMID,KDMID)
            WRITE(2,*)' date', IDAYMID,IMNTHMID,IYEARMID,'time',SECMID
            WRITE(*,*)' date', IDAYMID,IMNTHMID,IYEARMID,'time',SECMID
            CLOSE(2)
            STOP
         ENDIF
      ENDDO
C
C-----------------------------------------------------------------------
C SOLVE FOR THE VERTICAL STRUCTURE OF THE 3-D VARIABLES
C
C Objectives: Assemble and solve tridiagonal momentum,
C             equations at the center of the time step, then convert to
C             values at the end of the time step.
C             => UZNEW, VZNEW    at time level K+1
C
C Compute vertical structure
C
      CALL VERTICALS5(NN,NNV,DELT,SV,G0xmat,G0ymat,IQ,JQ,
     & D1,D2,D3,DZ, ZETAMID,UZMID,VZMID,ZETANEW,UZNEW,VZNEW)
C
C
C-----------------------------------------------------------------------
C INCREMENT TIMING PARAMETERS AND UPDATE/JUGGLE ARRAYS 
C
C Objectives: Increment time such that, for the next time step:
C               *OLD arrays contain information at time level K-1
C               *MID arrays contain information at time level K
C
C *OLD<=*MID:
C
      DO I=1,NN
         ZETAOLD(I)=ZETAMID(I)              
         UOLD(I)=UMID(I)              
         VOLD(I)=VMID(I)              
      ENDDO
C
C *MID<=*NEW:  Thus *MID arrays will contain information for the 
C              current time => (KDMID,SECMID), which is time level K 
C              for the next time step.
C
      KDMID=KDNEW
      SECMID=SECNEW
      DO J=1,NNV
         DO I=1,NN
            UZMID(I,J)=UZNEW(I,J)
            VZMID(I,J)=VZNEW(I,J)
         ENDDO
      ENDDO
      DO I=1,NN
         ZETAMID(I)=ZETANEW(I)       
      ENDDO
c     ii = 563
c     write(*,*)'after update  ',uzmid(ii,1),uznew(ii,1)
c    &  ,vzmid(ii,1),vzmid(ii,1),enzm(ii,1)
C
C Result: Arrays at the end of the current time step 
C                                      are now stored as *MID 
C
C-----------------------------------------------------------------------
C COMPUTE VERTICALLY AVERAGED VELOCITIES AT CURRENT TIME (*MID)
C
      CALL VERTAVG(UZMID,UMID,ZMID,NN,NNV,NNDIM)
      CALL VERTAVG(VZMID,VMID,ZMID,NN,NNV,NNDIM)
C
C-----------------------------------------------------------------------
C
C TIME STEP IS COMPLETE  
C
C Compute the difference between the model and the observations.
C Write the error field to arrays defined in observes.h
C
      CALL SACO_UVERROR(ITER,IN,NNV,UZMID,VZMID)
      CALL SACO_EERROR(ITER,IN,ZETAMID)
C
C CALL OUTPUT WITH RESULTS AT THE CURRENT TIME (*MID)
C
c     write(*,*)'saco_run, kdmid, secmid ', KDMID,SECMID
      CALL OUTPUTS2(ITER,KDMID,SECMID,DELT,MESHNAME,FILINQ,FILICQ,
     & NN,NNV,NE,X,Y,IN,HDOWN, AK, ENZM, ZETAMID,UMID,VMID,
     & ZMID,UZMID,VZMID,ZETAOLD,UOLD,VOLD,casco_iter)

C
C-----------------------------------------------------------------------
C  ELAPSED TIME CHECK, RETURN FOR ANOTHER TIME STEP OR STOP 
C  (AS APPROPRIATE)
C
      IF(KDMID.LT.KDL)GO TO 123
      IF((KDMID.EQ.KDL).AND.(SECMID.LT.SECL))GO TO 123
C
      WRITE(2,'(/A)')'TIME STEPPING SECTION SUCCESSFULLY COMPLETED'
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C END OF SIMULATION, COMPLETE .echo FILE AND CALL OUTPUTQ4 ONE LAST TIME
C (WITH ITER SET TO -1)
C
      ITER=-1
      write(*,*)'calling outputs2. final time'
      CALL OUTPUTS2(ITER,KDMID,SECMID,DELT,MESHNAME,FILINQ,FILICQ,
     & NN,NNV,NE,X,Y,IN,HDOWN, AK, ENZM, ZETAMID,UMID,VMID,
     & ZMID,UZMID,VZMID,ZETAOLD,UOLD,VOLD,casco_iter)
      
C
      CALL DMY(IDAYMID,IMNTHMID,IYEARMID,KDMID)
      WRITE(2,'(/A)')'SIMULATION SUCCESSFULLY COMPLETED'
      WRITE(2,'(/5X,A35,1X,I2,1X,I2,3X,I4,1X,F18.8)')
     & 'Simulation end date and time       ',
     & IDAYMID,IMNTHMID,IYEARMID,SECMID
c     CLOSE(2)
C
      RETURN
      END        
C***********************************************************************
C***********************************************************************
