model/src/ini_eos.F

a37a13034c Mart*0001 #include "CPP_OPTIONS.h"
57061cb432 Jean*0002 #undef INCLUDE_EOS_CHECK
a37a13034c Mart*0003
57061cb432 Jean*0004 C--  File ini_eos.F: Routines to initialise Equation of State
                0005 C--   Contents
                0006 C--   o INI_EOS
                0007 C--   o EOS_CHECK
                0008
                0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
a37a13034c Mart*0010 CBOP
                0011 C     !ROUTINE: INI_EOS
                0012 C     !INTERFACE:
06bb0cec77 Jean*0013       SUBROUTINE INI_EOS( myThid )
a37a13034c Mart*0014 C     !DESCRIPTION: \bv
                0015 C     *==========================================================*
06bb0cec77 Jean*0016 C     | SUBROUTINE INI_EOS
a37a13034c Mart*0017 C     | o Initialise coefficients of equation of state.
                0018 C     *==========================================================*
                0019 C     \ev
                0020
                0021 C     !USES:
                0022
06bb0cec77 Jean*0023       IMPLICIT NONE
a37a13034c Mart*0024 C     == Global variables ==
                0025 #include "SIZE.h"
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "EOS.h"
                0029
                0030 C     !INPUT/OUTPUT PARAMETERS:
                0031 C     == Routine arguments ==
                0032 C     myThid -  Number of this instance of INI_CORI
                0033       INTEGER myThid
                0034
533b1afce7 Jean*0035 C     !FUNCTIONS
                0036       INTEGER  ILNBLNK
                0037       EXTERNAL ILNBLNK
                0038
a37a13034c Mart*0039 C     !LOCAL VARIABLES:
                0040 C     == Local variables ==
06bb0cec77 Jean*0041 C     i,k    :: Loop counters
533b1afce7 Jean*0042       INTEGER  i, k, iLen
a37a13034c Mart*0043       CHARACTER*(MAX_LEN_MBUF) msgBuf
06bb0cec77 Jean*0044
9669509dca Jean*0045       IF ( .NOT.fluidIsWater ) RETURN
a37a13034c Mart*0046
06bb0cec77 Jean*0047       _BARRIER
                0048       _BEGIN_MASTER(myThid)
                0049
a37a13034c Mart*0050       equationOfState = eosType
533b1afce7 Jean*0051       iLen = ILNBLNK(equationOfState)
                0052       iLen = MAX(iLen,1)
a37a13034c Mart*0053
06bb0cec77 Jean*0054       DO k = 1,6
31a5949e20 Mart*0055          eosJMDCFw(k) = 0. _d 0
06bb0cec77 Jean*0056       ENDDO
                0057       DO k = 1,9
31a5949e20 Mart*0058          eosJMDCSw(k) = 0. _d 0
06bb0cec77 Jean*0059       ENDDO
                0060       DO k = 1,5
31a5949e20 Mart*0061          eosJMDCKFw(k) = 0. _d 0
06bb0cec77 Jean*0062       ENDDO
                0063       DO k = 1,7
31a5949e20 Mart*0064          eosJMDCKSw(k) = 0. _d 0
06bb0cec77 Jean*0065       ENDDO
                0066       DO k = 1,14
31a5949e20 Mart*0067          eosJMDCKP(k) = 0. _d 0
06bb0cec77 Jean*0068       ENDDO
                0069       DO k = 0,11
31a5949e20 Mart*0070          eosMDJWFnum(k) = 0. _d 0
06bb0cec77 Jean*0071       ENDDO
                0072       DO k = 0,12
31a5949e20 Mart*0073          eosMDJWFden(k) = 0. _d 0
06bb0cec77 Jean*0074       ENDDO
ac1c5b24e7 Mart*0075       DO k = 1,48
                0076          teos(k)        = 0. _d 0
                0077       ENDDO
31a5949e20 Mart*0078
533b1afce7 Jean*0079 C     For all current EOS that explicitly depend on pressure ('UNESCO',
                0080 C     'JMD95','MDJWF' & 'TEOS10' ), the input pressure is assumed to be
                0081 C     relative to reference atmospheric pressure of 1.atm = 101325 Pa
                0082       eosRefP0 = 101325. _d 0
                0083
06bb0cec77 Jean*0084       IF ( equationOfState .EQ. 'LINEAR' ) THEN
                0085          IF ( tAlpha .EQ. UNSET_RL ) tAlpha = 2.  _d -4
                0086          IF ( sBeta  .EQ. UNSET_RL ) sBeta  = 7.4 _d -4
                0087       ELSEIF ( equationOfState .EQ. 'POLY3' ) THEN
