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
701e10a905 Mart*0313        Sprac_Sref = 35.0 _d 0 / 35.16504 _d 0
                0314        I_S0       = 0.025 _d 0 *Sprac_Sref
                0315        I_Ts       = 0.025 _d 0
                0316        I_cp0      = 1.0 _d 0 / 3991.86795711963 _d 0
                0317
                0318        H00 = 61.01362420681071 _d 0   *  I_cp0
                0319        H01 = 168776.46138048015 _d 0  * (I_cp0*I_Ts)
                0320        H02 = -2735.2785605119625 _d 0 * (I_cp0*I_Ts**2)
                0321        H03 = 2574.2164453821433 _d 0  * (I_cp0*I_Ts**3)
                0322        H04 = -1536.6644434977543 _d 0 * (I_cp0*I_Ts**4)
                0323        H05 = 545.7340497931629 _d 0   * (I_cp0*I_Ts**5)
                0324        H06 = -50.91091728474331 _d 0  * (I_cp0*I_Ts**6)
                0325        H07 = -18.30489878927802 _d 0  * (I_cp0*I_Ts**7)
                0326        H20 = 268.5520265845071 _d 0   *  I_cp0
                0327        H21 = -12019.028203559312 _d 0 * (I_cp0*I_Ts)
                0328        H22 = 3734.858026725145 _d 0   * (I_cp0*I_Ts**2)
                0329        H23 = -2046.7671145057618 _d 0 * (I_cp0*I_Ts**3)
                0330        H24 = 465.28655623826234 _d 0  * (I_cp0*I_Ts**4)
                0331        H25 = -0.6370820302376359 _d 0 * (I_cp0*I_Ts**5)
                0332        H26 = -10.650848542359153 _d 0 * (I_cp0*I_Ts**6)
                0333        H30 = 937.2099110620707 _d 0   *  I_cp0
                0334        H31 = 588.1802812170108 _d 0   * (I_cp0*I_Ts)
                0335        H32 = 248.39476522971285 _d 0  * (I_cp0*I_Ts**2)
                0336        H33 = -3.871557904936333 _d 0  * (I_cp0*I_Ts**3)
                0337        H34 = -2.6268019854268356 _d 0 * (I_cp0*I_Ts**4)
                0338        H40 = -1687.914374187449 _d 0  *  I_cp0
                0339        H41 = 936.3206544460336 _d 0   * (I_cp0*I_Ts)
                0340        H42 = -942.7827304544439 _d 0  * (I_cp0*I_Ts**2)
                0341        H43 = 369.4389437509002 _d 0   * (I_cp0*I_Ts**3)
                0342        H44 = -33.83664947895248 _d 0  * (I_cp0*I_Ts**4)
                0343        H45 = -9.987880382780322 _d 0  * (I_cp0*I_Ts**5)
                0344        H50 = 246.9598888781377 _d 0   *  I_cp0
                0345        H60 = 123.59576582457964 _d 0  *  I_cp0
                0346        H70 = -48.5891069025409 _d 0   *  I_cp0
                0347
                0348        TPN00 = -1.446013646344788 _d -2
                0349        TPN10 = -3.305308995852924 _d -3 * Sprac_Sref
                0350        TPN20 =  1.062415929128982 _d -4 * Sprac_Sref**2
                0351        TPN01 =  9.477566673794488 _d -1
                0352        TPN11 =  2.166591947736613 _d -3 * Sprac_Sref
                0353        TPN02 =  3.828842955039902 _d -3
                0354        TPD10 =  6.506097115635800 _d -4 * Sprac_Sref
                0355        TPD01 =  3.830289486850898 _d -3
                0356        TPD02 =  1.247811760368034 _d -6
                0357
06bb0cec77 Jean*0358       ELSEIF( equationOfState .EQ. 'IDEALG' ) THEN
                0359
                0360       ELSE
                0361
533b1afce7 Jean*0362          WRITE(msgBuf,'(3A)') 'INI_EOS: eosType= "',
                0363      &                     equationOfState(1:iLen), '" not valid'
