C***********************************************************************
C***********************************************************************
      SUBROUTINE INITIALIZES3(MESHNAME,FILICQ4,NN,NNV,HDOWN,DELT,
     &             IDD0,IMM0,IYEAR0,SEC0,ZETAOLD,UOLD,VOLD,ZETAMID,
     &             UMID,VMID,ZMID,UZMID,VZMID)
C-----------------------------------------------------------------------
C
CASCO  For now we only allow COLD-START
C Objectives: Initialize dependent variable arrays using either
C             *COLD-START* or *HOT-START* information.  
C             For a *HOT-START*, conduct the following:
C               1. perform temporal interpolation if the time step for
C                  the .icq4 file does not match DELT for the current
C                  simulation
C               2. ensure the horizontal mesh matches that begin used
C                  for the current simulation
C               3. perform vertical interpolation from the .icq4 grid
C                  to the grid being utilized for the current simulation
C
C Complete Array Documentation: See main program QUODDY4_1.1
C
C Index Range: I=1,NN; J=1,NNV
C               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C
C Inputs: MESHNAME - the name of the geometric mesh
C         FILICQ4 - the name of the .icq4 hot-start file 
C                   (set to 'COLD-START' if a cold-start is desired).
C         NN,NNV,HDOWN - geometric mesh information
C         DELT - size of the time step
C         IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                                 of the simulation, under cold-start
C                                 conditions (these will be overwritten
C                                 if a hot-start is begin conducted)
C
C Outputs:IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                                 of the simulation (the input
C                                 values are only overwritten
C                                 if a hot-start is begin conducted)
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) - total depth
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) - total depth
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),WZMID(I,J) - x,y,z components of 
C                                            velocity
C
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'CASCO.DIM'

      INTEGER  NN, NNV
      INTEGER  IDD0,IMM0,IYEAR0
      REAL     SEC0
      REAL     DELT
C
C Global dimensioning and assignments
C
      REAL HDOWN(NNDIM)
      REAL ZETAOLD(NNDIM),UOLD(NNDIM),VOLD(NNDIM)
      REAL ZETAMID(NNDIM),UMID(NNDIM),VMID(NNDIM),ZMID(NNDIM,NNVDIM)
      REAL UZMID(NNDIM,NNVDIM),VZMID(NNDIM,NNVDIM)
      CHARACTER*72 MESHNAME,FILICQ4
C
C Local dimensioning and assignments
C
c     REAL ZICQ4I(NNVDIM),UI(NNVDIM),VI(NNVDIM),WI(NNVDIM),
c    &     Q2I(NNVDIM),Q2LI(NNVDIM),TI(NNVDIM),SI(NNVDIM)
c      CHARACTER*72 MESHICQ4

      INTEGER INDX, I, J
C
C START OF EXECUTABLE CODE
C
C Initialize 2-D horizontal nodal arrays, assuming a COLD-START
C
      DO I=1,NN
         ZETAMID(I)=0.0               
         UMID(I)=0.0          
         VMID(I)=0.0          
         ZETAOLD(I)=0.0
         UOLD(I)=0.0              
         VOLD(I)=0.0
      ENDDO
C
C Initialize 3-D nodal arrays, assuming a COLD-START
C
      DO J=1,NNV
         DO I=1,NN
            UZMID(I,J)=0.0
            VZMID(I,J)=0.0
         ENDDO
      ENDDO
C
C Initialize vertical grid, assuming a COLD-START
C
      CALL VERTGRIDS2(NN,NNV,HDOWN,ZETAMID,ZMID)
C
C Return if this is a COLD-START
C
      INDX=INDEX(FILICQ4,' ')-1
      IF(FILICQ4(1:INDX).EQ.'COLD-START')THEN
         RETURN
      ELSE
	 write(*,*)'Only allow cold starts'
	 write(2,*)'Only allow cold starts'
	 stop
      ENDIF
