model/src/find_rho.F

1eadaea85b Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
e34b4d724c Mart*0003 #ifdef TARGET_NEC_SX
                0004 C     make sure that the factorized EOS is used on NEC vector computers
                0005 # define USE_FACTORIZED_EOS
                0006 #endif
                0007 C     this was the default, so we keep it that way
9835224081 Alis*0008 #define USE_FACTORIZED_POLY
924557e60a Chri*0009
de7f345328 Jean*0010 C--  File find_rho.F: Routines to compute density
                0011 C--   Contents
                0012 C--   o FIND_RHO_2D
                0013 C--   o FIND_RHOP0
                0014 C--   o FIND_BULKMOD
                0015 C--   o FIND_RHONUM
                0016 C--   o FIND_RHODEN
701e10a905 Mart*0017 C--   o FIND_RHOTEOS
22cb42e69a Jean*0018 C--   o FIND_RHO_SCALAR: in-situ density for individual points
de7f345328 Jean*0019 C--   o LOOK_FOR_NEG_SALINITY
                0020
                0021 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9366854e02 Chri*0022 CBOP
de7f345328 Jean*0023 C     !ROUTINE: FIND_RHO_2D
                0024 C     !INTERFACE:
                0025       SUBROUTINE FIND_RHO_2D(
                0026      I                iMin, iMax, jMin, jMax, kRef,
                0027      I                tFld, sFld,
d4493909b3 Jean*0028      O                rhoLoc,
de7f345328 Jean*0029      I                k, bi, bj, myThid )
                0030
                0031 C     !DESCRIPTION: \bv
                0032 C     *==========================================================*
                0033 C     | o SUBROUTINE FIND_RHO_2D
                0034 C     |   Calculates [rho(S,T,z)-rhoConst] of a 2-D slice
                0035 C     *==========================================================*
                0036 C     | kRef - determines pressure reference level
                0037 C     |        (not used in 'LINEAR' mode)
                0038 C     | Note:  k is not used ; keep it for debugging.
                0039 C     *==========================================================*
                0040 C     \ev
                0041
                0042 C     !USES:
                0043       IMPLICIT NONE
                0044 C     == Global variables ==
                0045 #include "SIZE.h"
                0046 #include "EEPARAMS.h"
                0047 #include "PARAMS.h"
                0048 #include "EOS.h"
                0049
                0050 C     !INPUT/OUTPUT PARAMETERS:
                0051 C     == Routine arguments ==
                0052 C     k    :: Level of Theta/Salt slice
                0053 C     kRef :: Pressure reference level
                0054       INTEGER iMin,iMax,jMin,jMax
                0055       INTEGER kRef
                0056       _RL tFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       _RL sFld  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0058       _RL rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
de7f345328 Jean*0059       INTEGER k, bi, bj
                0060       INTEGER myThid
                0061
                0062 C     !LOCAL VARIABLES:
                0063 C     == Local variables ==
09f39a4963 Jean*0064       INTEGER i,j
533b1afce7 Jean*0065       _RL refTemp, refSalt, sigRef, tP, sP, deltaSig, dRho
                0066       _RL dp0, oneMinusKap, facPres, theta_v
1eadaea85b Jean*0067       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0068       _RL rhoP0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0069       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0070       _RL rhoNum (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0071       _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
bc6d61d764 Jean*0072       CHARACTER*(MAX_LEN_MBUF) msgBuf
9366854e02 Chri*0073 CEOP
924557e60a Chri*0074
02d90fb24c Jean*0075 #ifdef ALLOW_AUTODIFF
8adf9f02ba Patr*0076       DO j=1-OLy,sNy+OLy
                0077        DO i=1-OLx,sNx+OLx
d4493909b3 Jean*0078         rhoLoc(i,j)  = 0. _d 0
a37a13034c Mart*0079         rhoP0(i,j)   = 0. _d 0
                0080         bulkMod(i,j) = 0. _d 0
8adf9f02ba Patr*0081        ENDDO
                0082       ENDDO
                0083 #endif
                0084
59e59b79d8 Mart*0085 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
de7f345328 Jean*0086       CALL LOOK_FOR_NEG_SALINITY(
                0087      I          iMin, iMax, jMin, jMax,
                0088      U          sFld,
                0089      I          k, bi, bj, myThid )
59e59b79d8 Mart*0090 #endif
                0091
6ac14281aa Jean*0092       IF (equationOfState.EQ.'LINEAR') THEN
047cf42561 Alis*0093
09f39a4963 Jean*0094 C ***NOTE***
                0095 C In the linear EOS, to make the static stability calculation meaningful
bc6d61d764 Jean*0096 C we use reference temp & salt from level kRef ;
09f39a4963 Jean*0097 C **********
                0098        refTemp=tRef(kRef)
                0099        refSalt=sRef(kRef)
                0100
e305438401 Mart*0101        dRho = rhoNil-rhoConst
                0102
6ac14281aa Jean*0103        DO j=jMin,jMax
                0104         DO i=iMin,iMax
d4493909b3 Jean*0105          rhoLoc(i,j)=rhoNil*(
de7f345328 Jean*0106      &     sBeta*(sFld(i,j)-refSalt)
                0107      &   -tAlpha*(tFld(i,j)-refTemp) )
e305438401 Mart*0108      &        + dRho
6ac14281aa Jean*0109         ENDDO
                0110        ENDDO
de7f345328 Jean*0111
6ac14281aa Jean*0112       ELSEIF (equationOfState.EQ.'POLY3') THEN
047cf42561 Alis*0113
                0114        refTemp=eosRefT(kRef)
                0115        refSalt=eosRefS(kRef)
e305438401 Mart*0116        sigRef=eosSig0(kRef) + (1000.-rhoConst)
047cf42561 Alis*0117
6ac14281aa Jean*0118        DO j=jMin,jMax
                0119         DO i=iMin,iMax
de7f345328 Jean*0120          tP=tFld(i,j)-refTemp
                0121          sP=sFld(i,j)-refSalt
9835224081 Alis*0122 #ifdef USE_FACTORIZED_POLY
047cf42561 Alis*0123          deltaSig=
9835224081 Alis*0124      &    (( eosC(9,kRef)*sP + eosC(5,kRef) )*sP + eosC(2,kRef) )*sP
                0125      &   + ( ( eosC(6,kRef)
                0126      &         *tP
                0127      &        +eosC(7,kRef)*sP + eosC(3,kRef)
                0128      &       )*tP
                0129      &      +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
                0130      &     )*tP
                0131 #else
                0132          deltaSig=
                0133      &     eosC(1,kRef)*tP
                0134      &    +eosC(2,kRef)         *sP
                0135      &    +eosC(3,kRef)*tP*tP
                0136      &    +eosC(4,kRef)*tP      *sP
                0137      &    +eosC(5,kRef)         *sP*sP
                0138      &    +eosC(6,kRef)*tP*tP*tP
                0139      &    +eosC(7,kRef)*tP*tP   *sP
                0140      &    +eosC(8,kRef)*tP      *sP*sP
                0141      &    +eosC(9,kRef)         *sP*sP*sP
                0142 #endif
d4493909b3 Jean*0143          rhoLoc(i,j)=sigRef+deltaSig
6ac14281aa Jean*0144         ENDDO
                0145        ENDDO
047cf42561 Alis*0146
6ac14281aa Jean*0147       ELSEIF ( equationOfState(1:5).EQ.'JMD95'
                0148      &      .OR. equationOfState.EQ.'UNESCO' ) THEN
a37a13034c Mart*0149 C     nonlinear equation of state in pressure coordinates
533b1afce7 Jean*0150          dp0 = surf_pRef - eosRefP0
a37a13034c Mart*0151
6ac14281aa Jean*0152          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0153      I             bi, bj, iMin, iMax, jMin, jMax, kRef, dp0,
de7f345328 Jean*0154      O             locPres,
                0155      I             myThid )
