Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26481)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26482)
@@ -3607,31 +3607,29 @@
 	const int numvertices = this->GetNumberOfVertices();
 	int         basinid;
-	IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars 
+	IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin;  
 	IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
    IssmDouble* smbspin     = xNew<IssmDouble>(numvertices*arorder);
 
-	/*Get Basin ID*/
+	/*Get Basin ID and Basin coefficients*/
 	this->GetInputValue(&basinid,SmbBasinsIdEnum);
-
-	for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii];
-
+	for(int i=0;i<arorder;i++) phi_basin[i] = phi[basinid*arorder+i];
 	beta0_basin   = beta0[basinid];
 	beta1_basin   = beta1[basinid];
-	for(int jj=0;jj<nspin;jj++){	
-		tspin = starttime-((nspin-jj)*tstep_ar); //current time in spin-up loop
-      noisespin_basin = noisespin[jj*numbasins+basinid];
+
+	/*Loop over number of spin-up steps for all vertices*/	
+	for(int j=0;j<nspin;j++){	
+		tspin = starttime-((nspin-j)*tstep_ar); 
+      noisespin_basin = noisespin[j*numbasins+basinid];
       IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
-      for(int ii=0;ii<numvertices*arorder;ii++) oldsmbspin[ii]=smbspin[ii]; //copy smbspin in oldsmbspin
+      for(int i=0;i<numvertices*arorder;i++) oldsmbspin[i]=smbspin[i]; 
 
       for(int v=0;v<numvertices;v++){
          IssmDouble autoregressionterm = 0.;
-         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];
+         for(int i=0;i<arorder;i++) autoregressionterm += phi_basin[i]*smbspin[v+i*numvertices];
          smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
       }
 
-		/*correct older values in smbspin*/
-      for(int ii=0;ii<(arorder-1)*numvertices;ii++){
-			smbspin[ii+numvertices]=oldsmbspin[ii];
-		} 
+		/*Correct older values in smbspin*/
+      for(int i=0;i<(arorder-1)*numvertices;i++) smbspin[i+numvertices]=oldsmbspin[i]; 
 
       xDelete<IssmDouble>(oldsmbspin); 
@@ -3644,22 +3642,24 @@
 }/*}}}*/
 void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
+	
 	const int numvertices = this->GetNumberOfVertices();
 	int         basinid,M,N;
-	IssmDouble  beta0_basin,beta1_basin,noise_basin; //initialize scalars 
+	IssmDouble  beta0_basin,beta1_basin,noise_basin;  
 	IssmDouble* phi_basin  = xNew<IssmDouble>(arorder);
    IssmDouble* smblist    = xNew<IssmDouble>(numvertices);
-	IssmDouble* smbvaluesautoregression = NULL; //array for past SMB values that we are about to retrieve
-
+	IssmDouble* smbvaluesautoregression = NULL;
+
+	/*Get Basin ID and Basin coefficients*/
 	this->GetInputValue(&basinid,SmbBasinsIdEnum);
-
-	for(int ii=0;ii<arorder;ii++){phi_basin[ii] = phi[basinid*arorder+ii];}
+	for(int ii=0;ii<arorder;ii++) phi_basin[ii] = phi[basinid*arorder+ii]; 
 	beta0_basin   = beta0[basinid];
 	beta1_basin   = beta1[basinid];
-	noise_basin   = noiseterms[basinid]; //note that noiseterms is computed at every timestep
-	this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M); //get array of past SMB values to compute AR model
+	noise_basin   = noiseterms[basinid]; 
+	this->inputs->GetArray(SmbValuesAutoregressionEnum,this->lid,&smbvaluesautoregression,&M);
+
 	/*If not AR model timestep: take the old SMB values*/
 	if(isstepforar==false){
-		//VV do something with the lapse rate here if needed (12Oct2021)
-      for(int ii=0;ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
+		/*VV do something with the lapse rate here if needed (12Oct2021)*/
+      for(int i=0;i<numvertices;i++) smblist[i]=smbvaluesautoregression[i]; 
 	}
 	/*If AR model timestep: compute SMB values on vertices from AR*/
@@ -3667,21 +3667,18 @@
 		for(int v=0;v<numvertices;v++){
 
-			/*compute autoregressive term*/
+			/*Compute autoregressive term*/
 			IssmDouble autoregressionterm=0.;
-			for(int ii=0;ii<arorder;ii++){
-				autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
-			}
-
-			/*stochastic SMB value*/
+			for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];
+
+			/*Stochastic SMB value*/
 			smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin;
 		}
 		/*Update autoregression smb values*/
 		IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