C
C Overwrite initialization arrays with those from a .icq4 file, 
C if a HOT-START is desired
C
C  First, Read .icq4 file
C
C     CALL ICQ4READ(FILICQ4,MESHICQ4,NNICQ4,NNVICQ4,
C    &DELTICQ4,IDD0,IMM0,IYEAR0,SEC0,HMID,UMID,VMID,HOLD,UOLD,VOLD,
C    &ZMID,ZOLD,UZMID,VZMID,WZMID,Q2MID,Q2LMID,TMID,SMID,NNDIM,NNVDIM)
C
C  Now the *0 timing parameters, the *OLD arrays, and the *MID arrays
C  contain information read from the .icq4 file.  Hence, temporal,
C  horizontal, and vertical interpolation/extrapolation may be
C  necessary to convert from the .icq4 file to the desired conditions
C  herein.
C
C  Temporal interpolation/extrapolation:
C
C   Compare dt passed from the main program to dt read from the .icq4
C   file and make adjustment on field variables at t-dt (i.e. the *OLD
C   arrays) if the two dt's are different (using linear
C   interpolation/extrapolation)
C
C     IF(ABS(DELT-DELTICQ4).GT.0.00001)THEN
C        WRITE(2,*)'TIME STEP (IN SECS) FOR ".inq" FILE  =',DELT
C        WRITE(2,*)'TIME STEP (IN SECS) FOR ".icq4" FILE =',DELTICQ4
C        RDELT=DELT/DELTICQ4
C        DO I=1,NN
C           HOLD(I)=HMID(I)-RDELT*(HMID(I)-HOLD(I))
C           UOLD(I)=UMID(I)-RDELT*(UMID(I)-UOLD(I))
C           VOLD(I)=VMID(I)-RDELT*(VMID(I)-VOLD(I))
C        ENDDO
C        DO J=1,NNVICQ4
C           DO I=1,NN
C              ZOLD(I,J)=ZMID(I,J)-RDELT*(ZMID(I,J)-ZOLD(I,J))
C           ENDDO
C        ENDDO
C        WRITE(2,*)'TEMPORAL INTERPOLATION OF THE INITIAL CONDITIONS SUC
C    &CESSFULLY COMPLETED'
C     ENDIF
C
C Horizontal interpolation/extrapolation:  Currently not implemented =>
C ensure MESHICQ4=MESHNAME and NNICQ4=NN.
C
C     ISTOP=0
C     INDX1=INDEX(MESHNAME,' ')-1
C     INDX2=INDEX(MESHICQ4,' ')-1
C     IF(MESHNAME(1:INDX1).NE.MESHICQ4(1:INDX2))THEN
C        ISTOP=1
C        WRITE(2,*)'INCONSISTENT MESH INPUTS:'
C        WRITE(2,*)' .inq FILE SPECIFIES '//MESHNAME(1:INDX1)
C        WRITE(2,*)' .icq4 FILE SPECIFIES '//MESHICQ4(1:INDX2)
C     ENDIF
C     IF(NNICQ4.NE.NN) THEN
C        ISTOP=1
C        WRITE(2,*)'INCONSISTENT HORIZONTAL NODE NUMBER INPUTS:'
C        WRITE(2,*)' .inq FILE SPECIFIES NN=',NN
C        WRITE(2,*)' .icq4 FILE SPECIFIES NN=',NNICQ4
C     ENDIF
C     IF(ISTOP.EQ.1)THEN
C        WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE ICQ4READ'
C        STOP
C     ENDIF
C
C Vertical interpolation/extrapolation:
C
C  Store desired grid at beginning of simulation in ZOLD
C
C     CALL VERTGRIDQ4(NN,NNV,HDOWN,HMID,ZOLD)
C
C  Interpolate quantities at beginning of simulation to the desired grid
C  and update *MID(I,J) values within a node loop
C
C     DO I=1,NN
C        DO JICQ4=1,NNVICQ4
C           ZICQ4I(JICQ4)=ZMID(I,JICQ4)
C        ENDDO
C        DO J=1,NNV
C           CALL BASIS1D1(NNVICQ4,ZICQ4I,ZOLD(I,J),b2,n2,b1,n1)
C
C    notes:   ZOLD(I,J)<ZICQ4I(1) => n1=nn,b1=0.0,n2=1,b2=1.0 
C             ZOLD(I,J)>ZICQ4I(NNVICQ4) => n1=nn,b1=1.0,n2=1,b2=0.0 
C
C           UI(J)=b1*UZMID(I,n1)+b2*UZMID(I,n2)
C           VI(J)=b1*VZMID(I,n1)+b2*VZMID(I,n2)
C           WI(J)=b1*WZMID(I,n1)+b2*WZMID(I,n2)
C           Q2I(J)=b1*Q2MID(I,n1)+b2*Q2MID(I,n2)
C           Q2LI(J)=b1*Q2LMID(I,n1)+b2*Q2LMID(I,n2)
C           TI(J)=b1*TMID(I,n1)+b2*TMID(I,n2)
C           SI(J)=b1*SMID(I,n1)+b2*SMID(I,n2)
C        ENDDO
C        DO J=1,NNV
C           ZMID(I,J)=ZOLD(I,J)
C           UZMID(I,J)=UI(J)
C           VZMID(I,J)=VI(J)
C           WZMID(I,J)=WI(J)
C           Q2MID(I,J)=Q2I(J)
C           Q2LMID(I,J)=Q2LI(J)
C           TMID(I,J)=TI(J)
C           SMID(I,J)=SI(J)
C        ENDDO
C     ENDDO
C
C Now, Set the desired z-coordinates for the time-step just prior
C to the beginning of this simulation => ZOLD
C
C     CALL VERTGRIDQ4(NN,NNV,HDOWN,HOLD,ZOLD)
C
C END OF ROUTINE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE INPUTDECS3(MESHNAME,FILBEL,THETA,
     &                     XFCTR,YFCTR,HFCTR,HMIN,
     &                     NN,NE,NNV,X,Y,IN,HDOWN,SUMX,SUMY,
     &                     DSBDRY,INBE,NBE,INGEO,NGEO,NORM,NVN,
     &                     KST,LRVS,NR,NPCN,LST,NHBW,BCFTR)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via an INCLUDE file
C
      INCLUDE 'CASCO.DIM'
C
C Global dimensioning and assignments
C
      INTEGER NN,NE,NNV
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM)
      INTEGER LRVS(NNDIM),NORM(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM)
      REAL DSBDRY(NNDIM),BCFTR(NNDIM)
      CHARACTER*72 MESHNAME,FILBEL
      REAL XFCTR,YFCTR,HFCTR,HMIN
      INTEGER KST,NR,NPCN,LST, NHBW 
      INTEGER NBE, NGEO, NVN
      REAL    THETA
      
C
C Local dimensioning and assignments
C
      CHARACTER*72 FILBAT,FILELE,FILNOD
      CHARACTER*72 CLINES
