CO***********************************************************************
C***********************************************************************
C Beginning of User-Specified Subroutines for QUODDY5 (RELEASE 1.1)
C Beginning of User-Specified Subroutines for QUODDY5 (RELEASE 1.1)
C***********************************************************************
C***********************************************************************
      SUBROUTINE BCQ5(KD,SECKD,ITER,DELT,DEGLAT,XNOD,YNOD,
     &VAR1,VAR2,NBCS,III)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL XNOD(NNDIM),YNOD(NNDIM),VAR1(NNDIM),VAR2(NNDIM)
C DJMFEW
      REAL UBARTD(NNDIM)
C End DJMFEW
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper
C   Boundary Conditions.  There are four sections which require
C   modification, beginning at statement labels 3,4,5, and 6.  
C Each section must begin with the appropriately labelled continue 
C   statement, and end with a return statement.
C This subroutine is called once at the beginning of the simulation 
C   (with ITER=0 and III=0) for initialization purposes; then it is
C   called as necessary during each time step to set the desired
C   horizontal boundary conditions.
C================Beginning-of-user-specified-instructions===============
C
C User-specified local assignments
C
C Beginning of executable code
C
C  Declarations for multiple tidal constituents and perturbations:
      PARAMETER(M=2,NPRTDIM=18)
      REAL F(M),VU(M),SIG(M),RDUM
      REAL AMP(M,NBEDIM),PHASE(M,NBEDIM)
      INTEGER KDPRT(NPRTDIM)
      REAL SECKDPRT(NPRTDIM),FPRT(NPRTDIM)
      REAL AMPPRT(NPRTDIM,NBEDIM),PHASEPRT(NPRTDIM,NBEDIM)
      REAL RAMPKDPRT(NPRTDIM),RAMPPRT(NPRTDIM)
C DJMFEW
      REAL RAMPTD,RAMPUPTD
      INTEGER IYEARTD,YDGMT0TD
      INTEGER KDSTARTTD,SECKDSTARTTD
C End DJMFEW
      INTEGER NNIBCS(2),IBCS(2,NBEDIM)
      CHARACTER*72 HEADER,BCSPATH,BCSFILE,INTERNAL,CKP
      CHARACTER*5 KON(M),CDUM
CVL DECLARATIONS FOR TDBCS FLAG 3/30/99
      integer ITDBCSFG
C BOF normal flow
      REAL AMPUnz(NBEDIM,NFDIM),PHASEUnz(NBEDIM,NFDIM)
      REAL AMPVnz(NBEDIM,NFDIM),PHASEVnz(NBEDIM,NFDIM)
      INTEGER NODEnz(NBEDIM),NNODESnz
C end BOF normal flow
      save ITDBCSFG
CVL END NEW DECLARATIONS
      SAVE NNIBCS,IBCS,AMP,PHASE,KD0,SECKD0,KDVUF,F,VU,SIG,
     &NPRT,AMPPRT,PHASEPRT,KDPRT,SECKDPRT,FPRT
C non-zero velocity condition
      SAVE NNODESnz,NODEnz,AMPUnz,PHASEUnz,AMPVnz,PHASEVnz
C end non-zero
      SAVE RAMPKDGMT,RAMPKDPRT
C DJMFEW
      SAVE IYEARTD,YDGMT0TD,RAMPUPTD
C End DJMFEW
CVL DATA STATEMENT MOVED FOR UNC COMPILATION 3/30/99
C      DATA (KON(K),K=1,M)/'Z0','M2','N2','S2','K1','O1'/
       DATA (KON(K),K=1,M)/'Z0','M2'/
C
C Beginning of executable code
CVL BEGIN NEW CODE
CVL CALL TIME DEPENDENT BCS

      if (ITER.eq.0) then
         write (*,*) 'SUBROUTINE BCSTD INITIALIZATION:'
         write(*,*) 'Use time dependent BCs (1/0)?'
         read(*,*) ITDBCSFG
C DJMFEW
         if (ITDBCSFG.eq.1) then
            WRITE(*,*)'Enter time to begin TDBCS
     &                 (year,ydgmt0) and the rampup period (in days):'
            READ(*,*)IYEARTD,YDGMT0TD,RAMPUPTD
          endif

      endif

      if (ITDBCSFG.eq.1) then

         CALL GMT0_GDAY(IYEARTD,YDGMT0TD,KDSTARTTD,SECKDSTARTTD)

         DELTAKDTD=REAL(KD-KDSTARTTD)+(SECKD-SECKDSTARTTD)/24.0/3600.0

         IF(RAMPUPTD.EQ.0)THEN
            RAMPTD=1.0
           ELSE
            RAMPTD=DELTAKD/RAMPUPTD
            RAMPTD=AMAX1(RAMPUPTD,0.0)
            RAMPTD=AMIN1(RAMPUPTD,1.0)
          ENDIF

         call BCSTD(KD,SECKD,ITER,UBARTD,VBAR,III)

       else

         DO I=1,NNIBCS(1)
            NBC=IBCS(1,I)
            UBARTD(NBC)=0.0
         enddo
      endif

C End DJMFEW

      IF(III.EQ.0)GO TO 3
      IF(III-2)4,5,6
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C This section may be used for initialization purposes.
C
    3 CONTINUE

      WRITE(*,'(/A)')'SUBROUTINE BCQ5 INITIALIZATION:'
c
c   Set initial timing parameters
      KD0=KD
      SECKD0=SECKD
c
c   Read elevation amplitudes and phase BCS for M tidal consitiuents:
c     Note:  The frequencies are set using IOS routines after all of
c            the nodal values for Greenwich amplitude and phase are
c            entered at each of the M frequencies under consideration.
c

c    Set parameters
      IF(M.EQ.0)GO TO 12
      BCSPATH='CLIM_BCS/BCS_'
      CALL CSTRING(BCSPATH,IS,IE)
      WRITE(*,*)'Enter the ramp up period (in days) for the climatologic
     &al boundary conditions:'
      READ(*,*)RAMPKDGMT
      WRITE(*,*)'GMT BC Files: Ramp up period (days) =',RAMPKDGMT
      DO K=1,M
         KM=2
         IF(K.EQ.1)KM=1
c
c    Read individual data files
         INDX=INDEX(KON(K),' ')-1
         BCSFILE=BCSPATH(IS:IE)//KON(K)(1:INDX)//'.s2c'
         OPEN(21,file=BCSFILE,STATUS='OLD')
         READ(21,'(A)')HEADER
         READ(21,'(A)')HEADER
         READ(21,'(A)')HEADER
         DO I=1,NBEDIM+1
            READ(21,*,END=11)IBCS(KM,I),AMP(K,I),PHASE(K,I)
            NNIBCS(KM)=I
         ENDDO
         WRITE(*,*)'ERROR READING BC FILE ... NUMBER EXCEEDS NBEDIM'
         STOP
 11      CLOSE(21)
