model/src/find_alpha.F

1eadaea85b Jean*0001 #include "PACKAGES_CONFIG.h"
fb3dc7d949 Alis*0002 #include "CPP_OPTIONS.h"
                0003 #define USE_FACTORIZED_POLY
                0004
1eadaea85b Jean*0005 C--  File find_alpha.F: Calculate expansion Coeff
                0006 C--   Contents
                0007 C--   o FIND_ALPHA
                0008 C--   o FIND_BETA
                0009
                0010 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
9366854e02 Chri*0011 CBOP
                0012 C     !ROUTINE: FIND_ALPHA
                0013 C     !INTERFACE:
                0014       SUBROUTINE FIND_ALPHA (
de7f345328 Jean*0015      I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
                0016      O     alphaLoc,
18a5ed860a Mart*0017      I     myThid )
fb3dc7d949 Alis*0018
9366854e02 Chri*0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
de7f345328 Jean*0021 C     | o SUBROUTINE FIND_ALPHA
                0022 C     |   Calculates [drho(S,T,z) / dT] of a horizontal slice
9366854e02 Chri*0023 C     *==========================================================*
1eadaea85b Jean*0024 C     | k        - is the Theta/Salt level
                0025 C     | kRef     - determines pressure reference level
                0026 C     |            (not used in 'LINEAR' mode)
de7f345328 Jean*0027 C     | alphaLoc - drho / dT (kg/m^3/C)
9366854e02 Chri*0028 C     *==========================================================*
                0029 C     \ev
                0030
                0031 C     !USES:
                0032       IMPLICIT NONE
6ac14281aa Jean*0033 C     === Global variables ===
fb3dc7d949 Alis*0034 #include "SIZE.h"
                0035 #include "DYNVARS.h"
                0036 #include "EEPARAMS.h"
                0037 #include "PARAMS.h"
a37a13034c Mart*0038 #include "EOS.h"
                0039 #include "GRID.h"
fb3dc7d949 Alis*0040
9366854e02 Chri*0041 C     !INPUT/OUTPUT PARAMETERS:
6ac14281aa Jean*0042 C     == Routine arguments ==
18a5ed860a Mart*0043 C     k      :: Level of Theta/Salt slice
                0044 C     kRef   :: Pressure reference level
                0045 c     myThid :: thread number for this instance of the routine
                0046       INTEGER myThid
6ac14281aa Jean*0047       INTEGER bi,bj,iMin,iMax,jMin,jMax
                0048       INTEGER k
                0049       INTEGER kRef
1eadaea85b Jean*0050       _RL alphaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
fb3dc7d949 Alis*0051
9366854e02 Chri*0052 C     !LOCAL VARIABLES:
de7f345328 Jean*0053 C     == Local variables ==
6ac14281aa Jean*0054       INTEGER i,j
533b1afce7 Jean*0055       _RL refTemp, refSalt, tP, sP, dp0
31a5949e20 Mart*0056       _RL t1, t2, t3, s1, s3o2, p1, p2, sp5, p1t1
ac1c5b24e7 Mart*0057       _RL ct, sa, sqrtsa, p
a37a13034c Mart*0058       _RL drhoP0dtheta, drhoP0dthetaFresh, drhoP0dthetaSalt
                0059       _RL dKdtheta, dKdthetaFresh, dKdthetaSalt, dKdthetaPres
1eadaea85b Jean*0060       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL rhoP0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0062       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31a5949e20 Mart*0063       _RL dnum_dtheta, dden_dtheta
1eadaea85b Jean*0064       _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0065       _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
9366854e02 Chri*0066 CEOP
fb3dc7d949 Alis*0067
59e59b79d8 Mart*0068 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
de7f345328 Jean*0069 c     CALL LOOK_FOR_NEG_SALINITY(
                0070 c    I          iMin, iMax, jMin, jMax,
                0071 c    U          sFld,
                0072 c    I          k, bi, bj, myThid )
59e59b79d8 Mart*0073 #endif
                0074
6ac14281aa Jean*0075       IF (equationOfState.EQ.'LINEAR') THEN
fb3dc7d949 Alis*0076
6ac14281aa Jean*0077          DO j=jMin,jMax
                0078             DO i=iMin,iMax
de7f345328 Jean*0079                alphaLoc(i,j) = -rhonil * tAlpha
6ac14281aa Jean*0080             ENDDO
                0081          ENDDO