6ac14281aa Jean*0156
a37a13034c Mart*0157          CALL FIND_RHOP0(
de7f345328 Jean*0158      I             iMin, iMax, jMin, jMax,
                0159      I             tFld, sFld,
                0160      O             rhoP0,
                0161      I             myThid )
                0162
a37a13034c Mart*0163          CALL FIND_BULKMOD(
de7f345328 Jean*0164      I             iMin, iMax, jMin, jMax,
                0165      I             locPres, tFld, sFld,
                0166      O             bulkMod,
                0167      I             myThid )
                0168
d4493909b3 Jean*0169 c#ifdef ALLOW_AUTODIFF_TAMC
04522666de Ed H*0170 cph can not DO storing here since find_rho is called multiple times;
7103bd8015 Patr*0171 cph additional recomp. should be acceptable
d4493909b3 Jean*0172 c#endif /* ALLOW_AUTODIFF_TAMC */
6ac14281aa Jean*0173          DO j=jMin,jMax
de7f345328 Jean*0174           DO i=iMin,iMax
a37a13034c Mart*0175
                0176 C     density of sea water at pressure p
d4493909b3 Jean*0177             rhoLoc(i,j) = rhoP0(i,j)
de7f345328 Jean*0178      &              /(1. _d 0 -
                0179      &              locPres(i,j)*SItoBar/bulkMod(i,j) )
e305438401 Mart*0180      &              - rhoConst
a37a13034c Mart*0181
de7f345328 Jean*0182           ENDDO
6ac14281aa Jean*0183          ENDDO
a37a13034c Mart*0184
6ac14281aa Jean*0185       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
533b1afce7 Jean*0186          dp0 = surf_pRef - eosRefP0
6ac14281aa Jean*0187
                0188          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0189      I             bi, bj, iMin, iMax, jMin, jMax, kRef, dp0,
de7f345328 Jean*0190      O             locPres,
                0191      I             myThid )
                0192
                0193          CALL FIND_RHONUM(
                0194      I             iMin, iMax, jMin, jMax,
                0195      I             locPres, tFld, sFld,
                0196      O             rhoNum,
                0197      I             myThid )
                0198
                0199          CALL FIND_RHODEN(
                0200      I             iMin, iMax, jMin, jMax,
                0201      I             locPres, tFld, sFld,
                0202      O             rhoDen,
                0203      I             myThid )
d4493909b3 Jean*0204
                0205 c#ifdef ALLOW_AUTODIFF_TAMC
04522666de Ed H*0206 cph can not DO storing here since find_rho is called multiple times;
7103bd8015 Patr*0207 cph additional recomp. should be acceptable
d4493909b3 Jean*0208 c#endif /* ALLOW_AUTODIFF_TAMC */
6ac14281aa Jean*0209          DO j=jMin,jMax
de7f345328 Jean*0210           DO i=iMin,iMax
d4493909b3 Jean*0211             rhoLoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
de7f345328 Jean*0212           ENDDO
6ac14281aa Jean*0213          ENDDO
31a5949e20 Mart*0214
ac1c5b24e7 Mart*0215       ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
533b1afce7 Jean*0216          dp0 = surf_pRef - eosRefP0
ac1c5b24e7 Mart*0217
                0218          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0219      I             bi, bj, iMin, iMax, jMin, jMax, kRef, dp0,
ac1c5b24e7 Mart*0220      O             locPres,
                0221      I             myThid )
                0222
                0223          CALL FIND_RHOTEOS(
                0224      I             iMin, iMax, jMin, jMax,
                0225      I             locPres, tFld, sFld,
                0226      O             rhoNum, rhoDen,
                0227      I             myThid )
                0228
                0229 c#ifdef ALLOW_AUTODIFF_TAMC
                0230 cph can not DO storing here since find_rho is called multiple times;
                0231 cph additional recomp. should be acceptable
                0232 c#endif /* ALLOW_AUTODIFF_TAMC */
                0233          DO j=jMin,jMax
                0234           DO i=iMin,iMax
                0235             rhoLoc(i,j) = rhoNum(i,j)*rhoDen(i,j) - rhoConst
                0236           ENDDO
                0237          ENDDO
                0238
6ac14281aa Jean*0239       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
4ee1a938b6 Jean*0240
                0241          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0242      I             bi, bj, iMin, iMax, jMin, jMax, kRef, zeroRL,
4ee1a938b6 Jean*0243      O             locPres,
                0244      I             myThid )
                0245
                0246          oneMinusKap = oneRL - atm_kappa
                0247 C        rho = P/(R*T_v) = Po/(R*theta_v)*(P/Po)^(1-kappa)
                0248          DO j=jMin,jMax
                0249           DO i=iMin,iMax
                0250            IF ( locPres(i,j).GT.zeroRL .AND. tFld(i,j).GT.zeroRL ) THEN
                0251             facPres = ( locPres(i,j)/atm_Po )**oneMinusKap
                0252             theta_v = tFld(i,j)*( sFld(i,j)*atm_Rq + oneRL )
                0253             rhoLoc(i,j) = atm_Po*facPres
                0254      &                 /( atm_Rd*theta_v ) - rhoConst
                0255            ELSE
                0256             rhoLoc(i,j) = zeroRL
                0257            ENDIF
                0258           ENDDO
                0259          ENDDO
                0260
6ac14281aa Jean*0261       ELSE
bc6d61d764 Jean*0262        WRITE(msgBuf,'(3a)')
de7f345328 Jean*0263      &      ' FIND_RHO_2D: equationOfState = "',equationOfState,'"'
bc6d61d764 Jean*0264        CALL PRINT_ERROR( msgBuf, myThid )
de7f345328 Jean*0265        STOP 'ABNORMAL END: S/R FIND_RHO_2D'
6ac14281aa Jean*0266       ENDIF
924557e60a Chri*0267
de7f345328 Jean*0268       RETURN
6ac14281aa Jean*0269       END
a37a13034c Mart*0270
de7f345328 Jean*0271 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
a37a13034c Mart*0272 CBOP
                0273 C     !ROUTINE: FIND_RHOP0
                0274 C     !INTERFACE:
6ac14281aa Jean*0275       SUBROUTINE FIND_RHOP0(
de7f345328 Jean*0276      I                iMin, iMax, jMin, jMax,
                0277      I                tFld, sFld,
                0278      O                rhoP0,
                0279      I                myThid )
a37a13034c Mart*0280
                0281 C     !DESCRIPTION: \bv
                0282 C     *==========================================================*
                0283 C     | o SUBROUTINE FIND_RHOP0
de7f345328 Jean*0284 C     |   Calculates rho(S,T,0) of a slice
a37a13034c Mart*0285 C     *==========================================================*
                0286 C     *==========================================================*
                0287 C     \ev
                0288
                0289 C     !USES:
6ac14281aa Jean*0290       IMPLICIT NONE
a37a13034c Mart*0291 C     == Global variables ==
                0292 #include "SIZE.h"
                0293 #include "EEPARAMS.h"
                0294 #include "PARAMS.h"
                0295 #include "EOS.h"
                0296
                0297 C     !INPUT/OUTPUT PARAMETERS:
                0298 C     == Routine arguments ==
