model/src/seawater.F

701e10a905 Mart*0001 #include "PACKAGES_CONFIG.h"
603addf83b Mart*0002 #include "CPP_OPTIONS.h"
701e10a905 Mart*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
7c28a79aa8 Jean*0006
                0007 C--  File seawater.F: routines that compute quantities related to seawater.
                0008 C--   Contents
623b648aa3 Mart*0009 C     Seawater (SW) librabry routines
7c28a79aa8 Jean*0010 C--   o SW_PTMP: function to compute potential temperature
22cb42e69a Jean*0011 C--   o SW_TEMP: function to compute in-situ temperature from pot. temp.
                0012 C--   o SW_ADTG: function to compute adiabatic temperature gradient
                0013 C--              (used by both SW_PTMP & SW_TEMP)
701e10a905 Mart*0014 C     TEOS10 routines (renamed and modified from MOM6 implementation)
                0015 C--   o CONVERT_CT2PT: S/R to convert conservative to potential temperature
                0016 C--   o CONVERT_PT2CT: S/R to convert potential to conservative temperature
                0017 C--                    (used by CONVERT_CT2PT)
7c28a79aa8 Jean*0018
22cb42e69a Jean*0019 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0020
                0021 CBOP
                0022 C     !ROUTINE: SW_PTMP
                0023 C     !INTERFACE:
                0024       _RL FUNCTION SW_PTMP  (S,T,P,PR)
603addf83b Mart*0025
                0026 C     !DESCRIPTION: \bv
22cb42e69a Jean*0027 C     *=============================================================*
                0028 C     | S/R  SW_PTMP
                0029 C     | o compute potential temperature as per UNESCO 1983 report.
                0030 C     *=============================================================*
701e10a905 Mart*0031 C     \ev
22cb42e69a Jean*0032 C     started:
                0033 C              Armin Koehl akoehl@ucsd.edu
                0034 C
                0035 C     ==================================================================
                0036 C     SUBROUTINE SW_PTMP
                0037 C     ==================================================================
ba0b047096 Mart*0038 C     S  :: salinity    [         (PSS-78) ]
22cb42e69a Jean*0039 C     T  :: temperature [degree C (IPTS-68)]
0e6a4460e5 Ivan*0040 C     P  :: pressure    [dbar]
                0041 C     PR :: Reference pressure  [dbar]
623b648aa3 Mart*0042 C     \ev
603addf83b Mart*0043
                0044 C     !USES:
                0045       IMPLICIT NONE
                0046
                0047 C     !INPUT/OUTPUT PARAMETERS:
                0048       _RL S,T,P,PR
                0049
22cb42e69a Jean*0050 C     !FUNCTIONS:
701e10a905 Mart*0051       _RL SW_ADTG
                0052       EXTERNAL SW_ADTG
22cb42e69a Jean*0053
                0054 C     !LOCAL VARIABLES
603addf83b Mart*0055       _RL del_P ,del_th, th, q
                0056       _RL onehalf, two, three
22cb42e69a Jean*0057       PARAMETER ( onehalf = 0.5 _d 0, two = 2. _d 0, three = 3. _d 0 )
                0058 CEOP
603addf83b Mart*0059
22cb42e69a Jean*0060 C theta1
603addf83b Mart*0061       del_P  = PR - P
701e10a905 Mart*0062       del_th = del_P*SW_ADTG(S,T,P)
603addf83b Mart*0063       th     = T + onehalf*del_th
                0064       q      = del_th
22cb42e69a Jean*0065 C theta2
701e10a905 Mart*0066       del_th = del_P*SW_ADTG(S,th,P+onehalf*del_P)
603addf83b Mart*0067
                0068       th     = th + (1 - 1/sqrt(two))*(del_th - q)
                0069       q      = (two-sqrt(two))*del_th + (-two+three/sqrt(two))*q
                0070
22cb42e69a Jean*0071 C theta3
701e10a905 Mart*0072       del_th = del_P*SW_ADTG(S,th,P+onehalf*del_P)
603addf83b Mart*0073       th     = th + (1 + 1/sqrt(two))*(del_th - q)
                0074       q      = (two + sqrt(two))*del_th + (-two-three/sqrt(two))*q
                0075
