C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.21 2014/04/04 20:56:32 jmc Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#define USE_FACTORIZED_POLY

C--  File find_alpha.F: Calculate expansion Coeff
C--   Contents
C--   o FIND_ALPHA
C--   o FIND_BETA

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: FIND_ALPHA
C     !INTERFACE:
      SUBROUTINE FIND_ALPHA (
     I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
     O     alphaLoc,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | o SUBROUTINE FIND_ALPHA
C     |   Calculates [drho(S,T,z) / dT] of a horizontal slice
C     *==========================================================*
C     | k        - is the Theta/Salt level
C     | kRef     - determines pressure reference level
C     |            (not used in 'LINEAR' mode)
C     | alphaLoc - drho / dT (kg/m^3/C)
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     k      :: Level of Theta/Salt slice
C     kRef   :: Pressure reference level
c     myThid :: thread number for this instance of the routine
      INTEGER myThid
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k
      INTEGER kRef
      _RL alphaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j
      _RL refTemp,refSalt,tP,sP
      _RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
      _RL ct, sa, sqrtsa, p
      _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
      _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
      _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rhoP0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dnum_dtheta, dden_dtheta
      _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
c     CALL LOOK_FOR_NEG_SALINITY(
c    I          iMin, iMax, jMin, jMax,
c    U          sFld,
c    I          k, bi, bj, myThid )
#endif

      IF (equationOfState.EQ.'LINEAR') THEN

         DO j=jMin,jMax
            DO i=iMin,iMax
               alphaLoc(i,j) = -rhonil * tAlpha
            ENDDO
         ENDDO

      ELSEIF (equationOfState.EQ.'POLY3') THEN

         refTemp=eosRefT(kRef)
         refSalt=eosRefS(kRef)

         DO j=jMin,jMax
            DO i=iMin,iMax
               tP=theta(i,j,k,bi,bj)-refTemp
               sP=salt(i,j,k,bi,bj)-refSalt
#ifdef USE_FACTORIZED_POLY
               alphaLoc(i,j) =
     &     ( eosC(6,kRef)
     &         *tP*3.
     &        +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
     &       )*tP
     &      +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
     &
#else
               alphaLoc(i,j) =
     &              eosC(1,kRef)                +
     &              eosC(3,kRef)*tP*2.          +
     &              eosC(4,kRef)         *sP    +
     &              eosC(6,kRef)*tP*tP*3.       +
     &              eosC(7,kRef)*tP*2.   *sP    +
     &              eosC(8,kRef)         *sP*sP
#endif
            ENDDO
         ENDDO

      ELSEIF ( equationOfState(1:5).EQ.'JMD95'
     &        .OR. equationOfState.EQ.'UNESCO' ) THEN
C     nonlinear equation of state in pressure coordinates

         CALL PRESSURE_FOR_EOS(
     I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O        locPres,
     I        myThid )

         CALL FIND_RHOP0(
     I        iMin, iMax, jMin, jMax,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoP0,
     I        myThid )

         CALL FIND_BULKMOD(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        bulkMod,
     I        myThid )

         DO j=jMin,jMax
            DO i=iMin,iMax

C     abbreviations
               t1 = theta(i,j,k,bi,bj)
               t2 = t1*t1
               t3 = t2*t1

#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
               s3o2 = 0. _d 0
#endif
               s1  = salt(i,j,k,bi,bj)
              IF ( s1 .GT. 0. _d 0 ) THEN
               s3o2 = SQRT(s1*s1*s1)
              ELSE
               s1   = 0. _d 0
               s3o2 = 0. _d 0
              ENDIF

               p1  = locPres(i,j)*SItoBar
               p2 = p1*p1

C     d(rho)/d(theta)
C     of fresh water at p = 0
               drhoP0dthetaFresh =
     &                eosJMDCFw(2)
     &           + 2.*eosJMDCFw(3)*t1
     &           + 3.*eosJMDCFw(4)*t2
     &           + 4.*eosJMDCFw(5)*t3
     &           + 5.*eosJMDCFw(6)*t3*t1
C     of salt water at p = 0
               drhoP0dthetaSalt =
     &        s1*(
     &                eosJMDCSw(2)
     &           + 2.*eosJMDCSw(3)*t1
     &           + 3.*eosJMDCSw(4)*t2
     &           + 4.*eosJMDCSw(5)*t3
     &           )
     &       + s3o2*(
     &           +    eosJMDCSw(7)
     &           + 2.*eosJMDCSw(8)*t1
     &           )
C     d(bulk modulus)/d(theta)
C     of fresh water at p = 0
               dKdthetaFresh =
     &                eosJMDCKFw(2)
     &           + 2.*eosJMDCKFw(3)*t1
     &           + 3.*eosJMDCKFw(4)*t2
     &           + 4.*eosJMDCKFw(5)*t3
C     of sea water at p = 0
               dKdthetaSalt =
     &        s1*(    eosJMDCKSw(2)
     &           + 2.*eosJMDCKSw(3)*t1
     &           + 3.*eosJMDCKSw(4)*t2
     &           )
     &    + s3o2*(    eosJMDCKSw(6)
     &           + 2.*eosJMDCKSw(7)*t1
     &           )
C     of sea water at p
               dKdthetaPres =
     &        p1*(    eosJMDCKP(2)
     &           + 2.*eosJMDCKP(3)*t1
     &           + 3.*eosJMDCKP(4)*t2
     &           )
     &   + p1*s1*(    eosJMDCKP(6)
     &           + 2.*eosJMDCKP(7)*t1
     &           )
     &      + p2*(    eosJMDCKP(10)
     &           + 2.*eosJMDCKP(11)*t1
     &           )
     &   + p2*s1*(    eosJMDCKP(13)
     &           + 2.*eosJMDCKP(14)*t1
     &           )

               drhoP0dtheta  = drhoP0dthetaFresh
     &                       + drhoP0dthetaSalt
               dKdtheta      = dKdthetaFresh
     &                       + dKdthetaSalt
     &                       + dKdthetaPres
               alphaLoc(i,j) =
     &              ( bulkmod(i,j)**2*drhoP0dtheta
     &              - bulkmod(i,j)*p1*drhoP0dtheta
     &              - rhoP0(i,j)*p1*dKdtheta )
     &              /( bulkmod(i,j) - p1 )**2

            ENDDO
         ENDDO
      ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN

         CALL PRESSURE_FOR_EOS(
     I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O        locPres,
     I        myThid )

         CALL FIND_RHONUM(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoLoc,
     I        myThid )

         CALL FIND_RHODEN(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoDen,
     I        myThid )

         DO j=jMin,jMax
            DO i=iMin,iMax
               t1  = theta(i,j,k,bi,bj)
               t2  = t1*t1
               s1  = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
               sp5 = 0. _d 0
#endif
              IF ( s1 .GT. 0. _d 0 ) THEN
               sp5 = SQRT(s1)
              ELSE
               s1  = 0. _d 0
               sp5 = 0. _d 0
              ENDIF

               p1   = locPres(i,j)*SItodBar
               p1t1 = p1*t1

               dnum_dtheta = eosMDJWFnum(1)
     &              + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
     &              + eosMDJWFnum(5)*s1
     &              + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)

               dden_dtheta = eosMDJWFden(1)
     &              + t1*(2.*eosMDJWFden(2)
     &              +     t1*(3.*eosMDJWFden(3)
     &              +         4.*eosMDJWFden(4)*t1 ) )
     &              + s1*(eosMDJWFden(6)
     &              +     t1*(3.*eosMDJWFden(7)*t1
     &              +         2.*eosMDJWFden(9)*sp5 ) )
     &              + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)

               alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
     &              - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)

         ENDDO
      ENDDO

      ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN

       CALL PRESSURE_FOR_EOS(
     I      bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O      locPres,
     I      myThid )

       CALL FIND_RHOTEOS(
     I      iMin, iMax, jMin, jMax, locPres,
     I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O      rhoLoc, rhoDen,
     I      myThid )

       DO j=jMin,jMax
        DO i=iMin,iMax
         ct      = theta(i,j,k,bi,bj)
         sa      = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
         sqrtsa  = 0. _d 0