de7f345328 Jean*0299       INTEGER iMin,iMax,jMin,jMax
                0300       _RL tFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0301       _RL sFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0302       _RL rhoP0(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6ac14281aa Jean*0303       INTEGER myThid
a37a13034c Mart*0304
                0305 C     !LOCAL VARIABLES:
                0306 C     == Local variables ==
6ac14281aa Jean*0307       INTEGER i,j
5ef7d4a8ca Jean*0308       _RL rfresh, rsalt
a37a13034c Mart*0309       _RL t, t2, t3, t4, s, s3o2
e34b4d724c Mart*0310 #ifdef USE_FACTORIZED_EOS
                0311       _RL eosF1, eosF2, eosF3, eosF4, eosF5, eosF6
                0312       _RL eosS1, eosS2, eosS3, eosS4, eosS5, eosS6
                0313       _RL eosS7, eosS8, eosS9
                0314 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0315 CEOP
e34b4d724c Mart*0316 #ifdef USE_FACTORIZED_EOS
1eadaea85b Jean*0317       eosF1 = eosJMDCFw(1)
e34b4d724c Mart*0318       eosF2 = eosJMDCFw(2)
                0319       eosF3 = eosJMDCFw(3)
                0320       eosF4 = eosJMDCFw(4)
                0321       eosF5 = eosJMDCFw(5)
                0322       eosF6 = eosJMDCFw(6)
                0323       eosS1 = eosJMDCSw(1)
                0324       eosS2 = eosJMDCSw(2)
                0325       eosS3 = eosJMDCSw(3)
                0326       eosS4 = eosJMDCSw(4)
                0327       eosS5 = eosJMDCSw(5)
                0328       eosS6 = eosJMDCSw(6)
                0329       eosS7 = eosJMDCSw(7)
                0330       eosS8 = eosJMDCSw(8)
                0331       eosS9 = eosJMDCSw(9)
                0332 #endif /* USE_FACTORIZED_EOS */
1eadaea85b Jean*0333
6ac14281aa Jean*0334       DO j=jMin,jMax
                0335          DO i=iMin,iMax
a37a13034c Mart*0336 C     abbreviations
de7f345328 Jean*0337             t  = tFld(i,j)
a37a13034c Mart*0338             t2 = t*t
                0339             t3 = t2*t
                0340             t4 = t3*t
                0341
02d90fb24c Jean*0342 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0343 C     this line makes TAF create vectorizable code
                0344             s3o2 = 0. _d 0
                0345 #endif
de7f345328 Jean*0346             s  = sFld(i,j)
59e581f4e0 Jean*0347            IF ( s .GT. 0. _d 0 ) THEN
6ac14281aa Jean*0348             s3o2 = s*SQRT(s)
59e581f4e0 Jean*0349            ELSE
                0350             s    = 0. _d 0
                0351             s3o2 = 0. _d 0
                0352            ENDIF
de7f345328 Jean*0353
a37a13034c Mart*0354 C     density of freshwater at the surface
e34b4d724c Mart*0355 #ifdef USE_FACTORIZED_EOS
1eadaea85b Jean*0356             rfresh =
e34b4d724c Mart*0357      &                    eosF1 +
                0358      &               t* ( eosF2 +
                0359      &               t* ( eosF3 +
                0360      &               t* ( eosF4 +
                0361      &               t* ( eosF5 +
                0362      &               t*   eosF6  ))))
                0363 #else
1eadaea85b Jean*0364             rfresh =
de7f345328 Jean*0365      &             eosJMDCFw(1)
a37a13034c Mart*0366      &           + eosJMDCFw(2)*t
                0367      &           + eosJMDCFw(3)*t2
                0368      &           + eosJMDCFw(4)*t3
                0369      &           + eosJMDCFw(5)*t4
                0370      &           + eosJMDCFw(6)*t4*t
e34b4d724c Mart*0371 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0372 C     density of sea water at the surface
e34b4d724c Mart*0373 #ifdef USE_FACTORIZED_EOS
                0374             rsalt =  s * ( eosS1 +
                0375      &               t * ( eosS2 +
                0376      &               t * ( eosS3 +
                0377      &               t * ( eosS4 +
                0378      &               t *   eosS5  ))))
                0379      &          + s3o2 * ( eosS6 +
1eadaea85b Jean*0380      &               t * ( eosS7 +
e34b4d724c Mart*0381      &               t *   eosS8  ))
                0382      &          + s*s  *   eosS9
                0383 #else
de7f345328 Jean*0384             rsalt =
a37a13034c Mart*0385      &         s*(
                0386      &             eosJMDCSw(1)
                0387      &           + eosJMDCSw(2)*t
                0388      &           + eosJMDCSw(3)*t2
                0389      &           + eosJMDCSw(4)*t3
                0390      &           + eosJMDCSw(5)*t4
                0391      &           )
                0392      &       + s3o2*(
                0393      &             eosJMDCSw(6)
                0394      &           + eosJMDCSw(7)*t
                0395      &           + eosJMDCSw(8)*t2
                0396      &           )
                0397      &           + eosJMDCSw(9)*s*s
e34b4d724c Mart*0398 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0399
                0400             rhoP0(i,j) = rfresh + rsalt
                0401
6ac14281aa Jean*0402          ENDDO
                0403       ENDDO
a37a13034c Mart*0404
de7f345328 Jean*0405       RETURN
6ac14281aa Jean*0406       END
a37a13034c Mart*0407
22cb42e69a Jean*0408 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0409 CBOP
a37a13034c Mart*0410 C     !ROUTINE: FIND_BULKMOD
                0411 C     !INTERFACE:
6ac14281aa Jean*0412       SUBROUTINE FIND_BULKMOD(
de7f345328 Jean*0413      I                iMin, iMax, jMin, jMax,
                0414      I                locPres, tFld, sFld,
                0415      O                bulkMod,
                0416      I                myThid )
a37a13034c Mart*0417
                0418 C     !DESCRIPTION: \bv
                0419 C     *==========================================================*
de7f345328 Jean*0420 C     | o SUBROUTINE FIND_BULKMOD
a37a13034c Mart*0421 C     |   Calculates the secant bulk modulus K(S,T,p) of a slice
                0422 C     *==========================================================*
6ac14281aa Jean*0423 C     | k    - is the level of Theta/Salt slice
a37a13034c Mart*0424 C     *==========================================================*
                0425 C     \ev
                0426
                0427 C     !USES:
6ac14281aa Jean*0428       IMPLICIT NONE
a37a13034c Mart*0429 C     == Global variables ==
                0430 #include "SIZE.h"
                0431 #include "EEPARAMS.h"
                0432 #include "PARAMS.h"
                0433 #include "EOS.h"
                0434
                0435 C     !INPUT/OUTPUT PARAMETERS:
                0436 C     == Routine arguments ==
de7f345328 Jean*0437       INTEGER iMin,iMax,jMin,jMax
1eadaea85b Jean*0438       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
de7f345328 Jean*0439       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0440       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0441       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6ac14281aa Jean*0442       INTEGER myThid
a37a13034c Mart*0443
                0444 C     !LOCAL VARIABLES:
                0445 C     == Local variables ==
6ac14281aa Jean*0446       INTEGER i,j
5ef7d4a8ca Jean*0447       _RL bMfresh, bMsalt, bMpres
a37a13034c Mart*0448       _RL t, t2, t3, t4, s, s3o2, p, p2
e34b4d724c Mart*0449 #ifdef USE_FACTORIZED_EOS
                0450       _RL eosF1, eosF2, eosF3, eosF4, eosF5
                0451       _RL eosS1, eosS2, eosS3, eosS4, eosS5, eosS6, eosS7
                0452       _RL eosP1, eosP2, eosP3, eosP4, eosP5, eosP6, eosP7
                0453       _RL eosP8, eosP9, eosP10, eosP11, eosP12, eosP13, eosP14
                0454 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0455 CEOP
e34b4d724c Mart*0456 #ifdef USE_FACTORIZED_EOS
                0457       eosF1 = eosJMDCKFw(1)
                0458       eosF2 = eosJMDCKFw(2)
                0459       eosF3 = eosJMDCKFw(3)
                0460       eosF4 = eosJMDCKFw(4)
                0461       eosF5 = eosJMDCKFw(5)
                0462
                0463       eosS1 = eosJMDCKSw(1)
                0464       eosS2 = eosJMDCKSw(2)
                0465       eosS3 = eosJMDCKSw(3)
                0466       eosS4 = eosJMDCKSw(4)
                0467       eosS5 = eosJMDCKSw(5)
                0468       eosS6 = eosJMDCKSw(6)
                0469       eosS7 = eosJMDCKSw(7)
                0470
                0471       eosP1  =eosJMDCKP(1)
                0472       eosP2  =eosJMDCKP(2)
                0473       eosP3  =eosJMDCKP(3)
                0474       eosP4  =eosJMDCKP(4)
                0475       eosP5  =eosJMDCKP(5)
                0476       eosP6  =eosJMDCKP(6)
                0477       eosP7  =eosJMDCKP(7)
                0478       eosP8  =eosJMDCKP(8)
                0479       eosP9  =eosJMDCKP(9)
                0480       eosP10 =eosJMDCKP(10)
                0481       eosP11 =eosJMDCKP(11)
                0482       eosP12 =eosJMDCKP(12)
                0483       eosP13 =eosJMDCKP(13)
                0484       eosP14 =eosJMDCKP(14)
                0485 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0486
6ac14281aa Jean*0487       DO j=jMin,jMax
                0488          DO i=iMin,iMax
a37a13034c Mart*0489 C     abbreviations
de7f345328 Jean*0490             t  = tFld(i,j)
a37a13034c Mart*0491             t2 = t*t
                0492             t3 = t2*t
                0493             t4 = t3*t
                0494
02d90fb24c Jean*0495 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0496 C     this line makes TAF create vectorizable code
                0497             s3o2 = 0. _d 0
                0498 #endif
de7f345328 Jean*0499             s  = sFld(i,j)
59e581f4e0 Jean*0500           IF ( s .GT. 0. _d 0 ) THEN
6ac14281aa Jean*0501             s3o2 = s*SQRT(s)
59e581f4e0 Jean*0502           ELSE
                0503             s    = 0. _d 0
                0504             s3o2 = 0. _d 0
                0505           ENDIF
a37a13034c Mart*0506 C
6ac14281aa Jean*0507             p = locPres(i,j)*SItoBar
a37a13034c Mart*0508             p2 = p*p
                0509 C     secant bulk modulus of fresh water at the surface
e34b4d724c Mart*0510 #ifdef USE_FACTORIZED_EOS
                0511             bMfresh =       eosF1 +
                0512      &                t * ( eosF2 +
                0513      &                t * ( eosF3 +
                0514      &                t * ( eosF4 +
                0515      &                t *   eosF5  )))
                0516 #else
de7f345328 Jean*0517             bMfresh =
a37a13034c Mart*0518      &             eosJMDCKFw(1)
                0519      &           + eosJMDCKFw(2)*t
                0520      &           + eosJMDCKFw(3)*t2
                0521      &           + eosJMDCKFw(4)*t3
                0522      &           + eosJMDCKFw(5)*t4
e34b4d724c Mart*0523 #endif /* USE_FACTORIZED_EOS */
                0524 C     secant bulk modulus of sea water at the surface
                0525 #ifdef USE_FACTORIZED_EOS
                0526             bMsalt =  s * ( eosS1 +
                0527      &                t * ( eosS2 +
                0528      &                t * ( eosS3 +
                0529      &                t *   eosS4  )))
                0530      &            + s3o2* ( eosS5 +
                0531      &                t * ( eosS6 +
                0532      &                t *   eosS7  ))
                0533 C     secant bulk modulus of sea water at pressure p
                0534             bMpres =  p * ( eosP1 +
                0535      &                t * ( eosP2 +
                0536      &                t * ( eosP3 +
                0537      &                t *   eosP4  )))
                0538      &             + p*s* ( eosP5 +
                0539      &                t * ( eosP6 +
                0540      &                t *   eosP7  ))
                0541      &             + p*s3o2*eosP8
                0542      &             + p2 * ( eosP9 +
                0543      &                t * ( eosP10 +
                0544      &                t *   eosP11  ))
                0545      &            + p2*s* ( eosP12 +
                0546      &                t * ( eosP13 +
                0547      &                t *   eosP14 ))
                0548 #else
a37a13034c Mart*0549 C     secant bulk modulus of sea water at the surface
                0550             bMsalt =
                0551      &         s*( eosJMDCKSw(1)
                0552      &           + eosJMDCKSw(2)*t
                0553      &           + eosJMDCKSw(3)*t2
                0554      &           + eosJMDCKSw(4)*t3
                0555      &           )
                0556      &    + s3o2*( eosJMDCKSw(5)
                0557      &           + eosJMDCKSw(6)*t
                0558      &           + eosJMDCKSw(7)*t2
                0559      &           )
                0560 C     secant bulk modulus of sea water at pressure p
de7f345328 Jean*0561             bMpres =
a37a13034c Mart*0562      &         p*( eosJMDCKP(1)
                0563      &           + eosJMDCKP(2)*t
                0564      &           + eosJMDCKP(3)*t2
                0565      &           + eosJMDCKP(4)*t3
                0566      &           )
                0567      &     + p*s*( eosJMDCKP(5)
                0568      &           + eosJMDCKP(6)*t
                0569      &           + eosJMDCKP(7)*t2
                0570      &           )
                0571      &      + p*s3o2*eosJMDCKP(8)
                0572      &      + p2*( eosJMDCKP(9)
                0573      &           + eosJMDCKP(10)*t
                0574      &           + eosJMDCKP(11)*t2
                0575      &           )
                0576      &    + p2*s*( eosJMDCKP(12)
                0577      &           + eosJMDCKP(13)*t
                0578      &           + eosJMDCKP(14)*t2
                0579      &           )
e34b4d724c Mart*0580 #endif /* USE_FACTORIZED_EOS */
a37a13034c Mart*0581
                0582             bulkMod(i,j) = bMfresh + bMsalt + bMpres
                0583
6ac14281aa Jean*0584          ENDDO
                0585       ENDDO
a37a13034c Mart*0586
de7f345328 Jean*0587       RETURN
6ac14281aa Jean*0588       END
a37a13034c Mart*0589
22cb42e69a Jean*0590 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
31a5949e20 Mart*0591 CBOP
                0592 C     !ROUTINE: FIND_RHONUM
                0593 C     !INTERFACE:
6ac14281aa Jean*0594       SUBROUTINE FIND_RHONUM(
de7f345328 Jean*0595      I                iMin, iMax, jMin, jMax,
                0596      I                locPres, tFld, sFld,
                0597      O                rhoNum,
                0598      I                myThid )
31a5949e20 Mart*0599
                0600 C     !DESCRIPTION: \bv
                0601 C     *==========================================================*
                0602 C     | o SUBROUTINE FIND_RHONUM
                0603 C     |   Calculates the numerator of the McDougall et al.
                0604 C     |   equation of state
                0605 C     |   - the code is more or less a copy of MOM4
                0606 C     *==========================================================*
                0607 C     *==========================================================*
                0608 C     \ev
                0609
                0610 C     !USES:
6ac14281aa Jean*0611       IMPLICIT NONE
31a5949e20 Mart*0612 C     == Global variables ==
                0613 #include "SIZE.h"
                0614 #include "EEPARAMS.h"
                0615 #include "PARAMS.h"
                0616 #include "EOS.h"
                0617
                0618 C     !INPUT/OUTPUT PARAMETERS:
                0619 C     == Routine arguments ==
de7f345328 Jean*0620       INTEGER iMin,iMax,jMin,jMax
1eadaea85b Jean*0621       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
de7f345328 Jean*0622       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0623       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0624       _RL rhoNum (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6ac14281aa Jean*0625       INTEGER myThid
31a5949e20 Mart*0626
                0627 C     !LOCAL VARIABLES:
                0628 C     == Local variables ==
6ac14281aa Jean*0629       INTEGER i,j
31a5949e20 Mart*0630       _RL t1, t2, s1, p1
                0631 CEOP
6ac14281aa Jean*0632       DO j=jMin,jMax
                0633          DO i=iMin,iMax
31a5949e20 Mart*0634 C     abbreviations
de7f345328 Jean*0635             t1  = tFld(i,j)
31a5949e20 Mart*0636             t2 = t1*t1
de7f345328 Jean*0637             s1  = sFld(i,j)
                0638
6ac14281aa Jean*0639             p1   = locPres(i,j)*SItodBar
de7f345328 Jean*0640
31a5949e20 Mart*0641             rhoNum(i,j) = eosMDJWFnum(0)
de7f345328 Jean*0642      &           + t1*(eosMDJWFnum(1)
                0643      &           +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
31a5949e20 Mart*0644      &           + s1*(eosMDJWFnum(4)
de7f345328 Jean*0645      &           +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)
                0646      &           + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
                0647      &           +     eosMDJWFnum(9)*s1
31a5949e20 Mart*0648      &           +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
                0649
6ac14281aa Jean*0650          ENDDO
                0651       ENDDO
31a5949e20 Mart*0652
de7f345328 Jean*0653       RETURN
                0654       END
31a5949e20 Mart*0655
22cb42e69a Jean*0656 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
31a5949e20 Mart*0657 CBOP
                0658 C     !ROUTINE: FIND_RHODEN
                0659 C     !INTERFACE:
6ac14281aa Jean*0660       SUBROUTINE FIND_RHODEN(
de7f345328 Jean*0661      I                iMin, iMax, jMin, jMax,
                0662      I                locPres, tFld, sFld,
                0663      O                rhoDen,
                0664      I                myThid )
31a5949e20 Mart*0665
                0666 C     !DESCRIPTION: \bv
                0667 C     *==========================================================*
                0668 C     | o SUBROUTINE FIND_RHODEN
                0669 C     |   Calculates the denominator of the McDougall et al.
                0670 C     |   equation of state
                0671 C     |   - the code is more or less a copy of MOM4
                0672 C     *==========================================================*
                0673 C     *==========================================================*
                0674 C     \ev
                0675
                0676 C     !USES:
6ac14281aa Jean*0677       IMPLICIT NONE
31a5949e20 Mart*0678 C     == Global variables ==
                0679 #include "SIZE.h"
                0680 #include "EEPARAMS.h"
                0681 #include "PARAMS.h"
                0682 #include "EOS.h"
                0683
                0684 C     !INPUT/OUTPUT PARAMETERS:
                0685 C     == Routine arguments ==
de7f345328 Jean*0686       INTEGER iMin,iMax,jMin,jMax
1eadaea85b Jean*0687       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
de7f345328 Jean*0688       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0689       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0690       _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
6ac14281aa Jean*0691       INTEGER myThid
31a5949e20 Mart*0692
                0693 C     !LOCAL VARIABLES:
                0694 C     == Local variables ==
6ac14281aa Jean*0695       INTEGER i,j
31a5949e20 Mart*0696       _RL t1, t2, s1, sp5, p1, p1t1
                0697       _RL den, epsln
                0698       parameter ( epsln = 0. _d 0 )
                0699 CEOP
6ac14281aa Jean*0700       DO j=jMin,jMax
                0701          DO i=iMin,iMax
31a5949e20 Mart*0702 C     abbreviations
de7f345328 Jean*0703             t1  = tFld(i,j)
31a5949e20 Mart*0704             t2 = t1*t1
de7f345328 Jean*0705             s1  = sFld(i,j)
02d90fb24c Jean*0706 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0707 C     this line makes TAF create vectorizable code
                0708             sp5 = 0. _d 0
                0709 #endif
59e581f4e0 Jean*0710            IF ( s1 .GT. 0. _d 0 ) THEN
de7f345328 Jean*0711             sp5 = SQRT(s1)
59e581f4e0 Jean*0712            ELSE
                0713             s1  = 0. _d 0
                0714             sp5 = 0. _d 0
                0715            ENDIF
de7f345328 Jean*0716
6ac14281aa Jean*0717             p1   = locPres(i,j)*SItodBar
31a5949e20 Mart*0718             p1t1 = p1*t1
de7f345328 Jean*0719
31a5949e20 Mart*0720             den = eosMDJWFden(0)
                0721      &           + t1*(eosMDJWFden(1)
                0722      &           +     t1*(eosMDJWFden(2)
                0723      &           +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
                0724      &           + s1*(eosMDJWFden(5)
de7f345328 Jean*0725      &           +     t1*(eosMDJWFden(6)
                0726      &           +         eosMDJWFden(7)*t2)
                0727      &           +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
31a5949e20 Mart*0728      &           + p1*(eosMDJWFden(10)
                0729      &           +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
de7f345328 Jean*0730
                0731             rhoDen(i,j) = 1.0/(epsln+den)
31a5949e20 Mart*0732
6ac14281aa Jean*0733          ENDDO
                0734       ENDDO
31a5949e20 Mart*0735
de7f345328 Jean*0736       RETURN
                0737       END
60c223928f Mart*0738
22cb42e69a Jean*0739 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0740 CBOP
ac1c5b24e7 Mart*0741 C     !ROUTINE: FIND_RHOTEOS
                0742 C     !INTERFACE:
                0743       SUBROUTINE FIND_RHOTEOS(
                0744      I                iMin, iMax, jMin, jMax,
                0745      I                locPres, tFld, sFld,
                0746      O                rhoNum, rhoDen,
                0747      I                myThid )
                0748
                0749 C     !DESCRIPTION: \bv
                0750 C     *==========================================================*
e34b4d724c Mart*0751 C     | o SUBROUTINE FIND_RHOTEOS
1eadaea85b Jean*0752 C     |   Calculates numerator and inverse of denominator of the
ac1c5b24e7 Mart*0753 C     |   TEOS-10 (McDougall et al. 2011) equation of state
                0754 C     |   - the code is more or less a copy the teos-10 toolbox
                0755 C     |     available at http://www.teos-10.org
                0756 C     *==========================================================*
                0757 C     \ev
                0758
                0759 C     !USES:
                0760       IMPLICIT NONE
                0761 C     == Global variables ==
                0762 #include "SIZE.h"
                0763 #include "EEPARAMS.h"
                0764 #include "PARAMS.h"
                0765 #include "EOS.h"
                0766
                0767 C     !INPUT/OUTPUT PARAMETERS:
                0768 C     == Routine arguments ==
                0769       INTEGER iMin,iMax,jMin,jMax
1eadaea85b Jean*0770       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ac1c5b24e7 Mart*0771       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0772       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1eadaea85b Jean*0773       _RL rhoNum (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0774       _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
ac1c5b24e7 Mart*0775       INTEGER myThid
                0776
                0777 C     !LOCAL VARIABLES:
                0778 C     == Local variables ==
                0779       INTEGER i,j
                0780       _RL ct, sa, sqrtsa, p
                0781       _RL den, epsln
                0782       parameter ( epsln = 0. _d 0 )
                0783 CEOP
                0784       DO j=jMin,jMax
                0785        DO i=iMin,iMax
                0786 C     abbreviations
                0787         ct      = tFld(i,j)
                0788         sa      = sFld(i,j)
02d90fb24c Jean*0789 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0790 C     this line makes TAF create vectorizable code
                0791         sqrtsa = 0. _d 0
                0792 #endif
ac1c5b24e7 Mart*0793         IF ( sa .GT. 0. _d 0 ) THEN
                0794          sqrtsa = SQRT(sa)
                0795         ELSE
                0796          sa     = 0. _d 0
                0797          sqrtsa = 0. _d 0
                0798         ENDIF
                0799         p       = locPres(i,j)*SItodBar
1eadaea85b Jean*0800
                0801         rhoNum(i,j) = teos(01)
                0802      &   + ct*(teos(02) + ct*(teos(03) + teos(04)*ct))
                0803      &   + sa*(teos(05) + ct*(teos(06) + teos(07)*ct)
                0804      &   + sqrtsa*(teos(08) + ct*(teos(09)
ac1c5b24e7 Mart*0805      &                      + ct*(teos(10) + teos(11)*ct))))
1eadaea85b Jean*0806      &   + p*(teos(12) + ct*(teos(13) + teos(14)*ct)
                0807      &                 + sa*(teos(15) + teos(16)*ct)
ac1c5b24e7 Mart*0808      &   + p*(teos(17) + ct*(teos(18) + teos(19)*ct) + teos(20)*sa))
1eadaea85b Jean*0809
                0810         den = teos(21)
                0811      &   + ct*(teos(22) + ct*(teos(23) + ct*(teos(24) + teos(25)*ct)))
                0812      &   + sa*(teos(26) + ct*(teos(27) + ct*(teos(28)
                0813      &   + ct*(teos(29) + teos(30)*ct)))
                0814      &   + teos(36)*sa
                0815      &   + sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
                0816      &            + ct*(teos(34) + teos(35)*ct)))))
                0817      &   + p*(teos(37) + ct*(teos(38) + ct*(teos(39) + teos(40)*ct))
                0818      &   + sa*(teos(41) + teos(42)*ct)
                0819      &   + p*(teos(43) + ct*(teos(44) + teos(45)*ct + teos(46)*sa)
ac1c5b24e7 Mart*0820      &   + p*(teos(47) + teos(48)*ct)))
1eadaea85b Jean*0821
ac1c5b24e7 Mart*0822         rhoDen(i,j) = 1.0/(epsln+den)
                0823
                0824        ENDDO
                0825       ENDDO
                0826
                0827       RETURN
                0828       END
                0829
                0830 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0831 CBOP
22cb42e69a Jean*0832 C     !ROUTINE: FIND_RHO_SCALAR
                0833 C     !INTERFACE:
                0834       SUBROUTINE FIND_RHO_SCALAR(
                0835      I     tLoc, sLoc, pLoc,
                0836      O     rhoLoc,
                0837      I     myThid )
                0838
                0839 C     !DESCRIPTION: \bv
                0840 C     *==========================================================*
                0841 C     | o SUBROUTINE FIND_RHO_SCALAR
                0842 C     |   Calculates rho(S,T,p)
533b1afce7 Jean*0843 C     | Note: for seawater EOS, p is pressure (Pa) relative to
                0844 C     |   reference surface pressure "surf_pRef"
22cb42e69a Jean*0845 C     *==========================================================*
                0846 C     \ev
                0847
                0848 C     !USES:
                0849       IMPLICIT NONE
                0850 C     == Global variables ==
                0851 #include "SIZE.h"
                0852 #include "EEPARAMS.h"
                0853 #include "PARAMS.h"
                0854 #include "EOS.h"
                0855
                0856 C     !INPUT/OUTPUT PARAMETERS:
                0857 C     == Routine arguments ==
                0858       _RL tLoc, sLoc, pLoc
                0859       _RL rhoLoc
                0860       INTEGER myThid
                0861
                0862 C     !LOCAL VARIABLES:
                0863 C     == Local variables ==
                0864       _RL t1, t2, t3, t4, s1, s3o2, p1, p2, sp5, p1t1
533b1afce7 Jean*0865       _RL dp0, ct, sa, sqrtsa, p
22cb42e69a Jean*0866       _RL rfresh, rsalt, rhoP0
                0867       _RL bMfresh, bMsalt, bMpres, BulkMod
                0868       _RL rhoNum, rhoDen, den, epsln
4ee1a938b6 Jean*0869       _RL oneMinusKap, facPres, theta_v
22cb42e69a Jean*0870       PARAMETER ( epsln = 0. _d 0 )
                0871       CHARACTER*(MAX_LEN_MBUF) msgBuf
e34b4d724c Mart*0872 #ifdef USE_FACTORIZED_EOS
                0873       _RL eosF1, eosF2, eosF3, eosF4, eosF5, eosF6
                0874       _RL eosS1, eosS2, eosS3, eosS4, eosS5, eosS6
                0875       _RL eosS7, eosS8, eosS9
                0876       _RL eosP1, eosP2, eosP3, eosP4, eosP5, eosP6, eosP7
                0877       _RL eosP8, eosP9, eosP10, eosP11, eosP12, eosP13, eosP14
                0878 #endif /* USE_FACTORIZED_EOS */
22cb42e69a Jean*0879 CEOP
e34b4d724c Mart*0880 #ifdef USE_FACTORIZED_EOS
1eadaea85b Jean*0881       eosF1 = eosJMDCFw(1)
e34b4d724c Mart*0882       eosF2 = eosJMDCFw(2)
                0883       eosF3 = eosJMDCFw(3)
                0884       eosF4 = eosJMDCFw(4)
                0885       eosF5 = eosJMDCFw(5)
                0886       eosF6 = eosJMDCFw(6)
                0887       eosS1 = eosJMDCSw(1)
                0888       eosS2 = eosJMDCSw(2)
                0889       eosS3 = eosJMDCSw(3)
                0890       eosS4 = eosJMDCSw(4)
                0891       eosS5 = eosJMDCSw(5)
                0892       eosS6 = eosJMDCSw(6)
                0893       eosS7 = eosJMDCSw(7)
                0894       eosS8 = eosJMDCSw(8)
                0895       eosS9 = eosJMDCSw(9)
                0896
                0897       eosP1  =eosJMDCKP(1)
                0898       eosP2  =eosJMDCKP(2)
                0899       eosP3  =eosJMDCKP(3)
                0900       eosP4  =eosJMDCKP(4)
                0901       eosP5  =eosJMDCKP(5)
                0902       eosP6  =eosJMDCKP(6)
                0903       eosP7  =eosJMDCKP(7)
                0904       eosP8  =eosJMDCKP(8)
                0905       eosP9  =eosJMDCKP(9)
                0906       eosP10 =eosJMDCKP(10)
                0907       eosP11 =eosJMDCKP(11)
                0908       eosP12 =eosJMDCKP(12)
                0909       eosP13 =eosJMDCKP(13)
                0910       eosP14 =eosJMDCKP(14)
                0911 #endif /* USE_FACTORIZED_EOS */
22cb42e69a Jean*0912
                0913       rhoLoc  = 0. _d 0
                0914       rhoP0   = 0. _d 0
                0915       bulkMod = 0. _d 0
                0916       rfresh  = 0. _d 0
                0917       rsalt   = 0. _d 0
                0918       bMfresh = 0. _d 0
                0919       bMsalt  = 0. _d 0
                0920       bMpres  = 0. _d 0
                0921       rhoNum  = 0. _d 0
                0922       rhoDen  = 0. _d 0
                0923       den     = 0. _d 0
533b1afce7 Jean*0924       dp0     = surf_pRef - eosRefP0
22cb42e69a Jean*0925
                0926       t1 = tLoc
                0927       t2 = t1*t1
                0928       t3 = t2*t1
                0929       t4 = t3*t1
                0930
                0931       s1  = sLoc
                0932       IF ( s1 .LT. 0. _d 0 ) THEN
                0933 C     issue a warning
                0934          WRITE(msgBuf,'(A,E13.5)')
                0935      &        ' FIND_RHO_SCALAR:   WARNING, salinity = ', s1
                0936          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0937      &                       SQUEEZE_RIGHT , myThid )
                0938          s1 = 0. _d 0
                0939       ENDIF
                0940
                0941       IF (equationOfState.EQ.'LINEAR') THEN
                0942
                0943          rholoc = rhoNil*(
                0944      &                      sBeta *(sLoc-sRef(1))
                0945      &                     -tAlpha*(tLoc-tRef(1))
                0946      &                   ) + rhoNil
                0947 c        rhoLoc = 0. _d  0
                0948
                0949       ELSEIF (equationOfState.EQ.'POLY3') THEN
                0950
                0951 C     this is not correct, there is a field eosSig0 which should be use here
                0952 C     but I DO not intent to include the reference level in this routine
                0953          WRITE(msgBuf,'(A)')
                0954      &        ' FIND_RHO_SCALAR: for POLY3, the density is not'
                0955          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0956      &                       SQUEEZE_RIGHT , myThid )
                0957          WRITE(msgBuf,'(A)')
                0958      &         '                 computed correctly in this routine'
                0959          CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
                0960      &                       SQUEEZE_RIGHT , myThid )
                0961          rhoLoc = 0. _d 0
                0962
                0963       ELSEIF ( equationOfState(1:5).EQ.'JMD95'
                0964      &      .OR. equationOfState.EQ.'UNESCO' ) THEN
                0965 C     nonlinear equation of state in pressure coordinates
                0966
                0967          s3o2 = s1*SQRT(s1)
                0968
533b1afce7 Jean*0969          p1 = ( pLoc + dp0 )*SItoBar
22cb42e69a Jean*0970          p2 = p1*p1
                0971
                0972 C     density of freshwater at the surface
e34b4d724c Mart*0973 #ifdef USE_FACTORIZED_EOS
1eadaea85b Jean*0974          rfresh =
e34b4d724c Mart*0975      &              eosF1 +
                0976      &        t1* ( eosF2 +
                0977      &        t1* ( eosF3 +
                0978      &        t1* ( eosF4 +
                0979      &        t1* ( eosF5 +
                0980      &        t1*   eosF6  ))))
                0981 #else
22cb42e69a Jean*0982          rfresh =
                0983      &          eosJMDCFw(1)
                0984      &        + eosJMDCFw(2)*t1
                0985      &        + eosJMDCFw(3)*t2
                0986      &        + eosJMDCFw(4)*t3
                0987      &        + eosJMDCFw(5)*t4
                0988      &        + eosJMDCFw(6)*t4*t1
e34b4d724c Mart*0989 #endif /* USE_FACTORIZED_EOS */
22cb42e69a Jean*0990 C     density of sea water at the surface
e34b4d724c Mart*0991 #ifdef USE_FACTORIZED_EOS
                0992          rsalt =   s1 * ( eosS1 +
                0993      &             t1 * ( eosS2 +
                0994      &             t1 * ( eosS3 +
                0995      &             t1 * ( eosS4 +
                0996      &             t1 *   eosS5  ))))
                0997      &         + s3o2 * ( eosS6 +
1eadaea85b Jean*0998      &             t1 * ( eosS7 +
e34b4d724c Mart*0999      &             t1 *   eosS8  ))
                1000      &         + s1*s1*   eosS9
                1001 #else
22cb42e69a Jean*1002          rsalt =
                1003      &        s1*(
                1004      &             eosJMDCSw(1)
                1005      &           + eosJMDCSw(2)*t1
                1006      &           + eosJMDCSw(3)*t2
                1007      &           + eosJMDCSw(4)*t3
                1008      &           + eosJMDCSw(5)*t4
                1009      &           )
                1010      &        + s3o2*(
                1011      &             eosJMDCSw(6)
                1012      &           + eosJMDCSw(7)*t1
                1013      &           + eosJMDCSw(8)*t2
                1014      &           )
                1015      &           + eosJMDCSw(9)*s1*s1
e34b4d724c Mart*1016 #endif /* USE_FACTORIZED_EOS */
22cb42e69a Jean*1017
                1018          rhoP0 = rfresh + rsalt
                1019
                1020 C     secant bulk modulus of fresh water at the surface
e34b4d724c Mart*1021 #ifdef USE_FACTORIZED_EOS
                1022          bMfresh =   eosF1 +
                1023      &        t1 * ( eosF2 +
                1024      &        t1 * ( eosF3 +
                1025      &        t1 * ( eosF4 +
                1026      &        t1 *   eosF5  )))
1eadaea85b Jean*1027 #else
22cb42e69a Jean*1028          bMfresh =
                1029      &             eosJMDCKFw(1)
                1030      &           + eosJMDCKFw(2)*t1
                1031      &           + eosJMDCKFw(3)*t2
                1032      &           + eosJMDCKFw(4)*t3
                1033      &           + eosJMDCKFw(5)*t4
e34b4d724c Mart*1034 #endif /* USE_FACTORIZED_EOS */
                1035 C     secant bulk modulus of sea water at the surface
                1036 #ifdef USE_FACTORIZED_EOS
                1037          bMsalt = s1 * ( eosS1 +
                1038      &            t1 * ( eosS2 +
                1039      &            t1 * ( eosS3 +
                1040      &            t1 *   eosS4  )))
                1041      &        + s3o2 * ( eosS5 +
                1042      &            t1 * ( eosS6 +
                1043      &            t1 *   eosS7  ))
                1044 C     secant bulk modulus of sea water at pressure p
                1045          bMpres =  p1 * ( eosP1 +
                1046      &             t1 * ( eosP2 +
                1047      &             t1 * ( eosP3 +
                1048      &             t1 *   eosP4  )))
                1049      &          + p1*s1*( eosP5 +
                1050      &             t1 * ( eosP6 +
                1051      &             t1 *   eosP7  ))
                1052      &          + p1*s3o2*eosP8
                1053      &          + p2  * ( eosP9 +
                1054      &             t1 * ( eosP10 +
                1055      &             t1 *   eosP11  ))
                1056      &         + p2*s1* ( eosP12 +
                1057      &             t1 * ( eosP13 +
                1058      &             t1 *   eosP14 ))
                1059 #else
22cb42e69a Jean*1060 C     secant bulk modulus of sea water at the surface
                1061          bMsalt =
                1062      &        s1*( eosJMDCKSw(1)
                1063      &           + eosJMDCKSw(2)*t1
                1064      &           + eosJMDCKSw(3)*t2
                1065      &           + eosJMDCKSw(4)*t3
                1066      &           )
                1067      &    + s3o2*( eosJMDCKSw(5)
                1068      &           + eosJMDCKSw(6)*t1
                1069      &           + eosJMDCKSw(7)*t2
                1070      &           )
                1071 C     secant bulk modulus of sea water at pressure p
                1072          bMpres =
                1073      &        p1*( eosJMDCKP(1)
                1074      &           + eosJMDCKP(2)*t1
                1075      &           + eosJMDCKP(3)*t2
                1076      &           + eosJMDCKP(4)*t3
                1077      &           )
                1078      &   + p1*s1*( eosJMDCKP(5)
                1079      &           + eosJMDCKP(6)*t1
                1080      &           + eosJMDCKP(7)*t2
                1081      &           )
                1082      &      + p1*s3o2*eosJMDCKP(8)
                1083      &      + p2*( eosJMDCKP(9)
                1084      &           + eosJMDCKP(10)*t1
                1085      &           + eosJMDCKP(11)*t2
                1086      &           )
                1087      &    + p2*s1*( eosJMDCKP(12)
                1088      &           + eosJMDCKP(13)*t1
                1089      &           + eosJMDCKP(14)*t2
                1090      &           )
e34b4d724c Mart*1091 #endif /* USE_FACTORIZED_EOS */
22cb42e69a Jean*1092
                1093          bulkMod = bMfresh + bMsalt + bMpres
                1094
                1095 C     density of sea water at pressure p
                1096          rhoLoc = rhoP0/(1. _d 0 - p1/bulkMod)
                1097
                1098       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
                1099
e34b4d724c Mart*1100        sp5 = SQRT(s1)
                1101
533b1afce7 Jean*1102        p1 = ( pLoc + dp0 )*SItodBar
e34b4d724c Mart*1103        p1t1 = p1*t1
                1104
                1105        rhoNum = eosMDJWFnum(0)
                1106      &      + t1*(eosMDJWFnum(1)
                1107      &      +     t1*(eosMDJWFnum(2) + eosMDJWFnum(3)*t1) )
                1108      &      + s1*(eosMDJWFnum(4)
                1109      &      +     eosMDJWFnum(5)*t1  + eosMDJWFnum(6)*s1)
                1110      &      + p1*(eosMDJWFnum(7) + eosMDJWFnum(8)*t2
                1111      &      +     eosMDJWFnum(9)*s1
                1112      &      +     p1*(eosMDJWFnum(10) + eosMDJWFnum(11)*t2) )
                1113
                1114        den = eosMDJWFden(0)
                1115      &      + t1*(eosMDJWFden(1)
                1116      &      +     t1*(eosMDJWFden(2)
                1117      &      +         t1*(eosMDJWFden(3) + t1*eosMDJWFden(4) ) ) )
                1118      &      + s1*(eosMDJWFden(5)
                1119      &      +     t1*(eosMDJWFden(6)
                1120      &      +         eosMDJWFden(7)*t2)
                1121      &      +     sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2) )
                1122      &      + p1*(eosMDJWFden(10)
                1123      &      +     p1t1*(eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1) )
                1124
                1125        rhoDen = 1.0/(epsln+den)
1eadaea85b Jean*1126
e34b4d724c Mart*1127        rhoLoc = rhoNum*rhoDen
22cb42e69a Jean*1128
ac1c5b24e7 Mart*1129       ELSEIF( equationOfState .EQ. 'TEOS10' ) THEN
                1130
                1131        ct      = tLoc
                1132        sa      = sLoc
                1133        IF ( sa .GT. 0. _d 0 ) THEN
                1134         sqrtsa = SQRT(sa)
                1135        ELSE
                1136         sa     = 0. _d 0
                1137         sqrtsa = 0. _d 0
                1138        ENDIF
533b1afce7 Jean*1139        p = ( pLoc + dp0 )*SItodBar
1eadaea85b Jean*1140
                1141        rhoNum = teos(01)
                1142      &   + ct*(teos(02) + ct*(teos(03) + teos(04)*ct))
                1143      &   + sa*(teos(05) + ct*(teos(06) + teos(07)*ct)
                1144      &   + sqrtsa*(teos(08) + ct*(teos(09)
ac1c5b24e7 Mart*1145      &            + ct*(teos(10) + teos(11)*ct))))
1eadaea85b Jean*1146      &   + p*(teos(12) + ct*(teos(13) + teos(14)*ct)
                1147      &                      + sa*(teos(15) + teos(16)*ct)
ac1c5b24e7 Mart*1148      &   + p*(teos(17) + ct*(teos(18) + teos(19)*ct) + teos(20)*sa))
1eadaea85b Jean*1149
                1150        den = teos(21)
                1151      &   + ct*(teos(22) + ct*(teos(23) + ct*(teos(24) + teos(25)*ct)))
                1152      &   + sa*(teos(26) + ct*(teos(27) + ct*(teos(28)
                1153      &   + ct*(teos(29) + teos(30)*ct)))
                1154      &   + teos(36)*sa
                1155      %   + sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
                1156      &            + ct*(teos(34) + teos(35)*ct)))))
                1157      %   + p*(teos(37) + ct*(teos(38) + ct*(teos(39) + teos(40)*ct))
                1158      %   + sa*(teos(41) + teos(42)*ct)
                1159      %   + p*(teos(43) + ct*(teos(44) + teos(45)*ct + teos(46)*sa)
ac1c5b24e7 Mart*1160      %   + p*(teos(47) + teos(48)*ct)))
1eadaea85b Jean*1161
ac1c5b24e7 Mart*1162        rhoDen = 1.0/(epsln+den)
1eadaea85b Jean*1163
ac1c5b24e7 Mart*1164        rhoLoc = rhoNum*rhoDen
                1165
