Index: issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
===================================================================
--- issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 26242)
+++ issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 26272)
@@ -80,4 +80,5 @@
 
 	int         nl;
+	int         ntimesteps;
 	IssmDouble* love_h=NULL;
 	IssmDouble* love_k=NULL;
@@ -86,4 +87,6 @@
 	IssmDouble* love_tk=NULL;
 	IssmDouble* love_tl=NULL;
+	IssmDouble* love_timefreq=NULL;
+	bool        love_istime=true;
 	int         externalnature=0;
 	int         isexternal=0;
@@ -91,17 +94,28 @@
 	IssmDouble* G_rigid = NULL;
 	IssmDouble* G_rigid_local = NULL;
-	IssmDouble* G_elastic = NULL;
-	IssmDouble* G_elastic_local = NULL;
-	IssmDouble* U_elastic = NULL;
-	IssmDouble* U_elastic_local = NULL;
-	IssmDouble* H_elastic = NULL;
-	IssmDouble* H_elastic_local = NULL;
+	IssmDouble* G_viscoelastic = NULL;
+	IssmDouble* G_viscoelastic_interpolated = NULL;
+	IssmDouble* G_viscoelastic_local = NULL;
+	IssmDouble* U_viscoelastic = NULL;
+	IssmDouble* U_viscoelastic_interpolated = NULL;
+	IssmDouble* U_viscoelastic_local = NULL;
+	IssmDouble* H_viscoelastic = NULL;
+	IssmDouble* H_viscoelastic_interpolated= NULL;
+	IssmDouble* H_viscoelastic_local = NULL;
 	int         M,m,lower_row,upper_row;
 	IssmDouble  degacc=.01;
+	IssmDouble  timeacc;
 	IssmDouble  planetradius=0;
 	IssmDouble  planetarea=0;
 	bool		rigid=false;
 	bool		elastic=false;
+	bool		viscous=false;
 	bool		rotation=false;
+	int         ndeg;
+	int         horiz;
+
+	bool istime=true;
+	IssmDouble start_time,final_time;
+	int  nt,precomputednt;
 
 	int     numoutputs;
@@ -120,5 +134,5 @@
 	IssmDouble*  bslcocean_partition=NULL;
 	int          npartice,nparthydro,npartocean,nel;
-
+	int          grdmodel;
 
 	/*some constant parameters: */
@@ -134,4 +148,5 @@
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.horiz",SolidearthSettingsHorizEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.elastic",SolidearthSettingsElasticEnum));
+	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.viscous",SolidearthSettingsViscousEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.rotation",SolidearthSettingsRotationEnum));
 	parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
@@ -211,192 +226,308 @@
 	}
 
