C***********************************************************************
C***********************************************************************

      SUBROUTINE OUTPUTS2(ITER,KD,SECKD,DELT,MESH,FILINQ,FILICQ,
     & NN,NNV,NE,X,Y,IN,HDOWN, AK, ENZM,  
     & ZETA,UBAR,VBAR,Z,U,V,ZETAOLD,UBAROLD,VBAROLD, casco_iter)

C Include file for parameter statements
      INCLUDE 'CASCO.DIM'
C Fixed global variables
      INTEGER ITER
      INTEGER KD
      REAL    SECKD
      REAL    DELT
      INTEGER NN,NNV,NE
      REAL    X(NNDIM),Y(NNDIM),HDOWN(NNDIM)
      REAL    AK(NNDIM)
      INTEGER IN(3,NEDIM)
      REAL ZETA(NNDIM),UBAR(NNDIM),VBAR(NNDIM)
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM)
      REAL ENZM(NNDIM,NNVDIM)
      REAL ZETAOLD(NNDIM),UBAROLD(NNDIM),VBAROLD(NNDIM)
      CHARACTER*72 MESH,FILINQ,FILICQ
      integer casco_iter  
C-----------------------------------------------------------------------
C This user subroutine is supplied to act as an interface for the
C   user to output the 2-D and 3-D results at a particular time in the
C   simulation.
C
C Index range:  I=1,NN, J=1,NNV, K=1,NE, L=1,NEV
C
C   NN - number of horizontal nodes
C   NE - number of horizontal elements
C   NNV - number of vertical nodes
C   NEV - number of vertical elements (NEV=NNV-1)
C
C Timing parameters:
C
C   KD,SECKD - Gregorian calendar day and seconds since the beginning
C     of that calendar day.
C   ITER - Iteration number for the current time step
C   DELT - Time step size in seconds
C 
C Constant information regarding the horizontal mesh
C
C   X(I),Y(I) - horizontal cartesian nodal coordinates
C   HDOWN(I) - horizintal nodal bathymetry
C   IN(I,K) - horizontal triangular element incidence list
C
C The following arrays contain the 2-D and 3-D solution at the end of
C   each iteration, at the time indicated by KD and SECKD:
C
C   ZETA(I) -  free surface elevation
C   UBAR(I),VBAR(I) - vertically averaged velocity
C   Z(I,J) - nodal coordinate locations in 3-D
C   U(I,J),V(I,J) - nodal values of the X and Y components
C     of velocity.
C   ENZM(I,L) elemental values of the vertical diffusivities 
c	      for momentum
C
C The following arrays contain the 2-D and 3-D solution at one time
C   step prior to the end of the simulation.  They are particularly 
C   useful in the construction of ".icq4" HOT-START files:
C
C   ZETAOLD(I) - free surface elevation
C   UBAROLD(I),VBAROLD(I) - vertically averaged velocity
C
C Z and J are positive upward:
C          Bottom:  Z(I,1)=-HDOWN(I); J=1
C          Surface: Z(I,NNV)= 0  
C
C The following are the parameters related to the simulation at hand;
C   which are particularly useful in the construction of ".icq4"
C   HOT-START files:
C
C   MESH - Name of the mesh
C   FILINQ - Name if the .inq file being used to control the simulation
C   FILICQ - Name of the .icq4 file used to initialize the simulation
C
C This subroutine is called at the beginning of the simulation 
C   (with ITER=0) for initialization purposes, after each time step
C   (with ITER>0) to facilitate any recording of results desired by the
C   user and finally at the successful completion of the simulation
C   (with ITER set to -1) to facilitate any remaining post-processing
C   or recording.
C-----------------------------------------------------------------------
C A standard routine 'SUBROUTINE ICQ4WRITE' writes all computed results
C   to a user-specified file.
C
C      CALL ICQ4WRITE(MESH,FILINQ,FILICQ,FILOUT,
C     &KD,SECKD,DELT,ZETA,UBAR,VBAR,ZETAOLD,UBAROLD,VBAROLD,
C     &Z,ZOLD,U,V,W,Q2,Q2L,T,S,NN,NNV,NNDIM,NNVDIM)
C
C================Beginning-of-user-specified-instructions===============
c
c Array assignments and save parameters for SUBROUTINE DECOMPOSE
c     INTEGER M, K 
c     PARAMETER(M=6)
c      PARAMETER(M=2)
c     REAL ELEV(NNDIM,700),AMPRES(NNDIM,M),PHASERES(NNDIM,M)
c     REAL F(M),VU(M),SIG(M)
c     CHARACTER*5 KON(M)
c     DATA (KON(K),K=1,M)/'Z0','M2','N2','S2','K1','O1'/
c      DATA (KON(K),K=1,M)/'Z0','M2'/

