gsw_gibbs_ice
Gibbs energy of ice and its derivatives
Contents
USAGE:
gibbs_ice = gsw_gibbs_ice(nt,np,t,p)
DESCRIPTION:
Calculates Ice specific Gibbs energy and derivatives up
to order 2. The Gibbs function for sea ice is that of TEOS-10
(IOC et al., 2010), being the sum of IAPWS-08 for the saline part and
IAPWS-09 for the pure water part.
INPUT:
nt = order of t derivative [ integers 0, 1 or 2 ]
np = order of p derivative [ integers 0, 1 or 2 ]
t = in-situ temperature (ITS-90) [ deg C ]
p = sea pressure [ dbar ]
(ie. absolute pressure - 10.1325 dbar)
t and p need to have the same dimensions.
OUTPUT:
gibbs = Specific Gibbs energy or its derivatives.
The Gibbs energy (when nt = np = 0) has units of:
[ J/kg ]
The temperature derivatives are output in units of:
[ (J/kg) (K)^(-nt) ]
The pressure derivatives are output in units of:
[ (J/kg) (Pa)^(-np) ]
The mixed derivatives are output in units of:
[ (J/kg) (K)^(-nt) (Pa)^(-np) ]
Note. The derivatives are taken with respect to pressure in Pa, not
withstanding that the pressure input into this routine is in dbar.
AUTHOR:
Trevor McDougall and Paul Barker [ help@teos-10.org ]
VERSION NUMBER:
3.05 (16th February, 2015)
REFERENCES:
IAPWS, 2009: Revised Release release on the Equation of State 2006 for
H2O Ice Ih. The International Association for the Properties of Water
and Steam. Doorwerth, The Netherlands, September 2009, available from
http://www.iapws.org.
IOC, SCOR and IAPSO, 2010: The international thermodynamic equation of
seawater - 2010: Calculation and use of thermodynamic properties.
Intergovernmental Oceanographic Commission, Manuals and Guides No. 56,
UNESCO (English), 196 pp. Available from http://www.TEOS-10.org
The software is available from http://www.TEOS-10.org
function gibbs_ice = gsw_gibbs_ice(nt,np,t,p)
rec_Pt = 1.634903221903779e-3;
Tt = 273.16;
rec_Tt = 3.660858105139845e-3;
T = t + 273.15;
tau = T.*rec_Tt;
db2Pa = 1e4;
dzi = db2Pa.*p.*rec_Pt;
g00 = -6.32020233335886e5;
g01 = 6.55022213658955e-1;
g02 = -1.89369929326131e-8;
g03 = 3.3974612327105304e-15;
g04 = -5.564648690589909e-22;
s0 = -3.32733756492168e3;
t1 = (3.68017112855051e-2 + 5.10878114959572e-2i);
t2 = (3.37315741065416e-1 + 3.35449415919309e-1i);
r1 = (4.47050716285388e1 + 6.56876847463481e1i);
r20 = (-7.25974574329220e1 - 7.81008427112870e1i);
r21 = (-5.57107698030123e-5 + 4.64578634580806e-5i);
r22 = ( 2.34801409215913e-11 - 2.85651142904972e-11i);
if nt == 0 & np == 0
tau_t1 = tau./t1;
sqtau_t1 = tau_t1.*tau_t1;
tau_t2 = tau./t2;
sqtau_t2 = tau_t2.*tau_t2;
g0 = g00 + dzi.*(g01 + dzi.*(g02 + dzi.*(g03 + g04.*dzi)));
r2 = r20 + dzi.*(r21 + r22.*dzi);
g = r1.*(tau.*log((1 + tau_t1)./(1 - tau_t1)) + t1.*(log(1 - sqtau_t1) - sqtau_t1)) ...
+ r2.*(tau.*log((1 + tau_t2)./(1 - tau_t2)) + t2.*(log(1 - sqtau_t2) - sqtau_t2));
gibbs_ice = g0 - Tt.*(s0.*tau - real(g));
elseif nt == 1 & np == 0
tau_t1 = tau./t1;
tau_t2 = tau./t2;
r2 = r20 + dzi.*(r21 + r22.*dzi);
g = r1.*(log((1 + tau_t1)./(1 - tau_t1)) - 2.*tau_t1) ...
+ r2.*(log((1 + tau_t2)./(1 - tau_t2)) - 2.*tau_t2);
gibbs_ice = -s0 + real(g);
elseif nt == 0 & np == 1
tau_t2 = tau./t2;
sqtau_t2 = tau_t2.*tau_t2;
g0p = rec_Pt.*(g01 + dzi.*(2.*g02 + dzi.*(3.*g03 + 4.*g04.*dzi)));
r2p = rec_Pt.*(r21 + 2.*r22.*dzi);
g = r2p.*(tau.*log((1 + tau_t2)./(1 - tau_t2)) + t2.*(log(1 - sqtau_t2) ...
- sqtau_t2));
gibbs_ice = g0p + Tt.*real(g);
elseif nt == 1 & np == 1
tau_t2 = tau./t2;
r2p = rec_Pt.*(r21 + 2.*r22.*dzi);
g = r2p.*(log((1 + tau_t2)./(1 - tau_t2)) - 2.*tau_t2);
gibbs_ice = real(g);
elseif nt == 2 & np == 0
r2 = r20 + dzi.*(r21 + r22.*dzi);
g = r1.*(1./(t1 - tau) + 1./(t1 + tau) - 2./t1) ...
+ r2.*(1./(t2 - tau) + 1./(t2 + tau) - 2./t2);
gibbs_ice = rec_Tt.*real(g);
elseif nt == 0 & np == 2
sqrec_Pt = rec_Pt.*rec_Pt;
tau_t2 = tau./t2;
sqtau_t2 = tau_t2.*tau_t2;
g0pp = sqrec_Pt.*(2.*g02 + dzi.*(6.*g03 + 12.*g04.*dzi));
r2pp = 2.*r22.*sqrec_Pt;
g = r2pp.*(tau.*log((1 + tau_t2)./(1 - tau_t2)) + t2.*(log(1 - sqtau_t2) ...
- sqtau_t2));
gibbs_ice = g0pp + Tt.*real(g);
end
end