model/src/find_alpha.F

C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.16 2007/05/01 04:09:36 mlosch Exp $
C $Name:  $

#include "CPP_OPTIONS.h"
#define USE_FACTORIZED_POLY

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     |
C     | k - is the Theta/Salt level
C     | kRef - determines pressure reference level
C     |        (not used in 'LINEAR' mode)
C     |
C     | alphaLoc - drho / dT (kg/m^3/C)
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 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

               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 ( 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

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

      RETURN
      END

      SUBROUTINE FIND_BETA (
     I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
     O     betaLoc,
     I     myThid )
C     /==========================================================\
C     | o SUBROUTINE FIND_BETA                                   |
C     |   Calculates [drho(S,T,z) / dS] of a horizontal slice    |
C     |==========================================================|
C     |                                                          |
C     | k - is the Theta/Salt level                              |
C     | kRef - determines pressure reference level               |
C     |        (not used in 'LINEAR' mode)                       |
C     |                                                          |
C     | betaLoc - drho / dS (kg/m^3/PSU)                         |
C     |                                                          |
C     \==========================================================/
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"

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 ==
      INTEGER i,j
      _RL refTemp,refSalt,tP,sP
      _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
      _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 ( 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 ( 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

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

      RETURN
      END
1	jmc	1.17	C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.16 2007/05/01 04:09:36 mlosch Exp $
2	cnh	1.5	C $Name: $
3	adcroft	1.1
4			#include "CPP_OPTIONS.h"
5			#define USE_FACTORIZED_POLY
6
7	cnh	1.5	CBOP
8			C !ROUTINE: FIND_ALPHA
9			C !INTERFACE:
10			SUBROUTINE FIND_ALPHA (
11	jmc	1.17	I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
12			O alphaLoc,
13	mlosch	1.16	I myThid )
14	adcroft	1.1
15	cnh	1.5	C !DESCRIPTION: \bv
16			C ==========================================================
17	jmc	1.17	C \| o SUBROUTINE FIND_ALPHA
18			C \| Calculates [drho(S,T,z) / dT] of a horizontal slice
19	cnh	1.5	C ==========================================================
20	jmc	1.17	C \|
21			C \| k - is the Theta/Salt level
22			C \| kRef - determines pressure reference level
23			C \| (not used in 'LINEAR' mode)
24			C \|
25			C \| alphaLoc - drho / dT (kg/m^3/C)
26			C \|
27	cnh	1.5	C ==========================================================
28			C \ev
29
30			C !USES:
31			IMPLICIT NONE
32	jmc	1.12	C === Global variables ===
33	adcroft	1.1	#include "SIZE.h"
34			#include "DYNVARS.h"
35			#include "EEPARAMS.h"
36			#include "PARAMS.h"
37	mlosch	1.6	#include "EOS.h"
38			#include "GRID.h"
39	adcroft	1.1
40	cnh	1.5	C !INPUT/OUTPUT PARAMETERS:
41	jmc	1.12	C == Routine arguments ==
42	mlosch	1.16	C k :: Level of Theta/Salt slice
43			C kRef :: Pressure reference level
44			c myThid :: thread number for this instance of the routine
45			INTEGER myThid
46	jmc	1.12	INTEGER bi,bj,iMin,iMax,jMin,jMax
47			INTEGER k
48			INTEGER kRef
49	jmc	1.17	_RL alphaLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50	adcroft	1.1
51	cnh	1.5	C !LOCAL VARIABLES:
52	jmc	1.17	C == Local variables ==
53	jmc	1.12	INTEGER i,j
54	adcroft	1.1	_RL refTemp,refSalt,tP,sP
55	mlosch	1.9	_RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
56	mlosch	1.6	_RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
57			_RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
58	jmc	1.12	_RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
59			_RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
60	mlosch	1.6	_RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61	mlosch	1.9	_RL dnum_dtheta, dden_dtheta
62	jmc	1.12	_RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
63			_RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
64	cnh	1.5	CEOP
65	adcroft	1.1
66	mlosch	1.11	#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
67	jmc	1.17	c CALL LOOK_FOR_NEG_SALINITY(
68			c I iMin, iMax, jMin, jMax,
69			c U sFld,
70			c I k, bi, bj, myThid )
71	mlosch	1.11	#endif
72
73	jmc	1.12	IF (equationOfState.EQ.'LINEAR') THEN
74	adcroft	1.1
75	jmc	1.12	DO j=jMin,jMax
76			DO i=iMin,iMax
77	jmc	1.17	alphaLoc(i,j) = -rhonil * tAlpha
78	jmc	1.12	ENDDO
79			ENDDO
80	jmc	1.17
81	jmc	1.12	ELSEIF (equationOfState.EQ.'POLY3') THEN
82	adcroft	1.1
83			refTemp=eosRefT(kRef)
84			refSalt=eosRefS(kRef)
85
86	jmc	1.12	DO j=jMin,jMax
87			DO i=iMin,iMax
88	adcroft	1.1	tP=theta(i,j,k,bi,bj)-refTemp
89			sP=salt(i,j,k,bi,bj)-refSalt
90			#ifdef USE_FACTORIZED_POLY
91	jmc	1.17	alphaLoc(i,j) =
92	adcroft	1.1	& ( eosC(6,kRef)
93			& tP3.
94			& +(eosC(7,kRef)sP + eosC(3,kRef))2.
95			& )*tP
96			& +(eosC(8,kRef)sP + eosC(4,kRef) )sP + eosC(1,kRef)
97	jmc	1.17	&
98	adcroft	1.1	#else
99	jmc	1.17	alphaLoc(i,j) =
100	adcroft	1.1	& eosC(1,kRef) +
101			& eosC(3,kRef)tP2. +
102			& eosC(4,kRef) *sP +
103			& eosC(6,kRef)tPtP*3. +
104			& eosC(7,kRef)tP2. *sP +
105			& eosC(8,kRef) sPsP
106			#endif
107	jmc	1.12	ENDDO
108			ENDDO
109	jmc	1.17
110	jmc	1.12	ELSEIF ( equationOfState(1:5).EQ.'JMD95'
111			& .OR. equationOfState.EQ.'UNESCO' ) THEN
112	mlosch	1.6	C nonlinear equation of state in pressure coordinates
113
114	jmc	1.12	CALL PRESSURE_FOR_EOS(
115			I bi, bj, iMin, iMax, jMin, jMax, kRef,
116			O locPres,
117	jmc	1.17	I myThid )
118	jmc	1.12
119	mlosch	1.6	CALL FIND_RHOP0(
120	jmc	1.17	I iMin, iMax, jMin, jMax,
121			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
122	mlosch	1.6	O rhoP0,
123			I myThid )
124	jmc	1.17
125	mlosch	1.6	CALL FIND_BULKMOD(
126	jmc	1.17	I iMin, iMax, jMin, jMax, locPres,
127			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
128	mlosch	1.6	O bulkMod,
129			I myThid )
130
131	jmc	1.12	DO j=jMin,jMax
132			DO i=iMin,iMax
133	mlosch	1.6
134			C abbreviations
135	jmc	1.17	t1 = theta(i,j,k,bi,bj)
136	mlosch	1.9	t2 = t1*t1
137			t3 = t2*t1
138	jmc	1.17
139	mlosch	1.9	s1 = salt(i,j,k,bi,bj)
140	jmc	1.14	IF ( s1 .GT. 0. _d 0 ) THEN
141	jmc	1.12	s3o2 = SQRT(s1s1s1)
142	jmc	1.14	ELSE
143			s1 = 0. _d 0
144			s3o2 = 0. _d 0
145			ENDIF
146	jmc	1.17
147	jmc	1.12	p1 = locPres(i,j)*SItoBar
148	mlosch	1.9	p2 = p1*p1
149	mlosch	1.6
150			C d(rho)/d(theta)
151			C of fresh water at p = 0
152	jmc	1.17	drhoP0dthetaFresh =
153	mlosch	1.6	& eosJMDCFw(2)
154	mlosch	1.9	& + 2.eosJMDCFw(3)t1
155	mlosch	1.6	& + 3.eosJMDCFw(4)t2
156			& + 4.eosJMDCFw(5)t3
157	mlosch	1.9	& + 5.eosJMDCFw(6)t3*t1
158	mlosch	1.6	C of salt water at p = 0
159	jmc	1.17	drhoP0dthetaSalt =
160	mlosch	1.9	& s1*(
161	mlosch	1.6	& eosJMDCSw(2)
162	mlosch	1.9	& + 2.eosJMDCSw(3)t1
163	mlosch	1.6	& + 3.eosJMDCSw(4)t2
164			& + 4.eosJMDCSw(5)t3
165			& )
166			& + s3o2*(
167			& + eosJMDCSw(7)
168	mlosch	1.9	& + 2.eosJMDCSw(8)t1
169	mlosch	1.6	& )
170			C d(bulk modulus)/d(theta)
171			C of fresh water at p = 0
172	jmc	1.17	dKdthetaFresh =
173	mlosch	1.6	& eosJMDCKFw(2)
174	mlosch	1.9	& + 2.eosJMDCKFw(3)t1
175	mlosch	1.6	& + 3.eosJMDCKFw(4)t2
176			& + 4.eosJMDCKFw(5)t3
177			C of sea water at p = 0
178			dKdthetaSalt =
179	mlosch	1.9	& s1*( eosJMDCKSw(2)
180			& + 2.eosJMDCKSw(3)t1
181	mlosch	1.6	& + 3.eosJMDCKSw(4)t2
182			& )
183			& + s3o2*( eosJMDCKSw(6)
184	mlosch	1.9	& + 2.eosJMDCKSw(7)t1
185	mlosch	1.6	& )
186			C of sea water at p
187			dKdthetaPres =
188	mlosch	1.9	& p1*( eosJMDCKP(2)
189			& + 2.eosJMDCKP(3)t1
190	mlosch	1.6	& + 3.eosJMDCKP(4)t2
191			& )
192	mlosch	1.9	& + p1s1( eosJMDCKP(6)
193			& + 2.eosJMDCKP(7)t1
194	mlosch	1.6	& )
195			& + p2*( eosJMDCKP(10)
196	mlosch	1.9	& + 2.eosJMDCKP(11)t1
197	mlosch	1.6	& )
198	mlosch	1.9	& + p2s1( eosJMDCKP(13)
199			& + 2.eosJMDCKP(14)t1
200	mlosch	1.6	& )
201
202	jmc	1.17	drhoP0dtheta = drhoP0dthetaFresh
203	mlosch	1.6	& + drhoP0dthetaSalt
204	jmc	1.17	dKdtheta = dKdthetaFresh
205			& + dKdthetaSalt
206	mlosch	1.6	& + dKdthetaPres
207	jmc	1.17	alphaLoc(i,j) =
208	mlosch	1.6	& ( bulkmod(i,j)*2drhoP0dtheta
209	mlosch	1.9	& - bulkmod(i,j)p1drhoP0dtheta
210			& - rhoP0(i,j)p1dKdtheta )
211			& /( bulkmod(i,j) - p1 )**2
212	jmc	1.17
213
214	jmc	1.12	ENDDO
215			ENDDO
216			ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
217
218			CALL PRESSURE_FOR_EOS(
219			I bi, bj, iMin, iMax, jMin, jMax, kRef,
220			O locPres,
221	jmc	1.17	I myThid )
222	jmc	1.12
223	jmc	1.17	CALL FIND_RHONUM(
224			I iMin, iMax, jMin, jMax, locPres,
225			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
226			O rhoLoc,
227			I myThid )
228	heimbach	1.15
229	jmc	1.17	CALL FIND_RHODEN(
230			I iMin, iMax, jMin, jMax, locPres,
231			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
232			O rhoDen,
233			I myThid )
234	mlosch	1.9
235	jmc	1.12	DO j=jMin,jMax
236			DO i=iMin,iMax
237	mlosch	1.9	t1 = theta(i,j,k,bi,bj)
238			t2 = t1*t1
239			s1 = salt(i,j,k,bi,bj)
240	jmc	1.14	IF ( s1 .GT. 0. _d 0 ) THEN
241	jmc	1.17	sp5 = SQRT(s1)
242	jmc	1.14	ELSE
243			s1 = 0. _d 0
244			sp5 = 0. _d 0
245			ENDIF
246	mlosch	1.9
247	jmc	1.12	p1 = locPres(i,j)*SItodBar
248	mlosch	1.9	p1t1 = p1*t1
249	jmc	1.17
250	mlosch	1.9	dnum_dtheta = eosMDJWFnum(1)
251	jmc	1.17	& + t1(2.eosMDJWFnum(2) + 3.eosMDJWFnum(3)t1)
252			& + eosMDJWFnum(5)*s1
253	jmc	1.14	& + p1t1(2.eosMDJWFnum(8) + 2.eosMDJWFnum(11)p1)
254	jmc	1.17
255			dden_dtheta = eosMDJWFden(1)
256			& + t1(2.eosMDJWFden(2)
257			& + t1(3.eosMDJWFden(3)
258	mlosch	1.9	& + 4.eosMDJWFden(4)t1 ) )
259	jmc	1.17	& + s1*(eosMDJWFden(6)
260			& + t1(3.eosMDJWFden(7)*t1
261	mlosch	1.9	& + 2.eosMDJWFden(9)sp5 ) )
262			& + p1p1(3.eosMDJWFden(11)t2 + eosMDJWFden(12)*p1)
263	jmc	1.17
264			alphaLoc(i,j) = rhoDen(i,j)*(dnum_dtheta
265	heimbach	1.15	& - (rhoLoc(i,j)rhoDen(i,j))dden_dtheta)
266	jmc	1.17
267	jmc	1.12	ENDDO
268			ENDDO
269	mlosch	1.9
270	jmc	1.12	ELSE
271			WRITE(,) 'FIND_ALPHA: equationOfState = ',equationOfState
272			STOP 'FIND_ALPHA: "equationOfState" has illegal value'
273			ENDIF
274	adcroft	1.1
275	jmc	1.17	RETURN
276	jmc	1.12	END
277	adcroft	1.1
278	jmc	1.12	SUBROUTINE FIND_BETA (
279	mlosch	1.10	I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
280	jmc	1.17	O betaLoc,
281	mlosch	1.16	I myThid )
282	adcroft	1.1	C /==========================================================\
283			C \| o SUBROUTINE FIND_BETA \|
284			C \| Calculates [drho(S,T,z) / dS] of a horizontal slice \|
285			C \|==========================================================\|
286			C \| \|
287			C \| k - is the Theta/Salt level \|
288			C \| kRef - determines pressure reference level \|
289			C \| (not used in 'LINEAR' mode) \|
290			C \| \|
291	jmc	1.17	C \| betaLoc - drho / dS (kg/m^3/PSU) \|
292	adcroft	1.1	C \| \|
293			C \==========================================================/
294	jmc	1.12	IMPLICIT NONE
295	adcroft	1.1
296	jmc	1.12	C === Global variables ===
297	adcroft	1.1	#include "SIZE.h"
298			#include "DYNVARS.h"
299			#include "EEPARAMS.h"
300			#include "PARAMS.h"
301	mlosch	1.6	#include "EOS.h"
302			#include "GRID.h"
303	adcroft	1.1
304	jmc	1.12	C == Routine arguments ==
305	mlosch	1.16	C k :: Level of Theta/Salt slice
306			C kRef :: Pressure reference level
307			c myThid :: thread number for this instance of the routine
308			INTEGER myThid
309	jmc	1.12	INTEGER bi,bj,iMin,iMax,jMin,jMax
310			INTEGER k
311			INTEGER kRef
312	jmc	1.17	_RL betaLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
313	adcroft	1.1
314	jmc	1.12	C == Local variables ==
315			INTEGER i,j
316	adcroft	1.1	_RL refTemp,refSalt,tP,sP
317	mlosch	1.9	_RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
318	mlosch	1.6	_RL drhoP0dS
319			_RL dKdS, dKdSSalt, dKdSPres
320	jmc	1.12	_RL locPres(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
321			_RL rhoP0 (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
322	mlosch	1.6	_RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
323	mlosch	1.9	_RL dnum_dsalt, dden_dsalt
324	jmc	1.12	_RL rhoDen (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
325			_RL rhoLoc (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
326	mlosch	1.6	CEOP
327
328	mlosch	1.11	#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
329	jmc	1.17	c CALL LOOK_FOR_NEG_SALINITY(
330			c I iMin, iMax, jMin, jMax,
331			c U sFld,
332			c I k, bi, bj, myThid )
333	mlosch	1.11	#endif
334
335	jmc	1.12	IF (equationOfState.EQ.'LINEAR') THEN
336	adcroft	1.1
337	jmc	1.12	DO j=jMin,jMax
338			DO i=iMin,iMax
339	jmc	1.17	betaLoc(i,j) = rhonil * sBeta
340	jmc	1.12	ENDDO
341			ENDDO
342	jmc	1.17
343	jmc	1.12	ELSEIF (equationOfState.EQ.'POLY3') THEN
344	adcroft	1.1
345			refTemp=eosRefT(kRef)
346			refSalt=eosRefS(kRef)
347
348	jmc	1.12	DO j=jMin,jMax
349			DO i=iMin,iMax
350	adcroft	1.1	tP=theta(i,j,k,bi,bj)-refTemp
351			sP=salt(i,j,k,bi,bj)-refSalt
352			#ifdef USE_FACTORIZED_POLY
353	jmc	1.17	betaLoc(i,j) =
354	adcroft	1.1	& ( eosC(9,kRef)sP3. + eosC(5,kRef)2. )sP + eosC(2,kRef)
355			& + ( eosC(7,kRef)*tP
356			& +eosC(8,kRef)sP2. + eosC(4,kRef)
357			& )*tP
358			#else
359	jmc	1.17	betaLoc(i,j) =
360	adcroft	1.1	& eosC(2,kRef) +
361			& eosC(4,kRef)*tP +
362			& eosC(5,kRef) sP2. +
363			& eosC(7,kRef)tPtP +
364			& eosC(8,kRef)tP sP*2. +
365			& eosC(9,kRef) sPsP*3.
366			#endif
367	jmc	1.12	ENDDO
368			ENDDO
369	adcroft	1.1
370	jmc	1.12	ELSEIF ( equationOfState(1:5).EQ.'JMD95'
371			& .OR. equationOfState.EQ.'UNESCO' ) THEN
372	mlosch	1.6	C nonlinear equation of state in pressure coordinates
373
374	jmc	1.12	CALL PRESSURE_FOR_EOS(
375			I bi, bj, iMin, iMax, jMin, jMax, kRef,
376			O locPres,
377	jmc	1.17	I myThid )
378	jmc	1.12
379	mlosch	1.6	CALL FIND_RHOP0(
380	jmc	1.17	I iMin, iMax, jMin, jMax,
381			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
382	mlosch	1.6	O rhoP0,
383			I myThid )
384	jmc	1.17
385	mlosch	1.6	CALL FIND_BULKMOD(
386	jmc	1.17	I iMin, iMax, jMin, jMax, locPres,
387			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
388	mlosch	1.6	O bulkMod,
389			I myThid )
390
391	jmc	1.12	DO j=jMin,jMax
392			DO i=iMin,iMax
393	mlosch	1.6
394			C abbreviations
395	jmc	1.17	t1 = theta(i,j,k,bi,bj)
396	mlosch	1.9	t2 = t1*t1
397			t3 = t2*t1
398	jmc	1.17
399	mlosch	1.9	s1 = salt(i,j,k,bi,bj)
400	jmc	1.14	IF ( s1 .GT. 0. _d 0 ) THEN
401	jmc	1.12	s3o2 = 1.5*SQRT(s1)
402	jmc	1.14	ELSE
403			s1 = 0. _d 0
404			s3o2 = 0. _d 0
405			ENDIF
406	mlosch	1.6
407	jmc	1.12	p1 = locPres(i,j)*SItoBar
408	mlosch	1.6
409			C d(rho)/d(S)
410			C of fresh water at p = 0
411			drhoP0dS = 0. _d 0
412			C of salt water at p = 0
413			drhoP0dS = drhoP0dS
414			& + eosJMDCSw(1)
415	mlosch	1.9	& + eosJMDCSw(2)*t1
416	mlosch	1.6	& + eosJMDCSw(3)*t2
417			& + eosJMDCSw(4)*t3
418	mlosch	1.9	& + eosJMDCSw(5)t3t1
419	mlosch	1.6	& + s3o2*(
420			& eosJMDCSw(6)
421	mlosch	1.9	& + eosJMDCSw(7)*t1
422	mlosch	1.6	& + eosJMDCSw(8)*t2
423			& )
424	mlosch	1.9	& + 2eosJMDCSw(9)s1
425	mlosch	1.6	C d(bulk modulus)/d(S)
426			C of fresh water at p = 0
427			dKdS = 0. _d 0
428			C of sea water at p = 0
429			dKdSSalt =
430			& eosJMDCKSw(1)
431	mlosch	1.9	& + eosJMDCKSw(2)*t1
432	mlosch	1.6	& + eosJMDCKSw(3)*t2
433			& + eosJMDCKSw(4)*t3
434			& + s3o2*( eosJMDCKSw(5)
435	mlosch	1.9	& + eosJMDCKSw(6)*t1
436	mlosch	1.6	& + eosJMDCKSw(7)*t2
437			& )
438
439			C of sea water at p
440	jmc	1.17	dKdSPres =
441	mlosch	1.9	& p1*( eosJMDCKP(5)
442			& + eosJMDCKP(6)*t1
443	mlosch	1.6	& + eosJMDCKP(7)*t2
444			& )
445	mlosch	1.9	& + s3o2p1eosJMDCKP(8)
446			& + p1p1( eosJMDCKP(12)
447			& + eosJMDCKP(13)*t1
448	mlosch	1.6	& + eosJMDCKP(14)*t2
449			& )
450
451			dKdS = dKdSSalt + dKdSPres
452
453	jmc	1.17	betaLoc(i,j) =
454	mlosch	1.6	& ( bulkmod(i,j)*2drhoP0dS
455	mlosch	1.9	& - bulkmod(i,j)p1drhoP0dS
456			& - rhoP0(i,j)p1dKdS )
457			& /( bulkmod(i,j) - p1 )**2
458	jmc	1.17
459
460	jmc	1.12	ENDDO
461			ENDDO
462			ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
463
464			CALL PRESSURE_FOR_EOS(
465			I bi, bj, iMin, iMax, jMin, jMax, kRef,
466			O locPres,
467	jmc	1.17	I myThid )
468	jmc	1.12
469	jmc	1.17	CALL FIND_RHONUM(
470			I iMin, iMax, jMin, jMax, locPres,
471			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
472			O rhoLoc,
473			I myThid )
474	heimbach	1.15
475	jmc	1.17	CALL FIND_RHODEN(
476			I iMin, iMax, jMin, jMax, locPres,
477			I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
478			O rhoDen,
479			I myThid )
480	heimbach	1.13
481	jmc	1.12	DO j=jMin,jMax
482			DO i=iMin,iMax
483	mlosch	1.9	t1 = theta(i,j,k,bi,bj)
484			t2 = t1*t1
485			s1 = salt(i,j,k,bi,bj)
486	jmc	1.14	IF ( s1 .GT. 0. _d 0 ) THEN
487	jmc	1.17	sp5 = SQRT(s1)
488	jmc	1.14	ELSE
489			s1 = 0. _d 0
490			sp5 = 0. _d 0
491			ENDIF
492	jmc	1.17
493	jmc	1.12	p1 = locPres(i,j)*SItodBar
494	mlosch	1.9	p1t1 = p1*t1
495	jmc	1.17
496			dnum_dsalt = eosMDJWFnum(4)
497	mlosch	1.9	& + eosMDJWFnum(5)*t1
498			& + 2.eosMDJWFnum(6)s1 + eosMDJWFnum(9)*p1
499	jmc	1.17	dden_dsalt = eosMDJWFden(5)
500			& + t1( eosMDJWFden(6) + eosMDJWFden(7)t2 )
501	mlosch	1.9	& + 1.5sp5(eosMDJWFden(8) + eosMDJWFden(9)*t2)
502
503	jmc	1.17	betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
504	heimbach	1.15	& - (rhoLoc(i,j)rhoDen(i,j))dden_dsalt )
505	mlosch	1.9
506	jmc	1.12	ENDDO
507			ENDDO
508	mlosch	1.9
509	jmc	1.12	ELSE
510			WRITE(,) 'FIND_BETA: equationOfState = ',equationOfState
511			STOP 'FIND_BETA: "equationOfState" has illegal value'
512			ENDIF
513	adcroft	1.1
514	jmc	1.17	RETURN
515	jmc	1.12	END