model/src/gsw_teos10.F

e2bb6576b3 Mart*0001 #include "CPP_OPTIONS.h"
                0002
                0003 C--  File gsw_teos10.F: routines that compute quantities related to seawater.
                0004 C--   Contents
                0005 C--   TEOS-10 routines (Gibbs seawater, GSW)
                0006 C--   o GSW_PT_FROM_CT: function to compute potential temperature
                0007 C--              from conservative temperature and absolute salinity
                0008 C--   o GSW_CT_FROM_PT: function to compute conservative temperature with
                0009 C--              from potential temperature and absolute salinity
                0010 C--   o GSW_GIBBS_PT0_PT0: function to compute specific Gibbs free energy
                0011 C--              from potential temperature and absolute salinity
                0012
                0013 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0014
                0015 CBOP
                0016 C     !ROUTINE: GSW_PT_FROM_CT
                0017 C     !INTERFACE:
                0018       _RL FUNCTION GSW_PT_FROM_CT(SA,CT)
                0019 C     !DESCRIPTION: \bv
                0020 C     *=============================================================*
                0021 C     | S/R  GSW_PT_FROM_CT
                0022 C     | o compute potential temperature at reference level 0 dbar
                0023 C     |   from conservative temperature (CT) and absolute
                0024 C     |    salinity (SA)
                0025 C     | o this is a more or less shameless copy of the teos-10 code
                0026 C     |   available at http://www.teos-10.org
                0027 C     *=============================================================*
                0028 C     \ev
                0029
                0030 C     !USES:
                0031       IMPLICIT NONE
                0032
                0033 C     !INPUT/OUTPUT PARAMETERS:
                0034 C     == Routine arguments ==
                0035 C     SA :: Absolute Salinity                               (g/kg)
                0036 C     CT :: Conservative Temperature                        (deg C)
                0037       _RL SA,CT
                0038
                0039 C     !FUNCTIONS:
                0040 C     == Functions ==
9930e5740e Mart*0041 CML      _RL gsw_gibbs
e2bb6576b3 Mart*0042       _RL gsw_gibbs_pt0_pt0
                0043       _RL gsw_ct_from_pt
9930e5740e Mart*0044 CML      EXTERNAL gsw_gibbs, gsw_gibbs_pt0_pt0, gsw_ct_from_pt
                0045       EXTERNAL gsw_gibbs_pt0_pt0, gsw_ct_from_pt
