Index: /issm/trunk/configs/config-greenplanet.sh
===================================================================
--- /issm/trunk/configs/config-greenplanet.sh	(revision 14066)
+++ /issm/trunk/configs/config-greenplanet.sh	(revision 14067)
@@ -22,4 +22,4 @@
  --with-blacs-dir=$ISSM_DIR/externalpackages/petsc/install/ \
  --with-graphics-lib=/usr/lib64/libX11.so \
- --with-cxxoptflags="-O3 -xS" \
+ --with-cxxoptflags="-O3" \
  --with-vendor=intel-linux
Index: /issm/trunk/configure.ac
===================================================================
--- /issm/trunk/configure.ac	(revision 14066)
+++ /issm/trunk/configure.ac	(revision 14067)
@@ -2,5 +2,5 @@
 
 #AUTOCONF
-AC_INIT([ISSM],[4.2.3],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
+AC_INIT([ISSM],[4.2.4],[issm@jpl.nasa.gov],[issm],[http://issm.jpl.nasa.gov]) #Initializing configure
 AC_CONFIG_AUX_DIR([./aux-config])         #Put config files in aux-config
 AC_CONFIG_MACRO_DIR([m4])                 #m4 macros are located in m4
Index: /issm/trunk/etc/environment.sh
===================================================================
--- /issm/trunk/etc/environment.sh	(revision 14066)
+++ /issm/trunk/etc/environment.sh	(revision 14067)
@@ -86,5 +86,5 @@
 
 DOXYGEN_DIR="$ISSM_DIR/externalpackages/doxygen/install"
-pathappend "$DOXYGEN_DIR/bin"
+pathprepend "$DOXYGEN_DIR/bin"
 
 AUTOTOOLS_DIR="$ISSM_DIR/externalpackages/autotools/install"
Index: /issm/trunk/externalpackages/doxygen/install.sh
===================================================================
--- /issm/trunk/externalpackages/doxygen/install.sh	(revision 14066)
+++ /issm/trunk/externalpackages/doxygen/install.sh	(revision 14067)
@@ -3,8 +3,8 @@
 
 #Some cleanup
-rm -rf install
+rm -rf install src
 
 #Download latest version
-svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk install
+svn co https://doxygen.svn.sourceforge.net/svnroot/doxygen/trunk src
 
 #Configure doxygen
Index: /issm/trunk/externalpackages/gsl/Makefile.am.patch
===================================================================
--- /issm/trunk/externalpackages/gsl/Makefile.am.patch	(revision 14067)
+++ /issm/trunk/externalpackages/gsl/Makefile.am.patch	(revision 14067)
@@ -0,0 +1,27 @@
+--- Makefile.am	2011-04-14 08:13:48.000000000 -0700
++++ Makefile.am_tmp	2012-11-12 14:54:02.000000000 -0800
+@@ -24,10 +24,22 @@
+ noinst_HEADERS = templates_on.h templates_off.h build.h
+ 
+ MINGW32_HOST = @MINGW32_HOST@
+-if MINGW32_HOST
++
++# Origional 'Makefile.am' sets 'libgslcblas.la' as a requirement for
++#'libgsl.ls' only if the host system is detected to be MingW32. This
++# is unfortunate as 'libgsl.la' has undefined 'cblas' symbols the
++# result is linking errors at run time. This patch sets the 'cblas'
++# library as a requirement for linking 'libgsl'.
++#
++# Origional script:
++#
++# if MINGW32_HOST
++# libgsl_la_LIBADD += cblas/libgslcblas.la
++# libgsl_la_LDFLAGS += -no-undefined
++# endif
++
+ libgsl_la_LIBADD += cblas/libgslcblas.la
+ libgsl_la_LDFLAGS += -no-undefined
+-endif
+ 
+ m4datadir = $(datadir)/aclocal
+ m4data_DATA = gsl.m4
Index: /issm/trunk/externalpackages/gsl/install-android.sh
===================================================================
--- /issm/trunk/externalpackages/gsl/install-android.sh	(revision 14066)
+++ /issm/trunk/externalpackages/gsl/install-android.sh	(revision 14067)
@@ -25,6 +25,7 @@
     cd src
 
+    mv ./../Makefile.am.patch ./
     patch Makefile.am < Makefile.am.patch
-    
+
     autoreconf -iv --force -I $ISSM_DIR/externalpackages/autotools/install/share/aclocal
 
@@ -38,7 +39,8 @@
 #Compile gsl
 if [[ $step == "3" || $step == "0" ]]; then
-	cd src
+	cd $ISSM_DIR/externalpackages/gsl/src
+
     if [ $# -eq 0 ]; then
-	    make $j 
+	    make 
     else
 	    make -j $j 
Index: /issm/trunk/externalpackages/triangle/configs/android/configure.make
===================================================================
--- /issm/trunk/externalpackages/triangle/configs/android/configure.make	(revision 14066)
+++ /issm/trunk/externalpackages/triangle/configs/android/configure.make	(revision 14067)
@@ -10,9 +10,7 @@
 # http://www.codesourcery.com/gnu_toolchains/arm/arm_gnu_linux_abi.pdf
 
-source $ANDROID_DIR/android_aux.sh
-
-CC=${toolchain_path}-gcc
-AR=${toolchain_path}-ar
-RANLIB=${toolchain_path}-ranlib
+CC=${host_triplet}-gcc
+AR=${host_triplet}-ar
+RANLIB=${host_triplet}-ranlib
 CSWITCHES = $(CFLAGS)
 TRILIBDEFS = -DTRILIBRARY
Index: /issm/trunk/externalpackages/triangle/install-android.sh
===================================================================
--- /issm/trunk/externalpackages/triangle/install-android.sh	(revision 14066)
+++ /issm/trunk/externalpackages/triangle/install-android.sh	(revision 14067)
@@ -1,4 +1,6 @@
 #!/bin/bash
 set -eu
+
+source $ANDROID_DIR/android_aux.sh
 
 #Some cleanup 
Index: /issm/trunk/src/android/ISSM/.classpath
===================================================================
--- /issm/trunk/src/android/ISSM/.classpath	(revision 14066)
+++ /issm/trunk/src/android/ISSM/.classpath	(revision 14067)
@@ -1,8 +1,8 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <classpath>
+	<classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>
+	<classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>
 	<classpathentry kind="src" path="src"/>
 	<classpathentry kind="src" path="gen"/>
-	<classpathentry kind="con" path="com.android.ide.eclipse.adt.ANDROID_FRAMEWORK"/>
-	<classpathentry kind="con" path="com.android.ide.eclipse.adt.LIBRARIES"/>
 	<classpathentry kind="output" path="bin/classes"/>
 </classpath>
Index: /issm/trunk/src/android/ISSM/AndroidManifest.xml
===================================================================
--- /issm/trunk/src/android/ISSM/AndroidManifest.xml	(revision 14066)
+++ /issm/trunk/src/android/ISSM/AndroidManifest.xml	(revision 14067)
@@ -5,19 +5,21 @@
 
     <uses-sdk
-        android:minSdkVersion="15"
+        android:minSdkVersion="10"
         android:targetSdkVersion="15" />
-
+	<uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
     <application
         android:icon="@drawable/ic_launcher"
         android:label="@string/app_name"
         android:theme="@style/AppTheme" >
-        <activity
+        <activity android:name=".MapSelection"
+            	  android:label="@string/title_activity_issm" >
+            <intent-filter>
+                <action android:name="android.intent.action.MAIN" />
+				<category android:name="android.intent.category.LAUNCHER" />
+            </intent-filter>
+            </activity>
+            <activity
             android:name=".ISSM"
             android:label="@string/title_activity_issm" >
-            <intent-filter>
-                <action android:name="android.intent.action.MAIN" />
-
-                <category android:name="android.intent.category.LAUNCHER" />
-            </intent-filter>
         </activity>
     </application>
Index: /issm/trunk/src/android/ISSM/bin/AndroidManifest.xml
===================================================================
--- /issm/trunk/src/android/ISSM/bin/AndroidManifest.xml	(revision 14066)
+++ /issm/trunk/src/android/ISSM/bin/AndroidManifest.xml	(revision 14067)
@@ -5,19 +5,21 @@
 
     <uses-sdk
-        android:minSdkVersion="15"
+        android:minSdkVersion="10"
         android:targetSdkVersion="15" />
-
+	<uses-permission android:name="android.permission.WRITE_EXTERNAL_STORAGE" />
     <application
         android:icon="@drawable/ic_launcher"
         android:label="@string/app_name"
         android:theme="@style/AppTheme" >
-        <activity
+        <activity android:name=".MapSelection"
+            	  android:label="@string/title_activity_issm" >
+            <intent-filter>
+                <action android:name="android.intent.action.MAIN" />
+				<category android:name="android.intent.category.LAUNCHER" />
+            </intent-filter>
+            </activity>
+            <activity
             android:name=".ISSM"
             android:label="@string/title_activity_issm" >
-            <intent-filter>
-                <action android:name="android.intent.action.MAIN" />
-
-                <category android:name="android.intent.category.LAUNCHER" />
-            </intent-filter>
         </activity>
     </application>
Index: /issm/trunk/src/android/ISSM/gen/com/example/issm/R.java
===================================================================
--- /issm/trunk/src/android/ISSM/gen/com/example/issm/R.java	(revision 14066)
+++ /issm/trunk/src/android/ISSM/gen/com/example/issm/R.java	(revision 14067)
@@ -22,10 +22,12 @@
     public static final class id {
         public static final int button1=0x7f080001;
+        public static final int button2=0x7f080003;
         public static final int input=0x7f080000;
-        public static final int menu_settings=0x7f080003;
+        public static final int menu_settings=0x7f080004;
         public static final int output=0x7f080002;
     }
     public static final class layout {
         public static final int activity_issm=0x7f030000;
+        public static final int activity_mapselection=0x7f030001;
     }
     public static final class menu {
Index: /issm/trunk/src/android/ISSM/jni/Android.mk
===================================================================
--- /issm/trunk/src/android/ISSM/jni/Android.mk	(revision 14066)
+++ /issm/trunk/src/android/ISSM/jni/Android.mk	(revision 14067)
@@ -11,5 +11,5 @@
 LOCAL_PATH := $(my_LOCAL_PATH)
 LOCAL_ALLOW_UNDEFINED_SYMBOLS := true
-LOCAL_CFLAGS := -Wno-psabi
+LOCAL_CFLAGS := -Wno-psabi -DHAVE_CONFIG_H
 LOCAL_LDLIBS := -llog -ldl
 LOCAL_MODULE := IssmJni
Index: /issm/trunk/src/android/ISSM/jni/Application.mk
===================================================================
--- /issm/trunk/src/android/ISSM/jni/Application.mk	(revision 14067)
+++ /issm/trunk/src/android/ISSM/jni/Application.mk	(revision 14067)
@@ -0,0 +1,3 @@
+APP_STL:=stlport_static
+APP_CPPFLAGS := -frtti
+APP_CPPFLAGS += -fexceptions
Index: /issm/trunk/src/android/ISSM/jni/Main.cpp
===================================================================
--- /issm/trunk/src/android/ISSM/jni/Main.cpp	(revision 14066)
+++ /issm/trunk/src/android/ISSM/jni/Main.cpp	(revision 14067)
@@ -1,19 +1,55 @@
 #include <jni.h>
 #include "../../../c/android/fac.h"
+#include "../../../c/issm.h"
+#include <cstddef>
+#include <stdio.h>
+///////////////////////////////////////////////////////////////////////////////////////////
 namespace com_example_issm
 {
-	static jlong factorial(JNIEnv *env, jclass clazz, jlong n)
+	fac* f;
+	FemModel *fm;
+//------------------------------------------------------------------------------------
+	jint initilize(JNIEnv *env, jclass clazz, jstring file)
 	{
-		fac *f = new fac();
-		jlong result = (jlong) (f->factorial(n));
-		delete(f);
-		return (jlong) result;
+		char *nativefile = (char*)env->GetStringUTFChars(file,0);
+
+		f = new fac();
+
+		//call Model constructor passing in infile as File Descriptor parameter.
+		fm  = new FemModel(nativefile);
+
+		env->ReleaseStringUTFChars(file, nativefile); //must realease the char*
+		//jint size = (jint) fm->bufferSize();
+		//return size;
+
+		return 0; //return 0 for now.
 	}
 //------------------------------------------------------------------------------------
-	static JNINativeMethod method_table[] = {
-			{"fac" ,"(J)J" , (void *) factorial},
+	//fill out the first two slots, extract the same way in java using get(int position)
+	void solve(JNIEnv *env, jclass clazz , jint alpha, jobject buf)
+	{
+		jdouble *dBuf = (jdouble *)env->GetDirectBufferAddress(buf);
+
+		//pass bBuff to fem model to allocate data
+		// fm -> solve(dBuf, alpha);
+
+	}
+//------------------------------------------------------------------------------------
+	jlong factorial(JNIEnv *env, jclass clazz, jlong n)
+	{
+		if( f != NULL)
+			return (jlong) (f->factorial(n));
+		return 0;
+	}
+
+//------------------------------------------------------------------------------------
+	static JNINativeMethod method_table[] =
+	{
+			{"fac"      	,     "(J)J" 	, (void *) factorial},
+			{"createISSMModel"   ,"(Ljava/lang/String;)I"  , (void *) initilize},
+			{"processBuffer", "(ILjava/nio/DoubleBuffer;)V", (void *) solve}
 	};
 }
-
+//------------------------------------------------------------------------------------
 using namespace com_example_issm;
 
@@ -29,5 +65,5 @@
     	if(clazz)
     	{
-    		env->RegisterNatives(clazz, method_table, 1);
+    		env->RegisterNatives(clazz, method_table, 3);
     		env->DeleteLocalRef(clazz);
     		return JNI_VERSION_1_6;
@@ -35,6 +71,4 @@
     	else return -1;
     }
-
-    // Get jclass with env->FindClass.
-    // Register methods with env->RegisterNatives.
 }
+///////////////////////////////////////////////////////////////////////////////////////////
Index: /issm/trunk/src/android/ISSM/jni/issmlib/Android.mk
===================================================================
--- /issm/trunk/src/android/ISSM/jni/issmlib/Android.mk	(revision 14066)
+++ /issm/trunk/src/android/ISSM/jni/issmlib/Android.mk	(revision 14067)
@@ -2,5 +2,5 @@
 include $(CLEAR_VARS)
 LOCAL_MODULE    := libISSMCore
-LOCAL_SRC_FILES := ../../../../c/libISSMCore.a
-LOCAL_EXPORT_C_INCLUDES := ../../../../c/android/
+LOCAL_SRC_FILES := libISSMCore.a
+LOCAL_EXPORT_C_INCLUDES := $(ISSM_DIR)
 include $(PREBUILT_STATIC_LIBRARY)
Index: /issm/trunk/src/android/ISSM/res/layout/activity_mapselection.xml
===================================================================
--- /issm/trunk/src/android/ISSM/res/layout/activity_mapselection.xml	(revision 14067)
+++ /issm/trunk/src/android/ISSM/res/layout/activity_mapselection.xml	(revision 14067)
@@ -0,0 +1,27 @@
+
+    <RelativeLayout xmlns:android="http://schemas.android.com/apk/res/android"
+        android:layout_width="fill_parent"
+        android:layout_height="fill_parent" >
+
+        <Button
+            android:id="@+id/button1"
+            style="@layout/activity_mapselection"
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:layout_alignParentLeft="true"
+            android:layout_alignParentTop="true"
+            android:layout_marginLeft="63dp"
+            android:layout_marginTop="97dp"
+            android:text="Greenland" />
+
+        <Button
+            android:id="@+id/button2"
+            android:layout_width="wrap_content"
+            android:layout_height="wrap_content"
+            android:layout_alignBaseline="@+id/button1"
+            android:layout_alignBottom="@+id/button1"
+            android:layout_marginLeft="64dp"
+            android:layout_toRightOf="@+id/button1"
+            android:text="Antarctica" />
+
+    </RelativeLayout>
Index: /issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java
===================================================================
--- /issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java	(revision 14066)
+++ /issm/trunk/src/android/ISSM/src/com/example/issm/ISSM.java	(revision 14067)
@@ -1,6 +1,13 @@
 package com.example.issm;
+
+import java.io.IOException;
+import java.io.InputStream;
+import java.nio.ByteBuffer;
+import java.nio.ByteOrder;
+import java.nio.DoubleBuffer;
 
 import android.os.Bundle;
 import android.app.Activity;
+import android.content.res.AssetManager;
 import android.text.TextUtils;
 import android.view.Menu;
@@ -12,10 +19,25 @@
 
 
-public class ISSM extends Activity implements OnClickListener {
+public class ISSM extends Activity implements OnClickListener 
+{
 	private EditText input;
 	private TextView output;
+	private DoubleBuffer buff;
+	private IssmJni issmNative;
+	private String mapName;
+	private String issmFolder;
+	private int size;
     @Override
+  //------------------------------------------------------------------------------------------------    
     public void onCreate(Bundle savedInstanceState) {
         super.onCreate(savedInstanceState);
+        Bundle map = getIntent().getExtras();
+        {
+        	if(map!= null)
+        	{
+        		mapName = map.getString("map");
+        		issmFolder = map.getString("pathToFile");
+        	}
+        }
         setContentView(R.layout.activity_issm);
         this.input  = (EditText) super.findViewById(R.id.input);
@@ -24,16 +46,44 @@
         button.setOnClickListener(this);
         
+        //load up the ISSM library and create double buffer in java
+        //which later on will be pass for native allocation.
+        issmNative = new IssmJni();
+        buff = ByteBuffer.allocateDirect(15*8).order(ByteOrder.nativeOrder()).asDoubleBuffer();
+        this.createModel();
     }
-    
-	public void onClick(View view) 
+//------------------------------------------------------------------------------------------------    
+    public void createModel()
+    {
+    	String file = "";
+    	if( mapName.equals("greenland"))
+		{
+    		file = "test102.petsc";
+		}
+    	else file = "test102.petsc";
+    	size = issmNative.createISSMModel(issmFolder + file);
+    }
+//------------------------------------------------------------------------------------------------
+    public void fillBuffer()
+    {
+    	issmNative.processBuffer(5,buff);
+    }
+ //------------------------------------------------------------------------------------------------   
+    public void onClick(View view) 
 	{
-		// TODO Auto-generated method stub
+		//factorial method
 		String input = this.input.getText().toString();
 		if(TextUtils.isEmpty(input))
 		{
 			return;
-		}
-		long result = IssmJni.facIterative(Long.parseLong(input));
-		this.output.setText("Result = " + result + "\n");
+		}		
+		long resultfromFac = issmNative.fac(Long.parseLong(input));
+		//example of how to fill buffer Native
+		this.fillBuffer();
+		
+		//print result from fac and the first two slot of filled buffer.
+		this.output.setText("Result = " + resultfromFac + "\n" 
+										+ "First slot from buffer:\n" 
+										+ buff.get(0) + " size " +  size);
+
 	}
 }
Index: /issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java
===================================================================
--- /issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java	(revision 14066)
+++ /issm/trunk/src/android/ISSM/src/com/example/issm/IssmJni.java	(revision 14067)
@@ -1,13 +1,13 @@
 package com.example.issm;
+import java.nio.DoubleBuffer;
 
-public class IssmJni
+class IssmJni
 {
-	public static native long fac(long n);
-	static {
+	public native long fac(long n);
+	public native void processBuffer(int alpha, DoubleBuffer buff);
+	public native int createISSMModel(String file);
+	static 
+	{
         System.loadLibrary("IssmJni");
     }
-	public static long facIterative(long n)
-	{
-		return fac(n);
-	}
 }
Index: /issm/trunk/src/android/ISSM/src/com/example/issm/MapSelection.java
===================================================================
--- /issm/trunk/src/android/ISSM/src/com/example/issm/MapSelection.java	(revision 14067)
+++ /issm/trunk/src/android/ISSM/src/com/example/issm/MapSelection.java	(revision 14067)
@@ -0,0 +1,143 @@
+package com.example.issm;
+
+import java.io.File;
+import java.io.FileOutputStream;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.OutputStream;
+
+import android.app.Activity;
+import android.os.Bundle;
+import android.os.Environment;
+import android.text.TextUtils;
+import android.util.Log;
+import android.view.Menu;
+import android.view.View;
+import android.widget.Button;
+import android.widget.Toast;
+import android.content.Intent;
+import android.content.SharedPreferences;
+import android.content.res.AssetManager;
+/////////////////////////////////////////////////////////////////////////// 
+public class MapSelection extends Activity
+{
+	ISSM issm = new ISSM();
+	private static final String PREFERENCE_FIRST_RUN = null;
+	private String extStorageDirectory;
+	private String issmFolder;
+//------------------------------------------------------------------------	
+	public void onCreate(Bundle icicle)
+	{
+		super.onCreate(icicle);
+		setContentView(R.layout.activity_mapselection);
+		
+		 SharedPreferences settings = this.getSharedPreferences("ISSM", 0);
+		    boolean firstrun = settings.getBoolean(PREFERENCE_FIRST_RUN, true);
+		    extStorageDirectory = Environment.getExternalStorageDirectory().toString();
+		    
+		    issmFolder = extStorageDirectory + "/ISSM/input_files/";
+		    if (firstrun) 
+		    { // Checks to see if we've ran the application b4
+		        SharedPreferences.Editor e = settings.edit();
+		        e.putBoolean(PREFERENCE_FIRST_RUN, false);
+		        
+		        e.commit();
+		        // If not, run these methods:
+		        SetDirectory(issmFolder);
+		    }
+		    
+	
+		//this button represents greenland and pass signal to issm
+		Button gl = (Button) findViewById(R.id.button1);		
+		gl.setOnClickListener(new View.OnClickListener() 
+		{
+			public void onClick(View v) 
+			{
+				Intent i = new Intent(MapSelection.this, ISSM.class);
+				i.putExtra("map", "greenland");
+				i.putExtra("pathToFile", issmFolder);
+		        startActivity(i);
+			}
+		});
+		
+		//this button represents artarctica and pass signal to issm
+		Button art = (Button) findViewById(R.id.button2);
+		art.setOnClickListener(new View.OnClickListener() 
+		{
+			public void onClick(View v) 
+			{
+				Intent i = new Intent(MapSelection.this, ISSM.class);
+				i.putExtra("map", "antarctica");
+				i.putExtra("pathToFile", issmFolder);
+		        startActivity(i);
+			}
+		});
+	}
+//------------------------------------------------------------------------
+	private void SetDirectory(String directory) 
+	{
+	    if (android.os.Environment.getExternalStorageState().equals(android.os.Environment.MEDIA_MOUNTED)) 
+	    {
+	    	System.out.println(issmFolder);
+	        File txtDirectory = new File(issmFolder);
+	        
+	        // Check and create directory in SDcard
+	        if(!txtDirectory.exists())
+	        {
+	        	txtDirectory.mkdirs();
+	        	System.out.println("making directory");
+	        }
+	        CopyAssets(); // Then run the method to copy the file.
+
+	    } 
+	    else if (android.os.Environment.getExternalStorageState().equals(android.os.Environment.MEDIA_MOUNTED_READ_ONLY)) 
+	    {
+	    	Toast toast = Toast.makeText(this.getApplicationContext(), "Memory is not mounted to device", Toast.LENGTH_LONG);
+	    	toast.show();
+	    }
+
+	}
+//----------------------------------------------------------------------------
+	/**
+	 * -- Copy the file from the assets folder to the sdCard
+	 * ===========================================================
+	 **/
+		private void CopyAssets() 
+		{
+		    AssetManager assetManager = getAssets();
+		    String[] files = null;
+		    try {
+		        files = assetManager.list("Map");
+		    } catch (IOException e) {
+		        Log.e("tag", e.getMessage());
+		    }
+		    for (int i = 0; i < files.length; i++) {
+		        InputStream in = null;
+		        OutputStream out = null;
+		        try {
+		            in = assetManager.open("Map/"+files[i]);
+		            out = new FileOutputStream(issmFolder + files[i]);
+		            copyFile(in, out);
+		            in.close();
+		            in = null;
+		            out.flush();
+		            out.close();
+		            out = null;
+		        } catch (Exception e) {
+		            Log.e("tag", e.getMessage());
+		        }
+		    }
+		    System.out.println("Done");
+		}
+//----------------------------------------------------------------------------
+		private void copyFile(InputStream in, OutputStream out) throws IOException 
+		{
+		    byte[] buffer = new byte[1024];
+		    int read;
+		    while ((read = in.read(buffer)) != -1) {
+		        out.write(buffer, 0, read);
+		    }
+		}
+//----------------------------------------------------------------------------
+}
+/////////////////////////////////////////////////////////////////////////// 
Index: /issm/trunk/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk/src/c/Container/Observations.cpp	(revision 14066)
+++ /issm/trunk/src/c/Container/Observations.cpp	(revision 14067)
@@ -424,4 +424,88 @@
 	*pprediction = obs;
 }/*}}}*/
+/*FUNCTION Observations::InterpolationV4{{{*/
+void Observations::InterpolationV4(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata){
+	/* Reference:  David T. Sandwell, Biharmonic spline interpolation of GEOS-3
+	 * and SEASAT altimeter data, Geophysical Research Letters, 2, 139-142,
+	 * 1987.  Describes interpolation using value or gradient of value in any
+	 * dimension.*/
+
+	/*Intermediaries*/
+	int         i,j,n_obs;
+	IssmPDouble prediction,h;
+	IssmPDouble *x       = NULL;
+	IssmPDouble *y       = NULL;
+	IssmPDouble *obs     = NULL;
+	IssmPDouble *Green   = NULL;
+	IssmPDouble *weights = NULL;
+	IssmPDouble *g       = NULL;
+
+	/*Some checks*/
+	_assert_(maxdata>0);
+	_assert_(pprediction);
+
+	/*If radius is not provided or is 0, return all observations*/
+	if(radius==0) radius=this->quadtree->root->length;
+
+	/*Get list of observations for current point*/
+	this->ObservationList(&x,&y,&obs,&n_obs,x_interp,y_interp,radius,maxdata);
+
+	/*If we have less observations than mindata, return UNDEF*/
+	if(n_obs<mindata || n_obs<2){
+		prediction = UNDEF; 
+	}
+	else{
+
+		/*Allocate intermediary matrix and vectors*/
+		Green = xNew<IssmPDouble>(n_obs*n_obs);
+		g     = xNew<IssmPDouble>(n_obs);
+
+		/*First: distance vector*/
+		for(i=0;i<n_obs;i++){
+			h = sqrt( (x[i]-x_interp)*(x[i]-x_interp) + (y[i]-y_interp)*(y[i]-y_interp) );
+			if(h>0){
+				g[i] = h*h*(log(h)-1.);
+			}
+			else{
+				g[i] = 0.;
+			}
+		}
+
+		/*Build Green function matrix*/
+		for(i=0;i<n_obs;i++){
+			for(j=0;j<=i;j++){
+				h = sqrt( (x[i]-x[j])*(x[i]-x[j]) + (y[i]-y[j])*(y[i]-y[j]) );
+				if(h>0){
+					Green[j*n_obs+i] = h*h*(log(h)-1.);
+				}
+				else{
+					Green[j*n_obs+i] = 0.;
+				}
+				Green[i*n_obs+j] = Green[j*n_obs+i];
+			}
+		}
+
+		/*Compute weights*/
+#if _HAVE_GSL_
+		SolverxSeq(&weights,Green,obs,n_obs); // Green^-1 obs
+#else
+		_error_("GSL is required");
+#endif
+
+		/*Interpolate*/
+		prediction = 0;
+		for(i=0;i<n_obs;i++) prediction += weights[i]*g[i];
+
+	}
+
+	/*clean-up*/
+	*pprediction = prediction;
+	xDelete<IssmPDouble>(x);
+	xDelete<IssmPDouble>(y);
+	xDelete<IssmPDouble>(obs);
+	xDelete<IssmPDouble>(Green);
+	xDelete<IssmPDouble>(g);
+	xDelete<IssmPDouble>(weights);
+}/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
 void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
Index: /issm/trunk/src/c/Container/Observations.h
===================================================================
--- /issm/trunk/src/c/Container/Observations.h	(revision 14066)
+++ /issm/trunk/src/c/Container/Observations.h	(revision 14067)
@@ -28,4 +28,5 @@
 		void ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
 		void InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power);
+		void InterpolationV4(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata);
 		void InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram);
 		void InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius);
Index: /issm/trunk/src/c/Container/Options.h
===================================================================
--- /issm/trunk/src/c/Container/Options.h	(revision 14066)
+++ /issm/trunk/src/c/Container/Options.h	(revision 14067)
@@ -64,4 +64,5 @@
 			}
 			else{
+				if(GetOption(name)) _printLine_("WARNING: option "<<name<<" found but fetched format not consistent, defaulting...");
 				*pvalue=default_value;
 			}
Index: /issm/trunk/src/c/android/fac.cpp
===================================================================
--- /issm/trunk/src/c/android/fac.cpp	(revision 14066)
+++ /issm/trunk/src/c/android/fac.cpp	(revision 14067)
@@ -13,12 +13,14 @@
 
 /*}}}*/
-
+#include "stdio.h"
 long fac::factorial(long n) {
 	long f = 1;
 	long i;
+	printf("ok1\n");
 	
 	for(i = 1; i <= n; i++){
 		f *= i;
 	}
+	printf("ok2\n");
 	
 	return f;
Index: /issm/trunk/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk/src/c/classes/FemModel.cpp	(revision 14066)
+++ /issm/trunk/src/c/classes/FemModel.cpp	(revision 14067)
@@ -20,4 +20,9 @@
 
 /*Object constructors and destructor*/
+/*FUNCTION FemModel::FemModel(char* filename) {{{*/
+FemModel::FemModel(char* filename){
+
+}
+/*}}}*/
 /*FUNCTION FemModel::FemModel(int argc,char** argv){{{*/
 FemModel::FemModel(int argc,char** argv,COMM incomm){
@@ -502,5 +507,7 @@
 
 			/*Loop over elements that hold node number i*/
-			_assert_(head_e[node->Sid()]!=-1 || head_l[node->Sid()]!=-1);
+			//if(head_e[node->Sid()]==-1 && head_l[node->Sid()]==-1){
+			//	printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Sid()+1);
+			//}
 			for(j=head_e[node->Sid()];j!=-1;j=next_e[j]){
 				offset=count2offset_e[j];
Index: /issm/trunk/src/c/classes/FemModel.h
===================================================================
--- /issm/trunk/src/c/classes/FemModel.h	(revision 14066)
+++ /issm/trunk/src/c/classes/FemModel.h	(revision 14067)
@@ -47,4 +47,5 @@
 
 		/*constructors, destructors: */
+		FemModel(char* filename);
 		FemModel(int argc,char** argv,COMM comm_init);
 		FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, const int solution_type,const int* analyses,const int nummodels);
Index: /issm/trunk/src/c/classes/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/classes/objects/Elements/Penta.cpp	(revision 14066)
+++ /issm/trunk/src/c/classes/objects/Elements/Penta.cpp	(revision 14067)
@@ -66,7 +66,7 @@
 
 	/*Build neighbors list*/
-	if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index])) penta_elements_ids[1]=this->id; //upper penta is the same penta
+	if (xIsNan<IssmDouble>(iomodel->Data(MeshUpperelementsEnum)[index]) || iomodel->Data(MeshUpperelementsEnum)[index]==-1.) penta_elements_ids[1]=this->id; //upper penta is the same penta
 	else                                    penta_elements_ids[1]=reCast<int,IssmDouble>((iomodel->Data(MeshUpperelementsEnum)[index]));