c     CHARACTER*72 UNITS
c     INTEGER LRT(20),NCD(NNDIM)
      INTEGER NCD(NNDIM)

      INTEGER I,J,K
      INTEGER N1,N2
      INTEGER NBW
      INTEGER IMESH, NK, JBAT, NB
      INTEGER NBC1, NCDNK

      REAL  DELX, DELY, SX,SY 


      REAL ZERO
      DATA ZERO/0.0/
C
C Define the LABLIB standard node coordinate, element incidence,
C bathymetric depth, and boundary element files
C
      IMESH=INDEX(MESHNAME,' ')-1
      FILNOD=MESHNAME(1:IMESH)//'.nod'
      FILELE=MESHNAME(1:IMESH)//'.ele'
      FILBAT=MESHNAME(1:IMESH)//'.bat'
C
C Read the node coordinate list in MKS units (the ".nod" file)
C
      OPEN(UNIT=3,FILE=FILNOD,STATUS='OLD')
      NN=0
 10   READ(3,*,END=11)J,X(J),Y(J)
      NN=NN+1
      X(J)=XFCTR*X(J)
      Y(J)=YFCTR*Y(J)
      GO TO 10
 11   CLOSE(3)
C
C Read the element incidence list (the ".ele" file)
C                
      OPEN(UNIT=3,FILE=FILELE,STATUS='OLD')
      NE=0
 20   READ(3,*,END=21)J,(IN(I,J),I=1,3)
      NE=NE+1
      GO TO 20
 21   CLOSE(3)
C
C
C Check dimensioning of NN, NE, NNV, and NLEV
C
      IF(NNDIM.GE.NN .AND. NEDIM.GE.NE .AND. NNVDIM.GE.NNV)GO TO 25
      WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECQ4 ', 
     &  'DUE TO IN SUFFICIENT DIMENSIONS:'
      WRITE(2,*)' NN   =',NN,   ',NNDIM   =',NNDIM
      WRITE(2,*)' NE   =',NE,   ',NEDIM   =',NEDIM
      WRITE(2,*)' NNV  =', NNV, ',NNVDIM  =',NNVDIM
      STOP
   25 CONTINUE
C
C Compute bandwidth from the element incidence list and terminate
C execution if dimensioning parameter is to small
C
      CALL HALFBW(IN,3,NE,NHBW)
      IF(THETA.EQ.ZERO)NHBW=0
      NBW=2*NHBW+1
      IF(NBW.GT.NBWE)THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECQ4 ',
     &   'DUE TO INSUFFICIENT DIMENSIONS:'
         WRITE(2,*)' NBW   =',NBW,   ',NBWE   =',NBWE
         STOP
      END IF
C
C Read the bathymetric depth in MKS units (the ".bat" file) and check
C to ensure the number of entries matches the number of nodes
C
      OPEN(UNIT=3,FILE=FILBAT,STATUS='OLD')
      JBAT=0
   30 READ(3,*,END=31)J,HDOWN(J)
      JBAT=JBAT+1
      HDOWN(J)=HFCTR*HDOWN(J)
      IF(HDOWN(J).LT.HMIN)HDOWN(J)=HMIN
      GO TO 30
   31 CLOSE (3)
      IF(JBAT.NE.NN)THEN
      WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECS2 ', 
     &  'DUE TO IN COMPATABILITY BETWEEN THE NUMBER OF NODES IN ',
     &  ' .nod AND .bat FILES:'
         WRITE(2,*)' '//MESHNAME(1:IMESH)//'.nod CONTAINS',NN,' NODES'
         WRITE(2,*)' '//MESHNAME(1:IMESH)//'.bat CONTAINS',JBAT,' NODES'
         STOP
      ENDIF
C
C Read the boundary element incidence list, the ".bel" file,
C load geostrophic boundary element incidence list,
C and check dimensioning size 
C                
      OPEN(UNIT=3,FILE=FILBEL,STATUS='OLD')
      NBE=0
      NGEO=0
C
C  Read and ignore the first two lines of header
C
      READ(3,'(A)')CLINES
      READ(3,'(A)')CLINES
C
C  Read the boundary element information
C
CGH Dec 7 '98  Replace with code which checks the order of the 
CGH            boundary element incidence list.
CGH            Taken from INPUTDECF5 by CEN 
c40   CONTINUE 
c     READ(3,*,END=41)J,INBE(1,J),INBE(2,J),INBE(3,J),INBE(4,J)
c     NBE=NBE+1
c     IF(INBE(4,J).EQ.4)THEN
c        NGEO=NGEO+1
c        INGEO(1,NGEO)=INBE(1,J)
c        INGEO(2,NGEO)=INBE(2,J)
c     ENDIF
c     GO TO 40
c41   CLOSE(3)
c     IF(NBE.GT.NBEDIM)THEN
c        WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECS2 ',
c    &             'DUE TO INSUFFICIENT DIMENSIONS:'
c        WRITE(2,*)' NBE   =',NBE,   ',NBEDIM   =',NBEDIM
c        STOP
c     END IF
C
C Read the boundary element information
C
 40   READ(3,*,END=41)J,INBE(1,J),INBE(2,J),INBE(3,J),INBE(4,J)
      NBE=NBE+1
      GO TO 40
 41   CLOSE(3)
      IF(NBE.GT.NBEDIM)THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECS3 DUE TO
     & INSUFFICIENT DIMENSIONS:'
         WRITE(2,*)' NBE   =',NBE,   ',NBEDIM   =',NBEDIM
         STOP
      END IF
