C $Header: /u/gcmpack/MITgcm/model/src/find_alpha.F,v 1.21 2014/04/04 20:56:32 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
#define USE_FACTORIZED_POLY
C-- File find_alpha.F: Calculate expansion Coeff
C-- Contents
C-- o FIND_ALPHA
C-- o FIND_BETA
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: FIND_ALPHA
C !INTERFACE:
SUBROUTINE FIND_ALPHA (
I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
O alphaLoc,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | o SUBROUTINE FIND_ALPHA
C | Calculates [drho(S,T,z) / dT] of a horizontal slice
C *==========================================================*
C | k - is the Theta/Salt level
C | kRef - determines pressure reference level
C | (not used in 'LINEAR' mode)
C | alphaLoc - drho / dT (kg/m^3/C)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C k :: Level of Theta/Salt slice
C kRef :: Pressure reference level
c myThid :: thread number for this instance of the routine
INTEGER myThid
INTEGER bi,bj,iMin,iMax,jMin,jMax
INTEGER k
INTEGER kRef
_RL alphaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j
_RL refTemp,refSalt,tP,sP
_RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
_RL ct, sa, sqrtsa, p
_RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
_RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
_RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rhoP0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dnum_dtheta, dden_dtheta
_RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
c CALL LOOK_FOR_NEG_SALINITY(
c I iMin, iMax, jMin, jMax,
c U sFld,
c I k, bi, bj, myThid )
#endif
IF (equationOfState.EQ.'LINEAR') THEN
DO j=jMin,jMax
DO i=iMin,iMax
alphaLoc(i,j) = -rhonil * tAlpha
ENDDO
ENDDO
ELSEIF (equationOfState.EQ.'POLY3') THEN
refTemp=eosRefT(kRef)
refSalt=eosRefS(kRef)
DO j=jMin,jMax
DO i=iMin,iMax
tP=theta(i,j,k,bi,bj)-refTemp
sP=salt(i,j,k,bi,bj)-refSalt
#ifdef USE_FACTORIZED_POLY
alphaLoc(i,j) =
& ( eosC(6,kRef)
& *tP*3.
& +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
& )*tP
& +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
&
#else
alphaLoc(i,j) =
& eosC(1,kRef) +
& eosC(3,kRef)*tP*2. +
& eosC(4,kRef) *sP +
& eosC(6,kRef)*tP*tP*3. +
& eosC(7,kRef)*tP*2. *sP +
& eosC(8,kRef) *sP*sP
#endif
ENDDO
ENDDO
ELSEIF ( equationOfState(1:5).EQ.'JMD95'
& .OR. equationOfState.EQ.'UNESCO' ) THEN
C nonlinear equation of state in pressure coordinates
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHOP0(
I iMin, iMax, jMin, jMax,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoP0,
I myThid )
CALL FIND_BULKMOD(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O bulkMod,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
C abbreviations
t1 = theta(i,j,k,bi,bj)
t2 = t1*t1
t3 = t2*t1
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
s3o2 = 0. _d 0
#endif
s1 = salt(i,j,k,bi,bj)
IF ( s1 .GT. 0. _d 0 ) THEN
s3o2 = SQRT(s1*s1*s1)
ELSE
s1 = 0. _d 0
s3o2 = 0. _d 0
ENDIF
p1 = locPres(i,j)*SItoBar
p2 = p1*p1
C d(rho)/d(theta)
C of fresh water at p = 0
drhoP0dthetaFresh =
& eosJMDCFw(2)
& + 2.*eosJMDCFw(3)*t1
& + 3.*eosJMDCFw(4)*t2
& + 4.*eosJMDCFw(5)*t3
& + 5.*eosJMDCFw(6)*t3*t1
C of salt water at p = 0
drhoP0dthetaSalt =
& s1*(
& eosJMDCSw(2)
& + 2.*eosJMDCSw(3)*t1
& + 3.*eosJMDCSw(4)*t2
& + 4.*eosJMDCSw(5)*t3
& )
& + s3o2*(
& + eosJMDCSw(7)
& + 2.*eosJMDCSw(8)*t1
& )
C d(bulk modulus)/d(theta)
C of fresh water at p = 0
dKdthetaFresh =
& eosJMDCKFw(2)
& + 2.*eosJMDCKFw(3)*t1
& + 3.*eosJMDCKFw(4)*t2
& + 4.*eosJMDCKFw(5)*t3
C of sea water at p = 0
dKdthetaSalt =
& s1*( eosJMDCKSw(2)
& + 2.*eosJMDCKSw(3)*t1
& + 3.*eosJMDCKSw(4)*t2
& )
& + s3o2*( eosJMDCKSw(6)
& + 2.*eosJMDCKSw(7)*t1
& )
C of sea water at p
dKdthetaPres =
& p1*( eosJMDCKP(2)
& + 2.*eosJMDCKP(3)*t1
& + 3.*eosJMDCKP(4)*t2
& )
& + p1*s1*( eosJMDCKP(6)
& + 2.*eosJMDCKP(7)*t1
& )
& + p2*( eosJMDCKP(10)
& + 2.*eosJMDCKP(11)*t1
& )
& + p2*s1*( eosJMDCKP(13)
& + 2.*eosJMDCKP(14)*t1
& )
drhoP0dtheta = drhoP0dthetaFresh
& + drhoP0dthetaSalt
dKdtheta = dKdthetaFresh
& + dKdthetaSalt
& + dKdthetaPres
alphaLoc(i,j) =
& ( bulkmod(i,j)**2*drhoP0dtheta
& - bulkmod(i,j)*p1*drhoP0dtheta
& - rhoP0(i,j)*p1*dKdtheta )
& /( bulkmod(i,j) - p1 )**2
ENDDO
ENDDO
ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHONUM(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoLoc,
I myThid )
CALL FIND_RHODEN(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoDen,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
t1 = theta(i,j,k,bi,bj)
t2 = t1*t1
s1 = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
sp5 = 0. _d 0
#endif
IF ( s1 .GT. 0. _d 0 ) THEN
sp5 = SQRT(s1)
ELSE
s1 = 0. _d 0
sp5 = 0. _d 0
ENDIF
p1 = locPres(i,j)*SItodBar
p1t1 = p1*t1
dnum_dtheta = eosMDJWFnum(1)
& + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
& + eosMDJWFnum(5)*s1
& + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)
dden_dtheta = eosMDJWFden(1)
& + t1*(2.*eosMDJWFden(2)
& + t1*(3.*eosMDJWFden(3)
& + 4.*eosMDJWFden(4)*t1 ) )
& + s1*(eosMDJWFden(6)
& + t1*(3.*eosMDJWFden(7)*t1
& + 2.*eosMDJWFden(9)*sp5 ) )
& + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)
alphaLoc(i,j) = rhoDen(i,j)*(dnum_dtheta
& - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)
ENDDO
ENDDO
ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHOTEOS(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoLoc, rhoDen,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
ct = theta(i,j,k,bi,bj)
sa = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
sqrtsa = 0. _d 0
#endif
IF ( sa .GT. 0. _d 0 ) THEN
sqrtsa = SQRT(sa)
ELSE
sa = 0. _d 0
sqrtsa = 0. _d 0
ENDIF
p = locPres(i,j)*SItodBar
dnum_dtheta = teos(02)
& + ct*(2.*teos(03) + 3.*teos(04)*ct)
& + sa*(teos(06) + 2.*teos(07)*ct
& + sqrtsa*(teos(09) + ct*(2.*teos(10) + 3.*teos(11)*ct)))
& + p*( teos(13) + 2.*teos(14)*ct + sa*2.*teos(16)
& + p*(teos(18) + 2.*teos(19)*ct))
dden_dtheta = teos(22)
& + ct*(2.*teos(23) + ct*(3.*teos(24) + 4.*teos(25)*ct))
& + sa*(teos(27)
& + ct*(2.*teos(28) + ct*(3.*teos(29) + 4.*teos(30)*ct))
& + sqrtsa*(teos(32)
& + ct*(2.*teos(33) + ct*(3.*teos(34) + 4.*teos(35)*ct))))
& + p*(teos(38) + ct*(2.*teos(39) + 3.*teos(40)*ct)
& + teos(42)
& + p*(teos(44) + 2.*teos(45)*ct + teos(46)*sa
& + p*teos(48) ))
alphaLoc(i,j) = rhoDen(i,j)*(dnum_dtheta
& - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)
ENDDO
ENDDO
ELSE
WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
STOP 'FIND_ALPHA: "equationOfState" has illegal value'
ENDIF
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: FIND_BETA
C !INTERFACE:
SUBROUTINE FIND_BETA (
I bi, bj, iMin, iMax, jMin, jMax, k, kRef,
O betaLoc,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | o SUBROUTINE FIND_BETA
C | Calculates [drho(S,T,z) / dS] of a horizontal slice
C |==========================================================*
C | k - is the Theta/Salt level
C | kRef - determines pressure reference level
C | (not used in 'LINEAR' mode)
C | betaLoc - drho / dS (kg/m^3/PSU)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "EOS.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C k :: Level of Theta/Salt slice
C kRef :: Pressure reference level
c myThid :: thread number for this instance of the routine
INTEGER myThid
INTEGER bi,bj,iMin,iMax,jMin,jMax
INTEGER k
INTEGER kRef
_RL betaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j
_RL refTemp,refSalt,tP,sP
_RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
_RL ct, sa, sqrtsa, p
_RL drhoP0dS
_RL dKdS, dKdSSalt, dKdSPres
_RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rhoP0 (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dnum_dsalt, dden_dsalt
_RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP
#ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
c CALL LOOK_FOR_NEG_SALINITY(
c I iMin, iMax, jMin, jMax,
c U sFld,
c I k, bi, bj, myThid )
#endif
IF (equationOfState.EQ.'LINEAR') THEN
DO j=jMin,jMax
DO i=iMin,iMax
betaLoc(i,j) = rhonil * sBeta
ENDDO
ENDDO
ELSEIF (equationOfState.EQ.'POLY3') THEN
refTemp=eosRefT(kRef)
refSalt=eosRefS(kRef)
DO j=jMin,jMax
DO i=iMin,iMax
tP=theta(i,j,k,bi,bj)-refTemp
sP=salt(i,j,k,bi,bj)-refSalt
#ifdef USE_FACTORIZED_POLY
betaLoc(i,j) =
& ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
& + ( eosC(7,kRef)*tP
& +eosC(8,kRef)*sP*2. + eosC(4,kRef)
& )*tP
#else
betaLoc(i,j) =
& eosC(2,kRef) +
& eosC(4,kRef)*tP +
& eosC(5,kRef) *sP*2. +
& eosC(7,kRef)*tP*tP +
& eosC(8,kRef)*tP *sP*2. +
& eosC(9,kRef) *sP*sP*3.
#endif
ENDDO
ENDDO
ELSEIF ( equationOfState(1:5).EQ.'JMD95'
& .OR. equationOfState.EQ.'UNESCO' ) THEN
C nonlinear equation of state in pressure coordinates
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHOP0(
I iMin, iMax, jMin, jMax,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoP0,
I myThid )
CALL FIND_BULKMOD(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O bulkMod,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
C abbreviations
t1 = theta(i,j,k,bi,bj)
t2 = t1*t1
t3 = t2*t1
s1 = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
s3o2 = 0. _d 0
#endif
IF ( s1 .GT. 0. _d 0 ) THEN
s3o2 = 1.5*SQRT(s1)
ELSE
s1 = 0. _d 0
s3o2 = 0. _d 0
ENDIF
p1 = locPres(i,j)*SItoBar
C d(rho)/d(S)
C of fresh water at p = 0
drhoP0dS = 0. _d 0
C of salt water at p = 0
drhoP0dS = drhoP0dS
& + eosJMDCSw(1)
& + eosJMDCSw(2)*t1
& + eosJMDCSw(3)*t2
& + eosJMDCSw(4)*t3
& + eosJMDCSw(5)*t3*t1
& + s3o2*(
& eosJMDCSw(6)
& + eosJMDCSw(7)*t1
& + eosJMDCSw(8)*t2
& )
& + 2*eosJMDCSw(9)*s1
C d(bulk modulus)/d(S)
C of fresh water at p = 0
dKdS = 0. _d 0
C of sea water at p = 0
dKdSSalt =
& eosJMDCKSw(1)
& + eosJMDCKSw(2)*t1
& + eosJMDCKSw(3)*t2
& + eosJMDCKSw(4)*t3
& + s3o2*( eosJMDCKSw(5)
& + eosJMDCKSw(6)*t1
& + eosJMDCKSw(7)*t2
& )
C of sea water at p
dKdSPres =
& p1*( eosJMDCKP(5)
& + eosJMDCKP(6)*t1
& + eosJMDCKP(7)*t2
& )
& + s3o2*p1*eosJMDCKP(8)
& + p1*p1*( eosJMDCKP(12)
& + eosJMDCKP(13)*t1
& + eosJMDCKP(14)*t2
& )
dKdS = dKdSSalt + dKdSPres
betaLoc(i,j) =
& ( bulkmod(i,j)**2*drhoP0dS
& - bulkmod(i,j)*p1*drhoP0dS
& - rhoP0(i,j)*p1*dKdS )
& /( bulkmod(i,j) - p1 )**2
ENDDO
ENDDO
ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHONUM(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoLoc,
I myThid )
CALL FIND_RHODEN(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoDen,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
t1 = theta(i,j,k,bi,bj)
t2 = t1*t1
s1 = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
sp5 = 0. _d 0
#endif
IF ( s1 .GT. 0. _d 0 ) THEN
sp5 = SQRT(s1)
ELSE
s1 = 0. _d 0
sp5 = 0. _d 0
ENDIF
p1 = locPres(i,j)*SItodBar
p1t1 = p1*t1
dnum_dsalt = eosMDJWFnum(4)
& + eosMDJWFnum(5)*t1
& + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
dden_dsalt = eosMDJWFden(5)
& + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
& + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)
betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
& - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )
ENDDO
ENDDO
ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
CALL PRESSURE_FOR_EOS(
I bi, bj, iMin, iMax, jMin, jMax, kRef,
O locPres,
I myThid )
CALL FIND_RHOTEOS(
I iMin, iMax, jMin, jMax, locPres,
I theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
O rhoLoc, rhoDen,
I myThid )
DO j=jMin,jMax
DO i=iMin,iMax
ct = theta(i,j,k,bi,bj)
sa = salt(i,j,k,bi,bj)
#if (defined ALLOW_AUTODIFF defined TARGET_NEC_SX)
C this line makes TAF create vectorizable code
sqrtsa = 0. _d 0
#endif
IF ( sa .GT. 0. _d 0 ) THEN
sqrtsa = SQRT(sa)
ELSE
sa = 0. _d 0
sqrtsa = 0. _d 0
ENDIF
p = locPres(i,j)*SItodBar
dnum_dsalt = teos(05) + ct*(teos(06) + teos(07)*ct)
& + 1.5*sqrtsa*(teos(08)
& + ct*(teos(09) + ct*(teos(10) + teos(11)*ct)))
& + p*(teos(15) + teos(16)*ct + p*teos(20))
dden_dsalt = teos(26)
& + ct*(teos(27) + ct*(teos(28) + ct*(teos(29) + teos(30)*ct)))
& + 2.*teos(36)*sa
& + 1.5*sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
& + ct*(teos(34) + teos(35)*ct))))
& + p*(teos(41) + teos(42)*ct + p*teos(46))
betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
& - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )
ENDDO
ENDDO
ELSE
WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
STOP 'FIND_BETA: "equationOfState" has illegal value'
ENDIF
RETURN
END