-	if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index])) penta_elements_ids[0]=this->id; //lower penta is the same penta
+	if (xIsNan<IssmDouble>(iomodel->Data(MeshLowerelementsEnum)[index]) || iomodel->Data(MeshLowerelementsEnum)[index]==-1.) penta_elements_ids[0]=this->id; //lower penta is the same penta
 	else                                    penta_elements_ids[0]=reCast<int,IssmDouble>((iomodel->Data(MeshLowerelementsEnum)[index]));
 	this->InitHookNeighbors(penta_elements_ids);
Index: /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp
===================================================================
--- /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp	(revision 14066)
+++ /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.cpp	(revision 14067)
@@ -15,7 +15,7 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-int hascommondedge(double* element1,double* element2);
+int hascommondedge(int* element1,int* element2);
 
-void	ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width){
+void ElementConnectivityx(int** pelementconnectivity,int* elements, int nels,int* nodeconnectivity, int nods, int width){
 
 	int i,j,k,n;
@@ -23,6 +23,5 @@
 	/*intermediary: */
 	int    maxels;
-	double element;
-	double connectedelement;
+	int  connectedelement;
 	int    connectedelementindex;
 	int    node;
@@ -30,44 +29,41 @@
 	int    num_elements;
 
-	/*output: */
-	double* elementconnectivity=NULL;
-
 	/*maxels: */
 	maxels=width-1;
+
 	/*Allocate connectivity: */
-	elementconnectivity=xNewZeroInit<double>(nel*3);
+	int* elementconnectivity=xNewZeroInit<int>(nels*3);
 
 	/*Go through all elements, and for each element, go through its nodes, to get the neighbouring elements. 
 	 * Once we get the neighbouring elements, figure out if they share a segment with the current element. If so, 
 	 * plug them in the connectivity, unless they are already there.: */
+	for(n=0;n<nels;n++){
 
-	for(n=0;n<nel;n++){
-
-		element=(double)(n+1); //matlab indexing
+		//element=n+1; //matlab indexing
 
 		for(i=0;i<3;i++){
 
-			node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.
+			node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace.
 			index=node-1;
 
-			num_elements=(int)*(nodeconnectivity+width*index+maxels); //retrieve number of elements already  plugged into the connectivity of this node.
+			num_elements=nodeconnectivity[width*index+maxels]; //retrieve number of elements already  plugged into the connectivity of this node.
 
 			for(j=0;j<num_elements;j++){
 
 				/*for each element connected to node, figure out if it has a commond edge with element: */
-				connectedelement=*(nodeconnectivity+width*index+j);
-				connectedelementindex=(int)(connectedelement-1); //go from matlab indexing to c indexing.
+				connectedelement=nodeconnectivity[width*index+j];
+				connectedelementindex=connectedelement-1; //go from matlab indexing to c indexing.
 
-				if(hascommondedge(elements+n*3+0,elements+connectedelementindex*3+0)){
+				if(hascommondedge(&elements[n*3+0],&elements[connectedelementindex*3+0])){
 					/*Ok, this connected element has a commond edge  with element, plug it into elementconnectivity, unless 
 					 *it is already there: */
 
 					for(k=0;k<3;k++){
-						if (*(elementconnectivity+3*n+k)==0){
-							*(elementconnectivity+3*n+k)=connectedelement;
+						if(elementconnectivity[3*n+k]==0){
+							elementconnectivity[3*n+k]=connectedelement;
 							break;
 						}
 						else{
-							if(connectedelement==*(elementconnectivity+3*n+k))break;
+							if(connectedelement==elementconnectivity[3*n+k]) break;
 						}
 					}
@@ -81,16 +77,15 @@
 }
 
-int hascommondedge(double* element1,double* element2){
+int hascommondedge(int* el1,int* el2){
 
-	int i,j;
-	int count;
-
-	count=0;
-	for(i=0;i<3;i++){
-		for(j=0;j<3;j++){
-			if (*(element1+i)==*(element2+j))count++;
+	int count=0;
+	for(int i=0;i<3;i++){
+		for(int j=0;j<3;j++){
+			if(el1[i]==el2[j]) count++;
 		}
 	}
-	if (count==2)return 1;
-	else return 0;
+	if(count==2)
+	 return 1;
+	else
+	 return 0;
 }
Index: /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h
===================================================================
--- /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h	(revision 14066)
+++ /issm/trunk/src/c/modules/ElementConnectivityx/ElementConnectivityx.h	(revision 14067)
@@ -7,5 +7,5 @@
 
 /* local prototypes: */
-void	ElementConnectivityx( double** pelementconnectivity, double* elements, int nel, double* nodeconnectivity, int nods, int width);
+void	ElementConnectivityx(int** pelementconnectivity,int* elements,int nels,int* nodeconnectivity, int nods, int width);
 
 #endif  /* _ELEMENTCONNECTIVITYX_H */
Index: /issm/trunk/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Krigingx/Krigingx.cpp	(revision 14066)
+++ /issm/trunk/src/c/modules/Krigingx/Krigingx.cpp	(revision 14067)
@@ -22,4 +22,5 @@
 	/*Intermediaries*/
 	int           mindata,maxdata;
+	double        dmindata,dmaxdata,dnumthreads; //FIXME (Options come as double but we want to retrive integers)
 	double        radius;
 	char         *output       = NULL;
@@ -37,6 +38,7 @@
 	ProcessVariogram(&variogram,options);
 	options->Get(&radius,"searchradius",0.);
-	options->Get(&mindata,"mindata",1);
-	options->Get(&maxdata,"maxdata",50);
+	options->Get(&dmindata,"mindata",1.);  mindata=int(dmindata);//FIXME (Options come as double but we want to retrive integers)
+	options->Get(&dmaxdata,"maxdata",50.); maxdata=int(dmaxdata);//FIXME (Options come as double but we want to retrive integers)
+	options->Get(&dnumthreads,"numthreads",double(num)); num=int(dnumthreads);//FIXME (Options come as double but we want to retrive integers)
 
 	/*Process observation dataset*/
@@ -88,10 +90,10 @@
 		gate.predictions  = predictions;
 		gate.error        = error;
-		gate.percent      = xNewZeroInit<double>(num);
+		gate.numdone      = xNewZeroInit<int>(num);
 
 		/*launch the thread manager with Krigingxt as a core: */
 		LaunchThread(NearestNeighbort,(void*)&gate,num);
 		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
-		xDelete<double>(gate.percent);
+		xDelete<int>(gate.numdone);
 	}
 	else if(strcmp(output,"idw")==0){ //Inverse distance weighting
@@ -109,5 +111,5 @@
 		gate.predictions  = predictions;
 		gate.error        = error;
-		gate.percent      = xNewZeroInit<double>(num);
+		gate.numdone      = xNewZeroInit<int>(num);
 		gate.power        = power;
 
@@ -115,8 +117,10 @@
 		LaunchThread(idwt,(void*)&gate,num);
 		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
-		xDelete<double>(gate.percent);
-	}
-	else if(strcmp(output,"prediction")==0){
-
+		xDelete<int>(gate.numdone);
+	}
+	else if(strcmp(output,"v4")==0){ //Inverse distance weighting
+#if !defined(_HAVE_GSL_)
+		_error_("GSL is required for v4 interpolation");
+#endif
 		/*initialize thread parameters: */
 		gate.n_interp     = n_interp;
@@ -130,10 +134,33 @@
 		gate.predictions  = predictions;
 		gate.error        = error;
-		gate.percent      = xNewZeroInit<double>(num);
+		gate.numdone      = xNewZeroInit<int>(num);
+
+		/*launch the thread manager with Krigingxt as a core: */
+		LaunchThread(v4t,(void*)&gate,num);
+		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
+		xDelete<int>(gate.numdone);
+	}
+	else if(strcmp(output,"prediction")==0){
+#if !defined(_HAVE_GSL_)
+		_error_("GSL is required for v4 interpolation");
+#endif
+
+		/*initialize thread parameters: */
+		gate.n_interp     = n_interp;
+		gate.x_interp     = x_interp;
+		gate.y_interp     = y_interp;
+		gate.radius       = radius;
+		gate.mindata      = mindata;
+		gate.maxdata      = maxdata;
+		gate.variogram    = variogram;
+		gate.observations = observations;
+		gate.predictions  = predictions;
+		gate.error        = error;
+		gate.numdone      = xNewZeroInit<int>(num);
 
 		/*launch the thread manager with Krigingxt as a core: */
 		LaunchThread(Krigingxt,(void*)&gate,num);
 		_printLine_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%");
-		xDelete<double>(gate.percent);
+		xDelete<int>(gate.numdone);
 	}
 	else{
@@ -176,8 +203,5 @@
 	double       *predictions  = gate->predictions;
 	double       *error        = gate->error;
-	double       *percent      = gate->percent;
-
-	/*Intermediaries*/
-	double        localpercent;
+	int          *numdone      = gate->numdone;
 
 	/*partition loop across threads: */
@@ -186,8 +210,10 @@
 
 		/*Print info*/
-		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
-		localpercent=percent[0];
-		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
-		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+		numdone[my_thread]=idx-i0;
+		if(my_thread==0){
+			int alldone=numdone[0];
+			for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+			_printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+		}
 
 		/*Kriging interpolation*/
@@ -224,9 +250,5 @@
 	double       *predictions  = gate->predictions;
 	double       *error        = gate->error;
-	double       *percent      = gate->percent;
-
-	/*Intermediaries*/
-	int           i;
-	double        localpercent;
+	int          *numdone      = gate->numdone;
 
 	/*partition loop across threads: */
@@ -235,8 +257,10 @@
 
 		/*Print info*/
-		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
-		localpercent=percent[0];
-		for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
-		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+		numdone[my_thread]=idx-i0;
+		if(my_thread==0){
+			int alldone=numdone[0];
+			for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+			_printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+		}
 
 		observations->InterpolationNearestNeighbor(&predictions[idx],x_interp[idx],y_interp[idx],radius);
@@ -272,9 +296,6 @@
 	double       *predictions  = gate->predictions;
 	double       *error        = gate->error;
-	double       *percent      = gate->percent;
+	int          *numdone      = gate->numdone;
 	double        power        = gate->power;
-
-	/*Intermediaries*/
-	double localpercent;
 
 	/*partition loop across threads: */
@@ -283,10 +304,57 @@
 
 		/*Print info*/
-		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
-		localpercent=percent[0];
-		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
-		if(my_thread==0) _printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<localpercent<<"%");
+		numdone[my_thread]=idx-i0;
+		if(my_thread==0){
+			int alldone=numdone[0];
+			for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+			_printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+		}
 
 		observations->InterpolationIDW(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,power);
+	}
+	return NULL;
+}/*}}}*/
+/*FUNCTION v4t{{{*/
+void* v4t(void* vpthread_handle){
+
+	/*gate variables :*/
+	KrigingxThreadStruct *gate        = NULL;
+	pthread_handle       *handle      = NULL;
+	int my_thread;
+	int num_threads;
+	int i0,i1;
+
+	/*recover handle and gate: */
+	handle      = (pthread_handle*)vpthread_handle;
+	gate        = (KrigingxThreadStruct*)handle->gate;
+	my_thread   = handle->id;
+	num_threads = handle->num;
+
+	/*recover parameters :*/
+	int           n_interp     = gate->n_interp;
+	double       *x_interp     = gate->x_interp;
+	double       *y_interp     = gate->y_interp;
+	double        radius       = gate->radius;
+	int           mindata      = gate->mindata;
+	int           maxdata      = gate->maxdata;
+	Variogram    *variogram    = gate->variogram;
+	Observations *observations = gate->observations;
+	double       *predictions  = gate->predictions;
+	double       *error        = gate->error;
+	int          *numdone      = gate->numdone;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+	for(int idx=i0;idx<i1;idx++){
+
+		/*Print info*/
+		numdone[my_thread]=idx-i0;
+		if(my_thread==0){
+			int alldone=numdone[0];
+			for(int i=1;i<num_threads;i++) alldone+=numdone[i];
+			_printString_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(alldone)/double(n_interp)*100.<<"%");
+		}
+
+		observations->InterpolationV4(&predictions[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata);
 	}
 	return NULL;
Index: /issm/trunk/src/c/modules/Krigingx/Krigingx.h
===================================================================
--- /issm/trunk/src/c/modules/Krigingx/Krigingx.h	(revision 14066)
+++ /issm/trunk/src/c/modules/Krigingx/Krigingx.h	(revision 14067)
@@ -28,5 +28,5 @@
 	double       *predictions;
 	double       *error;
-	double       *percent;
+	int          *numdone;
 	double        power;//for idw
 }KrigingxThreadStruct;
@@ -35,3 +35,4 @@
 void* NearestNeighbort(void*);
 void* idwt(void*);
+void* v4t(void*);
 #endif /* _KRIGINGX_H */
Index: /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp
===================================================================
--- /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp	(revision 14066)
+++ /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp	(revision 14067)
@@ -18,5 +18,5 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void	NodeConnectivityx( double** pconnectivity, int* pwidth, double* elements, int nel, int nods){
+void	NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements, int nels, int nods){
 
 	int i,j,n;
@@ -29,24 +29,21 @@
 	int     num_elements;
 	int     already_plugged=0;
-	double  element;
-
-	/*output: */
-	double* connectivity=NULL;
+	int element;
 
 	/*Allocate connectivity: */
-	connectivity=xNewZeroInit<double>(nods*width);
+	int* connectivity=xNewZeroInit<int>(nods*width);
 
 	/*Go through all elements, and for each elements, plug into the connectivity, all the nodes. 
 	 * If nodes are already plugged into the connectivity, skip them.: */
-	for(n=0;n<nel;n++){
+	for(n=0;n<nels;n++){
 
-		element=(double)(n+1); //matlab indexing
+		element=n+1; //matlab indexing
 
 		for(i=0;i<3;i++){
 
-			node=(int)*(elements+n*3+i); //already matlab indexed, elements comes directly from the workspace.
+			node=elements[n*3+i]; //already matlab indexed, elements comes directly from the workspace.
 			index=node-1;
 
-			num_elements=(int)*(connectivity+width*index+maxels); //retrieve number of elements already  plugged into the connectivity of this node.
+			num_elements=connectivity[width*index+maxels]; //retrieve number of elements already  plugged into the connectivity of this node.
 
 			already_plugged=0;
@@ -60,6 +57,6 @@
 
 			/*this elements is not yet plugged  into the connectivity for this node, do it, and increase counter: */
-			*(connectivity+width*index+num_elements)=element;
-			*(connectivity+width*index+maxels)=(double)(num_elements+1);
+			connectivity[width*index+num_elements]=element;
+			connectivity[width*index+maxels]=num_elements+1;
 
 		}
@@ -69,5 +66,6 @@
 	 * warn the user to increase the connectivity width: */
 	for(i=0;i<nods;i++){
-		if (*(connectivity+width*i+maxels)>maxels)_error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table");
+		if (*(connectivity+width*i+maxels)>maxels)
+		 _error_("max connectivity width reached (" << *(connectivity+width*i+maxels) << ")! increase width of connectivity table");
 	}
 
Index: /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h
===================================================================
--- /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h	(revision 14066)
+++ /issm/trunk/src/c/modules/NodeConnectivityx/NodeConnectivityx.h	(revision 14067)
@@ -7,5 +7,5 @@
 
 /* local prototypes: */
-void	NodeConnectivityx( double** pconnectivity, int* pwidth,double* elements, int nel, int nods);
+void	NodeConnectivityx(int** pconnectivity,int* pwidth,int* elements,int nels, int nods);
 
 #endif  /* _NODECONNECTIVITYX_H */
Index: /issm/trunk/src/c/shared/Threads/LaunchThread.cpp
===================================================================
--- /issm/trunk/src/c/shared/Threads/LaunchThread.cpp	(revision 14066)
+++ /issm/trunk/src/c/shared/Threads/LaunchThread.cpp	(revision 14067)
@@ -31,4 +31,8 @@
 	pthread_t      *threads = NULL;
 	pthread_handle *handles = NULL;
+
+	/*check number of threads*/
+	if(num_threads<1)    _error_("number of threads must be at least 1");
+	if(num_threads>2000) _error_("number of threads seems to be too high ("<<num_threads<<">2000)");
 
 	/*dynamically allocate: */
Index: /issm/trunk/src/dox/issm.dox
===================================================================
--- /issm/trunk/src/dox/issm.dox	(revision 14066)
+++ /issm/trunk/src/dox/issm.dox	(revision 14067)
@@ -46,23 +46,23 @@
 </th>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">505</td><td  bgcolor=#FFFFFF style="text-align:right;">14316</td><td  bgcolor=#FFFFFF style="text-align:right;">16564</td><td  bgcolor=#FFFFFF style="text-align:right;">56522</td><td  bgcolor=#FFFFFF style="text-align:right;">87402</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> C++ </th><td  bgcolor=#FFFFFF style="text-align:right;">469</td><td  bgcolor=#FFFFFF style="text-align:right;">14087</td><td  bgcolor=#FFFFFF style="text-align:right;">16669</td><td  bgcolor=#FFFFFF style="text-align:right;">56706</td><td  bgcolor=#FFFFFF style="text-align:right;">87462</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">913</td><td  bgcolor=#C6E2FF style="text-align:right;">6625</td><td  bgcolor=#C6E2FF style="text-align:right;">13054</td><td  bgcolor=#C6E2FF style="text-align:right;">30160</td><td  bgcolor=#C6E2FF style="text-align:right;">49839</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> MATLAB </th><td  bgcolor=#C6E2FF style="text-align:right;">945</td><td  bgcolor=#C6E2FF style="text-align:right;">6583</td><td  bgcolor=#C6E2FF style="text-align:right;">13346</td><td  bgcolor=#C6E2FF style="text-align:right;">30647</td><td  bgcolor=#C6E2FF style="text-align:right;">50576</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">388</td><td  bgcolor=#FFFFFF style="text-align:right;">3069</td><td  bgcolor=#FFFFFF style="text-align:right;">2983</td><td  bgcolor=#FFFFFF style="text-align:right;">11855</td><td  bgcolor=#FFFFFF style="text-align:right;">17907</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> C/C++  Header </th><td  bgcolor=#FFFFFF style="text-align:right;">362</td><td  bgcolor=#FFFFFF style="text-align:right;">2772</td><td  bgcolor=#FFFFFF style="text-align:right;">2894</td><td  bgcolor=#FFFFFF style="text-align:right;">11914</td><td  bgcolor=#FFFFFF style="text-align:right;">17580</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> m4 </th><td  bgcolor=#C6E2FF style="text-align:right;">6</td><td  bgcolor=#C6E2FF style="text-align:right;">1052</td><td  bgcolor=#C6E2FF style="text-align:right;">79</td><td  bgcolor=#C6E2FF style="text-align:right;">8993</td><td  bgcolor=#C6E2FF style="text-align:right;">10124</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> Python </th><td  bgcolor=#C6E2FF style="text-align:right;">109</td><td  bgcolor=#C6E2FF style="text-align:right;">3361</td><td  bgcolor=#C6E2FF style="text-align:right;">4763</td><td  bgcolor=#C6E2FF style="text-align:right;">7050</td><td  bgcolor=#C6E2FF style="text-align:right;">15174</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Python </th><td  bgcolor=#FFFFFF style="text-align:right;">81</td><td  bgcolor=#FFFFFF style="text-align:right;">2677</td><td  bgcolor=#FFFFFF style="text-align:right;">3936</td><td  bgcolor=#FFFFFF style="text-align:right;">5169</td><td  bgcolor=#FFFFFF style="text-align:right;">11782</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> m4 </th><td  bgcolor=#FFFFFF style="text-align:right;">1</td><td  bgcolor=#FFFFFF style="text-align:right;">187</td><td  bgcolor=#FFFFFF style="text-align:right;">5</td><td  bgcolor=#FFFFFF style="text-align:right;">1247</td><td  bgcolor=#FFFFFF style="text-align:right;">1439</td>
 </tr>
 <tr>
-<th  bgcolor=#C6E2FF style="text-align:left;"> Objective  C </th><td  bgcolor=#C6E2FF style="text-align:right;">9</td><td  bgcolor=#C6E2FF style="text-align:right;">98</td><td  bgcolor=#C6E2FF style="text-align:right;">0</td><td  bgcolor=#C6E2FF style="text-align:right;">370</td><td  bgcolor=#C6E2FF style="text-align:right;">468</td>
+<th  bgcolor=#C6E2FF style="text-align:left;"> Objective  C </th><td  bgcolor=#C6E2FF style="text-align:right;">9</td><td  bgcolor=#C6E2FF style="text-align:right;">96</td><td  bgcolor=#C6E2FF style="text-align:right;">0</td><td  bgcolor=#C6E2FF style="text-align:right;">370</td><td  bgcolor=#C6E2FF style="text-align:right;">466</td>
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">5</td><td  bgcolor=#FFFFFF style="text-align:right;">58</td><td  bgcolor=#FFFFFF style="text-align:right;">81</td><td  bgcolor=#FFFFFF style="text-align:right;">268</td><td  bgcolor=#FFFFFF style="text-align:right;">407</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> Bourne  Shell </th><td  bgcolor=#FFFFFF style="text-align:right;">5</td><td  bgcolor=#FFFFFF style="text-align:right;">57</td><td  bgcolor=#FFFFFF style="text-align:right;">79</td><td  bgcolor=#FFFFFF style="text-align:right;">256</td><td  bgcolor=#FFFFFF style="text-align:right;">392</td>
 </tr>
 <tr>
@@ -76,5 +76,5 @@
 </tr>
 <tr>
-<th  bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td  bgcolor=#FFFFFF style="text-align:right;">1912</td><td  bgcolor=#FFFFFF style="text-align:right;">27936</td><td  bgcolor=#FFFFFF style="text-align:right;">36727</td><td  bgcolor=#FFFFFF style="text-align:right;">113707</td><td  bgcolor=#FFFFFF style="text-align:right;">178370</td>
+<th  bgcolor=#FFFFFF style="text-align:left;"> SUM: </th><td  bgcolor=#FFFFFF style="text-align:right;">1905</td><td  bgcolor=#FFFFFF style="text-align:right;">27184</td><td  bgcolor=#FFFFFF style="text-align:right;">37786</td><td  bgcolor=#FFFFFF style="text-align:right;">108560</td><td  bgcolor=#FFFFFF style="text-align:right;">173530</td>
 </tr>
 </table>
Index: /issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py
===================================================================
--- /issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py	(revision 14066)
+++ /issm/trunk/src/m/boundaryconditions/SetIceShelfBC.py	(revision 14067)
@@ -26,7 +26,7 @@
 			raise IOError("SetIceShelfBC error message: ice front file '%s' not found." % icefrontfile)
 		[nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
-		nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
+		nodeonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
 	else:
-		nodeonicefront=numpy.zeros((md.mesh.numberofvertices))
+		nodeonicefront=numpy.zeros((md.mesh.numberofvertices),bool)
 
 #	pos=find(md.mesh.vertexonboundary & ~nodeonicefront);
@@ -56,5 +56,5 @@
 	#segment on Neumann (Ice Front)
 #	pos=find(nodeonicefront(md.mesh.segments(:,1)) | nodeonicefront(md.mesh.segments(:,2)));
-	pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0].astype(int)-1],nodeonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
+	pos=numpy.nonzero(numpy.logical_or(nodeonicefront[md.mesh.segments[:,0]-1],nodeonicefront[md.mesh.segments[:,1]-1]))[0]
 	if   md.mesh.dimension==2:
 		pressureload=md.mesh.segments[pos,:]
@@ -62,5 +62,5 @@
 #		pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
 		pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
-		pressureload=numpy.zeros((0,5))
+		pressureload=numpy.zeros((0,5),int)
 		for i in xrange(1,md.mesh.numberoflayers):
 #			pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
@@ -69,5 +69,5 @@
 	#Add water or air enum depending on the element
 #	pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))];
-	pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
+	pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)))
 
 	#plug onto model
Index: /issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py
===================================================================
--- /issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 14066)
+++ /issm/trunk/src/m/boundaryconditions/SetMarineIceSheetBC.py	(revision 14067)
@@ -28,10 +28,10 @@
 			raise IOError("SetMarineIceSheetBC error message: ice front file '%s' not found." % icefrontfile)
 		[nodeinsideicefront,dum]=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),icefrontfile,'node',2)
-		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1)).astype(float)
+		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,nodeinsideicefront.reshape(-1))
 	else:
 		#Guess where the ice front is
-		vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices))
-		vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:].astype(int)-1]=1
-		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice).astype(float)
+		vertexonfloatingice=numpy.zeros((md.mesh.numberofvertices),bool)
+		vertexonfloatingice[md.mesh.elements[numpy.nonzero(md.mask.elementonfloatingice),:]-1]=True
+		vertexonicefront=numpy.logical_and(md.mesh.vertexonboundary,vertexonfloatingice)
 
 #	pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