a37a13034c Mart*0088          OPEN(37,FILE='POLY3.COEFFS',STATUS='OLD',FORM='FORMATTED')
06bb0cec77 Jean*0089          READ(37,*) k
                0090          IF (k.NE.Nr) THEN
a37a13034c Mart*0091             WRITE(msgBuf,'(A)')
                0092      &           'ini_eos: attempt to read POLY3.COEFFS failed'
533b1afce7 Jean*0093             CALL PRINT_ERROR( msgBuf, myThid )
a37a13034c Mart*0094             WRITE(msgBuf,'(A)')
                0095      &           '           because bad # of levels in data'
533b1afce7 Jean*0096             CALL PRINT_ERROR( msgBuf, myThid )
a37a13034c Mart*0097             STOP 'Bad data in POLY3.COEFFS'
                0098          ENDIF
06bb0cec77 Jean*0099          READ(37,*) (eosRefT(k),eosRefS(k),eosSig0(k),k=1,Nr)
                0100          DO k=1,Nr
                0101             READ(37,*) (eosC(i,k),i=1,9)
a37a13034c Mart*0102          ENDDO
                0103          CLOSE(37)
9b1fb10894 Mart*0104
9fb9047cf4 Jean*0105       ELSEIF ( equationOfState .EQ. 'JMD95Z'
                0106      &    .OR. equationOfState .EQ. 'JMD95P'
                0107      &    .OR. equationOfState .EQ. 'UNESCO' ) THEN
533b1afce7 Jean*0108
                0109 C     Jackett & McDougall (1995, JAOT) equation of state
a37a13034c Mart*0110 C     rho = R(salinity, potential temperature, pressure)
06bb0cec77 Jean*0111 C     pressure needs to be available (from the previous
a37a13034c Mart*0112 C     time step to linearize the problem)
533b1afce7 Jean*0113
06bb0cec77 Jean*0114          IF ( equationOfState .EQ. 'JMD95Z' .AND. usingPCoords ) THEN
                0115             WRITE(msgBuf,'(A)')
9b1fb10894 Mart*0116      &      'ini_eos: equation of state ''JMD95Z'' should not'
533b1afce7 Jean*0117             CALL PRINT_ERROR( msgBuf, myThid )
06bb0cec77 Jean*0118             WRITE(msgBuf,'(A)')
9b1fb10894 Mart*0119      &      '         be used together with pressure coordinates.'
533b1afce7 Jean*0120             CALL PRINT_ERROR( msgBuf, myThid )
06bb0cec77 Jean*0121             WRITE(msgBuf,'(A)')
9b1fb10894 Mart*0122      &      '         Use only ''JMD95P'' with ''OCEANICP''.'
533b1afce7 Jean*0123             CALL PRINT_ERROR( msgBuf, myThid )
9b1fb10894 Mart*0124             STOP 'ABNORMAL END: S/R INI_EOS'
06bb0cec77 Jean*0125          ENDIF
a37a13034c Mart*0126
                0127 C     coefficients nonlinear equation of state in pressure coordinates for
                0128 C     1. density of fresh water at p = 0
9c56ac6cb9 Jean*0129          eosJMDCFw(1) =  999.842594 _d +00
a37a13034c Mart*0130          eosJMDCFw(2) =    6.793952 _d -02
                0131          eosJMDCFw(3) = -  9.095290 _d -03
                0132          eosJMDCFw(4) =    1.001685 _d -04
                0133          eosJMDCFw(5) = -  1.120083 _d -06
                0134          eosJMDCFw(6) =    6.536332 _d -09
                0135 C     2. density of sea water at p = 0
                0136          eosJMDCSw(1) =    8.24493  _d -01
                0137          eosJMDCSw(2) = -  4.0899   _d -03
9c56ac6cb9 Jean*0138          eosJMDCSw(3) =    7.6438   _d -05
a37a13034c Mart*0139          eosJMDCSw(4) = -  8.2467   _d -07
9c56ac6cb9 Jean*0140          eosJMDCSw(5) =    5.3875   _d -09
                0141          eosJMDCSw(6) = -  5.72466  _d -03
                0142          eosJMDCSw(7) =    1.0227   _d -04
a37a13034c Mart*0143          eosJMDCSw(8) = -  1.6546   _d -06
                0144          eosJMDCSw(9) =    4.8314   _d -04