de7f345328 Jean*0082
6ac14281aa Jean*0083       ELSEIF (equationOfState.EQ.'POLY3') THEN
fb3dc7d949 Alis*0084
                0085          refTemp=eosRefT(kRef)
                0086          refSalt=eosRefS(kRef)
                0087
6ac14281aa Jean*0088          DO j=jMin,jMax
                0089             DO i=iMin,iMax
fb3dc7d949 Alis*0090                tP=theta(i,j,k,bi,bj)-refTemp
                0091                sP=salt(i,j,k,bi,bj)-refSalt
                0092 #ifdef USE_FACTORIZED_POLY
de7f345328 Jean*0093                alphaLoc(i,j) =
fb3dc7d949 Alis*0094      &     ( eosC(6,kRef)
                0095      &         *tP*3.
                0096      &        +(eosC(7,kRef)*sP + eosC(3,kRef))*2.
                0097      &       )*tP
                0098      &      +(eosC(8,kRef)*sP + eosC(4,kRef) )*sP + eosC(1,kRef)
de7f345328 Jean*0099      &
fb3dc7d949 Alis*0100 #else
de7f345328 Jean*0101                alphaLoc(i,j) =
fb3dc7d949 Alis*0102      &              eosC(1,kRef)                +
                0103      &              eosC(3,kRef)*tP*2.          +
                0104      &              eosC(4,kRef)         *sP    +
                0105      &              eosC(6,kRef)*tP*tP*3.       +
                0106      &              eosC(7,kRef)*tP*2.   *sP    +
                0107      &              eosC(8,kRef)         *sP*sP
                0108 #endif
6ac14281aa Jean*0109             ENDDO
                0110          ENDDO
de7f345328 Jean*0111
6ac14281aa Jean*0112       ELSEIF ( equationOfState(1:5).EQ.'JMD95'
                0113      &        .OR. equationOfState.EQ.'UNESCO' ) THEN
a37a13034c Mart*0114 C     nonlinear equation of state in pressure coordinates
533b1afce7 Jean*0115          dp0 = surf_pRef - eosRefP0
a37a13034c Mart*0116
6ac14281aa Jean*0117          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0118      I        bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
6ac14281aa Jean*0119      O        locPres,
de7f345328 Jean*0120      I        myThid )
6ac14281aa Jean*0121
a37a13034c Mart*0122          CALL FIND_RHOP0(
de7f345328 Jean*0123      I        iMin, iMax, jMin, jMax,
                0124      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
a37a13034c Mart*0125      O        rhoP0,
                0126      I        myThid )
de7f345328 Jean*0127
a37a13034c Mart*0128          CALL FIND_BULKMOD(
de7f345328 Jean*0129      I        iMin, iMax, jMin, jMax, locPres,
                0130      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
a37a13034c Mart*0131      O        bulkMod,
                0132      I        myThid )
                0133
6ac14281aa Jean*0134          DO j=jMin,jMax
                0135             DO i=iMin,iMax
a37a13034c Mart*0136
                0137 C     abbreviations
de7f345328 Jean*0138                t1 = theta(i,j,k,bi,bj)
31a5949e20 Mart*0139                t2 = t1*t1
                0140                t3 = t2*t1
de7f345328 Jean*0141
02d90fb24c Jean*0142 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0143 C     this line makes TAF create vectorizable code
                0144                s3o2 = 0. _d 0
                0145 #endif
31a5949e20 Mart*0146                s1  = salt(i,j,k,bi,bj)
59e581f4e0 Jean*0147               IF ( s1 .GT. 0. _d 0 ) THEN
6ac14281aa Jean*0148                s3o2 = SQRT(s1*s1*s1)
59e581f4e0 Jean*0149               ELSE
                0150                s1   = 0. _d 0
                0151                s3o2 = 0. _d 0
                0152               ENDIF
de7f345328 Jean*0153
6ac14281aa Jean*0154                p1  = locPres(i,j)*SItoBar
31a5949e20 Mart*0155                p2 = p1*p1
a37a13034c Mart*0156
                0157 C     d(rho)/d(theta)
                0158 C     of fresh water at p = 0
de7f345328 Jean*0159                drhoP0dthetaFresh =
a37a13034c Mart*0160      &                eosJMDCFw(2)
31a5949e20 Mart*0161      &           + 2.*eosJMDCFw(3)*t1
a37a13034c Mart*0162      &           + 3.*eosJMDCFw(4)*t2
                0163      &           + 4.*eosJMDCFw(5)*t3