@@ -62,5 +62,5 @@
 	#segment on Neumann (Ice Front)
 #	pos=find(vertexonicefront(md.mesh.segments(:,1)) | vertexonicefront(md.mesh.segments(:,2)));
-	pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0].astype(int)-1],vertexonicefront[md.mesh.segments[:,1].astype(int)-1]))[0]
+	pos=numpy.nonzero(numpy.logical_or(vertexonicefront[md.mesh.segments[:,0]-1],vertexonicefront[md.mesh.segments[:,1]-1]))[0]
 	if   md.mesh.dimension==2:
 		pressureload=md.mesh.segments[pos,:]
@@ -68,5 +68,5 @@
 #		pressureload_layer1=[md.mesh.segments(pos,1:2)  md.mesh.segments(pos,2)+md.mesh.numberofvertices2d  md.mesh.segments(pos,1)+md.mesh.numberofvertices2d  md.mesh.segments(pos,3)];
 		pressureload_layer1=numpy.hstack((md.mesh.segments[pos,0:2],md.mesh.segments[pos,1]+md.mesh.numberofvertices2d,md.mesh.segments[pos,0]+md.mesh.numberofvertices2d,md.mesh.segments[pos,2]))
-		pressureload=numpy.zeros((0,5))
+		pressureload=numpy.zeros((0,5),int)
 		for i in xrange(1,md.mesh.numberoflayers):
 #			pressureload=[pressureload ;pressureload_layer1(:,1:4)+(i-1)*md.mesh.numberofvertices2d pressureload_layer1(:,5)+(i-1)*md.mesh.numberofelements2d ];
@@ -75,5 +75,5 @@
 	#Add water or air enum depending on the element
 #	pressureload=[pressureload 1*md.mask.elementonfloatingice(pressureload(:,end))+ 0*md.mask.elementongroundedice(pressureload(:,end))];
-	pressureload=numpy.hstack((pressureload,1.*md.mask.elementonfloatingice[pressureload[:,-1].astype('int')-1].reshape(-1,1)+0.*md.mask.elementongroundedice[pressureload[:,-1].astype('int')-1].reshape(-1,1)))
+	pressureload=numpy.hstack((pressureload,1*md.mask.elementonfloatingice[pressureload[:,-1]-1].reshape(-1,1)+0*md.mask.elementongroundedice[pressureload[:,-1]-1].reshape(-1,1)))
 
 	#plug onto model
Index: /issm/trunk/src/m/classes/clusters/greenplanet.m
===================================================================
--- /issm/trunk/src/m/classes/clusters/greenplanet.m	(revision 14066)
+++ /issm/trunk/src/m/classes/clusters/greenplanet.m	(revision 14067)
@@ -14,5 +14,5 @@
 		 cpuspernode=8; 
 		 port=8000;
-		 queue='rignot';
+		 queue='c6145';
 		 codepath='';
 		 executionpath='';
@@ -50,5 +50,5 @@
 		 function md = checkconsistency(cluster,md,solution,analyses) % {{{
 
-			 available_queues={'rignot','default'};
+			 available_queues={'c6145','default'};
 			 queue_requirements_time=[Inf Inf];
 			 queue_requirements_np=[80 80];
Index: /issm/trunk/src/m/classes/flowequation.py
===================================================================
--- /issm/trunk/src/m/classes/flowequation.py	(revision 14066)
+++ /issm/trunk/src/m/classes/flowequation.py	(revision 14067)
@@ -19,8 +19,8 @@
 		# {{{ Properties
 		
-		self.ismacayealpattyn     = 0;
-		self.ishutter             = 0;
-		self.isl1l2             = 0;
-		self.isstokes             = 0;
+		self.ismacayealpattyn     = 0
+		self.ishutter             = 0
+		self.isl1l2               = 0
+		self.isstokes             = 0
 		self.vertex_equation      = float('NaN')
 		self.element_equation     = float('NaN')
@@ -37,13 +37,13 @@
 		string='   flow equation parameters:'
 
-		string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn','is the macayeal or pattyn approximation used ?'))
-		string="%s\n%s"%(string,fielddisplay(self,'ishutter','is the shallow ice approximation used ?'))
-		string="%s\n%s"%(string,fielddisplay(self,'isl1l2','are l1l2 equations used ?'))
-		string="%s\n%s"%(string,fielddisplay(self,'isstokes','are the Full-Stokes equations used ?'))
-		string="%s\n%s"%(string,fielddisplay(self,'vertex_equation','flow equation for each vertex'))
-		string="%s\n%s"%(string,fielddisplay(self,'element_equation','flow equation for each element'))
-		string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal','vertices on MacAyeal''s border (for tiling)'))
-		string="%s\n%s"%(string,fielddisplay(self,'borderpattyn','vertices on Pattyn''s border (for tiling)'))
-		string="%s\n%s"%(string,fielddisplay(self,'borderstokes','vertices on Stokes'' border (for tiling)'))
+		string="%s\n\n%s"%(string,fielddisplay(self,'ismacayealpattyn',"is the macayeal or pattyn approximation used ?"))
+		string="%s\n%s"%(string,fielddisplay(self,'ishutter',"is the shallow ice approximation used ?"))
+		string="%s\n%s"%(string,fielddisplay(self,'isl1l2',"are l1l2 equations used ?"))
+		string="%s\n%s"%(string,fielddisplay(self,'isstokes',"are the Full-Stokes equations used ?"))
+		string="%s\n%s"%(string,fielddisplay(self,'vertex_equation',"flow equation for each vertex"))
+		string="%s\n%s"%(string,fielddisplay(self,'element_equation',"flow equation for each element"))
+		string="%s\n%s"%(string,fielddisplay(self,'bordermacayeal',"vertices on MacAyeal's border (for tiling)"))
+		string="%s\n%s"%(string,fielddisplay(self,'borderpattyn',"vertices on Pattyn's border (for tiling)"))
+		string="%s\n%s"%(string,fielddisplay(self,'borderstokes',"vertices on Stokes' border (for tiling)"))
 		return string
 		#}}}
Index: /issm/trunk/src/m/classes/initialization.py
===================================================================
--- /issm/trunk/src/m/classes/initialization.py	(revision 14066)
+++ /issm/trunk/src/m/classes/initialization.py	(revision 14067)
@@ -70,6 +70,6 @@
 			md = checkfield(md,'initialization.vy','NaN',1,'size',[md.mesh.numberofvertices])
 			#Triangle with zero velocity
-			if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements.astype(int)-1]),axis=1)==0,\
-			                               numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements.astype(int)-1]),axis=1)==0)):
+			if numpy.any(numpy.logical_and(numpy.sum(numpy.abs(md.initialization.vx[md.mesh.elements-1]),axis=1)==0,\
+			                               numpy.sum(numpy.abs(md.initialization.vy[md.mesh.elements-1]),axis=1)==0)):
 				md.checkmessage("at least one triangle has all its vertices with a zero velocity")
 		if ThermalAnalysisEnum() in analyses:
Index: /issm/trunk/src/m/classes/mask.py
===================================================================
--- /issm/trunk/src/m/classes/mask.py	(revision 14066)
+++ /issm/trunk/src/m/classes/mask.py	(revision 14067)
@@ -33,5 +33,5 @@
 		string="%s\n%s"%(string,fielddisplay(self,"elementonfloatingice","element on floating ice flags list"))
 		string="%s\n%s"%(string,fielddisplay(self,"vertexonfloatingice","vertex on floating ice flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice  list"))
+		string="%s\n%s"%(string,fielddisplay(self,"elementongroundedice","element on grounded ice list"))
 		string="%s\n%s"%(string,fielddisplay(self,"vertexongroundedice","vertex on grounded ice flags list"))
 		string="%s\n%s"%(string,fielddisplay(self,"elementonwater","element on water flags list"))
Index: /issm/trunk/src/m/classes/mesh.py
===================================================================
--- /issm/trunk/src/m/classes/mesh.py	(revision 14066)
+++ /issm/trunk/src/m/classes/mesh.py	(revision 14067)
@@ -92,8 +92,8 @@
 		string="%s\n%s"%(string,fielddisplay(self,"vertexonsurface","upper vertices flags list"))
 		string="%s\n%s"%(string,fielddisplay(self,"elementonsurface","upper elements flags list"))
-		string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (NaN for vertex on the upper surface)"))
-		string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (NaN for element on the upper layer)"))
-		string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (NaN for vertex on the lower surface)"))
-		string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (NaN for element on the lower layer"))
+		string="%s\n%s"%(string,fielddisplay(self,"uppervertex","upper vertex list (-1 for vertex on the upper surface)"))
+		string="%s\n%s"%(string,fielddisplay(self,"upperelements","upper element list (-1 for element on the upper layer)"))
+		string="%s\n%s"%(string,fielddisplay(self,"lowervertex","lower vertex list (-1 for vertex on the lower surface)"))
+		string="%s\n%s"%(string,fielddisplay(self,"lowerelements","lower element list (-1 for element on the lower layer)"))
 		string="%s\n%s"%(string,fielddisplay(self,"vertexonboundary","vertices on the boundary of the domain flag list"))
 		string="%s\n%s"%(string,fielddisplay(self,"segments","edges on domain boundary (vertex1 vertex2 element)"))
Index: /issm/trunk/src/m/classes/model/model.py
===================================================================
--- /issm/trunk/src/m/classes/model/model.py	(revision 14066)
+++ /issm/trunk/src/m/classes/model/model.py	(revision 14067)
@@ -2,4 +2,5 @@
 import numpy
 import copy
+import sys
 from mesh import mesh
 from mask import mask
@@ -207,13 +208,13 @@
 		#kick out all elements with 3 dirichlets
 		spc_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
-		spc_node=numpy.unique(md1.mesh.elements[spc_elem,:]).astype(int)-1
+		spc_node=numpy.unique(md1.mesh.elements[spc_elem,:])-1
 		flag=numpy.ones(md1.mesh.numberofvertices)
 		flag[spc_node]=0
-		pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements.astype(int)-1],axis=1)))[0]
+		pos=numpy.nonzero(numpy.logical_not(numpy.sum(flag[md1.mesh.elements-1],axis=1)))[0]
 		flag_elem[pos]=0
 
 		#extracted elements and nodes lists
 		pos_elem=numpy.nonzero(flag_elem)[0]
-		pos_node=numpy.unique(md1.mesh.elements[pos_elem,:]).astype(int)-1
+		pos_node=numpy.unique(md1.mesh.elements[pos_elem,:])-1
 
 		#keep track of some fields
@@ -226,7 +227,7 @@
 
 		#Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
-		Pelem=numpy.zeros(numberofelements1)
+		Pelem=numpy.zeros(numberofelements1,int)
 		Pelem[pos_elem]=numpy.arange(1,numberofelements2+1)
-		Pnode=numpy.zeros(numberofvertices1)
+		Pnode=numpy.zeros(numberofvertices1,int)
 		Pnode[pos_node]=numpy.arange(1,numberofvertices2+1)
 
@@ -234,11 +235,11 @@
 		elements_1=copy.deepcopy(md1.mesh.elements)
 		elements_2=elements_1[pos_elem,:]
-		elements_2[:,0]=Pnode[elements_2[:,0].astype(int)-1]
-		elements_2[:,1]=Pnode[elements_2[:,1].astype(int)-1]
-		elements_2[:,2]=Pnode[elements_2[:,2].astype(int)-1]
+		elements_2[:,0]=Pnode[elements_2[:,0]-1]
+		elements_2[:,1]=Pnode[elements_2[:,1]-1]
+		elements_2[:,2]=Pnode[elements_2[:,2]-1]
 		if md1.mesh.dimension==3:
-			elements_2[:,3]=Pnode[elements_2[:,3].astype(int)-1]
-			elements_2[:,4]=Pnode[elements_2[:,4].astype(int)-1]
-			elements_2[:,5]=Pnode[elements_2[:,5].astype(int)-1]
+			elements_2[:,3]=Pnode[elements_2[:,3]-1]
+			elements_2[:,4]=Pnode[elements_2[:,4]-1]
+			elements_2[:,5]=Pnode[elements_2[:,5]-1]
 
 		#OK, now create the new model!
@@ -291,18 +292,18 @@
 		if md1.mesh.dimension==3:
 			md2.mesh.uppervertex=md1.mesh.uppervertex[pos_node]
-			pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.uppervertex)))[0]
-			md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos].astype(int)-1]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.uppervertex==-1))[0]
+			md2.mesh.uppervertex[pos]=Pnode[md2.mesh.uppervertex[pos]-1]
 
 			md2.mesh.lowervertex=md1.mesh.lowervertex[pos_node]
-			pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowervertex)))[0]
-			md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos].astype(int)-1]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowervertex==-1))[0]
+			md2.mesh.lowervertex[pos]=Pnode[md2.mesh.lowervertex[pos]-1]
 
 			md2.mesh.upperelements=md1.mesh.upperelements[pos_elem]
-			pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.upperelements)))[0]
-			md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos].astype(int)-1]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.upperelements==-1))[0]
+			md2.mesh.upperelements[pos]=Pelem[md2.mesh.upperelements[pos]-1]
 
 			md2.mesh.lowerelements=md1.mesh.lowerelements[pos_elem]
-			pos=numpy.nonzero(numpy.logical_not(numpy.isnan(md2.mesh.lowerelements)))[0]
-			md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos].astype(int)-1]
+			pos=numpy.nonzero(numpy.logical_not(md2.mesh.lowerelements==-1))[0]
+			md2.mesh.lowerelements[pos]=Pelem[md2.mesh.lowerelements[pos]-1]
 
 		#Initial 2d mesh 
@@ -316,7 +317,7 @@
 			md2.mesh.numberofvertices2d=numpy.size(pos_node_2d)
 			md2.mesh.elements2d=md1.mesh.elements2d[pos_elem_2d,:]
-			md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0].astype(int)-1]
-			md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1].astype(int)-1]
-			md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2].astype(int)-1]
+			md2.mesh.elements2d[:,0]=Pnode[md2.mesh.elements2d[:,0]-1]
+			md2.mesh.elements2d[:,1]=Pnode[md2.mesh.elements2d[:,1]-1]
+			md2.mesh.elements2d[:,2]=Pnode[md2.mesh.elements2d[:,2]-1]
 
 			md2.mesh.x2d=md1.mesh.x[pos_node_2d]
@@ -324,11 +325,11 @@
 
 		#Edges
-		if len(numpy.shape(md2.mesh.edges))>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
+		if numpy.ndim(md2.mesh.edges)>1 and numpy.size(md2.mesh.edges,axis=1)>1:    #do not use ~isnan because there are some NaNs...
 			#renumber first two columns
 			pos=numpy.nonzero(md2.mesh.edges[:,3]!=-1)[0]
-			md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0].astype(int)-1]
-			md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1].astype(int)-1]
-			md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2].astype(int)-1]
-			md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3].astype(int)-1]
+			md2.mesh.edges[:  ,0]=Pnode[md2.mesh.edges[:,0]-1]
+			md2.mesh.edges[:  ,1]=Pnode[md2.mesh.edges[:,1]-1]
+			md2.mesh.edges[:  ,2]=Pelem[md2.mesh.edges[:,2]-1]
+			md2.mesh.edges[pos,3]=Pelem[md2.mesh.edges[pos,3]-1]
 			#remove edges when the 2 vertices are not in the domain.
 			md2.mesh.edges=md2.mesh.edges[numpy.nonzero(numpy.logical_and(md2.mesh.edges[:,0],md2.mesh.edges[:,1]))[0],:]
@@ -364,6 +365,6 @@
 			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements,md2.mesh.vertexconnectivity)
 			md2.mesh.segments=contourenvelope(md2)
-			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2)
-			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2,bool)
+			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
 		else:
 			#First do the connectivity for the contourenvelope in 2d
@@ -371,6 +372,6 @@
 			[md2.mesh.elementconnectivity]=ElementConnectivity(md2.mesh.elements2d,md2.mesh.vertexconnectivity)
 			md2.mesh.segments=contourenvelope(md2)
-			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers)
-			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2].astype(int)-1]=1
+			md2.mesh.vertexonboundary=numpy.zeros(numberofvertices2/md2.mesh.numberoflayers,bool)
+			md2.mesh.vertexonboundary[md2.mesh.segments[:,0:2]-1]=True
 			md2.mesh.vertexonboundary=numpy.tile(md2.mesh.vertexonboundary,md2.mesh.numberoflayers)
 			#Then do it for 3d as usual
@@ -381,5 +382,5 @@
 		#Catch the elements that have not been extracted
 		orphans_elem=numpy.nonzero(numpy.logical_not(flag_elem))[0]
-		orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:]).astype(int)-1
+		orphans_node=numpy.unique(md1.mesh.elements[orphans_elem,:])-1
 		#Figure out which node are on the boundary between md2 and md1
 		nodestoflag1=numpy.intersect1d(orphans_node,pos_node)
@@ -443,6 +444,6 @@
 
 		#Keep track of pos_node and pos_elem
-		md2.mesh.extractedvertices=pos_node.astype(float)+1
-		md2.mesh.extractedelements=pos_elem.astype(float)+1
+		md2.mesh.extractedvertices=pos_node+1
+		md2.mesh.extractedelements=pos_elem+1
 
 		return md2
@@ -526,5 +527,5 @@
 
 		#Extrude elements 
-		elements3d=numpy.empty((0,6))
+		elements3d=numpy.empty((0,6),int)
 		for i in xrange(numlayers-1):
 			elements3d=numpy.vstack((elements3d,numpy.hstack((md.mesh.elements+i*md.mesh.numberofvertices,md.mesh.elements+(i+1)*md.mesh.numberofvertices))))    #Create the elements of the 3d mesh for the non extruded part
@@ -532,6 +533,6 @@
 
 		#Keep a trace of lower and upper nodes
-		mesh.lowervertex=float('NaN')*numpy.ones(number_nodes3d)
-		mesh.uppervertex=float('NaN')*numpy.ones(number_nodes3d)
+		mesh.lowervertex=-1*numpy.ones(number_nodes3d,int)
+		mesh.uppervertex=-1*numpy.ones(number_nodes3d,int)
 		mesh.lowervertex[md.mesh.numberofvertices:]=numpy.arange(1,(numlayers-1)*md.mesh.numberofvertices+1)
 		mesh.uppervertex[:(numlayers-1)*md.mesh.numberofvertices]=numpy.arange(md.mesh.numberofvertices+1,number_nodes3d+1)
@@ -540,6 +541,6 @@
 
 		#same for lower and upper elements
-		mesh.lowerelements=float('NaN')*numpy.ones(number_el3d)
-		mesh.upperelements=float('NaN')*numpy.ones(number_el3d)
+		mesh.lowerelements=-1*numpy.ones(number_el3d,int)
+		mesh.upperelements=-1*numpy.ones(number_el3d,int)
 		mesh.lowerelements[md.mesh.numberofelements:]=numpy.arange(1,(numlayers-2)*md.mesh.numberofelements+1)
 		mesh.upperelements[:(numlayers-2)*md.mesh.numberofelements]=numpy.arange(md.mesh.numberofelements+1,(numlayers-1)*md.mesh.numberofelements+1)
@@ -603,13 +604,13 @@
 
 		#bedinfo and surface info
-		md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',1)
-		md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d),'type','element','layer',md.mesh.numberoflayers-1)
-		md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',1)
-		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d),'type','node','layer',md.mesh.numberoflayers)
+		md.mesh.elementonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',1)
+		md.mesh.elementonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofelements2d,bool),'type','element','layer',md.mesh.numberoflayers-1)
+		md.mesh.vertexonbed=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',1)
+		md.mesh.vertexonsurface=project3d(md,'vector',numpy.ones(md.mesh.numberofvertices2d,bool),'type','node','layer',md.mesh.numberoflayers)
 
 		#elementstype
 		if not numpy.any(numpy.isnan(md.flowequation.element_equation)):
 			oldelements_type=md.flowequation.element_equation
-			md.flowequation.element_equation=numpy.zeros(number_el3d)
+			md.flowequation.element_equation=numpy.zeros(number_el3d,int)
 			md.flowequation.element_equation=project3d(md,'vector',oldelements_type,'type','element')
 
@@ -617,5 +618,5 @@
 		if not numpy.any(numpy.isnan(md.flowequation.vertex_equation)):
 			oldvertices_type=md.flowequation.vertex_equation
-			md.flowequation.vertex_equation=numpy.zeros(number_nodes3d)
+			md.flowequation.vertex_equation=numpy.zeros(number_nodes3d,int)
 			md.flowequation.vertex_equation=project3d(md,'vector',oldvertices_type,'type','node')
 
@@ -635,5 +636,5 @@
 		#in 3d, pressureload: [node1 node2 node3 node4 element]
 		pressureload_layer1=numpy.hstack((md.diagnostic.icefront[:,0:2],md.diagnostic.icefront[:,1:2]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,0:1]+md.mesh.numberofvertices2d,md.diagnostic.icefront[:,2:4]))    #Add two columns on the first layer 
-		pressureload=numpy.empty((0,6))
+		pressureload=numpy.empty((0,6),int)
 		for i in xrange(numlayers-1):
 			pressureload=numpy.vstack((pressureload,numpy.hstack((pressureload_layer1[:,0:4]+i*md.mesh.numberofvertices2d,pressureload_layer1[:,4:5]+i*md.mesh.numberofelements2d,pressureload_layer1[:,5:6]))))
@@ -642,9 +643,9 @@
 		#connectivity
 		md.mesh.elementconnectivity=numpy.tile(md.mesh.elementconnectivity,(numlayers-1,1))
-		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=float('NaN')
+		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity==0)]=-sys.maxint-1
 		for i in xrange(1,numlayers-1):
 			md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:] \
 				=md.mesh.elementconnectivity[i*md.mesh.numberofelements2d:(i+1)*md.mesh.numberofelements2d,:]+md.mesh.numberofelements2d
-		md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
+		md.mesh.elementconnectivity[numpy.nonzero(md.mesh.elementconnectivity<0)]=0
 
 		#materials
Index: /issm/trunk/src/m/classes/transient.m
===================================================================
--- /issm/trunk/src/m/classes/transient.m	(revision 14066)
+++ /issm/trunk/src/m/classes/transient.m	(revision 14067)
@@ -46,6 +46,6 @@
 
 			fielddisplay(obj,'isprognostic','indicates if a prognostic solution is used in the transient');
+			fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient');
 			fielddisplay(obj,'isthermal','indicates if a thermal solution is used in the transient');
-			fielddisplay(obj,'isdiagnostic','indicates if a diagnostic solution is used in the transient');
 			fielddisplay(obj,'isgroundingline','indicates if a groundingline migration is used in the transient');
 			fielddisplay(obj,'requested_outputs','list of additional outputs requested');
Index: /issm/trunk/src/m/classes/transient.py
===================================================================
--- /issm/trunk/src/m/classes/transient.py	(revision 14066)
+++ /issm/trunk/src/m/classes/transient.py	(revision 14067)
@@ -16,8 +16,8 @@
 	def __init__(self):
 		# {{{ Properties
-		self.isprognostic      = 0
-		self.isdiagnostic      = 0
-		self.isthermal         = 0
-		self.isgroundingline   = 0
+		self.isprognostic      = False
+		self.isdiagnostic      = False
+		self.isthermal         = False
+		self.isgroundingline   = False
 		self.requested_outputs = float('NaN')
 
@@ -30,6 +30,6 @@
 		string='   transient solution parameters:'
 		string="%s\n%s"%(string,fielddisplay(self,'isprognostic','indicates if a prognostic solution is used in the transient'))
+		string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient'))
 		string="%s\n%s"%(string,fielddisplay(self,'isthermal','indicates if a thermal solution is used in the transient'))
-		string="%s\n%s"%(string,fielddisplay(self,'isdiagnostic','indicates if a diagnostic solution is used in the transient'))
 		string="%s\n%s"%(string,fielddisplay(self,'isgroundingline','indicates if a groundingline migration is used in the transient'))
 		string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','list of additional outputs requested'))
@@ -41,8 +41,8 @@
 		
 		#full analysis: Diagnostic, Prognostic and Thermal but no groundingline migration for now
-		self.isprognostic=1
-		self.isdiagnostic=1
-		self.isthermal=1
-		self.isgroundingline=0
+		self.isprognostic=True
+		self.isdiagnostic=True
+		self.isthermal=True
+		self.isgroundingline=False
 
 		return self
Index: /issm/trunk/src/m/consistency/checkfield.py
===================================================================
--- /issm/trunk/src/m/consistency/checkfield.py	(revision 14066)
+++ /issm/trunk/src/m/consistency/checkfield.py	(revision 14067)
@@ -143,5 +143,5 @@
 	if options.getfieldvalue('forcing',0):
 		if   numpy.size(field,0)==md.mesh.numberofvertices:
-			if len(numpy.shape(field))>1 and not numpy.size(field,1)==1:
+			if numpy.ndim(field)>1 and not numpy.size(field,1)==1:
 				md = md.checkmessage(options.getfieldvalue('message',\
 					"field '%s' should have only one column as there are md.mesh.numberofvertices lines" % fieldname))
Index: /issm/trunk/src/m/exp/expcoarsen.m
===================================================================
--- /issm/trunk/src/m/exp/expcoarsen.m	(revision 14066)
+++ /issm/trunk/src/m/exp/expcoarsen.m	(revision 14067)
@@ -25,5 +25,5 @@
 
 %Get exp oldfile
-[path root ext ver]=fileparts(oldfile);
+[path root ext]=fileparts(oldfile);
 A=expread(oldfile);
 numprofiles=size(A,2);
Index: /issm/trunk/src/m/extrusion/project3d.py
===================================================================
--- /issm/trunk/src/m/extrusion/project3d.py	(revision 14066)
+++ /issm/trunk/src/m/extrusion/project3d.py	(revision 14067)
@@ -38,5 +38,5 @@
 
 	vector1d=False
-	if isinstance(vector2d,numpy.ndarray) and len(vector2d.shape)==1:
+	if isinstance(vector2d,numpy.ndarray) and numpy.ndim(vector2d)==1:
 		vector1d=True
 		vector2d=vector2d.reshape(numpy.size(vector2d),1)
@@ -49,7 +49,7 @@
 		#Initialize 3d vector
 		if   numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d:
-			projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices,  numpy.size(vector2d,axis=1)))
+			projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices,  numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
 		elif numpy.size(vector2d,axis=0)==md.mesh.numberofvertices2d+1:
-			projected_vector=paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))
+			projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofvertices+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
 			projected_vector[-1,:]=vector2d[-1,:]
 			vector2d=vector2d[:-1,:]
@@ -68,7 +68,7 @@
 		#Initialize 3d vector
 		if   numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d:
-			projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements,  numpy.size(vector2d,axis=1)))
+			projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements,  numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
 		elif numpy.size(vector2d,axis=0)==md.mesh.numberofelements2d+1:
-			projected_vector=paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))
+			projected_vector=(paddingvalue*numpy.ones((md.mesh.numberofelements+1,numpy.size(vector2d,axis=1)))).astype(vector2d.dtype)
 			projected_vector[-1,:]=vector2d[-1,:]
 			vector2d=vector2d[:-1,:]
Index: /issm/trunk/src/m/geometry/FlagElements.py
===================================================================
--- /issm/trunk/src/m/geometry/FlagElements.py	(revision 14066)
+++ /issm/trunk/src/m/geometry/FlagElements.py	(revision 14067)
@@ -24,8 +24,8 @@
 	if   isinstance(region,(str,unicode)):
 		if   not region:
-			flag=numpy.zeros(md.mesh.numberofelements,'bool')
+			flag=numpy.zeros(md.mesh.numberofelements,bool)
 			invert=0
 		elif strcmpi(region,'all'):
-			flag=numpy.ones(md.mesh.numberofelements,'bool')
+			flag=numpy.ones(md.mesh.numberofelements,bool)
 			invert=0
 		else:
@@ -43,9 +43,10 @@
 				raise RuntimeError("FlagElements.py calling basinzoom.py is not complete.")
 				xlim,ylim=basinzoom('basin',region)
-				flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0])).astype(float)
-				flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1)
+				flag_nodes=numpy.logical_and(numpy.logical_and(md.mesh.x<xlim[1],md.mesh.x>xlim[0]),numpy.logical_and(md.mesh.y<ylim[1],md.mesh.y>ylim[0]))
+				flag=numpy.prod(flag_nodes[md.mesh.elements],axis=1).astype(bool)
 			else:
 				#ok, flag elements
 				[flag,dum]=ContourToMesh(md.mesh.elements[:,0:3].copy(),md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),region,'element',1)
+				flag=flag.astype(bool)
 
 		if invert:
Index: /issm/trunk/src/m/geometry/GetAreas.py
===================================================================
--- /issm/trunk/src/m/geometry/GetAreas.py	(revision 14066)
+++ /issm/trunk/src/m/geometry/GetAreas.py	(revision 14067)
@@ -33,10 +33,10 @@
 	#initialization
 	areas=numpy.zeros(nels)
-	x1=x[index[:,0].astype(int)-1]
-	x2=x[index[:,1].astype(int)-1]
-	x3=x[index[:,2].astype(int)-1]
-	y1=y[index[:,0].astype(int)-1]
-	y2=y[index[:,1].astype(int)-1]
-	y3=y[index[:,2].astype(int)-1]
+	x1=x[index[:,0]-1]
+	x2=x[index[:,1]-1]
+	x3=x[index[:,2]-1]
+	y1=y[index[:,0]-1]
+	y2=y[index[:,1]-1]
+	y3=y[index[:,2]-1]
 
 	#compute the volume of each element
@@ -46,5 +46,5 @@
 	else:
 		#V=area(triangle)*1/3(z1+z2+z3)
-		thickness=numpy.mean(z[index[:,3:6].astype(int)-1])-numpy.mean(z[index[:,0:3].astype(int)-1])
+		thickness=numpy.mean(z[index[:,3:6]-1])-numpy.mean(z[index[:,0:3]-1])
 		areas=(0.5*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)))*thickness
 
Index: /issm/trunk/src/m/mesh/ComputeHessian.py
===================================================================
--- /issm/trunk/src/m/mesh/ComputeHessian.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/ComputeHessian.py	(revision 14067)
@@ -45,6 +45,6 @@
 
 	#Compute gradient for each element
-	grad_elx=numpy.sum(field[index.astype(int)-1,0]*alpha,axis=1) 
-	grad_ely=numpy.sum(field[index.astype(int)-1,0]*beta,axis=1)
+	grad_elx=numpy.sum(field[index-1,0]*alpha,axis=1) 
+	grad_ely=numpy.sum(field[index-1,0]*beta,axis=1)
 
 	#Compute gradient for each node (average of the elements around)
@@ -55,5 +55,5 @@
 
 	#Compute hessian for each element
-	hessian=numpy.hstack((numpy.sum(gradx[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index.astype(int)-1,0]*beta,axis=1).reshape(-1,1)))
+	hessian=numpy.hstack((numpy.sum(gradx[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*alpha,axis=1).reshape(-1,1),numpy.sum(grady[index-1,0]*beta,axis=1).reshape(-1,1)))
 
 	if strcmpi(type,'node'):
Index: /issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py
===================================================================
--- /issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/GetNodalFunctionsCoeff.py	(revision 14067)
@@ -39,10 +39,10 @@
 
 	#compute nodal functions coefficients N(x,y)=alpha x + beta y +gamma
-	x1=x[index[:,0].astype(int)-1]
-	x2=x[index[:,1].astype(int)-1]
-	x3=x[index[:,2].astype(int)-1]
-	y1=y[index[:,0].astype(int)-1]
-	y2=y[index[:,1].astype(int)-1]
-	y3=y[index[:,2].astype(int)-1]
+	x1=x[index[:,0]-1]
+	x2=x[index[:,1]-1]
+	x3=x[index[:,2]-1]
+	y1=y[index[:,0]-1]
+	y2=y[index[:,1]-1]
+	y3=y[index[:,2]-1]
 	invdet=1./(x1*(y2-y3)-x2*(y1-y3)+x3*(y1-y2))
 
Index: /issm/trunk/src/m/mesh/bamg.py
===================================================================
--- /issm/trunk/src/m/mesh/bamg.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/bamg.py	(revision 14067)
@@ -320,8 +320,8 @@
 	md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
 	md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
-	md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
-	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
-	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
+	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+	md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
 
 	#Fill in rest of fields:
@@ -331,13 +331,14 @@
 	md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
 	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
-	md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
-	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
 	md.mesh.elementconnectivity=md.private.bamg['mesh'].ElementConnectivity
 	md.mesh.elementconnectivity[numpy.nonzero(numpy.isnan(md.mesh.elementconnectivity))]=0
+	md.mesh.elementconnectivity=md.mesh.elementconnectivity.astype(int)
 
 	#Check for orphan
Index: /issm/trunk/src/m/mesh/meshconvert.py
===================================================================
--- /issm/trunk/src/m/mesh/meshconvert.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/meshconvert.py	(revision 14067)
@@ -35,8 +35,8 @@
 	md.mesh.x=bamgmesh_out['Vertices'][:,0].copy()
 	md.mesh.y=bamgmesh_out['Vertices'][:,1].copy()
-	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].copy()
-	md.mesh.edges=bamgmesh_out['IssmEdges'].copy()
-	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].copy()
-	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].copy()
+	md.mesh.elements=bamgmesh_out['Triangles'][:,0:3].astype(int)
+	md.mesh.edges=bamgmesh_out['IssmEdges'].astype(int)
+	md.mesh.segments=bamgmesh_out['IssmSegments'][:,0:3].astype(int)
+	md.mesh.segmentmarkers=bamgmesh_out['IssmSegments'][:,3].astype(int)
 
 	#Fill in rest of fields:
@@ -46,11 +46,11 @@
 	md.mesh.numberofedges=numpy.size(md.mesh.edges,axis=0)
 	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
-	md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
-	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mask.vertexonwater=numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.vertexonboundary=numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
 
 	return md
Index: /issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py
===================================================================
--- /issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/rifts/meshprocessoutsiderifts.py	(revision 14067)
@@ -41,5 +41,5 @@
 			B=node_connected_to_tip
 
