pkg/openad/seawater.F

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