31a5949e20 Mart*0164      &           + 5.*eosJMDCFw(6)*t3*t1
a37a13034c Mart*0165 C     of salt water at p = 0
de7f345328 Jean*0166                drhoP0dthetaSalt =
31a5949e20 Mart*0167      &        s1*(
a37a13034c Mart*0168      &                eosJMDCSw(2)
31a5949e20 Mart*0169      &           + 2.*eosJMDCSw(3)*t1
a37a13034c Mart*0170      &           + 3.*eosJMDCSw(4)*t2
                0171      &           + 4.*eosJMDCSw(5)*t3
                0172      &           )
                0173      &       + s3o2*(
                0174      &           +    eosJMDCSw(7)
31a5949e20 Mart*0175      &           + 2.*eosJMDCSw(8)*t1
a37a13034c Mart*0176      &           )
                0177 C     d(bulk modulus)/d(theta)
                0178 C     of fresh water at p = 0
de7f345328 Jean*0179                dKdthetaFresh =
a37a13034c Mart*0180      &                eosJMDCKFw(2)
31a5949e20 Mart*0181      &           + 2.*eosJMDCKFw(3)*t1
a37a13034c Mart*0182      &           + 3.*eosJMDCKFw(4)*t2
                0183      &           + 4.*eosJMDCKFw(5)*t3
                0184 C     of sea water at p = 0
                0185                dKdthetaSalt =
31a5949e20 Mart*0186      &        s1*(    eosJMDCKSw(2)
                0187      &           + 2.*eosJMDCKSw(3)*t1
a37a13034c Mart*0188      &           + 3.*eosJMDCKSw(4)*t2
                0189      &           )
                0190      &    + s3o2*(    eosJMDCKSw(6)
31a5949e20 Mart*0191      &           + 2.*eosJMDCKSw(7)*t1
a37a13034c Mart*0192      &           )
                0193 C     of sea water at p
                0194                dKdthetaPres =
31a5949e20 Mart*0195      &        p1*(    eosJMDCKP(2)
                0196      &           + 2.*eosJMDCKP(3)*t1
a37a13034c Mart*0197      &           + 3.*eosJMDCKP(4)*t2
                0198      &           )
31a5949e20 Mart*0199      &   + p1*s1*(    eosJMDCKP(6)
                0200      &           + 2.*eosJMDCKP(7)*t1
a37a13034c Mart*0201      &           )
                0202      &      + p2*(    eosJMDCKP(10)
31a5949e20 Mart*0203      &           + 2.*eosJMDCKP(11)*t1
a37a13034c Mart*0204      &           )
31a5949e20 Mart*0205      &   + p2*s1*(    eosJMDCKP(13)
                0206      &           + 2.*eosJMDCKP(14)*t1
a37a13034c Mart*0207      &           )
                0208
de7f345328 Jean*0209                drhoP0dtheta  = drhoP0dthetaFresh
a37a13034c Mart*0210      &                       + drhoP0dthetaSalt
de7f345328 Jean*0211                dKdtheta      = dKdthetaFresh
                0212      &                       + dKdthetaSalt
a37a13034c Mart*0213      &                       + dKdthetaPres
de7f345328 Jean*0214                alphaLoc(i,j) =
a37a13034c Mart*0215      &              ( bulkmod(i,j)**2*drhoP0dtheta
31a5949e20 Mart*0216      &              - bulkmod(i,j)*p1*drhoP0dtheta
                0217      &              - rhoP0(i,j)*p1*dKdtheta )
                0218      &              /( bulkmod(i,j) - p1 )**2
de7f345328 Jean*0219
6ac14281aa Jean*0220             ENDDO
                0221          ENDDO
                0222       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
533b1afce7 Jean*0223          dp0 = surf_pRef - eosRefP0
6ac14281aa Jean*0224
                0225          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0226      I        bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
6ac14281aa Jean*0227      O        locPres,
de7f345328 Jean*0228      I        myThid )
edd57506ae Patr*0229
de7f345328 Jean*0230          CALL FIND_RHONUM(
                0231      I        iMin, iMax, jMin, jMax, locPres,
                0232      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0233      O        rhoLoc,
                0234      I        myThid )
24de0d2c47 Patr*0235
de7f345328 Jean*0236          CALL FIND_RHODEN(
                0237      I        iMin, iMax, jMin, jMax, locPres,
                0238      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0239      O        rhoDen,
                0240      I        myThid )
6ac14281aa Jean*0241
                0242          DO j=jMin,jMax
                0243             DO i=iMin,iMax