-			elements=numpy.empty(0)
+			elements=numpy.empty(0,int)
 
 			while flags(B):    #as long as B does not belong to the domain outline, keep looking.
@@ -62,8 +62,8 @@
 		
 			#replace tip in elements
-			newelements=md.mesh.elements[elements.astype(int)-1,:]
+			newelements=md.mesh.elements[elements-1,:]
 			pos=numpy.nonzero(newelements==tip)
 			newelements[pos]=num
-			md.mesh.elements[elements.astype(int)-1,:]=newelements
+			md.mesh.elements[elements-1,:]=newelements
 			rift.tips=numpy.concatenate((rift.tips,num))
 
@@ -81,12 +81,12 @@
 	md.mesh.numberofvertices=numpy.size(md.mesh.x)
 	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
+	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
+	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
 	md.rifts.numrifts=length(md.rifts.riftstruct)
-	md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements)
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
+	md.flowequation.element_equation=3*numpy.ones(md.mesh.numberofelements,int)
+	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
 
 	return md
Index: /issm/trunk/src/m/mesh/rifts/meshprocessrifts.py
===================================================================
--- /issm/trunk/src/m/mesh/rifts/meshprocessrifts.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/rifts/meshprocessrifts.py	(revision 14067)
@@ -23,6 +23,9 @@
 	#Call MEX file
 	[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers)
+	md.mesh.elements=md.mesh.elements.astype(int)
 	md.mesh.x=md.mesh.x.reshape(-1)
 	md.mesh.y=md.mesh.y.reshape(-1)
+	md.mesh.segments=md.mesh.segments.astype(int)
+	md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
 	if not isinstance(md.rifts.riftstruct,list) or not md.rifts.riftstruct:
 		raise RuntimeError("TriMeshProcessRifts did not find any rift")
@@ -33,10 +36,10 @@
 	md.mesh.numberofvertices=numpy.size(md.mesh.x)
 	md.mesh.z=numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x))
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1
-	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices)
-	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements)
-	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements)
+	md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x),bool)
+	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
+	md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements,bool)
 
 	#get coordinates of rift tips
Index: /issm/trunk/src/m/mesh/triangle.py
===================================================================
--- /issm/trunk/src/m/mesh/triangle.py	(revision 14066)
+++ /issm/trunk/src/m/mesh/triangle.py	(revision 14067)
@@ -44,4 +44,7 @@
 	#Mesh using TriMesh
 	[md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers]=TriMesh(domainname,riftname,area)
+	md.mesh.elements=md.mesh.elements.astype(int)
+	md.mesh.segments=md.mesh.segments.astype(int)
+	md.mesh.segmentmarkers=md.mesh.segmentmarkers.astype(int)
 
 	#Fill in rest of fields:
@@ -49,13 +52,13 @@
 	md.mesh.numberofvertices = numpy.size(md.mesh.x)
 	md.mesh.z = numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices)
-	md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1] = 1.
-	md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices)
-	md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices)
-	md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements)
-	md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements)
+	md.mesh.vertexonboundary = numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1] = True
+	md.mesh.vertexonbed = numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.vertexonsurface = numpy.ones(md.mesh.numberofvertices,bool)
+	md.mesh.elementonbed = numpy.ones(md.mesh.numberofelements,bool)
+	md.mesh.elementonsurface = numpy.ones(md.mesh.numberofelements,bool)
 
 	#Now, build the connectivity tables for this mesh.
-	[md.mesh.vertexconnectivity]= NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
+	[md.mesh.vertexconnectivity] = NodeConnectivity(md.mesh.elements, md.mesh.numberofvertices)
 	[md.mesh.elementconnectivity] = ElementConnectivity(md.mesh.elements, md.mesh.vertexconnectivity)
 
Index: /issm/trunk/src/m/parameterization/contourenvelope.py
===================================================================
--- /issm/trunk/src/m/parameterization/contourenvelope.py	(revision 14066)
+++ /issm/trunk/src/m/parameterization/contourenvelope.py	(revision 14067)
@@ -76,5 +76,5 @@
 			pos=numpy.nonzero(flags)
 			elemin[pos]=1
-			nodein[mesh.elements[pos,:].astype(int)-1]=1
+			nodein[mesh.elements[pos,:]-1]=1
 
 			#modify element connectivity
@@ -88,14 +88,14 @@
 	if len(args)==1:
 		flag[numpy.nonzero(flag)]=elemin[flag[numpy.nonzero(flag)]]
-	elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0).astype(float)
+	elementonboundary=numpy.logical_and(numpy.prod(flag,axis=1)==0,numpy.sum(flag,axis=1)>0)
 
 	#Find segments on boundary
 	pos=numpy.nonzero(elementonboundary)[0]
 	num_segments=numpy.size(pos)
-	segments=numpy.zeros((num_segments*3,3))
+	segments=numpy.zeros((num_segments*3,3),int)
 	count=0
 
 	for el1 in pos:
-		els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]].astype(int)-1
+		els2=mesh.elementconnectivity[el1,numpy.nonzero(mesh.elementconnectivity[el1,:])[0]]-1
 		if numpy.size(els2)>1:
 			flag=numpy.intersect1d(mesh.elements[els2[0],:],mesh.elements[els2[1],:])
Index: /issm/trunk/src/m/parameterization/setflowequation.py
===================================================================
--- /issm/trunk/src/m/parameterization/setflowequation.py	(revision 14066)
+++ /issm/trunk/src/m/parameterization/setflowequation.py	(revision 14067)
@@ -50,9 +50,9 @@
 	#Flag the elements that have not been flagged as filltype
 	if   strcmpi(filltype,'hutter'):
-		hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=1
+		hutterflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(macayealflag,pattynflag)))]=True
 	elif strcmpi(filltype,'macayeal'):
-		macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=1
+		macayealflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(pattynflag,stokesflag))))]=True
 	elif strcmpi(filltype,'pattyn'):
-		pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=1
+		pattynflag[numpy.nonzero(numpy.logical_not(numpy.logical_or(hutterflag,numpy.logical_or(macayealflag,stokesflag))))]=True
 
 	#check that each element has at least one flag
@@ -63,7 +63,7 @@
 	if any(hutterflag+macayealflag+l1l2flag+pattynflag+stokesflag>1):
 		print "setflowequation warning message: some elements have several types, higher order type is used for them"
-		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=0
-		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=0
-		macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=0
+		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,macayealflag))]=False
+		hutterflag[numpy.nonzero(numpy.logical_and(hutterflag,pattynflag))]=False
+		macayealflag[numpy.nonzero(numpy.logical_and(macayealflag,pattynflag))]=False
 
 	#Check that no pattyn or stokes for 2d mesh
@@ -77,14 +77,14 @@
 
 	#Initialize node fields
-	nodeonhutter=numpy.zeros(md.mesh.numberofvertices)
-	nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:].astype(int)-1]=1
-	nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices)
-	nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
-	nodeonl1l2=numpy.zeros(md.mesh.numberofvertices)
-	nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:].astype(int)-1]=1
-	nodeonpattyn=numpy.zeros(md.mesh.numberofvertices)
-	nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
-	nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
-	noneflag=numpy.zeros(md.mesh.numberofelements)
+	nodeonhutter=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonhutter[md.mesh.elements[numpy.nonzero(hutterflag),:]-1]=True
+	nodeonmacayeal=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+	nodeonl1l2=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonl1l2[md.mesh.elements[numpy.nonzero(l1l2flag),:]-1]=True
+	nodeonpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+	nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+	noneflag=numpy.zeros(md.mesh.numberofelements,bool)
 
 	#First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
@@ -96,27 +96,27 @@
 		                              numpy.logical_and(nodeonpattyn,nodeonstokes).reshape(-1,1)).astype(int)    #find all the nodes on the boundary of the domain without icefront
 #		fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
-		fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements.astype(int)-1],axis=1)==6).astype(int)    #find all the nodes on the boundary of the domain without icefront
-		stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=0
-		nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
+		fullspcelems=(numpy.sum(fullspcnodes[md.mesh.elements-1],axis=1)==6).astype(int)    #find all the nodes on the boundary of the domain without icefront
+		stokesflag[numpy.nonzero(fullspcelems.reshape(-1))]=False
+		nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
 
 	#Then complete with NoneApproximation or the other model used if there is no stokes
 	if any(stokesflag): 
 		if   any(pattynflag):    #fill with pattyn
-			pattynflag[numpy.logical_not(stokesflag)]=1
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
+			pattynflag[numpy.logical_not(stokesflag)]=True
+			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
 		elif any(macayealflag):    #fill with macayeal
-			macayealflag[numpy.logical_not(stokesflag)]=1
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
+			macayealflag[numpy.logical_not(stokesflag)]=True
+			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
 		else:    #fill with none 
-			noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=1
+			noneflag[numpy.nonzero(numpy.logical_not(stokesflag))]=True
 
 	#Now take care of the coupling between MacAyeal and Pattyn
 	md.diagnostic.vertex_pairing=numpy.array([])
-	nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices)
-	nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices)
-	nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices)
-	macayealpattynflag=numpy.zeros(md.mesh.numberofelements)
-	macayealstokesflag=numpy.zeros(md.mesh.numberofelements)
-	pattynstokesflag=numpy.zeros(md.mesh.numberofelements)
+	nodeonmacayealpattyn=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonpattynstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+	nodeonmacayealstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+	macayealpattynflag=numpy.zeros(md.mesh.numberofelements,bool)
+	macayealstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
+	pattynstokesflag=numpy.zeros(md.mesh.numberofelements,bool)
 	if   strcmpi(coupling_method,'penalties'):
 		#Create the border nodes between Pattyn and MacAyeal and extrude them
@@ -135,103 +135,103 @@
 		if   any(macayealflag) and any(pattynflag):    #coupling macayeal pattyn
 			#Find node at the border
-			nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=1
+			nodeonmacayealpattyn[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonpattyn))]=True
 			#Macayeal elements in contact with this layer become MacAyealPattyn elements
-			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealpattyn)[0])
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealpattyn)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
-			macayealflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealpattynelements
-			macayealpattynflag[numpy.nonzero(commonelements)]=1
-			nodeonmacayeal[:]=0
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
+			commonelements[numpy.nonzero(pattynflag)]=False    #only one layer: the elements previously in macayeal
+			macayealflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealpattynelements
+			macayealpattynflag[numpy.nonzero(commonelements)]=True
+			nodeonmacayeal[:]=False
+			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
 			pos=numpy.nonzero(macayealpattynflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			macayealflag[pos[pos1]]=1
-			macayealpattynflag[pos[pos1]]=0
+			macayealflag[pos[pos1]]=True
+			macayealpattynflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			pattynflag[pos[pos2]]=1
-			macayealpattynflag[pos[pos2]]=0
+			pattynflag[pos[pos2]]=True
+			macayealpattynflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonmacayeal[:]=0
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
-			nodeonpattyn[:]=0
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
-			nodeonmacayealpattyn[:]=0
-			nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:].astype(int)-1]=1
+			nodeonmacayeal[:]=False
+			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+			nodeonpattyn[:]=False
+			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+			nodeonmacayealpattyn[:]=False
+			nodeonmacayealpattyn[md.mesh.elements[numpy.nonzero(macayealpattynflag),:]-1]=True
 
 		elif any(pattynflag) and any(stokesflag):    #coupling pattyn stokes
 			#Find node at the border
-			nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=1
+			nodeonpattynstokes[numpy.nonzero(numpy.logical_and(nodeonpattyn,nodeonstokes))]=True
 			#Stokes elements in contact with this layer become PattynStokes elements
-			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonpattynstokes)[0])
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonpattynstokes)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(pattynflag)]=0    #only one layer: the elements previously in macayeal
-			stokesflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealpattynelements
-			pattynstokesflag[numpy.nonzero(commonelements)]=1
-			nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
+			commonelements[numpy.nonzero(pattynflag)]=False    #only one layer: the elements previously in macayeal
+			stokesflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealpattynelements
+			pattynstokesflag[numpy.nonzero(commonelements)]=True
+			nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
 			pos=numpy.nonzero(pattynstokesflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonpattyn[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			stokesflag[pos[pos1]]=1
-			pattynstokesflag[pos[pos1]]=0
+			stokesflag[pos[pos1]]=True
+			pattynstokesflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			pattynflag[pos[pos2]]=1
-			pattynstokesflag[pos[pos2]]=0
+			pattynflag[pos[pos2]]=True
+			pattynstokesflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonstokes[:]=0
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
-			nodeonpattyn[:]=0
-			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:].astype(int)-1]=1
-			nodeonpattynstokes[:]=0
-			nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:].astype(int)-1]=1
+			nodeonstokes[:]=False
+			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+			nodeonpattyn[:]=False
+			nodeonpattyn[md.mesh.elements[numpy.nonzero(pattynflag),:]-1]=True
+			nodeonpattynstokes[:]=False
+			nodeonpattynstokes[md.mesh.elements[numpy.nonzero(pattynstokesflag),:]-1]=True
 
 		elif any(stokesflag) and any(macayealflag):
 			#Find node at the border
-			nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=1
+			nodeonmacayealstokes[numpy.nonzero(numpy.logical_and(nodeonmacayeal,nodeonstokes))]=True
 			#Stokes elements in contact with this layer become MacAyealStokes elements
-			matrixelements=ismember(md.mesh.elements.astype(int)-1,numpy.nonzero(nodeonmacayealstokes)[0])
+			matrixelements=ismember(md.mesh.elements-1,numpy.nonzero(nodeonmacayealstokes)[0])
 			commonelements=numpy.sum(matrixelements,axis=1)!=0
-			commonelements[numpy.nonzero(macayealflag)]=0    #only one layer: the elements previously in macayeal
-			stokesflag[numpy.nonzero(commonelements)]=0    #these elements are now macayealmacayealelements
-			macayealstokesflag[numpy.nonzero(commonelements)]=1
-			nodeonstokes=numpy.zeros(md.mesh.numberofvertices)
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
+			commonelements[numpy.nonzero(macayealflag)]=False    #only one layer: the elements previously in macayeal
+			stokesflag[numpy.nonzero(commonelements)]=False    #these elements are now macayealmacayealelements
+			macayealstokesflag[numpy.nonzero(commonelements)]=True
+			nodeonstokes=numpy.zeros(md.mesh.numberofvertices,bool)
+			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
 
 			#rule out elements that don't touch the 2 boundaries
 			pos=numpy.nonzero(macayealstokesflag)[0]
 			elist=numpy.zeros(numpy.size(pos),dtype=int)
-			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:].astype(int)-1],axis=1).astype(bool)
-			elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:].astype(int)-1]  ,axis=1).astype(bool)
+			elist = elist + numpy.sum(nodeonmacayeal[md.mesh.elements[pos,:]-1],axis=1).astype(bool)
+			elist = elist - numpy.sum(nodeonstokes[md.mesh.elements[pos,:]-1]  ,axis=1).astype(bool)
 			pos1=numpy.nonzero(elist==1)[0]
-			macayealflag[pos[pos1]]=1
-			macayealstokesflag[pos[pos1]]=0
+			macayealflag[pos[pos1]]=True
+			macayealstokesflag[pos[pos1]]=False
 			pos2=numpy.nonzero(elist==-1)[0]
-			stokesflag[pos[pos2]]=1
-			macayealstokesflag[pos[pos2]]=0
+			stokesflag[pos[pos2]]=True
+			macayealstokesflag[pos[pos2]]=False
 
 			#Recompute nodes associated to these elements
-			nodeonmacayeal[:]=0
-			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:].astype(int)-1]=1
-			nodeonstokes[:]=0
-			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:].astype(int)-1]=1
-			nodeonmacayealstokes[:]=0
-			nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:].astype(int)-1]=1
+			nodeonmacayeal[:]=False
+			nodeonmacayeal[md.mesh.elements[numpy.nonzero(macayealflag),:]-1]=True
+			nodeonstokes[:]=False
+			nodeonstokes[md.mesh.elements[numpy.nonzero(stokesflag),:]-1]=True
+			nodeonmacayealstokes[:]=False
+			nodeonmacayealstokes[md.mesh.elements[numpy.nonzero(macayealstokesflag),:]-1]=True
 
 		elif any(stokesflag) and any(hutterflag):
 			raise TypeError("type of coupling not supported yet")
 
-	#Create MacaAyealPattynApproximation where needed
-	md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements)
+	#Create MacAyealPattynApproximation where needed
+	md.flowequation.element_equation=numpy.zeros(md.mesh.numberofelements,int)
 	md.flowequation.element_equation[numpy.nonzero(noneflag)]=0
 	md.flowequation.element_equation[numpy.nonzero(hutterflag)]=1
@@ -250,5 +250,5 @@
 
 	#Create vertices_type
-	md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices)
+	md.flowequation.vertex_equation=numpy.zeros(md.mesh.numberofvertices,int)
 	pos=numpy.nonzero(nodeonmacayeal)
 	md.flowequation.vertex_equation[pos]=2
@@ -273,8 +273,8 @@
 
 	#figure out solution types
-	md.flowequation.ishutter=float(any(md.flowequation.element_equation==1))
-	md.flowequation.ismacayealpattyn=float(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3)))
-	md.flowequation.isl1l2=float(any(md.flowequation.element_equation==8))
-	md.flowequation.isstokes=float(any(md.flowequation.element_equation==4))
+	md.flowequation.ishutter=any(md.flowequation.element_equation==1)
+	md.flowequation.ismacayealpattyn=bool(numpy.any(numpy.logical_or(md.flowequation.element_equation==2,md.flowequation.element_equation==3)))
+	md.flowequation.isl1l2=any(md.flowequation.element_equation==8)
+	md.flowequation.isstokes=any(md.flowequation.element_equation==4)
 
 	return md
Index: /issm/trunk/src/m/parameterization/setmask.py
===================================================================
--- /issm/trunk/src/m/parameterization/setmask.py	(revision 14066)
+++ /issm/trunk/src/m/parameterization/setmask.py	(revision 14067)
@@ -37,15 +37,15 @@
 	vertexonfloatingice = numpy.zeros(md.mesh.numberofvertices,'bool')
 	vertexongroundedice = numpy.zeros(md.mesh.numberofvertices,'bool')
-	vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:].astype(int)-1]=True
+	vertexongroundedice[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=True
 	vertexonfloatingice[numpy.nonzero(numpy.logical_not(vertexongroundedice))]=True
 	#}}}
 
 	#Return: 
-	md.mask.elementonfloatingice = elementonfloatingice.astype(float)
-	md.mask.vertexonfloatingice = vertexonfloatingice.astype(float)
-	md.mask.elementongroundedice = elementongroundedice.astype(float)
-	md.mask.vertexongroundedice = vertexongroundedice.astype(float)
-	md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices)
-	md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements)
+	md.mask.elementonfloatingice = elementonfloatingice
+	md.mask.vertexonfloatingice = vertexonfloatingice
+	md.mask.elementongroundedice = elementongroundedice
+	md.mask.vertexongroundedice = vertexongroundedice
+	md.mask.vertexonwater = numpy.zeros(md.mesh.numberofvertices,bool)
+	md.mask.elementonwater = numpy.zeros(md.mesh.numberofelements,bool)
 
 	return md
Index: /issm/trunk/src/m/plot/applyoptions.m
===================================================================
--- /issm/trunk/src/m/plot/applyoptions.m	(revision 14066)
+++ /issm/trunk/src/m/plot/applyoptions.m	(revision 14067)
@@ -205,10 +205,17 @@
 			if exist(options,'log')
 				logval= getfieldvalue(options,'log');
-				for i= 1:size(tick_vals,1)
+				for i= 1:numel(tick_vals)
 					tick_vals(i)= logval^(tick_vals(i));
 				end
-			elseif size(tick_vals,1) == 3
+			elseif numel(tick_vals) == 3
 				tick_vals=[tick_vals(1); mean(tick_vals(1:2)); tick_vals(2); ...
 					mean(tick_vals(2:3)); tick_vals(3)];
+				set(c,'YTick',tick_vals);
+			end
+		else
+			if exist(options,'log')
+				logvalue=getfieldvalue(options,'log');
+				set(c,'YTick',log(tick_vals)./log(logvalue));
+			else
 				set(c,'YTick',tick_vals);
 			end
Index: /issm/trunk/src/m/plot/colormaps/getcolormap.m
===================================================================
--- /issm/trunk/src/m/plot/colormaps/getcolormap.m	(revision 14066)
+++ /issm/trunk/src/m/plot/colormaps/getcolormap.m	(revision 14067)
@@ -24,5 +24,5 @@
 	map = map(128:end,:);
 elseif strcmpi(map,'redblue'),
-	map = hsv;
+	map = hsv(256);
 	map = rgb2hsv(map);
 	map(:,2)       = max(min( abs(map(:,1)-0.5)/0.5 ,1),0);
Index: /issm/trunk/src/m/plot/radarpower.m
===================================================================
--- /issm/trunk/src/m/plot/radarpower.m	(revision 14066)
+++ /issm/trunk/src/m/plot/radarpower.m	(revision 14067)
@@ -20,4 +20,5 @@
 end
 
+geotiff_name = getfieldvalue(options,'geotiff_name','');
 highres = getfieldvalue(options,'highres',0);
 xlim    = getfieldvalue(options,'xlim',[min(md.mesh.x) max(md.mesh.x)]);
@@ -64,14 +65,16 @@
 		%md.radaroverlay.x=(xmin:(xmax-xmin)/(size(md.radaroverlay.pwr,2)-1):xmax);
 		%md.radaroverlay.y=(ymin:(ymax-ymin)/(size(md.radaroverlay.pwr,1)-1):ymax);
-		if highres,
-			if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']),
-				error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']);
+		if ~exist(options,'geotiff_name'),
+			if highres,
+				if ~exist(['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif']),
+					error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif not found.']);
+				end
+				geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
+			else
+				if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),
+					error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);
+				end
+				geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
 			end
-			geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog100_r2_hp1.tif'];
-		else
-			if ~exist(['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif']),
-				error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif not found.']);
-			end
-			geotiff_name=['/u/astrid-r1b/ModelData/MOG/mog500_r2_hp1.tif'];
 		end
 
@@ -91,14 +94,16 @@
 
 	elseif strcmpi(md.mesh.hemisphere,'s'),
-		if highres,
-			if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
-				error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
+		if ~exist(options,'geotiff_name'),
+			if highres,
+				if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif']),
+					error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif not found.']);
+				end
+				geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
+			else
+				if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
+					error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
+				end
+				geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
 			end
-			geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_200m.tif'];
-		else
-			if ~exist(['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif']),
-				error(['radarpower error message: file ' '/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif not found.']);
-			end
-			geotiff_name=['/u/astrid-r1b/ModelData/MosaicTiffRsat/amm125m_v2_1km.tif'];
 		end
 