c
c    Check for consistency with previously read files
         IF(KM.EQ.2)THEN
            IF(NNIBCS(1).NE.NNIBCS(2))THEN
               WRITE(*,*)'INCONSISTENT NUMBER OF BCS IN FILES'
               STOP
            ENDIF
            DO I=1,NNIBCS(1)
               IF(IBCS(1,I).NE.IBCS(2,I))THEN
                  WRITE(*,*)'INDEXING MISMATCH FOR BC FILES:'
                  STOP
               ENDIF
            ENDDO
         ENDIF
         WRITE(*,'(2X,I5,1X,A,1X,A,A)')NNIBCS(KM),KON(K)(1:INDX),
     &      'Elev. BCs successfully read from ',BCSFILE
      ENDDO

c
c   Use IOS files to load F(K) and VU(K) for the beginning of the
c    the current Gregorian Day (Day=KD,hour=KD*24):
      KDVUF=KD
      CALL DMY(ID,IM,IYEAR,KDVUF)
      KHVUF=KDVUF*24
      CALL OPNVUF(KHVUF,CDUM,DEGLAT,RDUM,RDUM,RDUM)
      CLOSE(8)
      DO K=1,M
         CALL VUF(KHVUF,KON(K),DEGLAT,F(K),VU(K),SIG(K))
      ENDDO
      WRITE(*,*)'  Tidal parameters initialized from IOS routines at:'
      WRITE(*,*)'   (day,month,year,sec)=',ID,IM,IYEAR,0.0
      WRITE(*,*)'   DEGLAT=',deglat
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      WRITE(*,'(3X,A)')
     &'KON     F            VU             SIG           FREQ'
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      DO K=1,M
         FREQK=2.0*ACOS(-1.0)*SIG(K)/3600.0
         WRITE(*,'(3X,A,1X,4e14.6)')KON(K),F(K),VU(K),SIG(K),FREQK
      ENDDO
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
c
c   Read elevation amplitudes and phase BCS for NPRT perturbations:
c     Notes: The frequencies are read directly from the files as are
c            the reference times for the phases of the BCS.
c
c    Set parameters

 12   WRITE(*,*)'Enter the number of perturbation frequencies'
      READ(*,*)NPRT
      IF(NPRT.EQ.0)GO TO 22
      IF(NPRT.GT.NPRTDIM)THEN
         WRITE(*,*)'EXECUTION TERMINATED:'
         WRITE(*,*)'TO MANY PERTURBATION BC FILES IN SUBROUTINE BCQ5'
         STOP
      ENDIF
      BCSPATH='CRUISE_T6.ALL/BCS_PRT'
      CALL CSTRING(BCSPATH,IS,IE)
      WRITE(*,*)'Enter the ramp up period (in days) for each set of pert
     &urbation boundary conditions:'
c      READ(*,*)(RAMPKDPRT(KP),KP=1,NPRT)
      READ(*,*)RAMPKDPRT(1)
      do k=2,NPRT
       RAMPKDPRT(k)=RAMPKDPRT(1)
      enddo
      DO KP=1,NPRT
         KM=2
         IF((M.EQ.0).AND.(KP.EQ.1))KM=1
c
c    Read individual data files
         WRITE(INTERNAL,'(I10)')KP
         READ(INTERNAL,'(A)')CKP
         CALL CSTRING(CKP,IS1,IE1)
         BCSFILE=BCSPATH(IS:IE)//CKP(IS1:IE1)//'.s2c'
         OPEN(21,file=BCSFILE,STATUS='OLD')
         READ(21,'(A)')HEADER
         PRINT*,HEADER
         READ(21,*)IYEAR,YDAYGMT0
         PRINT*,IYEAR, YDAYGMT0
         CALL GMT0_GDAY(IYEAR,YDAYGMT0,KDPRT(KP),SECKDPRT(KP))
         READ(21,*)FPRT(KP)
         PRINT*,FPRT(KP)
         DO I=1,NBEDIM+1
            READ(21,*,END=21)IBCS(KM,I),AMPPRT(KP,I),PHASEPRT(KP,I)
            NNIBCS(KM)=I
            PRINT*,IBCS(KM,I),AMPPRT(KP,I),PHASEPRT(KP,I)
         ENDDO
         WRITE(*,*)'ERROR READING BC FILE ... NUMBER EXCEEDS NBEDIM'
         STOP
 21      CLOSE(21)
c
c    Check for consistency with previously read files
         IF(KM.EQ.2)THEN
            IF(NNIBCS(1).NE.NNIBCS(2))THEN
               WRITE(*,*)'INCONSISTENT NUMBER OF BCS IN FILES'
               STOP
            ENDIF
            DO I=1,NNIBCS(1)
               IF(IBCS(1,I).NE.IBCS(2,I))THEN
                  WRITE(*,*)'INDEXING MISMATCH FOR BC FILES:'
                  STOP
               ENDIF
            ENDDO
         ENDIF
         WRITE(*,'(2X,I5,1X,A,1X,A,A)')NNIBCS(KM),KON(KP)(1:INDX),
     &      'Elev. BCs successfully read from ',BCSFILE
         WRITE(*,'(6X,A,1X,I4,1X,F10.4)')
     &      'reference time (IYEAR,YDAYGMT0) =',IYEAR,YDAYGMT0
         WRITE(*,'(6X,A,1X,E14.6)')
     &      'frequency (rad/sec) =',FPRT(KP)
         WRITE(*,'(6X,A,1X,E14.6)')
     &      'ramp up period (days) =',RAMPKDPRT(KP)
      ENDDO

c hardwire in option to use non-zero velocity conditions
  22   NFnz = 1    ! 1 RHE 01/07/2004
      BCSPATH='CLIM_BCS/BCS_'
      CALL CSTRING(BCSPATH,IS,IE)
       if (NFnz.eq.0)goto 23
       DO k=1,M
c load in the v2c files for normal velocity BC
         INDX=INDEX(KON(K),' ')-1
         BCSFILE=BCSPATH(IS:IE)//KON(K)(1:INDX)//'.v2c'
         OPEN(31,file=BCSFILE,STATUS='OLD')
         READ(31,'(A)')HEADER
         READ(31,'(A)')HEADER
         READ(31,'(A)')HEADER
         DO I=1,NBEDIM+1
           READ(31,*,END=77)NODEnz(I),AMPUnz(I,k),PHASEUnz(I,k),
     &                              AMPVnz(I,k),PHASEVnz(I,k)
             NNODESnz=I
         ENDDO
         WRITE(*,*)'ERROR READING BC FILE ... NUMBER EXCEEDS NBEDIM'
         STOP
 77      CLOSE(31)
       ENDDO

 23    RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Surface Elevation BC's are assigned in this section.  One statement of