C
C Reorder boundary element list and check for validity
C
      CALL REORDERINBE(NNDIM,NBEDIM,X,Y,INBE,NBE)
C
C Load geostrophic boundary element incidence list
C
      DO J=1,NBE
         IF(INBE(4,J).EQ.4)THEN
            NGEO=NGEO+1
            INGEO(1,NGEO)=INBE(1,J)
            INGEO(2,NGEO)=INBE(2,J)
         ENDIF
      ENDDO
C
C
CGH end of Dec 7 '98 change ^^^
C
C Determine node codes from boundary element incidence list and
C establish LRVS and NORM with their respective pointers.
C
C    NORM(I) - boundary node numbers where normal velocity is set
C                I =1,NVN       are V.n .NE. 0.0    (NCD=3)
C                I =KST,NNDIM   are V.n=0.0         (NCD=1,2)
C    LRVS(I) - node numbers of set depth and set velocity nodes
C                I=1,NR         are set depth nodes (NCD=5,6,9)
C                I=NR+1,NR+NPCN are V=0.0           (NCD=10)
C                I=LST,NNDIM    are V .NE. 0.0      (NCD=8,11)
C
      NVN=0
      KST=NNDIM+1
      NR=0
      NPCN=0
      LST=NNDIM+1
      DO 50 I=1,NNDIM
         NCD(I)=0
         NORM(I)=0
         LRVS(I)=0
 50   CONTINUE
      DO 60 J=1,NBE
         do 60 K=1,2
            NK=INBE(K,J)
            NCDNK=NCD(NK)
            NB=INBE(4,J)
            if(NCDNK.EQ.0)then
               NCD(NK)=NB
            else if(NCDNK.EQ.1)then
               if(NB.EQ.1)then
                  NCD(NK)=1
                  KST=KST-1
                  NORM(KST)=NK                 
               else if(NB.EQ.3)then
                  NCD(NK)=8
                  LST=LST-1
                  LRVS(LST)=NK
               else if(NB.EQ.4)then
                  NCD(NK)=7
                  KST=KST-1
                  NORM(KST)=NK
               else if(NB.EQ.5)then
                  NCD(NK)=6
                  KST=KST-1
                  NORM(KST)=NK
                  NR=NR+1
                  LRVS(NR)=NK
               else
               endif
            else if(NCDNK.EQ.2)then
               if(NB.EQ.2)then
                  NCD(NK)=2
                  KST=KST-1
                  NORM(KST)=NK                 
               else if(NB.EQ.3)then
                  NCD(NK)=8
                  LST=LST-1
                  LRVS(LST)=NK
               else if(NB.EQ.4)then
                  NCD(NK)=7
                  KST=KST-1
                  NORM(KST)=NK
               else if(NB.EQ.5)then
                  NCD(NK)=6
                  KST=KST-1
                  NORM(KST)=NK
                  NR=NR+1
                  LRVS(NR)=NK
               else
               endif
            else if(NCDNK.EQ.3)then
               if(NB.EQ.1.OR.NB.EQ.2)then
                  NCD(NK)=8
                  LST=LST-1
                  LRVS(LST)=NK
               else if(NB.EQ.3)then
                  NCD(NK)=3
                  NVN=NVN+1
                  NORM(NVN)=NK
               else
               endif
            else if(NCDNK.EQ.4)then
               if(NB.EQ.1.OR.NB.EQ.2)then
                  NCD(NK)=7
                  KST=KST-1
                  NORM(KST)=NK
               else if(NB.EQ.4)then
                  NCD(NK)=4
               else if(NB.EQ.5)then
                  NCD(NK)=9
                  NR=NR+1
                  LRVS(NR)=NK
               else
               endif
            else if(NCDNK.EQ.5)then
               if(NB.EQ.1.OR.NB.EQ.2)then
                  NCD(NK)=6
                  KST=KST-1
                  NORM(KST)=NK
                  NR=NR+1
                  LRVS(NR)=NK
               else if(NB.EQ.4)then
                  NCD(NK)=9
                  NR=NR+1
                  LRVS(NR)=NK
               else if(NB.EQ.5)then
                  NCD(NK)=5
                  NR=NR+1
                  LRVS(NR)=NK
               else
               endif
            else
            endif
 60   CONTINUE
C
C Check dimensioning size for NR
C
      IF(NR.GT.NRDIM)THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECQ4 DUE ',
     &             'TO INSUFFICIENT DIMENSIONS:'
         WRITE(2,*)' NR   =',NR,   ',NRDIM   =',NRDIM
         STOP
      END IF