06bb0cec77 Jean*0145          IF ( equationOfState(1:5) .EQ. 'JMD95' ) THEN
a37a13034c Mart*0146 C     3. secant bulk modulus K of fresh water at p = 0
                0147             eosJMDCKFw(1) =   1.965933 _d +04
                0148             eosJMDCKFw(2) =   1.444304 _d +02
                0149             eosJMDCKFw(3) = - 1.706103 _d +00
                0150             eosJMDCKFw(4) =   9.648704 _d -03
                0151             eosJMDCKFw(5) = - 4.190253 _d -05
                0152 C     4. secant bulk modulus K of sea water at p = 0
                0153             eosJMDCKSw(1) =   5.284855 _d +01
                0154             eosJMDCKSw(2) = - 3.101089 _d -01
                0155             eosJMDCKSw(3) =   6.283263 _d -03
                0156             eosJMDCKSw(4) = - 5.084188 _d -05
                0157             eosJMDCKSw(5) =   3.886640 _d -01
                0158             eosJMDCKSw(6) =   9.085835 _d -03
                0159             eosJMDCKSw(7) = - 4.619924 _d -04
                0160 C     5. secant bulk modulus K of sea water at p
                0161             eosJMDCKP( 1) =   3.186519 _d +00
                0162             eosJMDCKP( 2) =   2.212276 _d -02
9c56ac6cb9 Jean*0163             eosJMDCKP( 3) = - 2.984642 _d -04
a37a13034c Mart*0164             eosJMDCKP( 4) =   1.956415 _d -06
                0165             eosJMDCKP( 5) =   6.704388 _d -03
                0166             eosJMDCKP( 6) = - 1.847318 _d -04
                0167             eosJMDCKP( 7) =   2.059331 _d -07
                0168             eosJMDCKP( 8) =   1.480266 _d -04
                0169             eosJMDCKP( 9) =   2.102898 _d -04
                0170             eosJMDCKP(10) = - 1.202016 _d -05
                0171             eosJMDCKP(11) =   1.394680 _d -07
                0172             eosJMDCKP(12) = - 2.040237 _d -06
                0173             eosJMDCKP(13) =   6.128773 _d -08
                0174             eosJMDCKP(14) =   6.207323 _d -10
                0175
06bb0cec77 Jean*0176          ELSEIF ( equationOfState .EQ. 'UNESCO' ) THEN
a37a13034c Mart*0177
06bb0cec77 Jean*0178             WRITE(msgBuf,'(a)')
a37a13034c Mart*0179      &           'WARNING WARNING WARNING WARNING WARNING WARNING '
06bb0cec77 Jean*0180             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0181      &                          SQUEEZE_RIGHT, myThid )
06bb0cec77 Jean*0182             WRITE(msgBuf,'(a,a)')
a37a13034c Mart*0183      &           'WARNING: using the UNESCO formula with potential ',
                0184      &           'temperature'
06bb0cec77 Jean*0185             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0186      &                          SQUEEZE_RIGHT, myThid )
06bb0cec77 Jean*0187             WRITE(msgBuf,'(a)')
a37a13034c Mart*0188      &           'WARNING: can result in density errors of up to 5%'
06bb0cec77 Jean*0189             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0190      &                          SQUEEZE_RIGHT, myThid )
06bb0cec77 Jean*0191             WRITE(msgBuf,'(a)')
60c223928f Mart*0192      &           'WARNING: (see Jackett and McDougall 1995, JAOT)'
06bb0cec77 Jean*0193             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0194      &                          SQUEEZE_RIGHT, myThid )
06bb0cec77 Jean*0195             WRITE(msgBuf,'(a)')
a37a13034c Mart*0196      &           'WARNING WARNING WARNING WARNING WARNING WARNING '
06bb0cec77 Jean*0197             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0198      &                          SQUEEZE_RIGHT, myThid )
a37a13034c Mart*0199
                0200 C     3. secant bulk modulus K of fresh water at p = 0
                0201             eosJMDCKFw(1) =   1.965221 _d +04
                0202             eosJMDCKFw(2) =   1.484206 _d +02
                0203             eosJMDCKFw(3) = - 2.327105 _d +00
                0204             eosJMDCKFw(4) =   1.360477 _d -02
                0205             eosJMDCKFw(5) = - 5.155288 _d -05
                0206 C     4. secant bulk modulus K of sea water at p = 0
                0207             eosJMDCKSw(1) =   5.46746  _d +01
