== Parameterization file == {{{ #!m %Paths to datasets modeldatapath='/u/astrid-r1b/larour/ModelData'; thicknesspath=[modeldatapath '/ASEThickness/ASE_thickness']; surfacepath=[modeldatapath '/BamberDEMAntarctica1km/surface_smooth30_lowslope']; firnpath=[modeldatapath '/BroekeFirn1km/firn.mat']; %Read Surface and interpolate onto model md.surface=InterpFromFile(md.x,md.y,surfacepath,10); %Correct surface using firn layerfrom Broeke md.firn_layer=InterpFromFile(md.x,md.y,firnpath,0); rho_ice=917; rho_firn=830; firn_layer_correction=md.firn_layer*(1-rho_firn/rho_ice); md.surface=md.surface-firn_layer_correction; %Load thickness and interpolate onto model load(thicknesspath); md.thickness=InterpFromGridToMesh(x_m,y_m,thickness,md.x,md.y,10); %Correct thickness using Bamber's dem md=ThicknessCorrection(md,15000); }}} == Thickness correction routine == {{{ #!m function md=ThicknessCorrection(md,distance) %THICKNESSCORRECTION - correct the thickness of the ice shelf near the grounding line % % This routine corrects the thickness and the bed on the transition zone % by forcing the hydrostatic equilibrium. % the thickness is modified as follows: % thickness = coeff * thickness_observation + (1-coeff) * thickness_hydrostatic % where: % coeff=(d/distance)^2; % distance=10km by default but can be specified % % Usage: % md=ThicknessCorrection(md,distance); % % Example: % md=ThicknessCorrection(md,15000); %initialize thickness with the observations, and get hydrostatic thickness from the dem thickness=md.thickness; thickness_hydro=md.surface/(1-md.rho_ice/md.rho_water); %get nodes on ice sheet and on ice shelf pos_shelf=find(~md.gridonicesheet); pos_GL=intersect(unique(md.elements(find(md.elementonicesheet),:)),unique(md.elements(find(md.elementoniceshelf),:))); %modify thickness for i=1:length(pos_shelf) %search the grid on ice sheet the closest to i [d posd]=min(sqrt((md.x(pos_shelf(i))-md.x(pos_GL)).^2+(md.y(pos_shelf(i))-md.y(pos_GL)).^2)); if d>distance, %if d > 15km, hydrostatique equilibrium thickness(pos_shelf(i))=thickness_hydro(pos_shelf(i)); else %else: quadratic combination of hydrostatic equilibrium and observations (this line might be commented if the original thickness is good close to the GL) coeff=(d/distance)^2; thickness(pos_shelf(i))=(1-coeff)*thickness(pos_shelf(i))+coeff*thickness_hydro(pos_shelf(i)); end end %check the computed thickness minth=1/(1-md.rho_ice/md.rho_water); pos=find(isnan(thickness) | (thickness<=0)); thickness(pos)=minth; %change bed to take into account the changes in thickness md.thickness=thickness; md.bed=md.surface-md.thickness; }}} == stereographic conversion part 1 == {{{ #!m function [x,y] = map2ll(lat,lon,sgn) % % NAME: map2ll %=============================================================================== % % USAGE: [x,y] = map2ll(lat,lon,sgn) % % DESCRIPTION: % Converts from geodetic latitude and longitude to Polar % Stereographic (X,Y) coordinates for the polar regions. % % INPUT: % Lat = Geodetic Latitude (degrees, +90 to -90) % Lon = Geodetic Longitude (degrees, 0 to 360) % SGN = Sign of latitude % +1 : north latitude % -1 : south latitude % % % OUTPUT: X = Polar Stereographic X Coordinate (km) % Y = Polar Stereographic Y Coordinate (km) % % % AUTHOR: The original equations are from Snyder, J. P., % 1982, Map Projections Used by the U.S. Geological Survey, % Geological Survey Bulletin 1532, U.S. Government Printing % Office. See JPL Technical Memorandum 3349-85-101 for % further details. % % Original FORTRAN program written by C. S. Morris, % April 1985, Jet Propulsion Laboratory, California % Institute of Technology % % IDL conversion by Helmut Schottmueller, September 1995, % Institute for environmental physics, University of Bremen % % Matlab conversion by Michael P. Schodlok, December 2003, % Alfred Wegener Institute for Polar and Marine Research, Bremerhaven % % Error correction for RHO: 1997-11-03 by LML %==================================================================================== %------------- % CHECK INPUTS %------------- if nargin ~= 3 error('map2ll.m: must pass 3 parameters lat, lon, hemisphere') end %if % CHECK lat,lon,hems dimensions and verify consistent [ma,na] = size(lat); [mb,nb] = size(lon); [mc,nc] = size(sgn); % CHECK THAT lat & lon HAVE SAME SHAPE if (ma~=mb) | (na~=nb) error('check dims: lat & lon must have same dimensions') end %if % CHECK THAT sgn IS A SCALAR if (mc~=1 & np~=1 ) error('check sgn: sgn is either 1 or -1') end %if %====== % BEGIN %====== % Conversion constant from degrees to radians cde = 57.29577951; % Standard latitude for the SSM/I grid with no distorsion slat = 70.; % Radius of the earth in kilometers re = 6378.273; % Eccentricity of the Hughes ellipsoid squared ex2 = .006693883; % Eccentricity of the Hughes ellipsoid ex = sqrt(ex2); if sgn == 1 delta = 45.; else delta = 0.0; end %if latitude = abs(lat) * pi/180.; longitude = (lon + delta) * pi/180.; % compute X and Y in grid coordinates. T = tan(pi/4-latitude/2) ./ ((1-ex*sin(latitude))./(1+ex*sin(latitude))).^(ex/2); if (90 - slat) < 1.e-5 rho = 2.*re*T/sqrt((1.+ex)^(1.+ex)*(1.-ex)^(1.-ex)); else sl = slat*pi/180.; tc = tan(pi/4.-sl/2.)/((1.-ex*sin(sl))/(1.+ex*sin(sl)))^(ex/2.); mc = cos(sl)/sqrt(1.0-ex2*(sin(sl)^2)); rho = re*mc*T/tc; end % if y = -rho .* sgn .* cos(sgn.*longitude); x = rho .* sgn .* sin(sgn.*longitude); [cnt1,cnt2] = find(latitude >= pi / 2.); if cnt1 x(cnt1,1) = 0.0 y(cnt1,1) = 0.0 end %if return %========================================================================= }}} == stereographic conversion part 2 == {{{ #!m function [lat,lon] = map2xy(x,y,sgn) % % NAME: map2xy %=============================================================================== % % USAGE: [lat,lon] = map2xy(x,y,sgn) % % DESCRIPTION: % Converts from geodetic latitude and longitude to Polar % Stereographic (X,Y) coordinates for the polar regions. % % INPUT: % X = Polar Stereographic X Coordinate (km) % Y = Polar Stereographic Y Coordinate (km) % SGN = Sign of latitude % +1 : north latitude % -1 : south latitude % % % OUTPUT: % Lat = Geodetic Latitude (degrees, +90 to -90) % Lon = Geodetic Longitude (degrees, 0 to 360) % % % AUTHOR: The original equations are from Snyder, J. P., % 1982, Map Projections Used by the U.S. Geological Survey, % Geological Survey Bulletin 1532, U.S. Government Printing % Office. See JPL Technical Memorandum 3349-85-101 for % further details. % % Original FORTRAN program written by C. S. Morris, % April 1985, Jet Propulsion Laboratory, California % Institute of Technology % % IDL conversion by Helmut Schottmueller, September 1995, % Institute for environmental physics, University of Bremen % % Matlab conversion by Michael P. Schodlok, December 2003, % Alfred Wegener Institute for Polar and Marine Research, Bremerhaven % % Error correction for RHO: 1997-11-03 by LML %==================================================================================== %------------- % CHECK INPUTS %------------- if nargin ~= 3 error('map2ll.m: must pass 3 parameters lat, lon, hemisphere') end %if % CHECK lat,lon,hems dimensions and verify consistent [ma,na] = size(x); [mb,nb] = size(y); [mc,nc] = size(sgn); % CHECK THAT lat & lon HAVE SAME SHAPE if (ma~=mb) | (na~=nb) error('check dims: lat & lon must have same dimensions') end %if % CHECK THAT sgn IS A SCALAR if (mc~=1 & np~=1 ) error('check sgn: sgn is either 1 or -1') end %if %====== % BEGIN %====== % Conversion constant from degrees to radians cde = 57.29577951; % Standard latitude for the SSM/I grid with no distorsion slat = 70.; % Radius of the earth in kilometers re = 6378.273; % Eccentricity of the Hughes ellipsoid squared ex2 = .006693883; % Eccentricity of the Hughes ellipsoid ex = sqrt(ex2); if sgn == 1 delta = 45.; else delta = 0.0; end %if sl = slat*pi/180.; rho = sqrt(x.^2 + y.^2); cm = cos(sl) / sqrt(1.0 - ex2 * (sin(sl)^2)); T = tan((pi / 4.0) - (sl / 2.0)) / ((1.0 - ex * sin(sl)) / (1.0 + ex * sin(sl)))^(ex / 2.0); if abs(slat-90.) < 1.e-5 T = rho * sqrt((1. + ex)^(1. + ex) * (1. - ex)^(1. - ex)) / 2. / re; else T = rho * T / (re * cm); end %if chi = (pi / 2.0) - 2.0 * atan(T); lat = chi + ((ex2 / 2.0) + (5.0 * ex2^2.0 / 24.0) + (ex2^3.0 / 12.0)) * ... sin(2 * chi) + ((7.0 * ex2^2.0 / 48.0) + (29.0 * ex2^3 / 240.0)) * ... sin(4.0 * chi) + (7.0 * ex2^3.0 / 120.0) * sin(6.0 * chi) ; lat = sgn * lat; lon = atan2(sgn * x,-sgn * y); lon = sgn * lon; [res1,res2] = find(rho <= 0.1); if res1 lat(res1,1) = 90. * sgn; lon(res1,1) = 0.0; end %if lon = lon * 180. / pi; lat = lat * 180. / pi; lon = lon - delta; return %========================================================================= }}}