Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26513)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_sampling.cpp	(revision 26514)
@@ -10,11 +10,43 @@
 #include <random>
 
+#undef M_PI
+#define M_PI 3.141592653589793238462643
+
+double lcg(long long* seed){/*{{{*/
+
+
+	int N = 5;
+	long long M = 1LL<<32;
+	long long a = 1812433253;
+	long long c = 1;
+
+	for(int i=0;i<N;i++){
+		*seed = a * *seed + c;
+		*seed = *seed%M;
+	}
+
+  return ((double) *seed) / ((double) M);
+
+}
+
+double BoxMuller(long long* pseed){
+
+	double u1 = lcg(pseed);
+	double u2 = lcg(pseed);
+
+	double R = sqrt(-2*log(u1));
+	double theta = 2*M_PI*u2;
+
+	return R*cos(theta);
+
+}
+
 void GaussianVector(Vector<IssmDouble>* ppf,int seed){/*{{{*/
+
 
 	/*Intermediaries*/
 	double      rdnumber;
 
-	/*Define random number generator*/
-	std::normal_distribution<double> distribution(0.0,1.0);  // Define probability distribution as the standard normal distribution
+	/*Define seed*/
 	if(seed<0){
 		std::random_device rd;
@@ -27,6 +59,7 @@
 		seed = seed + 783728*my_rank; // change default seed for parallel simulations (by considering an arbitrary shif based on the rank number)
 	}
-	//std::default_random_engine generator(seed);
-	std::mt19937 generator(seed);
+
+	long long seedbis;
+ 	seedbis  = (long long) seed;
 
 	int        *local_indices = NULL;
@@ -37,5 +70,5 @@
 	ppf->GetLocalSize(&M);
 	for(int i=0;i<M;i++){
-		rdnumber = distribution(generator);
+		rdnumber = BoxMuller(&seedbis);
 		ppf->SetValue(local_indices[i],rdnumber,INS_VAL);
 	}