Changeset 11053

Timestamp:

12/13/11 15:48:20 (13 years ago)

Author:

schlegel

Message:

give them the answers

File:

• issm/workshop/17_SeaRISE/runme.m (modified) (9 diffs)

issm/workshop/17_SeaRISE/runme.m ¶

@@ -1 @@
-steps=[1];
+steps=[1:5];
 
 %Hard coded parameters
 cluster=generic('name',oshostname,'np',2);
-%What directory is you nc file in?
-modeldatapath=['./'];
-
+modeldatapath=['/Users/issm/Desktop/issm/workshop/14_JKS/Greenland5km_v1.2'];
 
 %Run Steps
+%Mesh{{{1
 if sum(find(steps==1))
 
@@ -22 +21 @@
 
         ncdata=[modeldatapath '/Greenland_5km_dev1.2.nc'];
+        %Use ncdisp('file') to see an ncdump
+        x1=double(ncread(ncdata,'x1'));
+        y1=double(ncread(ncdata,'y1'));
 
-        %Get x1 and y1 from ncdata
-        %Use ncdisp('file') to see an ncdump
+        velx=double(ncread(ncdata,'surfvelx'));
+        vely=double(ncread(ncdata,'surfvely'));
 
+        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);
 
-
-
-        %Get velx and vely from ncdata
-
-
-
-
-
-        %Interpolate grid to mesh velx and vely into md.inversion vx_obs, vy_obs, and
-        %calculate vel_obs
-
-
-
-
-        %Read thickness and Interpolate into md.geometry.thickness
-
-
-
+        thickness=double(ncread(ncdata,'thk'));
+        md.geometry.thickness=InterpFromGridToMesh(x1,y1,thickness',...
+                md.mesh.x,md.mesh.y,0);
 
         disp('   Mesh adaptation'); 
@@ -54 +44 @@
         name='./Models/Greenland.Mesh_generation';
         save(name,'md');
-end
+end% }}}
+%Parameterization{{{1
 if sum(find(steps==2)) 
-
-        %Set up your Parameter file.  Start by copying your Jks.par
 
         md=loadmodel('./Models/Greenland.Mesh_generation');
@@ -66 +55 @@
         name='./Models/Greenland.Parameterization'; 
         save(name,'md');
-end
+end%}}}
+%Control_drag{{{1
 if sum(find(steps==3)) 
 
@@ -79 +69 @@
         md.inversion.gradient_scaling(:)=100;
 
-
-        %Set up inversion control parameters
-
-
-
-
-
-
+        md.inversion.cost_functions=[];
+        md.inversion.cost_functions(1:floor(md.inversion.nsteps),1)=103;
+        md.inversion.cost_functions=[103*ones(md.inversion.nsteps,1)  ...
+                501*ones(md.inversion.nsteps,1)];
+        md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2);
+        md.inversion.cost_functions_coefficients(:,2)=8*10^-15;
+        md.inversion.control_parameters={'FrictionCoefficient'};
+        md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
 
         md=setcluster(md,cluster);
@@ -93 +84 @@
         md=solve(md,DiagnosticSolutionEnum);
 
-        %Save results in md.friction.coefficient, md.initialization.vx and md.initialization.vy
-        %before save
-
-
-
+        md.friction.coefficient=md.results.DiagnosticSolution.FrictionCoefficient;
+        md.initialization.vx=md.results.DiagnosticSolution.Vx;
+        md.initialization.vy=md.results.DiagnosticSolution.Vy;
 
         name='./Models/Greenland.Control_drag'; 
         save(name,'md');
-end
+end%}}}
+%Transient{{{1
 if sum(find(steps==4))
 
@@ -113 +103 @@
         y1=double(ncread(ncdata,'y1'));
 
-
-        %Get the surface mass balance from your nc file (Hint, your par file
-        %can help).  Interpolate and then Make sure you convert to m ice/year.
-
-
-        
-        %Set up forcing (forcing time in final row)
-
-
+        smb=double(ncread(ncdata,'smb'));
+        smb=InterpFromGridToMesh(x1,y1,smb',md.mesh.x,md.mesh.y,0)*...
+                1000/md.materials.rho_ice;
+        smb=[smb smb smb-1.0];
+        md.surfaceforcings.mass_balance=[smb; [1 10 20]];
 
         md.thermal.spctemperature=[md.initialization.temperature;1]; 
 
-        %Step of 0.2, for 20 years, output every year
-
-
-
+        md.timestepping.time_step=0.2;
+        md.timestepping.final_time=20;
+        md.settings.output_frequency=1;
 
         md.inversion.iscontrol=0;
@@ -140 +125 @@
         name='./Models/Greenland.Transient'; 
         save(name,'md');
-end% 
+end% }}}
+%Plot{{{1
 if sum(find(steps==5))
 
@@ -165 +151 @@
 
         %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');
 
         %Plot Volume
-
-
-
-end% 
+        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');
+        xlabel('years')
+end% }}}