C
C Determine unit normals for boundary nodes with at least one dry 
C boundary segment
C
      DO 70 I=1,NN
         SUMX(I)=0.0
         SUMY(I)=0.0
         DSBDRY(I)=0.0
 70   CONTINUE
      DO 71 J=1,NBE
         N1=INBE(1,J)
         N2=INBE(2,J)
         IF(INBE(4,J).EQ.1.OR.INBE(4,J).EQ.2.OR.INBE(4,J).EQ.3)THEN
            DELX=X(N2)-X(N1)
            DELY=Y(N2)-Y(N1)
            SUMX(N1)=SUMX(N1)+DELX
            SUMY(N1)=SUMY(N1)+DELY
            SUMX(N2)=SUMX(N2)+DELX
            SUMY(N2)=SUMY(N2)+DELY
         ENDIF
 71   CONTINUE
      DO 72 J=1,NBE
         N1=INBE(1,J)
         IF(NCD(N1).EQ.1.OR.NCD(N1).EQ.2.OR.NCD(N1).EQ.3
     &  .OR.NCD(N1).EQ.6.OR.NCD(N1).EQ.7.OR.NCD(N1).EQ.8)THEN
            SX=SUMX(N1)
            SY=SUMY(N1)
            DSBDRY(N1)=SQRT(SX*SX+SY*SY)
            SUMX(N1)=-SX/DSBDRY(N1)
            SUMY(N1)=SY/DSBDRY(N1)
         ENDIF
 72   CONTINUE
C                
C Load the BCFTR array to give the correct Coriolis factor of either
C   BCFTR(I) = 0 for solid boundary node with V.n = 0,
C   BCFTR(I) = 1 otherwise;for the (u + j*v) solution calculations
C
      DO 80 I=1,NNDIM
         BCFTR(I)=1.
 80   CONTINUE
      DO 82 I=KST,NNDIM
         NBC1=NORM(I)
         BCFTR(NBC1)=0.
 82   CONTINUE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE INPUTS2(FILINQ,
     &FILICQ,IDD0,IMM0,IYEAR0,SEC0,IDDL,IMML,IYEARL,SECL,
     &DEG,DELT,NNV,THETA,TAU0,NLBS,AKMIN,EKMIN,EPSN,
     &MESHNAME,NN,NE,NHBW,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,
     &             BCFTR)
C-----------------------------------------------------------------------
C Objectives: Read .ins file
C             Open .echo file
C             Read mesh information from files indicated in .ins file
C             Echo mesh information to the .echo file
C
C Complete Array Documentation: See main program QUODDY4_1.1
C
C Index Range: I=1,NN; J=1,NNV; JL=1,NLEV; K1=1,NBE; K2=1,NGEO; L=1,NE; 
C              M=1,NR; NB1=1,NHBW+1; NB2=1,2*NHBW+1
C               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C               NBE - number of boundary elements
C               NGEO - number of geostrophic boundary elements
C               NE - number of horizontal elements
C               NR - number of elevation boundary conditions
C               NHBW - half bandwidth of horizontal mesh
C
C Inputs: FILINQ - the name of the .inq file for this simulation
C
C Outputs:
C
C       Simulation Control Parameters:
C
C         FILICQ - name of hot start .icq4 file
C                  (set to 'COLD-START' if a cold start is desired)
C         IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                                 of the simulation (only set set if 
C                                 a cold start simulation is being
C                                 conducted)
C         IDDL,IMML,IYEARL,SECL - timing parameters for the end
C                                 of the simulation
C         DEG - degrees latitude to be used for coriolis terms (f-plane)
C         DELT - time step
C         NNV - number of vertical nodes
C         THETA - implicit gravity wave factor for the wave equation:
C                    THETA=0.0 => explicit
C                    THETA=1.0 => implicit
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         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 => read AK from a file 
C         AKMIN - minimum linearized bottom stress coefficient
C         EKMIN - minimum vertical eddy diffusivity for momentum, heat,
C                 and salinity
C         EPSN - implicit time weighting factor for vertical diffusion
C                of momentum:
C                    EPSN=0.0 => explicit
C                    EPSN=0.5 => Crank-Nicholson
C                    EPSN=1.0 => implicit
C
C      Horizontal Mesh Information:
C
C         MESHNAME - name of the horizontal mesh
C         NN - number of horizontal nodes
C         NE - the number of horizontal 2-D triangular elements
C         NHBW - half bandwidth of horizontal finite element mesh
C         X(I),Y(I) - cartesian nodal coordinates
C         IN(3,L) - LABLIB STANDARD triangular element incidence list
C         HDOWN(I) - Bathymetric depth at nodes
C         SUMX(I) - cosine value for boundary nodes where V.n is 
C                 specified
C         SUMY(I) - sine value for boundary nodes where V.n is 
C                 specified
C         DSBDRY(I) - surface integral information
C         INBE(2,K) - LABLIB STANDARD boundary element incidence list
C         NBE - number of boundary elements
C         INGEO(2,K) - LABLIB STANDARD geostrophic boundary element
C                 incidence list
C         NGEO - number of geostrophic boundary elements
C         NORM(I) - boundary node numbers where normal velocity is set
C                 I =1,NVN       are V.n .NE. 0.0    (NCD=3)
C                 I =KST,NNDIM   are V.n=0.0         (NCD=1,2)
C         LRVS(I) - node numbers of set depth and set velocity nodes
C                 I=1,NR         are set depth nodes (NCD=5,6,9)
C                 I=NR+1,NR+NPCN are V=0.0           (NCD=10)
C                 I=LST,NNDIM    are V .NE. 0.0      (NCD=8,11)
C
C
C       Other Information:
C
C         BCFTR(I)  - Real binary multiplicative factor which deletes
C                     the Coriolis term in momentum equation at solid
C                     boundaries
C                 BCFTR(I)=0.0 for land boundary nodes
C                 BCFTR(I)=1.0 for all other nodes
C
C History:
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'CASCO.DIM'

      integer NE,NN,NNV, NHBW
      integer NBE, NGEO, NVN, LST,  KST, NR, NPCN

      INTEGER  IDD0,IMM0,IYEAR0
      INTEGER  IDDL,IMML,IYEARL
      REAL     SEC0, SECL
      REAL     DELT
      REAL     DEG 