#endif
         IF ( sa .GT. 0. _d 0 ) THEN
          sqrtsa = SQRT(sa)
         ELSE
          sa     = 0. _d 0
          sqrtsa = 0. _d 0
         ENDIF
         p       = locPres(i,j)*SItodBar

         dnum_dtheta = teos(02)
     &   + ct*(2.*teos(03) + 3.*teos(04)*ct)
     &   + sa*(teos(06) + 2.*teos(07)*ct
     &   + sqrtsa*(teos(09) + ct*(2.*teos(10) + 3.*teos(11)*ct)))
     &   + p*(     teos(13) + 2.*teos(14)*ct  + sa*2.*teos(16)
     &        + p*(teos(18) + 2.*teos(19)*ct))

         dden_dtheta = teos(22)
     &   + ct*(2.*teos(23) + ct*(3.*teos(24) + 4.*teos(25)*ct))
     &   + sa*(teos(27)
     &   + ct*(2.*teos(28) + ct*(3.*teos(29) + 4.*teos(30)*ct))
     &   + sqrtsa*(teos(32)
     &   + ct*(2.*teos(33) + ct*(3.*teos(34) + 4.*teos(35)*ct))))
     &   + p*(teos(38) + ct*(2.*teos(39) + 3.*teos(40)*ct)
     &   + teos(42)
     &   + p*(teos(44) + 2.*teos(45)*ct + teos(46)*sa
     &   + p*teos(48) ))

         alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
     &        - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)

        ENDDO
       ENDDO

      ELSE
         WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
         STOP 'FIND_ALPHA: "equationOfState" has illegal value'
      ENDIF

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: FIND_BETA
C     !INTERFACE:
      SUBROUTINE FIND_BETA (
     I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
     O     betaLoc,
     I     myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | o SUBROUTINE FIND_BETA
C     |   Calculates [drho(S,T,z) / dS] of a horizontal slice
C     |==========================================================*
C     | k       - is the Theta/Salt level
C     | kRef    - determines pressure reference level
C     |           (not used in 'LINEAR' mode)
C     | betaLoc - drho / dS (kg/m^3/PSU)
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     k      :: Level of Theta/Salt slice
C     kRef   :: Pressure reference level
c     myThid :: thread number for this instance of the routine
      INTEGER myThid
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER k
      INTEGER kRef
      _RL betaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)

