source: issm/oecreview/Archive/20545-21336/ISSM-20742-20743.diff

Index: ../trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par
===================================================================
--- ../trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par    (revision 20742)
+++ ../trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par    (revision 20743)
@@ -1,7 +1,7 @@
-%Name and Coordinate system
+%Name and hemisphere
 md.miscellaneous.name='SeaRISEgreenland';
-md.mesh.epsg=3413;
 
+
 disp('   Loading SeaRISE data from NetCDF');
 ncdata='../Data/Greenland_5km_dev1.2.nc';
 x1    = ncread(ncdata,'x1');
@@ -43,7 +43,7 @@
 bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
 surf_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),surf,xi,yi,NaN);
 in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
-       './Jak_grounded.exp','node',1);
+       'Jak_grounded.exp','node',1);
 bed_jks(~in)=NaN;
 surf_jks(~in)=NaN;
 pos=find(~isnan(bed_jks));
@@ -65,8 +65,8 @@
 md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15;
 
 disp('   Interpolating surface mass balance');
-md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
-md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
+md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
+md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
 
 disp('   Construct basal friction parameters');
 md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
@@ -83,8 +83,7 @@
 
 disp('   Set other boundary conditions');
 md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
-md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
-md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
+md.basalforcings.melting_rate = zeros(md.mesh.numberofvertices,1);
 md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
 md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
 
Index: ../trunk-jpl/examples/IceBridge/Greenland.par
===================================================================
--- ../trunk-jpl/examples/IceBridge/Greenland.par       (revision 20742)
+++ ../trunk-jpl/examples/IceBridge/Greenland.par       (revision 20743)
@@ -26,28 +26,6 @@
 md.geometry.thickness(pos0)=1;
 md.geometry.surface=md.geometry.thickness+md.geometry.base;
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%Reading IceBridge data for Jakobshavn
-disp('      reading IceBridge Jakobshavn bedrock');
-fid  = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
-titles = fgets(fid); data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
-fclose(fid);
-
-[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
-bed  = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
-bedy = flipud(reshape(data(:,1),[360 740]));
-bedx = flipud(reshape(data(:,2),[360 740]));
-
-%Insert Icebridge bed and recalculate thickness
-bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
-in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
-       './Jak_grounded.exp','node',1);
-bed_jks(~in)=NaN;
-pos=find(~isnan(bed_jks));
-md.geometry.base(pos)=bed_jks(pos);
-md.geometry.thickness=md.geometry.surface-md.geometry.base;
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 disp('   Interpolating velocities ');
 md.inversion.vx_obs  = InterpFromGridToMesh(x1,y1,velx,md.mesh.x,md.mesh.y,0);
 md.inversion.vy_obs  = InterpFromGridToMesh(x1,y1,vely,md.mesh.x,md.mesh.y,0);
@@ -79,8 +57,8 @@
 
 disp('   Set other boundary conditions');
 md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
-md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
 md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
+md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
 md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
 md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
 
@@ -90,15 +68,10 @@
 disp('   Set Pressure');
 md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
 
-disp('   Single point constraint for continental model');
-%Set at least one vertex to velocity 0 so as to not get a singular problem (point on the wet peninsula)
+disp('   Single point constraints');
+%Initialize single point constraint arrarys
 md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
 md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices,1);
 md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices,1);
 md.stressbalance.spcvz = NaN*ones(md.mesh.numberofvertices,1);
-location = 1.0e+06 *[.32011 -2.2039];
-[dist pos]=min(sqrt((md.mesh.x - location(1)).^2 + (md.mesh.y - location(2)).^2));
-md.stressbalance.spcvx(pos) = 0;
-md.stressbalance.spcvy(pos) = 0;
-md.stressbalance.spcvz(pos) = 0;
 