Index: /issm/trunk/src/m/solve/WriteData.py
===================================================================
--- /issm/trunk/src/m/solve/WriteData.py	(revision 14066)
+++ /issm/trunk/src/m/solve/WriteData.py	(revision 14067)
@@ -105,5 +105,5 @@
 		elif isinstance(data,(list,tuple)):
 			data=numpy.array(data).reshape(-1,1)
-		if len(data.shape) == 1:
+		if numpy.ndim(data) == 1:
 			if numpy.size(data):
 				data=data.reshape(numpy.size(data),1)
@@ -138,5 +138,5 @@
 		elif isinstance(data,(list,tuple)):
 			data=numpy.array(data).reshape(-1,1)
-		if len(data.shape) == 1:
+		if numpy.ndim(data) == 1:
 			if numpy.size(data):
 				data=data.reshape(numpy.size(data),1)
@@ -171,5 +171,5 @@
 		elif isinstance(data,(list,tuple)):
 			data=numpy.array(data).reshape(-1,1)
-		if len(data.shape) == 1:
+		if numpy.ndim(data) == 1:
 			if numpy.size(data):
 				data=data.reshape(numpy.size(data),1)
@@ -207,5 +207,5 @@
 			elif isinstance(matrix,(list,tuple)):
 				matrix=numpy.array(matrix).reshape(-1,1)
-			if len(matrix.shape) == 1:
+			if numpy.ndim(matrix) == 1:
 				if numpy.size(matrix):
 					matrix=matrix.reshape(numpy.size(matrix),1)
@@ -231,5 +231,5 @@
 			elif isinstance(matrix,(list,tuple)):
 				matrix=numpy.array(matrix).reshape(-1,1)
-			if len(matrix.shape) == 1:
+			if numpy.ndim(matrix) == 1:
 				matrix=matrix.reshape(numpy.size(matrix),1)
 
Index: /issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp
===================================================================
--- /issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/ElementConnectivity/ElementConnectivity.cpp	(revision 14067)
@@ -13,11 +13,11 @@
 
 	/*inputs: */
-	double* elements=NULL;
-	double* nodeconnectivity=NULL;
-	int     nel,nods;
-	int     width;
+	int* elements=NULL;
+	int* nodeconnectivity=NULL;
+	int  nels,nods;
+	int  width;
 
 	/*outputs: */
-	double* elementconnectivity=NULL;
+	int* elementconnectivity=NULL;
 
 	/*Boot module: */
@@ -28,12 +28,12 @@
         
 	/*Input datasets: */
-	FetchData(&elements,&nel,NULL,ELEMENTS);
+	FetchData(&elements,&nels,NULL,ELEMENTS);
 	FetchData(&nodeconnectivity,&nods,&width,NODECONNECTIVITY);
 
 	/*!Generate internal degree of freedom numbers: */
-	ElementConnectivityx(&elementconnectivity, elements,nel, nodeconnectivity, nods, width);
+	ElementConnectivityx(&elementconnectivity,elements,nels,nodeconnectivity,nods,width);
 
 	/*write output datasets: */
-	WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nel,3);
+	WriteData(ELEMENTCONNECTIVITY,elementconnectivity,nels,3);
 
 	/*end module: */
Index: /issm/trunk/src/wrappers/Kriging/Kriging.cpp
===================================================================
--- /issm/trunk/src/wrappers/Kriging/Kriging.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/Kriging/Kriging.cpp	(revision 14067)
@@ -5,4 +5,8 @@
 
 void KrigingUsage(void){/*{{{*/
+	int num=1;
+#ifdef _MULTITHREADING_
+	num=_NUMTHREADS_;
+#endif
 	_pprintLine_("");
 	_pprintLine_("   usage: predictions=" << __FUNCT__ << "(x,y,observations,x_interp,y_interp,'options');");
@@ -21,4 +25,5 @@
 	_pprintLine_("      -'mintrimming':  minimum trimming value (default is +1.e+21)");
 	_pprintLine_("      -'minspacing':   minimum distance between observation (default is 0.01)");
+	_pprintLine_("      -'numthreads':   number of threads, default is "<<num);
 	_pprintLine_("");
 }/*}}}*/
Index: /issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp
===================================================================
--- /issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/NodeConnectivity/NodeConnectivity.cpp	(revision 14067)
@@ -13,11 +13,11 @@
 
 	/*inputs: */
-	double* elements=NULL;
-	int     nel;
-	int     nods;
+	int* elements=NULL;
+	int  nels;
+	int  nods;
 
 	/*outputs: */
-	double* connectivity=NULL;
-	int     width;
+	int* connectivity=NULL;
+	int  width;
 
 	/*Boot module: */
@@ -28,9 +28,9 @@
         
 	/*Input datasets: */
-	FetchData(&elements,&nel,NULL,ELEMENTS);
+	FetchData(&elements,&nels,NULL,ELEMENTS);
 	FetchData(&nods,NUMNODES);
 
 	/*!Generate internal degree of freedom numbers: */
-	NodeConnectivityx(&connectivity, &width,elements,nel, nods);
+	NodeConnectivityx(&connectivity,&width,elements,nels,nods);
 
 	/*write output datasets: */
Index: /issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp
===================================================================
--- /issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/matlab/io/FetchMatlabData.cpp	(revision 14067)
@@ -264,25 +264,17 @@
 
 	double* outvector=NULL;
-	int outvector_rows;
-
-	if(mxIsEmpty(dataref)){
-		/*Nothing to pick up. Just initialize matrix pointer to NULL: */
-		outvector_rows=0;
-		outvector=NULL;
-	}
-	else if (mxIsClass(dataref,"double") ){
-
-		/*Convert matlab vector to double*  vector: */
-		MatlabVectorToDoubleVector(&outvector,&outvector_rows,dataref);
-
-	}
-	else{
-		/*This is an error: we don't have the correct input!: */
-		_error_("Input parameter of class " << mxGetClassName(dataref) << " not supported yet");
+	int M,N;
+
+	/*Use Fetch matrix*/
+	FetchData(&outvector,&M,&N,dataref) ;
+
+	/*Check that it is a vector*/
+	if(M>0){
+		if(N!=1) _error_("input vector of size " << M << "x" << N << " should have only one column");
 	}
 
 	/*Assign output pointers:*/
 	*pvector=outvector;
-	if (pM)*pM=outvector_rows;
+	if(pM)*pM=M;
 }
 /*}}}*/
Index: /issm/trunk/src/wrappers/python/io/FetchPythonData.cpp
===================================================================
--- /issm/trunk/src/wrappers/python/io/FetchPythonData.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/python/io/FetchPythonData.cpp	(revision 14067)
@@ -67,5 +67,10 @@
 	npy_intp*  dims=NULL;
 
-	/*retrive dimensions: */
+	/*intermediary:*/
+	long* lmatrix=NULL;
+	bool* bmatrix=NULL;
+	int i;
+
+	/*retrieve dimensions: */
 	ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
 	if(ndim!=2)_error_("expecting an MxN matrix in input!");
@@ -74,10 +79,33 @@
 
 	if (M && N) {
-		/*retrieve internal value: */
-		dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
-
-		/*copy matrix: */
-		matrix=xNew<double>(M*N);
-		memcpy(matrix,dmatrix,(M*N)*sizeof(double));
+		if      (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) {
+			/*retrieve internal value: */
+			dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*copy matrix: */
+			matrix=xNew<double>(M*N);
+			memcpy(matrix,dmatrix,(M*N)*sizeof(double));
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) {
+			/*retrieve internal value: */
+			lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*transform into double matrix: */
+			matrix=xNew<double>(M*N);
+			for(i=0;i<M*N;i++)matrix[i]=(double)lmatrix[i];
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) {
+			/*retrieve internal value: */
+			bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*transform into double matrix: */
+			matrix=xNew<double>(M*N);
+			for(i=0;i<M*N;i++)matrix[i]=(double)bmatrix[i];
+		}
+
+		else
+			_error_("unrecognized matrix type in input!");
 	}
 	else
@@ -94,14 +122,16 @@
 
 	/*output: */
-	double* dmatrix=NULL;
 	int* matrix=NULL;
 	int M,N;
-
-	/*intermediary:*/
-	int i;
 	int ndim;
 	npy_intp*  dims=NULL;
 
-	/*retrive dimensions: */
+	/*intermediary:*/
+	double* dmatrix=NULL;
+	long* lmatrix=NULL;
+	bool* bmatrix=NULL;
+	int i;
+
+	/*retrieve dimensions: */
 	ndim=PyArray_NDIM((const PyArrayObject*)py_matrix);
 	if(ndim!=2)_error_("expecting an MxN matrix in input!");
@@ -110,10 +140,33 @@
 
 	if (M && N) {
-		/*retrieve internal value: */
-		dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
-
-		/*transform into integer matrix: */
-		matrix=xNew<int>(M*N);
-		for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i];
+		if      (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_DOUBLE) {
+			/*retrieve internal value: */
+			dmatrix=(double*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*transform into integer matrix: */
+			matrix=xNew<int>(M*N);
+			for(i=0;i<M*N;i++)matrix[i]=(int)dmatrix[i];
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_INT64) {
+			/*retrieve internal value: */
+			lmatrix=(long*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*transform into integer matrix: */
+			matrix=xNew<int>(M*N);
+			for(i=0;i<M*N;i++)matrix[i]=(int)lmatrix[i];
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_matrix) == NPY_BOOL) {
+			/*retrieve internal value: */
+			bmatrix=(bool*)PyArray_DATA((PyArrayObject*)py_matrix);
+
+			/*transform into integer matrix: */
+			matrix=xNew<int>(M*N);
+			for(i=0;i<M*N;i++)matrix[i]=(int)bmatrix[i];
+		}
+
+		else
+			_error_("unrecognized matrix type in input!");
 	}
 	else
@@ -136,5 +189,10 @@
 	npy_intp*  dims=NULL;
 
-	/*retrive dimensions: */
+	/*intermediary:*/
+	long* lvector=NULL;
+	bool* bvector=NULL;
+	int i;
+
+	/*retrieve dimensions: */
 	ndim=PyArray_NDIM((const PyArrayObject*)py_vector);
 	if(ndim!=1)_error_("expecting an Mx1 vector in input!");
@@ -143,10 +201,33 @@
 
 	if (M) {
-		/*retrieve internal value: */
-		dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector);
-
-		/*copy vector: */
-		vector=xNew<double>(M);
-		memcpy(vector,dvector,(M)*sizeof(double));
+		if      (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_DOUBLE) {
+			/*retrieve internal value: */
+			dvector=(double*)PyArray_DATA((PyArrayObject*)py_vector);
+
+			/*copy vector: */
+			vector=xNew<double>(M);
+			memcpy(vector,dvector,(M)*sizeof(double));
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_INT64) {
+			/*retrieve internal value: */
+			lvector=(long*)PyArray_DATA((PyArrayObject*)py_vector);
+
+			/*transform into double vector: */
+			vector=xNew<double>(M);
+			for(i=0;i<M;i++)vector[i]=(double)lvector[i];
+		}
+
+		else if (PyArray_TYPE((PyArrayObject*)py_vector) == NPY_BOOL) {
+			/*retrieve internal value: */
+			bvector=(bool*)PyArray_DATA((PyArrayObject*)py_vector);
+
+			/*transform into double vector: */
+			vector=xNew<double>(M);
+			for(i=0;i<M;i++)vector[i]=(double)bvector[i];
+		}
+
+		else
+			_error_("unrecognized vector type in input!");
 	}
 	else
Index: /issm/trunk/src/wrappers/python/io/WritePythonData.cpp
===================================================================
--- /issm/trunk/src/wrappers/python/io/WritePythonData.cpp	(revision 14066)
+++ /issm/trunk/src/wrappers/python/io/WritePythonData.cpp	(revision 14067)
@@ -40,4 +40,36 @@
 	dims[1]=(npy_intp)N;
 	array=PyArray_SimpleNewFromData(2,dims,NPY_DOUBLE,matrix);