9c56ac6cb9 Jean*0208             eosJMDCKSw(2) = - 0.603459 _d +00
a37a13034c Mart*0209             eosJMDCKSw(3) =   1.09987  _d -02
                0210             eosJMDCKSw(4) = - 6.1670   _d -05
                0211             eosJMDCKSw(5) =   7.944    _d -02
                0212             eosJMDCKSw(6) =   1.6483   _d -02
                0213             eosJMDCKSw(7) = - 5.3009   _d -04
                0214 C     5. secant bulk modulus K of sea water at p
                0215             eosJMDCKP( 1) =   3.239908 _d +00
                0216             eosJMDCKP( 2) =   1.43713  _d -03
                0217             eosJMDCKP( 3) =   1.16092  _d -04
                0218             eosJMDCKP( 4) = - 5.77905  _d -07
                0219             eosJMDCKP( 5) =   2.2838   _d -03
                0220             eosJMDCKP( 6) = - 1.0981   _d -05
                0221             eosJMDCKP( 7) = - 1.6078   _d -06
                0222             eosJMDCKP( 8) =   1.91075  _d -04
                0223             eosJMDCKP( 9) =   8.50935  _d -05
                0224             eosJMDCKP(10) = - 6.12293  _d -06
                0225             eosJMDCKP(11) =   5.2787   _d -08
9c56ac6cb9 Jean*0226             eosJMDCKP(12) = - 9.9348   _d -07
a37a13034c Mart*0227             eosJMDCKP(13) =   2.0816   _d -08
                0228             eosJMDCKP(14) =   9.1697   _d -10
06bb0cec77 Jean*0229          ELSE
a37a13034c Mart*0230             STOP 'INI_EOS: We should never reach this point!'
06bb0cec77 Jean*0231          ENDIF
a37a13034c Mart*0232
06bb0cec77 Jean*0233       ELSEIF ( equationOfState .EQ. 'MDJWF' ) THEN
31a5949e20 Mart*0234
71a3231ad9 Mart*0235          eosMDJWFnum( 0) =  9.99843699 _d +02
                0236          eosMDJWFnum( 1) =  7.35212840 _d +00
                0237          eosMDJWFnum( 2) = -5.45928211 _d -02
                0238          eosMDJWFnum( 3) =  3.98476704 _d -04
                0239          eosMDJWFnum( 4) =  2.96938239 _d +00
                0240          eosMDJWFnum( 5) = -7.23268813 _d -03
                0241          eosMDJWFnum( 6) =  2.12382341 _d -03
                0242          eosMDJWFnum( 7) =  1.04004591 _d -02
                0243          eosMDJWFnum( 8) =  1.03970529 _d -07
                0244          eosMDJWFnum( 9) =  5.18761880 _d -06
                0245          eosMDJWFnum(10) = -3.24041825 _d -08
                0246          eosMDJWFnum(11) = -1.23869360 _d -11
06bb0cec77 Jean*0247
9c56ac6cb9 Jean*0248          eosMDJWFden( 0) =  1.00000000 _d +00
71a3231ad9 Mart*0249          eosMDJWFden( 1) =  7.28606739 _d -03
9c56ac6cb9 Jean*0250          eosMDJWFden( 2) = -4.60835542 _d -05
71a3231ad9 Mart*0251          eosMDJWFden( 3) =  3.68390573 _d -07
                0252          eosMDJWFden( 4) =  1.80809186 _d -10
                0253          eosMDJWFden( 5) =  2.14691708 _d -03
                0254          eosMDJWFden( 6) = -9.27062484 _d -06
                0255          eosMDJWFden( 7) = -1.78343643 _d -10
                0256          eosMDJWFden( 8) =  4.76534122 _d -06
                0257          eosMDJWFden( 9) =  1.63410736 _d -09
                0258          eosMDJWFden(10) =  5.30848875 _d -06
                0259          eosMDJWFden(11) = -3.03175128 _d -16
                0260          eosMDJWFden(12) = -1.27934137 _d -17