Index: ../trunk-jpl/examples/IceBridge/runme.m
===================================================================
--- ../trunk-jpl/examples/IceBridge/runme.m     (revision 20742)
+++ ../trunk-jpl/examples/IceBridge/runme.m     (revision 20743)
@@ -1,4 +1,5 @@
 clear all;
+close all;
 steps=[1];
 
 %Location of SeaRISE dataset
@@ -20,16 +21,13 @@
        vel  = sqrt(vx.^2+vy.^2);
 
        %Mesh greenland without refinement in Jak basin
-       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8);
-       save ./Models/Greenland.Mesh_generation_old md;
        return;
-       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
+       
        %Refine mesh in the region of Jakobshavn (resolution = 3000 m)
        hmaxVertices=NaN*ones(md.mesh.numberofvertices,1);
        in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
-               './Jak_outline.exp','node',1);
+               'Jak_outline.exp','node',1);
        hmaxVertices(find(in))=3000;
        md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,...
                'err',8,'hmaxVertices',hmaxVertices);
@@ -45,8 +43,7 @@
        md = loadmodel('./Models/Greenland.Mesh_generation');
 
        md = setmask(md,'','');
-       md = parameterize(md,'./Greenland.par');
-       %md = parameterize(md,'./Greenland_noOIB.par');
+       md = parameterize(md,'Greenland.par');
        md = setflowequation(md,'SSA','all');
 
        save ./Models/Greenland.Parameterization2 md; 
@@ -76,7 +73,7 @@
        md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
        md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
        in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
-               './data_gaps.exp','node',1);
+               'data_gaps.exp','node',1);
        md.inversion.cost_functions_coefficients(find(in),1)=0.0;
        md.inversion.cost_functions_coefficients(find(in),2)=0.0;
 
@@ -142,18 +139,19 @@
 
        %Line Plots
        figure
+       time_plot=md.results.TransientSolution(1).time:md.timestepping.time_step:md.timestepping.final_time;
 
        %Plot surface mass balance
        surfmb=[]; for i=1:100; surfmb=[surfmb ...
                md.results.TransientSolution(i).SmbMassBalance]; end
-       subplot(3,1,1); plot([0.2:0.2:20],mean(surfmb)); title('Mean Surface mass balance');
+       subplot(3,1,1); plot(time_plot,mean(surfmb)); title('Mean Surface mass balance');
 
        %Plot velocity
        vel=[]; for i=1:100; vel=[vel md.results.TransientSolution(i).Vel]; end
-       subplot(3,1,2); plot([0.2:0.2:20],mean(vel)); title('Mean Velocity');
+       subplot(3,1,2); plot(time_plot,mean(vel)); title('Mean Velocity');
 
        %Plot Volume
        volume=[]; for i=1:100; volume=[volume md.results.TransientSolution(i).IceVolume]; end
-       subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume');
+       subplot(3,1,3); plot(time_plot,volume); title('Ice Volume');
        xlabel('years')
 end
