Index: /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumDefinitions.h (revision 12917) @@ -362,4 +362,6 @@ TemperatureOldEnum, TemperaturePicardEnum, + TemperatureSurfaceEnum, + TemperatureBasalEnum, ThicknessAbsMisfitEnum, TypeEnum, Index: /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12917) @@ -353,4 +353,6 @@ case TemperatureOldEnum : return "TemperatureOld"; case TemperaturePicardEnum : return "TemperaturePicard"; + case TemperatureSurfaceEnum : return "TemperatureSurface"; + case TemperatureBasalEnum : return "TemperatureBasal"; case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit"; case TypeEnum : return "Type"; Index: /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12917) @@ -360,4 +360,6 @@ 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,10 +384,10 @@ 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; Index: /issm/branches/trunk-jpl-damage/src/m/classes/initialization.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/classes/initialization.m (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/m/classes/initialization.m (revision 12917) @@ -12,4 +12,6 @@ pressure = NaN; temperature = NaN; + surfacetemp = NaN; + basaltemp = NaN; watercolumn = NaN; waterfraction = NaN; @@ -68,4 +70,6 @@ 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,4 +82,6 @@ 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); Index: /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/m/enum/EnumDefinitions.py (revision 12917) @@ -3367,4 +3367,24 @@ 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,4 +4615,4 @@ """ - return 458 - + return 460 + Index: /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/m/enum/MaximumNumberOfEnums.m (revision 12917) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=458; +macro=460; Index: /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureBasalEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureBasalEnum.m (revision 12917) +++ /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureBasalEnum.m (revision 12917) @@ -0,0 +1,11 @@ +function macro=TemperatureBasalEnum() +%TEMPERATUREBASALENUM - Enum of TemperatureBasal +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=TemperatureBasalEnum() + +macro=StringToEnum('TemperatureBasal'); Index: /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureSurfaceEnum.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureSurfaceEnum.m (revision 12917) +++ /issm/branches/trunk-jpl-damage/src/m/enum/TemperatureSurfaceEnum.m (revision 12917) @@ -0,0 +1,11 @@ +function macro=TemperatureSurfaceEnum() +%TEMPERATURESURFACEENUM - Enum of TemperatureSurface +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/EnumDefinitions/Synchronize.sh +% Please read src/c/EnumDefinitions/README for more information +% +% Usage: +% macro=TemperatureSurfaceEnum() + +macro=StringToEnum('TemperatureSurface'); Index: /issm/branches/trunk-jpl-damage/src/m/model/plot/applyoptions.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/plot/applyoptions.m (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/m/model/plot/applyoptions.m (revision 12917) @@ -121,49 +121,4 @@ 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 Index: /issm/branches/trunk-jpl-damage/src/m/model/steadystateiceshelftemp.m =================================================================== --- /issm/branches/trunk-jpl-damage/src/m/model/steadystateiceshelftemp.m (revision 12916) +++ /issm/branches/trunk-jpl-damage/src/m/model/steadystateiceshelftemp.m (revision 12917) @@ -8,6 +8,5 @@ % % 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,5 +15,5 @@ % Usage: -% temperature=steadystateiceshelftemp(md) +% temperature=steadystateiceshelftemp(md,surfacetemp,basaltemp) if (length(md.geometry.thickness)~=md.mesh.numberofvertices) @@ -23,32 +22,28 @@ %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));