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 ~strcmpi(md.type,'3d'),
	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 grids that are on the bedrock.
md.drag=project2d(md,md.drag,1);

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

%observations
md.vx_obs=project2d(md,md.vx_obs,md.numlayers);
md.vy_obs=project2d(md,md.vy_obs,md.numlayers);
md.vel_obs=project2d(md,md.vel_obs,md.numlayers);
md.accumulation=project2d(md,md.accumulation,md.numlayers);
md.firn_layer=project2d(md,md.firn_layer,md.numlayers);

%results
if ~isnan(md.vx),md.vx=DepthAverage(md,md.vx);end;
if ~isnan(md.vy),md.vy=DepthAverage(md,md.vy);end;
if ~isnan(md.vz),md.vz=DepthAverage(md,md.vz);end;
if ~isnan(md.vel),md.vel=DepthAverage(md,md.vel);end;
if ~isnan(md.surface_slopex),md.surface_slopex=project2d(md,md.surface_slopex,md.numlayers);end;
if ~isnan(md.surface_slopey),md.surface_slopey=project2d(md,md.surface_slopey,md.numlayers);end;
if ~isnan(md.bed_slopex),md.bed_slopex=project2d(md,md.bed_slopex,1);end;
if ~isnan(md.bed_slopey),md.bed_slopey=project2d(md,md.bed_slopey,1);end;

%bedinfo and surface info
md.elementonbed=ones(md.numberofelements2d,1);
md.elementonsurface=ones(md.numberofelements2d,1);
md.gridonbed=ones(md.numberofgrids2d,1);
md.gridonsurface=ones(md.numberofgrids2d,1);

%elementstype
if ~isnan(md.elements_type2d)
	md.elements_type=md.elements_type2d; 
else
	md.elements_type2d(:,1)=project2d(md,md.elements_type(:,1),1);
	md.elements_type2d(:,2)=project2d(md,md.elements_type(:,2),1);
end	
md.gridonhutter=project2d(md,md.gridonhutter,1);
md.gridonmacayeal=project2d(md,md.gridonmacayeal,1);
md.gridonpattyn=project2d(md,md.gridonpattyn,1);
md.gridonstokes=project2d(md,md.gridonstokes,1);

%boundary conditions
md.gridondirichlet_diag=project2d(md,md.gridondirichlet_diag,md.numlayers);
dirichletvalues_diag2d(:,1)=project2d(md,md.dirichletvalues_diag(:,1),md.numlayers); 
dirichletvalues_diag2d(:,2)=project2d(md,md.dirichletvalues_diag(:,2),md.numlayers); 
md.dirichletvalues_diag=dirichletvalues_diag2d;

%Extrusion of Neumann BC
%in 2d, segmentonnumann is: [grid1 grid2 element]
numberofneumann2d=size(md.segmentonneumann_diag,1)/md.numlayers;
md.segmentonneumann_diag=[md.segmentonneumann_diag(1:numberofneumann2d,1:2) md.segmentonneumann_diag(1:numberofneumann2d,5)]; %Add two columns on the first layer 

%Prognostic
md.gridondirichlet_prog=project2d(md,md.gridondirichlet_prog,md.numlayers);
md.dirichletvalues_prog=project2d(md,md.dirichletvalues_prog,md.numlayers);
%md.segmentonneumann_prog=[tproj(md.segmentonneumann_prog(:,1)) tproj(md.segmentonneumann_prog(:,2)) tproj2d_el(md.segmentonneumann_prog(:,5))];
%md.segmentonneumann_prog2=[tproj(md.segmentonneumann_prog2(:,1)) tproj(md.segmentonneumann_prog2(:,2)) tproj2d_el(md.segmentonneumann_prog2(:,5))];

%materials
md.B=DepthAverage(md,md.B);
md.n=project2d(md,md.n,1);

%special for thermal modeling:
md.melting=project2d(md,md.melting,1); 
md.observed_temperature=DepthAverage(md,md.observed_temperature); 
md.geothermalflux=project2d(md,md.geothermalflux,1); %bedrock only gets geothermal flux
md.gridondirichlet_thermal=project2d(md,md.gridondirichlet_thermal,md.numlayers); %surface temperature
md.dirichletvalues_thermal=project2d(md,md.dirichletvalues_thermal,md.numlayers); %surface temperature

%NaN the values that are not on an spc'd temperature.
pos=find(~md.gridondirichlet_thermal);
md.dirichletvalues_thermal(pos)=NaN;

%update of connectivity matrix
md.connectivity=25;

%Collapse the mesh
grids2d=md.numberofgrids2d;
elements2d=md.numberofelements2d;

%parameters
md.surface=project2d(md,md.surface,1);
md.thickness=project2d(md,md.thickness,1);
md.bed=project2d(md,md.bed,1);
md.gridonboundary=project2d(md,md.gridonboundary,1);
md.elementoniceshelf=project2d(md,md.elementoniceshelf,1);
md.gridoniceshelf=project2d(md,md.gridoniceshelf,1);
md.elementonicesheet=project2d(md,md.elementonicesheet,1);
md.gridonicesheet=project2d(md,md.gridonicesheet,1);

%Initialize with the 2d mesh
md.x=md.x2d;
md.y=md.y2d;
md.z=md.z2d;
md.numberofgrids=md.numberofgrids2d;
md.numberofelements=md.numberofelements2d;
md.elements=md.elements2d;

%Keep a trace of lower and upper grids
md.lowergrids=NaN;
md.uppergrids=NaN;

%Remove old mesh 
md.x2d=NaN;
md.y2d=NaN;
md.z2d=NaN;
md.elements2d=NaN;
md.elements_type2d=md.elements_type;
md.numberofelements2d=md.numberofelements;
md.numberofgrids2d=md.numberofgrids;
md.numlayers=0;

%Update mesh type
md.type='2d';

%Tie in with Cielo if requested
cieloend;
end
