Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26478)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 26479)
@@ -3604,29 +3604,42 @@
 /*}}}*/
 void       Element::SmbautoregressionInit(int numbasins,int arorder,int nspin,IssmDouble starttime,IssmDouble tstep_ar,IssmDouble tinit_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noisespin){/*{{{*/
+
 	const int numvertices = this->GetNumberOfVertices();
 	int         basinid;
 	IssmDouble  tspin,beta0_basin,beta1_basin,noisespin_basin; //initialize scalars 
-	IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
-   IssmDouble* smbspin                 = xNew<IssmDouble>(numvertices*arorder);
+	IssmDouble* phi_basin   = xNew<IssmDouble>(arorder);
+   IssmDouble* smbspin     = xNew<IssmDouble>(numvertices*arorder);
+
+	/*Get Basin ID*/
 	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];
-	for(int jj{0};jj<nspin;jj++){	
+	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];
-      IssmDouble* oldsmbspin = xNewZeroInit<IssmDouble>(numvertices*arorder);
-      for(int ii{0};ii<numvertices*arorder;ii++){oldsmbspin[ii]=smbspin[ii];} //copy smbspin in oldsmbspin
+      IssmDouble* oldsmbspin = xNew<IssmDouble>(numvertices*arorder);
+      for(int ii=0;ii<numvertices*arorder;ii++) oldsmbspin[ii]=smbspin[ii]; //copy smbspin in oldsmbspin
+
       for(int v=0;v<numvertices;v++){
-         double autoregressionterm{0.0};
-         for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];}
+         IssmDouble autoregressionterm = 0.;
+         for(int ii=0;ii<arorder;ii++) autoregressionterm += phi_basin[ii]*smbspin[v+ii*numvertices];
          smbspin[v] = beta0_basin+beta1_basin*(tspin-tinit_ar)+autoregressionterm+noisespin_basin; //compute newest values in smbspin
       }
-      for(int ii{0};ii<(arorder-1)*numvertices;ii++){smbspin[ii+numvertices]=oldsmbspin[ii];} //correct older values in smbspin
-      xDelete<IssmDouble>(oldsmbspin); //cleanup
-	}
-	this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);//save spin-up autoregression values
-   xDelete<IssmDouble>(smbspin); //cleanup
-   xDelete<IssmDouble>(phi_basin); //cleanup
+
+		/*correct older values in smbspin*/
+      for(int ii=0;ii<(arorder-1)*numvertices;ii++){
+			smbspin[ii+numvertices]=oldsmbspin[ii];
+		} 
+
+      xDelete<IssmDouble>(oldsmbspin); 
+	}
+	this->inputs->SetArrayInput(SmbValuesAutoregressionEnum,this->lid,smbspin,numvertices*arorder);
+
+	/*Cleanup and return*/
+   xDelete<IssmDouble>(smbspin); 
+   xDelete<IssmDouble>(phi_basin);
 }/*}}}*/
 void       Element::Smbautoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble* beta0,IssmDouble* beta1,IssmDouble* phi,IssmDouble* noiseterms){/*{{{*/
@@ -3634,9 +3647,11 @@
 	int         basinid,M,N;
 	IssmDouble  beta0_basin,beta1_basin,noise_basin; //initialize scalars 
-	IssmDouble* phi_basin               = xNew<IssmDouble>(arorder);
-   IssmDouble* smblist                 = xNew<IssmDouble>(numvertices);
+	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
+
 	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];
@@ -3646,17 +3661,25 @@
 	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];}