a37a13034c Mart*0261
ac1c5b24e7 Mart*0262       ELSEIF ( equationOfState .EQ. 'TEOS10' ) THEN
                0263
                0264        teos(01) =  9.998420897506056 _d +02
                0265        teos(02) =  2.839940833161907 _d  00
                0266        teos(03) = -3.147759265588511 _d -02
                0267        teos(04) =  1.181805545074306 _d -03
                0268        teos(05) = -6.698001071123802 _d  00
                0269        teos(06) = -2.986498947203215 _d -02
                0270        teos(07) =  2.327859407479162 _d -04
                0271        teos(08) = -3.988822378968490 _d -02
                0272        teos(09) =  5.095422573880500 _d -04
                0273        teos(10) = -1.426984671633621 _d -05
                0274        teos(11) =  1.645039373682922 _d -07
                0275        teos(12) = -2.233269627352527 _d -02
                0276        teos(13) = -3.436090079851880 _d -04
                0277        teos(14) =  3.726050720345733 _d -06
                0278        teos(15) = -1.806789763745328 _d -04
                0279        teos(16) =  6.876837219536232 _d -07
                0280        teos(17) = -3.087032500374211 _d -07
                0281        teos(18) = -1.988366587925593 _d -08
                0282        teos(19) = -1.061519070296458 _d -11
                0283        teos(20) =  1.550932729220080 _d -10
                0284        teos(21) =  1.000000000000000 _d  00
                0285        teos(22) =  2.775927747785646 _d -03
                0286        teos(23) = -2.349607444135925 _d -05
                0287        teos(24) =  1.119513357486743 _d -06
                0288        teos(25) =  6.743689325042773 _d -10
                0289        teos(26) = -7.521448093615448 _d -03
                0290        teos(27) = -2.764306979894411 _d -05
                0291        teos(28) =  1.262937315098546 _d -07
                0292        teos(29) =  9.527875081696435 _d -10
                0293        teos(30) = -1.811147201949891 _d -11
                0294        teos(31) = -3.303308871386421 _d -05
                0295        teos(32) =  3.801564588876298 _d -07
                0296        teos(33) = -7.672876869259043 _d -09
                0297        teos(34) = -4.634182341116144 _d -11
                0298        teos(35) =  2.681097235569143 _d -12
                0299        teos(36) =  5.419326551148740 _d -06
                0300        teos(37) = -2.742185394906099 _d -05
                0301        teos(38) = -3.212746477974189 _d -07
                0302        teos(39) =  3.191413910561627 _d -09
                0303        teos(40) = -1.931012931541776 _d -12
                0304        teos(41) = -1.105097577149576 _d -07
                0305        teos(42) =  6.211426728363857 _d -10
                0306        teos(43) = -1.119011592875110 _d -10
                0307        teos(44) = -1.941660213148725 _d -11
                0308        teos(45) = -1.864826425365600 _d -14
                0309        teos(46) =  1.119522344879478 _d -14
                0310        teos(47) = -1.200507748551599 _d -15
533b1afce7 Jean*0311        teos(48) =  6.057902487546866 _d -17
ac1c5b24e7 Mart*0312
06bb0cec77 Jean*0313       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
                0314
                0315       ELSE
                0316
533b1afce7 Jean*0317          WRITE(msgBuf,'(3A)') 'INI_EOS: eosType= "',
                0318      &                     equationOfState(1:iLen), '" not valid'
06bb0cec77 Jean*0319          CALL PRINT_ERROR( msgbuf, myThid )
                0320          STOP 'ABNORMAL END: S/R INI_EOS'
                0321
                0322       ENDIF
a37a13034c Mart*0323
e4e7364636 Jean*0324 C--   Check EOS initialisation:
                0325
57061cb432 Jean*0326       CALL EOS_CHECK( myThid )
e4e7364636 Jean*0327
a922e25940 Jean*0328       _END_MASTER( myThid )
06bb0cec77 Jean*0329       _BARRIER
a37a13034c Mart*0330
bbff648e65 Jean*0331       RETURN
                0332       END
a37a13034c Mart*0333
                0334 CBOP
57061cb432 Jean*0335 C     !ROUTINE: EOS_CHECK
a37a13034c Mart*0336 C     !INTERFACE:
57061cb432 Jean*0337       SUBROUTINE EOS_CHECK( myThid )
a37a13034c Mart*0338 C     !DESCRIPTION: \bv
                0339 C     *==========================================================*
57061cb432 Jean*0340 C     | SUBROUTINE EOS_CHECK
a37a13034c Mart*0341 C     | o check the equation of state.
                0342 C     *==========================================================*
                0343 C     \ev
                0344
                0345 C     !USES:
                0346
06bb0cec77 Jean*0347       IMPLICIT NONE
a37a13034c Mart*0348 #include "SIZE.h"
                0349 #include "EEPARAMS.h"
                0350 #include "PARAMS.h"
                0351 #include "EOS.h"
57061cb432 Jean*0352 #include "GRID.h"
a37a13034c Mart*0353
                0354 C     !INPUT/OUTPUT PARAMETERS:
                0355 C     == Routine arguments ==
57061cb432 Jean*0356 C     myThid ::  Number of this instance of EOS_CHECK
a37a13034c Mart*0357       INTEGER myThid
                0358
57061cb432 Jean*0359 #ifdef INCLUDE_EOS_CHECK
a37a13034c Mart*0360 C     !LOCAL VARIABLES:
                0361 C     == Local variables ==
                0362 C     bi,bj  - Loop counters