06bb0cec77 Jean*0364          CALL PRINT_ERROR( msgbuf, myThid )
                0365          STOP 'ABNORMAL END: S/R INI_EOS'
                0366
                0367       ENDIF
a37a13034c Mart*0368
e4e7364636 Jean*0369 C--   Check EOS initialisation:
                0370
57061cb432 Jean*0371       CALL EOS_CHECK( myThid )
e4e7364636 Jean*0372
a922e25940 Jean*0373       _END_MASTER( myThid )
06bb0cec77 Jean*0374       _BARRIER
a37a13034c Mart*0375
bbff648e65 Jean*0376       RETURN
                0377       END
a37a13034c Mart*0378
                0379 CBOP
57061cb432 Jean*0380 C     !ROUTINE: EOS_CHECK
a37a13034c Mart*0381 C     !INTERFACE:
57061cb432 Jean*0382       SUBROUTINE EOS_CHECK( myThid )
a37a13034c Mart*0383 C     !DESCRIPTION: \bv
                0384 C     *==========================================================*
57061cb432 Jean*0385 C     | SUBROUTINE EOS_CHECK
a37a13034c Mart*0386 C     | o check the equation of state.
                0387 C     *==========================================================*
                0388 C     \ev
                0389
                0390 C     !USES:
                0391
06bb0cec77 Jean*0392       IMPLICIT NONE
a37a13034c Mart*0393 #include "SIZE.h"
                0394 #include "EEPARAMS.h"
                0395 #include "PARAMS.h"
                0396 #include "EOS.h"
57061cb432 Jean*0397 #include "GRID.h"
a37a13034c Mart*0398
                0399 C     !INPUT/OUTPUT PARAMETERS:
                0400 C     == Routine arguments ==
57061cb432 Jean*0401 C     myThid ::  Number of this instance of EOS_CHECK
a37a13034c Mart*0402       INTEGER myThid
                0403
57061cb432 Jean*0404 #ifdef INCLUDE_EOS_CHECK
a37a13034c Mart*0405 C     !LOCAL VARIABLES:
                0406 C     == Local variables ==
                0407 C     bi,bj  - Loop counters
57061cb432 Jean*0408 C     i,j,k
a37a13034c Mart*0409       INTEGER bi, bj
57061cb432 Jean*0410       INTEGER iMin, iMax, jMin, jMax
06bb0cec77 Jean*0411       INTEGER  i, j, k
57061cb432 Jean*0412       _RL tFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0413       _RL sFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0414       _RL pFld   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
533b1afce7 Jean*0415       _RL rhoLoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0416       _RL bulkMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
a37a13034c Mart*0417
533b1afce7 Jean*0418       INTEGER ncheck, kcheck, k1c, k2c, iLen
ac1c5b24e7 Mart*0419       PARAMETER ( ncheck = 20 )
57061cb432 Jean*0420       _RL tLoc(ncheck), ptLoc(ncheck), sLoc(ncheck), pLoc(ncheck)
                0421       _RL rLoc(ncheck), bLoc(ncheck)
533b1afce7 Jean*0422       _RL surf_pRef_save, rhoDif
a37a13034c Mart*0423
                0424       CHARACTER*(MAX_LEN_MBUF) msgBuf
533b1afce7 Jean*0425       CHARACTER*(13) blkWrd
a37a13034c Mart*0426
57061cb432 Jean*0427       DATA tLoc
a37a13034c Mart*0428      &     /3.25905152915860 _d 0, 20.38687090048638 _d 0,
31a5949e20 Mart*0429      &     25.44820830309568 _d 0, 20.17368557065936 _d 0,
                0430      &     13.43397459640398 _d 0,
a37a13034c Mart*0431      &      5.               _d 0, 25.               _d 0,
                0432      &      5.               _d 0, 25.               _d 0,
                0433      &      5.               _d 0, 25.               _d 0,
ac1c5b24e7 Mart*0434      &      5.               _d 0, 25.               _d 0,
                0435      &      5.               _d 0, 25.               _d 0,
                0436      &      5.               _d 0, 25.               _d 0,
                0437      &      5.               _d 0, 25.               _d 0,
                0438      &      5.               _d 0/
