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