+      for(int ii=0;ii<numvertices;ii++){smblist[ii]=smbvaluesautoregression[ii];}
 	}
 	/*If AR model timestep: compute SMB values on vertices from AR*/
 	else{
 		for(int v=0;v<numvertices;v++){
-			double autoregressionterm{0.0};
-			for(int ii{0};ii<arorder;ii++){autoregressionterm += phi_basin[ii]*smbvaluesautoregression[v+ii*numvertices];} //compute autoregressive term
-			smblist[v] = beta0_basin+beta1_basin*telapsed_ar+autoregressionterm+noise_basin; //stochastic SMB value
+
+			/*compute autoregressive term*/
+			IssmDouble autoregressionterm=0.;
+			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
+		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
@@ -3664,4 +3687,5 @@
 	/*Add input to element*/
 	this->AddInput(SmbMassBalanceEnum,smblist,P1Enum);
+
 	/*Cleanup*/
 	xDelete<IssmDouble>(phi_basin);
Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26478)
+++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp	(revision 26479)
@@ -3,4 +3,5 @@
  */
 
+#include <config.h>
 #include "./SurfaceMassBalancex.h"
 #include "../../shared/shared.h"
@@ -195,5 +196,11 @@
 	/*Determine if this is a time step for the AR model*/
 	bool isstepforar{false};
-   if((std::fmod(time,tstep_ar)<std::fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+
+	#ifndef _HAVE_AD_
+   if((fmod(time,tstep_ar)<fmod((time-dt),tstep_ar)) || (time<=starttime+dt) || tstep_ar==dt){isstepforar = true;}
+	#else
+	_error_("not implemented yet");
+	#endif
+
 	/*Load parameters*/
 	int M,N,Nphi,arorder,numbasins;
Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26478)
+++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp	(revision 26479)
@@ -677,14 +677,21 @@
    Follows the Cholesky–Banachiewicz algorithm
    Lchol should point to an IssmDouble* of same dimensions as A*/
-   for(int ii{0};ii<(ndim*ndim);ii++){Lchol[ii]=0;}; //ensure zero-initialization
-   for(int ii{0};ii<ndim;ii++){
-                for(int jj{0};jj<=ii;jj++){
-                        double sum{0.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);}
-                }
-        }
+
+	/*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++){
+			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);
+			}
+		}
+	}
 } /*}}}*/
Index: /issm/trunk-jpl/src/c/shared/Random/random.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26478)
+++ /issm/trunk-jpl/src/c/shared/Random/random.cpp	(revision 26479)
@@ -20,40 +20,49 @@
 /*}}}*/
 
-void univariateNormal(IssmDouble* prand, IssmDouble mean, IssmDouble sdev) { /*{{{*/
+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<double> normdistri(mean,sdev); //Normal probability distribution
-	double tfunc;
+   std::normal_distribution<IssmPDouble> normdistri(mean,sdev); //Normal probability distribution
 	*prand = normdistri(generator);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble mean, IssmDouble* covariancematrix) { /*{{{*/
-   IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(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 ii=0;ii<dim;ii++){univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);}
    CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-   for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
-      double sum{0.0};
-      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
-      sampleMultivariateNormal[ii] = mean+sum; //assign value
+   for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+      /*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;
    }
+
+   /*Assign output pointer and cleanup*/
    *prand = sampleMultivariateNormal;
-   xDelete<IssmDouble>(sampleStandardNormal);
+   xDelete<IssmPDouble>(sampleStandardNormal);
    xDelete<IssmDouble>(Lchol);
 } /*}}}*/
 void multivariateNormal(IssmDouble** prand, int dim, IssmDouble* mean, IssmDouble* covariancematrix) { /*{{{*/
-	IssmDouble* sampleStandardNormal     = xNew<IssmDouble>(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 ii=0;ii<dim;ii++) univariateNormal(&(sampleStandardNormal[ii]),0.0,1.0);
+
 	CholeskyRealPositiveDefinite(Lchol,covariancematrix,dim);
-	for(int ii{0};ii<dim;ii++){ //matrix by vector multiplication
-		double sum{0.0};
-      for(int jj{0};jj<dim;jj++){sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];} //entry-by-entry multiplication along matrix row
+
+	for(int ii=0;ii<dim;ii++){ //matrix by vector multiplication
+		IssmDouble sum = 0.;
+      for(int jj=0.;jj<dim;jj++){
+			sum += sampleStandardNormal[jj]*Lchol[ii*dim+jj];
+		}
       sampleMultivariateNormal[ii] = mean[ii]+sum; //assign value
 	}
 	*prand = sampleMultivariateNormal;
-	xDelete<IssmDouble>(sampleStandardNormal);
+	xDelete<IssmPDouble>(sampleStandardNormal);
 	xDelete<IssmDouble>(Lchol);
 } /*}}}*/