22cb42e69a Jean*0076 C theta4
701e10a905 Mart*0077       del_th = del_P*SW_ADTG(S,th,P+del_P)
603addf83b Mart*0078       SW_PTMP     = th + (del_th - two*q)/(two*three)
                0079
22cb42e69a Jean*0080       RETURN
                0081       END
                0082
                0083 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
603addf83b Mart*0084
                0085 CBOP
                0086 C     !ROUTINE: SW_TEMP
                0087 C     !INTERFACE:
22cb42e69a Jean*0088       _RL FUNCTION SW_TEMP( S, T, P, PR )
603addf83b Mart*0089 C     !DESCRIPTION: \bv
                0090 C     *=============================================================*
                0091 C     | S/R  SW_TEMP
7c28a79aa8 Jean*0092 C     | o compute in-situ temperature from potential temperature
603addf83b Mart*0093 C     *=============================================================*
                0094 C
                0095 C     REFERENCES:
                0096 C     Fofonoff, P. and Millard, R.C. Jr
                0097 C     Unesco 1983. Algorithms for computation of fundamental properties of
                0098 C     seawater, 1983. _Unesco Tech. Pap. in Mar. Sci._, No. 44, 53 pp.
                0099 C     Eqn.(31) p.39
7c28a79aa8 Jean*0100 C
603addf83b Mart*0101 C     Bryden, H. 1973.
ff02675122 Jean*0102 C     New Polynomials for thermal expansion, adiabatic temperature gradient
                0103 C     and potential temperature of sea water.
603addf83b Mart*0104 C     DEEP-SEA RES., 1973, Vol20,401-408.
623b648aa3 Mart*0105 C     \ev
603addf83b Mart*0106
                0107 C     !USES:
                0108       IMPLICIT NONE
                0109 C     === Global variables ===
7c28a79aa8 Jean*0110
603addf83b Mart*0111 C     !INPUT/OUTPUT PARAMETERS:
                0112 C     === Routine arguments ===
22cb42e69a Jean*0113 C     S      :: salinity
                0114 C     T      :: potential temperature
                0115 C     P      :: pressure
                0116 C     PR     :: reference pressure
                0117       _RL  S, T, P, PR
603addf83b Mart*0118 CEOP
                0119
22cb42e69a Jean*0120 C     !FUNCTIONS:
701e10a905 Mart*0121       _RL SW_PTMP
                0122       EXTERNAL SW_PTMP
603addf83b Mart*0123
701e10a905 Mart*0124       SW_TEMP = SW_PTMP (S,T,PR,P)
603addf83b Mart*0125
                0126       RETURN
                0127       END
                0128
22cb42e69a Jean*0129 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
603addf83b Mart*0130
22cb42e69a Jean*0131 CBOP
                0132 C     !ROUTINE: SW_ADTG
                0133 C     !INTERFACE:
7c28a79aa8 Jean*0134       _RL FUNCTION SW_ADTG  (S,T,P)
603addf83b Mart*0135
22cb42e69a Jean*0136 C     !DESCRIPTION: \bv
                0137 C     *=============================================================*
                0138 C     | S/R  SW_ADTG
                0139 C     | o compute adiabatic temperature gradient as per UNESCO 1983 routines.
                0140 C     *=============================================================*
                0141 C
                0142 C     started:
                0143 C              Armin Koehl akoehl@ucsd.edu
623b648aa3 Mart*0144 C     \ev
22cb42e69a Jean*0145
                0146 C     !USES:
                0147       IMPLICIT NONE
                0148
                0149 C     !INPUT/OUTPUT PARAMETERS:
603addf83b Mart*0150       _RL S,T,P
22cb42e69a Jean*0151
                0152 C     !LOCAL VARIABLES:
                0153       _RL a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,d0,d1,e0,e1,e2
