Index: /issm/trunk-jpl/jenkins/javascript/karma/karma.conf.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/karma.conf.js	(revision 20882)
+++ /issm/trunk-jpl/jenkins/javascript/karma/karma.conf.js	(revision 20883)
@@ -40,8 +40,11 @@
       'lib/bin/qmu.js',
       'lib/bin/basalforcings.js',
+      'lib/bin/DepthAverage.js',
+      'lib/bin/project2d.js',
       'lib/bin/project3d.js',
       'lib/bin/model.js',
       //'scripts/specs/issm.spec.js'
-      'scripts/specs/temp.spec.js'
+      //'scripts/specs/temp.spec.js'
+      'scripts/test102.js'
       //'scripts/specs/3.spec.js'
     ],
Index: /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/DepthAverage.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/DepthAverage.js	(revision 20883)
+++ /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/DepthAverage.js	(revision 20883)
@@ -0,0 +1,1 @@
+link /home/andy/Programming/Research/trunk-jpl/src/m/extrusion/DepthAverage.js
Index: /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/contourenvelope.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/contourenvelope.js	(revision 20883)
+++ /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/contourenvelope.js	(revision 20883)
@@ -0,0 +1,1 @@
+link /home/andy/Programming/Research/trunk-jpl/src/m/parameterization/contourenvelope.js
Index: /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/issm-binaries.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/issm-binaries.js	(revision 20882)
+++ /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/issm-binaries.js	(revision 20883)
@@ -11182,5 +11182,4 @@
 	//Marshall into a binary array (fid) all the fields of model.
 	var fid = marshall(md);                                          // bin file
-	
 	//deal with toolkits options: 
 	toolkitsstring= md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits'); // toolkits file
Index: /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/project2d.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/project2d.js	(revision 20883)
+++ /issm/trunk-jpl/jenkins/javascript/karma/lib/bin/project2d.js	(revision 20883)
@@ -0,0 +1,1 @@
+link /home/andy/Programming/Research/trunk-jpl/src/m/extrusion/project2d.js
Index: /issm/trunk-jpl/jenkins/javascript/karma/scripts/test104.js
===================================================================
--- /issm/trunk-jpl/jenkins/javascript/karma/scripts/test104.js	(revision 20882)
+++ /issm/trunk-jpl/jenkins/javascript/karma/scripts/test104.js	(revision 20883)
@@ -7,5 +7,5 @@
 setflowequation(md,'FS','all');
 //md.cluster=generic('name',oshostname(),'np',3);
-md=solve(md,StressbalanceSolutionEnum());
+md=solve(md,StressbalanceSolutionEnum(),'checkconsistency','no');
 
 //Fields and tolerances to track changes
