Parameterization file

%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

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

  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

   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
%=========================================================================