Changeset 12917

Timestamp:

08/06/12 16:36:15 (13 years ago)

Author:

cborstad

Message:

CHG: added surface and basal temp enums, updated function to calculate depth-average ice shelf temp

Location:

issm/branches/trunk-jpl-damage/src

Files:

• c/EnumDefinitions/EnumDefinitions.h (modified) (1 diff)
• c/modules/EnumToStringx/EnumToStringx.cpp (modified) (1 diff)
• c/modules/StringToEnumx/StringToEnumx.cpp (modified) (2 diffs)
• m/classes/initialization.m (modified) (3 diffs)
• m/enum/EnumDefinitions.py (modified) (2 diffs)
• m/enum/MaximumNumberOfEnums.m (modified) (1 diff)
• m/enum/TemperatureBasalEnum.m (added)
• m/enum/TemperatureSurfaceEnum.m (added)
• m/model/plot/applyoptions.m (modified) (1 diff)
• m/model/steadystateiceshelftemp.m (modified) (3 diffs)

issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h

@@ -362 +362 @@
         TemperatureOldEnum,
         TemperaturePicardEnum,
+        TemperatureSurfaceEnum,
+        TemperatureBasalEnum,
         ThicknessAbsMisfitEnum,
         TypeEnum,

issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -353 +353 @@
                 case TemperatureOldEnum : return "TemperatureOld";
                 case TemperaturePicardEnum : return "TemperaturePicard";
+                case TemperatureSurfaceEnum : return "TemperatureSurface";
+                case TemperatureBasalEnum : return "TemperatureBasal";
                 case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit";
                 case TypeEnum : return "Type";

issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -360 +360 @@
               else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
+              else if (strcmp(name,"TemperatureSurface")==0) return TemperatureSurfaceEnum;
+              else if (strcmp(name,"TemperatureBasal")==0) return TemperatureBasalEnum;
               else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
               else if (strcmp(name,"Type")==0) return TypeEnum;
@@ -382 +384 @@
               else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
               else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
-              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
-              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
+              if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
+              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
+              else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
               else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
               else if (strcmp(name,"Outputfilename")==0) return OutputfilenameEnum;

issm/branches/trunk-jpl-damage/src/m/classes/initialization.m

@@ -12 +12 @@
                 pressure      = NaN;
                 temperature   = NaN;
+                surfacetemp   = NaN;
+                basaltemp     = NaN;
                 watercolumn   = NaN;
                 waterfraction = NaN;
@@ -68 +70 @@
                         fielddisplay(obj,'pressure','pressure field');
                         fielddisplay(obj,'temperature','temperature in Kelvins');
+                        fielddisplay(obj,'surfacetemp','surface temperature in Kelvins');
+                        fielddisplay(obj,'basaltemp','basal temperature in Kelvins');
                         fielddisplay(obj,'watercolumn','thickness of subglacial water');
                         fielddisplay(obj,'waterfraction','fraction of water in the ice');
@@ -78 +82 @@
                         WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum);
                         WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum);
+                        WriteData(fid,'data',obj.surfacetemp,'format','DoubleMat','mattype',1,'enum',TemperatureSurfaceEnum);
+                        WriteData(fid,'data',obj.basaltemp,'format','DoubleMat','mattype',1,'enum',TemperatureBasalEnum);
                         WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);
                         WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);

issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py

@@ -3367 +3367 @@
         return StringToEnum('TemperaturePicard')
 
+def TemperatureSurfaceEnum():
+        """
+        TEMPERATURESURFACEENUM - Enum of TemperatureSurface
+
+           Usage:
+              macro=TemperatureSurfaceEnum()
+        """
+
+        return StringToEnum('TemperatureSurface')
+
+def TemperatureBasalEnum():
+        """
+        TEMPERATUREBASALENUM - Enum of TemperatureBasal
+
+           Usage:
+              macro=TemperatureBasalEnum()
+        """
+
+        return StringToEnum('TemperatureBasal')
+
 def ThicknessAbsMisfitEnum():
         """
@@ -4595 +4615 @@
         """
 
