model/src/find_alpha.F

C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.14 2003/09/10 22:21:22 jmc 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 )

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
      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)
      INTEGER myThid
CEOP

#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
      CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,
     &     sFld, 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        bi, bj, iMin, iMax, jMin, jMax, k,
     I        theta, salt,
     O        rhoP0,
     I        myThid )
         
         CALL FIND_BULKMOD(
     I        bi, bj, iMin, iMax, jMin, jMax,  k,
     I        locPres, theta, salt,
     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 ) 

c$$$         CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k,
c$$$     &        kRef, theta, salt, rhoLoc, myThid )

         CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k,
     &      locPres, theta, salt, rhoLoc, myThid )

         CALL FIND_RHODEN( bi, bj, iMin, iMax, jMin, jMax, k,
     &        locPres, theta, salt, rhoDen, 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 )
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
      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)
      INTEGER myThid
CEOP

#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
      CALL LOOK_FOR_NEG_SALINITY( bi, bj, iMin, iMax, jMin, jMax,  k,
     &     sFld, 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        bi, bj, iMin, iMax, jMin, jMax, k,
     I        theta, salt,
     O        rhoP0,
     I        myThid )
         
         CALL FIND_BULKMOD(
     I        bi, bj, iMin, iMax, jMin, jMax,  k,
     I        locPres, theta, salt,
     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 ) 

c$$$         CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax,  k,
c$$$     &        kRef, theta, salt, rhoLoc, myThid )

         CALL FIND_RHONUM( bi, bj, iMin, iMax, jMin, jMax, k,
     &      locPres, theta, salt, rhoLoc, myThid )

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