source: issm/oecreview/Archive/15392-16133/ISSM-15831-15832.diff

Index: ../trunk-jpl/src/m/mech/steadystateiceshelftemp.py
===================================================================
--- ../trunk-jpl/src/m/mech/steadystateiceshelftemp.py  (revision 0)
+++ ../trunk-jpl/src/m/mech/steadystateiceshelftemp.py  (revision 15832)
@@ -0,0 +1,65 @@
+import numpy as npy
+
+def steadystateiceshelftemp(md,surfacetemp,basaltemp):
+       """
+       Compute the depth-averaged steady-state temperature of an ice shelf 
+       This routine computes the depth-averaged temperature accounting for vertical advection 
+       and diffusion of heat into the base of the ice shelf as a function of surface and basal 
+       temperature and the basal melting rate.  Horizontal advection is ignored.
+   The solution is a depth-averaged version of Equation 25 in Holland and Jenkins (1999).
+
+       In addition to supplying md, the surface and basal temperatures of the ice shelf must be supplied in degrees Kelvin.
+
+       The model md must also contain the fields: 
+       md.geometry.thickness
+       md.basalforcings.melting_rate (positive for melting, negative for freezing)
+
+   Usage:
+      temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp)
+       """
+
+       if len(md.geometry.thickness)!=md.mesh.numberofvertices:
+               raise ValueError('steadystateiceshelftemp error message: thickness should have a length of ' + md.mesh.numberofvertices)
+       
+       #surface and basal temperatures in degrees C
+       if len(surfacetemp)!=md.mesh.numberofvertices:
+               raise ValueError('steadystateiceshelftemp error message: surfacetemp should have a length of ' + md.mesh.numberofvertices)
+       
+       if len(basaltemp)!=md.mesh.numberofvertices:
+               raise ValueError('steadystateiceshelftemp error message: basaltemp should have a length of ' +md.mesh.numberofvertices)
+       
+       # Convert temps to Celsius for Holland and Jenkins (1999) equation
+       Ts=-273.15+surfacetemp
+       Tb=-273.15+basaltemp
+       
+       Hi=md.geometry.thickness
+       ki=1.14e-6*md.constants.yts # ice shelf thermal diffusivity from Holland and Jenkins (1999) converted to m^2/yr 
+       
+       #vertical velocity of ice shelf, calculated from melting rate 
+       wi=md.materials.rho_water/md.materials.rho_ice*md.basalforcings.melting_rate 
+       
+       #temperature profile is linear if melting rate is zero, depth-averaged temp is simple average in this case
+       temperature=(Ts+Tb)/2  # where wi~=0
+       
+       pos=npy.nonzero(abs(wi)>=1e-4) # to avoid division by zero
+
+       npy.seterr(over='raise',divide='raise') # raise errors if floating point exceptions are encountered in following calculation
+       #calculate depth-averaged temperature (in Celsius)
+       try:
+               temperature[pos]=-( (Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos] - (Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])*npy.exp(Hi[pos]*wi[pos]/ki) )/( Hi[pos]*(npy.exp(Hi[pos]*wi[pos]/ki)-1))
+       except FloatingPointError:
+               print 'steadystateiceshelf warning: overflow encountered in multipy/divide/exp, trying another formulation.' 
+               temperature[pos]=-( ((Tb[pos]-Ts[pos])*ki/wi[pos] + Hi[pos]*Tb[pos])/npy.exp(Hi[pos]*wi[pos]/ki) - Hi[pos]*Ts[pos] + (Tb[pos]-Ts[pos])*ki/wi[pos])/( Hi[pos]*(1-npy.exp(-Hi[pos]*wi[pos]/ki)))
+       
+       #temperature should not be less than surface temp
+       pos=npy.nonzero(temperature<Ts)
+       temperature[pos]=Ts[pos]
+       
+       # NaN where melt rates are too high (infinity/infinity in exponential)
+       pos=npy.nonzero(npy.isnan(temperature))
+       temperature[pos]=Ts[pos]
+       
+       #convert to Kelvin
+       temperature=temperature+273.15
+
+       return temperature
Index: ../trunk-jpl/src/m/mech/steadystateiceshelftemp.m
===================================================================
--- ../trunk-jpl/src/m/mech/steadystateiceshelftemp.m   (revision 15831)
+++ ../trunk-jpl/src/m/mech/steadystateiceshelftemp.m   (revision 15832)
@@ -45,8 +45,8 @@
 pos=find(abs(wi)>=1e-4); % to avoid division by zero
 
 %calculate depth-averaged temperature (in Celsius)
-temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1));
-%temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki)));
+%temperature(pos)=-( (Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos) - (Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos)).*exp(Hi(pos).*wi(pos)/ki) )./( Hi(pos).*(exp(Hi(pos).*wi(pos)/ki)-1));
+temperature(pos)=-( ((Tb(pos)-Ts(pos))*ki./wi(pos) + Hi(pos).*Tb(pos))./exp(Hi(pos).*wi(pos)/ki) - Hi(pos).*Ts(pos) + (Tb(pos)-Ts(pos))*ki./wi(pos))./( Hi(pos).*(1-exp(-Hi(pos).*wi(pos)/ki)));
 
 %temperature should not be less than surface temp
 pos=find(temperature<Ts);