C***********************************************************************
C***********************************************************************
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C END OF MAIN ROUTINE.  EXTERNALLY FIXED SUBROUTINES FOLLOW
C***********************************************************************
C***********************************************************************
      SUBROUTINE INPUTF6(FILINQ,
     &DEG,NNV,
     &MESHNAME,NN,NE,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &INBE,NBE,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST,NLEV,ZL)
C-----------------------------------------------------------------------
C Objectives: Read .inq file
C             Open .echo file
C             Read mesh information from files indicated in 
C               the .inq file
C             Set horizontal mesh boundary condition arrays
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               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C               NLEV - number of vertical level surfaces for
C                      baroclinic pressure gradient calculations
C               NBE - number of boundary elements
C               NGEO - number of geostrophic boundary elements
C               NE - number of horizontal elements
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
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         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,K1) - LABLIB STANDARD boundary element incidence list
C         NBE - number of boundary elements
C         INGEO(2,K2) - 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,6)
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       Vertical Level Mesh Information:
C
C         NLEV - number of level surfaces for baroclinic pressure
C                calculation
C         ZL(JL) - vertical position of level surfaces for baroclinic
C              pressure calculation
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via INCLUDE file
C
      INCLUDE 'FUNDY6.DIM'
C
C Global dimensioning and assignments
C
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM),
     &NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM),
     &DSBDRY(NNDIM),ZL(NLEVDIM)
      CHARACTER*72 FILINQ,MESHNAME
C
C Local dimensioning and assignments
C
      CHARACTER*72 INQLABEL,COMINQ,FILBEL,FILICQ,FILECH,UNITS
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
      READ(1,'(A)')INQLABEL
      READ(1,'(A)')FILECH
      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,*)NNV
      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 INPUTF5'
         WRITE(2,*)'INVALID UNITS, MODEL ONLY ALLOWS SI UNITS'
         STOP
      ENDIF
C
C Load geometry and boundary condition data
C
      CALL INPUTDECF6(MESHNAME,FILBEL,
     &                     XFCTR,YFCTR,HFCTR,HMIN,
     &                     NN,NE,NNV,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &                     NLEV,ZL,INBE,NBE,
     &                     INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST)
C
C Compute half bandwidth from the element incidence list and
C  perform check to ensure dimensioning parameter NBWE is of sufficient
C  size
C
      CALL HALFBW(IN,3,NE,NHBW)
      NBW=2*NHBW+1
      IF(NBW.GT.NBWE)THEN
        WRITE(2,*)'EXECUTION TERMINATED DUE TO INSUFFICIENT DIMENSIONS:'
        WRITE(2,*)' NBW   =',NBW,   ',NBWE   =',NBWE
        STOP
      END IF
C
C Echo input information to echo file
C
      CALL ECHOINPUTF6(COMINQ,MESHNAME,FILBEL,FILINQ,FILICQ,FILECH,
     &          XFCTR,YFCTR,HFCTR,DEG,HMIN,DELT,X,Y,HDOWN,
     &          NNV,NLEV,ZL,
     &          NN,NE,NBE,NHBW,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE INPUTDECF6(MESHNAME,FILBEL,
     &                     XFCTR,YFCTR,HFCTR,HMIN,
     &                     NN,NE,NNV,X,Y,IN,HDOWN,SUMX,SUMY,DSBDRY,
     &                     NLEV,ZL,INBE,NBE,
     &                     INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST)
C-----------------------------------------------------------------------
C Objectives: Read mesh information from files indicated in .inq file
C             Set horizontal boundary condition arrays
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               NN - number of horizontal nodes
C               NNV - number of vertical nodes
C               NLEV - number of vertical level surfaces for
C                      baroclinic pressure gradient calculations
C               NBE - number of boundary elements
C               NGEO - number of geostrophic boundary elements
C               NE - number of horizontal elements
C
C Inputs: MESHNAME -  name of the horizontal mesh => prefix for names of
C                     the input LABLIB STANDARD .nod, .ele, .bat, 
C                     and .lev files
C         FILBEL - name of the input LABLIB STANDARD .bel file
C         XFCTR,YFCTR - scale factors which convert (X,Y) locations in
C                       the meshname.nod file to MKS units (i.e. meters)
C         HFCTR - scale factor which converts bathymetries listed in
C                 meshname.bat file to MKS units (i.e. meters)
C         HMIN - minimum depth in MKS units (i.e. meters)
C         NNV -  number of vertical nodes
C
C Outputs:
C
C      Horizontal Mesh Information:
C
C         NN - number of horizontal nodes
C         NE - the number of horizontal 2-D triangular elements
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,K1) - LABLIB STANDARD boundary element incidence list
C         NBE - number of boundary elements
C         INGEO(2,K2) - 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,6)
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       Vertical Level Mesh Information:
C
C         NLEV - number of level surfaces for baroclinic pressure
C                calculation
C         ZL(JL) - vertical position of level surfaces for baroclinic
C              pressure calculation
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via an INCLUDE file
C
      INCLUDE 'FUNDY6.DIM'