31a5949e20 Mart*0244                t1  = theta(i,j,k,bi,bj)
                0245                t2  = t1*t1
                0246                s1  = salt(i,j,k,bi,bj)
02d90fb24c Jean*0247 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0248 C     this line makes TAF create vectorizable code
                0249                sp5 = 0. _d 0
                0250 #endif
59e581f4e0 Jean*0251               IF ( s1 .GT. 0. _d 0 ) THEN
de7f345328 Jean*0252                sp5 = SQRT(s1)
59e581f4e0 Jean*0253               ELSE
                0254                s1  = 0. _d 0
                0255                sp5 = 0. _d 0
                0256               ENDIF
31a5949e20 Mart*0257
6ac14281aa Jean*0258                p1   = locPres(i,j)*SItodBar
31a5949e20 Mart*0259                p1t1 = p1*t1
de7f345328 Jean*0260
31a5949e20 Mart*0261                dnum_dtheta = eosMDJWFnum(1)
de7f345328 Jean*0262      &              + t1*(2.*eosMDJWFnum(2) + 3.*eosMDJWFnum(3)*t1)
                0263      &              + eosMDJWFnum(5)*s1
59e581f4e0 Jean*0264      &              + p1t1*(2.*eosMDJWFnum(8) + 2.*eosMDJWFnum(11)*p1)
de7f345328 Jean*0265
                0266                dden_dtheta = eosMDJWFden(1)
                0267      &              + t1*(2.*eosMDJWFden(2)
                0268      &              +     t1*(3.*eosMDJWFden(3)
31a5949e20 Mart*0269      &              +         4.*eosMDJWFden(4)*t1 ) )
de7f345328 Jean*0270      &              + s1*(eosMDJWFden(6)
                0271      &              +     t1*(3.*eosMDJWFden(7)*t1
31a5949e20 Mart*0272      &              +         2.*eosMDJWFden(9)*sp5 ) )
                0273      &              + p1*p1*(3.*eosMDJWFden(11)*t2 + eosMDJWFden(12)*p1)
de7f345328 Jean*0274
                0275                alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
edd57506ae Patr*0276      &              - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)
de7f345328 Jean*0277
6ac14281aa Jean*0278          ENDDO
                0279       ENDDO
31a5949e20 Mart*0280
ac1c5b24e7 Mart*0281       ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
533b1afce7 Jean*0282        dp0 = surf_pRef - eosRefP0
ac1c5b24e7 Mart*0283
                0284        CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0285      I      bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
ac1c5b24e7 Mart*0286      O      locPres,
                0287      I      myThid )
                0288
                0289        CALL FIND_RHOTEOS(
                0290      I      iMin, iMax, jMin, jMax, locPres,
                0291      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0292      O      rhoLoc, rhoDen,
                0293      I      myThid )
                0294
                0295        DO j=jMin,jMax
                0296         DO i=iMin,iMax
                0297          ct      = theta(i,j,k,bi,bj)
                0298          sa      = salt(i,j,k,bi,bj)
02d90fb24c Jean*0299 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0300 C     this line makes TAF create vectorizable code
                0301          sqrtsa  = 0. _d 0
                0302 #endif
ac1c5b24e7 Mart*0303          IF ( sa .GT. 0. _d 0 ) THEN
                0304           sqrtsa = SQRT(sa)
                0305          ELSE
                0306           sa     = 0. _d 0
                0307           sqrtsa = 0. _d 0
                0308          ENDIF
                0309          p       = locPres(i,j)*SItodBar
                0310
                0311          dnum_dtheta = teos(02)
1eadaea85b Jean*0312      &   + ct*(2.*teos(03) + 3.*teos(04)*ct)
ac1c5b24e7 Mart*0313      &   + sa*(teos(06) + 2.*teos(07)*ct
                0314      &   + sqrtsa*(teos(09) + ct*(2.*teos(10) + 3.*teos(11)*ct)))
                0315      &   + p*(     teos(13) + 2.*teos(14)*ct  + sa*2.*teos(16)
                0316      &        + p*(teos(18) + 2.*teos(19)*ct))
1eadaea85b Jean*0317
                0318          dden_dtheta = teos(22)
ac1c5b24e7 Mart*0319      &   + ct*(2.*teos(23) + ct*(3.*teos(24) + 4.*teos(25)*ct))
1eadaea85b Jean*0320      &   + sa*(teos(27)
ac1c5b24e7 Mart*0321      &   + ct*(2.*teos(28) + ct*(3.*teos(29) + 4.*teos(30)*ct))
1eadaea85b Jean*0322      &   + sqrtsa*(teos(32)
                0323      &   + ct*(2.*teos(33) + ct*(3.*teos(34) + 4.*teos(35)*ct))))