22cb42e69a Jean*1166       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
4ee1a938b6 Jean*1167
                1168          oneMinusKap = oneRL - atm_kappa
                1169 C        rho = P/(R*T_v) = Po/(R*theta_v)*(P/Po)^(1-kappa)
                1170          facPres = ( pLoc/atm_Po )**oneMinusKap
                1171          theta_v = tLoc*( sLoc*atm_Rq + oneRL )
                1172          rhoLoc = atm_Po*facPres /( atm_Rd*theta_v )
                1173
22cb42e69a Jean*1174       ELSE
                1175        WRITE(msgBuf,'(3A)')
                1176      &        ' FIND_RHO_SCALAR : equationOfState = "',
                1177      &        equationOfState,'"'
                1178        CALL PRINT_ERROR( msgBuf, myThid )
                1179        STOP 'ABNORMAL END: S/R FIND_RHO_SCALAR'
                1180       ENDIF
                1181
                1182       RETURN
                1183       END
                1184
                1185 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
59e59b79d8 Mart*1186 CBOP
                1187 C     !ROUTINE: LOOK_FOR_NEG_SALINITY
                1188 C     !INTERFACE:
6ac14281aa Jean*1189       SUBROUTINE LOOK_FOR_NEG_SALINITY(
de7f345328 Jean*1190      I                iMin, iMax, jMin, jMax,
                1191      U                sFld,
                1192      I                k, bi, bj, myThid )