Index: ../trunk-jpl/examples/IceBridge/Greenland_solution.par
===================================================================
--- ../trunk-jpl/examples/IceBridge/Greenland_solution.par      (revision 0)
+++ ../trunk-jpl/examples/IceBridge/Greenland_solution.par      (revision 20743)
@@ -0,0 +1,107 @@
+%Name and hemisphere
+md.miscellaneous.name='SeaRISEgreenland';
+
+
+disp('   Loading SeaRISE data from NetCDF');
+ncdata='../Data/Greenland_5km_dev1.2.nc';
+x1    = ncread(ncdata,'x1');
+y1    = ncread(ncdata,'y1');
+usrf  = ncread(ncdata,'usrf')';
+topg  = ncread(ncdata,'topg')';
+velx  = ncread(ncdata,'surfvelx')';
+vely  = ncread(ncdata,'surfvely')';
+temp  = ncread(ncdata,'airtemp2m')';
+smb   = ncread(ncdata,'smb')';
+gflux = ncread(ncdata,'bheatflx')';
+
+disp('   Interpolating surface and bedrock');
+md.geometry.base     = InterpFromGridToMesh(x1,y1,topg,md.mesh.x,md.mesh.y,0);
+md.geometry.surface = InterpFromGridToMesh(x1,y1,usrf,md.mesh.x,md.mesh.y,0);
+
+disp('   Constructing thickness');
+md.geometry.thickness=md.geometry.surface-md.geometry.base;
+
+%Set min thickness to 1 meter
+pos0=find(md.geometry.thickness<=0);
+md.geometry.thickness(pos0)=1;
+md.geometry.surface=md.geometry.thickness+md.geometry.base;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Reading IceBridge data for Jakobshavn
+disp('      reading IceBridge Jakobshavn bedrock');
+fid  = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
+titles = fgets(fid); data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
+fclose(fid);
+
+[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
+bed  = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
+surf  = flipud(reshape(data(:,4),[360 740])); surf(find(surf==-9999))=NaN;
+bedy = flipud(reshape(data(:,1),[360 740]));
+bedx = flipud(reshape(data(:,2),[360 740]));
+
+%Insert Icebridge bed and recalculate thickness
+bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
+surf_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),surf,xi,yi,NaN);
+in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
+       'Jak_grounded.exp','node',1);
+bed_jks(~in)=NaN;
+surf_jks(~in)=NaN;
+pos=find(~isnan(bed_jks));
+md.geometry.base(pos)=bed_jks(pos);
+md.geometry.surface(pos)=surf_jks(pos);
+md.geometry.thickness=md.geometry.surface-md.geometry.base;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+disp('   Interpolating velocities ');
+md.inversion.vx_obs  = InterpFromGridToMesh(x1,y1,velx,md.mesh.x,md.mesh.y,0);
+md.inversion.vy_obs  = InterpFromGridToMesh(x1,y1,vely,md.mesh.x,md.mesh.y,0);
+md.inversion.vel_obs = sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
+md.initialization.vx = md.inversion.vx_obs;
+md.initialization.vy = md.inversion.vy_obs;
+md.initialization.vz = zeros(md.mesh.numberofvertices,1);
+md.initialization.vel= md.inversion.vel_obs;
+
+disp('   Interpolating temperatures');
+md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15;
+
+disp('   Interpolating surface mass balance');
+md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
+md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
+
+disp('   Construct basal friction parameters');
+md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
+pos=find(md.mask.groundedice_levelset<0);
+md.friction.coefficient(pos)=0; %no friction applied on floating ice
+md.friction.p=ones(md.mesh.numberofelements,1);
+md.friction.q=ones(md.mesh.numberofelements,1);
+
+disp('   Construct ice rheological properties');
+md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_B=paterson(md.initialization.temperature);
+md.friction.q=ones(md.mesh.numberofelements,1);
+md.friction.p=ones(md.mesh.numberofelements,1);
+
+disp('   Set other boundary conditions');
+md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
+md.basalforcings.melting_rate = zeros(md.mesh.numberofvertices,1);
+md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
+md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
+
+disp('   Set geothermal heat flux');
+md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md.mesh.x,md.mesh.y,0);
+
+disp('   Set Pressure');
+md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
+
+disp('   Single point constraint for continental model');
+%Set at least one vertex to velocity 0 so as to not get a singular problem (point on the wet peninsula)
+md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
+md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices,1);
+md.stressbalance.spcvz = NaN*ones(md.mesh.numberofvertices,1);
+location = 1.0e+06 *[.32011 -2.2039];
+[dist pos]=min(sqrt((md.mesh.x - location(1)).^2 + (md.mesh.y - location(2)).^2));
+md.stressbalance.spcvx(pos) = 0;
+md.stressbalance.spcvy(pos) = 0;
+md.stressbalance.spcvz(pos) = 0;
+

Property changes on: ../trunk-jpl/examples/IceBridge/Greenland_solution.par
___________________________________________________________________
Added: svn:executable
   + *