Index: sm/trunk-jpl/src/m/meca/initialdamage.m
===================================================================
--- /issm/trunk-jpl/src/m/meca/initialdamage.m	(revision 13123)
+++ 	(revision )
@@ -1,41 +1,0 @@
-function damage=initialdamage(md)
-%INITIALDAMAGE - compute initial damage guess for ICE SHELF inverse control method 
-%
-%	 This routine computes the initial damage guess as a function of
-%	 material properties, ice thickness, strain rate, and ice 
-%	 rigidity.  The model must contain computed strain rates.
-%
-%   Usage:
-%      md=initialdamage(md)
-%
-%   Example:
-%      md=initialdamage(md);
-
-%some checks
-if isempty(fieldnames(md.results)),
-	error(['md.results.strainrate is not present.  Calculate using md=mechanicalproperties(md,vx,vy)'])
-end
-if ~(md.mesh.dimension==2)
-	error('only 2d model supported currently');
-end
-if any(md.flowequation.element_equation~=2),
-	disp('Warning: the model has some non macayeal elements. These will be treated like MacAyeal''s elements');
-end
-
-%average results onto vertices
-eps1=averaging(md,md.results.strainrate.principalvalue1,0)/md.constants.yts; % s^-1 for strain rates
-epsxx=averaging(md,md.results.strainrate.xx,0)/md.constants.yts;
-epsyy=averaging(md,md.results.strainrate.yy,0)/md.constants.yts;
-epseff=averaging(md,md.results.strainrate.effectivevalue,0)/md.constants.yts;
-n=averaging(md,md.materials.rheology_n,0);
-
-%lump material constants together
-const=md.materials.rho_ice*md.constants.g*(1-md.materials.rho_ice/md.materials.rho_water)/2;
-
-damage=1-0.5*const.*md.geometry.thickness./md.materials.rheology_B./eps1.^(1./n);
-%damage=1-const.*md.geometry.thickness.*epseff.^((n-1)./n)./md.materials.rheology_B./(eps1+epsxx+epsyy);
-
-pos=find(damage>1);
-damage(pos)=1;
-pos=find(damage<0);
-damage(pos)=0;
Index: sm/trunk-jpl/src/m/meca/steadystateiceshelftemp.m
===================================================================
--- /issm/trunk-jpl/src/m/meca/steadystateiceshelftemp.m	(revision 13123)
+++ 	(revision )
@@ -1,51 +1,0 @@
-function temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp)
-%STEADYSTATEICESHELFTEMP - compute 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
-
-%   Usage:
-%      temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp)
-
-if (length(md.geometry.thickness)~=md.mesh.numberofvertices)
-	error(['steadystateiceshelftemp error message: thickness should have a length of ' num2str(md.mesh.numberofvertices)])
-end
-
-%surface and basal temperatures in degrees C
-if (length(surfacetemp)~=md.mesh.numberofvertices)
-	error(['steadystateiceshelftemp error message: surfacetemp should have a length of ' num2str(md.mesh.numberofvertices)])
-end
-
-if (length(basaltemp)~=md.mesh.numberofvertices)
-	error(['steadystateiceshelftemp error message: basaltemp should have a length of ' num2str(md.mesh.numberofvertices)])
-end
-
-% Convert temps to Celsius for Holland and Jenkins equation
-Ts=-273+surfacetemp;
-Tb=-273+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=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));
-
-%convert to Kelvin
-temperature=temperature+273;