57061cb432 Jean*0439      &     ptLoc
a37a13034c Mart*0440      &     /3.               _d 0, 20.               _d 0,
31a5949e20 Mart*0441      &     25.               _d 0, 20.               _d 0,
                0442      &     12.               _d 0,
a37a13034c Mart*0443      &      5.               _d 0, 25.               _d 0,
                0444      &      5.               _d 0, 25.               _d 0,
06bb0cec77 Jean*0445      &      4.03692566635316 _d 0, 22.84661726775120 _d 0,
ac1c5b24e7 Mart*0446      &      3.62720389416752 _d 0, 22.62420229124846 _d 0,
                0447      &     19.5130           _d 0,  3.5588           _d 0,
                0448      &      1.8157           _d 0,  1.3434           _d 0,
                0449      &      1.1583           _d 0,  1.0518           _d 0,
                0450      &      1.0073           _d 0/
57061cb432 Jean*0451      &     sLoc
a37a13034c Mart*0452      &     /35.5 _d 0, 35. _d 0,
31a5949e20 Mart*0453      &      35.0 _d 0, 20. _d 0,
                0454      &      40.0 _d 0,
a37a13034c Mart*0455      &       0.  _d 0,  0. _d 0,
                0456      &      35.  _d 0, 35. _d 0,
                0457      &       0.  _d 0,  0. _d 0,
ac1c5b24e7 Mart*0458      &      35.  _d 0, 35. _d 0,
                0459      &      34.7392 _d 0, 34.4652 _d 0,
533b1afce7 Jean*0460      &      34.7738 _d 0, 34.8435 _d 0,
ac1c5b24e7 Mart*0461      &      34.8637 _d 0, 34.8739 _d 0, 34.8776 _d 0/
57061cb432 Jean*0462      &     pLoc
a37a13034c Mart*0463      &     /300. _d 5,  200. _d 5,
31a5949e20 Mart*0464      &      200. _d 5,  100. _d 5,
                0465      &      800. _d 5,
a37a13034c Mart*0466      &        0. _d 0,    0. _d 0,
                0467      &        0. _d 0,    0. _d 0,
                0468      &     1000. _d 5, 1000. _d 5,
ac1c5b24e7 Mart*0469      &     1000. _d 5, 1000. _d 5,
                0470      &        0. _d 0, 1010. _d 4,
                0471      &     2025. _d 4, 3045. _d 4,
                0472      &     4069. _d 4, 5098. _d 4,
                0473      &     6131. _d 4/
57061cb432 Jean*0474       DATA rLoc
31a5949e20 Mart*0475      &     /1041.83267  _d 0, 1033.213387 _d 0,
                0476      &      1031.654229 _d 0, 1017.726743 _d 0,
                0477      &      1062.928258 _d 0,
                0478      &       999.96675  _d 0,  997.04796  _d 0,
                0479      &      1027.67547  _d 0, 1023.34306  _d 0,
                0480      &      1044.12802  _d 0, 1037.90204  _d 0,
ac1c5b24e7 Mart*0481      &      1069.48914  _d 0, 1062.53817  _d 0,
                0482      &      1024.571477039354932 _d 0, 1031.937207908966911 _d 0,
                0483      &      1037.002326523349893 _d 0, 1041.668901630106348 _d 0,
                0484      &      1046.179673478751511 _d 0, 1050.591667083023594 _d 0,
                0485      &      1054.901835905235885 _d 0/
57061cb432 Jean*0486      &     bLoc
a37a13034c Mart*0487      &     /   -1.00000 _d 0,    -1.00000 _d 0,
31a5949e20 Mart*0488      &         -1.00000 _d 0,    -1.00000 _d 0,
                0489      &         -1.00000 _d 0,
a37a13034c Mart*0490      &      20337.80375 _d 0, 22100.72106 _d 0,
                0491      &      22185.93358 _d 0, 23726.34949 _d 0,
                0492      &      23643.52599 _d 0, 25405.09717 _d 0,