C
C Global dimensioning and assignments
C
      INTEGER IN(3,NEDIM),INBE(4,NBEDIM),INGEO(2,NBEDIM),
     &NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),SUMX(NNDIM),SUMY(NNDIM),
     &DSBDRY(NNDIM),ZL(NLEVDIM)
      CHARACTER*72 MESHNAME,FILBEL
C
C Local dimensioning and assignments
C
      CHARACTER*72 FILBAT,FILELE,FILNOD,FILLEV
      CHARACTER*72 CLINES
      INTEGER NCD(NNDIM)
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'
      FILLEV=MESHNAME(1:IMESH)//'.lev'
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=4,FILE=FILELE,STATUS='OLD')
      NE=0
 20   READ(4,*,END=21)J,(IN(I,J),I=1,3)
      NE=NE+1
      GO TO 20
 21   CLOSE(4)
C
C Read the level surface list in MKS units (the ".lev" file)
C
      OPEN(UNIT=3,FILE=FILLEV,STATUS='OLD')
      READ(3,'(A)')CLINES
      READ(3,'(A)')CLINES
      NLEV=0
 15   READ(3,*,END=16)J,ZL(J)
      NLEV=NLEV+1
      GO TO 15
 16   CLOSE(3)
C
C Check dimensioning of NN, NE, NNV, and NLEV
C
      IF(NNDIM.GE.NN.AND.NEDIM.GE.NE
     &.AND.NNVDIM.GE.NNV.AND.NLEVDIM.GE.NLEV)GO TO 25
      WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECF5 DUE TO IN
     &SUFFICIENT DIMENSIONS:'
      WRITE(2,*)' NN   =',NN,   ',NNDIM   =',NNDIM
      WRITE(2,*)' NE   =',NE,   ',NEDIM   =',NEDIM
      WRITE(2,*)' NNV  =', NNV, ',NNVDIM  =',NNVDIM
      WRITE(2,*)' NLEV =',NLEV, ',NLEVDIM =',NLEVDIM
      STOP
   25 CONTINUE
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=7,FILE=FILBAT,STATUS='OLD')
      JBAT=0
   30 READ(7,*,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 (7)
      IF(JBAT.NE.NN)THEN
      WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECF5 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=8,FILE=FILBEL,STATUS='OLD')
      NBE=0
      NGEO=0
C
C  Ensure mesh indicated in .bel file matches that from .inq file
C
      READ(8,'(A)')CLINES
      CALL CROPSTRING(MESHNAME,IS1,IB1,IE1)
      CALL CROPSTRING(CLINES,IS2,IB2,IE2)
      IF(MESHNAME(IB1:IE1).NE.CLINES(IB2:IE2))THEN
      WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECF5 DUE TO IN
     &COMPATABILITY BETWEEN THE MESH NAME IN THE .inq AND .bel FILES:'
         WRITE(2,*)' .inq FILE INDICATES '//MESHNAME(IB1:IE1)//' MESH'
         WRITE(2,*)' .bel FILE INDICATES '//CLINES(IB2:IE2)//' MESH'
         STOP
      ENDIF
C
C  Read and ignore header line
C
      READ(8,'(A)')CLINES
