C***********************************************************************
C***********************************************************************
C Beginning of Fixed Subroutines for QUODDY4 (RELEASE 2.1)
C Beginning of Fixed Subroutines for QUODDY4 (RELEASE 2.1)
C***********************************************************************
C***********************************************************************
C
C Equations of state:
C
C***********************************************************************
C***********************************************************************
	SUBROUTINE EQSTATE2(T,S,R)
C-----------------------------------------------------------------------
C purpose: This subroutine computes the perturbation density (R; in 
C          sigmat units) as a function of two variables; temperature 
C          (T; 0 to 40 degrees C) and salinity (S; 0 to 42 psu).
C
C reference: UNESCO equation of state as published in Gill(1982)
C
C Contributed by Douglas T. Morgan
C-----------------------------------------------------------------------
C
C...Ensure T and S are within Range
C
      IF((T.GE.0).AND.(T.LE.40))GO TO 10
      WRITE(2,'(///A)')'EXECUTATION TERMINATED IN EQSTATE2 AS T LIES OUT
     &SIDE PERMISSIBLE RANGE (0 to 40 degrees C).'
      STOP
 10   IF((S.GE.0).AND.(S.LE.42))GO TO 20
      WRITE(2,'(///A)')'EXECUTATION TERMINATED IN EQSTATE2 AS S LIES OUT
     &SIDE PERMISSIBLE RANGE (0 to 42 psu).'
      STOP      
C
C...Density of pure water
C
 20   R=-0.157406
     &  +6.793952E-2*T-9.095290E-3*T**2.0
     &  +1.001685E-4*T**3.0-1.120083E-6*T**4.0
     &  +6.536332E-9*T**5.0
C
C...Density at one standard atmosphere
C
      R=R+S*(0.824493-4.0899E-3*T+7.6438E-5*T**2.0
     &          -8.2467E-7*T**3.0+5.3875E-9*T**4.0)
     &   +S**1.5*(-5.72466E-3+1.0227E-4*T-1.6546E-6*T**2.0)
     &   +4.8314E-4*S**2.0
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE EQSTATE1(T,R)
C-----------------------------------------------------------------------
C purpose: This subroutine computes the perturbation density (R; in 
C          sigmat units) as a linear function of one variable; 
C          temperature (in degrees C).
C-----------------------------------------------------------------------
      R=1028.0*(1.0-1.7E-4*(T-10.0))-1000.0
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Echo file writing:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE ECHOINPUTQ4(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 'QUODDY.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=' Q U O D D Y 4 '
      FILDIM='QUODDY.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,605)'**', 'WITH IMPLICIT LINEAR TRIANGULAR '//
     &                 'ELEMENTS', '**'
      WRITE(2,602)'**', '**'
      WRITE(2,606)'**', '  RELEASE 2.1; JUNE 20, 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(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)'NLEVDIM', NLEVDIM
      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 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 (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
      WRITE(2,611)'Number of level surfs. for BPG calc.',NLEV
      WRITE(2,610)'Level surf. locations (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 ECHOSETQ4(T0,S0,G,OMEGA,DEG,COR,NHBW,
     &NONLIN,NLBS,Cd,ISMAG,AH,AHMIN,PRESSURE,MASSVAR,IHHBC,CLOSURE,
     &Q20,Q2L0,ZLOGBOT,ZLOGTOP,IQADVDIF,IQ2TBC,IQ2BBC,IQ2LTBC,IQ2LBBC,
     &TAU0,EPSN,EPSH,EPSQ,THETA,EKMIN,AKMIN,Q2MIN,Q2LMIN,ELLMIN,NEV)
C-----------------------------------------------------------------------
C purpose: This subroutine writes the settings for various run 
C          parameters and switch settings to the .echo file.
C-----------------------------------------------------------------------
      CHARACTER*72 PRESSURE,MASSVAR,CLOSURE
      REAL ZERO
      DATA ZERO/0.0/
C
C
C
      WRITE(2,600)'Physical constants:'
      WRITE(2,601)'Gravitational Acceleration         ',G
      WRITE(2,601)'Rotation rate for celestial body   ',OMEGA
      WRITE(2,601)'Degrees latitude                   ',DEG
      WRITE(2,601)'Coriolis parameter                 ',COR
C
C
C
      WRITE(2,600)'Numerical weighting parameters:'
      WRITE(2,601)'TAU0                              ',TAU0
      WRITE(2,601)'EPSN                              ',EPSN
      WRITE(2,601)'EPSH                              ',EPSH
      WRITE(2,601)'EPSQ                              ',EPSQ
      WRITE(2,601)'THETA                             ',THETA
      IF(THETA.EQ.ZERO)THEN
      WRITE(2,604)'     =>Explicit Solution; the minimum'
      ELSE
      WRITE(2,604)'     =>Implicit Solution; the minimum'
      ENDIF
      WRITE(2,602)'      allowable value for NBWE is ',2*NHBW+1
C
C
C
      WRITE(2,600)'Application of nonlinearities:'
      IF(NONLIN.EQ.0)THEN
      WRITE(2,604)'Advection of momentum NOT included'
      ELSE IF(NONLIN.EQ.1)THEN
      WRITE(2,604)'Advection of momentum included'
      ELSE
      ENDIF
      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)'Quadratic bottom stress coefficient (AK) applied with
     &:'
      WRITE(2,601)'     Bottom stress drag coefficient',Cd
      WRITE(2,601)'     Minimum value of AK           ',AKMIN
      ELSE
      ENDIF
      IF(ISMAG.EQ.0)THEN
      WRITE(2,601)'Uniform hor. viscosity/diffusivity ',AH
      ELSE IF(ISMAG.EQ.1)THEN
      WRITE(2,604)'Local shear and grid dependent horizontal viscosity/d
     &iffusivity:'
      WRITE(2,601)'     Multiplicative constant set to   ',AH
      WRITE(2,601)'     Minimum viscosity/diffusivity    ',AHMIN
      ELSE
      ENDIF
C
C
C
      WRITE(2,600)'Application of mass variables:'
      INDX=INDEX(PRESSURE,' ')-1
      WRITE(2,604)PRESSURE(1:INDX)//' pressure calculation'
      INDX=INDEX(MASSVAR,' ')-1
      WRITE(2,604)MASSVAR(1:INDX)//' mass variables:'
      IF(MASSVAR.EQ.'TWO-PROG')THEN
      WRITE(2,603)'     Equation of state subroutine  ','    EQSTATE2'
      WRITE(2,601)'     Reference Temperature         ',T0
      WRITE(2,601)'     Reference Salinity            ',S0
      ELSE IF(MASSVAR.EQ.'TWO-DIAG')THEN
      WRITE(2,603)'     Equation of state subroutine  ','    EQSTATE2'
      WRITE(2,601)'     Reference Temperature         ',T0
      WRITE(2,601)'     Reference Salinity            ',S0
      ELSE IF((MASSVAR.EQ.'ONE-PROG'))THEN
      WRITE(2,603)'     Equation of state subroutine  ','    EQSTATE1'
      WRITE(2,601)'     Reference Temperature         ',T0
      ELSE IF((MASSVAR.EQ.'ONE-DIAG'))THEN
      WRITE(2,603)'     Equation of state subroutine  ','    EQSTATE1'
      WRITE(2,601)'     Reference Temperature         ',T0
      ELSE IF(MASSVAR.EQ.'ZERO')THEN
      WRITE(2,603)'     Equation of state subroutine  ','NONE'
      ELSE
      ENDIF
      WRITE(2,602)'Horizontal BC type at open bdrys   ',IHHBC
C
C
C
      INDX=INDEX(CLOSURE,' ')-1
      WRITE(2,600)CLOSURE(1:INDX)//' turbulence closure:'
      IF(CLOSURE.EQ.'MY25')THEN
      WRITE(2,601)'Reference value for Q2             ',Q20
      WRITE(2,601)'Reference value for Q2L            ',Q2L0
      WRITE(2,601)'Dist. of bot node above true bot   ',ZLOGBOT
      WRITE(2,601)'Dist. of top node below true surf  ',ZLOGTOP
      IF(IQADVDIF.EQ.0)THEN
      WRITE(2,603)'Include adv. and horiz. diffus.?   ','          NO'
      ELSE IF(IQADVDIF.EQ.1)THEN
      WRITE(2,603)'Include adv. and horiz. diffus.?   ','         YES'
      ELSE
      ENDIF
      WRITE(2,602)'Bottom BC type for Q2              ',IQ2BBC
      WRITE(2,602)'Surface BC type for Q2             ',IQ2TBC
      WRITE(2,602)'Bottom BC type for Q2L             ',IQ2LBBC
      WRITE(2,602)'Surface BC type for Q2L            ',IQ2LTBC
      WRITE(2,601)'Minimum vert. viscosity/diffusivity',EKMIN
      WRITE(2,601)'Minimum value for Q2               ',Q2MIN
      WRITE(2,601)'Minimum value for Q2L              ',Q2LMIN
      WRITE(2,601)'Minimum value for mixing length    ',ELLMIN
      ELSE IF(CLOSURE.EQ.'CONSTANT')THEN
      WRITE(2,601)'All vertical diffusivities set to ',EKMIN
      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
C***********************************************************************
C***********************************************************************
C
C Sparse matrix manipulation:
C
C***********************************************************************
C***********************************************************************
      FUNCTION KAY(I,J,IQ,JQ)