C      REAL  COR, G, RAD, OMEGA2
      REAL  AKMIN, EKMIN
      REAL  EPSN, THETA, TAU0
      INTEGER  NLBS


C
C Global dimensioning and assignments
C
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM)
      REAL SUMX(NNDIM),SUMY(NNDIM),DSBDRY(NNDIM),BCFTR(NNDIM)
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM)
      INTEGER NORM(NNDIM),LRVS(NNDIM)
      CHARACTER*72 FILINQ,FILICQ,MESHNAME
C
C Local dimensioning and assignments
C
      CHARACTER*72 INQLABEL,COMINQ,CHARINQ,FILBEL,FILECH,UNITS
      INTEGER      IECH, IICQ

      REAL  HMIN, XFCTR,YFCTR, HFCTR

C
C Open and read input file: FILINQ
C
      OPEN(UNIT=1,FILE=FILINQ,STATUS='OLD')
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')COMINQ
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')MESHNAME
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')FILBEL
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')FILICQ
C
C If a COLD START is desired, backspace one line and read starting 
C time parameters. 
C
      IICQ=INDEX(FILICQ,' ')-1
      if(FILICQ(1:IICQ).EQ.'COLD-START')then
         backspace(1)
         READ(1,'(a,1x,i2,1x,i2,1x,i4,1x,f21.12)')
     &       CHARINQ(1:IICQ),IDD0,IMM0,IYEAR0,SEC0
      endif
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')FILECH
      IECH=INDEX(FILECH,' ')-1
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')UNITS
      READ(1,*)XFCTR,YFCTR,HFCTR
      READ(1,*)DEG
      READ(1,*)HMIN
      read(1,*)IDDL,IMML,IYEARL,SECL
      READ(1,*)DELT
      READ(1,*)THETA
      READ(1,*)TAU0
      READ(1,*)NNV
      READ(1,*)NLBS
      READ(1,*)EKMIN
      READ(1,*)AKMIN
      READ(1,*)EPSN
      CLOSE(1)
C
C Open echo file
C
      OPEN(UNIT=2,FILE=FILECH,STATUS='UNKNOWN')
      REWIND(2)
C
C Make sure SI UNITS are used in the simulation
C
      IF(UNITS(1:8).NE.'SI UNITS')THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTQ4'
         WRITE(2,*)'INVALID UNITS, MODEL ONLY ALLOWS SI UNITS'
         STOP
      ENDIF
C
C Load geometry and boundary condition data
C
      CALL INPUTDECS3(MESHNAME,FILBEL,THETA,
     &                     XFCTR,YFCTR,HFCTR,HMIN,
     &                     NN,NE,NNV,X,Y,IN,HDOWN,SUMX,SUMY,
     &                     DSBDRY,INBE,NBE,INGEO,NGEO,NORM,NVN,
     &                     KST,LRVS,NR,NPCN,LST,NHBW,BCFTR)
C
C Echo input information to echo file
C
      CALL ECHOINPUTS2(COMINQ,MESHNAME,FILBEL,FILINQ,FILICQ,FILECH,
     &          XFCTR,YFCTR,HFCTR,DEG,HMIN,DELT,THETA,X,Y,HDOWN,
     &          TAU0,NNV,NLBS,EKMIN,AKMIN,EPSN,
     &          NN,NE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NHBW)

      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Echo file writing:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE ECHOINPUTS2(COMINQ,MESHNAME,FILBEL,
     &          FILINQ,FILICQ,FILECH,
     &          XFCTR,YFCTR,HFCTR,DEG,HMIN,DELT,THETA,X,Y,HDOWN,
     &          TAU0,NNV,NLBS,EKMIN,AKMIN,EPSN,
     &          NN,NE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NHBW)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via an INCLUDE file
C
      INCLUDE 'CASCO.DIM'
C
C Global dimensioning and assignments
C
      REAL XFCTR,YFCTR,HFCTR,DEG,HMIN,DELT,THETA
      REAL TAU0,EKMIN,AKMIN,EPSN
      INTEGER NNV, NN,NE,NBE,NGEO,NVN,KST,NR,NPCN,NHBW,LST, NLBS 
      INTEGER INGEO(2,NBEDIM),NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM)
      CHARACTER*72 COMINQ,MESHNAME,FILBEL,FILINQ,FILICQ,FILECH
C
C Declaration of local variables
C
      CHARACTER*15 CMODEL,FILDIM
      CHARACTER*15 finfo 

      REAL YMAX, YMIN, XMAX, XMIN
      REAL BMIN, BMAX
      INTEGER INDX, I, IMESH 
C
C Define model title, driver program, and user provided subroutine
C
      CMODEL=' CASCO '
      FILDIM='CASCO.DIM'
C
C Write output title to the echo file CASENAME.ech (unit 2)
C
      call casco_info(2) 
C
C
C Echo the one line comment from FININQ
C
      INDX=INDEX(FILINQ,' ')- 1
      WRITE(2,600)' '
      WRITE(2,600)'Comment from '//FILINQ(1:INDX)//':'
      WRITE(2,600)COMINQ
C
C Determine the exact length of the casename
C
      IMESH=INDEX(MESHNAME,' ')- 1
