Index: /issm/trunk-jpl/examples/IceBridge/Greenland.par =================================================================== --- /issm/trunk-jpl/examples/IceBridge/Greenland.par (revision 20742) +++ /issm/trunk-jpl/examples/IceBridge/Greenland.par (revision 20743) @@ -26,26 +26,4 @@ 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 '); @@ -80,6 +58,6 @@ disp(' Set other boundary conditions'); md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0; +md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1); md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1); -md.basalforcings.floatingice_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); @@ -91,14 +69,9 @@ 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: /issm/trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par =================================================================== --- /issm/trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par (revision 20742) +++ /issm/trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par (revision 20743) @@ -1,5 +1,5 @@ -%Name and Coordinate system +%Name and hemisphere md.miscellaneous.name='SeaRISEgreenland'; -md.mesh.epsg=3413; + disp(' Loading SeaRISE data from NetCDF'); @@ -44,5 +44,5 @@ 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; @@ -66,6 +66,6 @@ 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'); @@ -84,6 +84,5 @@ 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: /issm/trunk-jpl/examples/IceBridge/Greenland_solution.par =================================================================== --- /issm/trunk-jpl/examples/IceBridge/Greenland_solution.par (revision 20743) +++ /issm/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; + Index: /issm/trunk-jpl/examples/IceBridge/runme.m =================================================================== --- /issm/trunk-jpl/examples/IceBridge/runme.m (revision 20742) +++ /issm/trunk-jpl/examples/IceBridge/runme.m (revision 20743) @@ -1,3 +1,4 @@ clear all; +close all; steps=[1]; @@ -21,14 +22,11 @@ %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,... @@ -46,6 +44,5 @@ md = setmask(md,'',''); - md = parameterize(md,'./Greenland.par'); - %md = parameterize(md,'./Greenland_noOIB.par'); + md = parameterize(md,'Greenland.par'); md = setflowequation(md,'SSA','all'); @@ -77,5 +74,5 @@ 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; @@ -143,17 +140,18 @@ %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