C $Header: /home/ubuntu/mnt/e9_copy/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_gibbs_pt0_pt0, 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