59e59b79d8 Mart*1193
                1194 C     !DESCRIPTION: \bv
                1195 C     *==========================================================*
                1196 C     | o SUBROUTINE LOOK_FOR_NEG_SALINITY
                1197 C     |   looks for and fixes negative salinity values
6ac14281aa Jean*1198 C     |   this is necessary IF the equation of state uses
59e59b79d8 Mart*1199 C     |   the square root of salinity
                1200 C     *==========================================================*
de7f345328 Jean*1201 C     | k - is the Salt level
59e59b79d8 Mart*1202 C     *==========================================================*
                1203 C     \ev
                1204
                1205 C     !USES:
6ac14281aa Jean*1206       IMPLICIT NONE
59e59b79d8 Mart*1207 C     == Global variables ==
                1208 #include "SIZE.h"
                1209 #include "EEPARAMS.h"
                1210 #include "PARAMS.h"
                1211
                1212 C     !INPUT/OUTPUT PARAMETERS:
                1213 C     == Routine arguments ==
de7f345328 Jean*1214 C     k    :: Level of Salt slice
22cb42e69a Jean*1215       INTEGER iMin,iMax,jMin,jMax
de7f345328 Jean*1216       _RL     sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
22cb42e69a Jean*1217       INTEGER k, bi, bj
6ac14281aa Jean*1218       INTEGER myThid
59e59b79d8 Mart*1219
                1220 C     !LOCAL VARIABLES:
                1221 C     == Local variables ==