Index: /issm/trunk-jpl/src/m/classes/model.js
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.js	(revision 20882)
+++ /issm/trunk-jpl/src/m/classes/model.js	(revision 20883)
@@ -154,5 +154,5 @@
                 numlayers=extrusionlist.length;
             } else if (argc==3) { //one polynomial law
-                if (arguments[2]<=0) {
+                if (arguments[1]<=0) {
                     //help extrude;
                     throw 'extrusionexponent must be >=0';
@@ -163,5 +163,5 @@
 
                 for (var i = 0; i < numlayers; i++) {
-                    extrusionlist.push(Math.pow(i/(numlayers-1), arguments[2]));
+                    extrusionlist.push(Math.pow(i/(numlayers-1), arguments[1]));
                 };
 
@@ -265,6 +265,6 @@
             //uppervertex=NaN*ones(number_nodes3d,1);
 
-            lowervertex = fillArray(null, number_nodes3d);
-            uppervertex = fillArray(null, number_nodes3d);
+            lowervertex = fillArray(NaN, number_nodes3d);
+            uppervertex = fillArray(NaN, number_nodes3d);
 
             //lowervertex(md.mesh.numberofvertices+1:end)=1:(numlayers-1)*md.mesh.numberofvertices;
@@ -283,6 +283,6 @@
 
             //same for lower and upper elements
-            lowerelements = fillArray(null, number_el3d);
-            upperelements = fillArray(null, number_el3d);
+            lowerelements = fillArray(NaN, number_el3d);
+            upperelements = fillArray(NaN, number_el3d);
             //lowerelements(md.mesh.numberofelements+1:end)=1:(numlayers-2)*md.mesh.numberofelements;
             //upperelements(1:(numlayers-2)*md.mesh.numberofelements)=md.mesh.numberofelements+1:(numlayers-1)*md.mesh.numberofelements;
@@ -326,62 +326,252 @@
             md.mesh.long=project3d(md,'vector',md.mesh.long,'type','node');
 
-                md.geometry=md.geometry.extrude(md);
-                md.friction=md.friction.extrude(md);
-                md.inversion=md.inversion.extrude(md);
-                md.smb=md.smb.extrude(md);
-                md.initialization=md.initialization.extrude(md);
-
-                md.flowequation=md.flowequation.extrude(md);
-                md.stressbalance=md.stressbalance.extrude(md);
-                md.thermal=md.thermal.extrude(md);
-                md.masstransport=md.masstransport.extrude(md);
-                md.levelset=md.levelset.extrude(md);
-                md.calving=md.calving.extrude(md);
-                md.hydrology = md.hydrology.extrude(md);
-
-                //connectivity
-                //if ~isnan(md.mesh.elementconnectivity)
-
-                if (md.mesh.elementconnectivity.every(function(e, i, arr) { return e !== null && typeof e !== 'undefined'; })) {
-                    var temparr = [];
-                    for (var i = 0; i < numlayers; i++) {
-                        temparr.push(md.mesh.elementconnectivity);
-                    };
-
-                    //md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
-
-                    md.mesh.elementconnectivity = temparr;
-
-                    //md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
-
-                    for (var i = 0; i < md.mesh.elementconnectivity.length; i++) {
-                        if (md.mesh.elementconnectivity[i] == 0) {
-                            md.mesh.elementconnectivity[i] = NaN;
-                        }
-                    };
-
-                    for (var i = 2; i < numlayers; i++) {
-                        for (var j = (i-1)*md.mesh.nuberofelements2d+1; i <= (i)*md.mesh.numberofelements2d; i++) {
-                            md.mesh.elementconnectivity[i] += md.mesh.numberofelements2d;
-                        };
-                        //md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
-                            //=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
-                    };
-
-                    md.mesh.elementconnectivity.map(function(x) { return (isNaN(x)) ? 0 : x; });
-                }
-
-                md.materials=md.materials.extrude(md);
-                md.damage=md.damage.extrude(md);
-                md.mask=md.mask.extrude(md);
-                md.qmu=md.qmu.extrude(md);
-                md.basalforcings=md.basalforcings.extrude(md);
-
-                //increase connectivity if less than 25:
-                if (md.mesh.average_vertex_connectivity<=25)
-                    md.mesh.average_vertex_connectivity=100;
-                // }}}
-
+                //md.geometry=md.geometry.extrude(md);
+                //md.friction=md.friction.extrude(md);
+                //md.inversion=md.inversion.extrude(md);
+                //md.smb=md.smb.extrude(md);
+                //md.initialization=md.initialization.extrude(md);
+
+                //md.flowequation=md.flowequation.extrude(md);
+                //md.stressbalance=md.stressbalance.extrude(md);
+                //md.thermal=md.thermal.extrude(md);
+                //md.masstransport=md.masstransport.extrude(md);
+                //md.levelset=md.levelset.extrude(md);
+                //md.calving=md.calving.extrude(md);
+                //md.hydrology = md.hydrology.extrude(md);
+
+                ////connectivity
+                ////if ~isnan(md.mesh.elementconnectivity)
+
+                //if (md.mesh.elementconnectivity.every(function(e, i, arr) { return e !== null && typeof e !== 'undefined'; })) {
+                    //var temparr = [];
+                    //for (var i = 0; i < numlayers; i++) {
+                        //temparr.push(md.mesh.elementconnectivity);
+                    //};
+
+                    ////md.mesh.elementconnectivity=repmat(md.mesh.elementconnectivity,numlayers-1,1);
+
+                    //md.mesh.elementconnectivity = temparr;
+
+                    ////md.mesh.elementconnectivity(find(md.mesh.elementconnectivity==0))=NaN;
+
+                    //for (var i = 0; i < md.mesh.elementconnectivity.length; i++) {
+                        //if (md.mesh.elementconnectivity[i] == 0) {
+                            //md.mesh.elementconnectivity[i] = NaN;
+                        //}
+                    //};
+
+                    //for (var i = 2; i < numlayers; i++) {
+                        //for (var j = (i-1)*md.mesh.nuberofelements2d+1; i <= (i)*md.mesh.numberofelements2d; i++) {
+                            //md.mesh.elementconnectivity[i] += md.mesh.numberofelements2d;
+                        //};
+                        ////md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)...
+                            ////=md.mesh.elementconnectivity((i-1)*md.mesh.numberofelements2d+1:(i)*md.mesh.numberofelements2d,:)+md.mesh.numberofelements2d;
+                    //};
+
+                    //md.mesh.elementconnectivity.map(function(x) { return (isNaN(x)) ? 0 : x; });
+                //}
+
+                //md.materials=md.materials.extrude(md);
+                //md.damage=md.damage.extrude(md);
+                //md.mask=md.mask.extrude(md);
+                //md.qmu=md.qmu.extrude(md);
+                //md.basalforcings=md.basalforcings.extrude(md);
+
+                ////increase connectivity if less than 25:
+                //if (md.mesh.average_vertex_connectivity<=25)
+                    //md.mesh.average_vertex_connectivity=100;
         } //}}}
+        this.collapse = function(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.elementtype() !== 'Penta') {
+                console.error('collapse error message: only 3d mesh can be collapsed')
+            }
+
+            // Start with changing all the fields from the 3d mesh 
+
+            // dealing with the friction law
+            // drag is limited to nodes that are on the bedrock.
+            if (md.friction.classname() === 'friction') {
+                md.friction.coefficient=project2d(md,md.friction.coefficient,1);
+                md.friction.p=project2d(md,md.friction.p,1);
+                md.friction.q=project2d(md,md.friction.q,1);
+            } else if (md.friction.classname() === 'frictionhydro') {
+                md.friction.q=project2d(md,md.friction.q,1);
+                md.friction.C=project2d(md,md.friction.C,1);
+                md.friction.As=project2d(md,md.friction.As,1);
+                md.friction.effective_pressure=project2d(md,md.friction.effective_pressure,1);
+            } else if (md.friction.classname() === 'frictionwaterlayer') {
+                md.friction.coefficient=project2d(md,md.friction.coefficient,1);
+                md.friction.p=project2d(md,md.friction.p,1);
+                md.friction.q=project2d(md,md.friction.q,1);
+                md.friction.water_layer=project2d(md,md.friction.water_layer,1);
+            } else if (md.friction.classname() === 'frictionweertman') {
+                md.friction.C=project2d(md,md.friction.C,1);
+                md.friction.m=project2d(md,md.friction.m,1);
+            } else if (md.friction.classname() === 'frictionweertmantemp') {
+                md.friction.C=project2d(md,md.friction.C,1);
+                md.friction.m=project2d(md,md.friction.m,1);
+            } else {
+                disp('friction type not supported');
+            }
+
+            // observations
+            if (!isNaN(md.inversion.vx_obs))
+                md.inversion.vx_obs=project2d(md,md.inversion.vx_obs,md.mesh.numberoflayers);
+
+            if (!isNaN(md.inversion.vy_obs))
+                md.inversion.vy_obs=project2d(md,md.inversion.vy_obs,md.mesh.numberoflayers);
+
+            if (!isNaN(md.inversion.vel_obs))
+                md.inversion.vel_obs=project2d(md,md.inversion.vel_obs,md.mesh.numberoflayers);
+
+            if (!isNaN(md.inversion.cost_functions_coefficients))
+                md.inversion.cost_functions_coefficients=project2d(md,md.inversion.cost_functions_coefficients,md.mesh.numberoflayers);
+
+            if (numel(md.inversion.min_parameters)>1)
+                md.inversion.min_parameters=project2d(md,md.inversion.min_parameters,md.mesh.numberoflayers);
+
+            if (numel(md.inversion.max_parameters)>1) 
+                md.inversion.max_parameters=project2d(md,md.inversion.max_parameters,md.mesh.numberoflayers);
+
+            if (md.smb.classname() === 'SMBforcing' && !isNaN(md.smb.mass_balance)) {
+                md.smb.mass_balance=project2d(md,md.smb.mass_balance,md.mesh.numberoflayers);
+            } else if (md.smb.classname() === 'SMBhenning' && !isNaN(md.smb.smbref)) {
+                md.smb.smbref=project2d(md,md.smb.smbref,md.mesh.numberoflayers);
+            }
+
+            // results
+            if (!isNaN(md.initialization.vx))
+                md.initialization.vx=DepthAverage(md,md.initialization.vx);
+            if (!isNaN(md.initialization.vy))
+                md.initialization.vy=DepthAverage(md,md.initialization.vy);
+            if (!isNaN(md.initialization.vz))
+                md.initialization.vz=DepthAverage(md,md.initialization.vz);
+            if (!isNaN(md.initialization.vel))
+                md.initialization.vel=DepthAverage(md,md.initialization.vel);
+            if (!isNaN(md.initialization.temperature))
+                md.initialization.temperature=DepthAverage(md,md.initialization.temperature);
+            if (!isNaN(md.initialization.pressure))
+                md.initialization.pressure=project2d(md,md.initialization.pressure,1);
+            if (!isNaN(md.initialization.sediment_head))
+                md.initialization.sediment_head=project2d(md,md.initialization.sediment_head,1);
+            if (!isNaN(md.initialization.epl_head))
+                md.initialization.epl_head=project2d(md,md.initialization.epl_head,1);
+            if (!isNaN(md.initialization.epl_thickness))
+                md.initialization.epl_thickness=project2d(md,md.initialization.epl_thickness,1);
+
+            // gia
+            if (!isNaN(md.gia.mantle_viscosity))
+                md.gia.mantle_viscosity=project2d(md,md.gia.mantle_viscosity,1);
+            if (!isNaN(md.gia.lithosphere_thickness))
+                md.gia.lithosphere_thickness=project2d(md,md.gia.lithosphere_thickness,1);
+
+            // elementstype
+            if (!isNaN(md.flowequation.element_equation)) {
+                md.flowequation.element_equation=project2d(md,md.flowequation.element_equation,1);
+                md.flowequation.vertex_equation=project2d(md,md.flowequation.vertex_equation,1);
+                md.flowequation.borderSSA=project2d(md,md.flowequation.borderSSA,1);
+                md.flowequation.borderHO=project2d(md,md.flowequation.borderHO,1);
+                md.flowequation.borderFS=project2d(md,md.flowequation.borderFS,1);
+            }
+
+            // boundary conditions
+            md.stressbalance.spcvx=project2d(md,md.stressbalance.spcvx,md.mesh.numberoflayers);
+            md.stressbalance.spcvy=project2d(md,md.stressbalance.spcvy,md.mesh.numberoflayers);
+            md.stressbalance.spcvz=project2d(md,md.stressbalance.spcvz,md.mesh.numberoflayers);
+            md.stressbalance.referential=project2d(md,md.stressbalance.referential,md.mesh.numberoflayers);
+            md.stressbalance.loadingforce=project2d(md,md.stressbalance.loadingforce,md.mesh.numberoflayers);
+            md.masstransport.spcthickness=project2d(md,md.masstransport.spcthickness,md.mesh.numberoflayers);
+            if (!isNaN(md.damage.spcdamage))
+                md.damage.spcdamage=project2d(md,md.damage.spcdamage,md.mesh.numberoflayers);
+            md.thermal.spctemperature=project2d(md,md.thermal.spctemperature,md.mesh.numberoflayers);
+
+            // Hydrologydc variables
+            if (md.hydrology.classname() === 'hydrologydc') {
+                md.hydrology.spcsediment_head=project2d(md,md.hydrology.spcsediment_head,1);
+                md.hydrology.mask_eplactive_node=project2d(md,md.hydrology.mask_eplactive_node,1);
+                md.hydrology.sediment_transmitivity=project2d(md,md.hydrology.sediment_transmitivity,1);
+                md.hydrology.basal_moulin_input=project2d(md,md.hydrology.basal_moulin_input,1);
+
+                if(md.hydrology.isefficientlayer==1)
+                    md.hydrology.spcepl_head=project2d(md,md.hydrology.spcepl_head,1);
+            }
+            
+            // materials
+            md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
+            md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
+            
+            // damage: 
+            if (md.damage.isdamage)
+                md.damage.D=DepthAverage(md,md.damage.D);
+
+            // special for thermal modeling:
+            if (!isNaN(md.basalforcings.groundedice_melting_rate))
+                md.basalforcings.groundedice_melting_rate=project2d(md,md.basalforcings.groundedice_melting_rate,1); 
+
+            if (!isNaN(md.basalforcings.floatingice_melting_rate))
+                md.basalforcings.floatingice_melting_rate=project2d(md,md.basalforcings.floatingice_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
+            var nodes2d=md.mesh.numberofvertices2d;
+            var 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.base=project2d(md,md.geometry.base,1);
+            if (!isNaN(md.geometry.bed))
+                md.geometry.bed=project2d(md,md.geometry.bed,1);
+
+            if (!isNaN(md.mask.groundedice_levelset))
+                md.mask.groundedice_levelset=project2d(md,md.mask.groundedice_levelset,1);
+
+            if (!isNaN(md.mask.ice_levelset))
+                md.mask.ice_levelset=project2d(md,md.mask.ice_levelset,1);
+
+            // lat long
+            if (numel(md.mesh.lat) === md.mesh.numberofvertices)
+                md.mesh.lat=project2d(md,md.mesh.lat,1);
+            if (numel(md.mesh.long) === md.mesh.numberofvertices)
+                md.mesh.long=project2d(md,md.mesh.long,1);
+
+            // Initialize with the 2d mesh
+            var mesh = new mesh2d();
+            mesh.x=md.mesh.x2d;
+            mesh.y=md.mesh.y2d;
+            mesh.numberofvertices=md.mesh.numberofvertices2d;
+            mesh.numberofelements=md.mesh.numberofelements2d;
+            mesh.elements=md.mesh.elements2d;
+
+            if (!isNaN(md.mesh.vertexonboundary))
+                mesh.vertexonboundary=project2d(md,md.mesh.vertexonboundary,1);
+            if (!isNaN(md.mesh.elementconnectivity))
+                mesh.elementconnectivity=project2d(md,md.mesh.elementconnectivity,1);
+
+            md.mesh=mesh;
+            md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
+            md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
+            md.mesh.segments=contourenvelope(md.mesh);
+
+            return md;
+        }
 	//properties
 	// {{{
Index: /issm/trunk-jpl/src/m/extrusion/DepthAverage.js
===================================================================
--- /issm/trunk-jpl/src/m/extrusion/DepthAverage.js	(revision 20883)
+++ /issm/trunk-jpl/src/m/extrusion/DepthAverage.js	(revision 20883)
@@ -0,0 +1,52 @@
+function DepthAverage(md,vector) {
+    // DEPTHAVERAGE - computes depth average of 3d vector using the trapezoidal rule, and returns the value on 2d mesh. 
+    // 
+    //    Usage:
+    //       vector_average=DepthAverage(md,vector);
+    // 
+    //    Example:
+    //       vel_bar=DepthAverage(md,md.initialization.vel);
+
+    // check that the model given in input is 3d
+    if (md.mesh.elementtype() !== 'Penta') {
+        console.error('DepthAverage error message: the model given in input must be 3d');
+    }
+
+    // nods data
+    if (vector.length === md.mesh.numberofvertices) {
+        var vector_average=zeros(md.mesh.numberofvertices2d,1);
+
+        for (var i = 1; i < md.mesh.numberoflayers-1; ++i) {
+            vector_average = vector_average.map(function(x) {
+                return x + (project2d(md, vector, i) + project2d(md,vector,i+1))/2;
+            }).map(function(y) {
+                return y * (project2d(md, md.mesh.z, i+1) - project2d(md, md.mesh.z, i));
+            });
+        }
+
+        vector_average = vector_average.map(function(z) {
+            return z / project2d(md, md.geometry.thickness, 1);
+        });
+
+        return vector_average;
+    }
+    // element data
+    else if (vector.length === md.mesh.numberofelements) {
+        var vector_average=zeros(md.mesh.numberofelements2d,1);
+        for (var i = 1; i < md.mesh.numberoflayers-1; ++i) {
+            vector_average = vector_average.map(function(x) {
+                return x + project2d(md, vector, i);
+            }).map(function(y) {
+                return y * (project2d(md, md.mesh.z, i+1) - project2d(md, md.mesh.z, i));
+            });
+        }
+
+        vector_average = vector_average.map(function(z) {
+            return z / project2d(md, md.geometry.thickness, 1);
+        });
+
+        return vector_average;
+    } else {
+        console.error('vector size not supported yet');
+    }
+}
Index: /issm/trunk-jpl/src/m/extrusion/project2d.js
===================================================================
--- /issm/trunk-jpl/src/m/extrusion/project2d.js	(revision 20883)
+++ /issm/trunk-jpl/src/m/extrusion/project2d.js	(revision 20883)
@@ -0,0 +1,46 @@
+function project2d(md3d,value,layer) {
+    // PROJECT2D - returns the value of a field for a given layer of the mesh
+    // 
+    //    returns the value of a vector for a given layer from extruded mesh onto the 2d mesh 
+    //    used to do the extrusion. This function is used to compare values between different
+    //    layers of a 3d mesh.
+    // 
+    //    Usage:
+    //       projection_value=project2d(md3d,value,layer)
+    // 
+    //    Example:
+    //       vel2=project2d(md3d,md3d.initialization.vel,2);
+    //       returns the velocity of the second layer (1 is the base)
+
+    // some checks on list of arguments
+    if (arguments.length !== 3) {
+        console.error('project2d error message');
+    }
+
+    if (md3d.mesh.domaintype() !== '3D') {
+        console.error("wrong model type ... should be ''3d''");
+    }
+
+    if (layer<1 || layer>md3d.mesh.numberoflayers) {
+        console.error(['layer must be between 1 and ' + num2str(md3d.mesh.numberoflayers)]);
+    }
+
+    // Return the projection value
+    var temp = [];
+    if (value.length === md3d.mesh.numberofvertices) {
+        for (var i = (layer-1)*md3d.mesh.numberofvertices2d; i <= layer*md3d.mesh.numberofvertices2d; ++i) {
+            temp.push(value[i]);
+        }
+    } else if (value.length === md3d.mesh.numberofvertices+1) {
+        for (var i = (layer-1)*md3d.mesh.numberofvertices2d; i <= layer*md3d.mesh.numberofvertices2d; ++i) {
+            temp.push(value[i]);
+        }
+        temp.push(value[value.length-1]);
+    } else {
+        for (var i = (layer-1)*md3d.mesh.numberofelements; i <= layer*md3d.mesh.numberofelements2d; ++i) {
+            temp.push(value[i]);
+        }
+    }
+
+    return temp;
+}
Index: /issm/trunk-jpl/src/m/mesh/meshconvert.js
===================================================================
--- /issm/trunk-jpl/src/m/mesh/meshconvert.js	(revision 20883)
+++ /issm/trunk-jpl/src/m/mesh/meshconvert.js	(revision 20883)
@@ -0,0 +1,49 @@
+this.meshconvert = function(md,varargin)) {
+function zeros(...args) {
+	var array = [];
+	for (var i = 0; i < args[0]; ++i) {
+		array.push(args.length == 1 ? 0 : zeros(args.slice(1)));
+	}
+	return array;
+}
+var md = new model();
+//CONVERTMESH - convert mesh to bamg mesh
+//
+//   Usage:
+//      md=meshconvert(md);
+//      md=meshconvert(md,index,x,y);
+
+if (nargin~=1 & nargin~=4,) {
+	help meshconvert
+	console.error('meshconvert error message: bad usage');
+}
+
+if (nargin==1,) {
+	index = md.mesh.elements;
+	x     = md.mesh.x;
+	y     = md.mesh.y;
+} else {
+	index = varargin[0];
+	x     = varargin[1];
+	y     = varargin[2];
+}
+
+//call Bamg
+[bamgmesh_out bamggeom_out]=BamgConvertMesh(index,x,y);
+
+// plug results onto model
+md.private.bamg          = struct();
+md.private.bamg.mesh     = bamgmesh(bamgmesh_out);
+md.private.bamg.geometry = bamggeom(bamggeom_out);
+md.mesh.x              = bamgmesh_out.Vertices(:,1);
+md.mesh.y              = bamgmesh_out.Vertices(:,2);
+md.mesh.elements       = bamgmesh_out.Triangles(:,1:3);
+md.mesh.edges          = bamgmesh_out.IssmEdges;
+md.mesh.segments       = bamgmesh_out.IssmSegments(:,1:3);
+md.mesh.segmentmarkers = bamgmesh_out.IssmSegments(:,4);
+
+//Fill in rest of fields:
+md.mesh.numberofelements = size(md.mesh.elements,1);
+md.mesh.numberofvertices = length(md.mesh.x);
+md.mesh.numberofedges    = size(md.mesh.edges,1);
+md.mesh.vertexonboundary = zeros(md.mesh.numberofvertices,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2)) = 1;