C     !LOCAL VARIABLES:
C     == Local variables ==
      INTEGER i,j
      _RL refTemp,refSalt,tP,sP
      _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
      _RL ct, sa, sqrtsa, p
      _RL drhoP0dS
      _RL dKdS, dKdSSalt, dKdSPres
      _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rhoP0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL dnum_dsalt, dden_dsalt
      _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
c     CALL LOOK_FOR_NEG_SALINITY(
c    I          iMin, iMax, jMin, jMax,
c    U          sFld,
c    I          k, bi, bj, myThid )
#endif

      IF (equationOfState.EQ.'LINEAR') THEN

         DO j=jMin,jMax
            DO i=iMin,iMax
               betaLoc(i,j) = rhonil * sBeta
            ENDDO
         ENDDO

      ELSEIF (equationOfState.EQ.'POLY3') THEN

         refTemp=eosRefT(kRef)
         refSalt=eosRefS(kRef)

         DO j=jMin,jMax
            DO i=iMin,iMax
               tP=theta(i,j,k,bi,bj)-refTemp
               sP=salt(i,j,k,bi,bj)-refSalt
#ifdef USE_FACTORIZED_POLY
               betaLoc(i,j) =
     &    ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
     &   + ( eosC(7,kRef)*tP
     &      +eosC(8,kRef)*sP*2. + eosC(4,kRef)
     &     )*tP
#else
               betaLoc(i,j) =
     &              eosC(2,kRef)                +
     &              eosC(4,kRef)*tP             +
     &              eosC(5,kRef)         *sP*2. +
     &              eosC(7,kRef)*tP*tP          +
     &              eosC(8,kRef)*tP      *sP*2. +
     &              eosC(9,kRef)         *sP*sP*3.
#endif
            ENDDO
         ENDDO

      ELSEIF ( equationOfState(1:5).EQ.'JMD95'
     &        .OR. equationOfState.EQ.'UNESCO' ) THEN
C     nonlinear equation of state in pressure coordinates

         CALL PRESSURE_FOR_EOS(
     I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O        locPres,
     I        myThid )

         CALL FIND_RHOP0(
     I        iMin, iMax, jMin, jMax,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoP0,
     I        myThid )

         CALL FIND_BULKMOD(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        bulkMod,
     I        myThid )

         DO j=jMin,jMax
            DO i=iMin,iMax

C     abbreviations
               t1 = theta(i,j,k,bi,bj)
               t2 = t1*t1
               t3 = t2*t1

               s1  = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
               s3o2 = 0. _d 0
#endif
              IF ( s1 .GT. 0. _d 0 ) THEN
               s3o2 = 1.5*SQRT(s1)
              ELSE
               s1   = 0. _d 0
               s3o2 = 0. _d 0
              ENDIF

               p1  = locPres(i,j)*SItoBar

C     d(rho)/d(S)
C     of fresh water at p = 0
               drhoP0dS = 0. _d 0