c     INTEGER KDOLD 
c     INTEGER NHOUR 
c     INTEGER IHH1,IDD1,IMM1,IYEAR1,KD1
c     INTEGER IHHF,IDDF,IMMF,IYEARF,KDF
c     REAL    SECKDOLD 
c     SAVE ELEV,KDOLD,SECKDOLD,NHOUR
c     SAVE IHH1,IDD1,IMM1,IYEAR1,KD1
c     SAVE IHHF,IDDF,IMMF,IYEARF,KDF

c
c Character assignment for SUBROUTINE ICQ4WRITE
c     CHARACTER*72 FILOUT
c
c Save parameters for SUBROUTINE DFTCENQ4
c     INTEGER NPRD
c     REAL    FREQ
c     SAVE NPRD,FREQ
c
c Assignments and saves for daily surface and bottom values
c     REAL TAVG(NNDIM),SAVG(NNDIM)
c     CHARACTER*72 INFILE,CKD,HEADER


c     REAL SECKD0,SECKDOLD1
c     INTEGER KD0 
c     INTEGER OUTSET
c     SAVE KD0,SECKD0,SECKDOLD1, OUTSET

c     INTEGER I 
c     INTEGER NFLAG 
c     REAL    B1, B2 
c     INTEGER NPERIOD
c     REAL DEGLAT 
c     INTEGER IHHKD, IHHKDOLD 
c     REAL    TIME

      character*120 fileroot
      integer istart1,iend1
      integer istart2,iend2
      character*12 citer
      character*120 fname

      integer i 

c
c Echo the iteration number to standard out
c     WRITE(*,*)'Output: iter, kd, seckd',ITER,kd,seckd
c
c     return 
c     WRITE(*,*)ITER
c
c Perform desired comparisons with observation locations
c     CALL INTOBSDT(KD,SECKD,ITER,DELT,IN,X,Y,HDOWN,NN,NNV,NE,
c    &UBAR,VBAR,Z,U,V,ENZM,MESH)

c     CALL INTOBSNEP(KD,SECKD,ITER,DELT,IN,X,Y,HDOWN,NN,NNV,NE,
c    &UBAR,VBAR,Z,U,V,ENZM,MESH)
c     CALL INTOBSST2(KD,SECKD,ITER,DELT,IN,X,Y,HDOWN,NN,NNV,NE,
c    &UBAR,VBAR,Z,U,V,ENZM,MESH)
c
c
c Write .icq4 file at the end of the simulation
*      IF(ITER.EQ.2473)THEN
*         FILOUT='1995_04_05.icq4'
*         CALL ICQ4WRITE(MESH,FILINQ,FILICQ,FILOUT,
*     &   KD,SECKD,DELT,ZETA,UBAR,VBAR,ZETAOLD,UBAROLD,VBAROLD,
*     &   Z,ZOLD,U,V,W,Q2,Q2L,T,S,NN,NNV,NNDIM,NNVDIM)
*      ENDIF
c

       call  ts1_saco_write(casco_iter,kd,seckd,iter,delt,nn,
     $                      zeta,ubar,vbar)


       if(casco_iter.lt.0) then 
       if(iter.lt.0) then 
	  write(*,*)'Outputs2: try and write final zeta and uvbar'
           fileroot = 'zeta' 
           if(casco_iter.gt.0) then 
             write(citer,'(i2.2)')casco_iter
             call cstring(citer,istart1,iend1)
             call cstring(fileroot,istart2,iend2)
             fname =  fileroot(istart2:iend2) // '.'
     $           // citer(istart1:iend1) // '.s2r'
	   elseif(casco_iter.lt.0) then 
             call cstring(fileroot,istart2,iend2)
             fname =  fileroot(istart2:iend2) // '.s2r'
	   endif
           write(*,*)'Output file name ',fname

           open(48,file=fname)
           write(48,'(a)')mesh
           write(48,*)'elevation at end of simulation'
           do i = 1, nn
              write(48,*)i,zeta(i)
           enddo
           close(48)
 
           fileroot = 'uvbar' 
           if(casco_iter.gt.0) then 
             write(citer,'(i2.2)')casco_iter
             call cstring(citer,istart1,iend1)
             call cstring(fileroot,istart2,iend2)
             fname =  fileroot(istart2:iend2) // '.'
     $           // citer(istart1:iend1) // '.v2r'
	   elseif(casco_iter.lt.0) then 
             call cstring(fileroot,istart2,iend2)
             fname =  fileroot(istart2:iend2) // '.v2r'
	   endif

           write(*,*)'Output file name ',fname
           open(48,file=fname)
           write(48,'(a)')mesh
           write(48,*)'ubar, vbar at end of simulation'
           do i = 1, nn
              write(48,*)i,ubar(i), vbar(i)
           enddo
           close(48)

       endif
       endif


      RETURN
