model/src/find_alpha.F

C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.7 2002/08/15 18:48:47 adcroft 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, eqn,
     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     | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'   
C     |                                                           
C     | alphaloc - drho / dT (kg/m^3/C)                           
C     |                                                           
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
c Common
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"

C     !INPUT/OUTPUT PARAMETERS:
c Arguments
      integer bi,bj,iMin,iMax,jMin,jMax
      integer k                 ! Level of Theta/Salt slice
      integer kRef              ! Pressure reference level
      character*(*) eqn
      _RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)

C     !LOCAL VARIABLES:
c Local
      integer i,j
      _RL refTemp,refSalt,tP,sP
      _RL t, t2, t3, s, s3o2, p, p2
      _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
      _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
      _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      integer myThid
CEOP

Cml      stop 'myThid is not properly defined in this subroutine'

      if (eqn.eq.'LINEAR') then

         do j=jMin,jMax
            do i=iMin,iMax
               alphaloc(i,j) = -rhonil * tAlpha
            enddo
         enddo
         
      elseif (eqn.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 ( eqn(1:5).eq.'JMD95' ) then
C     nonlinear equation of state in pressure coordinates

         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, kRef,
     I        theta, salt,
     O        bulkMod,
     I        myThid )

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

C     abbreviations
               t  = theta(i,j,k,bi,bj) 
               t2 = t*t
               t3 = t2*t
               
               s  = salt(i,j,k,bi,bj)
               if ( s .lt. 0. _d 0 ) then
C     issue a warning
                  write(*,'(a,i3,a,i3,a,i3,a,e13.5)') 
     &             ' FIND_ALPHA: WARNING, salinity(',
     &                 i,',',j,',',k,') = ', s
                  s = 0. _d 0
               end if
               s3o2 = sqrt(s*s*s)
               
               p  = pressure(i,j,kRef,bi,bj)*SItoBar
               p2 = p*p

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

               drhoP0dtheta  = drhoP0dthetaFresh 
     &                       + drhoP0dthetaSalt
               dKdtheta      = dKdthetaFresh 
     &                       + dKdthetaSalt 
     &                       + dKdthetaPres
               alphaloc(i,j) = 
     &              ( bulkmod(i,j)**2*drhoP0dtheta
     &              - bulkmod(i,j)*p*drhoP0dtheta
     &              - rhoP0(i,j)*p*dKdtheta )
     &              /( bulkmod(i,j) - p )**2
               
              
            end do
         end do
      else
         write(*,*) 'FIND_ALPHA: eqn = ',eqn
         stop 'FIND_ALPHA: argument "eqn" has illegal value'
      endif

      return
      end

      subroutine FIND_BETA (
     I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef, eqn,
     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     | eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'  |
C     |                                                          |
C     | betaloc - drho / dS (kg/m^3/PSU)                         |
C     |                                                          |
C     \==========================================================/
      implicit none

c Common
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"

c Arguments
      integer bi,bj,iMin,iMax,jMin,jMax
      integer k                 ! Level of Theta/Salt slice
      integer kRef              ! Pressure reference level
      character*(*) eqn
      _RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)

c Local
      integer i,j
      _RL refTemp,refSalt,tP,sP
      _RL t, t2, t3, s, s3o2, p
      _RL drhoP0dS
      _RL dKdS, dKdSSalt, dKdSPres
      _RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      integer myThid
CEOP

Cml      stop 'myThid is not properly defined in this subroutine'

      if (eqn.eq.'LINEAR') then

         do j=jMin,jMax
            do i=iMin,iMax
               betaloc(i,j) = rhonil * sBeta
            enddo
         enddo
         
      elseif (eqn.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 ( eqn(1:5).eq.'JMD95' ) then
C     nonlinear equation of state in pressure coordinates

         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, kRef,
     I        theta, salt,
     O        bulkMod,
     I        myThid )

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

C     abbreviations
               t  = theta(i,j,k,bi,bj) 
               t2 = t*t
               t3 = t2*t
               
               s  = salt(i,j,k,bi,bj)
               if ( s .lt. 0. _d 0 ) then
