Changeset 20743


Ignore:
Timestamp:
06/15/16 17:56:54 (9 years ago)
Author:
schlegel
Message:

CHG: update for workshop

Location:
issm/trunk-jpl/examples/IceBridge
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/examples/IceBridge/Greenland.par

    r20532 r20743  
    2626md.geometry.thickness(pos0)=1;
    2727md.geometry.surface=md.geometry.thickness+md.geometry.base;
    28 
    29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    30 %Reading IceBridge data for Jakobshavn
    31 disp('      reading IceBridge Jakobshavn bedrock');
    32 fid  = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
    33 titles = fgets(fid); data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
    34 fclose(fid);
    35 
    36 [xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
    37 bed  = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
    38 bedy = flipud(reshape(data(:,1),[360 740]));
    39 bedx = flipud(reshape(data(:,2),[360 740]));
    40 
    41 %Insert Icebridge bed and recalculate thickness
    42 bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
    43 in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
    44         './Jak_grounded.exp','node',1);
    45 bed_jks(~in)=NaN;
    46 pos=find(~isnan(bed_jks));
    47 md.geometry.base(pos)=bed_jks(pos);
    48 md.geometry.thickness=md.geometry.surface-md.geometry.base;
    49 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    5028
    5129disp('   Interpolating velocities ');
    … …  
    8058disp('   Set other boundary conditions');
    8159md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
     60md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
    8261md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
    83 md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
    8462md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
    8563md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
    … …  
    9169md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
    9270
    93 disp('   Single point constraint for continental model');
    94 %Set at least one vertex to velocity 0 so as to not get a singular problem (point on the wet peninsula)
     71disp('   Single point constraints');
     72%Initialize single point constraint arrarys
    9573md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
    9674md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices,1);
    9775md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices,1);
    9876md.stressbalance.spcvz = NaN*ones(md.mesh.numberofvertices,1);
    99 location = 1.0e+06 *[.32011 -2.2039];
    100 [dist pos]=min(sqrt((md.mesh.x - location(1)).^2 + (md.mesh.y - location(2)).^2));
    101 md.stressbalance.spcvx(pos) = 0;
    102 md.stressbalance.spcvy(pos) = 0;
    103 md.stressbalance.spcvz(pos) = 0;
    10477

  • issm/trunk-jpl/examples/IceBridge/Greenland_cheatsheet.par

    r20532 r20743  
    1 %Name and Coordinate system
     1%Name and hemisphere
    22md.miscellaneous.name='SeaRISEgreenland';
    3 md.mesh.epsg=3413;
     3
    44
    55disp('   Loading SeaRISE data from NetCDF');
    … …  
    4444surf_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),surf,xi,yi,NaN);
    4545in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
    46         './Jak_grounded.exp','node',1);
     46        'Jak_grounded.exp','node',1);
    4747bed_jks(~in)=NaN;
    4848surf_jks(~in)=NaN;
    … …  
    6666
    6767disp('   Interpolating surface mass balance');
    68 md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
    69 md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
     68md.surfaceforcings.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
     69md.surfaceforcings.mass_balance=md.surfaceforcings.mass_balance*md.materials.rho_water/md.materials.rho_ice;
    7070
    7171disp('   Construct basal friction parameters');
    … …  
    8484disp('   Set other boundary conditions');
    8585md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
    86 md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
    87 md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
     86md.basalforcings.melting_rate = zeros(md.mesh.numberofvertices,1);
    8887md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
    8988md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);

  • issm/trunk-jpl/examples/IceBridge/runme.m

    r20532 r20743  
    11clear all;
     2close all;
    23steps=[1];
    34
    … …  
    2122
    2223        %Mesh greenland without refinement in Jak basin
    23         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    2424        md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,'err',8);
    25         save ./Models/Greenland.Mesh_generation_old md;
    2625        return;
    27         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    28 
     26       
    2927        %Refine mesh in the region of Jakobshavn (resolution = 3000 m)
    3028        hmaxVertices=NaN*ones(md.mesh.numberofvertices,1);
    3129        in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
    32                 './Jak_outline.exp','node',1);
     30                'Jak_outline.exp','node',1);
    3331        hmaxVertices(find(in))=3000;
    3432        md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel,...
    … …  
    4644
    4745        md = setmask(md,'','');
    48         md = parameterize(md,'./Greenland.par');
    49         %md = parameterize(md,'./Greenland_noOIB.par');
     46        md = parameterize(md,'Greenland.par');
    5047        md = setflowequation(md,'SSA','all');
    5148
    … …  
    7774        md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
    7875        in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
    79                 './data_gaps.exp','node',1);
     76                'data_gaps.exp','node',1);
    8077        md.inversion.cost_functions_coefficients(find(in),1)=0.0;
    8178        md.inversion.cost_functions_coefficients(find(in),2)=0.0;
    … …  
    143140        %Line Plots
    144141        figure
     142        time_plot=md.results.TransientSolution(1).time:md.timestepping.time_step:md.timestepping.final_time;
    145143
    146144        %Plot surface mass balance
    147145        surfmb=[]; for i=1:100; surfmb=[surfmb ...
    148146                md.results.TransientSolution(i).SmbMassBalance]; end
    149         subplot(3,1,1); plot([0.2:0.2:20],mean(surfmb)); title('Mean Surface mass balance');
     147        subplot(3,1,1); plot(time_plot,mean(surfmb)); title('Mean Surface mass balance');
    150148
    151149        %Plot velocity
    152150        vel=[]; for i=1:100; vel=[vel md.results.TransientSolution(i).Vel]; end
    153         subplot(3,1,2); plot([0.2:0.2:20],mean(vel)); title('Mean Velocity');
     151        subplot(3,1,2); plot(time_plot,mean(vel)); title('Mean Velocity');
    154152
    155153        %Plot Volume
    156154        volume=[]; for i=1:100; volume=[volume md.results.TransientSolution(i).IceVolume]; end
    157         subplot(3,1,3); plot([0.2:0.2:20],volume); title('Ice Volume');
     155        subplot(3,1,3); plot(time_plot,volume); title('Ice Volume');
    158156        xlabel('years')
    159157end
Note: See TracChangeset for help on using the changeset viewer.