C
C Echo all the input filenames
C
      WRITE(2,600)' '
      WRITE(2,600)'Case name:          '//MESHNAME(1:IMESH)
      WRITE(2,600)'Input file:         '//FILINQ(1:INDX)
      INDX=INDEX(FILDIM,' ')- 1
      WRITE(2,600)'Dimensioning file:  '//FILDIM(1:INDX)
      INDX=INDEX(FILICQ,' ')- 1
      WRITE(2,600)'Initial condition:  '//FILICQ(1:INDX)
      INDX=INDEX(FILECH,' ')- 1
      WRITE(2,600)'Echo file:          '//FILECH(1:INDX)
      WRITE(2,600)'Node list:          '//
     &                 MESHNAME(1:IMESH)//'.nod'
      WRITE(2,600)'Element list:       '//
     &                 MESHNAME(1:IMESH)//'.ele'
      WRITE(2,600)'Bathymetry file:    '//
     &                 MESHNAME(1:IMESH)//'.bat'
      INDX=INDEX(FILBEL,' ')- 1
      WRITE(2,600)'Bndy element file:  '//FILBEL(1:INDX)
C
C Write information as directed in SUBROUTINE SUBTITLE
C
      WRITE(2,600)' '
      CALL SUBTITLE(2)
C
C Write the units used in the simulation
C
      WRITE(2,600)' '
      WRITE(2,600)'Units: SI UNITS'
C
C Echo program parameters
C
      WRITE(2,600)' '
      WRITE(2,600)'Program parameters:'
      WRITE(2,609)'NNDIM  ', NNDIM
      WRITE(2,609)'NEDIM  ', NEDIM
      WRITE(2,609)'NNVDIM ', NNVDIM
      WRITE(2,609)'NFTRDIM', NFTRDIM
      WRITE(2,609)'NBEDIM ', NBEDIM
      WRITE(2,609)'NRDIM  ', NRDIM
      WRITE(2,609)'NBWE   ', NBWE
      WRITE(2,609)'NB1DIM ', NB1DIM
C
C Echo geometric parameters
C
      WRITE(2,600)' '
      WRITE(2,600)'Geometric parameters:'
      WRITE(2,610)'X scaling factor                    ', XFCTR
      WRITE(2,610)'Y scaling factor                    ', YFCTR
      WRITE(2,610)'Z scaling factor                    ', HFCTR
      WRITE(2,610)'Degrees latitude                    ', DEG
      WRITE(2,610)'Minimum depth                       ', HMIN
C
C Echo horizontal finite element data
C
      WRITE(2,600)' '
      WRITE(2,600)'Horizontal finite element mesh data:'
      WRITE(2,611)'Total number of nodes              ', NN
      WRITE(2,611)'Total number of elements           ', NE
      WRITE(2,611)'Total number of boundary elements  ', NBE
      WRITE(2,611)'Half bandwidth required            ', NHBW
      WRITE(2,600)' '
      XMIN=X(1)
      XMAX=X(1)
      YMIN=Y(1)
      YMAX=Y(1)
      BMIN=HDOWN(1)
      BMAX=HDOWN(1)
      DO I=2,NN
         XMIN=AMIN1(XMIN,X(I))
         XMAX=AMAX1(XMAX,X(I))
         YMIN=AMIN1(YMIN,Y(I))
         YMAX=AMAX1(YMAX,Y(I))
         BMIN=AMIN1(BMIN,HDOWN(I))
         BMAX=AMAX1(BMAX,HDOWN(I))
      ENDDO
      WRITE(2,616)'Range of X (MIN,MAX)               ',XMIN,XMAX
      WRITE(2,616)'Range of Y (MIN,MAX)               ',YMIN,YMAX
      WRITE(2,616)'Range of bathymetry (MIN,MAX)      ',BMIN,BMAX
C
C Output vertical mesh data
C
      WRITE(2,600)' '
      WRITE(2,600)'Vertical finite element mesh data:'
      WRITE(2,611)'Number of vertical nodes            ', NNV
C
C Output horizontal boundary conditions
C
      WRITE(2,600)' '
      WRITE(2,600)'Wave equation and momentum equation horizontal '//
     &             'BC totals:'
      WRITE(2,615)'Totals for boundary conditions to '//
     &              'be set in subroutine BC:'
      WRITE(2,612)'Set depth nodes                    ', NR
      WRITE(2,612)'Set nonzero normal velocity        ', NVN
      WRITE(2,612)'Set nonzero velocity               ', NNDIM+1-LST
      WRITE(2,615)'Totals for boundary conditions computed '//
     &              'by CASCO:'
      WRITE(2,612)'Zero normal velocity nodes         ', NNDIM+1-KST
      WRITE(2,612)'Zero velocity nodes                ', NPCN
      WRITE(2,612)'Geostrophic boundary elements      ', NGEO
      IF (NR .GT. 0)THEN
         WRITE(2,3070)
         WRITE(2,3090)(LRVS(I),I=1,NR)
      ENDIF
      IF (NVN .GT. 0)THEN
         WRITE(2,3120)
         WRITE(2,3090)(NORM(I),I=1,NVN)
      ENDIF
      IF (NNDIM+1-KST .GT. 0)THEN
        WRITE(2,3130)
        WRITE(2,3090)(NORM(I),I=KST,NNDIM)
      ENDIF
      IF (NPCN .GT. 0)THEN
         WRITE(2,3140)
         WRITE(2,3090)(LRVS(I),I=NR+1,NR+NPCN)
      ENDIF
      IF (NNDIM+1-LST .GT. 0)THEN
         WRITE(2,3150)
         WRITE(2,3090)(LRVS(I),I=LST,NNDIM)
      ENDIF
      IF (NGEO .GT. 0)THEN
         WRITE(2,3160)
         DO I=1,NGEO
            WRITE(2,3170)INGEO(1,I),INGEO(2,I)
         ENDDO
      ENDIF