6ac14281aa Jean*1222       INTEGER i,j, localWarning
22cb42e69a Jean*1223       CHARACTER*(MAX_LEN_MBUF) msgBuf
59e59b79d8 Mart*1224 CEOP
                1225
                1226       localWarning = 0
6ac14281aa Jean*1227       DO j=jMin,jMax
                1228        DO i=iMin,iMax
59e59b79d8 Mart*1229 C     abbreviations
de7f345328 Jean*1230         IF ( sFld(i,j) .LT. 0. _d 0 ) THEN
59e59b79d8 Mart*1231          localWarning = localWarning + 1
de7f345328 Jean*1232          sFld(i,j) = 0. _d 0
6ac14281aa Jean*1233         ENDIF
                1234        ENDDO
                1235       ENDDO
59e59b79d8 Mart*1236 C     issue a warning
6ac14281aa Jean*1237       IF ( localWarning .GT. 0 ) THEN
22cb42e69a Jean*1238         WRITE(msgBuf,'(2A,I5,A,2I4)') 'S/R LOOK_FOR_NEG_SALINITY:',
                1239      &      ' from level k =', k, ' ; bi,bj =', bi, bj
                1240         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                1241      &                      SQUEEZE_RIGHT , myThid )
                1242         WRITE(msgBuf,'(2A,I6,A)') 'S/R LOOK_FOR_NEG_SALINITY:',
                1243      &      ' reset to zero', localWarning, ' negative salinity.'
                1244         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                1245      &                      SQUEEZE_RIGHT , myThid )
6ac14281aa Jean*1246       ENDIF
59e59b79d8 Mart*1247
de7f345328 Jean*1248       RETURN
6ac14281aa Jean*1249       END