+
+	PyTuple_SetItem(tuple, index, array);
+}/*}}}*/
+/*FUNCTION WriteData(PyObject* py_tuple,int index, int* matrix, int M, int N){{{*/
+void WriteData(PyObject* tuple, int index, int* matrix, int M,int N){
+
+	npy_intp dims[2]={0,0};
+	PyObject* array=NULL;
+
+	/*intermediary:*/
+	long* lmatrix=NULL;
+	int i;
+
+	/*transform into long matrix: */
+	lmatrix=xNew<long>(M*N);
+	for(i=0;i<M*N;i++)lmatrix[i]=(long)matrix[i];
+
+	dims[0]=(npy_intp)M;
+	dims[1]=(npy_intp)N;
+	array=PyArray_SimpleNewFromData(2,dims,NPY_INT64,lmatrix);
+
+	PyTuple_SetItem(tuple, index, array);
+}/*}}}*/
+/*FUNCTION WriteData(PyObject* py_tuple,int index, bool* matrix, int M, int N){{{*/
+void WriteData(PyObject* tuple, int index, bool* matrix, int M,int N){
+
+	npy_intp dims[2]={0,0};
+	PyObject* array=NULL;
+
+	dims[0]=(npy_intp)M;
+	dims[1]=(npy_intp)N;
+	array=PyArray_SimpleNewFromData(2,dims,NPY_BOOL,matrix);
 
 	PyTuple_SetItem(tuple, index, array);
Index: /issm/trunk/src/wrappers/python/io/pythonio.h
===================================================================
--- /issm/trunk/src/wrappers/python/io/pythonio.h	(revision 14066)
+++ /issm/trunk/src/wrappers/python/io/pythonio.h	(revision 14067)
@@ -18,4 +18,6 @@
 
 void WriteData(PyObject* py_tuple,int index, double* matrix, int M,int N);
+void WriteData(PyObject* py_tuple,int index, int* matrix, int M,int N);
+void WriteData(PyObject* py_tuple,int index, bool* matrix, int M,int N);
 void WriteData(PyObject* py_tuple,int index, int integer);
 void WriteData(PyObject* py_tuple,int index, char* string);
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test1501.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test1501.py	(revision 14066)
+++ 	(revision )
@@ -1,267 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1501.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-printingflag = false
-
-
-md=triangle(model(),'../Exp/Square.exp',350000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelf.py')
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isthermal=0
-
-
-md.timestepping.time_step=1
-md.settings.output_frequency=1
-md.timestepping.final_time=2000
-
-
-# Solve for thinning rate -> -1 * surface mass balance
-
-smb= 2*ones(md.mesh.numberofvertices,1)   
-md.surfaceforcings.mass_balance= smb
-md.basalforcings.melting_rate= smb
-
-
-md=solve(md,PrognosticSolutionEnum())
-
-
-for i=1:10
-	 md=solve(md,PrognosticSolutionEnum())
-	 md.surfaceforcings.mass_balance= md.surfaceforcings.mass_balance - (md.results['PrognosticSolution'][1]['Thickness']-md.geometry.thickness)
-end
-
-
-# Set up transient
-
-smb = md.surfaceforcings.mass_balance
-
-
-tooth= [ [ones(400,1)*(smb') - 10]' [ones(400,1)*(smb')]' ]
-smb=[ [ones(399,1)*(smb')]' smb  tooth tooth]
-
-
-md.surfaceforcings.mass_balance= smb
-md.surfaceforcings.mass_balance(end+1,:)=[1:2000]
-
-
-md=solve(md,TransientSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names=['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SurfaceforcingsMassBalance1', \
-	'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SurfaceforcingsMassBalance2', \
-	'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SurfaceforcingsMassBalance3', \
-	'Vx4','Vy4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SurfaceforcingsMassBalance4', \
-	'Vx5','Vy5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SurfaceforcingsMassBalance5']
-field_tolerances=[1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
-	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
-	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
-	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
-	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10]
-field_values=[\
-	md.results['TransientSolution'][400]['Vx'],\
-	md.results['TransientSolution'][400]['Vy'],\
-	md.results['TransientSolution'][400]['Vel'],\
-	md.results['TransientSolution'][400]['Pressure'],\
-	md.results['TransientSolution'][400]['Bed'],\
-	md.results['TransientSolution'][400]['Surface'],\
-	md.results['TransientSolution'][400]['Thickness'],\
-	md.results['TransientSolution'][400]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][800]['Vx'],\
-	md.results['TransientSolution'][800]['Vy'],\
-	md.results['TransientSolution'][800]['Vel'],\
-	md.results['TransientSolution'][800]['Pressure'],\
-	md.results['TransientSolution'][800]['Bed'],\
-	md.results['TransientSolution'][800]['Surface'],\
-	md.results['TransientSolution'][800]['Thickness'],\
-	md.results['TransientSolution'][800]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][1200]['Vx'],\
-	md.results['TransientSolution'][1200]['Vy'],\
-	md.results['TransientSolution'][1200]['Vel'],\
-	md.results['TransientSolution'][1200]['Pressure'],\
-	md.results['TransientSolution'][1200]['Bed'],\
-	md.results['TransientSolution'][1200]['Surface'],\
-	md.results['TransientSolution'][1200]['Thickness'],\
-	md.results['TransientSolution'][1200]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][1600]['Vx'],\
-	md.results['TransientSolution'][1600]['Vy'],\
-	md.results['TransientSolution'][1600]['Vel'],\
-	md.results['TransientSolution'][1600]['Pressure'],\
-	md.results['TransientSolution'][1600]['Bed'],\
-	md.results['TransientSolution'][1600]['Surface'],\
-	md.results['TransientSolution'][1600]['Thickness'],\
-	md.results['TransientSolution'][1600]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][2000]['Vx'],\
-	md.results['TransientSolution'][2000]['Vy'],\
-	md.results['TransientSolution'][2000]['Vel'],\
-	md.results['TransientSolution'][2000]['Pressure'],\
-	md.results['TransientSolution'][2000]['Bed'],\
-	md.results['TransientSolution'][2000]['Surface'],\
-	md.results['TransientSolution'][2000]['Thickness'],\
-	md.results['TransientSolution'][2000]['SurfaceforcingsMassBalance'],\
-	]
-
-
-if printingflag,
-
-
-	starttime = 360
-	endtime = 2000
-	res = 40
-	ts = [starttime:res:endtime]
-
-
-	index = md.mesh.elements
-	x1=md.mesh.x(index(:,1)) x2=md.mesh.x(index(:,2)) x3=md.mesh.x(index(:,3))
-	y1=md.mesh.y(index(:,1)) y2=md.mesh.y(index(:,2)) y3=md.mesh.y(index(:,3))
-	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)))
-
-
-	thickness = []
-	volume = []
-	massbal = []
-	velocity = []
-	for t=starttime:endtime
-		thickness = [thickness md.results['TransientSolution'][t]['Thickness']]
-		volume = [volume meanmd.results['TransientSolution'][t]['Thicknessvalue,2'].*areas]
-		massbal = [massbal md.results['TransientSolution'][t]['SurfaceforcingsMassBalance']]
-		velocity = [velocity md.results['TransientSolution'][t]['Vel']]
-	end
-
-
-	figure('Position', [0 0 860 932])
-
-
-	options = plotoptions('data','transient_movie','unit','km')
-	options = options.list[1]
-	options = checkplotoptions(md,options)
-
-
-# 	loop over the time steps
-
-	results=md.results.TransientSolution
-	count = 1
-	for i=ts
-
-
-		subplot(5,9,[28:31 37:40])
-		set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
-		field = 'Thickness'
-
-
-# 		process data
-
-		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options)
-		[data datatype]=processdata(md,results(i).(field),options)
-
-
-		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year']
-		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-		options=changefieldvalue(options,'title',titlestring)
-		options=addfielddefault(options,'colorbar',1)
-		options=changefieldvalue(options,'caxis',[0 max(max(thickness))])
-		applyoptions(md,[],options)
-
-
-		subplot(5,9,[33:36 42:45])
-		set(gca,'pos',get(gca,'pos')+[-0.00 -0.08 0.07 0.08])
-		field = 'Vel'
-
-
-# 		process data
-
-		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options)
-		[data datatype]=processdata(md,results(i).(field),options)
-
-
-		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year']
-		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-		options=changefieldvalue(options,'title',titlestring)
-		options=addfielddefault(options,'colorbar',1)
-		options=changefieldvalue(options,'caxis',[0 max(max(velocity))])
-		applyoptions(md,[],options)
-
-
-		subplot(5,4,1:4)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
-		plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Surface Mass Balance','FontSize',14)
-		ylabel('m/year','FontSize',14)
-
-
-		subplot(5,4,5:8)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
-		plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Ice Volume','FontSize',14)
-		ylabel('km^3','FontSize',14)
-
-
-		subplot(5,4,9:12)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
-		plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Mean Velocity','FontSize', 14)
-		ylabel('km/year','FontSize', 14)
-		xlabel('year','FontSize', 14)
-
-
-		set(gcf,'Renderer','zbuffer','color','white') %fixes a bug on Mac OS X (not needed in future Matlab version)
-		if i==starttime,
-# 			initialize images and frame
-
-			frame=getframe(gcf)
-			[images,map]=rgb2ind(frame.cdata,256,'nodither')
-			images(1,1,1,length(ts))=0
-		else
-			frame=getframe(gcf)
-			images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither')
-		end
-
-
-		count = count+1
-
-
-	end
-
-
-	filename='transawtooth2d.gif'
-	imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
-
-
-end %printingflag
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test1502.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test1502.py	(revision 14066)
+++ 	(revision )
@@ -1,273 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1502.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-printingflag = false
-
-
-md=triangle(model(),'../Exp/Square.exp',450000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelf.py')
-md=setflowequation(md,'macayeal','all')
-md.extrude(3,1.)
-md.cluster=generic('name',oshostname(),'np',2)
-md.transient.isthermal=0
-
-
-md.timestepping.time_step=1
-md.settings.output_frequency=1
-md.timestepping.final_time=2000
-
-
-# Solve for thinning rate -> -1 * surface mass balance
-
-smb= 2*ones(md.mesh.numberofvertices,1)   
-md.surfaceforcings.mass_balance= smb
-md.basalforcings.melting_rate= smb
-
-
-md=solve(md,PrognosticSolutionEnum())
-
-
-for i=1:10
-	 md=solve(md,PrognosticSolutionEnum())
-	 md.surfaceforcings.mass_balance= md.surfaceforcings.mass_balance - (md.results['PrognosticSolution'][1]['Thickness']-md.geometry.thickness)
-end
-
-
-# Set up transient
-
-smb = md.surfaceforcings.mass_balance
-
-
-tooth= [ [ones(400,1)*(smb') - 10]' [ones(400,1)*(smb')]' ]
-smb=[ [ones(399,1)*(smb')]' smb  tooth tooth]
-
-
-md.surfaceforcings.mass_balance= smb
-md.surfaceforcings.mass_balance(end+1,:)=[1:2000]
-
-
-md=solve(md,TransientSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names=['Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SurfaceforcingsMassBalance1', \
-	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SurfaceforcingsMassBalance2', \
-	'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SurfaceforcingsMassBalance3', \
-	'Vx4','Vy4','Vz4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SurfaceforcingsMassBalance4', \
-	'Vx5','Vy5','Vz5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SurfaceforcingsMassBalance5']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
-field_values=[\
-	md.results['TransientSolution'][400]['Vx'],\
-	md.results['TransientSolution'][400]['Vy'],\
-	md.results['TransientSolution'][400]['Vz'],\
-	md.results['TransientSolution'][400]['Vel'],\
-	md.results['TransientSolution'][400]['Pressure'],\
-	md.results['TransientSolution'][400]['Bed'],\
-	md.results['TransientSolution'][400]['Surface'],\
-	md.results['TransientSolution'][400]['Thickness'],\
-	md.results['TransientSolution'][400]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][800]['Vx'],\
-	md.results['TransientSolution'][800]['Vy'],\
-	md.results['TransientSolution'][800]['Vz'],\
-	md.results['TransientSolution'][800]['Vel'],\
-	md.results['TransientSolution'][800]['Pressure'],\
-	md.results['TransientSolution'][800]['Bed'],\
-	md.results['TransientSolution'][800]['Surface'],\
-	md.results['TransientSolution'][800]['Thickness'],\
-	md.results['TransientSolution'][800]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][1200]['Vx'],\
-	md.results['TransientSolution'][1200]['Vy'],\
-	md.results['TransientSolution'][1200]['Vz'],\
-	md.results['TransientSolution'][1200]['Vel'],\
-	md.results['TransientSolution'][1200]['Pressure'],\
-	md.results['TransientSolution'][1200]['Bed'],\
-	md.results['TransientSolution'][1200]['Surface'],\
-	md.results['TransientSolution'][1200]['Thickness'],\
-	md.results['TransientSolution'][1200]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][1600]['Vx'],\
-	md.results['TransientSolution'][1600]['Vy'],\
-	md.results['TransientSolution'][1600]['Vz'],\
-	md.results['TransientSolution'][1600]['Vel'],\
-	md.results['TransientSolution'][1600]['Pressure'],\
-	md.results['TransientSolution'][1600]['Bed'],\
-	md.results['TransientSolution'][1600]['Surface'],\
-	md.results['TransientSolution'][1600]['Thickness'],\
-	md.results['TransientSolution'][1600]['SurfaceforcingsMassBalance'],\
-	md.results['TransientSolution'][2000]['Vx'],\
-	md.results['TransientSolution'][2000]['Vy'],\
-	md.results['TransientSolution'][2000]['Vz'],\
-	md.results['TransientSolution'][2000]['Vel'],\
-	md.results['TransientSolution'][2000]['Pressure'],\
-	md.results['TransientSolution'][2000]['Bed'],\
-	md.results['TransientSolution'][2000]['Surface'],\
-	md.results['TransientSolution'][2000]['Thickness'],\
-	md.results['TransientSolution'][2000]['SurfaceforcingsMassBalance'],\
-	]
-
-
-if printingflag,
-
-
-	starttime = 360
-	endtime = 2000
-	res = 40
-	ts = [starttime:res:endtime]
-
-
-	index = md.mesh.elements
-	x1=md.mesh.x(index(:,1)) x2=md.mesh.x(index(:,2)) x3=md.mesh.x(index(:,3))
-	y1=md.mesh.y(index(:,1)) y2=md.mesh.y(index(:,2)) y3=md.mesh.y(index(:,3))
-	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)))
-
-
-	thickness = []
-	volume = []
-	massbal = []
-	velocity = []
-	for t=starttime:endtime
-		thickness = [thickness md.results['TransientSolution'][t]['Thickness']]
-		volume = [volume meanmd.results['TransientSolution'][t]['Thicknessvalue,2'].*areas]
-		massbal = [massbal md.results['TransientSolution'][t]['SurfaceforcingsMassBalance']]
-		velocity = [velocity md.results['TransientSolution'][t]['Vel']]
-	end
-
-
-	figure('Position', [0 0 1060 1060])
-
-
-	options = plotoptions('data','transient_movie','unit','km')
-	options = options.list[1]
-	options = checkplotoptions(md,options)
-
-
-# 	loop over the time steps
-
-	results=md.results.TransientSolution
-	count = 1
-	for i=ts
-
-
-		subplot(5,9,[28:31 37:40])
-		set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
-		field = 'Thickness'
-
-
-# 		process data
-
-		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options)
-		[data datatype]=processdata(md,results(i).(field),options)
-
-
-		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year']
-		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-		options=changefieldvalue(options,'title',titlestring)
-		options=addfielddefault(options,'colorbar',1)
-		options=changefieldvalue(options,'caxis',[0 max(max(thickness))])
-		applyoptions(md,[],options)
-
-
-		subplot(5,9,[33:36 42:45])
-		set(gca,'pos',get(gca,'pos')+[-0.01 -0.08 0.07 0.08])
-		field = 'Vel'
-
-
-# 		process data
-
-		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options)
-		[data datatype]=processdata(md,results(i).(field),options)
-
-
-		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year']
-		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
-		options=changefieldvalue(options,'title',titlestring)
-		options=addfielddefault(options,'colorbar',1)
-		options=changefieldvalue(options,'caxis',[0 max(max(velocity))])
-		applyoptions(md,[],options)
-
-
-		subplot(5,4,1:4)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
-		plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Surface Mass Balance','FontSize',14)
-		ylabel('m/year','FontSize',14)
-
-
-		subplot(5,4,5:8)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
-		plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Ice Volume','FontSize',14)
-		ylabel('km^3','FontSize',14)
-
-
-		subplot(5,4,9:12)
-		cla
-		set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
-		plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
-		hold on
-		ya = ylim
-		plot([i i], ya, 'r', 'LineWidth',6)
-		ylim(ya) xlim([starttime endtime])
-		title('Mean Velocity','FontSize', 14)
-		ylabel('km/year','FontSize', 14)
-		xlabel('year','FontSize', 14)
-
-
-		set(gcf,'Renderer','zbuffer','color','white') %fixes a bug on Mac OS X (not needed in future Matlab version)
-		if i==starttime,
-# 			initialize images and frame
-
-			frame=getframe(gcf)
-			[images,map]=rgb2ind(frame.cdata,256,'nodither')
-			images(1,1,1,length(ts))=0
-		else
-			frame=getframe(gcf)
-			images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither')
-		end
-
-
-		count = count+1
-
-
-	end
-
-
-	filename='transawtooth3d.gif'
-	imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
-
-
-end %printingflag
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test1601.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test1601.py	(revision 14066)
+++ 	(revision )
@@ -1,76 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1601.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelf.py')
-md=setflowequation(md,'macayeal','all')
-md.diagnostic.spcvx(numpy.nonzero(md.mesh.y>0))=NaN
-md.initialization.vx(:)=0
-md.initialization.vy(:)=0
-md.initialization.vel(:)=0
-
-
-md.cluster=generic('name',oshostname,'np',2)
-md=solve(md,DiagnosticSolutionEnum())
-vel0=md.results.DiagnosticSolution.Vel
-
-
-theta=30*pi/180
-x=md.mesh.x
-y=md.mesh.y
-md.mesh.x=cos(theta)*x-sin(theta)*y
-md.mesh.y=sin(theta)*x+cos(theta)*y
-
-
-md.diagnostic.referential(:,1:3)=repmat([cos(theta),sin(theta),0],md.mesh.numberofvertices,1)
-md.diagnostic.referential(:,4:6)=repmat([0,0,1],md.mesh.numberofvertices,1)
-md=solve(md,DiagnosticSolutionEnum())
-vel1=md.results.DiagnosticSolution.Vel
-
-
-plotmodel()(md,'data',vel0,'data',vel1,'data',vel1-vel0,'title','Cartesian CS','title','Rotated CS','title','difference')
-disp(['Error between Cartesian and rotated CS: ' num2str(max(abs(vel0-vel1))/(max(abs(vel0))+eps)) ])
-
-
-# Now, put CS back to normal except on the side where the spc are applied
-
-pos=find(x==0 | x==1000000)
-md.diagnostic.referential(:)=NaN
-md.diagnostic.referential(pos,1:3)=repmat([cos(theta),sin(theta),0],size(pos,1),1)
-md.diagnostic.referential(pos,4:6)=repmat([0,0,1],size(pos,1),1)
-md=solve(md,DiagnosticSolutionEnum())
-vel2=md.results.DiagnosticSolution.Vel
-
-
-plotmodel()(md,'data',vel0,'data',vel2,'data',vel2-vel0,'title','Cartesian CS','title','Rotated CS','title','difference')
-disp(['Error between Cartesian and rotated CS: ' num2str(max(abs(vel0-vel2))/(max(abs(vel0))+eps)) ])
-
-
-# Fields and tolerances to track changes
-
-field_names     =['vel1','vel2']
-field_tolerances=[1e-11,1e-11]
-field_values=[\
-	vel1, \
-	vel2, \
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test1602.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test1602.py	(revision 14066)
+++ 	(revision )
@@ -1,62 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test1602.m
-Created on 2012-09-27 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-from numpy import *
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=setmask(md,'../Exp/SquareShelf.exp','')
-md=parameterize(md,'../Par/SquareSheetShelf.py')
-md.extrude(5,1.)
-md=setflowequation(md,'pattyn','all')
-md.diagnostic.spcvx(numpy.nonzero(md.mesh.y>0))=NaN
-md.initialization.vx(:)=0
-md.initialization.vy(:)=0
-md.initialization.vel(:)=0
-
-
-md.cluster=generic('name',oshostname,'np',3)
-md=solve(md,DiagnosticSolutionEnum())
-vel0=md.results.DiagnosticSolution.Vel
-
-
-theta=30*pi/180
-x=md.mesh.x
-y=md.mesh.y
-md.mesh.x=cos(theta)*x-sin(theta)*y
-md.mesh.y=sin(theta)*x+cos(theta)*y
-
-
-md.diagnostic.referential(:,1:3)=repmat([cos(theta),sin(theta),0],md.mesh.numberofvertices,1)
-md.diagnostic.referential(:,4:6)=repmat([0,0,1],md.mesh.numberofvertices,1)
-md=solve(md,DiagnosticSolutionEnum())
-vel1=md.results.DiagnosticSolution.Vel
-
-
-plotmodel()(md,'data',vel0,'data',vel1,'data',vel1-vel0,'title','Cartesian CS','title','Rotated CS','title','difference','view#all',2)
-disp(['Error between Cartesian and rotated CS: ' num2str(max(abs(vel0-vel1))/(max(abs(vel0))+eps)) ])
-
-
-# Fields and tolerances to track changes
-
-field_names     =['vel1']
-field_tolerances=[1e-9]
-field_values=[\
-	vel1, \
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3001.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3001.py	(revision 14066)
+++ 	(revision )
@@ -1,47 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3001.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',50000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.diagnostic.requested_outputs=StressTensorEnum()
-md.autodiff.isautodiff=true
-md=solve(md,DiagnosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Vx','Vy','Vel','Pressure',\
-	'StressTensorxx','StressTensoryy','StressTensorxy']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,\
-	1e-13,1e-13,1e-13]
-field_values=[\
-	md.results['DiagnosticSolution'][1]['Vx'],\
-	md.results['DiagnosticSolution'][1]['Vy'],\
-	md.results['DiagnosticSolution'][1]['Vel'],\
-	md.results['DiagnosticSolution'][1]['Pressure'],\
-	md.results['DiagnosticSolution'][1]['StressTensorxx'],\
-	md.results['DiagnosticSolution'][1]['StressTensoryy'],\
-	md.results['DiagnosticSolution'][1]['StressTensorxy'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3002.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3002.py	(revision 14066)
+++ 	(revision )
@@ -1,43 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3002.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',180000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md.extrude(3,2.)
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.autodiff.isautodiff=true
-md=solve(md,DiagnosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Vx','Vy','Vz','Vel','Pressure']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13]
-field_values=[\
-	md.results['DiagnosticSolution'][1]['Vx'],\
-	md.results['DiagnosticSolution'][1]['Vy'],\
-	md.results['DiagnosticSolution'][1]['Vz'],\
-	md.results['DiagnosticSolution'][1]['Vel'],\
-	md.results['DiagnosticSolution'][1]['Pressure'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3003.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3003.py	(revision 14066)
+++ 	(revision )
@@ -1,52 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3003.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',180000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md.extrude(3,2.)
-md=setflowequation(md,'pattyn','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.diagnostic.requested_outputs=StressTensorEnum()
-md.autodiff.isautodiff=true
-md=solve(md,DiagnosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Vx','Vy','Vz','Vel','Pressure',\
-	'StressTensorxx','StressTensoryy','StressTensorzz','StressTensorxy','StressTensorxz','StressTensoryz']
-field_tolerances=[1e-09,1e-09,1e-09,1e-09,1e-09,\
-	1e-09,1e-09,1e-09,1e-09,1e-09,1e-09]
-field_values=[\
-	md.results['DiagnosticSolution'][1]['Vx'],\
-	md.results['DiagnosticSolution'][1]['Vy'],\
-	md.results['DiagnosticSolution'][1]['Vz'],\
-	md.results['DiagnosticSolution'][1]['Vel'],\
-	md.results['DiagnosticSolution'][1]['Pressure'],\
-	md.results['DiagnosticSolution'][1]['StressTensorxx'],\
-	md.results['DiagnosticSolution'][1]['StressTensoryy'],\
-	md.results['DiagnosticSolution'][1]['StressTensorzz'],\
-	md.results['DiagnosticSolution'][1]['StressTensorxy'],\
-	md.results['DiagnosticSolution'][1]['StressTensorxz'],\
-	md.results['DiagnosticSolution'][1]['StressTensoryz'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3004.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3004.py	(revision 14066)
+++ 	(revision )
@@ -1,43 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3004.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',180000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md.extrude(3,2.)
-md=setflowequation(md,'stokes','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.autodiff.isautodiff=true
-md=solve(md,DiagnosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Vx','Vy','Vz','Vel','Pressure']
-field_tolerances=[1e-08,1e-08,1e-07,1e-08,1e-08]
-field_values=[\
-	md.results['DiagnosticSolution'][1]['Vx'],\
-	md.results['DiagnosticSolution'][1]['Vy'],\
-	md.results['DiagnosticSolution'][1]['Vz'],\
-	md.results['DiagnosticSolution'][1]['Vel'],\
-	md.results['DiagnosticSolution'][1]['Pressure'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3005.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3005.py	(revision 14066)
+++ 	(revision )
@@ -1,38 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3005.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.autodiff.isautodiff=true
-md=solve(md,PrognosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Thickness']
-field_tolerances=[1e-13]
-field_values=[\
-	md.results['PrognosticSolution'][1]['Thickness'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3006.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3006.py	(revision 14066)
+++ 	(revision )
@@ -1,42 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3006.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from meshconvert import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=meshconvert(md)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.prognostic.stabilization=3
-md.prognostic.spcthickness=md.geometry.thickness
-md.autodiff.isautodiff=true
-md=solve(md,PrognosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Thickness']
-field_tolerances=[1e-13]
-field_values=[\
-	md.results['PrognosticSolution'][1]['Thickness'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3007.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3007.py	(revision 14066)
+++ 	(revision )
@@ -1,39 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3007.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md=setflowequation(md,'macayeal','all')
-md.extrude(5,3.)
-md.cluster=generic('name',oshostname(),'np',3)
-md.autodiff.isautodiff=true
-md=solve(md,PrognosticSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Thickness']
-field_tolerances=[1e-13]
-field_values=[\
-	md.results['PrognosticSolution'][1]['Thickness'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3008.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3008.py	(revision 14066)
+++ 	(revision )
@@ -1,41 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3008.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',180000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md.extrude(3,1.)
-md=setflowequation(md,'macayeal','all')
-md.timestepping.time_step=0
-md.cluster=generic('name',oshostname(),'np',3)
-md.autodiff.isautodiff=true
-md=solve(md,ThermalSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Temperature','BasalforcingsMeltingRate']
-field_tolerances=[1e-13,1e-13]
-field_values=[\
-	md.results['ThermalSolution'][1]['Temperature'],\
-	md.results['ThermalSolution'][1]['BasalforcingsMeltingRate'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3009.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3009.py	(revision 14066)
+++ 	(revision )
@@ -1,44 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3009.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',180000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md.extrude(3,1.)
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
-md.autodiff.isautodiff=true
-md=solve(md,TransientSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Temperature','BasalforcingsMeltingRate']
-field_tolerances=[1e-13,1e-13]
-field_values=[\
-	md.results['TransientSolution'][1]['Temperature'],\
-	md.results['TransientSolution'][1]['BasalforcingsMeltingRate'],\
-	]
Index: sm/trunk/test/NightlyRun/InNeedOfDebugging/test3010.py
===================================================================
--- /issm/trunk/test/NightlyRun/InNeedOfDebugging/test3010.py	(revision 14066)
+++ 	(revision )
@@ -1,66 +1,0 @@
-"""
-== == == == == == == == == == == == == == == == == == ==
-Auto generated python script for ISSM:   test3010.m
-Created on 2012-09-25 via translateToPy.py Ver 1.0 by mikep
-== == == == == == == == == == == == == == == == == == ==
-
-Matlab script conversion into python
-translateToPy.py Author: Michael Pellegrin
-translateToPy.py Date: 09/24/12
-== == == == == == == == == == == == == == == == == == ==
-"""
-
-from MatlabFuncs import *
-from model import *
-from EnumDefinitions import *
-import numpy
-from triangle import *
-from setmask import *
-from parameterize import *
-from setflowequation import *
-from solve import *
-
-md=triangle(model(),'../Exp/Square.exp',150000)
-md=setmask(md,'all','')
-md=parameterize(md,'../Par/SquareShelfConstrained.py')
-md=setflowequation(md,'macayeal','all')
-md.cluster=generic('name',oshostname(),'np',3)
-md.transient.requested_outputs=IceVolumeEnum()
-
-
-md.autodiff.isautodiff=true
-md=solve(md,TransientSolutionEnum())
-
-
-# Fields and tolerances to track changes
-
-field_names     =['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','Volume1','Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','Volume2','Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','Volume3']
-field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-						1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
-						1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
-field_values=[\
-	md.results['TransientSolution'][1]['Vx'],\
-	md.results['TransientSolution'][1]['Vy'],\
-	md.results['TransientSolution'][1]['Vel'],\
-	md.results['TransientSolution'][1]['Pressure'],\
-	md.results['TransientSolution'][1]['Bed'],\
-	md.results['TransientSolution'][1]['Surface'],\
-	md.results['TransientSolution'][1]['Thickness'],\
-	md.results['TransientSolution'][1]['IceVolume'],\
-	md.results['TransientSolution'][2]['Vx'],\
-	md.results['TransientSolution'][2]['Vy'],\
-	md.results['TransientSolution'][2]['Vel'],\
-	md.results['TransientSolution'][2]['Pressure'],\
-	md.results['TransientSolution'][2]['Bed'],\
-	md.results['TransientSolution'][2]['Surface'],\
-	md.results['TransientSolution'][2]['Thickness'],\
-	md.results['TransientSolution'][2]['IceVolume'],\
-	md.results['TransientSolution'][3]['Vx'],\
-	md.results['TransientSolution'][3]['Vy'],\
-	md.results['TransientSolution'][3]['Vel'],\
-	md.results['TransientSolution'][3]['Pressure'],\
-	md.results['TransientSolution'][3]['Bed'],\
-	md.results['TransientSolution'][3]['Surface'],\
-	md.results['TransientSolution'][3]['Thickness'],\
-	md.results['TransientSolution'][3]['IceVolume'],\
-	]
Index: /issm/trunk/test/NightlyRun/runme.m
===================================================================
--- /issm/trunk/test/NightlyRun/runme.m	(revision 14066)
+++ /issm/trunk/test/NightlyRun/runme.m	(revision 14067)
@@ -97,4 +97,6 @@
 if strcmpi(benchmark,'nightly'),
 	test_ids=intersect(test_ids,[1:999]);
+elseif strcmpi(benchmark,'validation'),
+	test_ids=intersect(test_ids,[1001:1999]);
 elseif strcmpi(benchmark,'ismip'),
 	test_ids=intersect(test_ids,[1101:1199]);
@@ -105,10 +107,10 @@
 elseif strcmpi(benchmark,'mesh'),
 	test_ids=intersect(test_ids,[1401:1499]);
+elseif strcmpi(benchmark,'tranforcing'),
+	test_ids=intersect(test_ids,[1501:1502]);
+elseif strcmpi(benchmark,'referential'),
+	test_ids=intersect(test_ids,[1601:1602]);
 elseif strcmpi(benchmark,'adolc'),
 	test_ids=intersect(test_ids,[3001:3020]);
-elseif strcmpi(benchmark,'validation'),
-	test_ids=intersect(test_ids,[1001:1999]);
-elseif strcmpi(benchmark,'tranforcing'),
-	test_ids=intersect(test_ids,[1501:1502]);
 end
 % }}}
Index: /issm/trunk/test/NightlyRun/runme.py
===================================================================
--- /issm/trunk/test/NightlyRun/runme.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/runme.py	(revision 14067)
@@ -99,4 +99,6 @@
 	if   strcmpi(benchmark,'nightly'):
 		test_ids=test_ids.intersection(set(range(1,1000)))
+	elif strcmpi(benchmark,'validation'):
+		test_ids=test_ids.intersection(set(range(1001,2000)))
 	elif strcmpi(benchmark,'ismip'):
 		test_ids=test_ids.intersection(set(range(1101,1200)))
@@ -107,10 +109,10 @@
 	elif strcmpi(benchmark,'mesh'):
 		test_ids=test_ids.intersection(set(range(1401,1500)))
+	elif strcmpi(benchmark,'tranforcing'):
+		test_ids=test_ids.intersection(set(range(1501,1503)))
+	elif strcmpi(benchmark,'referential'):
+		test_ids=test_ids.intersection(set(range(1601,1603)))
 	elif strcmpi(benchmark,'adolc'):
 		test_ids=test_ids.intersection(set(range(3001,3020)))
-	elif strcmpi(benchmark,'validation'):
-		test_ids=test_ids.intersection(set(range(1001,2000)))
-	elif strcmpi(benchmark,'tranforcing'):
-		test_ids=test_ids.intersection(set(range(1501,1503)))
 	#print 'test_ids after benchmark =',test_ids
 	test_ids=list(test_ids)
Index: /issm/trunk/test/NightlyRun/test109.py
===================================================================
--- /issm/trunk/test/NightlyRun/test109.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test109.py	(revision 14067)
@@ -14,8 +14,8 @@
 md=setflowequation(md,'macayeal','all')
 md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())
 
Index: /issm/trunk/test/NightlyRun/test121.py
===================================================================
--- /issm/trunk/test/NightlyRun/test121.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test121.py	(revision 14067)
@@ -15,8 +15,8 @@
 md.cluster=generic('name',oshostname(),'np',3);
 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md.thermal.isenthalpy=1
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test1501.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1501.m	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test1501.m	(revision 14067)
@@ -13,5 +13,5 @@
 
 %Solve for thinning rate -> -1 * surface mass balance
-smb= 2.*ones(md.mesh.numberofvertices,1);   
+smb= 2.*ones(md.mesh.numberofvertices,1);
 md.surfaceforcings.mass_balance= smb;
 md.basalforcings.melting_rate= smb;
Index: /issm/trunk/test/NightlyRun/test1501.py
===================================================================
--- /issm/trunk/test/NightlyRun/test1501.py	(revision 14067)
+++ /issm/trunk/test/NightlyRun/test1501.py	(revision 14067)
@@ -0,0 +1,222 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+printingflag = False
+
+md=triangle(model(),'../Exp/Square.exp',350000.)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelf.py')
+md=setflowequation(md,'macayeal','all')
+md.cluster=generic('name',oshostname(),'np',3)
+md.transient.isthermal=False
+
+md.timestepping.time_step=1.
+md.settings.output_frequency=1
+md.timestepping.final_time=2000.
+
+#Solve for thinning rate -> -1 * surface mass balance
+smb= 2.*numpy.ones((md.mesh.numberofvertices,1))
+md.surfaceforcings.mass_balance= smb
+md.basalforcings.melting_rate= smb
+
+md=solve(md,PrognosticSolutionEnum())
+
+for i in xrange(1,11):
+	 md=solve(md,PrognosticSolutionEnum())
+	 md.surfaceforcings.mass_balance= md.surfaceforcings.mass_balance - ((md.results.PrognosticSolution.Thickness)-md.geometry.thickness)
+
+#Set up transient
+smb = md.surfaceforcings.mass_balance
+
+#tooth= [ [ones(400,1)*(smb') - 10.]' [ones(400,1)*(smb')]' ];
+tooth=numpy.hstack((numpy.tile(smb-10.,(1,400)),numpy.tile(smb,(1,400))))
+#smb=[ [ones(399,1)*(smb')]' smb  tooth tooth];
+smb=numpy.hstack((numpy.tile(smb,(1,399)),smb,tooth,tooth))
+
+#md.surfaceforcings.mass_balance= smb;
+#md.surfaceforcings.mass_balance(end+1,:)=[1.:2000.];
+md.surfaceforcings.mass_balance=numpy.vstack((smb,numpy.arange(1,2001)))
+
+md=solve(md,TransientSolutionEnum())
+
+#Fields and tolerances to track changes
+field_names=['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SurfaceforcingsMassBalance1', \
+	'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SurfaceforcingsMassBalance2', \
+	'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SurfaceforcingsMassBalance3', \
+	'Vx4','Vy4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SurfaceforcingsMassBalance4', \
+	'Vx5','Vy5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SurfaceforcingsMassBalance5']
+field_tolerances=[1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
+	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
+	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
+	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
+	1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10]
+field_values=[\
+	md.results.TransientSolution[400-1].Vx,\
+	md.results.TransientSolution[400-1].Vy,\
+	md.results.TransientSolution[400-1].Vel,\
+	md.results.TransientSolution[400-1].Pressure,\
+	md.results.TransientSolution[400-1].Bed,\
+	md.results.TransientSolution[400-1].Surface,\
+	md.results.TransientSolution[400-1].Thickness,\
+	md.results.TransientSolution[400-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[800-1].Vx,\
+	md.results.TransientSolution[800-1].Vy,\
+	md.results.TransientSolution[800-1].Vel,\
+	md.results.TransientSolution[800-1].Pressure,\
+	md.results.TransientSolution[800-1].Bed,\
+	md.results.TransientSolution[800-1].Surface,\
+	md.results.TransientSolution[800-1].Thickness,\
+	md.results.TransientSolution[800-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[1200-1].Vx,\
+	md.results.TransientSolution[1200-1].Vy,\
+	md.results.TransientSolution[1200-1].Vel,\
+	md.results.TransientSolution[1200-1].Pressure,\
+	md.results.TransientSolution[1200-1].Bed,\
+	md.results.TransientSolution[1200-1].Surface,\
+	md.results.TransientSolution[1200-1].Thickness,\
+	md.results.TransientSolution[1200-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[1600-1].Vx,\
+	md.results.TransientSolution[1600-1].Vy,\
+	md.results.TransientSolution[1600-1].Vel,\
+	md.results.TransientSolution[1600-1].Pressure,\
+	md.results.TransientSolution[1600-1].Bed,\
+	md.results.TransientSolution[1600-1].Surface,\
+	md.results.TransientSolution[1600-1].Thickness,\
+	md.results.TransientSolution[1600-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[2000-1].Vx,\
+	md.results.TransientSolution[2000-1].Vy,\
+	md.results.TransientSolution[2000-1].Vel,\
+	md.results.TransientSolution[2000-1].Pressure,\
+	md.results.TransientSolution[2000-1].Bed,\
+	md.results.TransientSolution[2000-1].Surface,\
+	md.results.TransientSolution[2000-1].Thickness,\
+	md.results.TransientSolution[2000-1].SurfaceforcingsMassBalance,\
+	]
+
+if printingflag:
+	pass
+
+	"""
+	starttime = 360;
+	endtime = 2000;
+	res = 40;
+	ts = [starttime:res:endtime];
+
+	index = md.mesh.elements;
+	x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3));
+	y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3));
+	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
+
+	thickness = [];
+	volume = [];
+	massbal = [];
+	velocity = [];
+	for t=starttime:endtime
+		thickness = [thickness (md.results.TransientSolution(t).Thickness)];
+		volume = [volume mean(md.results.TransientSolution(t).Thickness.value,2).*areas];
+		massbal = [massbal (md.results.TransientSolution(t).SurfaceforcingsMassBalance)];
+		velocity = [velocity (md.results.TransientSolution(t).Vel)];
+	end
+
+	figure('Position', [0 0 860 932])
+
+	options = plotoptions('data','transient_movie','unit','km');
+	options = options.list{1};
+	options = checkplotoptions(md,options);
+
+	%loop over the time steps
+	results=md.results.TransientSolution;
+	count = 1;
+	for i=ts
+
+		subplot(5,9,[28:31 37:40])
+		set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
+		field = 'Thickness';
+
+		%process data
+		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
+		[data datatype]=processdata(md,results(i).(field),options);
+
+		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
+		options=changefieldvalue(options,'title',titlestring);
+		options=addfielddefault(options,'colorbar',1);
+		options=changefieldvalue(options,'caxis',[0 max(max(thickness))]);
+		applyoptions(md,[],options);
+
+		subplot(5,9,[33:36 42:45])
+		set(gca,'pos',get(gca,'pos')+[-0.00 -0.08 0.07 0.08])
+		field = 'Vel';
+
+		%process data
+		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
+		[data datatype]=processdata(md,results(i).(field),options);
+
+		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
+		options=changefieldvalue(options,'title',titlestring);
+		options=addfielddefault(options,'colorbar',1);
+		options=changefieldvalue(options,'caxis',[0 max(max(velocity))]);
+		applyoptions(md,[],options);
+
+		subplot(5,4,1:4)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
+		plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Surface Mass Balance','FontSize',14)
+		ylabel('m/year','FontSize',14)
+
+		subplot(5,4,5:8)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
+		plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Ice Volume','FontSize',14)
+		ylabel('km^3','FontSize',14)
+
+		subplot(5,4,9:12)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
+		plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Mean Velocity','FontSize', 14)
+		ylabel('km/year','FontSize', 14)
+		xlabel('year','FontSize', 14)
+
+		set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
+		if i==starttime,
+			%initialize images and frame
+			frame=getframe(gcf);
+			[images,map]=rgb2ind(frame.cdata,256,'nodither');
+			images(1,1,1,length(ts))=0;
+		else
+			frame=getframe(gcf);
+			images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither');
+		end
+
+		count = count+1;
+
+	end
+
+	filename='transawtooth2d.gif';
+	imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
+	"""
+
Index: /issm/trunk/test/NightlyRun/test1502.m
===================================================================
--- /issm/trunk/test/NightlyRun/test1502.m	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test1502.m	(revision 14067)
@@ -14,5 +14,5 @@
 
 %Solve for thinning rate -> -1 * surface mass balance
-smb= 2.*ones(md.mesh.numberofvertices,1);   
+smb= 2.*ones(md.mesh.numberofvertices,1);
 md.surfaceforcings.mass_balance= smb;
 md.basalforcings.melting_rate= smb;
Index: /issm/trunk/test/NightlyRun/test1502.py
===================================================================
--- /issm/trunk/test/NightlyRun/test1502.py	(revision 14067)
+++ /issm/trunk/test/NightlyRun/test1502.py	(revision 14067)
@@ -0,0 +1,228 @@
+import numpy
+import sys
+from model import *
+from triangle import *
+from setmask import *
+from parameterize import *
+from setflowequation import *
+from EnumDefinitions import *
+from solve import *
+from MatlabFuncs import *
+
+printingflag = False
+
+md=triangle(model(),'../Exp/Square.exp',450000.)
+md=setmask(md,'all','')
+md=parameterize(md,'../Par/SquareShelf.py')
+md=setflowequation(md,'macayeal','all')
+md.extrude(3,1.)
+md.cluster=generic('name',oshostname(),'np',2)
+md.transient.isthermal=False
+
+md.timestepping.time_step=1.
+md.settings.output_frequency=1
+md.timestepping.final_time=2000.
+
+#Solve for thinning rate -> -1 * surface mass balance
+smb= 2.*numpy.ones((md.mesh.numberofvertices,1))
+md.surfaceforcings.mass_balance= smb
+md.basalforcings.melting_rate= smb
+
+md=solve(md,PrognosticSolutionEnum())
+
+for i in xrange(1,11):
+	 md=solve(md,PrognosticSolutionEnum())
+	 md.surfaceforcings.mass_balance= md.surfaceforcings.mass_balance - ((md.results.PrognosticSolution.Thickness)-md.geometry.thickness)
+
+#Set up transient
+smb = md.surfaceforcings.mass_balance
+
+#tooth= [ [ones(400,1)*(smb') - 10.]' [ones(400,1)*(smb')]' ];
+tooth=numpy.hstack((numpy.tile(smb-10.,(1,400)),numpy.tile(smb,(1,400))))
+#smb=[ [ones(399,1)*(smb')]' smb  tooth tooth];
+smb=numpy.hstack((numpy.tile(smb,(1,399)),smb,tooth,tooth))
+
+#md.surfaceforcings.mass_balance= smb;
+#md.surfaceforcings.mass_balance(end+1,:)=[1.:2000.];
+md.surfaceforcings.mass_balance=numpy.vstack((smb,numpy.arange(1,2001)))
+
+md=solve(md,TransientSolutionEnum())
+
+#Fields and tolerances to track changes
+field_names=['Vx1','Vy1','Vz1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SurfaceforcingsMassBalance1', \
+	'Vx2','Vy2','Vz2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SurfaceforcingsMassBalance2', \
+	'Vx3','Vy3','Vz3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SurfaceforcingsMassBalance3', \
+	'Vx4','Vy4','Vz4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SurfaceforcingsMassBalance4', \
+	'Vx5','Vy5','Vz5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SurfaceforcingsMassBalance5']
+field_tolerances=[1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
+	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
+	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
+	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,\
+	1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13,1e-13]
+field_values=[\
+	md.results.TransientSolution[400-1].Vx,\
+	md.results.TransientSolution[400-1].Vy,\
+	md.results.TransientSolution[400-1].Vz,\
+	md.results.TransientSolution[400-1].Vel,\
+	md.results.TransientSolution[400-1].Pressure,\
+	md.results.TransientSolution[400-1].Bed,\
+	md.results.TransientSolution[400-1].Surface,\
+	md.results.TransientSolution[400-1].Thickness,\
+	md.results.TransientSolution[400-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[800-1].Vx,\
+	md.results.TransientSolution[800-1].Vy,\
+	md.results.TransientSolution[800-1].Vz,\
+	md.results.TransientSolution[800-1].Vel,\
+	md.results.TransientSolution[800-1].Pressure,\
+	md.results.TransientSolution[800-1].Bed,\
+	md.results.TransientSolution[800-1].Surface,\
+	md.results.TransientSolution[800-1].Thickness,\
+	md.results.TransientSolution[800-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[1200-1].Vx,\
+	md.results.TransientSolution[1200-1].Vy,\
+	md.results.TransientSolution[1200-1].Vz,\
+	md.results.TransientSolution[1200-1].Vel,\
+	md.results.TransientSolution[1200-1].Pressure,\
+	md.results.TransientSolution[1200-1].Bed,\
+	md.results.TransientSolution[1200-1].Surface,\
+	md.results.TransientSolution[1200-1].Thickness,\
+	md.results.TransientSolution[1200-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[1600-1].Vx,\
+	md.results.TransientSolution[1600-1].Vy,\
+	md.results.TransientSolution[1600-1].Vz,\
+	md.results.TransientSolution[1600-1].Vel,\
+	md.results.TransientSolution[1600-1].Pressure,\
+	md.results.TransientSolution[1600-1].Bed,\
+	md.results.TransientSolution[1600-1].Surface,\
+	md.results.TransientSolution[1600-1].Thickness,\
+	md.results.TransientSolution[1600-1].SurfaceforcingsMassBalance,\
+	md.results.TransientSolution[2000-1].Vx,\
+	md.results.TransientSolution[2000-1].Vy,\
+	md.results.TransientSolution[2000-1].Vz,\
+	md.results.TransientSolution[2000-1].Vel,\
+	md.results.TransientSolution[2000-1].Pressure,\
+	md.results.TransientSolution[2000-1].Bed,\
+	md.results.TransientSolution[2000-1].Surface,\
+	md.results.TransientSolution[2000-1].Thickness,\
+	md.results.TransientSolution[2000-1].SurfaceforcingsMassBalance,\
+	]
+
+if printingflag:
+	pass
+	"""
+
+	starttime = 360;
+	endtime = 2000;
+	res = 40;
+	ts = [starttime:res:endtime];
+
+	index = md.mesh.elements;
+	x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3));
+	y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3));
+	areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
+
+	thickness = [];
+	volume = [];
+	massbal = [];
+	velocity = [];
+	for t=starttime:endtime
+		thickness = [thickness (md.results.TransientSolution(t).Thickness)];
+		volume = [volume mean(md.results.TransientSolution(t).Thickness.value,2).*areas];
+		massbal = [massbal (md.results.TransientSolution(t).SurfaceforcingsMassBalance)];
+		velocity = [velocity (md.results.TransientSolution(t).Vel)];
+	end
+
+	figure('Position', [0 0 1060 1060])
+
+	options = plotoptions('data','transient_movie','unit','km');
+	options = options.list{1};
+	options = checkplotoptions(md,options);
+
+	%loop over the time steps
+	results=md.results.TransientSolution;
+	count = 1;
+	for i=ts
+
+		subplot(5,9,[28:31 37:40])
+		set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
+		field = 'Thickness';
+
+		%process data
+		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
+		[data datatype]=processdata(md,results(i).(field),options);
+
+		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
+		options=changefieldvalue(options,'title',titlestring);
+		options=addfielddefault(options,'colorbar',1);
+		options=changefieldvalue(options,'caxis',[0 max(max(thickness))]);
+		applyoptions(md,[],options);
+
+		subplot(5,9,[33:36 42:45])
+		set(gca,'pos',get(gca,'pos')+[-0.01 -0.08 0.07 0.08])
+		field = 'Vel';
+
+		%process data
+		[x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
+		[data datatype]=processdata(md,results(i).(field),options);
+
+		titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
+		plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
+		options=changefieldvalue(options,'title',titlestring);
+		options=addfielddefault(options,'colorbar',1);
+		options=changefieldvalue(options,'caxis',[0 max(max(velocity))]);
+		applyoptions(md,[],options);
+
+		subplot(5,4,1:4)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
+		plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Surface Mass Balance','FontSize',14)
+		ylabel('m/year','FontSize',14)
+
+		subplot(5,4,5:8)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
+		plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Ice Volume','FontSize',14)
+		ylabel('km^3','FontSize',14)
+
+		subplot(5,4,9:12)
+		cla
+		set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
+		plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
+		hold on
+		ya = ylim;
+		plot([i i], ya, 'r', 'LineWidth',6)
+		ylim(ya); xlim([starttime endtime]);
+		title('Mean Velocity','FontSize', 14)
+		ylabel('km/year','FontSize', 14)
+		xlabel('year','FontSize', 14)
+
+		set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
+		if i==starttime,
+			%initialize images and frame
+			frame=getframe(gcf);
+			[images,map]=rgb2ind(frame.cdata,256,'nodither');
+			images(1,1,1,length(ts))=0;
+		else
+			frame=getframe(gcf);
+			images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither');
+		end
+
+		count = count+1;
+
+	end
+
+	filename='transawtooth3d.gif';
+	imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
+	"""
+
Index: /issm/trunk/test/NightlyRun/test207.py
===================================================================
--- /issm/trunk/test/NightlyRun/test207.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test207.py	(revision 14067)
@@ -15,8 +15,8 @@
 md=setflowequation(md,'macayeal','all')
 md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())
 
Index: /issm/trunk/test/NightlyRun/test228.py
===================================================================
--- /issm/trunk/test/NightlyRun/test228.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test228.py	(revision 14067)
@@ -24,5 +24,5 @@
 
 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
 
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test229.py
===================================================================
--- /issm/trunk/test/NightlyRun/test229.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test229.py	(revision 14067)
@@ -24,5 +24,5 @@
 
 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
 
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test230.py
===================================================================
--- /issm/trunk/test/NightlyRun/test230.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test230.py	(revision 14067)
@@ -25,5 +25,5 @@
 
 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
 
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test231.py
===================================================================
--- /issm/trunk/test/NightlyRun/test231.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test231.py	(revision 14067)
@@ -25,5 +25,5 @@
 
 md.surfaceforcings.mass_balance=numpy.vstack((smb,[1.5,3.]))
-md.transient.isthermal=0
+md.transient.isthermal=False
 
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test232.py
===================================================================
--- /issm/trunk/test/NightlyRun/test232.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test232.py	(revision 14067)
@@ -18,8 +18,8 @@
 md.timestepping.time_step=1.
 md.timestepping.final_time=4.
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())
 
Index: /issm/trunk/test/NightlyRun/test3009.py
===================================================================
--- /issm/trunk/test/NightlyRun/test3009.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test3009.py	(revision 14067)
@@ -14,8 +14,8 @@
 md=setflowequation(md,'macayeal','all')
 md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md.autodiff.isautodiff=True
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test313.py
===================================================================
--- /issm/trunk/test/NightlyRun/test313.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test313.py	(revision 14067)
@@ -15,8 +15,8 @@
 md.cluster=generic('name',oshostname(),'np',3)
 md.verbose=verbose('convergence',True,'solution',True)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())
 
Index: /issm/trunk/test/NightlyRun/test326.py
===================================================================
--- /issm/trunk/test/NightlyRun/test326.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test326.py	(revision 14067)
@@ -16,8 +16,8 @@
 md.cluster=generic('name',oshostname(),'np',3)
 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1))
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md.thermal.isenthalpy=1
 md=solve(md,TransientSolutionEnum())
Index: /issm/trunk/test/NightlyRun/test407.py
===================================================================
--- /issm/trunk/test/NightlyRun/test407.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test407.py	(revision 14067)
@@ -15,8 +15,8 @@
 md=setflowequation(md,'pattyn','all')
 md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())
 
Index: /issm/trunk/test/NightlyRun/test423.py
===================================================================
--- /issm/trunk/test/NightlyRun/test423.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test423.py	(revision 14067)
@@ -29,8 +29,8 @@
 md.cluster=generic('name',oshostname(),'np',3)
 
-md.transient.isthermal=0
-md.transient.isprognostic=0
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isthermal=False
+md.transient.isprognostic=False
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
 
 #test different grounding line dynamics.
Index: /issm/trunk/test/NightlyRun/test424.py
===================================================================
--- /issm/trunk/test/NightlyRun/test424.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test424.py	(revision 14067)
@@ -20,6 +20,6 @@
 md.geometry.surface=md.geometry.bed+md.geometry.thickness
 md.surfaceforcings.mass_balance[:]=100.
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
 md.groundingline.migration='AgressiveMigration'
 
Index: /issm/trunk/test/NightlyRun/test425.py
===================================================================
--- /issm/trunk/test/NightlyRun/test425.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test425.py	(revision 14067)
@@ -20,6 +20,6 @@
 md.geometry.surface=md.geometry.bed+md.geometry.thickness
 md.surfaceforcings.mass_balance[:]=-150.
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
 md.groundingline.migration='SoftMigration'
 
Index: /issm/trunk/test/NightlyRun/test426.py
===================================================================
--- /issm/trunk/test/NightlyRun/test426.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test426.py	(revision 14067)
@@ -21,6 +21,6 @@
 md.extrude(3,1.);
 md=setflowequation(md,'macayeal','all');
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
 md.groundingline.migration='AgressiveMigration'
 md.cluster=generic('name',oshostname(),'np',3)
Index: /issm/trunk/test/NightlyRun/test427.py
===================================================================
--- /issm/trunk/test/NightlyRun/test427.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test427.py	(revision 14067)
@@ -22,6 +22,6 @@
 
 md.surfaceforcings.mass_balance[:]=-150
-md.transient.isdiagnostic=0
-md.transient.isgroundingline=1
+md.transient.isdiagnostic=False
+md.transient.isgroundingline=True
 md.groundingline.migration='SoftMigration'
 md.cluster=generic('name',oshostname(),'np',3)
Index: /issm/trunk/test/NightlyRun/test515.py
===================================================================
--- /issm/trunk/test/NightlyRun/test515.py	(revision 14066)
+++ /issm/trunk/test/NightlyRun/test515.py	(revision 14067)
@@ -15,8 +15,8 @@
 md.thermal.stabilization=2
 md.cluster=generic('name',oshostname(),'np',3)
-md.transient.isdiagnostic=0
-md.transient.isprognostic=0
-md.transient.isthermal=1
-md.transient.isgroundingline=0
+md.transient.isdiagnostic=False
+md.transient.isprognostic=False
+md.transient.isthermal=True
+md.transient.isgroundingline=False
 md=solve(md,TransientSolutionEnum())