603addf83b Mart*0154       _RL sref
22cb42e69a Jean*0155 CEOP
603addf83b Mart*0156
                0157       sref = 35. _d 0
                0158       a0 =  3.5803 _d -5
                0159       a1 = +8.5258 _d -6
                0160       a2 = -6.836 _d -8
                0161       a3 =  6.6228 _d -10
                0162
                0163       b0 = +1.8932 _d -6
                0164       b1 = -4.2393 _d -8
                0165
                0166       c0 = +1.8741 _d -8
                0167       c1 = -6.7795 _d -10
                0168       c2 = +8.733 _d -12
                0169       c3 = -5.4481 _d -14
                0170
                0171       d0 = -1.1351 _d -10
                0172       d1 =  2.7759 _d -12
                0173
                0174       e0 = -4.6206 _d -13
                0175       e1 = +1.8676 _d -14
                0176       e2 = -2.1687 _d -16
                0177
                0178       SW_ADTG =      a0 + (a1 + (a2 + a3*T)*T)*T
                0179      &     + (b0 + b1*T)*(S-sref)
                0180      &     + ( (c0 + (c1 + (c2 + c3*T)*T)*T) + (d0 + d1*T)*(S-sref) )*P
                0181      &     + (  e0 + (e1 + e2*T)*T )*P*P
22cb42e69a Jean*0182
                0183       RETURN
                0184       END