C   the form VAR1(I)=ELEV must appear for each specified surface
C   elevation node, where I is the horizontal node number and ELEV is
C   the elevation for node I.  
C A return statement must follow the last elevation BC assignment.
C
    4 CONTINUE
c
c  Set ramp up factors for perturbation and GMT bcs with respect to the
c   beginning of the simulation.
      DELTAKD=REAL(KD-KD0)+(SECKD-SECKD0)/86400.0
      IF(RAMPKDGMT.EQ.0)THEN
         RAMPGMT=1.0
      ELSE
         RAMPGMT=DELTAKD/RAMPKDGMT
         RAMPGMT=AMIN1(RAMPGMT,1.0)
      ENDIF
      DO KP=1,NPRT
         IF(RAMPKDPRT(KP).EQ.0)THEN
            RAMPPRT(KP)=1.0
         ELSE
            RAMPPRT(KP)=DELTAKD/RAMPKDPRT(KP)
            RAMPPRT(KP)=AMIN1(RAMPPRT(KP),1.0)
         ENDIF
      ENDDO

c
c  Use IOS files to load F(K) and VU(K) for the beginning of the
c   the current Gregorian Day (Day=KD,hour=KD*24), if this is a new day.
      IF(KD.EQ.KDVUF)GO TO 31
      KDVUF=KD
      CALL DMY(ID,IM,IYEAR,KDVUF)
      KHVUF=KDVUF*24
      CALL OPNVUF(KHVUF,CDUM,DEGLAT,RDUM,RDUM,RDUM)
      CLOSE(8)
      DO K=1,M
         CALL VUF(KHVUF,KON(K),DEGLAT,F(K),VU(K),SIG(K))
      ENDDO
      WRITE(*,*)'  Tidal parameters reset from IOS routines at:'
      WRITE(*,*)'   (day,month,year,sec)=',ID,IM,IYEAR,0.0
      WRITE(*,*)'   DEGLAT=',deglat
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      WRITE(*,'(3X,A)')
     &'KON     F            VU             SIG           FREQ'
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      DO K=1,M
         FREQK=2.0*ACOS(-1.0)*SIG(K)/3600.0
         WRITE(*,'(3X,A,1X,4e14.6)')KON(K),F(K),VU(K),SIG(K),FREQK
      ENDDO
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
c
c  Determine the number of hours since F(M) and VU(M) were set
 31   DTHRSGMT=REAL(KD-KDVUF)*24.0+SECKD/3600.0
      TWOPI=2.0*ACOS(-1.0)
c
c  Synthesize the tides
      DO I=1,NNIBCS(1)
         NBC=IBCS(1,I)
         VAR1(NBC)=0.0
         DO K=1,M
            REVGMT=SIG(K)*DTHRSGMT+VU(K)-PHASE(K,I)/360.0
            RADGMT=TWOPI*REVGMT
            VAR1(NBC)=VAR1(NBC)+RAMPGMT*F(K)*AMP(K,I)*COS(RADGMT)
         ENDDO
         DO KP=1,NPRT
            DTRADKD=AMOD(FPRT(KP)*REAL(KD-KDPRT(KP))*86400.0,TWOPI)
            DTRADSEC=FPRT(KP)*(SECKD-SECKDPRT(KP))
            OMEGAT=DTRADKD+DTRADSEC
            RADPRT=OMEGAT-PHASEPRT(KP,I)*TWOPI/360.0
            VAR1(NBC)=VAR1(NBC)+RAMPPRT(KP)*AMPPRT(KP,I)*COS(RADPRT)
         ENDDO
C DJMFEW
C  Add in time dependent BCS
         VAR1(NBC)=VAR1(NBC)+RAMPTD*UBARTD(NBC)
C End DJMFEW
      ENDDO
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero normal velocity BC's are assigned in this section.  One 
C   statement of the form VAR1(I)=VNORM must appear for each set nonzero
C   normal velocity node, where I is the horizontal node number and 
C   VNORM is the normal velocity for node I.  
C   Positive normal velocity indicates flow out of the system.  
C A return statement must follow the last normal velocity BC assignment.
C
    5 CONTINUE

      PI=ACOS(-1.)
      DELTAKD=REAL(KD-KD0)+(SECKD-SECKD0)/86400.0
      IF(RAMPKDGMT.EQ.0)THEN
         RAMPGMT=1.0
      ELSE
         RAMPGMT=DELTAKD/RAMPKDGMT
         RAMPGMT=AMIN1(RAMPGMT,1.0)
      ENDIF

      IF(KD.EQ.KDVUF)GO TO 41
      KDVUF=KD
      CALL DMY(ID,IM,IYEAR,KDVUF)
      KHVUF=KDVUF*24
      CALL OPNVUF(KHVUF,CDUM,DEGLAT,RDUM,RDUM,RDUM)
      CLOSE(8)
      DO K=1,M
         CALL VUF(KHVUF,KON(K),DEGLAT,F(K),VU(K),SIG(K))
      ENDDO
      WRITE(*,*)'  Tidal parameters reset from IOS routines at:'
      WRITE(*,*)'   (day,month,year,sec)=',ID,IM,IYEAR,0.0
      WRITE(*,*)'   DEGLAT=',deglat
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      WRITE(*,'(3X,A)')
     &'KON     F            VU             SIG           FREQ'
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
      DO K=1,M
         FREQK=2.0*ACOS(-1.0)*SIG(K)/3600.0
         WRITE(*,'(3X,A,1X,4e14.6)')KON(K),F(K),VU(K),SIG(K),FREQK
      ENDDO
      WRITE(*,'(3X,A)')
     &'------- ------------ -------------  ------------  ------------'
c
c  Determine the number of hours since F(M) and VU(M) were set
 41   DTHRSGMT=REAL(KD-KDVUF)*24.0+SECKD/3600.0
      TWOPI=2.0*ACOS(-1.0)
      
      DO I=1,NNODESnz
         NBC=NODEnz(I)
         VAR1(NBC)=0.0
         VAR2(NBC)=0.0
         DO K = 1,M 
            REVGMT=SIG(K)*DTHRSGMT+VU(K)-PHASEUnz(I,K)/360.0
            RADGMT=TWOPI*REVGMT
            VAR1(NBC)=VAR1(NBC)+RAMPGMT*F(K)*AMPUnz(I,K)*COS(RADGMT)

            REVGMT=SIG(K)*DTHRSGMT+VU(K)-PHASEVnz(I,K)/360.0
            RADGMT=TWOPI*REVGMT
            VAR2(NBC)=VAR2(NBC)+RAMPGMT*F(K)*AMPVnz(I,K)*COS(RADGMT)
         ENDDO
       ENDDO
      RETURN
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C Nonzero velocity BC's are assigned in this section.  Two statements
C   of the following form must appear for each nonzero velocity node:
C      VAR1(I)=XVEL
C      VAR2(I)=YVEL
C   I is the horizontal node number, (XVEL,YVEL) are the
C   values of the (X,Y) velocities for node I.  
C A return statement must follow the last velocity BC assignment.
C
    6 CONTINUE
      WRITE(*,*)'SUBROUTINE BC5 called in error with III=3'
      STOP
      RETURN
