Index: /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp	(revision 26611)
+++ /issm/trunk-jpl/src/c/shared/Random/randomgenerator.cpp	(revision 26612)
@@ -13,38 +13,67 @@
 namespace rnd{
 
+  /* Linear congruential engine */
+
+  linear_congruential_engine::linear_congruential_engine(){/*{{{*/
+
+      pseed = new unsigned int;
+			*pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1);
+      a = 1103515245;		       // BSD Formula
+      c = 12345;					     // BSD Formula
+      m = 2147483648;			     // BSD Formula
+			return;
+	}/*}}}*/
+
+  linear_congruential_engine::linear_congruential_engine(unsigned int _a, unsigned int _b, unsigned int _m){/*{{{*/
+
+      pseed = new unsigned int;
+			*pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1);
+      a = _a;
+      c = _b;
+      m = _m;
+			return;
+	}/*}}}*/
+
+  linear_congruential_engine::~linear_congruential_engine(){}
+
+  unsigned int linear_congruential_engine::get_m() { return m; }
+
+  void linear_congruential_engine::seed( int s ) {
+      if(s<0)
+        *pseed = std::chrono::steady_clock::now().time_since_epoch()/std::chrono::milliseconds(1);
+      else
+        *pseed = (unsigned) s;
+  }
+
+  unsigned int linear_congruential_engine::generator(){
+    *pseed = ( a * *pseed + c ) % m ;
+    return *pseed;
+  }
+
+  /* Uniform distribution */
+
 	uniform_distribution::uniform_distribution(){/*{{{*/
-
-			a   = 1103515245; 	// BSD Formula
-			c  = 12345;					// BSD Formula
-			m = 2147483648;			// BSD Formula
-			_seed = 0;
-			lbound = 0.0;
-			ubound = 1.0;
+			a = 0.0;
+			b = 1.0;
 			return;
 	}
 	/*}}}*/
-	uniform_distribution::uniform_distribution(double lower,double upper){/*{{{*/
 
-			a   = 1103515245;		// BSD Formula
-			c  = 12345;					// BSD Formula
-			m = 2147483648;			// BSD Formula
-			_seed = 0;
-			lbound = lower;
-			ubound = upper;
+	uniform_distribution::uniform_distribution(double _a,double _b){/*{{{*/
+			a = _a;
+			b = _b;
 			return;
 	}
 	/*}}}*/
+
 	uniform_distribution::~uniform_distribution(){}
-	void uniform_distribution::seed( unsigned int s ) { _seed = s; }
-	unsigned int uniform_distribution::get_seed() { return _seed; }
-	double uniform_distribution::generator() {
-			_seed = ( a * _seed + c ) % m ;
-			return (ubound-lbound)*(double) _seed/ m + lbound;
+
+	double uniform_distribution::generator(rnd::linear_congruential_engine random_engine) {
+			return (b-a)*(double) random_engine.generator()/ random_engine.get_m() + a;
 	}
 
+  /* Normal distribution */
 
 	normal_distribution::normal_distribution(){/*{{{*/
-
-			_seed = 0;
 			mean   = 0;
 			sdev  = 1.0;
@@ -52,26 +81,23 @@
 	}
 	/*}}}*/
+
 	normal_distribution::normal_distribution(double m,double s){/*{{{*/
-
-			_seed = 0;
 			mean   = m;
 			sdev  = s;
 			return;
 	}
-		/*}}}*/
+	/*}}}*/
+
 	normal_distribution::~normal_distribution(){}
-	void normal_distribution::seed( unsigned int s ) { _seed = s; }
-	double normal_distribution::generator(){/*{{{*/
+
+	double normal_distribution::generator(rnd::linear_congruential_engine random_engine){/*{{{*/
 
 			rnd::uniform_distribution	unifdistri;
-			unifdistri.seed(_seed);
 
-			double u1 = unifdistri.generator();
-			double u2 = unifdistri.generator();
+			double u1 = unifdistri.generator(random_engine);
+			double u2 = unifdistri.generator(random_engine);
 
 			double R = sqrt(-2*log(u1));
 			double theta = 2*M_PI*u2;
-
-			seed(unifdistri.get_seed());
 
 			return mean + sdev * (R*cos(theta));
Index: /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h	(revision 26611)
+++ /issm/trunk-jpl/src/c/shared/Random/randomgenerator.h	(revision 26612)
@@ -11,14 +11,31 @@
 namespace rnd{
 
+  class linear_congruential_engine
+  {
+    private:
+      unsigned int a;
+      unsigned int c;
+      unsigned int m;
+      unsigned int *pseed;
+
+    public:
+
+      /*constructors, destructors: */
+      linear_congruential_engine();
+      linear_congruential_engine(unsigned int _a, unsigned int _b, unsigned int _m);
+      ~linear_congruential_engine();
+
+      unsigned int get_m();
+      void seed( int s );
+      unsigned int generator();
+
+  };
+
   class uniform_distribution
   {
 
     private:
-      int a;
-      int c;
-      unsigned int m;
-      unsigned _seed;
-      double lbound;
-      double ubound;
+      double a;  // lower bound of range
+      double b;  // upper bound of range
 
     public:
@@ -26,10 +43,8 @@
       /*constructors, destructors: */
       uniform_distribution();
-      uniform_distribution(double a_1,double a_2);
+      uniform_distribution(double _a,double _b);
       ~uniform_distribution();
 
-      void seed( unsigned int s );
-      unsigned int get_seed();
-      double generator();
+      double generator(rnd::linear_congruential_engine random_engine);
 
   };
@@ -39,5 +54,4 @@
 
     private:
-      unsigned _seed;
       double mean;
       double sdev;
@@ -50,10 +64,10 @@
       ~normal_distribution();
 
-      void seed( unsigned int s );
-      double generator();
+      double generator(rnd::linear_congruential_engine random_engine);
 
   };
 
+
 }
 
-#endif //ifndef _RANDOMGENERATOR_H_
+#endif //* _RANDOMGENERATOR_H_ */
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26611)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26612)
@@ -9,5 +9,5 @@
 #include "../analyses/analyses.h"
 #include "../shared/Random/randomgenerator.h"
-#include <random>
+//#include <random>
 
 void GaussianVector(Vector<IssmDouble>* ppf,int seed){/*{{{*/
@@ -18,10 +18,6 @@
 
 	/*Define seed*/
-	rnd::normal_distribution distribution;
-	if(seed<0){
-		std::random_device rd;
-		seed = rd();
-	}
-	else
+	rnd::linear_congruential_engine random_engine;
+	if(seed>=0)
 	{
 		int my_rank;
@@ -29,5 +25,9 @@
 		seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
 	}
-	distribution.seed(seed);
+	random_engine.seed(seed);
+
+	/* Define univariate distribution */
+
+	rnd::normal_distribution distribution(0.0,1.0);
 
 	int        *local_indices = NULL;
@@ -38,5 +38,5 @@
 	ppf->GetLocalSize(&M);
 	for(int i=0;i<M;i++){
-		rdnumber = distribution.generator();
+		rdnumber = distribution.generator(random_engine);
 		ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
 	}