c
c Save hourly elevations for decomposition
c
c  Initialize parameters for decomp when ITER=0
c     IF(ITER.EQ.0)THEN
c
c   set hour/day/month/year/century for beginning analysis
c     *1=>beginning of storage for analysis 
c         Notes: SUBROUTINE BCQ4 is ramping up the BC's for the
c                  first 24 hours => they are at full strength 
c                  1 day into the simulation).
c                I will wait 4 more days before saving data for
c                  decomposition to allow for convergence of the model.
c                Recall that KD and SECKD are starting time for the
c                  simulation when ITER=0.
c                If SECKD does not correspond to an integer hour, then
c                  storage for decomposition will begin at a time which
c                  is not exactly 4 days after the start of the
c                  simulation.
c
c        KD1=KD+1+4
c        IHH1=INT(SECKD)/3600
c        CALL DMY(IDD1,IMM1,IYEAR1,KD1)
c
c   initialize local timing parameters needed for temporal interpolation
c        KDOLD=KD
c        SECKDOLD=SECKD
c        NHOUR=0
c        WRITE(*,*)'Initializing for decomposition interval:'
c        WRITE(*,*)'  ITER,NHOUR =',ITER,NHOUR
c        WRITE(*,*)'  IHH1,IDD1,IMM1,IYEAR1=',IHH1,IDD1,IMM1,IYEAR1
c        WRITE(*,*)'  KD,HOURKD  =',KD,INT(SECKD)/3600
c        WRITE(*,*)'  KD1,IHH1   =',KD1,IHH1
c        WRITE(99,*)'Initializing for decomposition interval:'
c        WRITE(99,*)'  ITER,NHOUR =',ITER,NHOUR
c        WRITE(99,*)'  IHH1,IDD1,IMM1,IYEAR1=',IHH1,IDD1,IMM1,IYEAR1
c        WRITE(99,*)'  KD,HOURKD  =',KD,INT(SECKD)/3600
c        WRITE(99,*)'  KD1,IHH1   =',KD1,IHH1
c     ELSE
c
c  Save the elevations every hour starting at hour IHH1 on KD1
c    (NFLAG will be set to 1 and elevations will be stored in
c     decomposition array for timesteps during which the minute hand 
c     passes the 12 at or after IHH1 on KD1)
c        NFLAG=0
c        IHHKD=INT(SECKD)/3600
c        IHHKDOLD=INT(SECKDOLD)/3600
c        IF((KD.GE.KD1).AND.(IHHKD.GE.IHH1))THEN
c     Check for the special case that midnight has been passed
c           IF(KD.EQ.(KDOLD+1))THEN
c              IF((IHHKDOLD.NE.23).OR.(IHHKD.NE.0))THEN
c                 WRITE(*,*)'Potential error in decomposition'
c                 STOP
c              ENDIF
c              NFLAG=1
c              B1=SECKD/DELT
c              B2=(3600.0*24.0-SECKDOLD)/DELT
c     Make sure KD equals KDOLD if midnight has not been passed
c           ELSE IF(KD.NE.KDOLD)THEN
c              WRITE(*,*)'Potential error in decomposition'
c              STOP
c     Check to see if an hour has ellapsed
c           ELSE IF(IHHKD.EQ.IHHKDOLD+1)THEN
c              NFLAG=1
c              B1=AMOD(SECKD,3600.0)/DELT
c              B2=(3600.0-AMOD(SECKDOLD,3600.0))/DELT
C     Make sure IHHKD equals IHHKDOLD if an hour has not ellapsed
c           ELSE IF(IHHKD.NE.IHHKDOLD)THEN
c              WRITE(*,*)'Potential error in decomposition'
c              STOP
c           ELSE
c           ENDIF
c     NFLAG=1
c           IF(NFLAG.EQ.1)THEN
c
c       Call ZTM3D (note that the first time NFLAG=1, NHOUR=0 here)
C              CALL ZTM3D(KD,SECKD,NHOUR,
C    &         NN,NNV,NE,X,Y,Z,IN,RHO,U,V,W,ENZM,Q2,Q2L)
c       Save elevations
c              NHOUR=NHOUR+1
c              DO I=1,NN
c                ELEV(I,NHOUR)=B1*ZETAOLD(I)+B2*ZETA(I)
c                ELEV(I,NHOUR)=B1*ZOLD(I,NNV)+B2*Z(I,NNV)
c              ENDDO
c       Save actual starting time for analysis when NHOUR equals 1
c              IF(NHOUR.EQ.1)THEN
c                 CALL DMY(IDDF,IMMF,IYEARF,KD)
c                 KDF=KD
c                 IHHF=IHHKD
c              ENDIF
c              WRITE(*,'(3i8,6f10.3)')ITER,NHOUR,KD,SECKD,B1,B2,
c    &                   ZETAOLD(1),ZETA(1),ELEV(1,NHOUR)
c              WRITE(99,'(3i8,6f10.3)')ITER,NHOUR,KD,SECKD,B1,B2,
c    &                   ZETAOLD(1),ZETA(1),ELEV(1,NHOUR)
c           ENDIF
c        ENDIF
c     ENDIF
c
c  Now update KDOLD and SECKDOLD for next pass through OUTPUT
c     KDOLD=KD
c     SECKDOLD=SECKD
c
c Perform tidal decomposition at the end of the simulation
c     IF(ITER.EQ.-1)THEN
c
c  Set latitude, echo constants, and perform decomposition:
c        DEGLAT=43.50
c        WRITE(*,*)'Calling DECOMPCEN'
c        WRITE(*,*)'  NNDIM,NN=',NNDIM,NN
c        WRITE(*,*)'  IHHF,IDDF,IMMF,IYEARF=',IHHF,IDDF,IMMF,IYEARF
c        WRITE(*,*)'  NHOUR=',NHOUR
c        WRITE(*,*)'  DEGLAT=',DEGLAT
c        WRITE(*,*)'  (KON(I),I=1,M)=',(KON(I),I=1,M)
c        WRITE(99,*)'Calling DECOMPCEN'
c        WRITE(99,*)'  NNDIM,NN=',NNDIM,NN
c        WRITE(99,*)'  IHHF,IDDF,IMMF,IYEARF=',IHHF,IDDF,IMMF,IYEARF
c        WRITE(99,*)'  NHOUR=',NHOUR
c        WRITE(99,*)'  DEGLAT=',DEGLAT
c        WRITE(99,*)'  (KON(I),I=1,M)=',(KON(I),I=1,M)
c        CALL DECOMPCEN(NNDIM,NN,ELEV,IHHF,IDDF,IMMF,IYEARF,NHOUR,
c    &                  DEGLAT,M,KON,F,VU,SIG,AMPRES,PHASERES)
c        CALL IOSTOS2RC(MESH,NNDIM,NN,IHHF,IDDF,IMMF,IYEARF,NHOUR,
c    &                         M,KON,F,VU,SIG,AMPRES,PHASERES)
c        WRITE(*,*)'Returned from DECOMPCEN'
c     ENDIF
C
C End of routine
C
      RETURN
      END
