C $Header: /u/gcmpack/MITgcm/model/src/gsw_teos10.F,v 1.2 2015/02/27 08:18:12 mlosch Exp $
C $Name: $
#include "CPP_OPTIONS.h"
C-- File gsw_teos10.F: routines that compute quantities related to seawater.
C-- Contents
C-- TEOS-10 routines (Gibbs seawater, GSW)
C-- o GSW_PT_FROM_CT: function to compute potential temperature
C-- from conservative temperature and absolute salinity
C-- o GSW_CT_FROM_PT: function to compute conservative temperature with
C-- from potential temperature and absolute salinity
C-- o GSW_GIBBS_PT0_PT0: function to compute specific Gibbs free energy
C-- from potential temperature and absolute salinity
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: GSW_PT_FROM_CT
C !INTERFACE:
_RL FUNCTION GSW_PT_FROM_CT(SA,CT)
C !DESCRIPTION: \bv
C *=============================================================*
C | S/R GSW_PT_FROM_CT
C | o compute potential temperature at reference level 0 dbar
C | from conservative temperature (CT) and absolute
C | salinity (SA)
C | o this is a more or less shameless copy of the teos-10 code
C | available at http://www.teos-10.org
C *=============================================================*
C \ev
C !USES:
IMPLICIT NONE
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C SA :: Absolute Salinity (g/kg)
C CT :: Conservative Temperature (deg C)
_RL SA,CT
C !FUNCTIONS:
C == Functions ==
CML _RL gsw_gibbs
_RL gsw_gibbs_pt0_pt0
_RL gsw_ct_from_pt
CML EXTERNAL gsw_gibbs, gsw_gibbs_pt0_pt0, gsw_ct_from_pt
EXTERNAL , gsw_ct_from_pt
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER n0, n2
_RL s1, p0, cp0
_RL a0, a1, a2, a3, a4, a5, b0, b1, b2, b3
_RL a5ct, b3ct, ct_factor, pt_num, pt_den, ct_diff
_RL pt, pt_old, ptm, dct_dpt
CEOP
cp0 = 3991.86795711963 _d 0
n0 = 0
n2 = 2
s1 = SA * 35. _d 0 / 35.16504 _d 0
p0 = 0. _d 0
a0 = -1.446013646344788 _d -2
a1 = -3.305308995852924 _d -3
a2 = 1.062415929128982 _d -4
a3 = 9.477566673794488 _d -1
a4 = 2.166591947736613 _d -3
a5 = 3.828842955039902 _d -3
b0 = 1.000000000000000 _d +0
b1 = 6.506097115635800 _d -4
b2 = 3.830289486850898 _d -3
b3 = 1.247811760368034 _d -6
a5ct = a5*CT
b3ct = b3*CT
ct_factor = (a3 + a4*s1 + a5ct)
pt_num = a0 + s1*(a1 + a2*s1) + CT*ct_factor
pt_den = b0 + b1*s1 + CT*(b2 + b3ct)
pt = (pt_num)/(pt_den)
dct_dpt = (pt_den)/(ct_factor + a5ct - (b2 + b3ct + b3ct)*pt)
C start the 1.5 iterations through the modified Newton-Rapshon
C iterative method.
ct_diff = gsw_ct_from_pt(sa,pt) - CT
pt_old = pt
pt = pt_old - (ct_diff)/dct_dpt
ptm = 0.5 _d 0*(pt + pt_old)
dct_dpt = -(ptm + 273.15 _d 0)*gsw_gibbs_pt0_pt0(sa,ptm)/cp0
pt = pt_old - (ct_diff)/dct_dpt
ct_diff = gsw_ct_from_pt(sa,pt) - CT
pt_old = pt
GSW_PT_FROM_CT = pt_old - (ct_diff)/dct_dpt
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: GSW_CT_FROM_PT
C !INTERFACE:
_RL FUNCTION GSW_CT_FROM_PT(SA,PT)
C !DESCRIPTION: \bv
C *=============================================================*
C | S/R GSW_CT_FROM_PT
C | o compute conservative temperature from potential
C | temperature (PT) at reference level 0 dbar and absolute
C | salinity (SA)
C | o this is a more or less shameless copy of the teos-10 code
C | available at http://www.teos-10.org
C *=============================================================*
C \ev
C !USES:
IMPLICIT NONE
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C SA :: Absolute Salinity (g/kg)
C PT :: Potential Temperature (deg C)
_RL sa, pt
C !FUNCTIONS:
C == Functions ==
C !LOCAL VARIABLES:
C == Local variables ==
_RL pot_enthalpy, sfac
_RL x2, x, y, cp0
CEOP
sfac = 0.0248826675584615 _d 0
x2 = sfac*sa
x = 0. _d 0
if (x2.gt.0. _d 0) x = sqrt(x2)
C normalize for F03 and F08
y = pt*0.025 _d 0
pot_enthalpy = 61.01362420681071 _d 0 +
& y*(168776.46138048015 _d 0 +
& y*(-2735.2785605119625 _d 0 + y*(2574.2164453821433 _d 0 +
& y*(-1536.6644434977543 _d 0 + y*(545.7340497931629 _d 0 +
& (-50.91091728474331 _d 0 - 18.30489878927802 _d 0*y)*y))))) +
& x2*(268.5520265845071 _d 0 + y*(-12019.028203559312 _d 0 +
& y *(3734.858026725145 _d 0 + y*(-2046.7671145057618 _d 0 +
& y*(465.28655623826234 _d 0 + (-0.6370820302376359 _d 0 -
& 10.650848542359153 _d 0*y)*y)))) +
& x*(937.2099110620707 _d 0 + y*(588.1802812170108 _d 0 +
& y*(248.39476522971285 _d 0 + (-3.871557904936333 _d 0 -
& 2.6268019854268356 _d 0*y)*y)) +
& x*(-1687.914374187449 _d 0 + x*(246.9598888781377 _d 0 +
& x*(123.59576582457964 _d 0 - 48.5891069025409 _d 0*x)) +
& y*( 936.3206544460336 _d 0 +
& y*(-942.7827304544439 _d 0 + y*(369.4389437509002 _d 0 +
& (-33.83664947895248 _d 0 - 9.987880382780322 _d 0*y)*y))))))
cp0 = 3991.86795711963 _d 0
gsw_ct_from_pt = pot_enthalpy/cp0
RETURN
END
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: GSW_GIBBS_PT0_PT0
C !INTERFACE:
_RL FUNCTION GSW_GIBBS_PT0_PT0(SA,PT0)
C !DESCRIPTION: \bv
C *=============================================================*
C | S/R GSW_GIBBS_PT0_PT0
C | o helper routine that computes the specific Gibbs free
C | energy from potential temperature (PT) and absolute
C | salinity (SA) for pressure 0 dbar
C | o this is a more or less shameless copy of the teos-10 code
C | available at http://www.teos-10.org
C *=============================================================*
C \ev
C !USES:
IMPLICIT NONE
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C SA :: Absolute Salinity (g/kg)
C PT :: Potential Temperature at p = 0 dbar (deg C)
_RL sa, pt0
C !FUNCTIONS:
C == Functions ==
C !LOCAL VARIABLES:
C == Local variables ==
_RL sfac, x2, x, y, g03, g08
CEOP
sfac = 0.0248826675584615
x2 = sfac*sa
x = 0. _d 0
if (x2.gt.0. _d 0) x = sqrt(x2)
y = pt0*0.025 _d 0
g03 = -24715.571866078 +
& y*(4420.4472249096725 +
& y*(-1778.231237203896 +
& y*(1160.5182516851419 +
& y*(-569.531539542516 + y*128.13429152494615))))
g08 = x2*(1760.062705994408 + x*(-86.1329351956084 +
& x*( -137.1145018408982 + y*(296.20061691375236 +
& y* (-205.67709290374563 + 49.9394019139016*y))) +
& y*( -60.136422517125 + y*10.50720794170734)) +
& y*(-1351.605895580406 + y*(1097.1125373015109 +
& y*( -433.20648175062206 + 63.905091254154904*y))))
gsw_gibbs_pt0_pt0 = (g03 + g08)*0.000625
RETURN
END