C     issue a warning
                  write(*,'(a,i3,a,i3,a,i3,a,e13.5)') 
     &             ' FIND_BETA: WARNING, salinity(',
     &                 i,',',j,',',k,') = ', s
                  s = 0. _d 0
               end if
               s3o2 = 1.5*sqrt(s)

               p  = pressure(i,j,kRef,bi,bj)*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)*t
     &              + eosJMDCSw(3)*t2
     &              + eosJMDCSw(4)*t3
     &              + eosJMDCSw(5)*t3*t
     &       + s3o2*(
     &                eosJMDCSw(6)
     &              + eosJMDCSw(7)*t
     &              + eosJMDCSw(8)*t2
     &              )
     &              + 2*eosJMDCSw(9)*s
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)*t
     &              + eosJMDCKSw(3)*t2
     &              + eosJMDCKSw(4)*t3
     &       + s3o2*( eosJMDCKSw(5)
     &              + eosJMDCKSw(6)*t
     &              + eosJMDCKSw(7)*t2
     &              )

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

               dKdS = dKdSSalt + dKdSPres

               betaloc(i,j) = 
     &              ( bulkmod(i,j)**2*drhoP0dS
     &              - bulkmod(i,j)*p*drhoP0dS
     &              - rhoP0(i,j)*p*dKdS )
     &              /( bulkmod(i,j) - p )**2
               
              
            end do
         end do
      else
         write(*,*) 'FIND_BETA: eqn = ',eqn
         stop 'FIND_BETA: argument "eqn" has illegal value'
      endif

      return
      end