-		for(int ii=0;ii<numvertices;ii++) temparray[ii] = smblist[ii]; //first store newly computed smb values
-		for(int ii=0;ii<(arorder-1)*numvertices;ii++){
-			temparray[ii+numvertices] = smbvaluesautoregression[ii];
-		} //second shift older smb values
-		this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); //save updated autoregression values
-		xDelete<IssmDouble>(temparray); //cleanup
+		/*Assign newest values and shift older values*/
+		for(int i=0;i<numvertices;i++) temparray[i] = smblist[i]; 
+		for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = smbvaluesautoregression[i];
+		this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,temparray,numvertices*arorder); 
+		xDelete<IssmDouble>(temparray);
 	}
 	/*Add input to element*/
Index: /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp	(revision 26481)
+++ /issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp	(revision 26482)
@@ -60,8 +60,8 @@
 		femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
 		femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
+		oceangridx=xNew<IssmDouble>(ngrids_ocean);
+		oceangridy=xNew<IssmDouble>(ngrids_ocean);
 		if(my_rank==0){
-			oceangridx = xNew<IssmDouble>(ngrids_ocean);
 			ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
-			oceangridy = xNew<IssmDouble>(ngrids_ocean);
 			ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
 
@@ -70,8 +70,5 @@
 			ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
 		}
-		if(my_rank!=0){
-			oceangridx=xNew<IssmDouble>(ngrids_ocean);
-			oceangridy=xNew<IssmDouble>(ngrids_ocean);
-		}
+		
 		ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 		ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26481)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26482)
@@ -148,7 +148,10 @@
 }/*}}}*/
 void SmbautoregressionInitx(FemModel* femmodel){/*{{{*/
+	
 	/*Initialization step of Smbautoregressionx*/
-	int M,N,Nphi,arorder,numbasins;
+	int M,N,Nphi,arorder,numbasins,my_rank;
 	IssmDouble starttime,tstep_ar,tinit_ar;
+	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
+	femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
 	IssmDouble* beta0    = xNew<IssmDouble>(numbasins);
 	IssmDouble* beta1    = xNew<IssmDouble>(numbasins);
@@ -158,19 +161,23 @@
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
 	femmodel->parameters->FindParam(&tinit_ar,SmbAutoregressionInitialTimeEnum);
-	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
-	femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
 	femmodel->parameters->FindParam(&beta0,&M,SmbBeta0Enum);    _assert_(M==numbasins);
 	femmodel->parameters->FindParam(&beta1,&M,SmbBeta1Enum);    _assert_(M==numbasins);
 	femmodel->parameters->FindParam(&phi,&M,&Nphi,SmbPhiEnum);  _assert_(M==numbasins); _assert_(Nphi==arorder);
 	femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
+	
 	/*AR model spin-up*/
-	int nspin{2*arorder+5}; //number of spin-up steps to be executed
-	IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); //sample of basin-specific noise at each spinup step
-	for(int ii{0};ii<nspin;ii++){
-		IssmDouble* temparray = xNew<IssmDouble>(numbasins);
-		multivariateNormal(&temparray,numbasins,0.0,covmat);
-		for(int jj{0};jj<numbasins;jj++){noisespin[ii*numbasins+jj]=temparray[jj];}
-		xDelete<IssmDouble>(temparray);
-	}
+	int nspin{2*arorder+5}; 
+	IssmDouble* noisespin = xNewZeroInit<IssmDouble>(numbasins*nspin); 
+	my_rank=IssmComm::GetRank();
+	if(my_rank==0){
+		for(int i=0;i<nspin;i++){
+			IssmDouble* temparray = xNew<IssmDouble>(numbasins);
+			multivariateNormal(&temparray,numbasins,0.0,covmat);
+			for(int j=0;j<numbasins;j++){noisespin[i*numbasins+j]=temparray[j];}
+			xDelete<IssmDouble>(temparray);
+		}
+	}
+	ISSM_MPI_Bcast(noisespin,numbasins*nspin,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+	
 	/*Initialize SmbValuesAutoregressionEnum*/
 	for(Object* &object:femmodel->elements->objects){
@@ -178,4 +185,5 @@
 		element->SmbautoregressionInit(numbasins,arorder,nspin,starttime,tstep_ar,tinit_ar,beta0,beta1,phi,noisespin);
 	}
+	
 	/*Cleanup*/
 	xDelete<IssmDouble>(beta0);
@@ -192,11 +200,12 @@
 	femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    femmodel->parameters->FindParam(&tstep_ar,SmbAutoregressionTimestepEnum);
+
 	/*Initialize module at first time step*/
 	if(time<=starttime+dt){SmbautoregressionInitx(femmodel);}
 	/*Determine if this is a time step for the AR model*/
-	bool isstepforar{false};
+	bool isstepforar = false;
 
 	#ifndef _HAVE_AD_
-   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt) isstepforar = true;
 	#else
 	_error_("not implemented yet");
@@ -204,5 +213,5 @@
 
 	/*Load parameters*/
-	int M,N,Nphi,arorder,numbasins;
+	int M,N,Nphi,arorder,numbasins,my_rank;
 	femmodel->parameters->FindParam(&numbasins,SmbNumBasinsEnum);
 	femmodel->parameters->FindParam(&arorder,SmbAutoregressiveOrderEnum);