e2bb6576b3 Mart*0046
                0047 C     !LOCAL VARIABLES:
                0048 C     == Local variables ==
                0049       INTEGER n0, n2
                0050       _RL s1, p0, cp0
                0051       _RL a0, a1, a2, a3, a4, a5, b0, b1, b2, b3
                0052       _RL a5ct, b3ct, ct_factor, pt_num, pt_den, ct_diff
                0053       _RL pt, pt_old, ptm, dct_dpt
                0054 CEOP
                0055
                0056       cp0 = 3991.86795711963 _d 0
                0057
                0058       n0 = 0
                0059       n2 = 2
                0060
                0061       s1 = SA * 35. _d 0 / 35.16504 _d 0
                0062       p0 = 0. _d 0
                0063
                0064       a0 = -1.446013646344788 _d -2
                0065       a1 = -3.305308995852924 _d -3
                0066       a2 =  1.062415929128982 _d -4
                0067       a3 =  9.477566673794488 _d -1
                0068       a4 =  2.166591947736613 _d -3
                0069       a5 =  3.828842955039902 _d -3
                0070
                0071       b0 =  1.000000000000000 _d +0
                0072       b1 =  6.506097115635800 _d -4
                0073       b2 =  3.830289486850898 _d -3
                0074       b3 =  1.247811760368034 _d -6
                0075
                0076       a5ct = a5*CT
                0077       b3ct = b3*CT
                0078
                0079       ct_factor = (a3 + a4*s1 + a5ct)
                0080       pt_num    = a0 + s1*(a1 + a2*s1) + CT*ct_factor
                0081       pt_den    = b0 + b1*s1 + CT*(b2 + b3ct)
                0082       pt        = (pt_num)/(pt_den)
                0083
                0084       dct_dpt   = (pt_den)/(ct_factor + a5ct - (b2 + b3ct + b3ct)*pt)
                0085
                0086 C     start the 1.5 iterations through the modified Newton-Rapshon
                0087 C     iterative method.
                0088
                0089       ct_diff = gsw_ct_from_pt(sa,pt) - CT
                0090       pt_old  = pt
                0091       pt      = pt_old - (ct_diff)/dct_dpt
                0092       ptm     = 0.5 _d 0*(pt + pt_old)
                0093
                0094       dct_dpt = -(ptm + 273.15 _d 0)*gsw_gibbs_pt0_pt0(sa,ptm)/cp0
                0095
                0096       pt             = pt_old - (ct_diff)/dct_dpt
                0097       ct_diff        = gsw_ct_from_pt(sa,pt) - CT
                0098       pt_old         = pt
                0099       GSW_PT_FROM_CT = pt_old - (ct_diff)/dct_dpt
                0100
                0101       RETURN
                0102       END
                0103
                0104 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0105
                0106 CBOP
                0107 C     !ROUTINE: GSW_CT_FROM_PT
                0108 C     !INTERFACE:
                0109       _RL FUNCTION GSW_CT_FROM_PT(SA,PT)
                0110 C     !DESCRIPTION: \bv
                0111 C     *=============================================================*
                0112 C     | S/R  GSW_CT_FROM_PT
                0113 C     | o compute conservative temperature from potential
                0114 C     |   temperature (PT)  at reference level 0 dbar and absolute
                0115 C     |   salinity (SA)
                0116 C     | o this is a more or less shameless copy of the teos-10 code
                0117 C     |   available at http://www.teos-10.org
                0118 C     *=============================================================*
                0119 C     \ev
                0120
                0121 C     !USES:
                0122       IMPLICIT NONE
                0123
                0124 C     !INPUT/OUTPUT PARAMETERS:
                0125 C     == Routine arguments ==
                0126 C     SA :: Absolute Salinity                               (g/kg)
                0127 C     PT :: Potential Temperature                          (deg C)
                0128       _RL sa, pt
                0129
                0130 C     !FUNCTIONS:
                0131 C     == Functions ==
                0132
                0133 C     !LOCAL VARIABLES:
                0134 C     == Local variables ==
                0135       _RL pot_enthalpy, sfac
                0136       _RL x2, x, y, cp0
                0137 CEOP
                0138
                0139       sfac = 0.0248826675584615 _d 0
                0140
                0141       x2   = sfac*sa
                0142       x    = 0. _d 0
                0143       if (x2.gt.0. _d 0) x = sqrt(x2)
                0144 C     normalize for F03 and F08
                0145       y    = pt*0.025 _d 0
                0146
                0147       pot_enthalpy =  61.01362420681071 _d 0 +
                0148      &     y*(168776.46138048015 _d 0 +
                0149      &     y*(-2735.2785605119625 _d 0 + y*(2574.2164453821433 _d 0 +
                0150      &     y*(-1536.6644434977543 _d 0 + y*(545.7340497931629 _d 0 +
                0151      &     (-50.91091728474331 _d 0 - 18.30489878927802 _d 0*y)*y))))) +
                0152      &     x2*(268.5520265845071 _d 0 + y*(-12019.028203559312 _d 0 +
                0153      &     y *(3734.858026725145 _d 0 + y*(-2046.7671145057618 _d 0 +
                0154      &     y*(465.28655623826234 _d 0 + (-0.6370820302376359 _d 0 -
                0155      &     10.650848542359153 _d 0*y)*y)))) +
                0156      &     x*(937.2099110620707 _d 0 + y*(588.1802812170108 _d 0 +
                0157      &     y*(248.39476522971285 _d 0 + (-3.871557904936333 _d 0 -
                0158      &     2.6268019854268356 _d 0*y)*y)) +
                0159      &     x*(-1687.914374187449 _d 0 + x*(246.9598888781377 _d 0 +
                0160      &     x*(123.59576582457964 _d 0 - 48.5891069025409 _d 0*x)) +
                0161      &     y*( 936.3206544460336 _d 0 +
                0162      &     y*(-942.7827304544439 _d 0 + y*(369.4389437509002 _d 0 +
                0163      &     (-33.83664947895248 _d 0 - 9.987880382780322 _d 0*y)*y))))))
                0164
                0165       cp0 = 3991.86795711963 _d 0
                0166
                0167       gsw_ct_from_pt = pot_enthalpy/cp0
                0168
                0169       RETURN
                0170       END
                0171
                0172 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0173
                0174 CBOP
                0175 C     !ROUTINE: GSW_GIBBS_PT0_PT0
                0176 C     !INTERFACE:
                0177       _RL FUNCTION GSW_GIBBS_PT0_PT0(SA,PT0)
                0178 C     !DESCRIPTION: \bv
                0179 C     *=============================================================*
                0180 C     | S/R GSW_GIBBS_PT0_PT0
                0181 C     | o helper routine that computes the specific Gibbs free
                0182 C     |   energy from potential temperature (PT) and absolute
                0183 C     |   salinity (SA) for pressure 0 dbar
                0184 C     | o this is a more or less shameless copy of the teos-10 code
                0185 C     |   available at http://www.teos-10.org
                0186 C     *=============================================================*
                0187 C     \ev
                0188
                0189 C     !USES:
                0190       IMPLICIT NONE
                0191
                0192 C     !INPUT/OUTPUT PARAMETERS:
                0193 C     == Routine arguments ==
                0194 C     SA :: Absolute Salinity                               (g/kg)
                0195 C     PT :: Potential Temperature at p = 0 dbar             (deg C)
                0196       _RL sa, pt0
                0197
                0198 C     !FUNCTIONS:
                0199 C     == Functions ==
                0200
                0201 C     !LOCAL VARIABLES:
                0202 C     == Local variables ==
                0203       _RL sfac, x2, x, y, g03, g08
                0204 CEOP
                0205
                0206       sfac = 0.0248826675584615
                0207
                0208       x2   = sfac*sa
                0209       x    = 0. _d 0
                0210       if (x2.gt.0. _d 0) x = sqrt(x2)
                0211       y    = pt0*0.025 _d 0
                0212
                0213       g03 = -24715.571866078 +
                0214      &     y*(4420.4472249096725 +
                0215      &     y*(-1778.231237203896 +
                0216      &     y*(1160.5182516851419 +
                0217      &     y*(-569.531539542516  + y*128.13429152494615))))
                0218
                0219       g08 = x2*(1760.062705994408  + x*(-86.1329351956084 +
                0220      &     x*( -137.1145018408982  + y*(296.20061691375236 +
                0221      &     y* (-205.67709290374563 + 49.9394019139016*y))) +
                0222      &     y*(  -60.136422517125   + y*10.50720794170734)) +
                0223      &     y*(-1351.605895580406   + y*(1097.1125373015109 +
                0224      &     y*( -433.20648175062206 + 63.905091254154904*y))))
                0225
                0226       gsw_gibbs_pt0_pt0 = (g03 + g08)*0.000625
                0227
                0228       RETURN
                0229       END