57061cb432 Jean*0363 C     i,j,k
a37a13034c Mart*0364       INTEGER bi, bj
57061cb432 Jean*0365       INTEGER iMin, iMax, jMin, jMax
06bb0cec77 Jean*0366       INTEGER  i, j, k
57061cb432 Jean*0367       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0368       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0369       _RL pFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
533b1afce7 Jean*0370       _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0371       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a37a13034c Mart*0372
533b1afce7 Jean*0373       INTEGER ncheck, kcheck, k1c, k2c, iLen
ac1c5b24e7 Mart*0374       PARAMETER ( ncheck = 20 )
57061cb432 Jean*0375       _RL tLoc(ncheck), ptLoc(ncheck), sLoc(ncheck), pLoc(ncheck)
                0376       _RL rLoc(ncheck), bLoc(ncheck)
533b1afce7 Jean*0377       _RL surf_pRef_save, rhoDif
a37a13034c Mart*0378
                0379       CHARACTER*(MAX_LEN_MBUF) msgBuf
533b1afce7 Jean*0380       CHARACTER*(13) blkWrd
a37a13034c Mart*0381
57061cb432 Jean*0382       DATA tLoc
a37a13034c Mart*0383      &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,
31a5949e20 Mart*0384      &     25.44820830309568 _d 0, 20.17368557065936 _d 0,
                0385      &     13.43397459640398 _d 0,
a37a13034c Mart*0386      &      5.               _d 0, 25.               _d 0,
                0387      &      5.               _d 0, 25.               _d 0,
                0388      &      5.               _d 0, 25.               _d 0,
ac1c5b24e7 Mart*0389      &      5.               _d 0, 25.               _d 0,
                0390      &      5.               _d 0, 25.               _d 0,
                0391      &      5.               _d 0, 25.               _d 0,
                0392      &      5.               _d 0, 25.               _d 0,
                0393      &      5.               _d 0/
57061cb432 Jean*0394      &     ptLoc
a37a13034c Mart*0395      &     /3.               _d 0, 20.               _d 0,
31a5949e20 Mart*0396      &     25.               _d 0, 20.               _d 0,
                0397      &     12.               _d 0,
a37a13034c Mart*0398      &      5.               _d 0, 25.               _d 0,
                0399      &      5.               _d 0, 25.               _d 0,
06bb0cec77 Jean*0400      &      4.03692566635316 _d 0, 22.84661726775120 _d 0,
ac1c5b24e7 Mart*0401      &      3.62720389416752 _d 0, 22.62420229124846 _d 0,
                0402      &     19.5130           _d 0,  3.5588           _d 0,
                0403      &      1.8157           _d 0,  1.3434           _d 0,
                0404      &      1.1583           _d 0,  1.0518           _d 0,
                0405      &      1.0073           _d 0/
57061cb432 Jean*0406      &     sLoc
a37a13034c Mart*0407      &     /35.5 _d 0, 35. _d 0,
31a5949e20 Mart*0408      &      35.0 _d 0, 20. _d 0,
                0409      &      40.0 _d 0,
a37a13034c Mart*0410      &       0.  _d 0,  0. _d 0,
                0411      &      35.  _d 0, 35. _d 0,
                0412      &       0.  _d 0,  0. _d 0,
ac1c5b24e7 Mart*0413      &      35.  _d 0, 35. _d 0,
                0414      &      34.7392 _d 0, 34.4652 _d 0,
533b1afce7 Jean*0415      &      34.7738 _d 0, 34.8435 _d 0,
ac1c5b24e7 Mart*0416      &      34.8637 _d 0, 34.8739 _d 0, 34.8776 _d 0/
57061cb432 Jean*0417      &     pLoc
a37a13034c Mart*0418      &     /300. _d 5,  200. _d 5,
31a5949e20 Mart*0419      &      200. _d 5,  100. _d 5,
                0420      &      800. _d 5,
a37a13034c Mart*0421      &        0. _d 0,    0. _d 0,
                0422      &        0. _d 0,    0. _d 0,
                0423      &     1000. _d 5, 1000. _d 5,
ac1c5b24e7 Mart*0424      &     1000. _d 5, 1000. _d 5,
                0425      &        0. _d 0, 1010. _d 4,
                0426      &     2025. _d 4, 3045. _d 4,
                0427      &     4069. _d 4, 5098. _d 4,
                0428      &     6131. _d 4/