ac1c5b24e7 Mart*0324      &   + p*(teos(38) + ct*(2.*teos(39) + 3.*teos(40)*ct)
                0325      &   + teos(42)
                0326      &   + p*(teos(44) + 2.*teos(45)*ct + teos(46)*sa
                0327      &   + p*teos(48) ))
                0328
                0329          alphaLoc(i,j)    = rhoDen(i,j)*(dnum_dtheta
                0330      &        - (rhoLoc(i,j)*rhoDen(i,j))*dden_dtheta)
1eadaea85b Jean*0331
ac1c5b24e7 Mart*0332         ENDDO
                0333        ENDDO
                0334
6ac14281aa Jean*0335       ELSE
                0336          WRITE(*,*) 'FIND_ALPHA: equationOfState = ',equationOfState
                0337          STOP 'FIND_ALPHA: "equationOfState" has illegal value'
                0338       ENDIF
fb3dc7d949 Alis*0339
de7f345328 Jean*0340       RETURN
6ac14281aa Jean*0341       END
fb3dc7d949 Alis*0342
1eadaea85b Jean*0343 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0344 CBOP
                0345 C     !ROUTINE: FIND_BETA
                0346 C     !INTERFACE:
6ac14281aa Jean*0347       SUBROUTINE FIND_BETA (
9b24f7ff2d Mart*0348      I     bi, bj, iMin, iMax, jMin, jMax,  k, kRef,
de7f345328 Jean*0349      O     betaLoc,
18a5ed860a Mart*0350      I     myThid )
fb3dc7d949 Alis*0351
1eadaea85b Jean*0352 C     !DESCRIPTION: \bv
                0353 C     *==========================================================*
                0354 C     | o SUBROUTINE FIND_BETA
                0355 C     |   Calculates [drho(S,T,z) / dS] of a horizontal slice
                0356 C     |==========================================================*
                0357 C     | k       - is the Theta/Salt level
                0358 C     | kRef    - determines pressure reference level
                0359 C     |           (not used in 'LINEAR' mode)
ba0b047096 Mart*0360 C     | betaLoc - drho / dS (kg/m^3/(g/kg))
1eadaea85b Jean*0361 C     *==========================================================*
                0362 C     \ev
                0363
                0364 C     !USES:
                0365       IMPLICIT NONE
6ac14281aa Jean*0366 C     === Global variables ===
fb3dc7d949 Alis*0367 #include "SIZE.h"
                0368 #include "DYNVARS.h"
                0369 #include "EEPARAMS.h"
                0370 #include "PARAMS.h"
a37a13034c Mart*0371 #include "EOS.h"
                0372 #include "GRID.h"
fb3dc7d949 Alis*0373
1eadaea85b Jean*0374 C     !INPUT/OUTPUT PARAMETERS:
6ac14281aa Jean*0375 C     == Routine arguments ==
18a5ed860a Mart*0376 C     k      :: Level of Theta/Salt slice
                0377 C     kRef   :: Pressure reference level
                0378 c     myThid :: thread number for this instance of the routine
                0379       INTEGER myThid
6ac14281aa Jean*0380       INTEGER bi,bj,iMin,iMax,jMin,jMax
                0381       INTEGER k
                0382       INTEGER kRef
1eadaea85b Jean*0383       _RL betaLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
fb3dc7d949 Alis*0384
1eadaea85b Jean*0385 C     !LOCAL VARIABLES:
6ac14281aa Jean*0386 C     == Local variables ==
                0387       INTEGER i,j
533b1afce7 Jean*0388       _RL refTemp, refSalt, tP, sP, dp0
31a5949e20 Mart*0389       _RL t1, t2, t3, s1, s3o2, p1, sp5, p1t1
ac1c5b24e7 Mart*0390       _RL ct, sa, sqrtsa, p
a37a13034c Mart*0391       _RL drhoP0dS
                0392       _RL dKdS, dKdSSalt, dKdSPres