-	/*Deal with elasticity {{{*/
-	iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
-	iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
-	iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
-
-	if(elastic | rigid){
-		/*compute green functions for a range of angles*/
-		iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
-		M=reCast<int,IssmDouble>(180./degacc+1.);
-	}
-
-	/*love numbers: */
-	if(elastic){
-		iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
-		iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
-		iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l");
-		iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th");
-		iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
-		iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
-
-		parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
-		parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1));
-		parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1));
-		parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
-
-		// AD performance is sensitive to calls to ensurecontiguous.
-		// // Providing "t" will cause ensurecontiguous to be called.
-		#ifdef _HAVE_AD_
-		G_elastic=xNew<IssmDouble>(M,"t");
-		U_elastic=xNew<IssmDouble>(M,"t");
-		H_elastic=xNew<IssmDouble>(M,"t");
-		#else
-		G_elastic=xNew<IssmDouble>(M);
-		U_elastic=xNew<IssmDouble>(M);
-		H_elastic=xNew<IssmDouble>(M);
-		#endif
-	}
-	if(rigid){
-		#ifdef _HAVE_AD_
-		G_rigid=xNew<IssmDouble>(M,"t");
-		#else
-		G_rigid=xNew<IssmDouble>(M);
-		#endif
-	}
-	
-	if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
-
-	if(rigid | elastic){
-
-		/*compute combined legendre + love number (elastic green function:*/
-		m=DetermineLocalSize(M,IssmComm::GetComm());
-		GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
-	}
-	if(elastic){
-		#ifdef _HAVE_AD_
-		G_elastic_local=xNew<IssmDouble>(m,"t");
-		U_elastic_local=xNew<IssmDouble>(m,"t");
-		H_elastic_local=xNew<IssmDouble>(m,"t");
-		#else
-		G_elastic_local=xNew<IssmDouble>(m);
-		U_elastic_local=xNew<IssmDouble>(m);
-		H_elastic_local=xNew<IssmDouble>(m);
-		#endif
-	}
-	if(rigid){
-		#ifdef _HAVE_AD_
-		G_rigid_local=xNew<IssmDouble>(m,"t");
-		#else
-		G_rigid_local=xNew<IssmDouble>(m);
-		#endif
-	}
-
-	if(rigid){
-		for(int i=lower_row;i<upper_row;i++){
-			IssmDouble alpha,x;
-			alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
-			G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
-		}
-	}
-	if(elastic){
-		for(int i=lower_row;i<upper_row;i++){
-			IssmDouble alpha,x;
-			alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
-
-			G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
-			U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
-			H_elastic_local[i-lower_row]= 0; 
-			IssmDouble Pn = 0.; 
-			IssmDouble Pn1 = 0.; 
-			IssmDouble Pn2 = 0.; 
-			IssmDouble Pn_p = 0.; 
-			IssmDouble Pn_p1 = 0.; 
-			IssmDouble Pn_p2 = 0.; 
-
-			for (int n=0;n<nl;n++) {
-				IssmDouble deltalove_G;
-				IssmDouble deltalove_U;
-
-				deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
-				deltalove_U = (love_h[n]-love_h[nl-1]);
-
-				/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
-				if(n==0){
-					Pn=1; 
-					Pn_p=0; 
+	parameters->FindParam(&grdmodel,GrdModelEnum);
+	if(grdmodel!=IvinsEnum){
+		/*Deal with elasticity {{{*/
+		iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
+		iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
+		iomodel->FetchData(&viscous,"md.solidearth.settings.viscous");
+		iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
+		iomodel->FetchData(&horiz,"md.solidearth.settings.horiz");
+
+		if(rigid){
+			/*compute green functions for a range of angles*/
+			iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
+			M=reCast<int,IssmDouble>(180./degacc+1.);
+		}
+
+		/*love numbers: */
+		if(viscous || elastic){
+			int dummy;
+
+			iomodel->FetchData(&timeacc,"md.solidearth.settings.timeacc");
+			iomodel->FetchData(&start_time,"md.timestepping.start_time");
+			iomodel->FetchData(&final_time,"md.timestepping.final_time");
+			iomodel->FetchData(&love_istime,"md.solidearth.lovenumbers.istime");
+			iomodel->FetchData(&love_timefreq,&precomputednt,&dummy,"md.solidearth.lovenumbers.timefreq");
+			iomodel->FetchData(&love_h,&ndeg,&precomputednt,"md.solidearth.lovenumbers.h");
+			iomodel->FetchData(&love_k,&ndeg,&precomputednt,"md.solidearth.lovenumbers.k");
+			iomodel->FetchData(&love_l,&ndeg,&precomputednt,"md.solidearth.lovenumbers.l");
+			iomodel->FetchData(&love_th,&ndeg,&precomputednt,"md.solidearth.lovenumbers.th");
+			iomodel->FetchData(&love_tk,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tk");
+			iomodel->FetchData(&love_tl,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tl");
+
+			parameters->AddObject(new DoubleParam(SolidearthSettingsTimeAccEnum,timeacc));
+			parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
+			parameters->AddObject(new DoubleMatParam(LoveTimeFreqEnum,love_timefreq,precomputednt,1));
+			parameters->AddObject(new BoolParam(LoveIsTimeEnum,love_istime));
+
+			// AD performance is sensitive to calls to ensurecontiguous.
+			// // Providing "t" will cause ensurecontiguous to be called.
+			if(viscous){
+				IssmDouble* stacktimes=NULL;
+				ntimesteps=precomputednt; 
+				nt=reCast<int>((final_time-start_time)/timeacc)+1;
+
+				parameters->AddObject(new IntParam(StackNumStepsEnum,nt));
+				/*Initialize stack times:*/
+				stacktimes=xNew<IssmDouble>(nt);
+				for(int t=0;t<nt;t++){
+					stacktimes[t]=start_time+timeacc*t;
 				}
-				else if(n==1){ 
-					Pn = cos(alpha); 
-					Pn_p = 1; 
+				parameters->AddObject(new DoubleVecParam(StackTimesEnum,stacktimes,nt));
+				parameters->AddObject(new IntParam(StackIndexEnum,0));
+				xDelete<IssmDouble>(stacktimes);
+			}
+			else {
+				ntimesteps=1;
+				nt=1;
+			}
+
+#ifdef _HAVE_AD_
+			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
+			U_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
+			if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
+#else
+			G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
+			U_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
+			if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
+#endif
+		}
+		if(rigid){
+#ifdef _HAVE_AD_
+			G_rigid=xNew<IssmDouble>(M,"t");
+#else
+			G_rigid=xNew<IssmDouble>(M);
+#endif
+		}
+
+		if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
+
+		if(rigid){
+
+			/*compute combined legendre + love number (elastic green function:*/
+			m=DetermineLocalSize(M,IssmComm::GetComm());
+			GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
+		}
+		if(viscous | elastic){
+#ifdef _HAVE_AD_
+			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
+			U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
+			if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
+#else
+			G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
+			U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
+			if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
+#endif
+		}
+		if(rigid){
+#ifdef _HAVE_AD_
+			G_rigid_local=xNew<IssmDouble>(m,"t");
+#else
+			G_rigid_local=xNew<IssmDouble>(m);
+#endif
+		}
+
+		if(rigid){
+			for(int i=lower_row;i<upper_row;i++){
+				IssmDouble alpha,x;
+				alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+				G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
+			}
+		}
+		if(viscous | elastic){
+			for(int i=lower_row;i<upper_row;i++){
+				IssmDouble alpha,x;
+				alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+
+				for(int t=0;t<ntimesteps;t++){
+					G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
+					U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_rigid_local[i-lower_row];
+					if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t]= 0; 
 				}
-				else{
-					Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
-					Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
+
+				IssmDouble Pn = 0.; 
+				IssmDouble Pn1 = 0.; 
+				IssmDouble Pn2 = 0.; 
+				IssmDouble Pn_p = 0.; 
+				IssmDouble Pn_p1 = 0.; 
+				IssmDouble Pn_p2 = 0.; 
+
+				for (int n=0;n<ndeg;n++) {
+
+					/*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
+					if(n==0){
+						Pn=1; 
+						Pn_p=0; 
+					}
+					else if(n==1){ 
+						Pn = cos(alpha); 
+						Pn_p = 1; 
+					}
+					else{
+						Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
+						Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
+					}
+					Pn2=Pn1; Pn1=Pn;
+					Pn_p2=Pn_p1; Pn_p1=Pn_p;
+
+					for(int t=0;t<ntimesteps;t++){
+						IssmDouble deltalove_G;
+						IssmDouble deltalove_U;
+
+						deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
+						deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
+
+						G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_G*Pn;		                // gravitational potential 
+						U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_U*Pn;		                // vertical (up) displacement 
+						if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += sin(alpha)*love_l[n*precomputednt+t]*Pn_p;		// horizontal displacements 
+					}
 				}
-				Pn2=Pn1; Pn1=Pn;
-				Pn_p2=Pn_p1; Pn_p1=Pn_p;
-
-				G_elastic_local[i-lower_row] += deltalove_G*Pn;		// gravitational potential 
-				U_elastic_local[i-lower_row] += deltalove_U*Pn;		// vertical (up) displacement 
-				H_elastic_local[i-lower_row] += sin(alpha)*love_l[n]*Pn_p;		// horizontal displacements 
-			}
-		}
-	}
-	if(rigid){
-
-		/*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
-		int* recvcounts=xNew<int>(IssmComm::GetSize());
-		int* displs=xNew<int>(IssmComm::GetSize());
-
-		//recvcounts:
-		ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
-
-		/*displs: */
-		ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
-
-		/*All gather:*/
-		ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		if(elastic){
-			ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-			ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-			ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-		}
-		
-		/*free resources: */
-		xDelete<int>(recvcounts);
-		xDelete<int>(displs);
-
-		/*Avoid singularity at 0: */
-		G_rigid[0]=G_rigid[1];
-		if(elastic){
-			G_elastic[0]=G_elastic[1];
-			U_elastic[0]=U_elastic[1];
-			H_elastic[0]=H_elastic[1];
-		}
-		
-
-		/*Save our precomputed tables into parameters*/
-		parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M));
-		if(elastic){
-			parameters->AddObject(new DoubleVecParam(SealevelchangeGElasticEnum,G_elastic,M));
-			parameters->AddObject(new DoubleVecParam(SealevelchangeUElasticEnum,U_elastic,M));
-			parameters->AddObject(new DoubleVecParam(SealevelchangeHElasticEnum,H_elastic,M));
-		}
-
-		/*free resources: */
-		xDelete<IssmDouble>(G_rigid);
-		xDelete<IssmDouble>(G_rigid_local);
-		if(elastic){
-			xDelete<IssmDouble>(love_h);
-			xDelete<IssmDouble>(love_k);
-			xDelete<IssmDouble>(love_l);
-			xDelete<IssmDouble>(love_th);
-			xDelete<IssmDouble>(love_tk);
-			xDelete<IssmDouble>(love_tl);
-			xDelete<IssmDouble>(G_elastic);
-			xDelete<IssmDouble>(G_elastic_local);
-			xDelete<IssmDouble>(U_elastic);
-			xDelete<IssmDouble>(U_elastic_local);
-			xDelete<IssmDouble>(H_elastic);
-			xDelete<IssmDouble>(H_elastic_local);
-		}
-	} /*}}}*/
-
-	/*Indicate we have not yet run the Geometry Core module: */
-	parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));
-	/*}}}*/
+			}
+		}
+		if(rigid){
+
+			/*merge G_viscoelastic_local into G_viscoelastic; U_viscoelastic_local into U_viscoelastic; H_viscoelastic_local to H_viscoelastic:{{{*/
+			int* recvcounts=xNew<int>(IssmComm::GetSize());
+			int* displs=xNew<int>(IssmComm::GetSize());
+			int  rc;
+			int  offset;
+
+			//deal with rigid first: 
+			ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
+
+			/*displs: */
+			ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
+
+			/*All gather:*/
+			ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+
+			if(elastic){
+				rc=m*ntimesteps;
+				offset=lower_row*ntimesteps;
+				ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
+				ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
+
+				ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+				ISSM_MPI_Allgatherv(U_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, U_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+				if(horiz)ISSM_MPI_Allgatherv(H_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, H_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+			}
+
+			/*free resources: */
+			xDelete<int>(recvcounts);
+			xDelete<int>(displs);
+
+			
+
+			/*Avoid singularity at 0: */
+			G_rigid[0]=G_rigid[1];
+			if(elastic){
+				for(int t=0;t<ntimesteps;t++){
+					G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
+					U_viscoelastic[t]=U_viscoelastic[ntimesteps+t];
+					if(horiz)H_viscoelastic[t]=H_viscoelastic[ntimesteps+t];
+				}
+			}
+
+			
+
+			/*Reinterpolate viscoelastic green kernels onto a regular gridded time 
+			 *with steps equal to timeacc:*/
+			if(viscous){
+				nt=reCast<int>((final_time-start_time)/timeacc)+1;
+#ifdef _HAVE_AD_
+				G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
+				U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
+				if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
+#else
+				G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
+				U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
+				if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
+#endif
+
+				for(int t=0;t<nt;t++){
+					IssmDouble lincoeff;
+					IssmDouble viscoelastic_time=t*timeacc;
+					int        timeindex2=-1;
+
+					/*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
+					if(t!=0){
+						for(int t2=0;t2<ntimesteps;t2++){
+							if (viscoelastic_time<love_timefreq[t2]){
+								timeindex2=t2-1;
+								if(timeindex2<0)_error_("Temporal Love numbers are computed  with a time accuracy superior to the requested solution time step!");
+								lincoeff=(viscoelastic_time-love_timefreq[t2-1])/(love_timefreq[t2]-love_timefreq[t2-1]);
+								break;
+							}
+						}
+						if(timeindex2==-1)_error_("Temporal love numbers should be extended in time to encompass the requested solution time interval!");
+					}
+					else{
+						timeindex2=0;
+						lincoeff=0;
+					}
+
+					for(int index=0;index<M;index++){
+
+						int timeindex=index*nt+t;
+						int timepreindex= index*ntimesteps+timeindex2;
+						G_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*G_viscoelastic[timepreindex]+lincoeff*G_viscoelastic[timepreindex+1];
+						U_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*U_viscoelastic[timepreindex]+lincoeff*U_viscoelastic[timepreindex+1];
+						if(horiz)H_viscoelastic_interpolated[timeindex]+=(1-lincoeff)*H_viscoelastic[timepreindex]+lincoeff*H_viscoelastic[timepreindex+1];
+					}
+				}
+						
+			}
+			else if(elastic){
+				nt=1; //in elastic, or if we run only rigid, we need only one step
+#ifdef _HAVE_AD_
+				G_viscoelastic_interpolated=G_viscoelastic;
+				U_viscoelastic_interpolated=U_viscoelastic;
+				if(horiz)H_viscoelastic_interpolated=H_viscoelastic;
+#else
+				G_viscoelastic_interpolated=G_viscoelastic;
+				U_viscoelastic_interpolated=U_viscoelastic;
+				if(horiz)H_viscoelastic_interpolated=H_viscoelastic;
+#endif
+			}	
+
+			/*Save our precomputed tables into parameters*/
+			parameters->AddObject(new DoubleVecParam(SealevelchangeGRigidEnum,G_rigid,M));
+			if(viscous || elastic){
+				parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
+				parameters->AddObject(new DoubleVecParam(SealevelchangeUViscoElasticEnum,U_viscoelastic_interpolated,M*nt));
+				if(horiz)parameters->AddObject(new DoubleVecParam(SealevelchangeHViscoElasticEnum,H_viscoelastic_interpolated,M*nt));
+			}
+
+			/*free resources: */
+			xDelete<IssmDouble>(G_rigid);
+			xDelete<IssmDouble>(G_rigid_local);
+			if(elastic){
+				xDelete<IssmDouble>(love_h);
+				xDelete<IssmDouble>(love_k);
+				xDelete<IssmDouble>(love_l);
+				xDelete<IssmDouble>(love_th);
+				xDelete<IssmDouble>(love_tk);
+				xDelete<IssmDouble>(love_tl);
+				xDelete<IssmDouble>(G_viscoelastic);
+				xDelete<IssmDouble>(G_viscoelastic_local);
+				xDelete<IssmDouble>(U_viscoelastic);
+				xDelete<IssmDouble>(U_viscoelastic_local);
+				if(horiz){
+					xDelete<IssmDouble>(H_viscoelastic);
+					xDelete<IssmDouble>(H_viscoelastic_local);
+				}
+			}
+		} /*}}}*/
+
+		/*Indicate we have not yet run the Geometry Core module: */
+		parameters->AddObject(new BoolParam(SealevelchangeGeometryDoneEnum,false));
+		/*}}}*/
+	}
 
 	/*Transitions:{{{ */
