Index: /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par =================================================================== --- /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par (revision 14144) +++ /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par (revision 14145) @@ -9,4 +9,7 @@ disp(' Loading SeaRISE data from NetCDF'); ncdata='Data/Antarctica_5km_withshelves_v0.75.nc'; +if ~exist(ncdata,'file'), + error('Download Antarctica_5km_withshelves_v0.75.nc first on http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica and put it in Data/'); +end x1 = ncread(ncdata,'x1'); y1 = ncread(ncdata,'y1'); @@ -19,4 +22,7 @@ disp(' Loading velocities data from NetCDF'); nsidc_vel='Data/Antarctica_ice_velocity.nc'; +if ~exist(nsidc_vel,'file'), + error('Download Antarctica_ice_velocity.nc first on http://nsidc.org/data/nsidc-0484.html and put it in Data/'); +end xmin = ncreadatt(nsidc_vel,'/','xmin'); Index: /issm/workshop/2012/Talks/09_PIG_Application/runme.m =================================================================== --- /issm/workshop/2012/Talks/09_PIG_Application/runme.m (revision 14144) +++ /issm/workshop/2012/Talks/09_PIG_Application/runme.m (revision 14145) @@ -5,11 +5,13 @@ % steps=[6 7]; -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% SeaRISe dataset for Antarctica. Can be downloaded at: -% http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica searise='./Data/Antarctica_5km_withshelves_v0.75.nc'; - -% Rignot et al, 2011. NSIDC available velocities (at http://nsidc.org/data/nsidc-0484.html) +if ~exist(searise,'file'), + error('Download Antarctica_5km_withshelves_v0.75.nc first on http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica and put it in Data/'); +end + nsidc_vel='./Data/Antarctica_ice_velocity.nc'; +if ~exist(nsidc_vel,'file'), + error('Download Antarctica_ice_velocity.nc first on http://nsidc.org/data/nsidc-0484.html and put it in Data/'); +end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Index: sm/workshop/2012/Talks/10_PIG_UQ/Exp_Par/DomainOutline.exp =================================================================== --- /issm/workshop/2012/Talks/10_PIG_UQ/Exp_Par/DomainOutline.exp (revision 14144) +++ (revision ) @@ -1,126 +1,0 @@ -## Name: -## Icon:0 -# Points Count Value -121 25000. -# X pos Y pos --1712113.01792813 -349656.020505649 --1711509.19177597 -349354.107429569 --1709240.65861084 -349354.107429569 --1707080.15083453 -348922.005874306 --1705459.7700023 -347193.599653257 --1703083.21144836 -345249.142654577 --1701138.75444968 -344168.888766421 --1698438.11972929 -343952.83798879 --1698006.01817402 -341900.355601294 --1693793.02801022 -340279.97476906 --1691092.39328983 -339631.822436167 --1688391.75856944 -339415.771658536 --1686447.30157076 -340063.923991429 --1683530.61607274 -341036.152490769 --1682558.3875734 -343304.685655896 --1682134.72094546 -344182.148765835 --1682126.28601813 -345681.244209839 --1681072.04449632 -345932.439387964 --1680634.47184078 -346682.563940305 --1680571.96146142 -348557.875321158 --1677446.44249333 -348370.344183073 --1675508.62073312 -349495.531011584 --1676616.99118854 -353567.097593376 --1672991.86018721 -353799.617575363 --1668839.16319382 -353999.199148638 --1667132.22989864 -355121.465154142 --1666069.55344949 -357496.859569889 --1662944.03448141 -358309.494501592 --1660756.17120374 -358747.067157124 --1657922.98360981 -361682.807423854 --1656579.222847 -364003.848741 --1653004.88416289 -364748.063575853 --1653442.45681842 -361622.544607765 --1651129.57278203 -355684.058568398 --1648929.292881 -350845.969201 --1646566.31508863 -351433.352771799 --1646628.82546799 -352558.53960031 --1645441.12826011 -352245.987703502 --1643753.34801735 -353058.622635204 --1643938.56948073 -347802.280478547 --1646260.76983598 -347802.280478547 --1644571.89685034 -341785.670467234 --1643503.3064999 -340869.098659662 --1643190.75460309 -336680.903242424 --1642690.6715682 -335243.164517104 --1641252.93284288 -335430.695655189 --1640440.29791117 -334993.122999657 --1640127.74601436 -334180.488067954 --1639023.7827689 -331719.786283476 --1637183.63212881 -330493.019190082 --1635241.2508976 -328652.868549991 --1631765.41079965 -327017.179092133 --1625836.03651491 -325994.873180972 --1619599.97045683 -325892.642589856 --1615204.05503883 -326097.103772088 --1611319.29257642 -325994.873180972 --1608456.83602517 -326199.334363204 --1603754.22883383 -327323.870865482 --1600789.54169146 -329675.174461153 --1596383.29436089 -331974.149040973 --1592314.06481444 -335534.724894122 --1589259.83915 -338300.084057 --1588627.08983285 -343572.995032775 --1588056.60000602 -345432.35635307 --1586993.92355687 -347432.688492646 --1589259.83915 -349315.983207 --1585493.67445219 -352996.112255843 --1585118.61217602 -354496.361360525 --1585861.91571639 -356240.334251376 --1586368.81976325 -358747.067157124 --1587889.09727008 -360996.600678795 --1586681.37166006 -363935.42864415 --1585806.226349 -365748.229645641 --1586743.88203942 -368373.665578835 --1586056.26786644 -370123.956200964 --1585493.67445219 -371124.122270752 --1586243.79900453 -373562.027065861 --1428620.24729973 -364425.259722667 --1414257.66839347 -327460.548254126 --1412220.43513781 -325972.894202318 --1404515.92405384 -296165.173727996 --1387984.93322852 -267002.010668232 --1375329.9306695 -287124.577303688 --1374139.7661228 -286331.134272554 --1347030.46255906 -322432.792189152 --1351262.15872511 -326267.766839633 --1381809.71542376 -297174.855698053 --1381677.47491858 -295852.450646162 --1382481.1222932 -294804.43614984 --1398663.57133811 -316073.230869517 --1409694.46978971 -363328.212840853 --1202495.37465859 -351317.761956718 --1191463.48094 -350233.974803 --1232161.108357 -292094.507065 --1231243.116761 -288422.540681 --1360067.937381 -107884.193494 --1423274.349742 -81217.925965 --1402446.798767 -45916.992107 --1478343.80656 60338.818803 --1513291.731079 4210.33397000001 --1523175.992559 11976.539419 --1503760.478937 36334.18378 --1528509.56197759 63667.6809196406 --1623055.49187657 62841.1506965373 --1648011.301285 32874.517873 --1638219.390929 20634.629928 --1650153.281676 31038.534681 --1725734.589736 -61372.619303 --1719614.645763 -69022.549269 --1759088.284386 -114004.137466 --1755416.318002 -117064.109453 --1761230.264776 -126856.019808 --1758170.29279 -130527.986192 --1762148.256372 -135423.94137 --1764902.231159 -131445.977788 --1769798.186337 -134199.952575 --1775000.138714 -131445.977788 --1794277.962227 -153477.776089 --1712576.710195 -225081.120566 --1729939.58914591 -335652.456760472 --1712113.01792813 -349656.020505649 Index: sm/workshop/2012/Talks/10_PIG_UQ/Exp_Par/Pig.par =================================================================== --- /issm/workshop/2012/Talks/10_PIG_UQ/Exp_Par/Pig.par (revision 14144) +++ (revision ) @@ -1,92 +1,0 @@ -%Name and hemisphere -md.miscellaneous.name='PIG'; -md.mesh.hemisphere='s'; - -disp(' Loading SeaRISE data from NetCDF'); -ncdata='Data/Antarctica_5km_withshelves_v0.75.nc'; -x1 = ncread(ncdata,'x1'); -y1 = ncread(ncdata,'y1'); -usrf = ncread(ncdata,'usrf')'; -topg = ncread(ncdata,'topg')'; -temp = ncread(ncdata,'presartm')'; -smb = ncread(ncdata,'presprcp')'; -gflux = ncread(ncdata,'bheatflx_fox')'; - -disp(' Loading velocities data from NetCDF'); -nsidc_vel='Data/Antarctica_ice_velocity.nc'; - -xmin = ncreadatt(nsidc_vel,'/','xmin'); -xmin = strtrim(xmin); %this is a string, and we need to recover the double value -xmin = xmin(1:end-2); %get rid of the unit -xmin = str2num(xmin); %convert to double -ymax = ncreadatt(nsidc_vel,'/','ymax'); ymax = strtrim(ymax); ymax = ymax(1:end-2); ymax = str2num(ymax); -nx = ncreadatt(nsidc_vel,'/','nx'); -ny = ncreadatt(nsidc_vel,'/','ny'); -spacing = ncreadatt(nsidc_vel,'/','spacing'); spacing = strtrim(spacing); spacing = spacing(1:end-2); spacing = str2num(spacing); -velx = double(ncread(nsidc_vel,'vx')); -vely = double(ncread(nsidc_vel,'vy')); -x2=xmin+(0:1:nx)'*spacing; x2=double(x2); -y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; y2=double(y2); - -disp(' Interpolating surface and bedrock'); -md.geometry.bed = 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.bed; - -%ensure hydrostatic equilibrium on ice shelf: -di=md.materials.rho_ice/md.materials.rho_water; -pos=find(md.mask.vertexonfloatingice); -md.geometry.thickness(pos)=1/(1-di)*md.geometry.surface(pos); -md.geometry.bed(pos)=md.geometry.surface(pos)-md.geometry.thickness(pos); -md.geometry.hydrostatic_ratio=ones(md.mesh.numberofvertices,1); - -%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.bed; - -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=10*ones(md.mesh.numberofvertices,1); -pos=find(md.mask.vertexonfloatingice); -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 observed velocities') -vx_obs=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0); -vy_obs=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0); -vel_obs=sqrt(vx_obs.^2+vy_obs.^2); -md.initialization.vx=vx_obs; -md.initialization.vy=vy_obs; -md.initialization.vz=zeros(md.mesh.numberofvertices,1); -md.initialization.vel=vel_obs; -md.inversion.vx_obs=vx_obs; -md.inversion.vy_obs=vy_obs; -md.inversion.vel_obs=vel_obs; - -disp(' Set boundary conditions'); -md.basalforcings.melting_rate = zeros(md.mesh.numberofvertices,1); -md.thermal.spctemperature = [md.initialization.temperature;1]; %impose observed temperature on surface -md.prognostic.spcthickness = NaN*ones(md.mesh.numberofvertices,1); -md=SetMarineIceSheetBC(md); - -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.surface-md.mesh.z); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.m =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.m (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.m (revision 14145) @@ -0,0 +1,40 @@ +%Create model +md=model; + +%Create mesh +... + +%Define mask +... + +%Parameterize model +... + +%Extrude to have a 3d model +... + +%Choose flow approximation +... + +%Create dirichlet conditions +md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); +md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); +md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); +%Zero velocity on bed +... + +%Create MPCs to have periodic boundary conditions +posx=find(md.mesh.x==0); +posx2=find(md.mesh.x==max(md.mesh.x)); +posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same grids two times +posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); +md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; + +%Set cluster +... + +%Print messages in verbose +... + +%Solve +... Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.par =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.par (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPA.par (revision 14145) @@ -0,0 +1,20 @@ +%Create geometry: bed, thickness and surface +disp(' creating thickness'); +md.geometry.surface= +md.geometry.bed= +md.geometry.thickness= + +%Create drag parameters (won't be used as ice is frozen on the bed) +disp(' creating drag'); +md.friction.coefficient= +md.friction.p=ones(md.mesh.numberofvertices,1); +md.friction.q=ones(md.mesh.numberofvertices,1); + +%Create materials parameters: rheology parameters +disp(' creating flow law paramter'); +md.materials.rheology_B= +md.materials.rheology_n= + +%Create generic boundary conditions +disp(' boundary conditions for diagnostic model'); +... Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.m =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.m (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.m (revision 14145) @@ -0,0 +1,44 @@ +L=80000; % 5000,10000,20000,40000,80000,160000 + +%Create model +md=model; + +%Create mesh +nx=20; %number of nodes in x direction +ny=20; +md=squaremesh(model,L,L,nx,ny); + +%Define mask +md=setmask(md,'',''); %ice sheet test + +%Parameterize model +md=parameterize(md,'./ISMIPAback.par'); + +%Extrude to have a 3d model +md=extrude(md,5,1); + +%Choose flow approximation +md=setflowequation(md,'pattyn','all'); + +%Create dirichlet on the bed only +md.diagnostic.spcvx=NaN*ones(md.mesh.numberofvertices,1); +md.diagnostic.spcvy=NaN*ones(md.mesh.numberofvertices,1); +md.diagnostic.spcvz=NaN*ones(md.mesh.numberofvertices,1); +pos=find(md.mesh.vertexonbed); +md.diagnostic.spcvx(pos)=0; +md.diagnostic.spcvy(pos)=0; +md.diagnostic.spcvz(pos)=0; + +%Create MPCs to have periodic boundary conditions +posx=find(md.mesh.x==0); +posx2=find(md.mesh.x==max(md.mesh.x)); +posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); %Don't take the same grids two times +posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x)); +md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; + +%Set cluster and print messages +md.cluster=generic('name',oshostname(),'np',2); +md.verbose=verbose('convergence',true,'solution',true); + +%Solve +md=solve(md,DiagnosticSolutionEnum); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.par =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.par (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPAback.par (revision 14145) @@ -0,0 +1,20 @@ +%Create geometry: bed, thickness and surface +disp(' creating thickness'); +md.geometry.surface=-md.mesh.x*tan(0.5*pi/180); +md.geometry.bed=md.geometry.surface-1000+500*sin(md.mesh.x*2*pi/max(md.mesh.x)).*sin(md.mesh.y*2*pi/max(md.mesh.x)); +md.geometry.thickness=md.geometry.surface-md.geometry.bed; + +%Create drag parameters (won't be used as ice is frozen on the bed) +disp(' creating drag'); +md.friction.coefficient=200*ones(md.mesh.numberofvertices,1); %q=1. +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=ones(md.mesh.numberofelements,1); + +%Create materials parameters: rheology parameters +disp(' creating flow law paramter'); +md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1); +md.materials.rheology_n=3*ones(md.mesh.numberofelements,1); + +%Create generic boundary conditions +disp(' boundary conditions for diagnostic model'); +md=SetIceSheetBC(md); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.m =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.m (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.m (revision 14145) @@ -0,0 +1,41 @@ +%Create model +md=model; + +%Create mesh +... + +%Set mask +... + +%Parameterize +... + +%Extrude +... + +%Set flow equation +... + +%Boundary conditions +md.diagnostic.spcvx(:)=NaN; +md.diagnostic.spcvy(:)=NaN; +md.diagnostic.spcvz(:)=NaN; +%Create dirichlet on the bed if no slip +... + +%Create MPCs to have periodic boundary conditions +posx=find(md.mesh.x==0); +posx2=find(md.mesh.x==max(md.mesh.x)); +posy=find(md.mesh.y==0); %Don't take the same grids two times +posy2=find(md.mesh.y==max(md.mesh.y)); +md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; +md.prognostic.vertex_pairing=[posx,posx2;posy,posy2]; + +%Transient parameters in md.timestepping and md.transient + +%Compute the diagnostic +md.cluster=generic('name',oshostname(),'np',2); +md.verbose=verbose('convergence',true,'solution',true); + +%Solve +... Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.par =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.par (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPF.par (revision 14145) @@ -0,0 +1,22 @@ +disp(' creating geometry'); +md.geometry.surface= +md.geometry.bed= +md.geometry.thickness= + +disp(' creating drag'); +md.friction.coefficient= +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=zeros(md.mesh.numberofelements,1); + +disp(' creating flow law paramter'); +md.materials.rheology_B=(1/(2.140373*10^-7/(365.25*24*3600)))*ones(md.mesh.numberofvertices,1); +md.materials.rheology_n=1*ones(md.mesh.numberofelements,1); + +disp(' boundary conditions for diagnostic model'); +md=SetIceSheetBC(md); + +%Field initialization +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.m =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.m (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.m (revision 14145) @@ -0,0 +1,41 @@ +%Create model +md=model; + +%Create mesh +L=100000; %in m +nx=20; %numberof nodes in x direction +ny=20; +md=squaremesh(md,L,L,nx,ny); + +md=setmask(md,'',''); %ice sheet test +md=parameterize(md,'./ISMIPFback.par'); +md=extrude(md,5,1); +md=setflowequation(md,'pattyn','all'); + +%Boundary conditions +md.diagnostic.spcvx(:)=NaN; +md.diagnostic.spcvy(:)=NaN; +md.diagnostic.spcvz(:)=NaN; +%Create dirichlet on the bed if no slip +pos=find(md.mesh.vertexonbed); +md.diagnostic.spcvx(pos)=0; +md.diagnostic.spcvy(pos)=0; +md.diagnostic.spcvz(pos)=0; + +%Create MPCs to have periodic boundary conditions +posx=find(md.mesh.x==0); +posx2=find(md.mesh.x==max(md.mesh.x)); +posy=find(md.mesh.y==0); %Don't take the same grids two times +posy2=find(md.mesh.y==max(md.mesh.y)); +md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; +md.prognostic.vertex_pairing=[posx,posx2;posy,posy2]; + +%Transient parameters +md.timestepping.time_step=4; +md.timestepping.final_time=4*20; +md.transient.isthermal=0; + +%Compute the diagnostic +md.cluster=generic('name',oshostname(),'np',2); +md.verbose=verbose('convergence',true,'solution',true); +md=solve(md,TransientSolutionEnum); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.par =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.par (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/ISMIPFback.par (revision 14145) @@ -0,0 +1,24 @@ +disp(' creating geometry'); +L=100000; +alpha=3*pi/180; +md.geometry.surface=-md.mesh.x*tan(alpha); +md.geometry.bed=md.geometry.surface-1000+100*exp(-((md.mesh.x-L/2).^2+(md.mesh.y-L/2).^2)/(10000^2)); +md.geometry.thickness=md.geometry.surface-md.geometry.bed; + +disp(' creating drag'); +md.friction.coefficient=sqrt(365.25*24*3600/(1000*2.140373*10^-7))*ones(md.mesh.numberofvertices,1); +md.friction.p=ones(md.mesh.numberofelements,1); +md.friction.q=zeros(md.mesh.numberofelements,1); + +disp(' creating flow law paramter'); +md.materials.rheology_B=(1/(2.140373*10^-7/(365.25*24*3600)))*ones(md.mesh.numberofvertices,1); +md.materials.rheology_n=1*ones(md.mesh.numberofelements,1); + +disp(' boundary conditions for diagnostic model'); +md=SetIceSheetBC(md); + +%Field initialization +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.vz=zeros(md.mesh.numberofvertices,1); +md.initialization.pressure=zeros(md.mesh.numberofvertices,1); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip.par =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip.par (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip.par (revision 14145) @@ -0,0 +1,53 @@ +%Ok, start defining model parameters here + +disp(' creating geometry'); +md.geometry.bathymetry= +md.geometry.bed= +md.geometry.thickness= +md.geometry.surface= + +disp(' creating drag'); +md.friction.coefficient= +md.friction.p= +md.friction.q= + +disp(' creating flow law paramter'); +md.materials.rheology_B= +md.materials.rheology_n= +md.materials.rheology_law='None'; + +disp(' boundary conditions for diagnostic model'); +%Create node on boundary fist (because we cannot use mesh) +md=SetIceShelfBC(md,'MismipFront.exp'); +md.diagnostic.icefront(:,end)=1; +md.diagnostic.spcvx(:)=NaN; +md.diagnostic.spcvy(:)=NaN; +pos=find((md.mesh.y<50000.1 & md.mesh.y>49999.9) | (md.mesh.y<0.1 & md.mesh.y>-0.1)); +md.diagnostic.spcvy(pos)=0; +pos2=find(md.mesh.x<0.1 & md.mesh.x>-0.1); +md.diagnostic.spcvx(pos2)=0; +md.diagnostic.spcvz(pos2)=0; + +disp(' forcing conditions'); +md.surfaceforcings.mass_balance= +md.basalforcings.geothermalflux= + +md.thermal.spctemperature=NaN*ones(md.mesh.numberofvertices,1); +md.groundingline.migration='AgressiveMigration'; + +%Parameters +md.materials.rho_ice= +md.materials.rho_water= +md.constants.g= +md.constants.yts= +md.transient.isthermal= +md.transient.isgroundingline= +md.diagnostic.isnewton= + +%Initialization +md.initialization.vx=ones(md.mesh.numberofvertices,1); +md.initialization.vy=ones(md.mesh.numberofvertices,1); +md.initialization.vz=ones(md.mesh.numberofvertices,1); +md.initialization.vel=sqrt(2)*ones(md.mesh.numberofvertices,1); +md.initialization.pressure=md.constants.g*md.materials.rho_ice*md.geometry.thickness; +md.initialization.temperature=273*ones(md.mesh.numberofvertices,1); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip3d.m =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip3d.m (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/Mismip3d.m (revision 14145) @@ -0,0 +1,8 @@ +md=triangle(model,'MismipDomain.exp',10000); +md=setmask(md,'',''); +md=parameterize(md,'Mismip.par'); +md=setflowequation(md,'macayeal','all'); +md.transient.isgroundingline=1; +md.timestepping.time_step=100; +md.timestepping.final_time=1000; +md=solve(md,TransientSolutionEnum); Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipDomain.exp =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipDomain.exp (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipDomain.exp (revision 14145) @@ -0,0 +1,10 @@ +## Name:domainoutline +## Icon:0 +# Points Count Value +5 1. +# X pos Y pos +0 0 +800000 0 +800000 50000 +0 50000 +0 0 Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipFront.exp =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipFront.exp (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/MismipFront.exp (revision 14145) @@ -0,0 +1,10 @@ +## Name:front +## Icon:0 +# Points Count Value +5 1. +# X pos Y pos +800001 1 +799999 1 +799999 49999 +800001 49999 +800001 1 Index: /issm/workshop/2012/Talks/11_ISMIP/Codes/Profile.exp =================================================================== --- /issm/workshop/2012/Talks/11_ISMIP/Codes/Profile.exp (revision 14145) +++ /issm/workshop/2012/Talks/11_ISMIP/Codes/Profile.exp (revision 14145) @@ -0,0 +1,7 @@ +## Name:icefront +## Icon:0 +# Points Count Value +2 1. +# X pos Y pos +0 50000 +100000 50000