1eadaea85b Jean*0393       _RL locPres(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0394       _RL rhoP0  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0395       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
31a5949e20 Mart*0396       _RL dnum_dsalt, dden_dsalt
1eadaea85b Jean*0397       _RL rhoDen (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0398       _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a37a13034c Mart*0399 CEOP
                0400
59e59b79d8 Mart*0401 #ifdef CHECK_SALINITY_FOR_NEGATIVE_VALUES
de7f345328 Jean*0402 c     CALL LOOK_FOR_NEG_SALINITY(
                0403 c    I          iMin, iMax, jMin, jMax,
                0404 c    U          sFld,
                0405 c    I          k, bi, bj, myThid )
59e59b79d8 Mart*0406 #endif
                0407
6ac14281aa Jean*0408       IF (equationOfState.EQ.'LINEAR') THEN
fb3dc7d949 Alis*0409
6ac14281aa Jean*0410          DO j=jMin,jMax
                0411             DO i=iMin,iMax
de7f345328 Jean*0412                betaLoc(i,j) = rhonil * sBeta
6ac14281aa Jean*0413             ENDDO
                0414          ENDDO
de7f345328 Jean*0415
6ac14281aa Jean*0416       ELSEIF (equationOfState.EQ.'POLY3') THEN
fb3dc7d949 Alis*0417
                0418          refTemp=eosRefT(kRef)
                0419          refSalt=eosRefS(kRef)
                0420
6ac14281aa Jean*0421          DO j=jMin,jMax
                0422             DO i=iMin,iMax
fb3dc7d949 Alis*0423                tP=theta(i,j,k,bi,bj)-refTemp
                0424                sP=salt(i,j,k,bi,bj)-refSalt
                0425 #ifdef USE_FACTORIZED_POLY
de7f345328 Jean*0426                betaLoc(i,j) =
fb3dc7d949 Alis*0427      &    ( eosC(9,kRef)*sP*3. + eosC(5,kRef)*2. )*sP + eosC(2,kRef)
                0428      &   + ( eosC(7,kRef)*tP
                0429      &      +eosC(8,kRef)*sP*2. + eosC(4,kRef)
                0430      &     )*tP
                0431 #else
de7f345328 Jean*0432                betaLoc(i,j) =
fb3dc7d949 Alis*0433      &              eosC(2,kRef)                +
                0434      &              eosC(4,kRef)*tP             +
                0435      &              eosC(5,kRef)         *sP*2. +
                0436      &              eosC(7,kRef)*tP*tP          +
                0437      &              eosC(8,kRef)*tP      *sP*2. +
                0438      &              eosC(9,kRef)         *sP*sP*3.
                0439 #endif
6ac14281aa Jean*0440             ENDDO
                0441          ENDDO
fb3dc7d949 Alis*0442
6ac14281aa Jean*0443       ELSEIF ( equationOfState(1:5).EQ.'JMD95'
                0444      &        .OR. equationOfState.EQ.'UNESCO' ) THEN
a37a13034c Mart*0445 C     nonlinear equation of state in pressure coordinates
533b1afce7 Jean*0446          dp0 = surf_pRef - eosRefP0
a37a13034c Mart*0447
6ac14281aa Jean*0448          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0449      I        bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
6ac14281aa Jean*0450      O        locPres,
de7f345328 Jean*0451      I        myThid )
6ac14281aa Jean*0452
a37a13034c Mart*0453          CALL FIND_RHOP0(
de7f345328 Jean*0454      I        iMin, iMax, jMin, jMax,
                0455      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
a37a13034c Mart*0456      O        rhoP0,
                0457      I        myThid )
de7f345328 Jean*0458
a37a13034c Mart*0459          CALL FIND_BULKMOD(
de7f345328 Jean*0460      I        iMin, iMax, jMin, jMax, locPres,
                0461      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
a37a13034c Mart*0462      O        bulkMod,
                0463      I        myThid )
                0464
6ac14281aa Jean*0465          DO j=jMin,jMax
                0466             DO i=iMin,iMax
a37a13034c Mart*0467
                0468 C     abbreviations
de7f345328 Jean*0469                t1 = theta(i,j,k,bi,bj)
31a5949e20 Mart*0470                t2 = t1*t1
                0471                t3 = t2*t1
de7f345328 Jean*0472
31a5949e20 Mart*0473                s1  = salt(i,j,k,bi,bj)
02d90fb24c Jean*0474 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0475 C     this line makes TAF create vectorizable code
                0476                s3o2 = 0. _d 0
                0477 #endif
59e581f4e0 Jean*0478               IF ( s1 .GT. 0. _d 0 ) THEN
6ac14281aa Jean*0479                s3o2 = 1.5*SQRT(s1)
59e581f4e0 Jean*0480               ELSE
                0481                s1   = 0. _d 0
                0482                s3o2 = 0. _d 0
                0483               ENDIF
a37a13034c Mart*0484
6ac14281aa Jean*0485                p1  = locPres(i,j)*SItoBar
a37a13034c Mart*0486
                0487 C     d(rho)/d(S)
                0488 C     of fresh water at p = 0
                0489                drhoP0dS = 0. _d 0
                0490 C     of salt water at p = 0
                0491                drhoP0dS = drhoP0dS
                0492      &              + eosJMDCSw(1)
31a5949e20 Mart*0493      &              + eosJMDCSw(2)*t1
a37a13034c Mart*0494      &              + eosJMDCSw(3)*t2
                0495      &              + eosJMDCSw(4)*t3
31a5949e20 Mart*0496      &              + eosJMDCSw(5)*t3*t1
a37a13034c Mart*0497      &       + s3o2*(
                0498      &                eosJMDCSw(6)
31a5949e20 Mart*0499      &              + eosJMDCSw(7)*t1
a37a13034c Mart*0500      &              + eosJMDCSw(8)*t2
                0501      &              )
31a5949e20 Mart*0502      &              + 2*eosJMDCSw(9)*s1
a37a13034c Mart*0503 C     d(bulk modulus)/d(S)
                0504 C     of fresh water at p = 0
                0505                dKdS = 0. _d 0
                0506 C     of sea water at p = 0
                0507                dKdSSalt =
                0508      &                eosJMDCKSw(1)
31a5949e20 Mart*0509      &              + eosJMDCKSw(2)*t1
a37a13034c Mart*0510      &              + eosJMDCKSw(3)*t2
                0511      &              + eosJMDCKSw(4)*t3
                0512      &       + s3o2*( eosJMDCKSw(5)
31a5949e20 Mart*0513      &              + eosJMDCKSw(6)*t1
a37a13034c Mart*0514      &              + eosJMDCKSw(7)*t2
                0515      &              )
                0516
                0517 C     of sea water at p
de7f345328 Jean*0518                dKdSPres =
31a5949e20 Mart*0519      &           p1*( eosJMDCKP(5)
                0520      &              + eosJMDCKP(6)*t1
a37a13034c Mart*0521      &              + eosJMDCKP(7)*t2
                0522      &              )
31a5949e20 Mart*0523      &        + s3o2*p1*eosJMDCKP(8)
                0524      &      + p1*p1*( eosJMDCKP(12)
                0525      &              + eosJMDCKP(13)*t1
a37a13034c Mart*0526      &              + eosJMDCKP(14)*t2
                0527      &              )
                0528
                0529                dKdS = dKdSSalt + dKdSPres
                0530
de7f345328 Jean*0531                betaLoc(i,j) =
a37a13034c Mart*0532      &              ( bulkmod(i,j)**2*drhoP0dS
31a5949e20 Mart*0533      &              - bulkmod(i,j)*p1*drhoP0dS
                0534      &              - rhoP0(i,j)*p1*dKdS )
                0535      &              /( bulkmod(i,j) - p1 )**2
de7f345328 Jean*0536
6ac14281aa Jean*0537             ENDDO
                0538          ENDDO
                0539       ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
533b1afce7 Jean*0540          dp0 = surf_pRef - eosRefP0
6ac14281aa Jean*0541
                0542          CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0543      I        bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
6ac14281aa Jean*0544      O        locPres,
de7f345328 Jean*0545      I        myThid )
6ac14281aa Jean*0546
de7f345328 Jean*0547          CALL FIND_RHONUM(
                0548      I        iMin, iMax, jMin, jMax, locPres,
                0549      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0550      O        rhoLoc,
                0551      I        myThid )
