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

Last change on this file since 21337 was 21337, checked in by Mathieu Morlighem, 8 years ago

CHG: added Archive/20545-21336

File size: 12.7 KB

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

     
    1 %Name and Coordinate system
     1%Name and hemisphere
    22md.miscellaneous.name='SeaRISEgreenland';
    3 md.mesh.epsg=3413;
    43
     4
    55disp('   Loading SeaRISE data from NetCDF');
    66ncdata='../Data/Greenland_5km_dev1.2.nc';
    77x1    = ncread(ncdata,'x1');
    … …  
    4343bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
    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;
    4949pos=find(~isnan(bed_jks));
    … …  
    6565md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15;
    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');
    7272md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
    … …  
    8383
    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);
    9089

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

     
    2626md.geometry.thickness(pos0)=1;
    2727md.geometry.surface=md.geometry.thickness+md.geometry.base;
    2828
    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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    50 
    5129disp('   Interpolating velocities ');
    5230md.inversion.vx_obs  = InterpFromGridToMesh(x1,y1,velx,md.mesh.x,md.mesh.y,0);
    5331md.inversion.vy_obs  = InterpFromGridToMesh(x1,y1,vely,md.mesh.x,md.mesh.y,0);
    … …  
    7957
    8058disp('   Set other boundary conditions');
    8159md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
    82 md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
    8360md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
     61md.basalforcings.groundedice_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);
    8664
    … …  
    9068disp('   Set Pressure');
    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

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

     
    11clear all;
     2close all;
    23steps=[1];
    34
    45%Location of SeaRISE dataset
    … …  
    2021        vel  = sqrt(vx.^2+vy.^2);
    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,...
    3533                'err',8,'hmaxVertices',hmaxVertices);
    … …  
    4543        md = loadmodel('./Models/Greenland.Mesh_generation');
    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
    5249        save ./Models/Greenland.Parameterization2 md;
    … …  
    7673        md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
    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;
    8279
    … …  
    142139
    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

  • ../trunk-jpl/examples/IceBridge/Greenland_solution.par

     
     1%Name and hemisphere
     2md.miscellaneous.name='SeaRISEgreenland';
     3
     4
     5disp('   Loading SeaRISE data from NetCDF');
     6ncdata='../Data/Greenland_5km_dev1.2.nc';
     7x1    = ncread(ncdata,'x1');
     8y1    = ncread(ncdata,'y1');
     9usrf  = ncread(ncdata,'usrf')';
     10topg  = ncread(ncdata,'topg')';
     11velx  = ncread(ncdata,'surfvelx')';
     12vely  = ncread(ncdata,'surfvely')';
     13temp  = ncread(ncdata,'airtemp2m')';
     14smb   = ncread(ncdata,'smb')';
     15gflux = ncread(ncdata,'bheatflx')';
     16
     17disp('   Interpolating surface and bedrock');
     18md.geometry.base     = InterpFromGridToMesh(x1,y1,topg,md.mesh.x,md.mesh.y,0);
     19md.geometry.surface = InterpFromGridToMesh(x1,y1,usrf,md.mesh.x,md.mesh.y,0);
     20
     21disp('   Constructing thickness');
     22md.geometry.thickness=md.geometry.surface-md.geometry.base;
     23
     24%Set min thickness to 1 meter
     25pos0=find(md.geometry.thickness<=0);
     26md.geometry.thickness(pos0)=1;
     27md.geometry.surface=md.geometry.thickness+md.geometry.base;
     28
     29%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     30%Reading IceBridge data for Jakobshavn
     31disp('      reading IceBridge Jakobshavn bedrock');
     32fid  = fopen('../Data/Jakobshavn_2008_2011_Composite_XYZGrid.txt');
     33titles = fgets(fid); data = fscanf(fid,'%g,%g,%g,%g,%g',[5 266400])';
     34fclose(fid);
     35
     36[xi,yi]= ll2xy(md.mesh.lat,md.mesh.long,+1,45,70);
     37bed  = flipud(reshape(data(:,5),[360 740])); bed(find(bed==-9999))=NaN;
     38surf  = flipud(reshape(data(:,4),[360 740])); surf(find(surf==-9999))=NaN;
     39bedy = flipud(reshape(data(:,1),[360 740]));
     40bedx = flipud(reshape(data(:,2),[360 740]));
     41
     42%Insert Icebridge bed and recalculate thickness
     43bed_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),bed,xi,yi,NaN);
     44surf_jks=InterpFromGridToMesh(bedx(1,:)',bedy(:,1),surf,xi,yi,NaN);
     45in=ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,...
     46        'Jak_grounded.exp','node',1);
     47bed_jks(~in)=NaN;
     48surf_jks(~in)=NaN;
     49pos=find(~isnan(bed_jks));
     50md.geometry.base(pos)=bed_jks(pos);
     51md.geometry.surface(pos)=surf_jks(pos);
     52md.geometry.thickness=md.geometry.surface-md.geometry.base;
     53%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     54
     55disp('   Interpolating velocities ');
     56md.inversion.vx_obs  = InterpFromGridToMesh(x1,y1,velx,md.mesh.x,md.mesh.y,0);
     57md.inversion.vy_obs  = InterpFromGridToMesh(x1,y1,vely,md.mesh.x,md.mesh.y,0);
     58md.inversion.vel_obs = sqrt(md.inversion.vx_obs.^2+md.inversion.vy_obs.^2);
     59md.initialization.vx = md.inversion.vx_obs;
     60md.initialization.vy = md.inversion.vy_obs;
     61md.initialization.vz = zeros(md.mesh.numberofvertices,1);
     62md.initialization.vel= md.inversion.vel_obs;
     63
     64disp('   Interpolating temperatures');
     65md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15;
     66
     67disp('   Interpolating surface mass balance');
     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;
     70
     71disp('   Construct basal friction parameters');
     72md.friction.coefficient=30*ones(md.mesh.numberofvertices,1);
     73pos=find(md.mask.groundedice_levelset<0);
     74md.friction.coefficient(pos)=0; %no friction applied on floating ice
     75md.friction.p=ones(md.mesh.numberofelements,1);
     76md.friction.q=ones(md.mesh.numberofelements,1);
     77
     78disp('   Construct ice rheological properties');
     79md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
     80md.materials.rheology_B=paterson(md.initialization.temperature);
     81md.friction.q=ones(md.mesh.numberofelements,1);
     82md.friction.p=ones(md.mesh.numberofelements,1);
     83
     84disp('   Set other boundary conditions');
     85md.mask.ice_levelset(md.mesh.vertexonboundary==1)=0;
     86md.basalforcings.melting_rate = zeros(md.mesh.numberofvertices,1);
     87md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
     88md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
     89
     90disp('   Set geothermal heat flux');
     91md.basalforcings.geothermalflux=InterpFromGridToMesh(x1,y1,gflux,md.mesh.x,md.mesh.y,0);
     92
     93disp('   Set Pressure');
     94md.initialization.pressure=md.materials.rho_ice*md.constants.g*md.geometry.thickness;
     95
     96disp('   Single point constraint for continental model');
     97%Set at least one vertex to velocity 0 so as to not get a singular problem (point on the wet peninsula)
     98md.stressbalance.referential=NaN*ones(md.mesh.numberofvertices,6);
     99md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices,1);
     100md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices,1);
     101md.stressbalance.spcvz = NaN*ones(md.mesh.numberofvertices,1);
     102location = 1.0e+06 *[.32011 -2.2039];
     103[dist pos]=min(sqrt((md.mesh.x - location(1)).^2 + (md.mesh.y - location(2)).^2));
     104md.stressbalance.spcvx(pos) = 0;
     105md.stressbalance.spcvy(pos) = 0;
     106md.stressbalance.spcvz(pos) = 0;
     107
Note: See TracBrowser for help on using the repository browser.