ac1c5b24e7 Mart*0493      &      25577.49819 _d 0, 27108.94504 _d 0,
                0494      &         -1.00000 _d 0,    -1.00000 _d 0,
                0495      &         -1.00000 _d 0,    -1.00000 _d 0,
                0496      &         -1.00000 _d 0,    -1.00000 _d 0,
                0497      &         -1.00000 _d 0/
06bb0cec77 Jean*0498
533b1afce7 Jean*0499       blkWrd = '             '
a37a13034c Mart*0500       bi   = 1
                0501       bj   = 1
                0502       k    = 1
57061cb432 Jean*0503       iMin = 1
                0504       iMax = 1
                0505       jMin = 1
                0506       jMax = 1
a37a13034c Mart*0507       i    = 1
                0508       j    = 1
01c57eed7d Mart*0509       IF (       equationOfState.NE.'LINEAR'
57061cb432 Jean*0510      &     .AND. equationOfState.NE.'POLY3'
                0511      &     .AND. equationOfState.NE.'IDEALG' ) THEN
a37a13034c Mart*0512 C     check nonlinear EOS
57061cb432 Jean*0513         WRITE(msgBuf,'(A,A)')
                0514      &        'EOS_CHECK: Check the equation of state: Type ',
a37a13034c Mart*0515      &        equationOfState
57061cb432 Jean*0516         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0517      &                      SQUEEZE_RIGHT, myThid )
                0518         IF ( equationOfState.EQ.'JMD95Z' .OR.
                0519      &       equationOfState.EQ.'JMD95P' ) THEN
                0520           WRITE(msgBuf,'(A)')
                0521      &         'EOS_CHECK: check Rho values for eosType=JMD95:'
                0522           k1c = 1
                0523           k2c = 1
                0524         ELSEIF ( equationOfState.EQ.'MDJWF' ) THEN
                0525           WRITE(msgBuf,'(A)')
                0526      &         'EOS_CHECK: check Rho values for eosType=MDJWF:'
                0527           k1c = 2
                0528           k2c = 5
                0529         ELSEIF ( equationOfState.EQ.'UNESCO' ) THEN
                0530           WRITE(msgBuf,'(A)')
                0531      &         'EOS_CHECK: check Rho & K values for eosType=UNESCO:'
                0532           k1c = 6
                0533           k2c = 13
                0534         ELSEIF ( equationOfState.EQ.'TEOS10' ) THEN
                0535           WRITE(msgBuf,'(A)')
                0536      &         'EOS_CHECK: check Rho values for eosType=TEOS10:'
                0537           k1c = 14
                0538           k2c = ncheck
                0539         ELSE
                0540           WRITE(msgBuf,'(3A)') 'EOS_CHECK: Invalid eosType= ',
                0541      &          equationOfState
                0542           CALL PRINT_ERROR( msgBuf, myThid )
                0543           STOP 'ABNORMAL END: S/R EOS_CHECK'
                0544         ENDIF
                0545         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0546      &                      SQUEEZE_RIGHT, myThid )
                0547 C-    Will use "surf_pRef" to get the right pressure out of S/R PRESSURE_FOR_EOS
                0548 C     --> save current value into "surf_pRef_save"
                0549         surf_pRef_save = surf_pRef
                0550         DO kcheck = k1c, k2c
                0551
                0552           surf_pRef = eosRefP0 + pLoc(kcheck)
                0553           IF ( usingPCoords ) THEN
                0554             surf_pRef = surf_pRef - rC(k)
57061cb432 Jean*0555           ENDIF
                0556           IF ( equationOfState.NE.'UNESCO' ) THEN
                0557              tFld(i,j) = ptLoc(kcheck)
                0558           ELSE
                0559              tFld(i,j) = tLoc(kcheck)
                0560           ENDIF
                0561           sFld(i,j)    = sLoc(kcheck)
                0562           pFld(i,j)    = pLoc(kcheck)
                0563           rhoLoc(i,j)  =  0. _d 0
                0564
                0565           CALL FIND_RHO_2D(
                0566      I           iMin, iMax, jMin, jMax, k,
                0567      I           tFld, sFld,
                0568      O           rhoLoc,
                0569      I           k, bi, bj, myThid )
