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