C
C Output format statements
C
  600 FORMAT(A)
  601 FORMAT(A72)
  602 FORMAT(A2,68X,A2)
  603 FORMAT(A2,26X,A16,26X,A2)
  604 FORMAT(A2,18X,A32,18X,A2)
  605 FORMAT(A2,14X,A40,14X,A2)
  606 FORMAT(A2,19X,A30,19X,A2)
  609 FORMAT(5X,A7,6X,I6)
  610 FORMAT(5X,A35,1X,E12.4)
  611 FORMAT(5X,A35,1X,I12)
  612 FORMAT(10X,A30,1X,I12)
  613 FORMAT(A,I7,A)
  614 FORMAT(A,I7,5X,A,I7,5X,A,I7)
  615 FORMAT(5X,A)
  616 FORMAT(5X,A35,1X,E12.4,1X,E12.4)
  618 FORMAT(5X,A35,1X,I2,1X,I2,3X,I4,1X,F18.8)
 2006 FORMAT(4A4,4X,2G10.0)
 2009 FORMAT(11X,'INSUFFICIENT DATA FOR NORMAL DIRECTIONS')
 2010 FORMAT(13I6,A2)
 2011 FORMAT(I5,G10.0,63X,A2)
 3030 FORMAT(/1X,10(1H*),1X,'ERROR IN DATA INPUT',1X,4A4)
 3070 FORMAT(/20X,'NODES AT WHICH TOTAL DEPTH IS SET',/20X,33(1H-))
 3090 FORMAT(12I6)
 3120 FORMAT(/12X,'BOUNDARY NODES WITH NORMAL VELOCITY SET (NONZERO)',
     &/12X,49(1H-))
 3130 FORMAT(/13X,'BOUNDARY NODES WITH NORMAL VELOCITY SET TO ZERO',
     &/13X,47(1H-))
 3140 FORMAT(/12X,'PHYSICAL CORNER NODES ... VELOCITY IS SET TO ZERO',
     &/12X,49(1H-))
 3150 FORMAT(/14X,'PHYSICAL CORNER NODES ... VELOCITY IS NONZERO',
     &/14X,45(1H-))
 3160 FORMAT(/22X,'GEOSTROPHIC BOUNDARY ELEMENTS'/22X,29(1H-))
 3170 FORMAT(30X,2I6)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ECHOSETS3(G,OMEGA2,RAD,DEG,COR,NLBS,
     &   CLOSURE,THETA,TAU0,EPSN,EKMIN,AKMIN,NEV)
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
      REAL    G, OMEGA2, RAD, DEG, COR
      INTEGER NLBS
      REAL    THETA, TAU0, EPSN, EKMIN, AKMIN
      INTEGER NEV

      CHARACTER*72 CLOSURE
C
C
      INTEGER INDX
C
      WRITE(2,600)'Physical constants:'
      WRITE(2,601)'Gravitational Acceleration         ',G
      WRITE(2,601)'2*OMEGA for Celestial Body         ',OMEGA2
      WRITE(2,601)'Radius of Celestial Body           ',RAD
      WRITE(2,601)'Degrees latitude                   ',DEG
      WRITE(2,601)'Coriolis parameter                 ',COR
C
C
C
      WRITE(2,600)'Numerical weighting parameters:'
      WRITE(2,601)'THETA                             ',THETA
      WRITE(2,601)'TAU0                              ',TAU0
      WRITE(2,601)'EPSN                              ',EPSN
C
C
C
      WRITE(2,600)'Application of nonlinearities:'
      IF(NLBS.EQ.0)THEN
      WRITE(2,604)'Linear bottom stress coefficient (AK) applied with:'
      WRITE(2,601)'     All AKs set to                ',AKMIN
      ELSE IF(NLBS.EQ.1)THEN
      WRITE(2,604)'get ak(i) from LINSTRESS ' 
      WRITE(2,601)'     Minimum value of AK           ',AKMIN
      ELSE
	 write(2,*)'illegal NLBS ', NLBS
	 stop
      ENDIF
C
C
C
      INDX=INDEX(CLOSURE,' ')-1
      WRITE(2,600)CLOSURE(1:INDX)//' turbulence closure:'
      IF(CLOSURE.EQ.'SPATIAL')THEN
         WRITE(2,604)' get ENZ(I,J) from LINENZ '
      ELSE IF(CLOSURE.EQ.'CONSTANT')THEN
         WRITE(2,601)'All vertical diffusivities set to ',EKMIN
      ELSE 
	 write(2,*)'illegal CLOSURE ', CLOSURE
	 stop
      ENDIF
C
C
C
 600  FORMAT(/A)
 601  FORMAT(5X,A35,1X,E12.4)
 602  FORMAT(5X,A35,1X,I12)
 603  FORMAT(5X,A35,1X,A)
 604  FORMAT(5X,A)
      RETURN
      END