@@ -219,9 +228,13 @@
 	femmodel->parameters->FindParam(&covmat,&M,&N,SmbCovmatEnum); _assert_(M==numbasins); _assert_(N==numbasins);
 
-	/*time elapsed with respect to AR model initial time*/
+	/*Time elapsed with respect to AR model initial time*/
 	IssmDouble telapsed_ar = time-tinit_ar; 
 
 	/*Before looping through elements: compute noise term specific to each basin from covmat*/
-	multivariateNormal(&noiseterms,numbasins,0.0,covmat);
+	my_rank=IssmComm::GetRank();
+	if(my_rank==0){
+		multivariateNormal(&noiseterms,numbasins,0.0,covmat);
+	}
+	ISSM_MPI_Bcast(noiseterms,numbasins,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
 
 	/*Loop over each element to compute SMB at vertices*/
Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26481)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26482)
@@ -679,18 +679,13 @@
 
 	/*ensure zero-initialization*/
-	for(int ii=0;ii<(ndim*ndim);ii++) Lchol[ii]=0;;
-
-	for(int ii=0;ii<ndim;ii++){
-		for(int jj{0};jj<=ii;jj++){
+	for(int i=0;i<(ndim*ndim);i++) Lchol[i]=0;;
+
+	for(int i=0;i<ndim;i++){
+		for(int j=0;j<=i;j++){
 			IssmDouble sum=0.;
-			for(int kk{0};kk<jj;kk++){
-				sum += Lchol[ii*ndim+kk]*Lchol[jj*ndim+kk];
-			}
-			if(ii==jj){
-				Lchol[ii*ndim+jj] = sqrt(A[ii*ndim+jj]-sum);
-			}
-			else{
-				Lchol[ii*ndim+jj] = 1./Lchol[jj*ndim+jj] * (A[ii*ndim+jj]-sum);
-			}
+			for(int k=0;k<j;k++) sum += Lchol[i*ndim+k]*Lchol[j*ndim+k];
+			
+			if(i==j) Lchol[i*ndim+j] = sqrt(A[i*ndim+j]-sum);
+			else Lchol[i*ndim+j]     = 1./Lchol[j*ndim+j] * (A[i*ndim+j]-sum);
 		}
 	}
Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26481)
+++ /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26482)
@@ -21,24 +21,27 @@
 
 void univariateNormal(IssmPDouble* prand, IssmPDouble mean, IssmPDouble sdev) { /*{{{*/
-	/*univariateNormal	generates a random value follwoing Normal distribution*/
-	unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); //random seed using time_since_epoch
-   std::default_random_engine generator(seed); //generator of random numbers
-   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); //Normal probability distribution
+
+	/*Random seed using time_since_epoch*/
+	unsigned seed = std::chrono::steady_clock::now().time_since_epoch().count(); 
+   std::default_random_engine generator(seed);
+	/*Normal Probability Distribution*/
+   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); 
 	*prand = normdistri(generator);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
-   IssmPDouble* sampleStandardNormal     = xNew<IssmPDouble>(dim);
+   
+	IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
    IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
    IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
 
-   for(int ii=0;ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
+   for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0); 
    CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-   for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
-      /*entry-by-entry multiplication along matrix row*/
+   
+	/*Matrix by vector multiplication*/
+	for(int i=0;i<dim;i++){ 
+      /*Entry-by-entry multiplication along matrix row*/
       IssmDouble sum=0.;
-      for(int jj{0};jj<dim;jj++){
-         sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
-      }
-      sampleMultivariateNormal[ii] = mean+sum;
+      for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j]; 
+      sampleMultivariateNormal[i] = mean+sum;
    }
 
@@ -49,18 +52,20 @@
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
-	IssmPDouble* sampleStandardNormal     = xNew<IssmPDouble>(dim);
+	
+	IssmPDouble* sampleStandardNormal    = xNew<IssmPDouble>(dim);
 	IssmDouble* sampleMultivariateNormal = xNew<IssmDouble>(dim);
 	IssmDouble* Lchol                    = xNewZeroInit<IssmDouble>(dim*dim);
-	for(int ii=0;ii<dim;ii++) univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);
+	for(int i=0;i<dim;i++) univariateNormal(&(sampleStandardNormal[i]),0.0,1.0);
 
 	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
 
-	for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+	/*Matrix by vector multiplication*/
+	for(int i=0;i<dim;i++){
 		IssmDouble sum = 0.;
-      for(int jj=0.;jj<dim;jj++){
-			sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
-		}
-      sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
+      for(int j=0;j<dim;j++) sum += sampleStandardNormal[j]*Lchol[i*dim+j]; 
+      sampleMultivariateNormal[i] = mean[i]+sum;
 	}
+   
+	/*Assign output pointer and cleanup*/
 	*prand = sampleMultivariateNormal;
 	xDelete<IssmPDouble>(sampleStandardNormal);