edd57506ae Patr*0552
de7f345328 Jean*0553          CALL FIND_RHODEN(
                0554      I        iMin, iMax, jMin, jMax, locPres,
                0555      I        theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0556      O        rhoDen,
                0557      I        myThid )
24de0d2c47 Patr*0558
6ac14281aa Jean*0559          DO j=jMin,jMax
                0560             DO i=iMin,iMax
31a5949e20 Mart*0561                t1  = theta(i,j,k,bi,bj)
                0562                t2  = t1*t1
                0563                s1  = salt(i,j,k,bi,bj)
02d90fb24c Jean*0564 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0565 C     this line makes TAF create vectorizable code
                0566                sp5 = 0. _d 0
                0567 #endif
59e581f4e0 Jean*0568               IF ( s1 .GT. 0. _d 0 ) THEN
de7f345328 Jean*0569                sp5 = SQRT(s1)
59e581f4e0 Jean*0570               ELSE
                0571                s1  = 0. _d 0
                0572                sp5 = 0. _d 0
                0573               ENDIF
de7f345328 Jean*0574
6ac14281aa Jean*0575                p1   = locPres(i,j)*SItodBar
31a5949e20 Mart*0576                p1t1 = p1*t1
de7f345328 Jean*0577
                0578                dnum_dsalt = eosMDJWFnum(4)
