C***********************************************************************
C***********************************************************************
C EQUATION OF STATE SUBROUTINES
C=======================================================================
      SUBROUTINE EQSTATE1_2D(T0,NN,NNV,T,RHO)
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
C history: SUBROUTINE EQSTATE1 - Released with Q4_2.1
C          SUBROUTINE EQSTATE1_2D - Revised to loop over 2-D arrays
C            internally and perform calculations without intermediate
C            sigmat step
C-----------------------------------------------------------------------
      INCLUDE 'DCMS.DIM'
      REAL T(NNDIM,NNVDIM),RHO(NNDIM,NNVDIM)
C
C Calculate reference density
      REFDENS=1028.0*(1.0-1.7E-4*(T0-10.0))
C
C Calculate density and RHO for all locations in 2-D array
      DO J=1,NNV
         DO I=1,NN
            DENSITY=1028.0*(1.0-1.7E-4*(T(I,J)-10.0))
            RHO(I,J)=(DENSITY-REFDENS)/REFDENS
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
      SUBROUTINE EQSTATE2_2D(T0,S0,NN,NNV,T,S,RHO)
C-----------------------------------------------------------------------
C purpose: This subroutine computes the perturbation density (R; in 
C          sigmat units) as a function of two variables; temperature 
C          (T; -10 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 Nested by B. Blanton
C Changed IF statements.  Jennifer Shore. December 1999.
C  Note "freezing" in Gill -> likely for salt-saturated state = 0 Fahr
C-----------------------------------------------------------------------
C
C Global Arrays
      INCLUDE 'DCMS.DIM'
      REAL T(NNDIM,NNVDIM),S(NNDIM,NNVDIM),RHO(NNDIM,NNVDIM)
C
C Equation of state coefficients
      PARAMETER (A0=-0.157406)
      PARAMETER (A1=6.793952E-2)
      PARAMETER (A2=-9.095290E-3)
      PARAMETER (A3=1.001685E-4)
      PARAMETER (A4=-1.120083E-6)
      PARAMETER (A5=6.536332E-9)
      PARAMETER (B0=0.824493)
      PARAMETER (B1=-4.0899E-3)
      PARAMETER (B2=7.6438E-5)
      PARAMETER (B3=-8.2467E-7)
      PARAMETER (B4=5.3875E-9)
      PARAMETER (C0=-5.72466E-3)
      PARAMETER (C1=1.0227E-4)
      PARAMETER (C2=-1.6546E-6)
      PARAMETER (D0=4.8314E-4)
C
C Ensure T and S are within Range, terminate execution if they are not
      NERR=0
      DO J=1,NNV
         DO I=1,NN
            TEMP=T(I,J)
            SALI=S(I,J)
            IF((TEMP.GE.-10.0).AND.(TEMP.LE.40))GO TO 10
            WRITE(2,'(A)')'T LIES OUTSIDE RANGE (-10 to 40 degrees C).'
            WRITE(2,*)'I,J,T(I,J) = ',I,J,T(I,J)
            NERR=1
 10         IF((SALI.GE.0.0).AND.(SALI.LE.42))GO TO 20
            WRITE(2,'(A)')'S LIES OUTSIDE RANGE (  0 to 42 psu).'
            WRITE(2,*)'I,J,S(I,J) = ',I,J,S(I,J)
            NERR=1
 20         CONTINUE
         ENDDO
      ENDDO
      IF(NERR.EQ.1)THEN
         WRITE(2,'(///A)')'EXECUTATION TERMINATED IN EQSTATE2_2D.'
c         STOP
      ENDIF
C
C Calculate reference density:
C   First, calculate sigmat of pure water
C   Then, adjust to one standard atmosphere
      TEMP=T0
      SALI=S0
      REFSIGT=A0+TEMP*(A1+TEMP*(A2+TEMP*(A3+TEMP*(A4+TEMP*A5))))
      REFSIGT=REFSIGT+SALI*(B0+TEMP*(B1+TEMP*(B2+TEMP*(B3+TEMP*B4))))
     &            +SALI**1.5*(C0+TEMP*(C1+TEMP*C2))+D0*SALI*SALI
      REFDENS=1000.0+REFSIGT
C
C Calculate density and RHO for all locations in 2-D array
C   First, calculate sigmat of pure water
C   Then, adjust to one standard atmosphere
C   Finally, convert to RHO
      DO J=1,NNV
         DO I=1,NN
            TEMP=T(I,J)
            SALI=S(I,J)
            SIGT=A0+TEMP*(A1+TEMP*(A2+TEMP*(A3+TEMP*(A4+TEMP*A5))))
            SIGT=SIGT+SALI*(B0+TEMP*(B1+TEMP*(B2+TEMP*(B3+TEMP*B4))))
     &            +SALI**1.5*(C0+TEMP*(C1+TEMP*C2))+D0*SALI*SALI
            RHO(I,J)=(SIGT-REFSIGT)/(REFDENS)
         ENDDO
      ENDDO
      RETURN
      END
C***********************************************************************
C***********************************************************************
C BAROCLINIC PRESSURE GRADIENT SUBROUTINES
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 'DCMS.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 'DCMS.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 'DCMS.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 = MIN(J+1,NNV)
            JM = MAX(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=MIN(J+1,NNV)
            JM=MAX(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,cs)
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 'DCMS.DIM'
C
C Global arrays
C  fem mesh
      integer in(3,nedim)
      real x(nndim),y(nndim),cs(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,CS)
      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)
         csl=(cs(i1)+cs(i2)+cs(i3))/3.0
         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=MAX(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/csl
            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
C       (base calculation on relative measure which will satisfy single
C       precision comparisions)
         deptol=0.00001
         if(abs(zdepthmin/z1-1.0).LE.deptol)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=MAX(-bathy(n1),-bathy(n2),-bathy(n3))
         DNEIGH(n1)=MAX(DNEIGH(n1),botmin)
         DNEIGH(n2)=MAX(DNEIGH(n2),botmin)
         DNEIGH(n3)=MAX(DNEIGH(n3),botmin)
      enddo
      return
      end
C***********************************************************************
C***********************************************************************
      SUBROUTINE BUILDPPX1PPY1(NN,NE,X,Y,IN,BATHY,
     &IQ,JQ,PPX1,PPY1,CS)
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 'DCMS.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),CS(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=MIN(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*CS(II)/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)*CS(II)/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***********************************************************************
