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