57061cb432 Jean*0429       DATA rLoc
31a5949e20 Mart*0430      &     /1041.83267  _d 0, 1033.213387 _d 0,
                0431      &      1031.654229 _d 0, 1017.726743 _d 0,
                0432      &      1062.928258 _d 0,
                0433      &       999.96675  _d 0,  997.04796  _d 0,
                0434      &      1027.67547  _d 0, 1023.34306  _d 0,
                0435      &      1044.12802  _d 0, 1037.90204  _d 0,
ac1c5b24e7 Mart*0436      &      1069.48914  _d 0, 1062.53817  _d 0,
                0437      &      1024.571477039354932 _d 0, 1031.937207908966911 _d 0,
                0438      &      1037.002326523349893 _d 0, 1041.668901630106348 _d 0,
                0439      &      1046.179673478751511 _d 0, 1050.591667083023594 _d 0,
                0440      &      1054.901835905235885 _d 0/
57061cb432 Jean*0441      &     bLoc
a37a13034c Mart*0442      &     /   -1.00000 _d 0,    -1.00000 _d 0,
31a5949e20 Mart*0443      &         -1.00000 _d 0,    -1.00000 _d 0,
                0444      &         -1.00000 _d 0,
a37a13034c Mart*0445      &      20337.80375 _d 0, 22100.72106 _d 0,
                0446      &      22185.93358 _d 0, 23726.34949 _d 0,
                0447      &      23643.52599 _d 0, 25405.09717 _d 0,
ac1c5b24e7 Mart*0448      &      25577.49819 _d 0, 27108.94504 _d 0,
                0449      &         -1.00000 _d 0,    -1.00000 _d 0,
                0450      &         -1.00000 _d 0,    -1.00000 _d 0,
                0451      &         -1.00000 _d 0,    -1.00000 _d 0,
                0452      &         -1.00000 _d 0/
06bb0cec77 Jean*0453
533b1afce7 Jean*0454       blkWrd = '             '
a37a13034c Mart*0455       bi   = 1
                0456       bj   = 1
                0457       k    = 1
57061cb432 Jean*0458       iMin = 1
                0459       iMax = 1
                0460       jMin = 1
                0461       jMax = 1
a37a13034c Mart*0462       i    = 1
                0463       j    = 1
01c57eed7d Mart*0464       IF (       equationOfState.NE.'LINEAR'
57061cb432 Jean*0465      &     .AND. equationOfState.NE.'POLY3'
                0466      &     .AND. equationOfState.NE.'IDEALG' ) THEN
a37a13034c Mart*0467 C     check nonlinear EOS
57061cb432 Jean*0468         WRITE(msgBuf,'(A,A)')
                0469      &        'EOS_CHECK: Check the equation of state: Type ',
a37a13034c Mart*0470      &        equationOfState
57061cb432 Jean*0471         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0472      &                      SQUEEZE_RIGHT, myThid )
                0473         IF ( equationOfState.EQ.'JMD95Z' .OR.
                0474      &       equationOfState.EQ.'JMD95P' ) THEN
                0475           WRITE(msgBuf,'(A)')
                0476      &         'EOS_CHECK: check Rho values for eosType=JMD95:'
                0477           k1c = 1
                0478           k2c = 1
                0479         ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
                0480           WRITE(msgBuf,'(A)')
                0481      &         'EOS_CHECK: check Rho values for eosType=MDJWF:'
                0482           k1c = 2
                0483           k2c = 5
                0484         ELSEIF ( equationOfState.EQ.'UNESCO' ) THEN
                0485           WRITE(msgBuf,'(A)')
                0486      &         'EOS_CHECK: check Rho & K values for eosType=UNESCO:'
                0487           k1c = 6
                0488           k2c = 13
                0489         ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
                0490           WRITE(msgBuf,'(A)')
                0491      &         'EOS_CHECK: check Rho values for eosType=TEOS10:'
                0492           k1c = 14
                0493           k2c = ncheck
                0494         ELSE
                0495           WRITE(msgBuf,'(3A)') 'EOS_CHECK: Invalid eosType= ',
                0496      &          equationOfState
                0497           CALL PRINT_ERROR( msgBuf, myThid )
                0498           STOP 'ABNORMAL END: S/R EOS_CHECK'
                0499         ENDIF
                0500         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0501      &                      SQUEEZE_RIGHT, myThid )
                0502 C-    Will use "surf_pRef" to get the right pressure out of S/R PRESSURE_FOR_EOS
                0503 C     --> save current value into "surf_pRef_save"
                0504         surf_pRef_save = surf_pRef
                0505         DO kcheck = k1c, k2c
                0506
                0507           surf_pRef = eosRefP0 + pLoc(kcheck)
                0508           IF ( usingPCoords ) THEN
                0509             surf_pRef = surf_pRef - rC(k)