C
C  Read the boundary element information
C
 40   READ(8,*,END=41)J,INBE(1,J),INBE(2,J),INBE(3,J),INBE(4,J)
      NBE=NBE+1
      GO TO 40
 41   CLOSE(8)
      IF(NBE.GT.NBEDIM)THEN
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE INPUTDECF5 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 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,8)
C                I=KST,NNDIM    are V.n=0.0         (NCD=1,2,6,8)
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=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 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
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE REORDERINBE(nndim,nbedim,x,y,inbe,nbe)
C-----------------------------------------------------------------------
C purpose: This subroutine determines the number of closed boundaries in
C           a 2-D triangular mesh
C
C restrictions:  Applicable only for Lablib Standard 2-D triangular mesh
C
C inputs:  nndim - dimension parameter for the number of nodes
C          nedim - dimension parameter for the number of elements
C          nbe - the number of boundary elements in the mesh
C          inbe - Lablib Standard incidence list for the boundary 
C                  elements
C
C outputs: The number of closed boundaries found is written to standard
C           out.  If an unclosed boundary is found, a message is written
C           to standard out.
C
C history: Written by Christopher E. Naimie
C          Dartmouth College
C          03 December 1998
C-----------------------------------------------------------------------
C
C ARGUMENTS
C
      real x(nndim),y(nndim)
      integer inbe(4,nbedim)
C
C START OF EXECUTABLE CODE
C
      nopen=0
      nclo=0
      ncflag=0
      istop=0
C
C Normalize all boundary elements such that the material to the left
C of an observer standing on inbe(1,i), looking toward inbe(2,i) is 0
C
      do i=1,nbe
         if(inbe(3,i).EQ.0)then
            if(inbe(4,i).EQ.0)then
               write(2,*)'INVALID BOUNDARY ELEMENT:',i
               write(2,*)'  (INBE(J,I),J=1,4)=',(INBE(J,I),J=1,4)
               istop=1
            endif
         elseif(inbe(4,i).EQ.0)then
            if(inbe(3,i).EQ.0)then
               write(2,*)'INVALID BOUNDARY ELEMENT:',i
               write(2,*)'  (INBE(J,I),J=1,4)=',(INBE(J,I),J=1,4)
               istop=1
            endif
            inbe1i=inbe(1,i)
            inbe(1,i)=inbe(2,i)
            inbe(2,i)=inbe1i
            inbe(3,i)=0
            inbe(4,i)=inbe(3,i)
         else
            write(2,*)'INVALID BOUNDARY ELEMENT:',i
            write(2,*)'  (INBE(J,I),J=1,4)=',(INBE(J,I),J=1,4)
            istop=1
         endif
      enddo
C
C Ensure that each boundary node occurs once and only once in the first
C two rows of inbe
C
      do i=1,nbe
         node=inbe(1,i)
         if(inbe(2,i).EQ.node)then
            write(2,*)'INVALID BOUNDARY ELEMENT:',i
            write(2,*)'  (INBE(J,I),J=1,4)=',(INBE(J,I),J=1,4)
            istop=1
         endif
         ncount=0
         do j=1,i-1
            if(inbe(2,j).EQ.node)ncount=ncount+1
         enddo
         do j=i+1,nbe
            if(inbe(2,j).EQ.node)ncount=ncount+1
         enddo
         if(ncount.NE.1)then
            write(2,*)'INVALID BOUNDARY ELEMENT LIST - NODE: ',node
            write(2,*)' APPEARS ',ncount,' TIMES IN 2ND ROW OF INBE'
            istop=1
         endif
         ncount=1
         do j=1,i-1
            if(inbe(1,j).EQ.node)ncount=ncount+1
         enddo
         do j=i+1,nbe
            if(inbe(1,j).EQ.node)ncount=ncount+1
         enddo
         if(ncount.NE.1)then
            write(2,*)'INVALID BOUNDARY ELEMENT LIST - NODE: ',node
            write(2,*)' APPEARS ',ncount,' TIMES IN 1ST ROW OF INBE'
            istop=1
         endif
      enddo
