function md=collapse(md)
%COLLAPSE - collapses a 3d mesh into a 2d mesh
%
%   This routine collapses a 3d model into a 2d model
%   and collapses all the fileds of the 3d model by
%   taking their depth-averaged values
%
%   Usage:
%      md=collapse(md)
%
%   See also: EXTRUDE, MODELEXTRACT

%Check that the model is really a 3d model
if ~md.mesh.dimension==3,
	error('collapse error message: only 3d mesh can be collapsed')
end

%Start with changing alle the fields from the 3d mesh 

%drag is limited to nodes that are on the bedrock.
md.friction.coefficient=project2d(md,md.friction.coefficient,1);

%p and q (same deal, except for element that are on the bedrock: )
md.friction.p=project2d(md,md.friction.p,1);
md.friction.q=project2d(md,md.friction.q,1);

%observations
if ~isnan(md.inversion.vx_obs), md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers); end;
if ~isnan(md.inversion.vy_obs), md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers); end;
if ~isnan(md.inversion.vel_obs), md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers); end;
if ~isnan(md.surfaceforcings.mass_balance),
	md.surfaceforcings.mass_balance=project2d(md,md.surfaceforcings.mass_balance,md.mesh.numberoflayers); 
end;
if ~isnan(md.balancethickness.thickening_rate), md.balancethickness.thickening_rate=project2d(md,md.balancethickness.thickening_rate,md.mesh.numberoflayers); end;

%results
if ~isnan(md.initialization.vx),md.initialization.vx=DepthAverage(md,md.initialization.vx);end;
if ~isnan(md.initialization.vy),md.initialization.vy=DepthAverage(md,md.initialization.vy);end;
if ~isnan(md.initialization.vz),md.initialization.vz=DepthAverage(md,md.initialization.vz);end;
if ~isnan(md.initialization.vel),md.initialization.vel=DepthAverage(md,md.initialization.vel);end;

%bedinfo and surface info
md.mesh.elementonbed=ones(md.mesh.numberofelements2d,1);
md.mesh.elementonsurface=ones(md.mesh.numberofelements2d,1);
md.mesh.vertexonbed=ones(md.mesh.numberofvertices2d,1);
md.mesh.vertexonsurface=ones(md.mesh.numberofvertices2d,1);

%elementstype
if ~isnan(md.flowequation.element_equation)
	md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
end	

%boundary conditions
md.diagnostic.spcvx=project2d(md,md.diagnostic.spcvx,md.mesh.numberoflayers);
md.diagnostic.spcvy=project2d(md,md.diagnostic.spcvy,md.mesh.numberoflayers);
md.diagnostic.spcvz=project2d(md,md.diagnostic.spcvz,md.mesh.numberoflayers);
md.prognostic.spcthickness=project2d(md,md.prognostic.spcthickness,md.mesh.numberoflayers);
md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);

%Extrusion of Neumann BC
if ~isnan(md.diagnostic.icefront),
	numberofneumann2d=size(md.diagnostic.icefront,1)/md.mesh.numberoflayers;
	md.diagnostic.icefront=[md.diagnostic.icefront(1:numberofneumann2d,1:2) md.diagnostic.icefront(1:numberofneumann2d,5:6)]; %Add two columns on the first layer 
end

%materials
md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);

%special for thermal modeling:
md.basalforcings.melting_rate=project2d(md,md.basalforcings.melting_rate,1); 
md.basalforcings.geothermalflux=project2d(md,md.basalforcings.geothermalflux,1); %bedrock only gets geothermal flux

%update of connectivity matrix
md.mesh.average_vertex_connectivity=25;

%Collapse the mesh
nodes2d=md.mesh.numberofvertices2d;
elements2d=md.mesh.numberofelements2d;

%parameters
md.geometry.surface=project2d(md,md.geometry.surface,1);
md.geometry.thickness=project2d(md,md.geometry.thickness,1);
md.geometry.bed=project2d(md,md.geometry.bed,1);
md.mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
md.mask.elementonfloatingice=project2d(md,md.mask.elementonfloatingice,1);
md.mask.vertexonfloatingice=project2d(md,md.mask.vertexonfloatingice,1);
md.mask.elementongroundedice=project2d(md,md.mask.elementongroundedice,1);
md.mask.vertexongroundedice=project2d(md,md.mask.vertexongroundedice,1);

%Initialize with the 2d mesh
md.x=md.x2d;
md.y=md.y2d;
md.z=zeros(size(md.x2d));
md.mesh.numberofvertices=md.mesh.numberofvertices2d;
md.mesh.numberofelements=md.mesh.numberofelements2d;
md.elements=md.elements2d;

%Keep a trace of lower and upper nodes
md.mesh.lowervertex=NaN;
md.mesh.uppervertex=NaN;

%Remove old mesh 
md.x2d=NaN;
md.y2d=NaN;
md.elements2d=NaN;
md.mesh.numberofelements2d=md.mesh.numberofelements;
md.mesh.numberofvertices2d=md.mesh.numberofvertices;
md.mesh.numberoflayers=0;

%Update mesh type
md.mesh.dimension=2;