533b1afce7 Jean*0570           rhoDif = rhoLoc(i,j) + rhoConst - rLoc(kcheck)
57061cb432 Jean*0571
533b1afce7 Jean*0572           IF ( equationOfState.EQ.'UNESCO' ) THEN
06bb0cec77 Jean*0573             CALL FIND_BULKMOD(
57061cb432 Jean*0574      I           iMin, iMax, jMin, jMax,
                0575      I           pFld, tFld, sFld,
                0576      O           bulkMod,
                0577      I           myThid )
533b1afce7 Jean*0578 C-    Reference value
                0579             WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,A,F11.5)')
ba0b047096 Mart*0580      &           'rho(', sFld(i,j), ' g/kg,',
533b1afce7 Jean*0581      &           tFld(i,j), ' degC,',
                0582      &           pLoc(kcheck)*SItoBar, ' bar)',
                0583      &           ' = ', rLoc(kcheck), ' ', bLoc(kcheck)
                0584             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0585      &                          SQUEEZE_RIGHT, myThid )
                0586             WRITE(msgBuf,'(A,F10.5,A,F11.5,A,1PE10.3)')
                0587      &           ' rho(find_rho_2d)                 = ',
                0588      &           rhoLoc(i,j)+rhoConst, ' ', bulkMod(i,j), ' ', rhoDif
                0589             iLen = LEN( blkWrd )
                0590           ELSE
                0591 C-    Reference value
                0592             WRITE(msgBuf,'(2(A,F4.1),A,F5.0,2A,F10.5,A,F11.5)')
ba0b047096 Mart*0593      &           'rho(', sFld(i,j), ' g/kg,',
533b1afce7 Jean*0594      &           tFld(i,j), ' degC,',
                0595      &           pLoc(kcheck)*SItoBar, ' bar)',
                0596      &           ' = ', rLoc(kcheck)
                0597             CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0598      &                          SQUEEZE_RIGHT, myThid )
                0599             WRITE(msgBuf,'(A,F10.5,A,1PE10.3)')
01c57eed7d Mart*0600      &         ' rho(find_rho_2d)                 = ',
533b1afce7 Jean*0601      &         rhoLoc(i,j)+rhoConst, ' ', rhoDif
                0602             iLen = 1
                0603           ENDIF
57061cb432 Jean*0604           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0605      &                        SQUEEZE_RIGHT, myThid )
57061cb432 Jean*0606
533b1afce7 Jean*0607           surf_pRef = eosRefP0
57061cb432 Jean*0608           CALL FIND_RHO_SCALAR( tFld(i,j), sLoc(kcheck),
                0609      &         pLoc(kcheck), rhoLoc(i,j), myThid )
533b1afce7 Jean*0610           rhoDif = rhoLoc(i,j) - rLoc(kcheck)
                0611           WRITE(msgBuf,'(A,F10.5,A,1PE10.3)')
01c57eed7d Mart*0612      &         ' rho(find_rho_scalar)             = ',
533b1afce7 Jean*0613      &         rhoLoc(i,j), blkWrd(1:iLen), rhoDif
57061cb432 Jean*0614           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0615      &                        SQUEEZE_RIGHT, myThid )
57061cb432 Jean*0616
                0617         ENDDO
533b1afce7 Jean*0618 C     end check nonlinear EOS ; restore surf_pRef value
                0619         surf_pRef = surf_pRef_save
a37a13034c Mart*0620
57061cb432 Jean*0621         WRITE(msgBuf,'(A)') 'EOS_CHECK: Done'
                0622         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
533b1afce7 Jean*0623      &                      SQUEEZE_RIGHT, myThid )
a37a13034c Mart*0624
06bb0cec77 Jean*0625       ENDIF
57061cb432 Jean*0626 #endif /* INCLUDE_EOS_CHECK */
a37a13034c Mart*0627
bbff648e65 Jean*0628       RETURN
                0629       END