C
C Find the entry in the boundary element list with the minimum x
C (this must be part of the outer boundary)
C
      nmin=inbe(1,1)
      do i=2,nbe
         if(x(inbe(1,i)).LT.x(nmin))nmin=inbe(1,i)
      enddo
C
C Find the boundary element which has its head at nmin
C
      do i=1,nbe
         if(inbe(2,i).EQ.nmin)then
            nstart=inbe(1,i)
            nhead=inbe(2,i)
            inbe(1,i)=inbe(1,1)
            inbe(2,i)=inbe(2,1)
            inbe(1,1)=nstart
            inbe(2,1)=nhead
            go to 10
         endif
      enddo
C
C Start with first boundary element in list and close boundaries
C
 10   do 20 i=2,nbe
         if(ncflag.EQ.1)then
            ncflag=0
            go to 20
         endif
         do 18 j=i,nbe
            if(inbe(1,j).EQ.nhead)then
               n1=inbe(1,i)
               n2=inbe(2,i)
               inbe(1,i)=inbe(1,j)
               inbe(2,i)=inbe(2,j)
               inbe(1,j)=n1
               inbe(2,j)=n2
               nhead=inbe(2,i)
               go to 19
            endif
 18      continue
         if(i.LT.(nbe))then
            nopen=nopen+1
            nstart=inbe(1,i)
            nhead=inbe(2,i)
            istop=1
         endif
         go to 20           
 19      if(inbe(2,i).EQ.nstart.AND.i.LT.nbe)then
            nclo=nclo+1
            ncflag=1
            nstart=inbe(1,i+1)
            nhead=inbe(2,i+1)
         endif
         if(inbe(2,i).EQ.nstart.AND.i.EQ.nbe)then
            nclo=nclo+1
            go to 21
         endif
 20   continue
      nopen=nopen+1
      istop=1
 21   continue
C
C Report statistical info and stop if errors were found
C
      if(istop.EQ.1)then
         WRITE(2,*)'EXECUTION TERMINATED IN SUBROUTINE REORDERINBE DUE',
     &             ' TO FATAL FLAW IN INPUT BOUNDARY ELEMENT FILE:'
         WRITE(2,*)NCLO,' CLOSED BOUNDARIES WERE FOUND.'
         WRITE(2,*)NOPEN,' OPEN BOUNDARIES WERE FOUND.'
         STOP
      endif
C
C END OF ROUTINE
C
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE ECHOINPUTF6(COMINQ,MESHNAME,FILBEL,
     &          FILINQ,FILICQ,FILECH,
     &          XFCTR,YFCTR,HFCTR,DEG,HMIN,DELT,X,Y,HDOWN,
     &          NNV,NLEV,ZL,
     &          NN,NE,NBE,NHBW,INGEO,NGEO,NORM,NVN,KST,LRVS,NR,NPCN,LST)
C-----------------------------------------------------------------------
C purpose: This subroutine opens the .echo file and then writes 
C          information from the .inq file and information computed
C          using data from the horizontal mesh files.
C-----------------------------------------------------------------------
C
C Declare all the integer program parameters via an INCLUDE file
C
      INCLUDE 'FUNDY6.DIM'
C
C Global dimensioning and assignments
C
      INTEGER INGEO(2,NBEDIM),NORM(NNDIM),LRVS(NNDIM)
      REAL X(NNDIM),Y(NNDIM),HDOWN(NNDIM),ZL(NLEVDIM)
      CHARACTER*72 COMINQ,MESHNAME,FILBEL,FILINQ,FILICQ,FILECH
C
C Declaration of local variables
C
      CHARACTER*15 CMODEL,FILDIM
C
C Define model title, driver program, and user provided subroutine
C
      CMODEL=' T R U X T O N 6 '
      FILDIM='FUNDY6.DIM'
