Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26620)
+++ /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26621)
@@ -9,8 +9,6 @@
 #include <math.h>
 #include <float.h>    /*  DBL_EPSILON  */
-#include <chrono>
 #include <cstdarg>
 #include <iostream>
-#include <random>
 
 #include "../Matrix/matrix.h"
@@ -18,19 +16,17 @@
 #include "../MemOps/MemOps.h"
 #include "../io/io.h"
+#include "./randomgenerator.h"
 /*}}}*/
 
-void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seedfixed=-1) { /*{{{*/
+void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev, int seed=-1) { /*{{{*/
 
-	unsigned seed;
-	/*Random seed using time_since_epoch*/
-	if(seedfixed==-1) seed = std::chrono::steady_clock::now().time_since_epoch().count(); 
-	/*Seed fixed by input argument*/
-	else seed = seedfixed;
-	std::default_random_engine generator(seed);
-	/*Normal Probability Distribution*/
-   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); 
-	*prand = normdistri(generator);
+	/*Seed the pseudo-random number generator*/
+	rnd::linear_congruential_engine randomengine;
+	randomengine.seed(seed);
+	/*Normal distribution*/
+	rnd::normal_distribution distriNormal(mean,sdev);
+	*prand = distriNormal.generator(randomengine);
 } /*}}}*/
-void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
    
 	IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
@@ -38,6 +34,15 @@
    IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
 
-	/*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
-	for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed); 
+	/*True randomness if seed<0, otherwise random seed is fixed at seed*/
+	/*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
+	rnd::linear_congruential_engine randomengine;
+	randomengine.seed(seed);
+	/*Normal distribution*/
+	rnd::normal_distribution distriNormal(0.0,1.0);
+	for(int i=0;i<dim;i++){
+		sampleStandardNormal[i] = distriNormal.generator(randomengine);
+	}
+
+	/*Cholsesky decomposition of the covariance matrix*/
 	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
    
@@ -55,12 +60,21 @@
    xDelete<IssmDouble>(Lchol);
 } /*}}}*/
-void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1) { /*{{{*/
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1) { /*{{{*/
 	
 	IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
 	IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
 	IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
+	
+	/*True randomness if seed<0, otherwise random seed is fixed at seed*/
+	/*Seed the pseudo-random number generator, repeatedly calling univariateNormal does not ensure randomness*/
+	rnd::linear_congruential_engine randomengine;
+	randomengine.seed(seed);
+	/*Normal distribution*/
+	rnd::normal_distribution distriNormal(0.0,1.0);
+	for(int i=0;i<dim;i++){
+		sampleStandardNormal[i] = distriNormal.generator(randomengine);
+	}
 
-	/*True randomness if seedfixed==-1, otherwise random seed is fixed at seedfixed*/
-	for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0,seedfixed); 
+	/*Cholsesky decomposition of the covariance matrix*/
 	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
 
Index: /issm/trunk-jpl/src/c/shared/Random/random.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.h	(revision 26620)
+++ /issm/trunk-jpl/src/c/shared/Random/random.h	(revision 26621)
@@ -6,7 +6,7 @@
 #define _RANDOM_H_
 
-void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev, int seedfixed=-1);
-void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seedfixed=-1);
-void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seedfixed=-1);
+void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev, int seed=-1);
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix, int seed=-1);
+void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix, int seed=-1);
 
 #endif //ifndef _RANDOM_H_