C
C End of routine
C
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ATMOSQ5(KD,SECKD,ITER,NN,NNV,XNOD,YNOD,TMID,SMID,
     &ATMX,ATMY,ATEMP,BTEMP,PME)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL XNOD(NNDIM),YNOD(NNDIM),TMID(NNDIM,NNVDIM),SMID(NNDIM,NNVDIM)
      REAL ATMX(NNDIM),ATMY(NNDIM),ATEMP(NNDIM),BTEMP(NNDIM),PME(NNDIM)
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper values
C   to arrays related to the atmospheric forcing.  Assignments to the
C   following arrays must appear for each node where the atmospheric 
C   forcing is nonzero:
C
C   Index range:  I=1,NN, where NN is the number of horizontal nodes
C
C   ATMX(I),ATMY(I) - X,Y directed wind stress in pascals
C   ATEMP(I),BTEMP(I) - Terms which set the type 3 boundary condition
C                     for temperature; which is set in accordance
C                     with the expression:
C                       NdT/dz=h*(Te-T)
C                     thus:
C                       h=boundary layer coeff
C                       Te=equilibruim value of T
C                       ATEMP=h                  
C                       BTEMP=h*Te
C   PME(I) - Net precipitation minus evaporation
C
C All of these arrays are set to zero by the main routine before
C   each call to this routine.  Therefore, if no value is assigned for a
C   node in this routine, the default value of zero is implemented.
C This subroutine is called once at the beginning of the simulation 
C   (with ITER=0) for initialization purposes; then it is called once
C   during each time step to set the desired atmospheric forcing.
C================Beginning-of-user-specified-instructions===============
C
C Local Assignments
C
      CHARACTER*72 OBSFILETAU,OBSFILEQ
      SAVE KDSTART,SECKDSTART,RAMPKDATM
c
c  Perform initialization operations if ITER=0
      IF(ITER.NE.0)GO TO 10
      WRITE(*,'(/A)')'SUBROUTINE ATMOSQ5 INITIALIZATION:'
      WRITE(*,*)'Enter time to begin atmospheric forcing (year,ydgmt0) a
     &nd the rampup period (in days):'
      READ(*,*)IYEARS,YDGMT0S,RAMPKDATM
      CALL GMT0_GDAY(IYEARS,YDGMT0S,KDSTART,SECKDSTART)
      OBSFILETAU='OBS/TAUXY.o1d'
      OBSFILEQ='OBS/QFLUX.o1d'
c
c Set wind stress by interpolating temporally between observations at
c a single site
 10   CALL TAU_OF_TIME(KD,SECKD,ITER,ATMXT,ATMYT,
     &KDSTART,SECKDSTART,RAMPKDATM,OBSFILETAU,93)
      DO I=1,NN
         ATMX(I)=ATMXT
         ATMY(I)=ATMYT
      ENDDO
c
c Set the heat flux by interpolating temporally between observations at
c a single site
      CALL QFLUX_OF_TIME(KD,SECKD,ITER,BTEMPT,
     &KDSTART,SECKDSTART,RAMPKDATM,OBSFILEQ,92)
      DO I=1,NN
         PME(I)=0.0
         ATEMP(I)=0.0
         BTEMP(I)=BTEMPT
      ENDDO
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE POINTSOURCEQ5(KD,SECKD,ITER,NN,NNV,Z,
     &SRCRATE,USRC,VSRC,TSRC,SSRC)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL Z(NNDIM,NNVDIM),USRC(NNDIM,NNVDIM),VSRC(NNDIM,NNVDIM),
     &TSRC(NNDIM,NNVDIM),SSRC(NNDIM,NNVDIM)
      REAL SRCRATE(NNDIM,NNVDIM)
C Local Variables
      REAL STRENGTH,RIVTEMP

      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE POINTSOURCEQ5_RIVERS(KD,SECKD,ITER,NN,NNV,Z,
     &SRCRATE,USRC,VSRC,TSRC,SSRC)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL Z(NNDIM,NNVDIM),USRC(NNDIM,NNVDIM),VSRC(NNDIM,NNVDIM),
     &TSRC(NNDIM,NNVDIM),SSRC(NNDIM,NNVDIM)
      REAL SRCRATE(NNDIM,NNVDIM)
C Local Variables
      REAL STRENGTH,RIVTEMP
      NEV=NNV-1
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign proper values
C   to arrays related to point sources.  The implementation of point
C   sources in the main program is node based in the horizontal and
C   element based in the vertical.  By default, the program assumes
C   that sources are added with momentum, temperature, and salinity
C   which are neutral with respect to the ambient conditions.
C
C Index range:  II - indicates the horizontal node(s) at which point 
C                  sources are desired.
C               KK - indicates the vertical elements(s) at which point 
C                  sources are desired.
C               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C               NEV - number of vertical elements (NEV=NNV-1)
C
C               II must be a member of the set given by the         
C               the horizontal node index: I=1,NN.
C
C               KK must be a member of the set given by the         
C               the vertical element index: K=1,NEV.
C
C 1. Assignments required whenever a point source is desired at 
C    horizontal node(s) II and vertical element(s) KK.
C
C    SRCRATE(II,KK) - Volume source strength (m^3/sec)
C                       at horizontal node II and vertical element KK.
C
C 2. Assignments required whenever a point source which is non-neutral
C    with respect to momentum is desired at horizontal node(s) II and
C    vertical element(s) KK.
C
C    USRC(II,NNV)=1.0 - Real flag to indicate that a source
C                       is being specified at horizontal node II
C                       which is non-neutral with respect to momentum.
C    USRC(II,KK)=UVAL - X-directed component of velocity for source
C                       at horizontal node II and vertical element KK.
C    VSRC(II,KK)=VVAL - Y-directed component of velocity for source
C                       at horizontal node II and vertical element KK.
C
C 3. Assignments required whenever a point source which is non-neutral
C    with respect to temperature is desired at horizontal node(s) II and
C    vertical element KK.
C
C    TSRC(II,NNV)=1.0 - Real flag to indicate that a source
C                       is being specified at horizontal node II
C                       which is non-neutral with respect to 
C                       temperature.
C    TSRC(II,KK)=TVAL - Temperature for source
C                       at horizontal node II and vertical element KK.
C
C 4. Assignments required whenever a point source which is non-neutral
C    with respect to salinity is desired at horizontal node(s) II and 
C    vertical element KK.
C
C    SSRC(II,NNV)=1.0 - Real flag to indicate that a source
C                       is being specified at horizontal node II
C                       which is non-neutral with respect to salinity.
C    SSRC(II,KK)=SVAL - Salinity for source
C                       at horizontal node II and vertical element KK.
C
C All related arrays are set to zero by the main routine before
C   each call to this routine.  Thus, the main program assumes that
C   there are no point sources within the domain.
C This subroutine is called once at the beginning of the simulation 
C   (with ITER=0) for initialization purposes; then it is called once
C   during each time step to set the desired point source information.
C================Beginning-of-user-specified-instructions===============