C     of salt water at p = 0
               drhoP0dS = drhoP0dS
     &              + eosJMDCSw(1)
     &              + eosJMDCSw(2)*t1
     &              + eosJMDCSw(3)*t2
     &              + eosJMDCSw(4)*t3
     &              + eosJMDCSw(5)*t3*t1
     &       + s3o2*(
     &                eosJMDCSw(6)
     &              + eosJMDCSw(7)*t1
     &              + eosJMDCSw(8)*t2
     &              )
     &              + 2*eosJMDCSw(9)*s1
C     d(bulk modulus)/d(S)
C     of fresh water at p = 0
               dKdS = 0. _d 0
C     of sea water at p = 0
               dKdSSalt =
     &                eosJMDCKSw(1)
     &              + eosJMDCKSw(2)*t1
     &              + eosJMDCKSw(3)*t2
     &              + eosJMDCKSw(4)*t3
     &       + s3o2*( eosJMDCKSw(5)
     &              + eosJMDCKSw(6)*t1
     &              + eosJMDCKSw(7)*t2
     &              )

C     of sea water at p
               dKdSPres =
     &           p1*( eosJMDCKP(5)
     &              + eosJMDCKP(6)*t1
     &              + eosJMDCKP(7)*t2
     &              )
     &        + s3o2*p1*eosJMDCKP(8)
     &      + p1*p1*( eosJMDCKP(12)
     &              + eosJMDCKP(13)*t1
     &              + eosJMDCKP(14)*t2
     &              )

               dKdS = dKdSSalt + dKdSPres

               betaLoc(i,j) =
     &              ( bulkmod(i,j)**2*drhoP0dS
     &              - bulkmod(i,j)*p1*drhoP0dS
     &              - rhoP0(i,j)*p1*dKdS )
     &              /( bulkmod(i,j) - p1 )**2

            ENDDO
         ENDDO
      ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN

         CALL PRESSURE_FOR_EOS(
     I        bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O        locPres,
     I        myThid )

         CALL FIND_RHONUM(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoLoc,
     I        myThid )

         CALL FIND_RHODEN(
     I        iMin, iMax, jMin, jMax, locPres,
     I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O        rhoDen,
     I        myThid )

         DO j=jMin,jMax
            DO i=iMin,iMax
               t1  = theta(i,j,k,bi,bj)
               t2  = t1*t1
               s1  = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
               sp5 = 0. _d 0
#endif
              IF ( s1 .GT. 0. _d 0 ) THEN
               sp5 = SQRT(s1)
              ELSE
               s1  = 0. _d 0
               sp5 = 0. _d 0
              ENDIF

               p1   = locPres(i,j)*SItodBar
               p1t1 = p1*t1

               dnum_dsalt = eosMDJWFnum(4)
     &              + eosMDJWFnum(5)*t1
     &              + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
               dden_dsalt = eosMDJWFden(5)
     &              + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
     &              + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)

               betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
     &              - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )

            ENDDO
         ENDDO

      ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN

       CALL PRESSURE_FOR_EOS(
     I      bi, bj, iMin, iMax, jMin, jMax,  kRef,
     O      locPres,
     I      myThid )

       CALL FIND_RHOTEOS(
     I      iMin, iMax, jMin, jMax, locPres,
     I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
     O      rhoLoc, rhoDen,
     I      myThid )

       DO j=jMin,jMax
        DO i=iMin,iMax
         ct      = theta(i,j,k,bi,bj)
         sa      = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
C     this line makes TAF create vectorizable code
         sqrtsa  = 0. _d 0
#endif
         IF ( sa .GT. 0. _d 0 ) THEN
          sqrtsa = SQRT(sa)
         ELSE
          sa     = 0. _d 0
          sqrtsa = 0. _d 0
         ENDIF
         p       = locPres(i,j)*SItodBar

         dnum_dsalt =  teos(05) + ct*(teos(06) + teos(07)*ct)
     &   + 1.5*sqrtsa*(teos(08)
     &           + ct*(teos(09) + ct*(teos(10) + teos(11)*ct)))
     &   + p*(teos(15) + teos(16)*ct + p*teos(20))

         dden_dsalt = teos(26)
     &   + ct*(teos(27) + ct*(teos(28) + ct*(teos(29) + teos(30)*ct)))
     &   + 2.*teos(36)*sa
     &   + 1.5*sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
     &                          + ct*(teos(34) + teos(35)*ct))))
     &   + p*(teos(41) + teos(42)*ct + p*teos(46))

         betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
     &        - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )

        ENDDO
       ENDDO

      ELSE
         WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
         STOP 'FIND_BETA: "equationOfState" has illegal value'
      ENDIF

      RETURN
      END