C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.15 2003/10/27 22:32:55 heimbach 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