C-----------------------------------------------------------------------
C  FINDS ADDRESS K IN QV FOR Q(I,J)
C
C  Notes: This version differs from the SPRSPAK version only in that 
C         (*)'s are used to dimension the agruments [instead of (1)'s].
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION IQ(*),JQ(*)
      KAY=IQ(I)+1
   10 KAY=KAY-1
      IF(JQ(KAY).NE.J)GO TO 10
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN(I,QV,IQ,JQ,X,B,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRIX X,
C     FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTOR B
C     (OR: OBTAINS ONLY Bi OF Q*X=B)
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION QV(*),IQ(*),JQ(*),X(NNDIM,*),B(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         QVK=QV(K)
         JQK=JQ(K)
         DO J=1,NNV
            B(J)=B(J)+QVK*X(JQK,J)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN2(I,QV,IQ,JQ,X1,B1,X2,B2,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV BY THE Jth COLUMN OF THE MATRICIES
C     X1 and X2, FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTORS B1 and B2
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X1 AND X2
C  THIS ROUTINE IS NONDESTRUCTIVE
C
C  This routine is a variation of the SPRSPAK SUBROUTINE SPRSMLTIN.  
C  SPRSMLTIN was converted to SPRSMLTIN2 by Christopher E. Naimie.
C-----------------------------------------------------------------------
      INTEGER IQ(*),JQ(*)
      REAL QV(*),X1(NNDIM,*),B1(*),X2(NNDIM,*),B2(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B1(J)=0.
         B2(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         JQK=JQ(K)
         QVK=QV(K)
         DO J=1,NNV
            B1(J)=B1(J)+QVK*X1(JQK,J)
            B2(J)=B2(J)+QVK*X2(JQK,J)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSMLTIN3(I,IQ,JQ,X,QV1,B1,QV2,B2,QV3,B3,NNV,NNDIM)
C-----------------------------------------------------------------------
C  MULTIPLIES THE Ith ROW OF QV1, QV2, AND QV3 BY THE Jth COLUMN OF THE 
C     MATRIX X,FOR ALL VALUES OF J FROM 1 TO NNV
C  ANSWER RETURNS IN THE VECTORS B1, B2, AND B3
C  IQ,JQ ARE THE STANDARD POINTERS FOR QV
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF X
C  THIS ROUTINE IS NONDESTRUCTIVE
C
C  This routine is a variation of the SPRSPAK SUBROUTINE SPRSMLTIN.  
C  SPRSMLTIN was converted to SPRSMLTIN2 by Christopher E. Naimie.
C-----------------------------------------------------------------------
      INTEGER IQ(*),JQ(*)
      REAL X(NNDIM,*),QV1(*),B1(*),QV2(*),B2(*),QV3(*),B3(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         B1(J)=0.
         B2(J)=0.
         B3(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         JQK=JQ(K)
         QVK1=QV1(K)
         QVK2=QV2(K)
         QVK3=QV3(K)
         DO J=1,NNV
            XJQKJ=X(JQK,J)
            B1(J)=B1(J)+QVK1*XJQKJ
            B2(J)=B2(J)+QVK2*XJQKJ
            B3(J)=B3(J)+QVK3*XJQKJ
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSCONV(I,PPX,PPY,IQ,JQ,U,V,Z,
     &  DUDX,DUDY,DVDX,DVDY,DZDX,DZDY,
     &  NNV,NNDIM)
C-----------------------------------------------------------------------
CCOMPUTES THE HORIZONTAL DERIVATIVES AT DEPTH, UNDER HORIZONTAL NODE #I.
C  PPX, PPY ARE THE SPRSPAK MATRICES FOR D/DX, D/DY
C  IQ,JQ ARE THE STANDARD POINTERS FOR PPX, PPY
C  U(I,J),V(I,J) ARE THE X- AND Y-VELOCITIES 
C      AT HORIZONTAL NODE I, VERTICAL NODE J
C  Z(I,J) IS THE VERTICAL POSITION
C      OF HORIZONTAL NODE I, VERTICAL NODE J
C  THE HORIZONTAL DERIVATIVE TERMS DUDX(J), ETC. ARE RETURNED 
C      FOR ALL NODES J UNDER THE GIVEN NODE I
C  NNV IS THE ACTUAL # OF NODES IN THE VERTICAL
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF U AND V 
C  THIS ROUTINE IS NONDESTRUCTIVE
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION PPX(*),PPY(*),IQ(*),JQ(*),
     &U(NNDIM,*),V(NNDIM,*),Z(NNDIM,*),
     &DUDX(*),DUDY(*),DVDX(*),DVDY(*),DZDX(*),DZDY(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         DUDX(J)=0.
         DUDY(J)=0.
         DVDX(J)=0.
         DVDY(J)=0.
         DZDX(J)=0.
         DZDY(J)=0.
      ENDDO
      DO K=KMIN,KMAX
         PPXK=PPX(K)
         PPYK=PPY(K)
         DO J=1,NNV
            UJQKJ=U(JQ(K),J)
            VJQKJ=V(JQ(K),J)
            ZJQKJ=Z(JQ(K),J)
            DUDX(J)=DUDX(J)+PPXK*UJQKJ
            DUDY(J)=DUDY(J)+PPYK*UJQKJ
            DVDX(J)=DVDX(J)+PPXK*VJQKJ
            DVDY(J)=DVDY(J)+PPYK*VJQKJ
            DZDX(J)=DZDX(J)+PPXK*ZJQKJ
            DZDY(J)=DZDY(J)+PPYK*ZJQKJ
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSGRAD(I,PPX,PPY,IQ,JQ,F,Z,
     &  DFDX,DFDY,DZDX,DZDY,
     &  NNV,NNDIM)
C-----------------------------------------------------------------------
C  A SIMPLE VARIANT OF SPRSCONV
C  COMPUTES THE HORIZONTAL GRADIENT OF F AT DEPTH, UNDER HORIZONTAL
C  NODE #I.
C  PPX, PPY ARE THE SPRSPAK MATRICES FOR D/DX, D/DY
C  IQ,JQ ARE THE STANDARD POINTERS FOR PPX, PPY
C  F(I,J) IS A SCALAR FUNCTION DEFINED 
C      AT HORIZONTAL NODE I, VERTICAL NODE J
C  Z(I,J) IS THE VERTICAL POSITION
C      OF HORIZONTAL NODE I, VERTICAL NODE J
C  THE HORIZONTAL DERIVATIVE TERMS DFDX(J), DFDY(J) ARE RETURNED
C      FOR ALL NODES J UNDER THE GIVEN NODE I
C  NNV IS THE ACTUAL # OF NODES IN THE VERTICAL
C  NNDIM IS THE FIRST PHYSICAL DIMENSION OF F
C  THIS ROUTINE IS NONDESTRUCTIVE
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION PPX(*),PPY(*),IQ(*),JQ(*),Z(NNDIM,*),DZDX(*),DZDY(*)
      DIMENSION F(NNDIM,*),DFDX(*),DFDY(*)
      KMIN=1
      IF(I.GT.1) KMIN=IQ(I-1)+1
      KMAX=IQ(I)
      DO J=1,NNV
         DFDX(J)=0.
         DFDY(J)=0.
         DZDX(J)=0.
         DZDY(J)=0.
         DO K=KMIN,KMAX
            DFDX(J)=DFDX(J)+PPX(K)*F(JQ(K),J)
            DFDY(J)=DFDY(J)+PPY(K)*F(JQ(K),J)
            DZDX(J)=DZDX(J)+PPX(K)*Z(JQ(K),J)
            DZDY(J)=DZDY(J)+PPY(K)*Z(JQ(K),J)
         ENDDO
      ENDDO
      RETURN
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE SPRSMLT(QV,IQ,JQ,X,B,NN)
C-----------------------------------------------------------------------
C  MULTIPLIES QV*X; ANSWER RETURNS IN B
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION QV(*),X(*),B(*),IQ(*),JQ(*)
      KMAX=0
      DO I=1,NN
         SUM=0.
         KMIN=KMAX+1
         KMAX=IQ(I)
         DO K=KMIN,KMAX
            SUM=SUM+QV(K)*X(JQ(K))
         ENDDO
         B(I)=SUM
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSBLD1(IN,NDIM,NNE,IQ,JQ,NN,NE,MN,NDF)
C-----------------------------------------------------------------------
C  NOVEMBER 30, 1990
C  D.LYNCH REVISION OF K.PAULSEN REVISION OF WEIPING VERSION.
C  SIMPLE ADAPTATION OF SPRSBLDW1 BY REVERSING ORDER OF IN(,) ARGUMENTS.
C
C  BUILDS SPRSPAK ARRAYS IQ,JQ FROM ELEMENT INCIDENCE LIST IN(K,L)
C     WITH L=GLOBAL ELEMENT NUMBER, K=LOCAL NODE NUMBER
C     FOR NDF DEGREES OF FREEDOM AT EACH NODE -- i.e. 
C     NDF EQUATIONS AND NDF UNKNOWNS.
C     NDF IS THE SAME AT EACH NODE.
C
C  SAME AS SPRSBLDW1 EXCEPT FOR ORDER OF ARGUMENTS FOR THE IN ARRAY
C 
C  NDIM IS THE FIRST PHYSICAL DIMENSION OF IN; NDIM.GE.NNE
C  NNE INDICATES THAT IN(I,L) FOR I=1 THRU NNE ARE NODE NUMBERS; 
C     THE BALANCE (IF ANY) OF IN WILL BE IGNORED.  
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C  JQ(K) IS THE COLUMN INDEX IN Q FOR QV(K)
C  NN IS THE ACTUAL NUMBER OF NODES
C  NE IS THE ACTUAL NUMBER OF ELEMENTS
C  MN EQUALS OR EXCEEDS THE MAXIMUM NUMBER OF NONZREO ENTRYS PER ROW, 
C     FOR A SCALAR (i.e. NDF=1) PROBLEM.
C  NDF IS THE DEGREES OF FREEDOM, ASSUMED THE SAME AT EACH NODE.
C
C  IT IS ASSUMED THAT IN MAIN PROGRAM, DIMENSIONS ARE AT LEAST 
C     IQ(NN*NDF), JQ((NN*NDF)*(MN*NDF)).  NOTE THAT JQ IS OVERSIZED; 
C     BUT THE ALGORITHM MAY USE THIS EXTRA SIZE.  
C
C  CALLS ARE MADE TO SPRSINSRTD AND SPRSBLD3DD
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION IN(NDIM,*),IQ(*),JQ(*)
C
C  INITIALIZE JQ
C
      DO I=1,NN*MN
         JQ(I)=0
      ENDDO
C
C  INITIALIZE JQ AS DIAGONAL MATRIX
C
      DO I=1,NN
         J=(I-1)*MN+1
         JQ(J)=I
      ENDDO
C
C  ELEMENT LOOP FOR JQ, ASSUMING NDF=1.
C  ALL ROWS ALLOWED MN ENTRIES.
C  TRAILING ZEROS WILL BE REMOVED BELOW.
C
      DO L=1,NE
         DO I=1,NNE
            II=IN(I,L)
            DO J=1,NNE
               JJ=IN(J,L)
               IF(I.NE.J)CALL SPRSINSRTD(II,JJ,JQ,MN)
            ENDDO
         ENDDO
      ENDDO
C
C FIND IQ FROM JQ
C
       K=0
       KK=0
       DO I=1,NN
          DO J=1,MN
             KK=KK+1
             IF(JQ(KK).NE.0)K=K+1
          ENDDO
          IQ(I)=K
       ENDDO      
C
C REARRANGE (CONDENSE) THE JQ -- REMOVE THE TRAILING ZEROS
C
      DO I=2,NN
         KMIN=IQ(I-1)+1
         KMAX=IQ(I)
         KJMIN=(I-1)*MN+1
         DO K=KMIN,KMAX
            JQ(K)=JQ(KJMIN)
            KJMIN=KJMIN+1
         ENDDO
      ENDDO
C
C IQ, JQ NOW CORRECT FOR NDF=1
C EXPAND FOR NDF.GT.1
C
      IF(NDF.GT.1)CALL SPRSBLD3DD(IQ,JQ,NN,NDF)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSBLD3DD(IQ,JQ,NN,NDF)
C-----------------------------------------------------------------------
C BUILD 3D IDEX IQQ AND JQQ
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION IQ(*),JQ(*)
      DO I=NN,1,-1
         IF (I.NE.1) THEN
            II=IQ(I)-IQ(I-1)
         ELSE
            II=IQ(I)
         ENDIF
         DO J=NDF,1,-1
            IQ(NDF*(I-1)+J)=IQ(I)*NDF*NDF-(NDF-J)*NDF*II
         ENDDO
         DO K=1,II
            JK=IQ(NDF*I)+1-(K-1)*NDF
            IF (I.NE.1) THEN
               JJ=JQ(IQ(I)-K+1)*NDF+1
            ELSE
               JJ=JQ(IQ(I)/NDF-K+1)*NDF+1
            END IF
            DO J=1,NDF
               JK=JK-1
               JJ=JJ-1
C
C THE FOLLOWING 3 LINES FROM WEIPING'S CODE
C	   JQ(JK)=JJ
C	   JQ(JK-NDF*II)=JJ
C	   IF(NDF.EQ.3)	   JQ(JK-6*II)=JJ
C WORKED FOR NDF=2 OR 3.  REPLACED BY INDUCTION  
C WITH THE DO 15 LOOP BELOW (DRL 1 DECEMBER 1990).   
C IT IS NOT OBVIOUS THAT JQ REMAINS UNCORRUPTED IN THE 
C OVERWRITING.  
C
               DO L=1,NDF
                  JQ(JK-(L-1)*NDF*II)=JJ
               ENDDO
            ENDDO
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSINSRTD(I,J,JQ,MN)
C-----------------------------------------------------------------------
C  INSERTS GLOBAL ADDRESS (I,J) INTO JQ(), 
C  ASSUMING EXACTLY MN ENTRIES IN EACH ROW.
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION JQ(*)
      KMIN=(I-1)*MN+1
      KMAX=I*MN
C
C  SEARCH FOR INSERTION POINT KK
C
      DO K=KMIN,KMAX
         KK=K
         IF(J.LT.JQ(K))GO TO 20
         IF(J.EQ.JQ(K))RETURN
         IF(JQ(K).EQ.0)GO TO 20
      ENDDO
C
C  INSERT J AT JQ(KK); BUMP JQ
C
C CHECK INSERTED BY K PAULSEN 11/29...NOT VALIDATED AT TIME OF
C INSERTION !!!!
C
   20 IF(JQ(KMAX).NE.0)THEN
         WRITE(*,*) 'MAX CONNECTIVITY MN=',MN,
     &   'VIOLATED AT I,J=',I,J,'IN SPRSINSRTD'
         STOP
      ENDIF
      DO K=KMAX-1,KK,-1
         JQ(K+1)=JQ(K)
      ENDDO
      JQ(KK)=J   
C
C  DONE
C
      RETURN     
      END
C*********************************************************************
C*********************************************************************
      SUBROUTINE SPRSINVMLT(QV,IQ,D,NN)
C-----------------------------------------------------------------------
C  REPLACES QV WITH INV(D)*QV
C  D IS A DIAGONAL MATRIX, LENGTH NN 
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C
C  Notes: This version differs from the SPRSPAK version in 2 ways:
C           1.(*)'s are used to dimension agruments [instead of (1)'s],
C           2.The do loops have been modified to a DO-ENDDO style.
C         Modified by Christopher E. Naimie
C-----------------------------------------------------------------------
      DIMENSION QV(*),D(*),IQ(*)
      KMAX=0
      DO I=1,NN
         KMIN=KMAX+1
         KMAX=IQ(I)
         FTR=1./D(I)
         DO K=KMIN,KMAX
            QV(K)=FTR*QV(K)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C...Horizontal viscosity/diffusivity:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE SMAGOR(AGPGP,U,V,AR,X,Y,IN,IQ,JQ,NE,SV,NN,AHI,
     &                  SMAGR,AHMIN,NFTR)
C-----------------------------------------------------------------------
C  Purpose:  This subroutine determines the node-based nonlinear 
C            horizontal diffusion matrix (AGPGP) and node-based values
C            of horizontal viscosity/diffusivity (AHI).  The formulation
C            is similar in spirit to Smagorinsky (1963) (i.e. horizontal 
C            viscosity/diffusivity is parameterized in terms of the 
C            local shear and grid scale).  The exact Laplacian form 
C            implemented herein is as presented in Lynch et al (1996) 
C            [CSR, pp 875-906].
C
C  Notes: This version differs from the earlier version in 2 ways:
C           1. A minimum horizontal viscosity/diffusivity (AHMIN) is 
C              implemented.
C           2. The do loops have been modified to a DO-ENDDO style.
C         As in Lynch et al (1996), the coefficient specified by
C              the user multiplies the square of the grid scale and
C              the local shear.
C         Original Version: Daniel R. Lynch and Justin T.C. Ip, 3 AUG 93
C         Modified by Christopher E. Naimie, 19 JUN 98
C-----------------------------------------------------------------------
C
C Global arrays
C
      REAL AGPGP(*),U(*),V(*),AR(*),X(*),Y(*),AHI(*),SV(*)
      INTEGER IN(3,*),IQ(*),JQ(*)
C
C Local dimensioning and assignments
C
      REAL DELY(3),DELX(3)
C
C Initialize AGPGP and AHI
C
      DO I=1,NFTR*NN
         AGPGP(I)=0.0
      ENDDO
      DO I=1,NN
         AHI(I)=0.0
      ENDDO
C
C  Laplacian matrix:
C     AGPGP = -<AH*grad(phi_i) . grad(phi_j)> 
C     AGPGP*V = + <DEL.(AH GRAD V), phi_i>
C
      DO L=1,NE       
         ARR=AR(L) 
         I1=IN(1,L)
         I2=IN(2,L)
         I3=IN(3,L)
         DELX(1)=X(I2)-X(I3)          
         DELX(2)=X(I3)-X(I1)          
         DELX(3)=X(I1)-X(I2)          
         DELY(1)=Y(I2)-Y(I3)          
         DELY(2)=Y(I3)-Y(I1)          
         DELY(3)=Y(I1)-Y(I2)          
         D1=(U(I1)*DELY(1)+U(I2)*DELY(2)+U(I3)*DELY(3)   
     &      +V(I1)*DELX(1)+V(I2)*DELX(2)+V(I3)*DELX(3))/2./ARR    
         D2=(V(I1)*DELY(1)+V(I2)*DELY(2)+V(I3)*DELY(3)   
     &      -U(I1)*DELX(1)-U(I2)*DELX(2)-U(I3)*DELX(3))/2./ARR
C
C  Notes: 2.*ARR is a measure of the local grid scale squared
C         SQRT(D1**2+D2**2) is a measure of the local shear
C
         AH=SMAGR*2.*ARR*SQRT(D1**2+D2**2)    
         DO I=1,3
            II=IN(I,L)
            AHI(II)=AHI(II)+AH*ARR/3.
            DO J=1,3
               IJ=IN(J,L)
               K=KAY(II,IJ,IQ,JQ)
               AGPGP(K)=AGPGP(K)-AH*(DELX(J)*DELX(I)+DELY(J)*DELY(I))
     &                  /4./ARR
            ENDDO
         ENDDO
      ENDDO
C
C Divide by mass matrix to compute AHI and compare to minimum value
C
      DO I=1,NN
         AHI(I)=AHI(I)/SV(I)
         AHI(I)=AMAX1(AHI(I),AHMIN)
      ENDDO
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C...Turbulence closure schemes:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE QUADSTRESS(AK,Cd,AKMIN,NN,UZ,VZ,NNDIM)
C-----------------------------------------------------------------------
C This subroutine outputs AK for quadratic bottom stress law
C-----------------------------------------------------------------------
      REAL AK(*),UZ(NNDIM,*),VZ(NNDIM,*)
      DO I=1,NN
         VEB=SQRT(UZ(I,1)**2+VZ(I,1)**2)
         AK(I)=AMAX1(Cd*VEB,AKMIN)
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE GALPERINQ4(UI,VI,Q2I,Q2LI,RHOI,ZI,
     &                G,NNV,NEV,ZLOGBOT,ZLOGTOP,EKMIN,Q2MIN,ELLMIN,
     &                RHSQ2,RHSQ2L,DKQ2,DKQ2L,EKM,EKH,EKQ)
C-----------------------------------------------------------------------
C purpose: This subroutine applies the level 2.5 turbulence closure to
C          determine the vertical mixing coefficients (EK*) along with
C          the source (RHS*) and decay (DK*) terms needed the next time
C          the transport equations for Q2 and Q2L will be solved.
C
C algorithmic design: This subroutine is applied for each node in the
C                     horizontal finite element mesh during each time
C                     step.
C
C index range: J=1,NNV, JL=1,NEV
C
C inputs:  UI(J),VI(J) - x,y components of horizontal velocity
C          Q2I(J),Q2LI(J) - turbulence transport variables
C          RHOI(J) - (density-reference density)/reference density 
C          ZI(J) - vertical grid
C          G - gravitational acceleration
C          NNV - number of vertical nodes
C          NEV - number of vertical elements (NEV=NNV-1)
C          ZLOGBOT - height of the bottom node in the computational 
C                    domain above the true bottom
C          ZLOGTOP - distance from the top node in the computational 
C                    domain to the true surface
C          EKMIN - minimum vertical eddy diffusivity for momentum, heat,
C          Q2MIN,ELLMIN - minimum values of Q2 and the mixing length.
C
C outputs: RHSQ2(J) - RHS source for Q2 transport equation
C          RHSQ2L(J) - RHS sources for Q2L transport equation
C          DKQ2(J) - LHS semi-implicit decay term for Q2 transport equ.
C          DKQ2(J) - LHS semi-implicit decay term for Q2L transport equ.
C          ENM(JL) - vertical eddy mixing diffusivity for momentum
C          EKH(JL) - vertical eddy mixing diffusivity for heat and salt
C          EKQ(JL) - vertical eddy mixing diffusivity for Q2 and Q2L
C
C history:
C
C SUBROUTINE MY25
C   subroutine to return right hand side and turbulent eddy 
C   viscosity terms based on the Mellor-Yamada 2.5 closure scheme
C   (Dan and Cisco, 14 April 1992)
C
C   modified 23 april by DRL 
C     to enforce locally-defined minima for ekm, ekq, ekh
C   modified 24 april by DRL 
C     cd removed from this subroutine; rely on main pgm for that.
C   modified 19 july 93 by DRL and CEN
C     two one-line errors removed (see below)
C   modified 4 august 93 by DRL and JTCI
C     sign of gh changed 
C     maximum values for gh and gm enforced as in mellor/yamada 82
C 
C SUBROUTINE GALPERIN
C   modified 24 september 93 by DRL:
C     stability fctns Sm, Sh redone according to Galperin et al, 
C     (1988).  Maxima for gh, gm removed.
C   modified 1 nov 93 by DRL and JTI:
C     maximum mixing length enforced
C     maximum for GH enforced
C   modified 21 april 1994 by JTI:
C     change limits on q2, l, Kmmin, Khmin, and Kqmin to be 1E-6
C     (the same as in Naimie's routine)
C     inplement wpf according to Blumberg et al., 1992
C
C SUBROUTINE GALPERIN2
C   DRL 2 September 94
C     From subr. galperin in /u/GLOBEC/quoddy2/sourcecode; 
C     Minor changes (noted) and general cleanup.
C
C   21 Oct. 94: added zlogbot to calling arguments; and adjusted WPF
C        at top.
C
C SUBROUTINE GALPERIN3
C   In QUODDY3 version 3, subroutine GALPERIN3 in fixsubs.f is used
C   instead of GALPERIN2 in all the preceding versions.  Subroutine
C   GALPERIN3 is revised from GALPERIN2 with an additional variable,
C   EKMIN representing minimum vertical viscosity and diffusivity
C   inputted from the quoddy3.inq file.  With EKMIN, users are able
C   to set the values for the minimum vertical viscosity and
C   diffusivity in the turbulence closure calculation.  The same
C   number EKMIN will be used as a minimum for viscosity and
C   diffusivity for temperature, salinity, and turbulence.
C   JTI May 31, 1995
C
C SUBROUTINE GALPERINQ4
C   Identical procedure to GALPERIN3, except that GALPERINQ4 is called
C   for every horizontal node (i.e. the node loop is outside of the
C   subroutine now).
C   CEN August 14, 1996.
C
C References: 
C   Mellor and Yamada (1982) Development of a turb. closure
C     model for geophys. fluid probs.  In: Rev. Geophys. and
C     Space Phys., vol. 20, 851-875.
C   Blumberg and Mellor (1987) A description of a 3D coastal ocean
C     ocean circulation model.  In: Three-Dimensional Coastal Ocean
C     Models; AGU Geophysical Monograph Series, Coastal and Estuarine 
C     Sciences 4, 1-16.
C   Galperin, Kantha, Hassid, Rosati (1988) A quasi-equilibrium 
C     turbulent energy model for geophysical flows.  
C     J. of the Atmospheric Sciences 45, 1, 55-62. (1988)
C-----------------------------------------------------------------------
C
C Declare integer program parameters via INCLUDE file
C
      INCLUDE 'QUODDY.DIM'
      PARAMETER(NEVDIM=NNVDIM-1)
C
C Global dimensioning and assignments
C
      DIMENSION UI(*),VI(*),ZI(*),RHOI(*)
      DIMENSION Q2I(*),Q2LI(*),EKM(*),EKH(*),EKQ(*)
      DIMENSION RHSQ2(*),RHSQ2L(*),DKQ2(*),DKQ2L(*)
C
C Local dimensioning and assignments
C
      DIMENSION ELL(NNVDIM),WPF(NNVDIM),GM(NEVDIM),GH(NEVDIM),
     *SM(NEVDIM),SH(NEVDIM)
      DIMENSION RHSQ2E(NEVDIM),RHSQ2LE(NEVDIM)
C
C  "Universal" constants from Mellor/Yamada
C
      DATA a1,a2,b1,b2,c1,e1,e2,e3,sq/0.92,0.74,16.6,10.1,0.08,
     &                              1.8,1.33,0.25,0.2/
C
C  von Karman constant
C
      data VK/0.4/
C
C  Minimum viscosity, diffusivities (JTI, 31 May, 1995)
C
      EKMMIN=EKMIN
      EKHMIN=EKMIN
      EKQMIN=EKMIN
C
C  Galperin constants (DRL notes, 24 Sept. 1993)
C
C         g0=(1.-6.*a1/b1) 
C         g1=6.*a1+b2        
C         g2=a1*(g0-3.*c1)
C         g3=3.*a1*a2*( (b2-3.*a2)*g0 - 3.*c1*g1 )
C         g4=3.*a2*g1
C         g5=9.*a1*a2
C         g6=a2*g0
C         g7=g4
C
      data g0,g1,g2,g3,g4,g5,g6,g7
     &/.66747,15.620,.39327,3.0858,34.676,6.1272,.49393,34.676/
C
C BEGINNING OF EXECUTABLE ROUTINE
C
C Constrain Q2L by limitation imposed in Galperin et al 1988, eqn. 22
C
      DO J=1,NNV
         JP=MIN0(J+1,NNV)
         JM=MAX0(J-1,1)
         DRHODZ=(RHOI(JP)-RHOI(JM))/(ZI(JP)-ZI(JM))
         IF(DRHODZ.LT.0.)THEN
            BVF=SQRT(-G*DRHODZ)
            FTR=(Q2I(J)**(1.5))*.53/BVF
            Q2LI(J)=AMIN1(Q2LI(J),FTR) 
         ENDIF
      ENDDO
C
C Compute nodal values of mixing length and wall-proximity-function
C via BGO92 and incorporating that:
C   The true bottom is (bathmetry + ZLOGBOT) 
C   The true top  is (elevation + ZLOGTOP)
C
      DO I=1,NNV
         Q2I(I)=AMAX1(Q2I(I),Q2MIN)
         ELL(I)=AMAX1(Q2LI(I)/Q2I(I),ELLMIN)
         TEMP=ELL(I)/VK
         TEMP2=TEMP/((ZI(I)-ZI(1))+ZLOGBOT)
         TEMP3=TEMP/(ZI(NNV)-ZI(I)+ZLOGTOP)
         WPF(I)=1.+E2*TEMP2*TEMP2+E3*TEMP3*TEMP3
      ENDDO
C
C Element based computation
C
      DO I=1,NEV
C
C  Production terms due to shear (GM) and buoyancy (GH)
C                   Note GM not used in Galperin
C
         ELL2TQ2=ELL(I+1)*ELL(I+1)/Q2I(I+1)
         ELL2BQ2=ELL(I)*ELL(I)/Q2I(I)
         FACTOR=0.5*(ELL2TQ2+ELL2BQ2)
         DZ=ZI(I+1)-ZI(I)
         DUDZ=(UI(I+1)-UI(I))/DZ
         DVDZ=(VI(I+1)-VI(I))/DZ
         SHEAR=DUDZ*DUDZ+DVDZ*DVDZ
         BUOY=g*(RHOI(I+1)-RHOI(I))/DZ
C         GM(I)=FACTOR*SHEAR
C
C  Upper limit on GH (Galperin '88, eqn. 29)
C
         GH(I)=AMIN1(FACTOR*BUOY,.0233)
C
C  Next compute SM and SH (Galperin et al 1988) 
C
         SM(I)=(g2-g3*GH(I))/((1.-g4*GH(I))*(1.-g5*GH(I)))
         SH(I)=g6/(1.-g7*GH(I))
C
C  Now compute the mixing coefficients and impose minimum values
C
         ELLQ=0.5*(ELL(I+1)*SQRT(Q2I(I+1))+ELL(I)*SQRT(Q2I(I)))
         EKM(I)=AMAX1(ELLQ*SM(I),EKMMIN)
         EKH(I)=AMAX1(ELLQ*SH(I),EKHMIN)
         EKQ(I)=AMAX1(ELLQ*sq   ,EKQMIN)
C
C  Finally, compute element-based right-hand side terms
C
         RHSQ2E(I)=2.*(EKM(I)*SHEAR+EKH(I)*BUOY)
         ELLe1=0.5*(ELL(I+1)+ELL(I))*e1
         RHSQ2LE(I)=ELLe1*(EKM(I)*SHEAR+EKH(I)*BUOY)
      ENDDO
C
C Compute node-based right hand side terms
C
      DO I=2,NNV-1
         DZT=ZI(I+1)-ZI(I)
         DZB=ZI(I)-ZI(I-1)
         DZJ=DZT+DZB
         RHSQ2(I)=(RHSQ2E(I)*DZT+RHSQ2E(I-1)*DZB)/DZJ
         RHSQ2L(I)=(RHSQ2LE(I)*DZT+RHSQ2LE(I-1)*DZB)/DZJ
      ENDDO
      RHSQ2(NNV)=RHSQ2E(NEV)
      RHSQ2(1)=RHSQ2E(1)
      RHSQ2L(NNV)=RHSQ2LE(NEV)
      RHSQ2L(1)=RHSQ2LE(1)
C
C Compute node=based dissipation/decay terms 
C      as quasilinear first-order decay for Q2, Q2L
C   
      DO I=1,NNV
         QOVERB1=SQRT(Q2I(I))/b1
         DKQ2(I)=2.*QOVERB1/ELL(I)
         DKQ2L(I)=QOVERB1*WPF(I)/ELL(I)
      ENDDO
C
C END OF ROUTINE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Vertical meshing:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE UNISIGMA (NNDIM, NN, NNV, BATHY, ELEVN, Z)
C-----------------------------------------------------------------------
C ASSIGNS UNIFORM VERTICAL MESH Z(I,J) UNDER EACH HORIZONRAL NODE
C  (UNIFORM DZ AMONG NNV NODES)
C TOP:     J=NNV,  Z(I,J) =  ELEVN(I)
C BOTTOM:  J=1,    Z(I,J) = -BATHY(I)
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF Z(I,J) IN MAIN.
C-----------------------------------------------------------------------
      DIMENSION Z(NNDIM,*),ELEVN(*),BATHY(*)
C
C Set uniform grid for FUNDY series
C
      DO I=1,NN
         DZ=BATHY(I)/(NNV-1.)
         ZJ=-BATHY(I)
         DO J=1,NNV
            Z(I,J)=ZJ
            ZJ=ZJ+DZ
         ENDDO
         Z(I,NNV)=0.
      ENDDO
C
C Adjust for QUODDY series
C
      DO I=1,NN
         ZETA=ELEVN(I)-BATHY(I)
         DO J=1,NNV
            ZFUNDY=Z(I,J)
            ZFACTR=(ZFUNDY+BATHY(I))/BATHY(I)
            Z(I,J)=ZFUNDY+ZETA*ZFACTR
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SINEGRID (NNDIM, NN, NNV, HDOWN, HMID, DZBL, Z)
C-----------------------------------------------------------------------
C Sinusoidal mesh and coefficients version of vertgrid
C-----------------------------------------------------------------------
C
C Global Arrays
C
      REAL Z(NNDIM,*), HDOWN(*), HMID(*)
C
C Set sinusoidal grid for the QUODDY series
C
      NEV=NNV-1
      PI=ACOS(-1.)
      DO I=1,NN
         ZETA=HMID(I)-HDOWN(I)
         A=(HDOWN(I)+ZETA-DZBL*NEV)/(NEV*SIN(2.*PI/NEV))
         IF(A .LT.0.)A=0.
         Z(I,1)=-HDOWN(I)
         DO J=2,NNV
            EPS=(J-1.0)/(NNV-1.0)
            Z(I,J)=-HDOWN(I)*(1.-EPS)-A*SIN(2.*PI*EPS)+EPS*ZETA
         ENDDO
         Z(I,NNV)=ZETA
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Convection calculation:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE CONVECTION(I,NNV,UZMID,VZMID,WZMID,ZMID,ZOLD,DELT,
     &           INCONV,DUZDX,DUZDY,DVZDX,DVZDY,DZDX,DZDY,VISCX,VISCY, 
     &           CONXZ,CONYZ,VERTFLX,CONX,CONY)
C-----------------------------------------------------------------------
C Computes RHS source terms for the wave equation and the momentum 
C   equation at node I.
C
C   CONXZ(I,J),CONYZ(I,J) are the 3-D advective plus diffusion terms
C     without FL*DV/DZ => CONXZ=(+UDUDX+VDUDY-Ah*DEL^2 U,ETC).
C     The vertical part is implicit in the tri-diagonal velocity matrix.
C   CONX(I),CONY(I) are the vertical integrals of the advective
C       plus diffusion terms for the wave equation
C       => CONX=INT[CONXZ+(FL+Wmesh)*DU/DZ]dz.
C   INCONV(I)=1 => CON* arrays include advection + horizontal diffusion
C   INCONV(I)=0 => CON* arrays include only horizontal diffusion
C                  (i.e. contributions from VISCX and VISCY)
C-----------------------------------------------------------------------
      INCLUDE 'QUODDY.DIM'
C
C Global arrays
C
      REAL UZMID(NNDIM,*),VZMID(NNDIM,*),WZMID(NNDIM,*),
     &ZMID(NNDIM,*),ZOLD(NNDIM,*),
     &DUZDX(*),DUZDY(*),DVZDX(*),DVZDY(*),DZDX(*),DZDY(*),
     &VISCX(*),VISCY(*),
     &CONXZ(NNDIM,*),CONYZ(NNDIM,*),VERTFLX(NNDIM,*),CONX(*),CONY(*)
      INTEGER INCONV(*)
C
C Local dimensioning and assignments
C
      REAL ZMIDI(NNVDIM),UZMIDI(NNVDIM),VZMIDI(NNVDIM),WZMIDI(NNVDIM)
C
C Initialize arrays
C
      CONX(I)=0.
      CONY(I)=0.
C
C Store *MID arrays in local 1-D arrays
C
      DO J=1,NNV
         ZMIDI(J)=ZMID(I,J)
         UZMIDI(J)=UZMID(I,J)
         VZMIDI(J)=VZMID(I,J)
         WZMIDI(J)=WZMID(I,J)
      ENDDO
C
C Perform computations
C
      DO J=1,NNV
         WMESH=(ZMIDI(J)-ZOLD(I,J))/DELT
         FL=WZMIDI(J)-WMESH-(UZMIDI(J)*DZDX(J)+VZMIDI(J)*DZDY(J))
         VERTFLX(I,J)=FL
         JP=MIN0(J+1,NNV)
         JM=MAX0(J-1,1)
         DELTAZ=0.5*(ZMIDI(JP)-ZMIDI(JM))
         DUDZ=(UZMIDI(JP)-UZMIDI(JM))/(ZMIDI(JP)-ZMIDI(JM))
         DVDZ=(VZMIDI(JP)-VZMIDI(JM))/(ZMIDI(JP)-ZMIDI(JM))
         CONVX=UZMIDI(J)*DUZDX(J)+VZMIDI(J)*DUZDY(J)
         CONVX=CONVX*INCONV(I)-VISCX(J)
         CONVY=UZMIDI(J)*DVZDX(J)+VZMIDI(J)*DVZDY(J)
         CONVY=CONVY*INCONV(I)-VISCY(J)
         CONXZ(I,J)=CONVX
         CONYZ(I,J)=CONVY
         CONX(I)=CONX(I)+DELTAZ*(CONVX+(FL+WMESH)*DUDZ*INCONV(I))
         CONY(I)=CONY(I)+DELTAZ*(CONVY+(FL+WMESH)*DVDZ*INCONV(I))
      ENDDO 
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTVEL3_2(I,NNV,UZMID,VZMID,WZMID,PME,SRC,SV,
     &                ZMID,ZOLD,DELT,DUZDX,DVZDY,DZDX,DZDY)
C-----------------------------------------------------------------------
c  Computes vertical velocity under a horizontal node I.
c  Procedure as in FUNDY5:
c     d/dz of continuity equation; 
c     Dirichlet BC's at top and bottom based on kinematic conditions.
C  Adapted from vertvel.f for Quoddy3.2
c     Adds sources as constants on vertical elements
c  Locally-defined tridiagonal matrix.
C-----------------------------------------------------------------------
      INCLUDE 'QUODDY.DIM'
C
C Global arrays
C
      REAL UZMID(NNDIM,*),VZMID(NNDIM,*),WZMID(NNDIM,*),
     &ZMID(NNDIM,*),ZOLD(NNDIM,*),
     &DUZDX(*),DVZDY(*),DZDX(*),DZDY(*)
      REAL PME(*),SRC(NNDIM,*),SV(*)
C
C Local dimensioning and assignments
C
      REAL ZMIDI(NNVDIM),UZMIDI(NNVDIM),VZMIDI(NNVDIM),SRCI(NNVDIM)
      REAL ACONT(NNVDIM),BCONT(NNVDIM), CCONT(NNVDIM),DCONT(NNVDIM)
C
C Store *MID arrays in local 1-D arrays
C
      DO J=1,NNV
         ZMIDI(J)=ZMID(I,J)
         UZMIDI(J)=UZMID(I,J)
         VZMIDI(J)=VZMID(I,J)
         SRCI(J)=SRC(I,J)
      ENDDO
C
C  Boundary Conditions tri-diagonal entries
C
      ACONT(1)=0.0
      BCONT(1)=1.0
      CCONT(1)=0.0
      DCONT(1)=UZMIDI(1)*DZDX(1)+VZMIDI(1)*DZDY(1)
      ACONT(NNV)=0.0
      BCONT(NNV)=1.0
      CCONT(NNV)=0.0
      DCONT(NNV)=(ZMIDI(NNV)-ZOLD(i,NNV))/DELT
     &          +UZMIDI(NNV)*DZDX(NNV)+VZMIDI(NNV)*DZDY(NNV)-PME(I)
C
C  Iternal tri-diagonal entries
C
      DO J=2,NNV-1
         DZTOP=ZMIDI(J+1)-ZMIDI(J)
         DU=UZMIDI(J+1)-UZMIDI(J)
         DV=VZMIDI(J+1)-VZMIDI(J)
         DIVTOP=DUZDX(J+1)-DZDX(J+1)*DU/DZTOP
     &         +DVZDY(J+1)-DZDY(J+1)*DV/DZTOP
         DZBOT=ZMIDI(J)-ZMIDI(J-1)
         DU=UZMIDI(J)-UZMIDI(J-1)
         DV=VZMIDI(J)-VZMIDI(J-1)
         DIVBOT=DUZDX(J-1)-DZDX(J-1)*DU/DZBOT
     &         +DVZDY(J-1)-DZDY(J-1)*DV/DZBOT
         ACONT(J)=1./DZBOT
         CCONT(J)=1./DZTOP
         BCONT(J)=-ACONT(J)-CCONT(J)
         DCONT(J)=0.5*(-DIVTOP+DIVBOT)
     &           +(SRCI(J)/DZTOP-SRCI(J-1)/DZBOT)/SV(I)
      ENDDO
C
C Solve and write results to 3-D array for W
C
      CALL THOMAS(3,ACONT,BCONT,CCONT,DCONT,NNV)
      DO J=1,NNV
         WZMID(I,J)=DCONT(J)
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Baroclinic forcing terms:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHOXYQ2(ITER,G,nn,ne,x,y,in,ppx,ppy,IQ,JQ,nnv,z,rho,
     &nlev,zl,rhox,rhoy,hrbarxe,hrbarye)
C-----------------------------------------------------------------------
C purpose: This routine computes baroclinic forcing terms from sigma 
C          surface values of rho
C          (pertubation density-reference density)/reference density.
C
C algorithmic design: This routine has the identical algorithmic design
C          as RHOXY2 (which was the routine implemented to compute
C          baroclinic pressure gradients for QUODDY2* and QUODDY3*).
C
C index range:  i=1,nn, j=1,nnv, l=1,nlev
C
C inputs:  ITER - iteration number (not used herein)
C          G - gravitational acceleration (m/sec^2)
C          nn - number of horizontal nodes
C          ne - number of horizontal elements
C          x(i),y(i) - x,y horizontal nodal coordinates
C          in(3,j) - triangular element incidence list
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          nnv - number of vertical nodes
C          z(i,j) - z coordinate of sigma mesh 3-D nodes
C          rho(i,j) - sigma mesh values of rho
C          nlev - the number of level surfaces
C          zl(l) - level surface z values
C
C outputs: rhox(i,j) = x-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhoy(i,j) = y-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhbarxe(k)= x-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C          rhbarye(k)= y-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           RHOXY2.f - October 6, 1993
C           RHOXYQ2.f - March 7, 1996 
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      include 'QUODDY.DIM'
C
C Global arrays
C  fem mesh
      integer in(3,nedim)
      real x(nndim),y(nndim)
C  sigma mesh
      real z(nndim,nnvdim),rho(nndim,nnvdim)
      real rhox(nndim,nnvdim),rhoy(nndim,nnvdim)
      real hrbarxe(nedim),hrbarye(nedim)
C  level surface mesh
      real zl(nlevdim)
C  Sparse matrix arrays 
      INTEGER IQ(nndim),JQ(nftrdim*nndim)
      REAL PPX(nftrdim*nndim),PPY(nftrdim*nndim)
C
C Local arrays
C
C  sigma mesh
      real zs(nnvdim),drdx(nnvdim),drdy(nnvdim)
C  level mesh 
      real rhol(nndim,nlevdim)
      real drdxl(nlevdim),drdyl(nlevdim),rhoxl(nlevdim),rhoyl(nlevdim)
C
C Interpolate sigma mesh rho data to level mesh
C
      do i=1,nn
         do j=1,nnv
            zs(j)=z(i,j)
         enddo
         do l=1,nlev
            CALL BASIS1D1(nnv,zs,zl(l),b2,n2,b1,n1)
C
C  notes:   zdes<zref(1)   => n1=nnv,b1=0.0,n2=1,b2=1.0 
C           zdes>zref(nnv) => n1=nnv,b1=1.0,n2=1,b2=0.0 
C
            rhol(i,l)=b1*rho(i,n1)+b2*rho(i,n2)
         enddo
      enddo
C
C Begin node loop to compute rhox and rhoy
C
      do i=1,nn
C
C  Compute level mesh nodal values of drdxl and drdyl, using 
C  appropriate parts of SUBROUTINE SPRSGRAD and SPRSPAK arrays
C
         KMIN=1
         IF(I.GT.1) KMIN=IQ(I-1)+1
         KMAX=IQ(I)
         DO J=1,nlev
            DRDXL(J)=0.0
            DRDYL(J)=0.0
            DO K=KMIN,KMAX
               DRDXL(J)=DRDXL(J)+PPX(K)*rhol(JQ(K),J)
               DRDYL(J)=DRDYL(J)+PPY(K)*rhol(JQ(K),J)
            ENDDO
         ENDDO
C
C  Interpolate drdxl and drdyl back to sigma mesh
C
         do j=1,nnv
            CALL BASIS1D1(nlev,zl,z(i,j),b2,n2,b1,n1)
C
C  notes:   zdes<zref(1)    => n1=nlev,b1=0.0,n2=1,b2=1.0 
C           zdes>zref(nlev) => n1=nlev,b1=1.0,n2=1,b2=0.0 
C
            drdx(j)=b1*drdxl(n1)+b2*drdxl(n2)
            drdy(j)=b1*drdyl(n1)+b2*drdyl(n2)
         enddo
C
C  Compute rhox and rhoy by direct integration on the sigma mesh
C
         rhox(i,nnv)=0.0
         rhoy(i,nnv)=0.0
         rhoxu=0.0
         rhoyu=0.0
         do j=nnv-1,1,-1
            dz=z(i,j+1)-z(i,j)
            fac=0.5*g*dz
            rhoxb=rhoxu+fac*(drdx(j+1)+drdx(j))
            rhoyb=rhoyu+fac*(drdy(j+1)+drdy(j))
            rhox(i,j)=rhoxb
            rhoy(i,j)=rhoyb
            rhoxu=rhoxb
            rhoyu=rhoyb
         enddo
      enddo
C
C Begin element loop to compute hrbarxe and hrbarye
C
      do k=1,ne
         i1=in(1,k)
         i2=in(2,k)
         i3=in(3,k)
         delx1=x(i2)-x(i3)
         delx2=x(i3)-x(i1)
         delx3=x(i1)-x(i2)
         dely1=y(i2)-y(i3)
         dely2=y(i3)-y(i1)
         dely3=y(i1)-y(i2)
         A2I=1.0/((x(i1)*dely1+x(i2)*dely2+x(i3)*dely3))
C        note:A2I=1/(2A)
C
C  Compute level mesh element values of drdxl and drdyl
C
         do l=1,nlev
            rho1=rhol(i1,l)
            rho2=rhol(i2,l)
            rho3=rhol(i3,l)
            drdxl(l)= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdyl(l)=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
         enddo
C
C  Interpolate drdxl and drdyl back to sigma mesh
C
         do j=1,nnv
            zs(j)=(z(i1,j)+z(i2,j)+z(i3,j))/3.0
            CALL BASIS1D1(nlev,zl,zs(j),b2,n2,b1,n1)
C
C  notes:   zdes<zref(1)    => n1=nlev,b1=0.0,n2=1,b2=1.0 
C           zdes>zref(nlev) => n1=nlev,b1=1.0,n2=1,b2=0.0 
C
            drdx(j)=b1*drdxl(n1)+b2*drdxl(n2)
            drdy(j)=b1*drdyl(n1)+b2*drdyl(n2)
         enddo
C
C  Compute element based hrbar by direct integration on the sigma mesh
C
         rhoxeu=0.0
         rhoyeu=0.0
         hrbarx=0.0
         hrbary=0.0
         do j=nnv-1,1,-1
            dz=zs(j+1)-zs(j)
            fac=0.5*g*dz
            rhoxeb=rhoxeu+fac*(drdx(j+1)+drdx(j))
            rhoyeb=rhoyeu+fac*(drdy(j+1)+drdy(j))
            hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
            hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
            rhoxeu=rhoxeb
            rhoyeu=rhoyeb
         enddo
         hrbarxe(k)=hrbarx
         hrbarye(k)=hrbary
      enddo
C
C End of routine
C
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHOXYQ2E(ITER,G,nn,ne,x,y,in,ppx,ppy,IQ,JQ,nnv,z,rho,
     &nlev,zl,rhox,rhoy,hrbarxe,hrbarye)
C-----------------------------------------------------------------------
C purpose: This routine computes baroclinic forcing terms from sigma 
C          surface values of rho
C          (pertubation density-reference density)/reference density.
C
C algorithmic design: This routine has the same basic algorithmic design
C          as RHOXY2 (which was the routine implemented to compute
C          baroclinic pressure gradients for QUODDY2* and QUODDY3*).
C          However, some efficiency measures have been implemented to
C          decrease the run time of this routine and the elemental
C          calculation of hrbar is conducted directly on the level
C          mesh.
C
C index range:  i=1,nn, j=1,nnv, l=1,nlev
C
C inputs:  ITER - iteration number (not used herein)
C          G - gravitational acceleration (m/sec^2)
C          nn - number of horizontal nodes
C          ne - number of horizontal elements
C          x(i),y(i) - x,y horizontal nodal coordinates
C          in(3,j) - triangular element incidence list
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          nnv - number of vertical nodes
C          z(i,j) - z coordinate of sigma mesh 3-D nodes
C          rho(i,j) - sigma mesh values of rho
C          nlev - the number of level surfaces
C          zl(l) - level surface z values
C
C outputs: rhox(i,j) = x-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhoy(i,j) = y-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhbarxe(k)= x-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C          rhbarye(k)= y-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           RHOXY2.f - October 6, 1993
C           RHOXYQ2E.f - March 7, 1996 
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      include 'QUODDY.DIM'
C
C Global arrays
C  fem mesh
      integer in(3,nedim)
      real x(nndim),y(nndim)
C  sigma mesh
      real z(nndim,nnvdim),rho(nndim,nnvdim)
      real rhox(nndim,nnvdim),rhoy(nndim,nnvdim)
      real hrbarxe(nedim),hrbarye(nedim)
C  level surface mesh
      real zl(nlevdim)
C  Sparse matrix arrays 
      INTEGER IQ(nndim),JQ(nftrdim*nndim)
      REAL PPX(nftrdim*nndim),PPY(nftrdim*nndim)
C
C Local arrays
C
C  sigma mesh
      real zs(nnvdim),drdx(nnvdim),drdy(nnvdim)
C  level mesh 
      real rhol(nndim,nlevdim)
      real drdxl(nlevdim),drdyl(nlevdim),rhoxl(nlevdim),rhoyl(nlevdim)
C
C Interpolate sigma mesh rho data to level mesh
C
      do 10 i=1,nn
         l=nlev
         do while(zl(l).GE.z(i,nnv))
            rhol(i,l)=rho(i,nnv)
            l=l-1
            if(l.EQ.0)go to 10
         enddo
         n1=nnv-1
         do while(n1.GE.1)
            n2=n1+1
            z2=z(i,n2)
            z1=z(i,n1)
            fac=1.0/(z2-z1)
            do while(zl(l).GE.z(i,n1))
               b1=fac*(z2-zl(l))
               b2=fac*(zl(l)-z1)
               rhol(i,l)=b1*rho(i,n1)+b2*rho(i,n2)
               l=l-1
               if(l.EQ.0)go to 10
            enddo
            n1=n1-1
         enddo
         do while(l.GE.1)
            rhol(i,l)=rho(i,1)
            l=l-1
         enddo
 10   continue
C
C Begin node loop to compute rhox and rhoy
C
      do i=1,nn
C
C  Compute level mesh nodal values of drdxl and drdyl, using 
C  appropriate parts of SUBROUTINE SPRSGRAD and SPRSPAK arrays
C
         KMIN=1
         IF(I.GT.1) KMIN=IQ(I-1)+1
         KMAX=IQ(I)
         l=nlev
         DO WHILE(zl(l).GE.z(i,1))
            DRDXL(l)=0.0
            DRDYL(l)=0.0
            DO K=KMIN,KMAX
               DRDXL(l)=DRDXL(l)+PPX(K)*rhol(JQ(K),l)
               DRDYL(l)=DRDYL(l)+PPY(K)*rhol(JQ(K),l)
            ENDDO
            l=l-1
            if(l.EQ.0)go to 20
         ENDDO
         DRDXL(l)=0.0
         DRDYL(l)=0.0
         DO K=KMIN,KMAX
            DRDXL(l)=DRDXL(l)+PPX(K)*rhol(JQ(K),l)
            DRDYL(l)=DRDYL(l)+PPY(K)*rhol(JQ(K),l)
         ENDDO
 20      CONTINUE

C
C  Interpolate level mesh drdxl and drdyl back to sigma mesh
C
         n2=2
         do j=1,nnv
            zhere=z(i,j)
            if(zhere.LE.zl(1))then
               drdx(j)=drdxl(1)
               drdy(j)=drdyl(1)
            else if(zhere.GE.zl(nlev))then
               drdx(j)=drdxl(nlev)
               drdy(j)=drdyl(nlev)
            else
 30            if(zhere.LE.zl(n2))then
                  n1=n2-1
                  z2=zl(n2)
                  z1=zl(n1)
                  fac=1.0/(z2-z1)
                  b1=fac*(z2-zhere)
                  b2=fac*(zhere-z1)
                  drdx(j)=b1*drdxl(n1)+b2*drdxl(n2)
                  drdy(j)=b1*drdyl(n1)+b2*drdyl(n2)
               else
                  n2=n2+1
                  go to 30
               endif
            endif
         enddo
C
C  Compute rhox and rhoy by direct integration on the sigma mesh
C
         rhox(i,nnv)=0.0
         rhoy(i,nnv)=0.0
         rhoxu=0.0
         rhoyu=0.0
         do j=nnv-1,1,-1
            dz=z(i,j+1)-z(i,j)
            fac=0.5*g*dz
            rhoxb=rhoxu+fac*(drdx(j+1)+drdx(j))
            rhoyb=rhoyu+fac*(drdy(j+1)+drdy(j))
            rhox(i,j)=rhoxb
            rhoy(i,j)=rhoyb
            rhoxu=rhoxb
            rhoyu=rhoyb
         enddo
      enddo
C
C Begin element loop to compute hrbarxe and hrbarye
C
      do k=1,ne
         i1=in(1,k)
         i2=in(2,k)
         i3=in(3,k)
         znnv=(z(i1,nnv)+z(i2,nnv)+z(i3,nnv))/3.0
         z1=(z(i1,1)+z(i2,1)+z(i3,1))/3.0
         delx1=x(i2)-x(i3)
         delx2=x(i3)-x(i1)
         delx3=x(i1)-x(i2)
         dely1=y(i2)-y(i3)
         dely2=y(i3)-y(i1)
         dely3=y(i1)-y(i2)
         A2I=1.0/((x(i1)*dely1+x(i2)*dely2+x(i3)*dely3))
C        note:A2I=1/(2A)
C
C  Compute level mesh element values of drdxl and drdyl from the top
C  level to the level just below the average depth for the element
C
         l=nlev
         DO WHILE(zl(l).GE.z1)
            rho1=rhol(i1,l)
            rho2=rhol(i2,l)
            rho3=rhol(i3,l)
            drdxl(l)= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdyl(l)=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
            l=l-1
            if(l.EQ.0)go to 40
         ENDDO 
         rho1=rhol(i1,l)
         rho2=rhol(i2,l)
         rho3=rhol(i3,l)
         drdxl(l)= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
         drdyl(l)=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
 40      CONTINUE
C
C  Compute element based hrbar by direct integration on level mesh
C
C   First find drdx,drdy at the top and bottom of sigma mesh
C
         CALL BASIS1D1(nlev,zl,znnv,b2,n2,b1,n1)
         maxlev=n1
         drdxnnv=b1*drdxl(n1)+b2*drdxl(n2)
         drdynnv=b1*drdyl(n1)+b2*drdyl(n2)
         CALL BASIS1D1(nlev,zl,z1,b2,n2,b1,n1)
         drdx1=b1*drdxl(n1)+b2*drdxl(n2)
         drdy1=b1*drdyl(n1)+b2*drdyl(n2)
         minlev=n2
C
C   Initialize with the portion of the integral between znnv and the
C   level surface just below znnv
C
         dz=znnv-zl(maxlev)
         fac=0.5*g*dz
         rhoxeb=fac*(drdxnnv+drdxl(maxlev))
         rhoyeb=fac*(drdynnv+drdyl(maxlev))
         hrbarx=.5*(rhoxeb)*dz
         hrbary=.5*(rhoyeb)*dz
         rhoxeu=rhoxeb
         rhoyeu=rhoyeb
C
C   Now add portion from level surface just below znnv to the level
C   surface just above z1
C
         do l=maxlev-1,minlev,-1
            dz=zl(l+1)-zl(l)
            fac=0.5*g*dz
            rhoxeb=rhoxeu+fac*(drdxl(l+1)+drdxl(l))
            rhoyeb=rhoyeu+fac*(drdyl(l+1)+drdyl(l))
            hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
            hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
            rhoxeu=rhoxeb
            rhoyeu=rhoyeb
         enddo
C
C   Finally, add the portion between level position zl(minlev) and z1 
C
         dz=zl(minlev)-z1
         fac=0.5*g*dz
         rhoxeb=rhoxeu+fac*(drdx1+drdxl(minlev))
         rhoyeb=rhoyeu+fac*(drdy1+drdyl(minlev))
         hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
         hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
         hrbarxe(k)=hrbarx
         hrbarye(k)=hrbary
      enddo
C
C End of routine
C
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHOXYQ2S(ITER,G,nn,ne,x,y,in,ppx,ppy,IQ,JQ,nnv,z,rho,
     &nlev,zl,rhox,rhoy,hrbarxe,hrbarye)
C-----------------------------------------------------------------------
C purpose: This routine computes baroclinic forcing terms from sigma 
C          surface values of rho
C          (pertubation density-reference density)/reference density.
C
C algorithmic design: This routine performs calculations directly on
C          the sigma mesh.
C
C index range:  i=1,nn, j=1,nnv, l=1,nlev
C
C inputs:  ITER - iteration number (not used herein)
C          G - gravitational acceleration (m/sec^2)
C          nn - number of horizontal nodes
C          ne - number of horizontal elements
C          x(i),y(i) - x,y horizontal nodal coordinates
C          in(3,j) - triangular element incidence list
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          nnv - number of vertical nodes
C          z(i,j) - z coordinate of sigma mesh 3-D nodes
C          rho(i,j) - sigma mesh values of rho
C          nlev - the number of level surfaces
C          zl(l) - level surface z values
C
C outputs: rhox(i,j) = x-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhoy(i,j) = y-directed sigma nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhbarxe(k)= x-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C          rhbarye(k)= y-directed sigma element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           RHOXY2S.f - October 6, 1993
C           RHOXYQ2S.f - August 15, 1996 
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      include 'QUODDY.DIM'
C
C Global arrays
C  fem mesh
      integer in(3,nedim)
      real x(nndim),y(nndim)
C  sigma mesh
      real z(nndim,nnvdim),rho(nndim,nnvdim)
      real rhox(nndim,nnvdim),rhoy(nndim,nnvdim)
      real hrbarxe(nedim),hrbarye(nedim)
C  level surface mesh
      real zl(nlevdim)
C  Sparse matrix arrays 
      INTEGER IQ(nndim),JQ(nftrdim*nndim)
      REAL PPX(nftrdim*nndim),PPY(nftrdim*nndim)
C
C Local arrays
C
C  sigma mesh
      real ze(nnvdim),rhoe(nnvdim),drdx(nnvdim),drdy(nnvdim)
      real dzdx(nnvdim),dzdy(nnvdim)
C
C Begin node loop to compute rhox and rhoy
C
      do 40 i=1,nn
C
C   Compute sigma mesh nodal values of drdx and drdy, using 
C      SUBROUTINE SPRSGRAD and SPRSPAK arrays
C
         CALL SPRSGRAD (I, PPX, PPY, IQ, JQ, RHO, Z, DRDX, DRDY,
     &                  DZDX, DZDY, NNV, NNDIM)
         DO 20 J=1,NNV
            JP = MIN0(J+1,NNV)
            JM = MAX0(J-1,1)
            DRDZ = (RHO(I,JP) - RHO(I,JM))/(Z(I,JP) - Z(I,JM))
            DRDX(J) = DRDX(J) - DZDX(J)*DRDZ
            DRDY(J) = DRDY(J) - DZDY(J)*DRDZ
   20    CONTINUE
C
C   Compute rhox and rhoy by direct integration on the sigma mesh
C
         rhox(i,nnv)=0.0
         rhoy(i,nnv)=0.0
         rhoxu=0.0
         rhoyu=0.0
         do 40 j=nnv-1,1,-1
            dz=z(i,j+1)-z(i,j)
            fac=0.5*g*dz
            rhoxb=rhoxu+fac*(drdx(j+1)+drdx(j))
            rhoyb=rhoyu+fac*(drdy(j+1)+drdy(j))
            rhox(i,j)=rhoxb
            rhoy(i,j)=rhoyb
            rhoxu=rhoxb
            rhoyu=rhoyb
 40   continue
C
C Begin element loop to compute hrbarxe and hrbarye
C
      do 80 k=1,ne
         i1=in(1,k)
         i2=in(2,k)
         i3=in(3,k)
         delx1=x(i2)-x(i3)
         delx2=x(i3)-x(i1)
         delx3=x(i1)-x(i2)
         dely1=y(i2)-y(i3)
         dely2=y(i3)-y(i1)
         dely3=y(i1)-y(i2)
         A2I=1.0/((x(i1)*dely1+x(i2)*dely2+x(i3)*dely3))
C        note:A2I=1/(2A)
C
C   Compute sigma mesh element values of drdx and drdy
C
         do 45 j=1,nnv
            ze(j)=(z(i1,j)+z(i2,j)+z(i3,j))/3.0
            rhoe(j)=(rho(i1,j)+rho(i2,j)+rho(i3,j))/3.0
 45      continue
         do 50 j=1,nnv
            rho1=rho(i1,j)
            rho2=rho(i2,j)
            rho3=rho(i3,j)
            z1=z(i1,j)
            z2=z(i2,j)
            z3=z(i3,j)
            drdx(j)= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdy(j)=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
C
C   Correct for tilt of sigma surfaces
C
            dzdxe= (z1*dely1+z2*dely2+z3*dely3)*A2I
            dzdye=-(z1*delx1+z2*delx2+z3*delx3)*A2I
            JP=MIN0(J+1,NNV)
            JM=MAX0(J-1,1)
            DRDZ= ( rhoe(JP)-rhoe(JM) ) / (  ze(JP) - ze(JM) )
            DRDX(J)=DRDX(J)-DZDXe*DRDZ
            DRDY(J)=DRDY(J)-DZDYe*DRDZ
 50      continue
C
C   Compute element based hrbar by direct integration on the sigma mesh
C
         rhoxeu=0.0
         rhoyeu=0.0
         hrbarx=0.0
         hrbary=0.0
         do 70 j=nnv-1,1,-1
            dz=ze(j+1)-ze(j)
            fac=0.5*g*dz
            rhoxeb=rhoxeu+fac*(drdx(j+1)+drdx(j))
            rhoyeb=rhoyeu+fac*(drdy(j+1)+drdy(j))
            hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
            hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
            rhoxeu=rhoxeb
            rhoyeu=rhoyeb
 70      continue
      hrbarxe(k)=hrbarx
      hrbarye(k)=hrbary
 80   continue
C
C End of routine
C
      return
      end
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHOXYQ4(ITER,G,nn,ne,x,y,in,ppx,ppy,IQ,JQ,nnv,z,rho,
     &nlev,zl,rhox,rhoy,hrbarxe,hrbarye)
C-----------------------------------------------------------------------
C purpose: This routine computes baroclinic forcing terms from sigma 
C          surface values of rho
C          (pertubation density-reference density)/reference density.
C
C index range:  i=1,nn, j=1,nnv, k=1,ne, l=1,nlev
C
C inputs:  ITER - iteration number (not used herein)
C          G - gravitational acceleration (m/sec^2)
C          nn - number of horizontal nodes
C          ne - number of horizontal elements
C          x(i),y(i) - x,y horizontal nodal coordinates
C          in(3,j) - triangular element incidence list
C          PPX(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix
C          PPY(I*NFTR) - SPRSPAK STANDARD Normalized Derivative matrix 
C          IQ(I),JQ(I*NFTR) - SPRSPAK STANDARD pointers for PPX and PPY
C          nnv - number of vertical nodes
C          z(i,j) - z coordinate of sigma mesh 3-D nodes
C          rho(i,j) - sigma mesh values of rho
C          nlev - the number of level surfaces
C          zl(l) - level surface z values
C
C outputs: rhox(i,j) = x-directed z nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhoy(i,j) = y-directed z nodal based baroclinic forcing
C                        positive on same side of momentum equ as dV/dt
C          rhbarxe(k)= x-directed z element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C          rhbarye(k)= y-directed z element based vertical average 
C                        of the depth times the baroclinic forcing,
C                        positive on same side of momentum equ as dV/dt
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           1993 OCT 06 - RHOXY2
C           1995 JUL 13 - RHOXY81
C           1996 AUG 15 - RHOXYQ4
C           1997 APR 28 - RHOXYQ4 REVISION:
C             Corrected minor bugs related to vertical interpolation. 
C             These bugs were discovered by QUG members studying new
C             domains.  The bugs do not effect earlier Gulf of Maine
C             simulations
C-----------------------------------------------------------------------
C
C Include file for parameter statements
      include 'QUODDY.DIM'
C
C Global arrays
C  fem mesh
      integer in(3,nedim)
      real x(nndim),y(nndim)
C  sigma mesh
      real z(nndim,nnvdim),rho(nndim,nnvdim)
      real rhox(nndim,nnvdim),rhoy(nndim,nnvdim)
      real hrbarxe(nedim),hrbarye(nedim)
C  level surface mesh
      real zl(nlevdim)
C  Sparse matrix arrays 
      INTEGER IQ(nndim),JQ(nftrdim*nndim)
      REAL PPX(nftrdim*nndim),PPY(nftrdim*nndim)
C
C Local arrays
C  fem mesh
      real bathy(nndim),dneigh(nndim)
C  z mesh
      real zs(nnvdim),drdx(nnvdim),drdy(nnvdim)
C  level mesh 
      real rhol(nndim,nlevdim)
      real drdxl(nlevdim),drdyl(nlevdim),rhoxl(nlevdim),rhoyl(nlevdim)
C  Sparse matrix arrays 
      REAL PPX1(NFTRDIM*NNDIM), PPY1(NFTRDIM*NNDIM)
C  Save necessary local information
      save PPX1,PPY1,dneigh
C
C Build nonstandard sprspak arrays on first call
      if(ITER.EQ.0)then
         do 5 i=1,nn
            bathy(i)=-z(i,1)
 5       continue
         CALL CALCDNEIGH(nndim,nedim,x,y,in,bathy,nn,ne,DNEIGH)
         CALL BUILDPPX1PPY1(NN,NE,X,Y,IN,BATHY,IQ,JQ,PPX1,PPY1)
      endif
C
C Interpolate z mesh rho data to level mesh from the top of the
C  level mesh to the level mesh node just above the bottom of the z
C  mesh.  Set level mesh values of rho equal to zero for level surfaces
C  below the bottom of the z mesh.
      do 10 i=1,nn
         l=nlev
         do 6 while(zl(l).GE.z(i,nnv))
            rhol(i,l)=rho(i,nnv)
            l=l-1
            if(l.EQ.0)go to 10
 6       continue
         n1=nnv-1
         do 8 while(n1.GE.1)
            n2=n1+1
            z2=z(i,n2)
            z1=z(i,n1)
            fac=1.0/(z2-z1)
            do 7 while(zl(l).GE.z(i,n1))
               b1=fac*(z2-zl(l))
               b2=fac*(zl(l)-z1)
               rhol(i,l)=b1*rho(i,n1)+b2*rho(i,n2)
               l=l-1
               if(l.EQ.0)go to 10
 7          continue
            n1=n1-1
 8       continue
         do 9 while(l.GE.1)
            rhol(i,l)=0.0
            l=l-1
 9       continue
 10   continue
C
C Begin node loop to compute rhox and rhoy
      do 40 i=1,nn
         dminnei=dneigh(i)
         zbot=z(i,1)
C
C  Compute level mesh nodal values of drdxl and drdyl, using 
C   appropriate parts of SUBROUTINE SPRSGRAD and SPRSPAK arrays
         KMIN=1
         IF(I.GT.1) KMIN=IQ(I-1)+1
         KMAX=IQ(I)
         l=nlev
C
C   From the surface to level surface just above dminnei
         do 20 WHILE(zl(l).GE.dminnei)
            DRDXL(l)=0.0
            DRDYL(l)=0.0
            do 19 K=KMIN,KMAX
               DRDXL(l)=DRDXL(l)+PPX(K)*rhol(JQ(K),l)
               DRDYL(l)=DRDYL(l)+PPY(K)*rhol(JQ(K),l)
 19         continue
            l=l-1
            IF(l.EQ.0)GO TO 24
 20      continue
C
C   From the level surface just above dminnei to level surface just
C    above zbot
         do 22 WHILE(zl(l).GE.zbot)
            DRDXL(l)=0.0
            DRDYL(l)=0.0
            do 21 K=KMIN,KMAX
               DRDXL(l)=DRDXL(l)+PPX1(K)*rhol(JQ(K),l)
               DRDYL(l)=DRDYL(l)+PPY1(K)*rhol(JQ(K),l)
 21         continue
            l=l-1
            IF(l.EQ.0)GO TO 24
 22      continue
 24      minlev=l+1
C
C  Compute drdx and drdy at the bottom node of the z mesh
         DRDX(1)=0.0
         DRDY(1)=0.0
         do 25 K=KMIN,KMAX
            ii=JQ(K)
            rhohere=0.0
            if(zbot.GE.z(ii,1))then
               do 26 j=2,nnv 
                  if(zbot.LE.z(ii,j))then
                      n2=j
                      n1=j-1
                      z1=z(ii,n1)
                      z2=z(ii,n2)
                      zdiff=z2-z1
                      b1=(z2-zbot)/zdiff
                      b2=(zbot-z1)/zdiff
                      rhohere=b1*rho(ii,n1)+b2*rho(ii,n2)
                      go to 27
                  endif
 26            continue
               rhohere=rho(ii,nnv)
            endif
 27         continue
            DRDX(1)=DRDX(1)+PPX1(K)*rhohere
            DRDY(1)=DRDY(1)+PPY1(K)*rhohere
 25      continue
C
C   Interpolate level mesh drdxl and drdyl back to z mesh
         j=2
         zhere=z(i,j)
         do 28 while(zhere.LE.zl(minlev))
            z2=zl(minlev)
            fac=1.0/(z2-zbot)
            b1=fac*(z2-zhere)
            b2=fac*(zhere-zbot)
            drdx(j)=b1*drdx(1)+b2*drdxl(minlev)
            drdy(j)=b1*drdy(1)+b2*drdyl(minlev)
            if(j.EQ.nnv)go to 33
            j=j+1
            zhere=z(i,j)
 28      continue
         n2=minlev+1
         do 29 while(zhere.LE.zl(nlev))
            do 31 while(zhere.LE.zl(n2))
               n1=n2-1
               z2=zl(n2)
               z1=zl(n1)
               fac=1.0/(z2-z1)
               b1=fac*(z2-zhere)
               b2=fac*(zhere-z1)
               drdx(j)=b1*drdxl(n1)+b2*drdxl(n2)
               drdy(j)=b1*drdyl(n1)+b2*drdyl(n2)
               if(j.EQ.nnv)go to 33
               j=j+1
               zhere=z(i,j)
 31         continue
            n2=n2+1
 29      continue
         do 32 jj=j,nnv
            drdx(jj)=drdxl(nlev)
            drdy(jj)=drdyl(nlev)
 32      continue
 33      continue
C
C   Compute rhox and rhoy by direct integration on the z mesh
         rhox(i,nnv)=0.0
         rhoy(i,nnv)=0.0
         rhoxu=0.0
         rhoyu=0.0
         do 40 j=nnv-1,1,-1
            dz=z(i,j+1)-z(i,j)
            fac=0.5*g*dz
            rhoxb=rhoxu+fac*(drdx(j+1)+drdx(j))
            rhoyb=rhoyu+fac*(drdy(j+1)+drdy(j))
            rhox(i,j)=rhoxb
            rhoy(i,j)=rhoyb
            rhoxu=rhoxb
            rhoyu=rhoyb
 40   continue
C
C Begin element loop to compute hrbarxe and hrbarye
      do 70 k=1,ne
         i1=in(1,k)
         i2=in(2,k)
         i3=in(3,k)
         znnv=(z(i1,nnv)+z(i2,nnv)+z(i3,nnv))/3.0
         z1=(z(i1,1)+z(i2,1)+z(i3,1))/3.0
         zdepthmin=AMAX1(z(i1,1),z(i2,1),z(i3,1))
         delx1=x(i2)-x(i3)
         delx2=x(i3)-x(i1)
         delx3=x(i1)-x(i2)
         dely1=y(i2)-y(i3)
         dely2=y(i3)-y(i1)
         dely3=y(i1)-y(i2)
         A2I=1.0/((x(i1)*dely1+x(i2)*dely2+x(i3)*dely3))
C        note:A2I=1/(2A)
C
C  Compute level mesh element values of drdxl and drdyl from the top
C   level to the level just above the shallowest node on the element
         l=nlev
         do 50 WHILE(zl(l).GE.zdepthmin)
            rho1=rhol(i1,l)
            rho2=rhol(i2,l)
            rho3=rhol(i3,l)
            drdxl(l)= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdyl(l)=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
            l=l-1
            if(l.EQ.0)go to 51
 50      continue 
 51      continue
         minlev=l+1
C
C  Compute element based hrbar by direct integration on level mesh
C
C   Find drdx,drdy at the top of the z mesh
         CALL BASIS1D1(nlev,zl,znnv,b2,n2,b1,n1)
         maxlev=n1
         drdxnnv=b1*drdxl(n1)+b2*drdxl(n2)
         drdynnv=b1*drdyl(n1)+b2*drdyl(n2)
C
C   Find drdx,drdy at the bottom of the z mesh
C
C    First find minimum, intermediate and maximum depth nodes on element
         ndepthmin=i1
         zdepthmin=z(i1,1)
         ndepthmax=i2
         zdepthmax=z(i2,1)
         if(z(i2,1).GT.zdepthmin)then
            ndepthmin=i2
            zdepthmin=z(i2,1)
            ndepthmax=i1
            zdepthmax=z(i1,1)
         endif
         if(z(i3,1).GT.zdepthmin)then
            ndepthmid=ndepthmin
            zdepthmid=zdepthmin
            ndepthmin=i3
            zdepthmin=z(i3,1)
         else if(z(i3,1).LT.zdepthmax)then
            ndepthmid=ndepthmax
            zdepthmid=zdepthmax
            ndepthmax=i3
            zdepthmax=z(i3,1)
         else
            ndepthmid=i3
            zdepthmid=z(i3,1)
         endif
C
C    Now compute drdx1,drdy1 based which of this elements nodes
C     have depths exceeding the average depth on the element
C
C     If all local nodes are wet, use standard element based calculation
         if(zdepthmin.EQ.z1)then
            rho1=rho(i1,1)
            rho2=rho(i2,1)
            rho3=rho(i3,1)
            drdx1= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdy1=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
C
C     If 2 local nodes are wet, use element based calculation which
C      includes wet local nodes plus center of element
         else if(zdepthmid.LE.z1)then
C      Center of element
            x1=(x(i1)+x(i2)+x(i3))/3.0
            y1=(y(i1)+y(i2)+y(i3))/3.0
            rho1=(rho(i1,1)+rho(i2,1)+rho(i3,1))/3.0
C      Intermediate depth node
            ii=ndepthmid
            do 54 j=2,nnv 
               if(z1.LT.z(ii,j))then
                  n2=j
                  n1=j-1
                  zdiff=z(ii,n2)-z(ii,n1)
                  b1=(z(ii,n2)-z1)/zdiff
                  b2=(z1-z(ii,n1))/zdiff  
                  go to 55
               endif
 54         continue
            b1=1.0
            n1=nnv
            b2=0.0
            n2=nnv-1
 55         continue
            x2=x(ii)
            y2=y(ii)
            rho2=b1*rho(ii,n1)+b2*rho(ii,n2)
C      Maximum depth node
            ii=ndepthmax
            do 52 j=2,nnv 
               if(z1.LT.z(ii,j))then
                  n2=j
                  n1=j-1
                  zdiff=z(ii,n2)-z(ii,n1)
                  b1=(z(ii,n2)-z1)/zdiff
                  b2=(z1-z(ii,n1))/zdiff  
                  go to 53
               endif
 52         continue
            b1=1.0
            n1=nnv
            b2=0.0
            n2=nnv-1
 53         continue
            x3=x(ii)
            y3=y(ii)
            rho3=b1*rho(ii,n1)+b2*rho(ii,n2)
C      Calculate drdx1,drdy1
            delx1=x2-x3
            delx2=x3-x1
            delx3=x1-x2
            dely1=y2-y3
            dely2=y3-y1
            dely3=y1-y2
            A2I=1.0/((x1*dely1+x2*dely2+x3*dely3))
            drdx1= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdy1=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
C
C     If only one local node is wet, find wet locations along element 
C      sides and use 2-D derivative based on the scaled down element 
C      which includes the wet node plus these new nodes
         else
C      Node between minimum and maximum depth nodes
            zfac=(z1-zdepthmax)/(zdepthmin-zdepthmax)
            x1=x(ndepthmax)+zfac*(x(ndepthmin)-x(ndepthmax))
            y1=y(ndepthmax)+zfac*(y(ndepthmin)-y(ndepthmax))
            rho1=rho(ndepthmax,1)
     &                     +zfac*(rho(ndepthmin,1)-rho(ndepthmax,1))
C      Node between middle and maximum depth nodes
            zfac=(z1-zdepthmax)/(zdepthmid-zdepthmax)
            x2=x(ndepthmax)+zfac*(x(ndepthmid)-x(ndepthmax))
            y2=y(ndepthmax)+zfac*(y(ndepthmid)-y(ndepthmax))
            rho2=rho(ndepthmax,1)
     &                     +zfac*(rho(ndepthmid,1)-rho(ndepthmax,1))
C      Maximum depth node
            ii=ndepthmax
            do 56 j=2,nnv 
               if(z1.LT.z(ii,j))then
                  n2=j
                  n1=j-1
                  zdiff=z(ii,n2)-z(ii,n1)
                  b1=(z(ii,n2)-z1)/zdiff
                  b2=(z1-z(ii,n1))/zdiff  
                  go to 57
               endif
 56         continue
            b1=1.0
            n1=nnv
            b2=0.0
            n2=nnv-1
 57         continue
            x3=x(ii)
            y3=y(ii)
            rho3=b1*rho(ii,n1)+b2*rho(ii,n2)
C      Calculate drdx1,drdy1
            delx1=x2-x3
            delx2=x3-x1
            delx3=x1-x2
            dely1=y2-y3
            dely2=y3-y1
            dely3=y1-y2
            A2I=1.0/((x1*dely1+x2*dely2+x3*dely3))
            drdx1= (rho1*dely1+rho2*dely2+rho3*dely3)*A2I
            drdy1=-(rho1*delx1+rho2*delx2+rho3*delx3)*A2I
         endif
C
C   Initialize with the portion of the integral between znnv and the
C    level surface just below znnv
         dz=znnv-zl(maxlev)
         fac=0.5*g*dz
         rhoxeb=fac*(drdxnnv+drdxl(maxlev))
         rhoyeb=fac*(drdynnv+drdyl(maxlev))
         hrbarx=.5*(rhoxeb)*dz
         hrbary=.5*(rhoyeb)*dz
         rhoxeu=rhoxeb
         rhoyeu=rhoyeb
C
C   Now add portion from level surface just below znnv to the level
C    surface just above zdepthmin
         do 60 l=maxlev-1,minlev,-1
            dz=zl(l+1)-zl(l)
            fac=0.5*g*dz
            rhoxeb=rhoxeu+fac*(drdxl(l+1)+drdxl(l))
            rhoyeb=rhoyeu+fac*(drdyl(l+1)+drdyl(l))
            hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
            hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
            rhoxeu=rhoxeb
            rhoyeu=rhoyeb
 60      continue
C
C   Now add the portion between level position zl(minlev) and z1 
         dz=zl(minlev)-z1
         fac=0.5*g*dz
         rhoxeb=rhoxeu+fac*(drdx1+drdxl(minlev))
         rhoyeb=rhoyeu+fac*(drdy1+drdyl(minlev))
         hrbarx=hrbarx+.5*(rhoxeb+rhoxeu)*dz
         hrbary=hrbary+.5*(rhoyeb+rhoyeu)*dz
         hrbarxe(k)=hrbarx
         hrbarye(k)=hrbary
 70   continue
C
C End of routine
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE CALCDNEIGH(nndim,nedim,x,y,in,bathy,nn,ne,DNEIGH)
C-----------------------------------------------------------------------
C purpose: This subroutine calculates a nodal array which contains
C          the minimum depth of local nodes (i.e. neighbors)
C-----------------------------------------------------------------------
C Global arrays
C
      integer in(3,nedim)
      real x(nndim),y(nndim),bathy(nndim),DNEIGH(nndim)
C
C Initialize DNEIGH
C
      do i=1,nn
         DNEIGH(i)=-bathy(i)
      enddo
C
C Fill DNEIGH
C
      do k=1,ne
         n1=in(1,k)
         n2=in(2,k)
         n3=in(3,k)
         botmin=AMAX1(-bathy(n1),-bathy(n2),-bathy(n3))
         DNEIGH(n1)=AMAX1(DNEIGH(n1),botmin)
         DNEIGH(n2)=AMAX1(DNEIGH(n2),botmin)
         DNEIGH(n3)=AMAX1(DNEIGH(n3),botmin)
      enddo
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE BUILDPPX1PPY1(NN,NE,X,Y,IN,BATHY,
     &IQ,JQ,PPX1,PPY1)
C-----------------------------------------------------------------------
C Purpose:  To compute SPRSPAK STANDARD Normalized Derivative matircies
C              for use under special conditions of RHOXYQ4 (i.e. only
C              use wet elements for the computation).
C
C Index range: I=1,NN, J=1,NNV, L=1,NE, K=1,NBE  
C
C inputs:  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
C outputs: IQ(I),JQ(I*NFTR)-SPRSPAK STANDARD pointers for PPX1 and PPY1
C          PPX1(I*NFTR)-SPRSPAK NONSTANDARD Normalized Derivative matrix
C          PPY1(I*NFTR)-SPRSPAK NONSTANDARD Normalized Derivative matrix 
C
C History: Developed by Daniel R. Lynch and Christopher E. Naimie
C          Dartmouth College
C          FUNDY4 - Part of Main Routine written by Daniel R. Lynch
C          FUNDY5 - Converted to a Subroutine by Christopher. E. Naimie
C-----------------------------------------------------------------------
C
C Include file for parameter statements
C
      INCLUDE 'QUODDY.DIM' 
C
C Global Arrays 
C
C   Sent to routine and returned unaltered
C
      INTEGER IN(3,NEDIM)
      REAL X(NNDIM),Y(NNDIM),BATHY(NNDIM)
C
C   Computed in routine and returned to calling routine
C
      INTEGER IQ(NNDIM),JQ(NFTRDIM*NNDIM)
      REAL PPX1(NFTRDIM*NNDIM),PPY1(NFTRDIM*NNDIM)
C
C Local arrays
C
      REAL DELX(3),DELY(3)
      REAL SV(NNDIM)
C
C Initialize the sprspak arrays PPX1, PPY1
C
      DO 10 I=1,NFTRDIM*NNDIM
         PPX1(I)=0.
         PPY1(I)=0.
 10   CONTINUE
      DO 20 I=1,NN
         SV(I)=0.0
 20   CONTINUE
C
C Compute diagonal mass matrix SV(I)=<Pi> and derivative matricies in 
C   sprspak form
C
C     PPX1 = <Pi*dP/dxj>
C     PPY1 = <Pi*dP/dyj>
C
      DO 30 L=1,NE
         I1=IN(1,L)
         I2=IN(2,L)
         I3=IN(3,L)
         DELX(1)=X(I2)-X(I3)
         DELX(2)=X(I3)-X(I1)
         DELX(3)=X(I1)-X(I2)
         DELY(1)=Y(I2)-Y(I3)
         DELY(2) =Y(I3)-Y(I1)
         DELY(3)=Y(I1)-Y(I2)
         ARR=0.5*(X(I1)*DELY(1) + X(I2)*DELY(2) + X(I3)*DELY(3))
         BMIN=AMIN1(BATHY(I1),BATHY(I2),BATHY(I3))
         DO 30 I=1,3
            II=IN(I,L)
            IF(BATHY(II).EQ.BMIN)THEN
               SV(II)=SV(II)+ARR/3.
               DO 29 J=1,3
                  IJ=IN(J,L)
                  K=KAY(II,IJ,IQ,JQ)
                  PPX1(K)=PPX1(K)+DELY(J)/6.
                  PPY1(K)=PPY1(K)-DELX(J)/6.
 29            CONTINUE
            ENDIF
 30   CONTINUE
C
C Premulitply PPX1 and PPY1 by inverse of SV so that now,
C 
C     PPX1 = <Pi*dP/dxj>/<Pi>
C     PPY1 = <Pi*dP/dyj>/<Pi>
C
      CALL SPRSINVMLT1(PPX1,IQ,SV,NN)
      CALL SPRSINVMLT1(PPY1,IQ,SV,NN)
C
C Normalize 
C
      CALL SPRSNORM1(PPX1,IQ,NN)
      CALL SPRSNORM1(PPY1,IQ,NN)
      RETURN
      END
C***********************************************************************
C*****************************************************************FUNDY5
      SUBROUTINE SPRSINVMLT1(QV,IQ,D,NN)
C-----------------------------------------------------------------------
C  REPLACES QV WITH INV(D)*QV
C  D IS A DIAGONAL MATRIX, LENGTH NN 
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C
C  This routine is a variation of the SPRSPAK routine SPRSINVMLT.  
C  SPRSINVMLT was converted to SPRSINVMLT1 by Christopher E. Naimie,
C  specifically for use with SUBROUTINE RHOXYQ4.  The modification
C  is to ignore those entries for which D equals zero.
C-----------------------------------------------------------------------
      DIMENSION QV(*),D(*),IQ(*)
      KMAX=0
      DO 10 I=1,NN
         KMIN=KMAX+1
         KMAX=IQ(I)
         IF(D(I).NE.0.0)THEN
            FTR=1./D(I)
            DO 20 K=KMIN,KMAX
               QV(K)=FTR*QV(K)
 20          CONTINUE
         ENDIF
 10   CONTINUE
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE SPRSNORM1(QV,IQ,NN)
C-----------------------------------------------------------------------
C  CORRECTS QV SO THAT ROW SUMS ARE ZERO
C  THE CORRECTION IS SPREAD EVENLY OVER THE WHOLE ROW
C  NN IS THE SIZE OF Q (# OF EQUATIONS)
C  IQ(I) IS THE ADDRESS IN QV OF THE LAST NONZERO ENTRY IN ROW I OF Q
C
C  This routine is a variation of the SPRSPAK routine SPRSNORM.  
C  SPRSNORM was converted to SPRSNORM1 by Christopher E. Naimie,
C  specifically for use with SUBROUTINE RHOXYQ4.  The modification
C  is to ignore those entries for which QV equals zero.
C-----------------------------------------------------------------------
      DIMENSION QV(*),IQ(*)
      KMAX=0
      DO 10 I=1,NN
      KMIN=KMAX+1
      KMAX=IQ(I)
      SUM=0.0
      NSUM=0
      DO 20 K=KMIN,KMAX
         IF(QV(K).NE.0.0)THEN
            SUM=SUM+QV(K)
            NSUM=NSUM+1
         ENDIF
   20 CONTINUE
      IF(NSUM.NE.0)THEN
         CORR=SUM/REAL(NSUM)
         DO 25 K=KMIN,KMAX
            IF(QV(K).NE.0.0)QV(K)=QV(K)-CORR
   25    CONTINUE
      ENDIF
   10 CONTINUE
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C .icq4 reading and writing
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE ICQ4READ(FILICQ4,MESHICQ4,NNICQ4,NNVICQ4,
     &DELTICQ4,IDD0,IMM0,IYEAR0,SEC0,HMID,UMID,VMID,HOLD,UOLD,VOLD,
     &ZMID,ZOLD,UZMID,VZMID,WZMID,Q2MID,Q2LMID,TMID,SMID,NNDIM,NNVDIM)
C-----------------------------------------------------------------------
C Objectives: Read a hot-start .icq4 file
C             Ensure the .icq4 file was written using QUODDY4
C             Ensure array sizes are within with dimensioning parameters
C
C Complete Array Documentation: See main program QUODDY4_1.1
C
C Index Range: I=1,NNICQ4; J=1,NNVICQ4
C               NNICQ4 - number of horizontal nodes in the .icq4 file
C               NNVICQ4 - the number of vertical nodes in the .icq4 file
C
C Inputs: FILICQ4 - the name of the .icq4 hot-start file
C         NNDIM - dimensioned array size for index I
C         NNVDIM - dimensioned array size for index J
C
C Outputs: All outputs contain information related specifications from
C          the .icq4 file
C
C         MESHICQ4 - the name of the geometric mesh
C         NNICQ4 - number of horizontal nodes
C         NNVICQ4 - number of vertical nodes
C         DELTICQ4 - time step size
C         IDD0,IMM0,IYEAR0,SEC0 - timing parameters for the beginning
C                                 of the simulation
C
C       *OLD arrays are initialized with dependent variable data at one
C       time step prior to the beginning of the simulation
C
C         HOLD(I) - total depth
C         UOLD(I),VOLD(I) - x,y components of vertically averaged
C                           horizontal velocity
C         ZOLD(I,J) - 3-D nodal grid
C
C       *MID arrays are initialized with dependent variable data at the
C       beginning of the simulation  
C
C         HMID(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         Q2MID(I,J),Q2LMID(I,J) - twice the turbulent kinetic energy
C                                  and the master length scale times
C                                  twice the turbulent kinetic energy
C         TMID(I,J),SMID(I,J) - temperature and salinity
C-----------------------------------------------------------------------
C
C Global dimensioning and assignments
C
      REAL HOLD(NNDIM),UOLD(NNDIM),VOLD(NNDIM),ZOLD(NNDIM,NNVDIM)
      REAL HMID(NNDIM),UMID(NNDIM),VMID(NNDIM),ZMID(NNDIM,NNVDIM),
     &UZMID(NNDIM,NNVDIM),VZMID(NNDIM,NNVDIM),WZMID(NNDIM,NNVDIM),
     &Q2MID(NNDIM,NNVDIM),Q2LMID(NNDIM,NNVDIM),
     &TMID(NNDIM,NNVDIM),SMID(NNDIM,NNVDIM)
      CHARACTER*72 FILICQ4,MESHICQ4
C
C Local assignments
C
      CHARACTER*72 MODICQ4,COMICQ4
C
C START OF EXECUTABLE CODE
C
C Open the .icq4 file for "HOT-START" Initialization
C
      INDX=INDEX(FILICQ4,' ')-1
      OPEN(UNIT=27,FILE=FILICQ4(1:INDX),FORM='FORMATTED',STATUS='OLD')
C
C Read the model code from the first line of the .icq4 file and make
C sure it is 'QUODDY4'
C
      READ(27,'(A)')MODICQ4
      INDX=INDEX(MODICQ4,' ')-1
      IF(MODICQ4(1:INDX).NE.'QUODDY4')THEN
         WRITE(2,'(A)')'EXECUTION TERMINATED IN SUBROUTINE ICQ4READ'
         WRITE(2,'(5A)')'QUODDY4 CANNOT BE HOT-STARTED WITH A ',
     &                   MODICQ4(1:INDX),' .icq* FILE.'
         STOP
      ENDIF
C
C Read the mesh name from the second line of the .icq4 file
C
      READ(27,'(A)')MESHICQ4
C
C Read and ignore the following three comment lines from the .icq4 file
C
      DO I=1,2
         READ(27,'(A)')COMICQ4
      ENDDO
C
C Read in the number of nodes in the horizontal and vertical
C dimensions and check whether they are compatiable with dimensioning
C parameters
C
      READ(27,*)NNICQ4,NNVICQ4
      ISTOP=0
      IF(NNICQ4.GT.NNDIM) THEN
         ISTOP=1
         WRITE(2,'(A)')'NNICQ4 = ',NNICQ4
         WRITE(2,'(A)')'THE NUMBER OF HORIZONTAL NODES IN .icq4 FILE EXC
     &EEDS NNDIM'
      ENDIF
      IF(NNVICQ4.GT.NNVDIM) THEN
         ISTOP=1
         WRITE(2,'(A)')'NNVICQ4 = ',NNVICQ4
         WRITE(2,'(A)')'THE NUMBER OF VERTICAL NODES IN .icq4 FILE EXCEE
     &DS NNVDIM'
      ENDIF
      IF(ISTOP.EQ.1)THEN
         WRITE(2,'(A)')'EXECUTION TERMINATED IN SUBROUTINE ICQ4READ'
         STOP
      ENDIF
C
C Read in time step label, current simulation time, and time-step size
C
      READ(27,*)IDD0,IMM0,IYEAR0,SEC0,DELTICQ4
C
C Read in the 2-D data: HMID, UMID, VMID, HOLD, UOLD, VOLD
C
      READ(27,*)(HMID(I),UMID(I),VMID(I),HOLD(I),UOLD(I),VOLD(I),
     &I=1,NNICQ4)
C
C Read in the 3-D data: ZMID,ZOLD,UZMID,VZMID,WZMID,Q2MID,
C                              Q2LMID,TMID,SMID
C
      READ(27,*)((ZMID(I,J),ZOLD(I,J),UZMID(I,J),VZMID(I,J),WZMID(I,J),
     &Q2MID(I,J),Q2LMID(I,J),TMID(I,J),SMID(I,J),I=1,NNICQ4),
     &J=1,NNVICQ4)
C
C The .icq4 file has been successfully read
C
      INDX=INDEX(FILICQ4,' ')-1
      WRITE(2,'(/A)')'     '//FILICQ4(1:INDX)//' successfully read'
      CLOSE(27)
C
C END OF ROUTINE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE ICQ4WRITE(CASNAM,FILINQ,FILICQ,FILOUT,
     &KDMID,SECMID,DELT,HMID,UMID,VMID,HOLD,UOLD,VOLD,ZMID,ZOLD,
     &UZMID,VZMID,WZMID,Q2MID,Q2LMID,TMID,SMID,NN,NNV,NNDIM,NNVDIM)
C-----------------------------------------------------------------------
C   Description: This subroutine is used to output all the required
C                data to a ".icq4" file for the "hot start" from the
C                current status.
C     Variables: [Argument list]
C                CASNAM -- Casename (mesh name) of the simulation;
C                        there should be mesh data files associated
C                        with this casename
C                FILINQ -- Filename of maximum 45 characters
C                FILICQ -- Filename of maximum 45 characters of the
C                        previous ".icq4" file
C                FILOUT -- Output ".icq4" file for the next "hot start"
C 	         KDMID is the day number of the current date based on
C 	            the Gregorian calendar (KDMID=1 is January 1, 0000)
C                   idd0,imm0,iyear0
C                SECMID is current time on day KDMID, in seconds after 
C                   midnight.
C                HMID(*), UMID(*), VMID(*) -- Sea surface and 2-D vel.
C                        at time t
C                ZMID(NNDIM,*) -- Vertical coordinate at current time
C                UZMID(NNDIM,*), VZMID(NNDIM,*), WZMID(NNDIM,*) -- 3-D
C                        velocity solution at the current time step
C                Q2MID(NNDIM,*), Q2LMID(NNDIM,*) -- 3-D turbulence
C                        kinetic energy and its product with the
C                        mixing length at the current time step
C                        velocity solution at the current time step
C                TMID(NNDIM,*), SMID(NNDIM,*) -- 3-D temperature
C                        and salinity solution
C                DELT -- Time step size, dt in hrs.
C                HOLD(*), UOLD(*), VOLD(*) -- Sea surface and 2-D vel.
C                        at time t-dt
C                ZOLD(NNDIM,*) -- Vertical coordinate at t-dt
C                HDOWN(*) -- Bathymetric depth at nodes
C                NN -- Actual number of nodes in the horizontal mesh
C                NNDIM -- Maximum number of nodes allowed in the
C                        horizontal mesh
C                NNV -- Total number of nodes along the vertical
C                [Local integer variables]
C                I, J -- Do loop indices
C      Accuracy: I*2, R*4 (single precision)
C         Input: All the variables in the argument list
C        Output: A ".icq4" files for the next "hot start" from a
C                previous simulation
C-----------------------------------------------------------------------
      CHARACTER FILOUT*(*)
      CHARACTER*72 CASNAM
      CHARACTER*72 FILICQ,FILINQ
      REAL DELT, SECMID
      REAL HMID(*),UMID(*),VMID(*),HOLD(*),UOLD(*),VOLD(*),
     &     ZMID(NNDIM,*),ZOLD(NNDIM,*),UZMID(NNDIM,*),
     &     VZMID(NNDIM,*),WZMID(NNDIM,*),Q2MID(NNDIM,*),
     &     Q2LMID(NNDIM,*),TMID(NNDIM,*),SMID(NNDIM,*)
C
C BEGINNING OF EXECUTABLE ROUTINE
C
C Open output ".icq4" file
C
      OPEN(UNIT=20,FILE=FILOUT,FORM='FORMATTED',STATUS='UNKNOWN')
      REWIND(20)
C
C Write the first five lines of the output, ".icq4" file
C
      WRITE(20,11)'QUODDY4'
      INDX=INDEX(CASNAM,' ')-1
      WRITE(20,11)CASNAM(1:INDX)
      INDX=INDEX(FILINQ,' ')-1
      WRITE(20,11)'Input file:         '//FILINQ(1:INDX)
      INDX=INDEX(FILICQ,' ')-1
      WRITE(20,11)'Initial condition:  '//FILICQ(1:INDX)
C
C Write the number of nodes in the horizontal and in the vertical
C dimensions
C
      WRITE(20,13)NN,NNV
C
C Write time step label, current simulation time, and previous time
C
      CALL DMY(IDD,IMM,IYEAR,KDMID)
      WRITE(20,12)IDD,IMM,IYEAR,SECMID,DELT
C
C Write the 2-D data: HMID, UMID, VMID, HOLD, UOLD, VOLD
C
      WRITE(20,16)(HMID(I),UMID(I),VMID(I),HOLD(I),UOLD(I),
     &              VOLD(I),I=1,NN)
C
C Write the 3-D data: ZMID, ZOLD, UZMID, VZMID, WZMID
C
      WRITE(20,17)((ZMID(I,J),ZOLD(I,J),UZMID(I,J),VZMID(I,J),
     &              WZMID(I,J),Q2MID(I,J),Q2LMID(I,J),TMID(I,J),
     &              SMID(I,J),I=1,NN),J=1,NNV)
      CLOSE(20)
C
C Output format statements
C
 11   FORMAT(A)
 12   FORMAT(1x,2(I2,1x),i4,2(1X,f21.12))
 13   FORMAT(I8,1X,I8)
 14   FORMAT(5X,A11,I7,A4,A45)
 15   FORMAT(27X,A)
 16   FORMAT(6E17.9)
 17   FORMAT(9E17.9)
C
C END OF ROUTINE
C
      RETURN
      END
C***********************************************************************
C***********************************************************************
C
C Miscellaneous:
C
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)
            ml=inbe(3,i)
            mr=inbe(4,i)
            inbe(1,i)=inbe(1,1)
            inbe(2,i)=inbe(2,1)
            inbe(3,i)=inbe(3,1)
            inbe(4,i)=inbe(4,1)
            inbe(1,1)=nstart
            inbe(2,1)=nhead
            inbe(3,1)=ml
            inbe(4,1)=mr
            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)
               ml=inbe(3,i)
               mr=inbe(4,i)
               inbe(1,i)=inbe(1,j)
               inbe(2,i)=inbe(2,j)
               inbe(3,i)=inbe(3,j)
               inbe(4,i)=inbe(4,j)
               inbe(1,j)=n1
               inbe(2,j)=n2
               inbe(3,j)=ml
               inbe(4,j)=mr
               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 BASIS1D1(nn,zref,zdes,b2,n2,b1,n1)
C-----------------------------------------------------------------------
C purpose: This subroutine evaluates the value of the basis functions 
C            alive on a 1-D linear element at a point
C
C restrictions: Applicable only for  1-D linear elements
C
C inputs:  nn is the number of entries in the reference z array
C          zref(nn) is the 1-D array containing the reference z array
C               NOTE: zref must increase from zref(1) to zref(nn)
C          zdes is the z value at which the basis functions are desired
C
C outputs: n1 is the reference node immediately below zdes
C          b1 is the value of the n1's basis function at zres
C          n2 is the reference node immediately above zdes
C          b2 is the value of the n2's basis function at zres
C
C notes:   zdes<zref(1)  => n1=nn,b1=0.0,n2=1,b2=1.0 
C          zdes>zref(nn) => n1=nn,b1=1.0,n2=1,b2=0.0 
C
C history:  Written by Christopher E. Naimie
C           Dartmouth College
C           26 AUGUST 1992
C-----------------------------------------------------------------------
C
C ARGUMENTS
      dimension zref(nn)
C
C START OF EXECUTABLE CODE
      do 10 k=1,nn 
         if(zdes.LT.zref(k))then
C
C if zdes<zref(1), zdes is below the reference data => set b1=0,n1=nn
            n2=k
            n1=k-1
            if(n2.EQ.1)then
                n1=nn
                b1=0.0
                b2=1.0
                go to 20
             endif
             zdiff=zref(n2)-zref(n1)
             b1=(zref(n2)-zdes)/zdiff
             b2=(zdes-zref(n1))/zdiff  
             go to 20
         endif
 10   continue
C
C if zdes>zref(nn), zdes is above the reference data => set b2=n2=0
      b1=1.0
      n1=nn
      b2=0.0
      n2=1
 20   continue
C
C END OF ROUTINE
      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 ELEMCOEFS(I,NNDIM,NNV,ENZMI,ENZHI,ENZQI,
     &                          ENZM,ENZH,ENZQ)
C-----------------------------------------------------------------------
C...Store the elemental based values of the diffusion coefficients
C   into 3-D arrays
C-----------------------------------------------------------------------
      REAL ENZMI(*),ENZHI(*),ENZQI(*),
     &     ENZM(NNDIM,*),ENZH(NNDIM,*),ENZQ(NNDIM,*)
      NNV1=NNV-1
      ENZM(I,1)=ENZMI(1)
      ENZH(I,1)=ENZHI(1)
      ENZQ(I,1)=ENZQI(1)
      DO J=2,NNV1
         ENZM(I,J)=ENZMI(J)
         ENZH(I,J)=ENZHI(J)
         ENZQ(I,J)=ENZQI(J)
      ENDDO
      ENZM(I,NNV)=0.
      ENZH(I,NNV)=0.
      ENZQ(I,NNV)=0.
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE NODECOEFS(I,NNDIM,NNV,ZMIDI,ENZMI,ENZHI,ENZQI,
     &                          ENZM,ENZH,ENZQ)
C-----------------------------------------------------------------------
C...Convert the elemental based values of Nz to nodal based using
C   Nz(i) = (dz(I)*Nz(I) + dz(II)*Nz(II))/(dz(I) + dz(II))
C-----------------------------------------------------------------------
      REAL ENZMI(*),ENZHI(*),ENZQI(*),ZMIDI(*),
     &    ENZM(NNDIM,*),ENZH(NNDIM,*),ENZQ(NNDIM,*)
      NNV1=NNV-1
      ENZM(I,1)=ENZMI(1)
      ENZH(I,1)=ENZHI(1)
      ENZQ(I,1)=ENZQI(1)
      DO J=2,NNV1
         DZ1=ZMIDI(J)-ZMIDI(J-1)
         DZ2=ZMIDI(J+1)-ZMIDI(J)
         DZTMP=DZ1+DZ2
         ENZM(I,J)=(DZ1*ENZMI(J-1)+DZ2*ENZMI(J))/DZTMP
         ENZH(I,J)=(DZ1*ENZHI(J-1)+DZ2*ENZHI(J))/DZTMP
         ENZQ(I,J)=(DZ1*ENZQI(J-1)+DZ2*ENZQI(J))/DZTMP
      ENDDO
      ENZM(I,NNV)=ENZMI(NNV1)
      ENZH(I,NNV)=ENZHI(NNV1)
      ENZQ(I,NNV)=ENZQI(NNV1)
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTAVG(F,FAVG,Z,NN,NNV,NNDIM)
      REAL F(NNDIM,*), FAVG(*), Z(NNDIM,*)
C-----------------------------------------------------------------------
C COMPUTES THE VERTICAL AVERAGE OF F, RETURNS IT IN FAVG.
C DOES THIS FOR ALL HORIZONTAL NODES
C DOES NOT REQUIRE UNIFORM MESH
C INTEGRATION IS PERFORMED FROM NODE 1 TO NODE NNV
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF F AND Z IN MAIN.
C-----------------------------------------------------------------------
      DO I=1,NN
         FINT=0.0
         DO J=2,NNV
            FINT=FINT+0.5*(F(I,J)+F(I,J-1))*(Z(I,J)-Z(I,J-1))
         ENDDO
         FAVG(I)=FINT/(Z(I,NNV)-Z(I,1))
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTINTI(F,FINT,Z,NNV)
      REAL F(*),Z(*)
C-----------------------------------------------------------------------
C COMPUTES THE VERTICAL INTEGRAL OF F, RETURNS IT IN FINT.
C DOES NOT REQUIRE UNIFORM MESH
C INTEGRATION IS PERFORMED FROM NODE 1 TO NODE NNV
C-----------------------------------------------------------------------
      FINT=0.0
      DO J=2,NNV
         FINT=FINT+0.5*(F(J)+F(J-1))*(Z(J)-Z(J-1))
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE VERTSUM(F,FSUM,NN,NEV,NNDIM)
C-----------------------------------------------------------------------
C COMPUTES THE VERTICAL SUM OF F, RETURNS IT IN FSUM.
C DOES THIS FOR ALL HORIZONTAL NODES
C SUMMATION IS PERFORMED FROM ELEMENT 1 TO ELEMENT NEV
C NNDIM IS THE FIRST PHYSICAL DIMENSION OF F IN MAIN.
C-----------------------------------------------------------------------
      REAL F(NNDIM,*), FSUM(*)
      DO I=1,NN
         FS=0.0
         DO J=1,NEV
            FS=FS+F(I,J)
         ENDDO
         FSUM(I)=FS
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE GDAY(IDD,IMM,IYear,KD)
C!
C!  GIVEN DAY,MONTH,(EACH 2 DIGITS) and year (four digits), GDAY RETURNS
C!  THE DAY#, KD BASED ON THE GREGORIAN CALENDAR.
C!  THE GREGORIAN CALENDAR, CURRENTLY 'UNIVERSALLY' IN USE WAS
C!  INITIATED IN EUROPE IN THE SIXTEENTH CENTURY. NOTE THAT GDAY
C!  IS VALID ONLY FOR GREGORIAN CALENDAR DATES.
C
C   KD=1 CORRESPONDS TO JANUARY 1, 0000
C	
c 	Note that the Gregorian reform of the Julian calendar 
c	omitted 10 days in 1582 in order to restore the date
c	of the vernal equinox to March 21 (the day after
c	Oct 4, 1582 became Oct 15, 1582), and revised the leap 
c	year rule so that centurial years not divisible by 400
c	were not leap years.
c
C   THIS ROUTINE WAS WRITTEN BY EUGENE NEUFELD, AT IOS, IN JUNE 1990.
C
      INTEGER NDP(13)
      INTEGER NDM(12)
	
      DATA NDP/0,31,59,90,120,151,181,212,243,273,304,334,365/
      DATA NDM/31,28,31,30,31,30,31,31,30,31,30,31/
C!
	lp=6
c make iyy and icc variables
        icc=iyear/100
        iyy=iyear-icc*100

C!  TEST FOR INVALID INPUT:
      IF(ICC.LT.0)THEN
	 WRITE(LP,5000)ICC
	 STOP
      ENDIF
      IF(IYY.LT.0.OR.IYY.GT.99)THEN
	 WRITE(LP,5010)IYY
	 STOP
      ENDIF
      IF(IMM.LE.0.OR.IMM.GT.12)THEN
	 WRITE(LP,5020)IMM
	 STOP
      ENDIF
      IF(IDD.LE.0)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.NE.2.AND.IDD.GT.NDM(IMM))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.29)THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
      IF(IMM.EQ.2.AND.IDD.GT.28.AND.((IYY/4)*4-IYY.NE.0.OR.(IYY.EQ.0.AND
     .    .(ICC/4)*4-ICC.NE.0)))THEN
	 WRITE(LP,5030)IDD
	 STOP
      ENDIF
5000  FORMAT(' INPUT ERROR. ICC = ',I7)
5010  FORMAT(' INPUT ERROR. IYY = ',I7)
5020  FORMAT(' INPUT ERROR. IMM = ',I7)
5030  FORMAT(' INPUT ERROR. IDD = ',I7)
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST CENTURY:
      KD = ICC*36524 + (ICC+3)/4
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST YEAR:
      KD = KD + IYY*365 + (IYY+3)/4
C!
C!  ADJUST FOR CENTURY RULE:
C!  (VIZ. NO LEAP-YEARS ON CENTURYS EXCEPT WHEN THE 2-DIGIT
C!  CENTURY IS DIVISIBLE BY 4.)
      IF(IYY.GT.0.AND.(ICC-(ICC/4)*4).NE.0) KD=KD-1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF LAST YEAR.
C!
C!  CALCULATE DAY# OF LAST DAY OF LAST MONTH:
      KD = KD + NDP(IMM)
C!
C!  ADJUST FOR LEAP YEARS:
      IF(IMM.GT.2.AND.((IYY/4)*4-IYY).EQ.0.AND.((IYY.NE.0).OR.
     .   (((ICC/4)*4-ICC).EQ.0)))   KD=KD+1
C!  KD NOW TRULY REPRESENTS THE DAY# OF THE LAST DAY OF THE LAST
C!  MONTH.
C!
C!  CALCULATE THE CURRENT DAY#:
      KD = KD + IDD
      RETURN
C!
C!
      ENTRY DMY(IDD,IMM,IYear,KD)
C!
C!  GIVEN THE (GREGORIAN) DAY#, KD, AS CALCULATED ABOVE IN THIS ROUTINE,
C!  ENTRY DMY RETURNS THE (GREGORIAN) DAY, MONTH, YEAR AND CENTURY.
C!
C!  TEST FOR VALID INPUT:
      IF(KD.LE.0) WRITE(LP,5040)KD
5040  FORMAT(' KD = ',I7,'  INVALID INPUT. DMY STOP.')
C!
C!  SAVE KD
      KKD=KD
C!  CALCULATE ICC AND SUBTRACT THE NUMBER OF DAYS REPRESENTED BY ICC
C!  FROM KKD
C!  JFH IS THE NUMBER OF 400 YEAR INTERVALS UP TO KKD
C!  JCC IS THE NUMBER OF ADDITIONAL CENTURIES UP TO KKD
      JFH = KKD/146097
      KKD = KKD - JFH*146097
      IF(KKD.LT.36525)THEN
	 JCC = 0
      ELSE
	 KKD = KKD - 36525
	 JCC = 1 + KKD/36524
	 KKD = KKD - (JCC-1)*36524
      END IF
      ICC = 4*JFH + JCC
      IF(KKD.EQ.0)THEN
	 ICC = ICC-1
	 IYY = 99
	 IMM = 12
	 IDD = 31
c	 RETURN
	go to 110
      ENDIF
C!
C!  CALCULATE IYY. JFY IS THE NUMBER OF FOUR YEAR INTERVALS IN THE
C!  CURRENT CENTURY. THE FIRST FOUR YEAR INTERVAL IS SHORT (1460 DAYS
C!  RATHER THAN 1461)IF THE CURRENT CENTURY IS NOT DIVISIBLE BY 4, AND
C!  IN THIS CASE JCC.NE.0 AS CALCULATED ABOVE.
C!
C!  CALCULATE JFY:
      JFY = 0
      IF(JCC.EQ.0)GOTO 10
      IF(KKD.LT.1460)GOTO 10
      JFY = 1
      KKD = KKD - 1460
10    KK = KKD/1461
      JFY = JFY + KK
      KKD = KKD - KK*1461
C!
C!  CALCULATE JYY, THE REMAINING YEARS OF THE CURRENT CENTURY UP TO THE
C!  CURRENT DAY:
      JYY = 0
C!  THE NEXT YEAR IS NOT A LEAP YEAR IF JFY=0 AND JCC.NE.0.
      IF(JFY.EQ.0.AND.JCC.NE.0)GOTO 20
      IF(KKD.LT.366)GOTO 30
      JYY = 1
      KKD = KKD - 366
20    JYYY = KKD/365
      JYY = JYY + JYYY
      KKD = KKD - JYYY*365
30    IYY = 4*JFY + JYY
      IF(KKD.EQ.0) THEN
	 IYY=IYY-1
	 IMM=12
	 IDD=31
c	 RETURN
	go to 110
      END IF
C!
C!  SET L=1 IF WE HAVE A LEAP YEAR.
      L=0
      IF(IYY-(IYY/4)*4.NE.0)GOTO 40
      IF(IYY.EQ.0.AND.(ICC-(ICC/4)*4).NE.0)GOTO 40
      L=1
C!
C!  CALCULATE IMM AND IDD
40    IF(KKD.GT.31) GOTO 50
      IMM=1
      IDD=KKD
c      RETURN
       go to 110
C!
50    IF(KKD.GT.59)GOTO 60
      IMM = 2
      IDD = KKD-31
c      RETURN
	go to 110
C!
60    IF(KKD.GT.60)GOTO 70
      IF(L.EQ.0)GOTO 70
      IMM = 2
      IDD = 29
c      RETURN
	go to 110
C!
70    IF(L.EQ.1) KKD=KKD-1
      DO 80 I=4,13
	 IF(KKD.GT.NDP(I))GOTO 80
	 IMM = I-1
	 IDD = KKD - NDP(I-1)
c	 RETURN
	go to 110
C!
80    CONTINUE

 110 	iyear=icc*100+iyy
	return

90    WRITE(LP,5050)
5050  FORMAT(' ERROR IN DMY.')
      STOP
      END

	Subroutine Up_date(kd,ssec)
C 
C this subroutine updates the current gregorian day number kd
C and resets the time, ssec, accordingly
C it inherently assumes that the variable is updated every day!

      integer kd
      real ssec
	
      kd=kd+aint(ssec/8.64E4)
      ssec=amod(ssec,8.64E4)
 
      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***********************************************************************
C
C   Symmetric banded matrix setup:
C  
C***********************************************************************
C***********************************************************************
      SUBROUTINE RHSBUILD(A,AA,IDBC,NN,NB,NDBC,NB2DIM,NB1DIM)
C-----------------------------------------------------------------------
C DRL, 9 March 1994
C
C Condensation of a symmetric banded matrix in BANSOLtr (transposed addressing) 
C   form for symmetry-preserving Dirichlet BC enforcement 
C
C Global Columns [IDBC(I), I=1,NDBC] of A are copied into AA and replaced with zeros in A.
C Global Rows    [IDBC(I), I=1,NDBC] of A are set to zero with diagonal set to unity.
C Corresponding rows in AA are set to zero.
C 
C A is a matrix stored in transpose, column-banded, symmetric form
C     A(column, row) is A(nb+1, nn) (BANSOLtr form)
C     A(1,*) is the diagonal
C AA is a matrix in column-banded, column-condensed form:
C     A(row, column) is A((2*nb+1), nx)
C     Note this is conventional, not transposed, addressing;
C          and AA is not symmetric
C
C Column condensation of AA:
C   IDBC(k) is a list of global node numbers for Dirichlet BC's
C   IDBC(k) is the global column number of column k in AA
C       AA(*,k) multiplies the BC data for the kth Dirichlet BC
C   NDBC is the number of Dirichlet BC's (number of columns in AA)
C
C Column-banding of AA:
C   Columns stored have total band height 2NB+1
C   The column is stored beginning in row 1
C   The main diagonal is at row NB+1
C
C NB is the half-bandheight (Full BW=2*NB+1)
C NDBC is the number of Dirichlet BC's
C            (number of columns to be copied into AA)
C NN number of nodes globally
C NB1DIM is the first physical dimension of A in the main program
C NB2DIM is the first physical dimension of AA in the main program
C    --- NB1DIM.ge.(  NB+1) ---
C    --- NB2DIM.ge.(2*NB+1) ---
C-----------------------------------------------------------------------
      DIMENSION A(NB1DIM,*),AA(NB2DIM,*),IDBC(*)
C
C Check dimensions passed
C
      NB2=2*NB+1
      IF (NB2DIM.LT.NB2) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX AA IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB2DIM =', NB2DIM, ' NB2 =', NB2
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
      NB1=NB+1
      IF (NB1DIM.LT.NB1) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX A IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB1DIM =', NB1DIM, ' NB1 =', NB1
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
C Zero the AA matrix
C
      DO 5 I=1,NB2
      DO 5 J=1,NDBC
    5 AA(I,J)=0.0
C
C Begin the copying/zeroing
C
      DO 30 I=1,NDBC
      II=IDBC(I)
C
C On and Above the Diagonal:
C
      DO 10 L=0,MIN0(NB,II-1)
   10 AA(NB+1-L,I)=A(L+1,II-L)
C 
C Below the Diagonal:(exploiting symmetry, use row of A)
C
      DO 20 J=2,NB1
   20 AA(NB+J,I)=A(J,II)
C
   30 CONTINUE
C
C Now the AA matrix is complete.  
C Note that AA transpose is the matrix of Galerkin equations for 
C  Dirichlet nodes.
C Begin zeroing the A matrix
C
      DO 60 I=1,NDBC
      II=IDBC(I)
C
C On and Above the Diagonal:
C  
      A(1,II)=1.0
      IF (II .GT. 1) THEN
      DO 40 L=1,MIN0(NB,II-1)
   40 A(L+1,II-L)=0.
      ENDIF
C 
C Below the Diagonal:(exploiting symmetry)
C
      DO 50 J=2,NB1
   50 A(J,II)=0.
C
   60 CONTINUE
C
C The A matrix is now ready with rows and columns zeroed and unity on 
C  diagonal.
C Now zero the Dirichlet rows of AA
C Note AA transpose will no longer hold the Dirichlet Galerkin eqns.
C
C V.N BC Modification-keep dirichlet galerkin equations(CEN&MJH 4/30/97)
*      DO 90 J=1,NDBC
*      JJ=IDBC(J)
*      DO 90 I=1,NDBC
*      II=IDBC(I)
*      IA=(NB+1)+(II-JJ)
*      IF(IA.GE.1 .AND. IA.LE.NB2) AA(IA,J)=0.0
*   90 CONTINUE
C End of V.N BC Modification
C
C Finished
C
      RETURN
      END
C
C In Quoddy:
C      
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***********************************************************************
      SUBROUTINE RHSMULT(A,X,R,JX,NN,NB,NX,NB2DIM)
C-----------------------------------------------------------------------
C DRL, 9 March 1994
C
C Performs the matrix algebra r=r-A*x for known r, A, and x
C A and x are undisturbed; r is overwritten as above
C
C r is a n-length vector
C A is a matrix in column-banded, column-condensed form:
C     A(row, column) is A((2*nb+1), nx)
C x is a n-length vector
C
C Column condensation:
C   JX(k) is a pointer from local to global node number
C   JX(k) is the global column number of column k in A
C       A(*,j) multiplies x(jx(j))
C
C Column-banding:
C   Columns stored have total band height 2NB+1
C   The column is stored beginning in row 1
C   The main diagonal is at row NB+1
C
C NB is the half-bandheight (Full BW=2*NB+1)
C NX is the number of columns in A (Number of Type I BC's)
C NN is the length of the r and x vectors (Number of nodes globally)
C NB2DIM is the first physical dimension of A in the main program
C    --- NB2DIM.ge.(2*NB+1) ---
C-----------------------------------------------------------------------
      DIMENSION A(NB2DIM,*), R(*), X(*),JX(*)
C
C Check dimensions passed
C
      NB2=2*NB+1
      IF (NB2DIM.LT.NB2) THEN
         WRITE(*,*) ' INSUFFICIENT DIMENSION FOR MATRIX A IN'//
     &              ' SUBROUTINE RHSMULT'
         WRITE(*,*) ' NB2DIM =', NB2DIM, ' NB2 =', NB2
         WRITE(*,*) ' EXECUTION TERMINATED'
         STOP
      ENDIF
C
      DO 10 J=1,NX
C
C IR is the global row number (Position in r)
C I is the local row number in A (column-banded)
C J  is the local column number in A (columns are condensed)
C JX is the global column number (Position in x)
C Row correspondence:  I = (NB+1) + (IR-JX(J))    
C
      IR1 = MAX0(1,JX(J)-NB)
      IR2 = MIN0(NN,JX(J)+NB)
      I  = NB+1 + (IR1-JX(J))
C
      DO 20 IR=IR1,IR2
      R(IR)=R(IR)-A(I,J)*X( JX(J) )
      I=I+1
   20 CONTINUE
C
   10 CONTINUE
      RETURN
      END
C
C In Quoddy:
C      
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***********************************************************************
C
C...Banded matrix solvers:
C
C***********************************************************************
C***********************************************************************
      SUBROUTINE THOMAS(KKK,A,B,C,R,N)            
C-----------------------------------------------------------------------
C THOMAS ALGORITHM FOR TRIDIAGONAL MATRICES  -  SYSTEM FORTRAN-77
C REAL VERSION
C A, B, AND C = LHS COLUMNS, DIMENSIONED (=>N) IN CALLING PROGRAM
C R = RIGHT-HAND-SIDE, DIMENSIONED (=>N) IN CALLING PROGRAM 
C N = # OF EQUATIONS             
C SOLUTION RETURNS IN R WHEN KKK => 2
C KKK = 1 PERFORMS LU DECOMPOSITION DESTRUCTIVELY            
C KKK = 2 PERFORMS BACK SUBSTITUTION
C KKK = 3 PERFORMS OPTIONS 1 AND 2
C-----------------------------------------------------------------------
      DIMENSION A(N),B(N),C(N),R(N)
C                                                                               C  STEP 1 - TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      IF(KKK.EQ.2)GO TO 45
      C(1)=C(1)/B(1)           
      DO 40 I=2,N              
      IM=I-1            
      C(I)=C(I)/(B(I)-A(I)*C(IM))            
   40 CONTINUE
      IF(KKK.EQ.1)RETURN
C                
C STEP 2 - MODIFY LOAD VECTOR R 
C
   45 CONTINUE
      R(1)=R(1)/B(1)           
      DO 41 I=2,N              
      IM=I-1            
      R(I)=(R(I)-A(I)*R(IM))/(B(I)-A(I)*C(IM))      
   41 CONTINUE          
C
C BACK SUBSTITUTE SOLUTION INTO R VECTOR
C                
      DO 50 I=1,N-1            
      J=N-I             
      JP=J+1            
      R(J)=R(J)-C(J)*R(JP)            
   50 CONTINUE          
      RETURN            
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE CTHOMAS(KKK,A,B,C,R,N)            
C-----------------------------------------------------------------------
C THOMAS ALGORITHM FOR TRIDIAGONAL MATRICES  -  SYSTEM FORTRAN-77
C COMPLEX VERSION
C A, B, AND C = LHS COLUMNS, DIMENSIONED (=>N) IN CALLING PROGRAM
C R = RIGHT-HAND-SIDE, DIMENSIONED (=>N) IN CALLING PROGRAM 
C N = # OF EQUATIONS             
C SOLUTION RETURNS IN R WHEN KKK => 2
C KKK = 1 PERFORMS LU DECOMPOSITION DESTRUCTIVELY            
C KKK = 2 PERFORMS BACK SUBSTITUTION
C KKK = 3 PERFORMS OPTIONS 1 AND 2
C-----------------------------------------------------------------------
      COMPLEX A(N),B(N),C(N),R(N)          
C                                                                               
C  STEP 1 - TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      IF(KKK.EQ.2)GO TO 45
      C(1)=C(1)/B(1)           
      DO 40 I=2,N              
      IM=I-1            
      C(I)=C(I)/(B(I)-A(I)*C(IM))            
   40 CONTINUE
      IF(KKK.EQ.1)RETURN
C                
C STEP 2 - MODIFY LOAD VECTOR R 
C
   45 CONTINUE
      R(1)=R(1)/B(1)           
      DO 41 I=2,N              
      IM=I-1            
      R(I)=(R(I)-A(I)*R(IM))/(B(I)-A(I)*C(IM))      
   41 CONTINUE          
C
C BACK SUBSTITUTE SOLUTION INTO R VECTOR
C                
      DO 50 I=1,N-1            
      J=N-I             
      JP=J+1            
      R(J)=R(J)-C(J)*R(JP)            
   50 CONTINUE          
      RETURN            
      END  
C***********************************************************************
C***********************************************************************
      SUBROUTINE SOLVETR(KKK,B,R,NEQ,IHALFB,NDIM,MDIM)        
C-----------------------------------------------------------------------
C   16 FEB 1992 -- SOLVES FOR B TRANSPOSED -- (COLUMN, ROW)
C                
C   ASYMMETRIC BAND MATRIX EQUATION SOLVER           
C   DOCTORED TO IGNORE ZEROS IN LU DECOMP. STEP      
C                
C   KKK=1 TRIANGULARIZES THE BAND MATRIX B           
C   KKK=2 SOLVES FOR RIGHT SIDE R, SOLUTION RETURNS IN R      
C-----------------------------------------------------------------------                
      DIMENSION B(MDIM,ndim),R(NDIM)        
      NRS=NEQ-1           
      IHBP=IHALFB+1       
      IF (KKK.EQ.2) GO TO 30       
C                
C  TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD     
C                
      DO 10 K=1,NRS       
      PIVOT=B(IHBP,k)              
      KK=K+1              
      KC=IHBP             
      DO 21 I=KK,NEQ      
      KC=KC-1             
      IF(KC.LE.0) GO TO 10         
      C=-B(KC,i)/PIVOT             
      IF (C.EQ.0.0) GO TO 21       
      B(KC,i)=C           
      KI=KC+1             
      LIM=KC+IHALFB       
      DO 20 J=KI,LIM      
      JC=IHBP+J-KC        
   20 B(J,i)=B(J,i)+C*B(JC,k)      
   21 CONTINUE            
   10 CONTINUE            
      GO TO 100           
C                
C  MODIFY LOAD VECTOR R            
C                
   30 NN=NEQ+1            
      IBAND=2*IHALFB+1             
      DO 40 I=2,NEQ       
      JC=IHBP-I+1         
      JI=1       
      IF (JC.LE.0) GO TO 50        
      GO TO 60            
   50 JC=1       
      JI=I-IHBP+1         
   60 SUM=0.0             
      DO 70 J=JC,IHALFB            
      SUM=SUM+B(J,i)*R(JI)         
   70 JI=JI+1             
   40 R(I)=R(I)+SUM       
C                
C   BACK SOLUTION         
C                
      R(NEQ)=R(NEQ)/B(IHBP,neq)             
      DO 80 IBACK=2,NEQ            
      I=NN-IBACK          
      JP=I       
      KR=IHBP+1           
      MR=MIN0(IBAND,IHALFB+IBACK)           
      SUM=0.0             
      DO 90 J=KR,MR       
      JP=JP+1             
   90 SUM=SUM+B(J,i)*R(JP)         
   80 R(I)=(R(I)-SUM)/B(IHBP,i)             
  100 RETURN              
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE BANSOLtr(KKK,B,R,NEQ,IBAND,NDIM,MDIM)
C-----------------------------------------------------------------------
C SYSTEM FORTRAN-77
C
C     TRANSPOSED VERSION of BANSOL (DRL 8 MARCH 94):
C            B(column, row)
C            Sequence of calling arguments unchanged from BANSOL
C
C B = LHS MATRIX , DIMENSIONED (mDIM,nDIM) IN CALLING PROGRAM
C R = RIGHT-HAND-SIDE, DIMENSIONED (NDIM) IN CALLING PROGRAM
C NEQ= # OF EQUATIONS
C IBAND = DIAGONAL(COLUMN 1) + HALF-BANDWIDTH
C SOLUTION RETURNS IN R
C KKK = 1 PERFORMS LU DECOMPOSITION DESTRUCTIVELY
C KKK = 2 PERFORMS BACK SUBSTITUTION
C KKK = 3 PERFORMS OPTIONS 1 AND 2
C
C     SYMMETRIC BAND MATRIX EQUATION SOLVER,   (REF, 2)
C     DOCTORED TO IGNORE ZEROS IN LU DECOMP. STEP (DRL 8 MARCH 94)
C     TRANSPOSED VERSION: B(column, row) (DRL 8 MARCH 94)
C-----------------------------------------------------------------------
      DIMENSION B(mDIM,nDIM), R(NDIM)
      NRS = NEQ - 1
      NR = NEQ
      IF(KKK .EQ. 2) GO TO 200
C
C  TRIANGULARIZE MATRIX A USING DOOLITTLE METHOD
C
      DO 120 N = 1,NRS
      M= N-1
      MR=MIN0(IBAND,NR-M)
      PIVOT = B(1,N)
      DO  120 L=2,MR
      CP= B(l,N)/PIVOT
      IF (CP.EQ.0.0) GO TO 120
      I = M+L
      J = 0
      DO 110 K=L,MR
      J = J + 1
  110 B(j,I)= B(j,I) -CP*B(k,N)
  120 B(l,N) = CP
      IF(KKK.EQ.1)GO TO 400
C
C  MODIFY LOAD VECTOR R
C
  200 DO 220 N=1,NRS
      M= N-1
      MR = MIN0(IBAND,NR-M)
      CP= R(N)
      R(N)=CP/B(1,N)
      DO 220  L=2,MR
      I = M + L
  220 R(I)=R(I)-B(l,N)*CP
C
C   BACK SOLUTION
C
      R(NR) = R(NR)/B(1,NR)
      DO 320 I = 1,NRS
      N=NR-I
      M= N-1
      MR=MIN0(IBAND,NR-M)
      DO 320 K = 2,MR
      L = M+K
C     STORE COMPUTED DISPLACEMENTS IN LOAD VECTOR R
  320 R(N)= R(N)- B(k,N)*R(L)
  400 RETURN
      END
C***********************************************************************
C***********************************************************************
C End of Fixed Subroutines for QUODDY4 (RELEASE 2.1)
C End of Fixed Subroutines for QUODDY4 (RELEASE 2.1)
C***********************************************************************
C***********************************************************************
