| 1 |
C $Header: /u/gcmpack/MITgcm/model/src/gsw_teos10.F,v 1.1 2015/02/26 17:19:52 mlosch Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "CPP_OPTIONS.h" |
| 5 |
|
| 6 |
C-- File gsw_teos10.F: routines that compute quantities related to seawater. |
| 7 |
C-- Contents |
| 8 |
C-- TEOS-10 routines (Gibbs seawater, GSW) |
| 9 |
C-- o GSW_PT_FROM_CT: function to compute potential temperature |
| 10 |
C-- from conservative temperature and absolute salinity |
| 11 |
C-- o GSW_CT_FROM_PT: function to compute conservative temperature with |
| 12 |
C-- from potential temperature and absolute salinity |
| 13 |
C-- o GSW_GIBBS_PT0_PT0: function to compute specific Gibbs free energy |
| 14 |
C-- from potential temperature and absolute salinity |
| 15 |
|
| 16 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 17 |
|
| 18 |
CBOP |
| 19 |
C !ROUTINE: GSW_PT_FROM_CT |
| 20 |
C !INTERFACE: |
| 21 |
_RL FUNCTION GSW_PT_FROM_CT(SA,CT) |
| 22 |
C !DESCRIPTION: \bv |
| 23 |
C *=============================================================* |
| 24 |
C | S/R GSW_PT_FROM_CT |
| 25 |
C | o compute potential temperature at reference level 0 dbar |
| 26 |
C | from conservative temperature (CT) and absolute |
| 27 |
C | salinity (SA) |
| 28 |
C | o this is a more or less shameless copy of the teos-10 code |
| 29 |
C | available at http://www.teos-10.org |
| 30 |
C *=============================================================* |
| 31 |
C \ev |
| 32 |
|
| 33 |
C !USES: |
| 34 |
IMPLICIT NONE |
| 35 |
|
| 36 |
C !INPUT/OUTPUT PARAMETERS: |
| 37 |
C == Routine arguments == |
| 38 |
C SA :: Absolute Salinity (g/kg) |
| 39 |
C CT :: Conservative Temperature (deg C) |
| 40 |
_RL SA,CT |
| 41 |
|
| 42 |
C !FUNCTIONS: |
| 43 |
C == Functions == |
| 44 |
CML _RL gsw_gibbs |
| 45 |
_RL gsw_gibbs_pt0_pt0 |
| 46 |
_RL gsw_ct_from_pt |
| 47 |
CML EXTERNAL gsw_gibbs, gsw_gibbs_pt0_pt0, gsw_ct_from_pt |
| 48 |
EXTERNAL gsw_gibbs_pt0_pt0, gsw_ct_from_pt |
| 49 |
|
| 50 |
C !LOCAL VARIABLES: |
| 51 |
C == Local variables == |
| 52 |
INTEGER n0, n2 |
| 53 |
_RL s1, p0, cp0 |
| 54 |
_RL a0, a1, a2, a3, a4, a5, b0, b1, b2, b3 |
| 55 |
_RL a5ct, b3ct, ct_factor, pt_num, pt_den, ct_diff |
| 56 |
_RL pt, pt_old, ptm, dct_dpt |
| 57 |
CEOP |
| 58 |
|
| 59 |
cp0 = 3991.86795711963 _d 0 |
| 60 |
|
| 61 |
n0 = 0 |
| 62 |
n2 = 2 |
| 63 |
|
| 64 |
s1 = SA * 35. _d 0 / 35.16504 _d 0 |
| 65 |
p0 = 0. _d 0 |
| 66 |
|
| 67 |
a0 = -1.446013646344788 _d -2 |
| 68 |
a1 = -3.305308995852924 _d -3 |
| 69 |
a2 = 1.062415929128982 _d -4 |
| 70 |
a3 = 9.477566673794488 _d -1 |
| 71 |
a4 = 2.166591947736613 _d -3 |
| 72 |
a5 = 3.828842955039902 _d -3 |
| 73 |
|
| 74 |
b0 = 1.000000000000000 _d +0 |
| 75 |
b1 = 6.506097115635800 _d -4 |
| 76 |
b2 = 3.830289486850898 _d -3 |
| 77 |
b3 = 1.247811760368034 _d -6 |
| 78 |
|
| 79 |
a5ct = a5*CT |
| 80 |
b3ct = b3*CT |
| 81 |
|
| 82 |
ct_factor = (a3 + a4*s1 + a5ct) |
| 83 |
pt_num = a0 + s1*(a1 + a2*s1) + CT*ct_factor |
| 84 |
pt_den = b0 + b1*s1 + CT*(b2 + b3ct) |
| 85 |
pt = (pt_num)/(pt_den) |
| 86 |
|
| 87 |
dct_dpt = (pt_den)/(ct_factor + a5ct - (b2 + b3ct + b3ct)*pt) |
| 88 |
|
| 89 |
C start the 1.5 iterations through the modified Newton-Rapshon |
| 90 |
C iterative method. |
| 91 |
|
| 92 |
ct_diff = gsw_ct_from_pt(sa,pt) - CT |
| 93 |
pt_old = pt |
| 94 |
pt = pt_old - (ct_diff)/dct_dpt |
| 95 |
ptm = 0.5 _d 0*(pt + pt_old) |
| 96 |
|
| 97 |
dct_dpt = -(ptm + 273.15 _d 0)*gsw_gibbs_pt0_pt0(sa,ptm)/cp0 |
| 98 |
|
| 99 |
pt = pt_old - (ct_diff)/dct_dpt |
| 100 |
ct_diff = gsw_ct_from_pt(sa,pt) - CT |
| 101 |
pt_old = pt |
| 102 |
GSW_PT_FROM_CT = pt_old - (ct_diff)/dct_dpt |
| 103 |
|
| 104 |
RETURN |
| 105 |
END |
| 106 |
|
| 107 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 108 |
|
| 109 |
CBOP |
| 110 |
C !ROUTINE: GSW_CT_FROM_PT |
| 111 |
C !INTERFACE: |
| 112 |
_RL FUNCTION GSW_CT_FROM_PT(SA,PT) |
| 113 |
C !DESCRIPTION: \bv |
| 114 |
C *=============================================================* |
| 115 |
C | S/R GSW_CT_FROM_PT |
| 116 |
C | o compute conservative temperature from potential |
| 117 |
C | temperature (PT) at reference level 0 dbar and absolute |
| 118 |
C | salinity (SA) |
| 119 |
C | o this is a more or less shameless copy of the teos-10 code |
| 120 |
C | available at http://www.teos-10.org |
| 121 |
C *=============================================================* |
| 122 |
C \ev |
| 123 |
|
| 124 |
C !USES: |
| 125 |
IMPLICIT NONE |
| 126 |
|
| 127 |
C !INPUT/OUTPUT PARAMETERS: |
| 128 |
C == Routine arguments == |
| 129 |
C SA :: Absolute Salinity (g/kg) |
| 130 |
C PT :: Potential Temperature (deg C) |
| 131 |
_RL sa, pt |
| 132 |
|
| 133 |
C !FUNCTIONS: |
| 134 |
C == Functions == |
| 135 |
|
| 136 |
C !LOCAL VARIABLES: |
| 137 |
C == Local variables == |
| 138 |
_RL pot_enthalpy, sfac |
| 139 |
_RL x2, x, y, cp0 |
| 140 |
CEOP |
| 141 |
|
| 142 |
sfac = 0.0248826675584615 _d 0 |
| 143 |
|
| 144 |
x2 = sfac*sa |
| 145 |
x = 0. _d 0 |
| 146 |
if (x2.gt.0. _d 0) x = sqrt(x2) |
| 147 |
C normalize for F03 and F08 |
| 148 |
y = pt*0.025 _d 0 |
| 149 |
|
| 150 |
pot_enthalpy = 61.01362420681071 _d 0 + |
| 151 |
& y*(168776.46138048015 _d 0 + |
| 152 |
& y*(-2735.2785605119625 _d 0 + y*(2574.2164453821433 _d 0 + |
| 153 |
& y*(-1536.6644434977543 _d 0 + y*(545.7340497931629 _d 0 + |
| 154 |
& (-50.91091728474331 _d 0 - 18.30489878927802 _d 0*y)*y))))) + |
| 155 |
& x2*(268.5520265845071 _d 0 + y*(-12019.028203559312 _d 0 + |
| 156 |
& y *(3734.858026725145 _d 0 + y*(-2046.7671145057618 _d 0 + |
| 157 |
& y*(465.28655623826234 _d 0 + (-0.6370820302376359 _d 0 - |
| 158 |
& 10.650848542359153 _d 0*y)*y)))) + |
| 159 |
& x*(937.2099110620707 _d 0 + y*(588.1802812170108 _d 0 + |
| 160 |
& y*(248.39476522971285 _d 0 + (-3.871557904936333 _d 0 - |
| 161 |
& 2.6268019854268356 _d 0*y)*y)) + |
| 162 |
& x*(-1687.914374187449 _d 0 + x*(246.9598888781377 _d 0 + |
| 163 |
& x*(123.59576582457964 _d 0 - 48.5891069025409 _d 0*x)) + |
| 164 |
& y*( 936.3206544460336 _d 0 + |
| 165 |
& y*(-942.7827304544439 _d 0 + y*(369.4389437509002 _d 0 + |
| 166 |
& (-33.83664947895248 _d 0 - 9.987880382780322 _d 0*y)*y)))))) |
| 167 |
|
| 168 |
cp0 = 3991.86795711963 _d 0 |
| 169 |
|
| 170 |
gsw_ct_from_pt = pot_enthalpy/cp0 |
| 171 |
|
| 172 |
RETURN |
| 173 |
END |
| 174 |
|
| 175 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
| 176 |
|
| 177 |
CBOP |
| 178 |
C !ROUTINE: GSW_GIBBS_PT0_PT0 |
| 179 |
C !INTERFACE: |
| 180 |
_RL FUNCTION GSW_GIBBS_PT0_PT0(SA,PT0) |
| 181 |
C !DESCRIPTION: \bv |
| 182 |
C *=============================================================* |
| 183 |
C | S/R GSW_GIBBS_PT0_PT0 |
| 184 |
C | o helper routine that computes the specific Gibbs free |
| 185 |
C | energy from potential temperature (PT) and absolute |
| 186 |
C | salinity (SA) for pressure 0 dbar |
| 187 |
C | o this is a more or less shameless copy of the teos-10 code |
| 188 |
C | available at http://www.teos-10.org |
| 189 |
C *=============================================================* |
| 190 |
C \ev |
| 191 |
|
| 192 |
C !USES: |
| 193 |
IMPLICIT NONE |
| 194 |
|
| 195 |
C !INPUT/OUTPUT PARAMETERS: |
| 196 |
C == Routine arguments == |
| 197 |
C SA :: Absolute Salinity (g/kg) |
| 198 |
C PT :: Potential Temperature at p = 0 dbar (deg C) |
| 199 |
_RL sa, pt0 |
| 200 |
|
| 201 |
C !FUNCTIONS: |
| 202 |
C == Functions == |
| 203 |
|
| 204 |
C !LOCAL VARIABLES: |
| 205 |
C == Local variables == |
| 206 |
_RL sfac, x2, x, y, g03, g08 |
| 207 |
CEOP |
| 208 |
|
| 209 |
sfac = 0.0248826675584615 |
| 210 |
|
| 211 |
x2 = sfac*sa |
| 212 |
x = 0. _d 0 |
| 213 |
if (x2.gt.0. _d 0) x = sqrt(x2) |
| 214 |
y = pt0*0.025 _d 0 |
| 215 |
|
| 216 |
g03 = -24715.571866078 + |
| 217 |
& y*(4420.4472249096725 + |
| 218 |
& y*(-1778.231237203896 + |
| 219 |
& y*(1160.5182516851419 + |
| 220 |
& y*(-569.531539542516 + y*128.13429152494615)))) |
| 221 |
|
| 222 |
g08 = x2*(1760.062705994408 + x*(-86.1329351956084 + |
| 223 |
& x*( -137.1145018408982 + y*(296.20061691375236 + |
| 224 |
& y* (-205.67709290374563 + 49.9394019139016*y))) + |
| 225 |
& y*( -60.136422517125 + y*10.50720794170734)) + |
| 226 |
& y*(-1351.605895580406 + y*(1097.1125373015109 + |
| 227 |
& y*( -433.20648175062206 + 63.905091254154904*y)))) |
| 228 |
|
| 229 |
gsw_gibbs_pt0_pt0 = (g03 + g08)*0.000625 |
| 230 |
|
| 231 |
RETURN |
| 232 |
END |
Parent Directory
| 
Revision Log
| 
Revision Graph