C...Merrimack River (node #1763)
      STRENGTH=313.2
      RIVTEMP= 9.7
      II=1763
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C...Saco River (node #2524)
      STRENGTH=116.6
      RIVTEMP=10.1
      II=2524
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C...Kennebec and Androscoggin Rivers (node #3614)
      STRENGTH=642.9
      RIVTEMP=9.2
      II=3614
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C...Penobscot River (node #6096)
      STRENGTH=649.8
      RIVTEMP=7.7
      II=6096
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C...St. Croix River (node #6534)
      STRENGTH=84.0
      RIVTEMP=9.0
      II=6534
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C...St. John River (node #6636)
      STRENGTH=1308.0
      RIVTEMP=5.9
      II=6636
      SRCRATE(II,NEV)=STRENGTH
      USRC(II,NNV)=1.0
      USRC(II,NEV)=0.0
      VSRC(II,NEV)=0.0
      TSRC(II,NNV)=1.0
      TSRC(II,NEV)=RIVTEMP
      SSRC(II,NNV)=1.0
      SSRC(II,NEV)=0.0

C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTGRIDQ5(NN,NNV,BATHY,ZETA,Z)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL BATHY(NNDIM),ZETA(NNDIM)
      REAL Z(NNDIM,NNVDIM)
C-----------------------------------------------------------------------
C This subroutine must be modified by the user to assign the vertical 
C   grid under each horizontal node.  There must be assignments to
C   specify the location of the nodes in the 3-D array Z(I,J):
C
C Index range:  I=1,NN, J=1,NNV
C          NN - number of horizontal nodes
C          NNV - number of vertical nodes
C
C Z and J are positive upward:
C          Bottom:  Z=-BATHY; J=1
C          Surface: Z= ZETA; J=NNV
C
C This subroutine is called multiple times during the initialization
C   phase for a simulation, then it is called once during each time step
C   to set the vertical grid as required by the solver.
C-----------------------------------------------------------------------
C Two standard subroutines are available to assign values to these
C   arrays:  Subroutines UNISIGMA and SINEGRID.
C
CCCCC
C
C   Subroutine UNISIGMAQ5 assigns uniform vertical spacing.
C   Assignment of NNV is required before this subroutine is called.
C
C      CALL UNISIGMAQ5(NNDIM,NN,NNV,BATHY,ZETA,Z)
C
CCCCC
C
C   Subroutine SINEGRIDQ5 assigns sinusiodal vertical spacing.
C   Assignment of the number of vertical nodes (NNV) and the minimum
C   node spacing at the vertical extremes of the water column (DZBL) is
C   required before this subroutine is called
C
C      CALL SINEGRIDQ5(NNDIM,NN,NNV,BATHY,ZETA,DZBL,Z)
C
C================Beginning-of-user-specified-instructions===============
C
      DZBL=1.0
      CALL SINEGRIDQ5(NNDIM,NN,NNV,BATHY,ZETA,DZBL,Z)
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE OUTPUTQ5(KD,SECKD,ITER,DELT,MESH,FILINQ,FILICQ,
     &Cd,AHI,NN,NNV,NE,XNOD,YNOD,IN,BATHY,akmin,
     &ZETA,UBAR,VBAR,Z,U,V,W,Q2,Q2L,RHO,T,S,ENZM,ENZH,ENZQ,
     &ZETAOLD,UBAROLD,VBAROLD,ZOLD,QP)
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C Fixed global variables
      REAL AHI(NNDIM),XNOD(NNDIM),YNOD(NNDIM),BATHY(NNDIM)
      INTEGER IN(3,NEDIM)
      REAL ZETA(NNDIM),UBAR(NNDIM),VBAR(NNDIM)
      REAL Z(NNDIM,NNVDIM),U(NNDIM,NNVDIM),V(NNDIM,NNVDIM),
     &W(NNDIM,NNVDIM),Q2(NNDIM,NNVDIM),Q2L(NNDIM,NNVDIM),
     &RHO(NNDIM,NNVDIM),T(NNDIM,NNVDIM),S(NNDIM,NNVDIM),
     &ENZM(NNDIM,NNVDIM),ENZH(NNDIM,NNVDIM),ENZQ(NNDIM,NNVDIM)
      REAL ZETAOLD(NNDIM),UBAROLD(NNDIM),VBAROLD(NNDIM),
     &ZOLD(NNDIM,NNVDIM),QP(NNDIM)

C cas implementing kestons tidal averaging
      integer iter, iterout 
      REAL Z_a(NNDIM,NNVDIM),ZOLD_a(NNDIM,NNVDIM),U_a(NNDIM,NNVDIM),
     &V_a(NNDIM,NNVDIM),W_a(NNDIM,NNVDIM),Q2_a(NNDIM,NNVDIM),
     &Q2L_a(NNDIM,NNVDIM),T_a(NNDIM,NNVDIM),S_a(NNDIM,NNVDIM)
      REAL ZETAOLD_a(NNDIM),UBAROLD_a(NNDIM),VBAROLD_a(NNDIM)
      REAL ZETA_a(NNDIM),UBAR_a(NNDIM),VBAR_a(NNDIM)

      save iterout
      save Z_a,ZOLD_a,U_a,V_a,W_a,Q2_a,Q2L_a,T_a,S_a
      save ZETAOLD_a,UBAROLD_a,VBAROLD_a
      save ZETA_a,UBAR_a,VBAR_a

      CHARACTER*72 MESH,FILINQ,FILICQ
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   BATHY(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) - sea 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),W(I,J) - nodal values of the X, Y, and Z components
C     of velocity.
C   Q2(I,J),Q2L(I,J) - nodal values of turbulent kinetic energy and the
C     turbulent kinetic energy times the master length scale.
C   RHO(I,J),T(I,J),S(I,J) - nodal values of density, temperature, and
C     salinity.
C   ENZM(I,L),ENZH(I,L),ENZQ(I,L) - elemental values of the vertical
C     diffusivities for momentum, mass variables, and the turbulent
C     variables.
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 ".icq5" HOT-START files:
C
C   ZETAOLD(I) - free surface elevation
C   UBAROLD(I),VBAROLD(I) - vertically averaged velocity
C   ZOLD(I,J) - nodal coordinate locations in 3-D
C
C Z and J are positive upward:
C          Bottom:  Z(I,1)=-BATHY(I)
C          Surface: Z(I,NNV)=ZETA(I)
C                   ZOLD(I,NNV)=ZETAOLD(I)
C
C The following are the parameters related to the simulation at hand;
C   which are particularly useful in the construction of ".icq5"
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 .icq5 file used to initialize the simulation
C   Cd - bottom slip coefficient; used with quadratic bottom slip
C        parameterization
C   AHI - nodal values of horizontal eddy viscosity/diffusivity
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 ICQ5WRITE' writes all computed results
C   to a user-specified file.
C
C      CALL ICQ5WRITE(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 Assignments for data input/output subroutines
      CHARACTER*72 OBSFILE,MESHL
      SAVE MESHL
c
c Assignments for SUBROUTINE DROG3DDT
      CHARACTER*72 INDFILE,GR2PATH
c
c Character assignment for SUBROUTINE ICQ5WRITE
      CHARACTER*72 FILOUT,INTERNAL,CITER,HEADER
cFEW24h added FILOUT24 (one line)
      CHARACTER*72 FILOUT24
      SAVE KDICQ5,SECKD0LD

C DJM
      integer idriftw
C End DJM

      IF(ITER.EQ.0)THEN
         WRITE(*,'(/A)')'SUBROUTINE OUTPUTQ5 INITIALIZATION:'
         OBSFILE='OBS/'
         WRITE(*,*)'Overwrite mass variables (1/0)?'
         READ(*,*)ANS
         IF(ANS.EQ.1)CALL OVERWRITEMASSVAR(OBSFILE,NN,NNV,Z,RHO,T,S,94)
         CALL CROPSTRING(MESH,istart,ibegin,iend)
         MESHL=MESH(ibegin:iend)
C DJM
         write(*,*) 'Are numerical drifters to reside at fixed depth (0),
     &               or move vertically (1)?'
         read(*,*) idriftw
C End DJM
      ENDIF
c
c Compare model results to those in current observation file and
c construct averages of friction and vertically integrated velocity over
c the period of the observations (i.e. output information required
c by TRUXTON6 for inversion purposes.
      OBSFILE='OBS/UV_INTEG.m3d'
      CALL NUMERICAL_UV_INTEG(KD,SECKD,ITER,DELT,
     &IN,XNOD,YNOD,BATHY,NN,NNV,NE,
     &UBAR,VBAR,Z,U,V,W,RHO,T,S,Q2,Q2L,Cd,ENZM,MESHL,OBSFILE,99,98)

      CALL CRUISE_AVERAGE(KD,SECKD,ITER,DELT,IN,XNOD,YNOD,NN,NNV,
     &NE,UBAR,VBAR,Z,U,V,Cd,ENZM,MESHL,OBSFILE,97)

c Write hourly .icq4 files [by commenting out only the third line]
c -OR-
c Write .icq4 file when ITER=11872 (i.e. 3 days into simulation when
c time-step=21.83203125sec) [by commenting out the first two lines]
c      IF (ITER.EQ.11872)THEN
       IF (ITER.EQ.1500) THEN
         FILOUT='ITER1500.icq5'
      CALL ICQ5WRITE(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

*      ITPERHR=NINT(3600.0/DELT)
*      IF(MOD(ITER,ITPERHR).EQ.0)THEN
*      IF(ITER.EQ.1440)THEN
*         WRITE(INTERNAL,'(I10)')100000+ITER
*         READ(INTERNAL,'(A)')CITER
*         CALL CSTRING(CITER,IS,IE)
*         FILOUT='ITER_'//CITER(IS+1:IE)//'.icq4'
*         CALL ICQ5WRITE(MESHL,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
c Track numerical drifters launched at times specified in INDFILE
      INDFILE='OBS/DRIFTERS.ind'
      GR2PATH=MESH
      CALL NUMERICAL_DRIFTERS(KD,SECKD,ITER,DELT,
     &NN,NE,NNV,Z,U,V,W,Q2,Q2L,RHO,T,S,ENZM,ENZH,ENZQ,INDFILE,GR2PATH,
     &idriftw)
c Output vertical structure of dependent variables at (x,y)
c locations specified in OBSFILE every hour
c      ITPERHR=NINT(3600.0/DELT)
c      IF(MOD(ITER,ITPERHR).EQ.0)THEN
c         OBSFILE='OBS/MOORING.nod'
c         CALL NUMERICAL_MOORINGS(NN,NNV,NE,XNOD,YNOD,Z,IN,U,V,W,
c     &   RHO,T,S,ENZM,ENZH,ENZQ,Q2,Q2L,KD,SECKD,ITER,OBSFILE,80)
c      ENDIF
c
c Interpolate model results to (t,x,y,z) locations specified in OBSFILE
*      OBSFILE='OBS/DRIFTERS.m4d'
*      CALL NUMERICAL_VPR(KD,SECKD,ITER,DELT,IN,XNOD,YNOD,BATHY,
*     &NN,NNV,NE,UBAR,VBAR,Z,U,V,W,RHO,T,S,Q2,Q2L,Cd,ENZM,MESH,
*     &OBSFILE,79)
c
c Output vertical transects of dependent variables for transects
c indicated in SUBROUTINE VQ5D every hour
*      ITPERHR=NINT(3600.0/DELT)
*      IF(MOD(ITER,ITPERHR).EQ.0)THEN
*         CALL VQ5D(NN,NNV,NE,XNOD,YNOD,Z,IN,U,V,W,RHO,T,S,
*     &   ENZM,ENZH,ENZQ,Q2,Q2L,KD,SECKD,ITER)
*      ENDIF
c
c Echo the iteration number and time to standard out every 10 timesteps
      IF(MOD(ITER,10).EQ.0)THEN
         CALL GDAY_GMT0(KD,SECKD,IYEAR,YDGMT0)
         WRITE(*,*)'ITER IYEAR YDGMT0 =',ITER,IYEAR,YDGMT0
      ENDIF

C cas - implemented Keston tidal cycle averaging to replace 24 hour
C averages
      IF (ITER.EQ.0) THEN
         do k=1,nn
          ZETA_a(k)=0.
          UBAR_a(k)=0.
          VBAR_a(k)=0.
          ZETAOLD_a(k)=0.
          UBAROLD_a(k)=0.
          VBAROLD_a(k)=0.
          do j=1,nnv
             U_a(k,j)=0.
             V_a(k,j)=0.
             W_a(k,j)=0.
             Z_a(k,j)=0.
             Zold_a(k,j)=0.
             Q2_a(k,j)=0.
             Q2L_a(k,j)=0.
             T_a(k,j)=0.
             S_a(k,j)=0.
          enddo
         enddo
       ENDIF
      iterout = 250
      RITEROUT=real(iterout)
      do k=1,nn
         ZETA_a(k)=ZETA_a(k)+ZETA(k)/RITEROUT
         UBAR_a(k)=UBAR_a(k)+UBAR(k)/RITEROUT
         VBAR_a(k)=VBAR_a(k)+VBAR(k)/RITEROUT
         ZETAOLD_a(k)=ZETAOLD_a(k)+ZETAOLD(k)/RITEROUT
         UBAROLD_a(k)=UBAROLD_a(k)+UBAROLD(k)/RITEROUT
         VBAROLD_a(k)=VBAROLD_a(k)+VBAROLD(k)/RITEROUT
         do j=1,nnv
             U_a(k,j)=U_a(k,j)+U(k,j)/RITEROUT
             V_a(k,j)=V_a(k,j)+V(k,j)/RITEROUT
             W_a(k,j)=W_a(k,j)+W(k,j)/RITEROUT
             Z_a(k,j)=Z_a(k,j)+Z(k,j)/RITEROUT
             Zold_a(k,j)=Zold_a(k,j)+Zold(k,j)/RITEROUT
             Q2_a(k,j)=Q2_a(k,j)+Q2(k,j)/RITEROUT
             Q2L_a(k,j)=Q2L_a(k,j)+Q2L(k,j)/RITEROUT
             T_a(k,j)=T_a(k,j)+T(k,j)/RITEROUT
             S_a(k,j)=S_a(k,j)+S(k,j)/RITEROUT
         enddo
      enddo

      IF (MOD(ITER,ITEROUT).EQ.0)THEN
          NSTF = ITER/ITEROUT 
          WRITE(INTERNAL,'(I10)')1000+NSTF
          READ(INTERNAL,'(A)')CITER
          CALL CSTRING(CITER,IS,IE)
          FILOUT='STF'//CITER(IS+1:IE)//'.icq5'
         CALL ICQ5WRITE(MESH,FILINQ,FILICQ,FILOUT,
     &   KD,SECKD,DELT,ZETA_a,UBAR_a,VBAR_a,ZETAOLD_a,
     &   UBAROLD_a,VBAROLD_a,Z_a,ZOLD_a,U_a,V_a,W_a,Q2_a,
     &   Q2L_a,T_a,S_a,NN,NNV,NNDIM,NNVDIM)
         do k=1,nn
          ZETA_a(k)=0.
          UBAR_a(k)=0.
          VBAR_a(k)=0.
          ZETAOLD_a(k)=0.
          UBAROLD_a(k)=0.
          VBAROLD_a(k)=0.
          do j=1,nnv
             U_a(k,j)=0.
             V_a(k,j)=0.
             W_a(k,j)=0.
             Z_a(k,j)=0.
             Zold_a(k,j)=0.
             Q2_a(k,j)=0.
             Q2L_a(k,j)=0.
             T_a(k,j)=0.
             S_a(k,j)=0.
          enddo
         enddo
      ENDIF

c Notify user that end of initialization has been reached if ITER=0
      IF(ITER.EQ.0)WRITE(*,'(A/)')'END USER SUBROUTINE INITIALIZATION:'
C
C End of routine
C
      RETURN
      END

C***********************************************************************
C***********************************************************************
      SUBROUTINE VQ5D(NN,NNV,NE,XNOD,YNOD,Z,IN,U,V,W,RHO,T,S,
     &ENZM,ENZH,ENZQ,Q2,Q2L,KD,SECKD,NQ5D)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C Include file for parameter statements
      INCLUDE 'DCMS.DIM'
C
C Global arrays
      real X(nndim),Y(nndim)
      real Z(nndim,nnvdim)
      real U(nndim,nnvdim),V(nndim,nnvdim),W(nndim,nnvdim)
      real RHO(nndim,nnvdim),T(nndim,nnvdim),S(nndim,nnvdim)
      real ENZM(nndim,nnvdim),ENZH(nndim,nnvdim),ENZQ(nndim,nnvdim)
      real Q2(nndim,nnvdim),Q2L(nndim,nnvdim)
      integer IN(3,nedim)
C
C Local declarations
      character*72 q4dfile,elefile,infile,cnq4d
C
C Determine character string
      write(infile,'(i10)')100000+nq4d
      read(infile,'(a)')cnq4d
      CALL CSTRING(cnq4d,istart,iend)
C
C Create v.q4d files
C  SF
*      px1=2.0779584e+05
*      py1=-1.7741785e+05
*      px2=2.1874687e+05
*      py2=-2.0986179e+05
*      nseg=20
*      q4dfile='ITER_'//cnq4d(istart+1:iend)//'_SF_xz.q4d'
*      elefile='ITER_SF_xz.ele'
*      CALL QUODDY2VSQ5(NN,NNV,NE,X,Y,Z,IN,U,V,W,RHO,T,S,
*     &ENZM,ENZH,ENZQ,Q2,Q2L,KD,SECKD,
*     &px1,py1,px2,py2,nseg,q4dfile,elefile)
C
C END OF ROUTINE
      return
      end

C***********************************************************************
C***********************************************************************
      SUBROUTINE SETQ5(T0,S0,G,OMEGA,NONLIN,NLBS,Cd,ISMAG,AH,AHMIN,
     &PRESSURE,MASSVAR,IHHBC,CLOSURE,Q20,Q2L0,
     &ZLOGBOT,ZLOGTOP,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,IUSTARQ2,
     &TAU0,EPSN,EPSH,EPSQ,THETA,AKMIN,Q2MIN,Q2LMIN,ELLMIN,
     &EKMMIN,EKHMIN,EKQMIN)
      CHARACTER*72 PRESSURE,MASSVAR,CLOSURE
C-----------------------------------------------------------------------
C purpose: This subroutine is provided to facilitate the overwriting of
C          parameters and switches set in the main program.
C
C algorithmic design: Any assignments made to passed arguments within 
C                     this subroutine will supersede the default 
C                     settings set in the main program.  Any changes
C                     will be reflected in the .echo file.
C                     
C passed arguments:
C   T0 - Reference temperature
C   S0 - Reference salinity
C   G - Gravitational acceleration
C   OMEGA - Radian frequency of celestial body (rad/sec)
C           (used to set the coriolis parameter)
C   NONLIN - integer flag to indicate whether or not nonlinear
C            advection is to be included in the momentum equation:
C              NONLIN=0 => no advection in momentum equation
C              NONLIN=1 => use advection in momemtum equation
C   NLBS - integer flag to indicate linearized bottom stress
C          coefficient (AK) parameterizaton in horizontal momentum
C          equations (bottom stress=AK*bottom velocity):
C              NLBS=0 => apply constant AK=AKMIN
C              NLBS=1 => apply quadratic stress condition
C                        => AK=Cd*bottom speed
C   Cd - bottom slip coefficient; used when NLBS=1
C   ISMAG - integer flag to indicate type of horizontal 
C           viscosity/diffusivity implemented within horizontal 
C           diffusion operator:
C              ISMAG=0 => apply uniform viscosity/diffusivity 
C                         (equal to AH)
C              ISMAG=1 => viscosity/diffusivity based on the horizontal 
C                         shear and grid scale (i.e. similar in spirit  
C                         to Smagorinsky(1963)), with multiplicative
C                         constant set to AH.
C   AH - real data utilized for horizontal viscosity/diffusivity closure
C        (see ISMAG for implementation)
C   AHMIN - minimum horizontal viscosity/diffusivity (only implemented 
C           when ISMAG=1)
C   PRESSURE - Character string to indicate pressure computation:
C               'BAROTROPIC' or 'BAROCLINIC'
C   MASSVAR - Character string to indicate treatment of mass variables:
C             'TWO-PROG' => Temperature and salinity transport equations
C                           are solved every time-step and 
C                           SUBROUTINE EQSTATE2_2D is used 
C                           to update RHO(T,S) after each time-step
C             'TWO-DIAG' => Initial values of temperature and salinity 
C                           are used by SUBROUTINE EQSTATE2_2D to set
C                           RHO(T,S), which is subsequently held
C                           constant throughout the simulation
C             'ONE-PROG' => Temperature transport equation is solved
C                           every time-step and SUBROUTINE EQSTATE1_2D 
C                           is used to update RHO(T) after each 
C                           time-step
C             'ONE-DIAG' => Initial values of temperature are used by
C                           SUBROUTINE EQSTATE1_2D to set RHO(T), 
C                           which is subsequently held
C                           constant throughout the simulation
C             'ZERO'=> No mass variables are implemented => temperature,
C                      salinity, and RHO are homogeneous and fixed
C   IHHBC - Binary switch to indicate BC type for temperature and
C           salinity at horizontal open boundaries:
C             1 => Fixed at the initial conditions
C             2 => No flux BC
C   CLOSURE - Character string to indicate turbulence closure:
C             'CONSTANT' => vertical eddy viscosity/diffusivities set
C                           to EKMIN
C             'MY25' => vertical eddy viscosity/diffusivities computed
C                       using Mellor and Yamada(1984,1982) turbulence
C                       closure with adjustments of Galperin et. al.
C                       (1988) and Blumberg et. al. (1992)
C   Q20 - Reference value of Q2 used to initialization
C   Q2L0 -  Reference value of Q2L used to initialization
C   ZLOGBOT - Height of the bottom node in the computational domain
C             above the true bottom
C   ZLOGTOP - Height of the top node in the computational domain
C             below the true surface
C   IQADVDIF - Binary switch to indicate treatment of advection and
C              horizontal diffusion in the Q2 and Q2L transport
C              equations:
C               0 => No advection and horizontal diffusion
C               1 => Include advection and horizontal diffusion
C   IQ2TBC  - Binary switch to indicate BC type for the Q2 transport
C             equation at the top of the water column:
C             1 => Type 1 BC given in Mellor and Yamada(1982)
C             2 => No flux BC
C   IQ2BBC  - Binary switch to indicate BC type for the Q2 transport
C             equation at the bottom of the water column:
C             1 => Type 1 BC given in Mellor and Yamada(1982)
C             2 => No flux BC
C   IQ2LTBC - Binary switch to indicate BC type for the Q2L transport
C             equation at the top of the water column:
C             1 => Type 1 BC given in Naimie(1995)[Ph.D. Thesis]
C             2 => No flux BC
C   IQ2LBBC - Binary switch to indicate BC type for the Q2L transport
C             equation at the bottom of the water column:
C             1 => Type 1 BC given in Naimie(1995)[Ph.D. Thesis]
C             2 => No flux BC
C   TAU0 - constant for weight applied to the continuity equation
C          when it is added to the vertically averaged momentum
C          equation to obtain the wave equation (typically .0002)
C   EPSN - implicit time weighting factor for vertical diffusion
C          of momentum:
C              EPSN=0.0 => explicit (at time level K)
C              EPSN=0.5 => implicit (at time level K+1/2)
C   EPSH - implicit time weighting factor for vertical diffusion
C          of temperature and salinity:
C              EPSH=0.0 => explicit (at time level K)
C              EPSH=0.5 => implicit (at time level K+1/2)
C   EPSQ - implicit time weighting factor for vertical diffusion
C          of Q2 and Q2L:
C              EPSQ=0.0 => explicit (at time level K)
C              EPSQ=0.5 => implicit (at time level K+1/2)
C   THETA - implicit gravity wave factor for the wave equation:
C              THETA=0.0 => explicit (at time level K-1)
C              THETA=1.0 => implicit (at time level K)
C   AKMIN - minimum linearized bottom stress coefficient
C   EKMIN - minimum vertical eddy diffusivity for momentum, heat,
C           and salinity
C   Q2MIN,Q2LMIN,ELLMIN - Minimum values of Q2, Q2L, and the mixing
C                         length.  Used only when CLOSURE='MY25'
C-----------------------------------------------------------------------
Ccas        EKHMIN = 1e-4
Ccas        EKMMIN = 1e-4
Ccas        EKQMIN = 1e-4
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SUBTITLEQ5(NFIL)
C-----------------------------------------------------------------------
C   Description: This user supplied subroutine is used to record
C                information on all the user provided subroutines
C                if so desired.  This routine can be left empty and
C                the execution of the QUODDY model is not affected.
C                If so desired, the provided information will be
C                written in the echo file (".echo") with the assigned
C                logical name NFIL.
C     Variables: [Argument list]
C                NFIL -- A pre-assigned logical file name of the echo
C                        file (".echo")
C                [Local character variables]
C                COMENT -- Comment on each user subroutine
C                FILSUB -- File contains all the user supplied routines
C      Accuracy: I*2, R*4 (single precision)
C         Input: NFIL
C        Output: Information on the user supplied subroutines is written
C                to the echo file (".echo")
C================Beginning-of-user-specified-instructions===============
C
C End of routine
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C End of User-Specified Subroutines for QUODDY5 (RELEASE 1.1)
C End of User-Specified Subroutines for QUODDY5 (RELEASE 1.1)
C***********************************************************************
C***********************************************************************