1	mlosch	1.8	C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.7 2002/08/15 18:48:47 adcroft 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	adcroft	1.1	I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
12			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 \| eqn - determines the eqn. of state: 'LINEAR' or 'POLY3'
24			C \|
25			C \| alphaloc - drho / dT (kg/m^3/C)
26			C \|
27			C ==========================================================
28			C \ev
29
30			C !USES:
31			IMPLICIT NONE
32	adcroft	1.1	c Common
33			#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	adcroft	1.1	c Arguments
42			integer bi,bj,iMin,iMax,jMin,jMax
43			integer k ! Level of Theta/Salt slice
44			integer kRef ! Pressure reference level
45			character() eqn
46			_RL alphaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
47
48	cnh	1.5	C !LOCAL VARIABLES:
49	adcroft	1.1	c Local
50			integer i,j
51			_RL refTemp,refSalt,tP,sP
52	mlosch	1.6	_RL t, t2, t3, s, s3o2, p, p2
53			_RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
54			_RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
55			_RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
56			_RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
57			integer myThid
58	cnh	1.5	CEOP
59	adcroft	1.1
60	mlosch	1.6	Cml stop 'myThid is not properly defined in this subroutine'
61
62	adcroft	1.1	if (eqn.eq.'LINEAR') then
63
64			do j=jMin,jMax
65			do i=iMin,iMax
66			alphaloc(i,j) = -rhonil * tAlpha
67			enddo
68			enddo
69
70			elseif (eqn.eq.'POLY3') then
71
72			refTemp=eosRefT(kRef)
73			refSalt=eosRefS(kRef)
74
75			do j=jMin,jMax
76			do i=iMin,iMax
77			tP=theta(i,j,k,bi,bj)-refTemp
78			sP=salt(i,j,k,bi,bj)-refSalt
79			#ifdef USE_FACTORIZED_POLY
80			alphaloc(i,j) =
81			& ( eosC(6,kRef)
82			& tP3.
83			& +(eosC(7,kRef)sP + eosC(3,kRef))2.
84			& )*tP
85			& +(eosC(8,kRef)sP + eosC(4,kRef) )sP + eosC(1,kRef)
86			&
87			#else
88			alphaloc(i,j) =
89			& eosC(1,kRef) +
90			& eosC(3,kRef)tP2. +
91			& eosC(4,kRef) *sP +
92			& eosC(6,kRef)tPtP*3. +
93			& eosC(7,kRef)tP2. *sP +
94			& eosC(8,kRef) sPsP
95			#endif
96			enddo
97			enddo
98
99	mlosch	1.6	elseif ( eqn(1:5).eq.'JMD95' ) then
100			C nonlinear equation of state in pressure coordinates
101
102			CALL FIND_RHOP0(
103			I bi, bj, iMin, iMax, jMin, jMax, k,
104			I theta, salt,
105			O rhoP0,
106			I myThid )
107
108			CALL FIND_BULKMOD(
109	adcroft	1.7	I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
110	mlosch	1.6	I theta, salt,
111			O bulkMod,
112			I myThid )
113
114			do j=jMin,jMax
115			do i=iMin,iMax
116
117			C abbreviations
118			t = theta(i,j,k,bi,bj)
119			t2 = t*t
120			t3 = t2*t
121
122	mlosch	1.8	s = salt(i,j,k,bi,bj)
123			if ( s .lt. 0. _d 0 ) then
124	mlosch	1.6	C issue a warning
125			write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
126	mlosch	1.8	& ' FIND_ALPHA: WARNING, salinity(',
127			& i,',',j,',',k,') = ', s
128	mlosch	1.6	s = 0. _d 0
129			end if
130			s3o2 = sqrt(sss)
131
132	adcroft	1.7	p = pressure(i,j,kRef,bi,bj)*SItoBar
133	mlosch	1.6	p2 = p*p
134
135			C d(rho)/d(theta)
136			C of fresh water at p = 0
137			drhoP0dthetaFresh =
138			& eosJMDCFw(2)
139			& + 2.eosJMDCFw(3)t
140			& + 3.eosJMDCFw(4)t2
141			& + 4.eosJMDCFw(5)t3
142			& + 5.eosJMDCFw(6)t3*t
143			C of salt water at p = 0
144			drhoP0dthetaSalt =
145			& s*(
146			& eosJMDCSw(2)
147			& + 2.eosJMDCSw(3)t
148			& + 3.eosJMDCSw(4)t2
149			& + 4.eosJMDCSw(5)t3
150			& )
151			& + s3o2*(
152			& + eosJMDCSw(7)
153			& + 2.eosJMDCSw(8)t
154			& )
155			C d(bulk modulus)/d(theta)
156			C of fresh water at p = 0
157			dKdthetaFresh =
158			& eosJMDCKFw(2)
159			& + 2.eosJMDCKFw(3)t
160			& + 3.eosJMDCKFw(4)t2
161			& + 4.eosJMDCKFw(5)t3
162			C of sea water at p = 0
163			dKdthetaSalt =
164			& s*( eosJMDCKSw(2)
165			& + 2.eosJMDCKSw(3)t
166			& + 3.eosJMDCKSw(4)t2
167			& )
168			& + s3o2*( eosJMDCKSw(6)
169			& + 2.eosJMDCKSw(7)t
170			& )
171			C of sea water at p
172			dKdthetaPres =
173			& p*( eosJMDCKP(2)
174			& + 2.eosJMDCKP(3)t
175			& + 3.eosJMDCKP(4)t2
176			& )
177			& + ps( eosJMDCKP(6)
178			& + 2.eosJMDCKP(7)t
179			& )
180			& + p2*( eosJMDCKP(10)
181			& + 2.eosJMDCKP(11)t
182			& )
183			& + p2s( eosJMDCKP(13)
184			& + 2.eosJMDCKP(14)t
185			& )
186
187			drhoP0dtheta = drhoP0dthetaFresh
188			& + drhoP0dthetaSalt
189			dKdtheta = dKdthetaFresh
190			& + dKdthetaSalt
191			& + dKdthetaPres
192			alphaloc(i,j) =
193			& ( bulkmod(i,j)*2drhoP0dtheta
194			& - bulkmod(i,j)pdrhoP0dtheta
195			& - rhoP0(i,j)pdKdtheta )
196			& /( bulkmod(i,j) - p )**2
197
198
199			end do
200			end do
201	adcroft	1.1	else
202	heimbach	1.4	write(,) 'FIND_ALPHA: eqn = ',eqn
203	adcroft	1.1	stop 'FIND_ALPHA: argument "eqn" has illegal value'
204			endif
205
206			return
207			end
208
209			subroutine FIND_BETA (
210			I bi, bj, iMin, iMax, jMin, jMax, k, kRef, eqn,
211			O betaloc )
212			C /==========================================================\
213			C \| o SUBROUTINE FIND_BETA \|
214			C \| Calculates [drho(S,T,z) / dS] of a horizontal slice \|
215			C \|==========================================================\|
216			C \| \|
217			C \| k - is the Theta/Salt level \|
218			C \| kRef - determines pressure reference level \|
219			C \| (not used in 'LINEAR' mode) \|
220			C \| eqn - determines the eqn. of state: 'LINEAR' or 'POLY3' \|
221			C \| \|
222			C \| betaloc - drho / dS (kg/m^3/PSU) \|
223			C \| \|
224			C \==========================================================/
225			implicit none
226
227			c Common
228			#include "SIZE.h"
229			#include "DYNVARS.h"
230			#include "EEPARAMS.h"
231			#include "PARAMS.h"
232	mlosch	1.6	#include "EOS.h"
233			#include "GRID.h"
234	adcroft	1.1
235			c Arguments
236			integer bi,bj,iMin,iMax,jMin,jMax
237			integer k ! Level of Theta/Salt slice
238			integer kRef ! Pressure reference level
239			character() eqn
240			_RL betaloc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
241
242			c Local
243			integer i,j
244			_RL refTemp,refSalt,tP,sP
245	mlosch	1.6	_RL t, t2, t3, s, s3o2, p
246			_RL drhoP0dS
247			_RL dKdS, dKdSSalt, dKdSPres
248			_RL rhoP0(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
249			_RL bulkMod(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
250			integer myThid
251			CEOP
252
253			Cml stop 'myThid is not properly defined in this subroutine'
254	adcroft	1.1
255			if (eqn.eq.'LINEAR') then
256
257			do j=jMin,jMax
258			do i=iMin,iMax
259			betaloc(i,j) = rhonil * sBeta
260			enddo
261			enddo
262
263			elseif (eqn.eq.'POLY3') then
264
265			refTemp=eosRefT(kRef)
266			refSalt=eosRefS(kRef)
267
268			do j=jMin,jMax
269			do i=iMin,iMax
270			tP=theta(i,j,k,bi,bj)-refTemp
271			sP=salt(i,j,k,bi,bj)-refSalt
272			#ifdef USE_FACTORIZED_POLY
273			betaloc(i,j) =
274			& ( eosC(9,kRef)sP3. + eosC(5,kRef)2. )sP + eosC(2,kRef)
275			& + ( eosC(7,kRef)*tP
276			& +eosC(8,kRef)sP2. + eosC(4,kRef)
277			& )*tP
278			#else
279			betaloc(i,j) =
280			& eosC(2,kRef) +
281			& eosC(4,kRef)*tP +
282			& eosC(5,kRef) sP2. +
283			& eosC(7,kRef)tPtP +
284			& eosC(8,kRef)tP sP*2. +
285			& eosC(9,kRef) sPsP*3.
286			#endif
287			enddo
288			enddo
289
290	mlosch	1.6	elseif ( eqn(1:5).eq.'JMD95' ) then
291			C nonlinear equation of state in pressure coordinates
292
293			CALL FIND_RHOP0(
294			I bi, bj, iMin, iMax, jMin, jMax, k,
295			I theta, salt,
296			O rhoP0,
297			I myThid )
298
299			CALL FIND_BULKMOD(
300	adcroft	1.7	I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
301	mlosch	1.6	I theta, salt,
302			O bulkMod,
303			I myThid )
304
305			do j=jMin,jMax
306			do i=iMin,iMax
307
308			C abbreviations
309			t = theta(i,j,k,bi,bj)
310			t2 = t*t
311			t3 = t2*t
312
313	mlosch	1.8	s = salt(i,j,k,bi,bj)
314			if ( s .lt. 0. _d 0 ) then
315	mlosch	1.6	C issue a warning
316			write(*,'(a,i3,a,i3,a,i3,a,e13.5)')
317	mlosch	1.8	& ' FIND_BETA: WARNING, salinity(',
318			& i,',',j,',',k,') = ', s
319	mlosch	1.6	s = 0. _d 0
320			end if
321			s3o2 = 1.5*sqrt(s)
322
323	adcroft	1.7	p = pressure(i,j,kRef,bi,bj)*SItoBar
324	mlosch	1.6
325			C d(rho)/d(S)
326			C of fresh water at p = 0
327			drhoP0dS = 0. _d 0
328			C of salt water at p = 0
329			drhoP0dS = drhoP0dS
330			& + eosJMDCSw(1)
331			& + eosJMDCSw(2)*t
332			& + eosJMDCSw(3)*t2
333			& + eosJMDCSw(4)*t3
334			& + eosJMDCSw(5)t3t
335			& + s3o2*(
336			& eosJMDCSw(6)
337			& + eosJMDCSw(7)*t
338			& + eosJMDCSw(8)*t2
339			& )
340			& + 2eosJMDCSw(9)s
341			C d(bulk modulus)/d(S)
342			C of fresh water at p = 0
343			dKdS = 0. _d 0
344			C of sea water at p = 0
345			dKdSSalt =
346			& eosJMDCKSw(1)
347			& + eosJMDCKSw(2)*t
348			& + eosJMDCKSw(3)*t2
349			& + eosJMDCKSw(4)*t3
350			& + s3o2*( eosJMDCKSw(5)
351			& + eosJMDCKSw(6)*t
352			& + eosJMDCKSw(7)*t2
353			& )
354
355			C of sea water at p
356			dKdSPres =
357			& p*( eosJMDCKP(5)
358			& + eosJMDCKP(6)*t
359			& + eosJMDCKP(7)*t2
360			& )
361			& + s3o2peosJMDCKP(8)
362			& + pp( eosJMDCKP(12)
363			& + eosJMDCKP(13)*t
364			& + eosJMDCKP(14)*t2
365			& )
366
367			dKdS = dKdSSalt + dKdSPres
368
369			betaloc(i,j) =
370			& ( bulkmod(i,j)*2drhoP0dS
371			& - bulkmod(i,j)pdrhoP0dS
372			& - rhoP0(i,j)pdKdS )
373			& /( bulkmod(i,j) - p )**2
374
375
376			end do
377			end do
378	adcroft	1.1	else
379	heimbach	1.4	write(,) 'FIND_BETA: eqn = ',eqn
380	adcroft	1.1	stop 'FIND_BETA: argument "eqn" has illegal value'
381			endif
382
383			return
384			end