701e10a905 Mart*0185
                0186 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0187
                0188 CBOP
                0189 C !ROUTINE: CONVERT_PT2CT
                0190
                0191 C !INTERFACE:
                0192       SUBROUTINE CONVERT_PT2CT(
                0193      I     Tp, Sa,
                0194      O     Tc,
                0195      I     myTime, myIter, myThid )
                0196 C     !DESCRIPTION:
                0197 C     Convert input potential temperature (degC) and absolute salinity
                0198 C     (g kg-1) to returned conservative temperature (degC) using the
                0199 C     polynomial expressions from TEOS-10.
                0200
                0201 C     !USES:
                0202       IMPLICIT NONE
                0203 C     == Global variables ===
                0204 #include "SIZE.h"
                0205 #include "EOS.h"
                0206
                0207 C     !INPUT/OUTPUT PARAMETERS:
                0208 C     Tc        :: Conservative temperature (degC)
                0209 C     Sa        :: Absolute salinity        (g kg-1)
                0210 C     T         :: potential temperature (degC)
                0211 C     myTime    :: Current time in simulation
                0212 C     myIter    :: Current iteration number
                0213 C     myThid    :: my Thread Id number
                0214       _RL Tc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0215       _RL Sa(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0216       _RL Tp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0217       _RL myTime
                0218       INTEGER myIter, myThid
                0219 CEOP
                0220
                0221 C     !LOCAL VARIABLES:
                0222       INTEGER i,j
                0223 C     Absolute salinity normalized by a plausible salinity range and its
                0224 C     square root (nondim)
                0225       _RL x2, x
                0226 C
                0227       _RL     zeRL
                0228       PARAMETER ( zeRL = 0.0 _d 0 )
                0229
                0230       DO j=1-OLy,sNy+OLy
                0231        DO i=1-OLx,sNx+OLx
                0232         x2 = MAX(I_S0 * Sa(i,j), zeRL)
                0233 #ifdef ALLOW_AUTODIFF
                0234         x = 0. _d 0
                0235         IF ( x2 .GT. 0. _d 0 ) x = SQRT(x2)
                0236 #else
                0237         x = SQRT(x2)
                0238 #endif
                0239         Tc(i,j) = H00 + (Tp(i,j)*(H01 + Tp(i,j)*(H02 +  Tp(i,j)*(H03
                0240      &       +  Tp(i,j)*(H04  + Tp(i,j)*(H05
                0241      &       + Tp(i,j)*(H06 + Tp(i,j)* H07))))))
                0242      &       + x2*(H20 + (Tp(i,j)*(H21 +  Tp(i,j)*(H22  + Tp(i,j)*(H23
                0243      &       + Tp(i,j)*(H24 + Tp(i,j)*(H25 + Tp(i,j)*H26)))))
                0244      &       +  x*(H30 + (Tp(i,j)*(H31  + Tp(i,j)*(H32
                0245      &       + Tp(i,j)*(H33 + Tp(i,j)* H34)))
                0246      &       + x*(H40 + (Tp(i,j)*(H41 + Tp(i,j)*(H42 + Tp(i,j)*(H43
                0247      &       + Tp(i,j)*(H44 + Tp(i,j)*H45))))
                0248      &       + x*(H50 + x*(H60 + x* H70)) )) )) )) )
                0249        ENDDO
                0250       ENDDO
                0251
                0252       RETURN
                0253       END
                0254
                0255 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0256
                0257 CBOP
                0258 C !ROUTINE: CONVERT_CT2PT
                0259
                0260 C !INTERFACE:
                0261       SUBROUTINE CONVERT_CT2PT(
                0262      I         Tc, Sa,
                0263      O         Tp,
                0264      I         myTime, myIter, myThid )
                0265
                0266 C     !DESCRIPTION:
                0267 C     Convert input conservative (degC) and absolute salinity
                0268 C     (g kg-1) to returned potential temperature (degC) by inverting
                0269 C     the polynomial expressions from TEOS-10.
                0270
                0271 C     !USES:
                0272       IMPLICIT NONE
                0273 C     == Global variables ===
                0274 #include "SIZE.h"
                0275 #include "EOS.h"
                0276
                0277 C     !INPUT/OUTPUT PARAMETERS:
                0278 C     Tc        :: Conservative temperature (degC)
                0279 C     Sa        :: Absolute salinity        (g kg-1)
                0280 C     Tp        :: potential temperature (degC)
                0281 C     bi,bj     :: Current tile indices
                0282 C     myTime    :: Current time in simulation
                0283 C     myIter    :: Current iteration number
                0284 C     myThid    :: my Thread Id number
                0285       _RL Tc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0286       _RL Sa(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0287       _RL Tp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0288       _RL myTime
                0289       INTEGER myIter, myThid
                0290 CEOP
                0291
                0292 C     !LOCAL VARIABLES:
                0293       INTEGER i,j
                0294 C     The numerator of a simple expression for potential temperature (degC)
                0295       _RL Tp_num
                0296 C     The inverse of the denominator of a simple expression for
                0297 C     potential temperature (nondim)
                0298       _RL I_Tp_den
                0299 C     The difference between an estimate of conservative temperature and
                0300 C     its target (degC)
                0301       _RL Tc_diff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0302 C     A previous estimate of the potential tempearture (degC)
                0303       _RL Tp_old(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0304 C     The partial derivative of potential temperature with conservative
                0305 C     temperature and inverse (nondim)
                0306       _RL dTp_dTc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0307       _RL dTc_dTp
                0308 C     Absolute salinity normalized by a plausible salinity range and its
                0309 C     square root (nondim)
                0310       _RL x2, x
                0311 C     intermediate temperature (degC)
                0312       _RL Tmp
                0313 C
                0314       _RL     zeRL
                0315       PARAMETER ( zeRL = 0.0 _d 0 )
                0316
                0317       DO j=1-OLy,sNy+OLy
                0318        DO i=1-OLx,sNx+OLx
                0319         Tp_num = TPN00 + (Sa(i,j)*(TPN10 + TPN20*Sa(i,j))
                0320      &       + Tc(i,j)*(TPN01 + (TPN11*Sa(i,j) + TPN02*Tc(i,j))))
                0321         I_Tp_den = 1.0 _d 0 / (1.0 _d 0
                0322      &       + (TPD10*Sa(i,j) + Tc(i,j)*(TPD01 + TPD02*Tc(i,j))))
                0323         Tp(i,j) = Tp_num*I_Tp_den
                0324         dTp_dTc(i,j) = ((TPN01 + (TPN11*Sa(i,j) + 2.*TPN02*Tc(i,j)))
                0325      &       - (TPD01 + 2.*TPD02*Tc(i,j))*Tp(i,j))*I_Tp_den
                0326        ENDDO
                0327       ENDDO
                0328
                0329 C--   Start the 1.5 iterations through the modified Newton-Raphson
                0330 C     iterative method, which is also known as the Newton-McDougall
                0331 C     method.  In this case 1.5 iterations converge to 64-bit machine
                0332 C     precision for oceanographically relevant temperatures and
                0333 C     salinities.
                0334       CALL CONVERT_PT2CT (
                0335      I     Tp, Sa,
                0336      O     Tc_diff,
                0337      I     myTime, myIter, myThid )
                0338       DO j=1-OLy,sNy+OLy
                0339        DO i=1-OLx,sNx+OLx
                0340         Tc_diff(i,j) = Tc_diff(i,j) - Tc(i,j)
                0341         Tp_old(i,j) = Tp(i,j)
                0342 #ifdef ALLOW_AUTODIFF
                0343        ENDDO
                0344       ENDDO
                0345 # ifdef ALLOW_AUTODIFF_TAMC
                0346 CADJ INIT loctape_cp2ct = COMMON, 1
                0347 CADJ STORE Tc_diff = loctape_cp2ct
                0348 CADJ STORE dTp_dTc = loctape_cp2ct
                0349 # endif
                0350       DO j=1-OLy,sNy+OLy
                0351        DO i=1-OLx,sNx+OLx
                0352 #endif
                0353         Tp(i,j) = Tp_old(i,j) - Tc_diff(i,j)*dTp_dTc(i,j)
                0354 C--   Estimate the potential temperature and its derivative from an C
                0355 C     approximate rational function fit.
                0356         x2 = MAX(I_S0 * Sa(i,j), zeRL)
                0357 #ifdef ALLOW_AUTODIFF
                0358         x = 0. _d 0
                0359         IF ( x2 .GT. 0. _d 0 ) x = SQRT(x2)
                0360 #else
                0361         x = SQRT(x2)
                0362 #endif
                0363         Tmp = 0.5 _d 0 *(Tp(i,j) + Tp_old(i,j))
                0364         dTc_dTp = (     H01 + Tmp*(2.*H02 + Tmp*(3.*H03 + Tmp*(4.*H04
                0365      &       + Tmp*(5.*H05 + Tmp*(6.*H06 + Tmp*(7.*H07)))))) )
                0366      &       + x2*(     (H21 + Tmp*(2.*H22 + Tmp*(3.*H23 + Tmp*(4.*H24
                0367      &       + Tmp*(5.*H25 + Tmp*(6.*H26))))))
                0368      &       +  x*(  (H31 + Tmp*(2.*H32 + Tmp*(3.*H33 + Tmp*(4.*H34))))
                0369      &       + x*(H41 + Tmp*(2.*H42 + Tmp*(3.*H43
                0370      &       + Tmp*(4.*H44 + Tmp*(5.*H45))))) ) )
                0371         dTp_dTc(i,j) = 1. _d 0 / dTc_dTp
                0372 #ifdef ALLOW_AUTODIFF
                0373        ENDDO
                0374       ENDDO
                0375 # ifdef ALLOW_AUTODIFF_TAMC
                0376 CADJ STORE dTp_dTc = loctape_cp2ct
                0377 # endif
                0378       DO j=1-OLy,sNy+OLy
                0379        DO i=1-OLx,sNx+OLx
                0380 #endif
                0381         Tp(i,j) = Tp_old(i,j) - Tc_diff(i,j)*dTp_dTc(i,j)
                0382        ENDDO
                0383       ENDDO
                0384       CALL CONVERT_PT2CT (
                0385      I     Tp, Sa,
                0386      O     Tc_diff,
                0387      I     myTime, myIter, myThid )
                0388       DO j=1-OLy,sNy+OLy
                0389        DO i=1-OLx,sNx+OLx
                0390         Tc_diff(i,j) = Tc_diff(i,j) - Tc(i,j)
                0391         Tp(i,j) = Tp(i,j) - Tc_diff(i,j)*dTp_dTc(i,j)
                0392        ENDDO
                0393       ENDDO
                0394
                0395       RETURN
                0396       END