31a5949e20 Mart*0579      &              + eosMDJWFnum(5)*t1
                0580      &              + 2.*eosMDJWFnum(6)*s1 + eosMDJWFnum(9)*p1
de7f345328 Jean*0581                dden_dsalt = eosMDJWFden(5)
                0582      &              + t1*( eosMDJWFden(6) + eosMDJWFden(7)*t2 )
31a5949e20 Mart*0583      &              + 1.5*sp5*(eosMDJWFden(8) + eosMDJWFden(9)*t2)
                0584
de7f345328 Jean*0585                betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
edd57506ae Patr*0586      &              - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )
31a5949e20 Mart*0587
6ac14281aa Jean*0588             ENDDO
                0589          ENDDO
31a5949e20 Mart*0590
ac1c5b24e7 Mart*0591       ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
533b1afce7 Jean*0592        dp0 = surf_pRef - eosRefP0
ac1c5b24e7 Mart*0593
                0594        CALL PRESSURE_FOR_EOS(
533b1afce7 Jean*0595      I      bi, bj, iMin, iMax, jMin, jMax,  kRef, dp0,
ac1c5b24e7 Mart*0596      O      locPres,
                0597      I      myThid )
                0598
                0599        CALL FIND_RHOTEOS(
                0600      I      iMin, iMax, jMin, jMax, locPres,
                0601      I      theta(1-OLx,1-OLy,k,bi,bj), salt(1-OLx,1-OLy,k,bi,bj),
                0602      O      rhoLoc, rhoDen,
                0603      I      myThid )
1eadaea85b Jean*0604
ac1c5b24e7 Mart*0605        DO j=jMin,jMax
                0606         DO i=iMin,iMax
                0607          ct      = theta(i,j,k,bi,bj)
                0608          sa      = salt(i,j,k,bi,bj)
02d90fb24c Jean*0609 #if (defined ALLOW_AUTODIFF && defined TARGET_NEC_SX)
6cedc7114a Mart*0610 C     this line makes TAF create vectorizable code
                0611          sqrtsa  = 0. _d 0
                0612 #endif
ac1c5b24e7 Mart*0613          IF ( sa .GT. 0. _d 0 ) THEN
                0614           sqrtsa = SQRT(sa)
                0615          ELSE
                0616           sa     = 0. _d 0
                0617           sqrtsa = 0. _d 0
                0618          ENDIF
                0619          p       = locPres(i,j)*SItodBar
                0620
                0621          dnum_dsalt =  teos(05) + ct*(teos(06) + teos(07)*ct)
                0622      &   + 1.5*sqrtsa*(teos(08)
                0623      &           + ct*(teos(09) + ct*(teos(10) + teos(11)*ct)))
                0624      &   + p*(teos(15) + teos(16)*ct + p*teos(20))
                0625
1eadaea85b Jean*0626          dden_dsalt = teos(26)
ac1c5b24e7 Mart*0627      &   + ct*(teos(27) + ct*(teos(28) + ct*(teos(29) + teos(30)*ct)))
1eadaea85b Jean*0628      &   + 2.*teos(36)*sa
                0629      &   + 1.5*sqrtsa*(teos(31) + ct*(teos(32) + ct*(teos(33)
ac1c5b24e7 Mart*0630      &                          + ct*(teos(34) + teos(35)*ct))))
                0631      &   + p*(teos(41) + teos(42)*ct + p*teos(46))
                0632
                0633          betaLoc(i,j) = rhoDen(i,j)*( dnum_dsalt
                0634      &        - (rhoLoc(i,j)*rhoDen(i,j))*dden_dsalt )
                0635
                0636         ENDDO
                0637        ENDDO
                0638
6ac14281aa Jean*0639       ELSE
                0640          WRITE(*,*) 'FIND_BETA: equationOfState = ',equationOfState
                0641          STOP 'FIND_BETA: "equationOfState" has illegal value'
                0642       ENDIF
fb3dc7d949 Alis*0643
de7f345328 Jean*0644       RETURN
6ac14281aa Jean*0645       END