57061cb432 Jean*0510           ENDIF
                0511           IF ( equationOfState.NE.'UNESCO' ) THEN
                0512              tFld(i,j) = ptLoc(kcheck)
                0513           ELSE
                0514              tFld(i,j) = tLoc(kcheck)
                0515           ENDIF
                0516           sFld(i,j)    = sLoc(kcheck)
                0517           pFld(i,j)    = pLoc(kcheck)
                0518           rhoLoc(i,j)  =  0. _d 0
                0519
                0520           CALL FIND_RHO_2D(
                0521      I           iMin, iMax, jMin, jMax, k,
                0522      I           tFld, sFld,
                0523      O           rhoLoc,
                0524      I           k, bi, bj, myThid )
533b1afce7 Jean*0525           rhoDif = rhoLoc(i,j) + rhoConst - rLoc(kcheck)
57061cb432 Jean*0526
533b1afce7 Jean*0527           IF ( equationOfState.EQ.'UNESCO' ) THEN
06bb0cec77 Jean*0528             CALL FIND_BULKMOD(
57061cb432 Jean*0529      I           iMin, iMax, jMin, jMax,
                0530      I           pFld, tFld, sFld,
                0531      O           bulkMod,
                0532      I           myThid )
533b1afce7 Jean*0533 C-    Reference value
                0534             WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,A,F11.5)')
ba0b047096 Mart*0535      &           'rho(', sFld(i,j), ' g/kg,',
533b1afce7 Jean*0536      &           tFld(i,j), ' degC,',
                0537      &           pLoc(kcheck)*SItoBar, ' bar)',
                0538      &           ' = ', rLoc(kcheck), ' ', bLoc(kcheck)
                0539             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0540      &                          SQUEEZE_RIGHT, myThid )
                0541             WRITE(msgBuf,'(A,F10.5,A,F11.5,A,1PE10.3)')
                0542      &           ' rho(find_rho_2d)                 = ',
                0543      &           rhoLoc(i,j)+rhoConst, ' ', bulkMod(i,j), ' ', rhoDif
                0544             iLen = LEN( blkWrd )
                0545           ELSE
                0546 C-    Reference value
                0547             WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,A,F11.5)')
ba0b047096 Mart*0548      &           'rho(', sFld(i,j), ' g/kg,',
533b1afce7 Jean*0549      &           tFld(i,j), ' degC,',
                0550      &           pLoc(kcheck)*SItoBar, ' bar)',
                0551      &           ' = ', rLoc(kcheck)
                0552             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0553      &                          SQUEEZE_RIGHT, myThid )
                0554             WRITE(msgBuf,'(A,F10.5,A,1PE10.3)')
01c57eed7d Mart*0555      &         ' rho(find_rho_2d)                 = ',
533b1afce7 Jean*0556      &         rhoLoc(i,j)+rhoConst, ' ', rhoDif
                0557             iLen = 1
                0558           ENDIF
57061cb432 Jean*0559           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0560      &                        SQUEEZE_RIGHT, myThid )
57061cb432 Jean*0561
533b1afce7 Jean*0562           surf_pRef = eosRefP0
57061cb432 Jean*0563           CALL FIND_RHO_SCALAR( tFld(i,j), sLoc(kcheck),
                0564      &         pLoc(kcheck), rhoLoc(i,j), myThid )
533b1afce7 Jean*0565           rhoDif = rhoLoc(i,j) - rLoc(kcheck)
                0566           WRITE(msgBuf,'(A,F10.5,A,1PE10.3)')
01c57eed7d Mart*0567      &         ' rho(find_rho_scalar)             = ',
533b1afce7 Jean*0568      &         rhoLoc(i,j), blkWrd(1:iLen), rhoDif
57061cb432 Jean*0569           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0570      &                        SQUEEZE_RIGHT, myThid )
57061cb432 Jean*0571
                0572         ENDDO
533b1afce7 Jean*0573 C     end check nonlinear EOS ; restore surf_pRef value
                0574         surf_pRef = surf_pRef_save
a37a13034c Mart*0575
57061cb432 Jean*0576         WRITE(msgBuf,'(A)') 'EOS_CHECK: Done'
                0577         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0578      &                      SQUEEZE_RIGHT, myThid )
a37a13034c Mart*0579
06bb0cec77 Jean*0580       ENDIF
57061cb432 Jean*0581 #endif /* INCLUDE_EOS_CHECK */
a37a13034c Mart*0582
bbff648e65 Jean*0583       RETURN
                0584       END