C
C Write output title to the echo file CASENAME.ech (unit 2)
C
      WRITE(2,601)'************************************'//
     &                 '************************************'
      WRITE(2,601)'************************************'//
     &                 '************************************'
      WRITE(2,602)'**', '**'
      WRITE(2,603)'**', '  '//CMODEL, '**'
      WRITE(2,602)'**', '**'
      WRITE(2,604)'**', ' 3-D SHALLOW WATER HYDRODYNAMICS', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,604)'**', '  LINEARIZED HARMONIC EQUATIONS ', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,604)'**', ' WITH LINEAR TRIANGULAR ELEMENTS', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,606)'**', '  RELEASE 1.0; DECEMBER 6, 1998 ', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,606)'**', '  Point of contact (e-mail):  ', '**'
      WRITE(2,606)'**', '  chris.naimie@dartmouth.edu  ', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,601)'************************************'//
     &                 '************************************'
      WRITE(2,601)'************************************'//
     &                 '************************************'
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(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 required 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)'NLEVDIM', NLEVDIM
      WRITE(2,609)'NFTRDIM', NFTRDIM
      WRITE(2,609)'NBEDIM ', NBEDIM
      WRITE(2,609)'NBWE   ', NBWE
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 of horizontal mesh  ', 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 in meters (MIN,MAX)     ',XMIN,XMAX
      WRITE(2,616)'Range of Y in meters (MIN,MAX)     ',YMIN,YMAX
      WRITE(2,616)'Range of depth in meters (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
      WRITE(2,611)'Number of level surfs. for BPG calc',NLEV
      WRITE(2,610)'Level surfs. in meters (top to bot)',ZL(NLEV)
      DO JL=NLEV-1,1,-1
      WRITE(2,610)'                                   ',ZL(JL)
      ENDDO
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 QUODDY:'
      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 ECHOSETF6(G,OMEGA,DEG,COR,FREQ)
C-----------------------------------------------------------------------
C purpose: This subroutine writes the settings for various run 
C          parameters and switch settings to the .echo file.
C-----------------------------------------------------------------------
C
C
C
      WRITE(2,600)'Physical parameters:'
      WRITE(2,601)'Gravitational acceleration (m/s^2):',G
      WRITE(2,601)'Celestial body rot. rate   (rad/s):',OMEGA
      WRITE(2,601)'Latitude                     (deg):',DEG
      WRITE(2,601)'Coriolis parameter         (rad/s):',COR      
      WRITE(2,601)'Tidal frequency            (rad/s):',FREQ
C
C
C
 600  FORMAT(/A)
 601  FORMAT(5X,A35,1X,E12.4)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE HALFBW(IN,NLOC,NE,NB)
C-----------------------------------------------------------------------
C  NB IS THE HALF-BANDWIDTH
C  NE IS THE NUMBER OF ELEMENTS
C  NOC IS THE FIRST DIMENSION OF IN
C  IN CONVENTION:IN(3,NEDIM)
C-----------------------------------------------------------------------
      DIMENSION IN(NLOC,*)      
      NB=0       
      DO L=1,NE       
         MI=IN(1,L)               
         MA=0       
         DO I=1,3        
            J=IN(I,L)           
            IF(J.LT.MI) MI=J             
            IF(J.GT.MA) MA=J             
         ENDDO 
         IDIFF=MA-MI         
         IF(IDIFF.GT.NB) NB=IDIFF     
      ENDDO 
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SUBTITLE (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***********************************************************************
      SUBROUTINE CROPSTRING(string,istart,ibeg,iend)
C-----------------------------------------------------------------------
C purpose: This subroutine finds 3 pointers for for a character string
C
C inputs:  string - a character string of length .LE. 72
C
C outputs: istart - first nonblank character in string
C          ibeg   - location of character in string which follows the 
C                   last "/" entry before iend
C          iend   - first blank character in string after istart
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           08 DEC 1998
C-----------------------------------------------------------------------
C
C ARGUMENTS
      character*72 string
C
C START OF EXECUTABLE CODE
C
C find beginning of string
      do 10 i=1,72
         istart=i
	 if(string(i:i).NE.' ')go to 11
 10   continue
C
C find end of string
 11   do 20 i=istart,72
         iend=i-1
	 if(string(i:i).EQ.' ')go to 21
 20   continue
C 
C find the last / before end of string
 21   do 30 i=iend,istart,-1
         ibeg=i+1
         if(string(i:i).EQ.'/')go to 31
 30   continue
      ibeg=istart
C
C END OF ROUTINE
 31   return
      end
C***********************************************************************
C***********************************************************************