-        return 458
-
+        return 460
+

issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m

@@ -9 +9 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=458;
+macro=460;

issm/branches/trunk-jpl-damage/src/m/model/plot/applyoptions.m

@@ -121 +121 @@
 c = getcolormap(options);
 h = colormap(c);
-
-                c = hsv(64);
-                c = rgb2hsv(c);
-                c(:,2) = max(min( abs(c(:,1)-0.5)/0.5 ,1),0);
-                c(1:32,1)   = 0.7;
-                c(33:end,1) = 1;
-                c = hsv2rgb(c);
-
-        elseif strcmpi(cname,'Rignot'),
-                c = hsv;
-
-                %adjust saturation
-                c = rgb2hsv(c);
-                alpha=getfieldvalue(options,'alpha',1);
-                c(:,2) = max(min( (0.1+c(:,1)).^(1/alpha) ,1),0);
-                c = hsv2rgb(c);
-
-        elseif strcmpi(cname,'Rignot2'),
-                c = hsv;
-
-                %adjust saturation
-                c = rgb2hsv(c);
-                alpha=getfieldvalue(options,'alpha',1);
-                c(:,2) = max(min( (0.1+c(:,1)).^(1/alpha) ,1),0);
-                c = hsv2rgb(c);
-
-                c=flipud(c);
-
-        elseif strcmpi(cname,'damage'),
-                c=hsv(64);
-                c(:,3)=1;
-                a=1:-(1/35):0.1;
-                c(1:32,2)=a;
-                c(33:end,2)=fliplr(a);
-                c(1:32,1)=2/3;
-                c(33:end,1)=0.07;
-                c=hsv2rgb(c);
-
-        else
-                c = cname;
-        end
-        h=colormap(c);
-else
-        h=colormap(jet(60));
-end
 
 %wrapping

issm/branches/trunk-jpl-damage/src/m/model/steadystateiceshelftemp.m

@@ -8 +8 @@
 %
 %        In addition to supplying md, the surface and basal temperatures of the ice shelf must
-%        be supplied to the function IN DEGREES CELSIUS.  These temperatures can be supplied at 
-%        each vertex of the mesh or as a single constant value to be used at every vertex.  
+%        be supplied in degrees Kelvin.
 %
 %        The model md must also contain the fields: 
@@ -16 +15 @@
 
 %   Usage:
-%      temperature=steadystateiceshelftemp(md)
+%      temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp)
 
 if (length(md.geometry.thickness)~=md.mesh.numberofvertices)
@@ -23 +22 @@
 
 %surface and basal temperatures in degrees C
-if (length(surfacetemp)==md.mesh.numberofvertices)
-        Ts=surfacetemp;
-elseif (length(surfacetemp)==1) 
-        Ts=surfacetemp*ones(md.mesh.numberofvertices,1);
-else
-        error(['steadystateiceshelftemp error message: surfacetemp should have a length of 1 or ' num2str(md.mesh.numberofvertices)])
+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)
-        Tb=basaltemp;
-elseif (length(basaltemp)==1) 
-        Tb=basaltemp*ones(md.mesh.numberofvertices,1);
-else
-        error(['steadystateiceshelftemp error message: basaltemp should have a length of 1 or ' num2str(md.mesh.numberofvertices)])
+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, positive for melting (thus minus sign)
-wi=-md.materials.rho_water/md.materials.rho_ice.*md.basalforcings.melting_rate; 
+%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
-md.initialization.temperature=(Ts+Tb)/2;  % where wi~=0
+temperature=(Ts+Tb)/2;  % where wi~=0
 
 pos=find(abs(wi)>=1e-4); % to avoid division by zero
 
-%calculate depth-averaged temperature
+%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));