Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp	(revision 26271)
+++ /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:{{{ */
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 26272)
@@ -401,4 +401,5 @@
 		virtual void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
 		virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
+		virtual void       SealevelchangeUpdateViscousFields()=0;
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 26272)
@@ -234,4 +234,5 @@
 		void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
+		void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
 		#endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 26272)
@@ -189,4 +189,5 @@
 		void       SealevelchangeDeformationConvolution(GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
+		void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 26272)
@@ -196,4 +196,5 @@
 		void       SealevelchangeDeformationConvolution(GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
 		void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
+		void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
 #endif
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 26272)
@@ -6100,4 +6100,10 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 
+	/*stacks:*/
+	IssmDouble* stackRSL = NULL;
+	IssmDouble* stackU = NULL;
+	IssmDouble* stackN = NULL;
+	IssmDouble* stackE = NULL;
+
 	#ifdef _HAVE_RESTRICT_
 	IssmDouble* __restrict__ G=NULL;
@@ -6105,8 +6111,8 @@
 	IssmDouble* __restrict__ GN=NULL;
 	IssmDouble* __restrict__ GE=NULL;
-	IssmDouble* __restrict__ G_elastic_precomputed=NULL;
+	IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
 	IssmDouble* __restrict__ G_rigid_precomputed=NULL;
-	IssmDouble* __restrict__ U_elastic_precomputed=NULL;
-	IssmDouble* __restrict__ H_elastic_precomputed=NULL;
+	IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
+	IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
 	#else
 	IssmDouble* G=NULL;
@@ -6114,11 +6120,11 @@
 	IssmDouble* GN=NULL;
 	IssmDouble* GE=NULL;
-	IssmDouble* G_elastic_precomputed=NULL;
+	IssmDouble* G_viscoelastic_precomputed=NULL;
 	IssmDouble* G_rigid_precomputed=NULL;
-	IssmDouble* U_elastic_precomputed=NULL;
-	IssmDouble* H_elastic_precomputed=NULL;
+	IssmDouble* U_viscoelastic_precomputed=NULL;
+	IssmDouble* H_viscoelastic_precomputed=NULL;
 	#endif
 	
-	/*elastic green function:*/
+	/*viscoelastic green function:*/
 	int index;
 	int         M;
@@ -6128,5 +6134,10 @@
 	bool computeelastic = false;
 	bool computerotation = false;
+	bool computeviscous = false;
 	int  horiz;
+	bool istime=true;
+	IssmDouble timeacc=0;
+	IssmDouble start_time,final_time;
+	int  nt,precomputednt;
 
 	/*Rotational:*/
@@ -6149,4 +6160,5 @@
 	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
 	this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum);
+	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
 	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
 	this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
@@ -6173,14 +6185,14 @@
 
 	if(computeelastic){
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
-
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
 		
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
+
+		if(horiz){
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
+		}
 	}
 	/*}}}*/
@@ -6194,13 +6206,24 @@
 	/*}}}*/
 	/*Compute green functions:{{{ */
+	if(computeviscous){
+		this->parameters->FindParam(&istime,LoveIsTimeEnum);
+		if(!istime)_error_("Frequency love numbers not supported yet!");
+		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
+		this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
+		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
+		nt=reCast<int>((final_time-start_time)/timeacc)+1;
+	}
+	else{
+		nt=1; //in elastic, or if we run only rigid, we need only one step
+	}
+
 	constant=3/rho_earth/planetarea;
-	
-	//Allocate: 
-	G=xNewZeroInit<IssmDouble>(3*nel);
+
+	G=xNewZeroInit<IssmDouble>(3*nel*nt);
 	if(computeelastic){
-		GU=xNewZeroInit<IssmDouble>(3*nel);
+		GU=xNewZeroInit<IssmDouble>(3*nel*nt);
 		if(horiz){
-			GN=xNewZeroInit<IssmDouble>(3*nel);
-			GE=xNewZeroInit<IssmDouble>(3*nel);
+			GN=xNewZeroInit<IssmDouble>(3*nel*nt);
+			GE=xNewZeroInit<IssmDouble>(3*nel*nt);
 		}
 	}
@@ -6224,43 +6247,49 @@
 			_assert_(index>0 && index<M);
 
-			/*Rigid earth gravitational perturbation: */
-			G[i*nel+e]+=G_rigid_precomputed[index];
-			
-			if(computeelastic){
-				G[i*nel+e]+=G_elastic_precomputed[index];
-			}
-			G[i*nel+e]=G[i*nel+e]*ratioe;
-
-			/*Elastic components:*/
-			if(computeelastic){
-				GU[i*nel+e] =  ratioe * U_elastic_precomputed[index];
-				if(horiz){
-					/*Compute azimuths, both north and east components: */
-					x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
-					if(lati==M_PI/2){
-						x=1e-12; y=1e-12;
+			if(horiz){
+				/*Compute azimuths, both north and east components: */
+				x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
+				if(lati==M_PI/2){
+					x=1e-12; y=1e-12;
+				}
+				if(lati==-M_PI/2){
+					x=1e-12; y=1e-12;
+				}
+				dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
+				N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+				E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+			}
+
+			for(int t=0;t<nt;t++){
+				int timeindex=i*nel*nt+e*nt+t; 
+
+				/*Rigid earth gravitational perturbation: */
+				G[timeindex]+=G_rigid_precomputed[index];
+
+				/*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
+				if(computeelastic){
+					G[timeindex]+=G_viscoelastic_precomputed[index*nt+t];
+				}
+				G[timeindex]=G[timeindex]*ratioe;
+
+				/*Elastic components:*/
+				if(computeelastic){
+					GU[timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
+					if(horiz){
+						GN[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
+						GE[timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
 					}
-					if(lati==-M_PI/2){
-						x=1e-12; y=1e-12;
-					}
-					dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
-					N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-					E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-
-					GN[i*nel+e] = ratioe*H_elastic_precomputed[index]*N_azim;
-					GE[i*nel+e] = ratioe*H_elastic_precomputed[index]*E_azim;
 				}
-			}
-		}
-	}
-
+			} //for(int t=0;t<nt;t++)
+		} //for (int i=0;i<3;i++)
+	} //for(int e=0;e<nel;e++)
 
 	/*Add in inputs:*/
-	this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3);
+	this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3*nt);
 	if(computeelastic){
-		this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3);
+		this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3*nt);
 		if(horiz){
-			this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3);
-			this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3);
+			this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3*nt);
+			this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3*nt);
 		}
 	}
@@ -6346,4 +6375,20 @@
 	}
 	/*}}}*/
+	/*Initialize stacks: {{{*/
+	if(computeviscous){
+		stackRSL=xNewZeroInit<IssmDouble>(3*nt);
+		stackU=xNewZeroInit<IssmDouble>(3*nt);
+
+		this->inputs->SetArrayInput(StackRSLEnum,this->lid,stackRSL,3*nt);
+		this->inputs->SetArrayInput(StackUEnum,this->lid,stackU,3*nt);
+		if(horiz){
+			stackN=xNewZeroInit<IssmDouble>(3*nt);
+			stackE=xNewZeroInit<IssmDouble>(3*nt);
+			this->inputs->SetArrayInput(StackNEnum,this->lid,stackRSL,3*nt);
+			this->inputs->SetArrayInput(StackEEnum,this->lid,stackU,3*nt);
+		}
+	}
+	/*}}}*/
+
 	/*Free allocations:{{{*/
 	#ifdef _HAVE_RESTRICT_
@@ -6692,8 +6737,8 @@
 
 	#ifdef _HAVE_RESTRICT_
-	IssmDouble* __restrict__ G_elastic_precomputed=NULL;
+	IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
 	IssmDouble* __restrict__ G_rigid_precomputed=NULL;
-	IssmDouble* __restrict__ U_elastic_precomputed=NULL;
-	IssmDouble* __restrict__ H_elastic_precomputed=NULL;
+	IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
+	IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
 	IssmDouble** __restrict__ Gsubel=NULL;
 	IssmDouble** __restrict__ GUsubel=NULL;
@@ -6702,8 +6747,8 @@
 
 	#else
-	IssmDouble* G_elastic_precomputed=NULL;
+	IssmDouble* G_viscoelastic_precomputed=NULL;
 	IssmDouble* G_rigid_precomputed=NULL;
-	IssmDouble* U_elastic_precomputed=NULL;
-	IssmDouble* H_elastic_precomputed=NULL;
+	IssmDouble* U_viscoelastic_precomputed=NULL;
+	IssmDouble* H_viscoelastic_precomputed=NULL;
 	IssmDouble** Gsubel=NULL;
 	IssmDouble** GUsubel=NULL;
@@ -6719,5 +6764,11 @@
 	bool computerigid = false;
 	bool computeelastic = false;
+	bool computeviscous = false;
 	int  horiz;
+
+	bool istime=true;
+	IssmDouble timeacc=0;
+	IssmDouble start_time,final_time;
+	int  nt,precomputednt;
 
 	/*}}}*/
@@ -6726,4 +6777,5 @@
 	this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
 	this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
+	this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
 	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
 	this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
@@ -6740,12 +6792,15 @@
 
 	if(computeelastic){
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
-
-		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUElasticEnum)); _assert_(parameter);
-		parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
+
+		parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
+		parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
+
+		if(horiz){
+			parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
+			parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
+
+		}
 	}
 	/*}}}*/
@@ -6760,4 +6815,15 @@
 	constant=3/rho_earth/planetarea;
 
+	if(computeviscous){
+		this->parameters->FindParam(&istime,LoveIsTimeEnum);
+		if(!istime)_error_("Frequency love numbers not supported yet!");
+		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
+		this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
+		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
+		nt=reCast<int>((final_time-start_time)/timeacc)+1;
+	}
+	else{
+		nt=1; //in elastic, or if we run only rigid, we need only one step
+	}
 	Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
 	if(computeelastic){
@@ -6772,10 +6838,10 @@
 	for(int l=0;l<SLGEOM_NUMLOADS;l++){
 		int nbar=slgeom->nbar[l];
-		Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
+		Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
 		if(computeelastic){
-			GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
+			GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
 			if(horiz){
-				GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
-				GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar);
+				GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
+				GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
 			}
 		}
@@ -6795,4 +6861,22 @@
 				longi=longitude[i];
 
+				if(horiz){
+					/*Compute azimuths, both north and east components: */
+					x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
+					if(lati==90){
+						x=1e-12; y=1e-12;
+					}
+					if(lati==-90){
+						x=1e-12; y=1e-12;
+					}
+					IssmDouble xbar=planetradius*cos(late)*cos(longe);
+					IssmDouble ybar=planetradius*cos(late)*sin(longe);
+					IssmDouble zbar=planetradius*sin(late);
+
+					dx = xbar-x; dy = ybar-y; dz = zbar-z;
+					N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+					E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+				}
+
 				/*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
 				delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
@@ -6800,34 +6884,24 @@
 				index=reCast<int,IssmDouble>( alpha/M_PI*reCast<IssmDouble,int>(M-1) );
 
-				/*Rigid earth gravitational perturbation: */
-				Gsubel[l][i*nbar+e]+=G_rigid_precomputed[index];
-
-				if(computeelastic){
-					Gsubel[l][i*nbar+e]+=G_elastic_precomputed[index];
-				}
-				Gsubel[l][i*nbar+e]*=ratioe;
-
-				/*Elastic components:*/
-				if(computeelastic){
-					GUsubel[l][i*nbar+e] =  ratioe * U_elastic_precomputed[index];
-					if(horiz){
-						/*Compute azimuths, both north and east components: */
-						x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
-						if(lati==90){
-							x=1e-12; y=1e-12;
+
+				for(int t=0;t<nt;t++){
+					int timeindex=i*nbar*nt+e*nt+t; 
+
+					/*Rigid earth gravitational perturbation: */
+					Gsubel[l][timeindex]+=G_rigid_precomputed[index];
+
+					if(computeelastic){
+						Gsubel[l][timeindex]+=G_viscoelastic_precomputed[index*nt+t];
+					}
+					Gsubel[l][timeindex]*=ratioe;
+
+					/*Elastic components:*/
+					if(computeelastic){
+						GUsubel[l][timeindex] =  ratioe * U_viscoelastic_precomputed[index*nt+t];
+
+						if(horiz){
+							GNsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*N_azim;
+							GEsubel[l][timeindex] = ratioe*H_viscoelastic_precomputed[index*nt+t]*E_azim;
 						}
-						if(lati==-90){
-							x=1e-12; y=1e-12;
-						}
-						IssmDouble xbar=planetradius*cos(late)*cos(longe);
-						IssmDouble ybar=planetradius*cos(late)*sin(longe);
-						IssmDouble zbar=planetradius*sin(late);
-
-						dx = xbar-x; dy = ybar-y; dz = zbar-z;
-						N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-						E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-
-						GNsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*N_azim;
-						GEsubel[l][i*nbar+e] = ratioe*H_elastic_precomputed[index]*E_azim;
 					}
 				}
@@ -6838,10 +6912,10 @@
 	/*Save all these arrayins for each element:*/
 	for (int l=0;l<SLGEOM_NUMLOADS;l++){
-		this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3);
+		this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3*nt);
 		if(computeelastic){
-			this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3);
+			this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3*nt);
 			if(horiz){
-				this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3);
-				this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3);
+				this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3*nt);
+				this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3*nt);
 			}
 		}
@@ -6869,4 +6943,69 @@
 	/*}}}*/
 	return;
+
+}
+/*}}}*/
+void       Tria::SealevelchangeUpdateViscousFields(){ /*{{{*/
+	
+	/*Inputs:*/
+	IssmDouble* stackRSL=NULL;
+	IssmDouble* stackU=NULL;
+	IssmDouble* stackN=NULL;
+	IssmDouble* stackE=NULL;
+	IssmDouble* stacktimes=NULL;
+	int         stacknumsteps;
+	int         stackindex; 
+	int         newindex; 
+	int         dummy;
+	bool        viscous=false;
+	IssmDouble  currenttime;
+	IssmDouble  lincoeff=0;
+	int horiz;
+		
+	this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
+	
+	if(viscous){
+		this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
+		
+		this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
+		this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
+		this->parameters->FindParam(&stackindex,StackIndexEnum);
+		this->parameters->FindParam(&currenttime,TimeEnum);
+
+		this->inputs->GetArrayPtr(StackRSLEnum,this->lid,&stackRSL,&dummy);
+		this->inputs->GetArrayPtr(StackUEnum,this->lid,&stackU,&dummy);
+		if(horiz){
+			this->inputs->GetArrayPtr(StackNEnum,this->lid,&stackN,&dummy);
+			this->inputs->GetArrayPtr(StackEEnum,this->lid,&stackE,&dummy);
+		}
+
+		lincoeff=0;
+		newindex=stacknumsteps-2;
+		for(int t=stackindex;t<stacknumsteps;t++){
+			if (stacktimes[t]>currenttime){
+				newindex=t-1;
+				lincoeff=(currenttime-stacktimes[newindex])/(stacktimes[t]-stacktimes[newindex]);
+				break;
+			}
+		}
+		if(newindex==(stacknumsteps-2))lincoeff=1;
+		stacktimes[newindex]=currenttime;
+		for(int i=0;i<NUMVERTICES;i++){
+			stackRSL[i*stacknumsteps+newindex]=(1-lincoeff)*stackRSL[i*stacknumsteps+newindex]+lincoeff*stackRSL[i*stacknumsteps+newindex+1];
+			stackU[i*stacknumsteps+newindex]=(1-lincoeff)*stackU[i*stacknumsteps+newindex]+lincoeff*stackU[i*stacknumsteps+newindex+1];
+			if(horiz){
+				stackN[i*stacknumsteps+newindex]=(1-lincoeff)*stackN[i*stacknumsteps+newindex]+lincoeff*stackN[i*stacknumsteps+newindex+1];
+				stackE[i*stacknumsteps+newindex]=(1-lincoeff)*stackE[i*stacknumsteps+newindex]+lincoeff*stackE[i*stacknumsteps+newindex+1];
+			}
+		}
+		stackindex=newindex;
+
+		this->parameters->SetParam(stackindex,StackIndexEnum);
+		this->parameters->SetParam(stacktimes,stacknumsteps,StackTimesEnum);
+
+		/*free allocations:*/
+		xDelete<IssmDouble>(stacktimes);
+	}
+
 
 }
@@ -6959,6 +7098,8 @@
 	IssmDouble oceanarea=0;
 	IssmDouble rho_water;
+	bool computefuture=false;
 	
 	bool sal = false;
+	bool viscous = false;
 	bool rotation= false;
 	bool percpu= false;
@@ -6968,6 +7109,7 @@
 	IssmDouble Grotm2[3];
 	IssmDouble Grotm3[3];
-	
+		
 	this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
+	this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
 	this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
 	this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
@@ -6980,5 +7122,5 @@
 		this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
 
-		this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true);
+		this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,StackRSLEnum,computefuture=false);
 	}
 
@@ -7051,4 +7193,5 @@
 	IssmDouble* GNsub[SLGEOM_NUMLOADS];
 	IssmDouble* GEsub[SLGEOM_NUMLOADS];
+	bool computefuture=false;
 
 	int horiz;
@@ -7102,12 +7245,11 @@
 			}
 		}
-
-		this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false);
+		this->SealevelchangeGxL(&SealevelRSL[0],G, Gsub, loads, slgeom, nel,percpu=false,StackRSLEnum,computefuture=true);
 
 		if(elastic){
-			this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false);
+			this->SealevelchangeGxL(&SealevelU[0],GU, GUsub, loads, slgeom, nel,percpu=false,StackUEnum,computefuture=true);
 			if(horiz ){
-				this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false);
-				this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false);
+				this->SealevelchangeGxL(&SealevelN[0],GN, GNsub, loads, slgeom, nel,percpu=false,StackNEnum,computefuture=true);
+				this->SealevelchangeGxL(&SealevelE[0],GE, GEsub, loads, slgeom, nel,percpu=false,StackEEnum,computefuture=true);
 			}
 		}
@@ -7166,23 +7308,49 @@
 
 } /*}}}*/
-void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu) { /*{{{*/
-
-	IssmDouble sealevel[3]={0};
-	int i,e,l,nbar;
+void       Tria::SealevelchangeGxL(IssmDouble* sealevelout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int stackenum, bool computefuture) { /*{{{*/
+
+	IssmDouble* sealevel=NULL;
+	int i,e,l,t,nbar;
+	bool viscous=false;
+	IssmDouble* stack=NULL;
+	int nt=1; //important
+	int stackindex=0; //important
+	int stacknumsteps=1; //important
+
+	this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
+	if(viscous){
+		this->parameters->FindParam(&stacknumsteps,StackNumStepsEnum);
+		if(computefuture) nt=stacknumsteps;
+		else nt=1;
+
+		//allocate
+		sealevel=xNewZeroInit<IssmDouble>(3*nt);
+	}
+	else sealevel=xNewZeroInit<IssmDouble>(3*nt);
 
 	if(loads->sealevelloads){
+
 		for(i=0;i<NUMVERTICES;i++) {
 			if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
 			for (e=0;e<nel;e++){
-				sealevel[i]+=G[i*nel+e]*(loads->sealevelloads[e]+loads->loads[e]);
+				for(t=0;t<nt;t++){
+					int index=i*nel*stacknumsteps+e*stacknumsteps+t;
+					sealevel[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
+				}
 			}
 			for(l=0;l<SLGEOM_NUMLOADS;l++){
 				nbar=slgeom->nbar[l];
 				for (e=0;e<nbar;e++){
-					sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
+					for(t=0;t<nt;t++){
+						int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
+						sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
+					}
 				}
 				if(l==SLGEOM_OCEAN){
 					for (e=0;e<nbar;e++){
-						sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subsealevelloads[e]);
+						for(t=0;t<nt;t++){
+							int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
+							sealevel[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
+						}
 					}
 				}
@@ -7191,16 +7359,73 @@
 	}
 	else{  //this is the initial convolution where only loads are provided
+
 		for(i=0;i<NUMVERTICES;i++) {
 			if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
 			for (e=0;e<nel;e++){
-				sealevel[i]+=G[i*nel+e]*(loads->loads[e]);
+				for(t=0;t<nt;t++){
+					int index=i*nel*stacknumsteps+e*stacknumsteps+t;
+					sealevel[i*nt+t]+=G[index]*(loads->loads[e]);
+				}
 			}
 			for(l=0;l<SLGEOM_NUMLOADS;l++){
 				nbar=slgeom->nbar[l];
 				for (e=0;e<nbar;e++){
-					sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subloads[l][e]);
+					for(t=0;t<nt;t++){
+						int index=i*nbar*stacknumsteps+e*stacknumsteps+t;
+						sealevel[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
+					}
 				}
 			}
 		}
+	}
+	
+	if(viscous){
+		IssmDouble* sealevelinterp=NULL;
+		IssmDouble* stacktimes=NULL;
+		IssmDouble  final_time;
+		IssmDouble  lincoeff;
+		IssmDouble  timeacc;
+
+		this->parameters->FindParam(&stackindex,StackIndexEnum);
+		this->parameters->FindParam(&stacktimes,NULL,StackTimesEnum);
+		this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
+		this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
+		this->inputs->GetArrayPtr(stackenum,this->lid,&stack,NULL);
+		if(computefuture){
+			sealevelinterp=xNew<IssmDouble>(3*nt); 
+			if(stacktimes[stackindex]<final_time){
+				lincoeff=1-(stacktimes[stackindex+1]-stacktimes[stackindex])/timeacc;
+				for(int t=stackindex;t<nt;t++){
+					if(t==stackindex){
+						for(int i=0;i<NUMVERTICES;i++){
+							sealevelinterp[i*nt+t]=  sealevel[i*nt+0];
+						}
+					}
+					else{
+						for(int i=0;i<NUMVERTICES;i++){
+							sealevelinterp[i*nt+t]=  (1-lincoeff)*sealevel[i*nt+(t-stackindex-1)]+lincoeff*sealevel[i*nt+(t-stackindex)];
+						}
+					}
+				}
+			}
+		}
+
+		/*update sealevel at present time using stack at present time: */
+		for(int i=0;i<NUMVERTICES;i++){
+			sealevel[i*nt+0]+=stack[i*stacknumsteps+stackindex]; 
+		}
+
+		if(computefuture){ /*update stack with future deformation from present load: */
+			for(int t=stackindex;t<nt;t++){
+				for(int i=0;i<NUMVERTICES;i++){
+					stack[i*stacknumsteps+t]+=sealevelinterp[i*nt+t];
+				}
+			}
+			/*Re-add stack now that we updated:*/
+			this->inputs->SetArrayInput(stackenum,this->lid,stack,3*stacknumsteps);
+		}
+
+		/*Free allocatoins:*/
+		xDelete<IssmDouble>(stacktimes);
 	}
 
@@ -7208,9 +7433,9 @@
 	if(percpu){
 		for(i=0;i<NUMVERTICES;i++){
-			if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i];
+			if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i*nt+0];
 		}
 	}
 	else{
-		for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i];
+		for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i*nt+0];
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 26272)
@@ -180,4 +180,5 @@
 		void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom);
 		void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom);
+		void       SealevelchangeUpdateViscousFields();
 		#endif
 		/*}}}*/
@@ -243,5 +244,5 @@
 		void           UpdateConstraintsExtrudeFromBase(void);
 		void           UpdateConstraintsExtrudeFromTop(void);
-		void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool optimize);
+		void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
 		/*}}}*/
 
Index: /issm/trunk-jpl/src/c/cores/love_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26271)
+++ /issm/trunk-jpl/src/c/cores/love_core.cpp	(revision 26272)
@@ -1019,4 +1019,5 @@
 	doubletype*  LoveHf = NULL;
 	doubletype*  LoveLf = NULL;
+	
 	doubletype*  LoveKernels= NULL;
 	LoveVariables* vars=NULL;
@@ -1104,14 +1105,46 @@
 	//Temporal love numbers
 	if (istemporal && !complex_computation){
+
+		IssmDouble*  LoveHtDouble=NULL;
+		IssmDouble*  LoveKtDouble=NULL;
+		IssmDouble*  LoveLtDouble=NULL;
+		doubletype*  LoveHt=NULL;
+		doubletype*  LoveLt=NULL;
+		doubletype*  LoveKt=NULL;
+
 		int NTit;
 		femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
 		int nt = nfreq/2/NTit;
 
-		doubletype* LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
-		doubletype* LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
-		doubletype* LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+		LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+		LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+		LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
 
 		love_freq_to_temporal<doubletype>(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
 
+		/*Downcast and add into parameters:*/
+		LoveHtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
+		LoveLtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
+		LoveKtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
+		for(int i=0;i<(sh_nmax+1)*nt;i++){
+			LoveHtDouble[i]=std::real(LoveHt[i]);
+			LoveLtDouble[i]=std::real(LoveLt[i]);
+			LoveKtDouble[i]=std::real(LoveKt[i]);
+		}
+		if(forcing_type==9){ //tidal loading
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
+		}
+		else if(forcing_type==11){ //surface loading
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
+		}
+		xDelete<IssmDouble>(LoveHtDouble);
+		xDelete<IssmDouble>(LoveKtDouble);
+		xDelete<IssmDouble>(LoveLtDouble);
+	
+		/*Add into external results, no need to downcast, we can handle complexes/quad precision: */
 		femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
 		femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
@@ -1122,5 +1155,35 @@
 		xDelete<doubletype>(LoveKt);
 	}
-
+	else{
+
+		IssmDouble*  LoveHfDouble=NULL;
+		IssmDouble*  LoveKfDouble=NULL;
+		IssmDouble*  LoveLfDouble=NULL;
+
+		/*Add into parameters:*/
+		LoveHfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
+		LoveLfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
+		LoveKfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
+		for(int i=0;i<(sh_nmax+1)*nfreq;i++){
+			LoveHfDouble[i]=std::real(LoveHf[i]);
+			LoveLfDouble[i]=std::real(LoveLf[i]);
+			LoveKfDouble[i]=std::real(LoveKf[i]);
+		}
+		if(forcing_type==9){ //tidal loading
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
+		}
+		else if(forcing_type==11){ //surface loading
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
+			femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
+		}
+		xDelete<IssmDouble>(LoveHfDouble);
+		xDelete<IssmDouble>(LoveKfDouble);
+		xDelete<IssmDouble>(LoveLfDouble);
+	}
+
+	/*Add into external results:*/
 	femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
 	femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));
Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26271)
+++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp	(revision 26272)
@@ -245,4 +245,5 @@
 	int  grdmodel; 
 	int  computesealevel=0;
+	bool viscous=false;
 	IssmDouble*           sealevelpercpu=NULL;
 
@@ -259,4 +260,5 @@
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
 	femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
+	femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
 	/*}}}*/
 
@@ -296,4 +298,10 @@
 
 	if(VerboseSolution()) _printf0_("	  starting  GRD convolutions\n");
+	
+	/*update viscous RSL:*/
+	for(Object* & object : femmodel->elements->objects){
+		Element* element = xDynamicCast<Element*>(object);
+		element->SealevelchangeUpdateViscousFields();
+	}
 
 	/*buildup loads: */
Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26271)
+++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim	(revision 26272)
@@ -364,4 +364,5 @@
 syn keyword cConstant RotationalAngularVelocityEnum
 syn keyword cConstant SolidearthSettingsElasticEnum
+syn keyword cConstant SolidearthSettingsViscousEnum
 syn keyword cConstant SealevelchangeGeometryDoneEnum
 syn keyword cConstant RotationalEquatorialMoiEnum
@@ -373,10 +374,13 @@
 syn keyword cConstant LoadLoveKEnum
 syn keyword cConstant LoadLoveLEnum
+syn keyword cConstant LoveTimeFreqEnum
+syn keyword cConstant LoveIsTimeEnum
 syn keyword cConstant SealevelchangeGRigidEnum
-syn keyword cConstant SealevelchangeGElasticEnum
+syn keyword cConstant SealevelchangeGViscoElasticEnum
 syn keyword cConstant SolidearthSettingsComputesealevelchangeEnum
 syn keyword cConstant SolidearthSettingsGRDEnum
 syn keyword cConstant SolidearthSettingsRunFrequencyEnum
-syn keyword cConstant SealevelchangeHElasticEnum
+syn keyword cConstant SolidearthSettingsTimeAccEnum
+syn keyword cConstant SealevelchangeHViscoElasticEnum
 syn keyword cConstant SolidearthSettingsHorizEnum
 syn keyword cConstant SolidearthSettingsMaxiterEnum
@@ -389,5 +393,5 @@
 syn keyword cConstant SealevelchangeRunCountEnum
 syn keyword cConstant SealevelchangeTransitionsEnum
-syn keyword cConstant SealevelchangeUElasticEnum
+syn keyword cConstant SealevelchangeUViscoElasticEnum
 syn keyword cConstant SettingsIoGatherEnum
 syn keyword cConstant SettingsNumResultsOnNodesEnum
@@ -398,4 +402,7 @@
 syn keyword cConstant SettingsSolverResidueThresholdEnum
 syn keyword cConstant SettingsWaitonlockEnum
+syn keyword cConstant StackNumStepsEnum
+syn keyword cConstant StackTimesEnum
+syn keyword cConstant StackIndexEnum
 syn keyword cConstant SmbAIceEnum
 syn keyword cConstant SmbAIdxEnum
@@ -949,4 +956,8 @@
 syn keyword cConstant SolidearthExternalGeoidRateEnum
 syn keyword cConstant SolidearthExternalBarystaticSeaLevelRateEnum
+syn keyword cConstant StackRSLEnum
+syn keyword cConstant StackUEnum
+syn keyword cConstant StackNEnum
+syn keyword cConstant StackEEnum
 syn keyword cConstant StrainRateeffectiveEnum
 syn keyword cConstant StrainRateparallelEnum
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26271)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h	(revision 26272)
@@ -358,4 +358,5 @@
 	RotationalAngularVelocityEnum,
 	SolidearthSettingsElasticEnum,
+	SolidearthSettingsViscousEnum,
 	SealevelchangeGeometryDoneEnum,
 	RotationalEquatorialMoiEnum,
@@ -367,10 +368,13 @@
 	LoadLoveKEnum,
 	LoadLoveLEnum,
+	LoveTimeFreqEnum,
+	LoveIsTimeEnum,
 	SealevelchangeGRigidEnum,
-	SealevelchangeGElasticEnum,
+	SealevelchangeGViscoElasticEnum,
 	SolidearthSettingsComputesealevelchangeEnum,
 	SolidearthSettingsGRDEnum,
 	SolidearthSettingsRunFrequencyEnum,
-	SealevelchangeHElasticEnum,
+	SolidearthSettingsTimeAccEnum,
+	SealevelchangeHViscoElasticEnum,
 	SolidearthSettingsHorizEnum,
 	SolidearthSettingsMaxiterEnum,
@@ -383,5 +387,5 @@
 	SealevelchangeRunCountEnum,
 	SealevelchangeTransitionsEnum,
-	SealevelchangeUElasticEnum,
+	SealevelchangeUViscoElasticEnum,
 	SettingsIoGatherEnum,
 	SettingsNumResultsOnNodesEnum,
@@ -392,4 +396,7 @@
 	SettingsSolverResidueThresholdEnum,
 	SettingsWaitonlockEnum,
+	StackNumStepsEnum,
+	StackTimesEnum,
+	StackIndexEnum,
 	SmbAIceEnum,
 	SmbAIdxEnum,
@@ -946,4 +953,8 @@
 	SolidearthExternalGeoidRateEnum,
 	SolidearthExternalBarystaticSeaLevelRateEnum,
+	StackRSLEnum,
+	StackUEnum,
+	StackNEnum,
+	StackEEnum,
 	StrainRateeffectiveEnum,
 	StrainRateparallelEnum,
Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26271)
+++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp	(revision 26272)
@@ -366,4 +366,5 @@
 		case RotationalAngularVelocityEnum : return "RotationalAngularVelocity";
 		case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
+		case SolidearthSettingsViscousEnum : return "SolidearthSettingsViscous";
 		case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone";
 		case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi";
@@ -375,10 +376,13 @@
 		case LoadLoveKEnum : return "LoadLoveK";
 		case LoadLoveLEnum : return "LoadLoveL";
+		case LoveTimeFreqEnum : return "LoveTimeFreq";
+		case LoveIsTimeEnum : return "LoveIsTime";
 		case SealevelchangeGRigidEnum : return "SealevelchangeGRigid";
-		case SealevelchangeGElasticEnum : return "SealevelchangeGElastic";
+		case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
 		case SolidearthSettingsComputesealevelchangeEnum : return "SolidearthSettingsComputesealevelchange";
 		case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
 		case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
-		case SealevelchangeHElasticEnum : return "SealevelchangeHElastic";
+		case SolidearthSettingsTimeAccEnum : return "SolidearthSettingsTimeAcc";
+		case SealevelchangeHViscoElasticEnum : return "SealevelchangeHViscoElastic";
 		case SolidearthSettingsHorizEnum : return "SolidearthSettingsHoriz";
 		case SolidearthSettingsMaxiterEnum : return "SolidearthSettingsMaxiter";
@@ -391,5 +395,5 @@
 		case SealevelchangeRunCountEnum : return "SealevelchangeRunCount";
 		case SealevelchangeTransitionsEnum : return "SealevelchangeTransitions";
-		case SealevelchangeUElasticEnum : return "SealevelchangeUElastic";
+		case SealevelchangeUViscoElasticEnum : return "SealevelchangeUViscoElastic";
 		case SettingsIoGatherEnum : return "SettingsIoGather";
 		case SettingsNumResultsOnNodesEnum : return "SettingsNumResultsOnNodes";
@@ -400,4 +404,7 @@
 		case SettingsSolverResidueThresholdEnum : return "SettingsSolverResidueThreshold";
 		case SettingsWaitonlockEnum : return "SettingsWaitonlock";
+		case StackNumStepsEnum : return "StackNumSteps";
+		case StackTimesEnum : return "StackTimes";
+		case StackIndexEnum : return "StackIndex";
 		case SmbAIceEnum : return "SmbAIce";
 		case SmbAIdxEnum : return "SmbAIdx";
@@ -951,4 +958,8 @@
 		case SolidearthExternalGeoidRateEnum : return "SolidearthExternalGeoidRate";
 		case SolidearthExternalBarystaticSeaLevelRateEnum : return "SolidearthExternalBarystaticSeaLevelRate";
+		case StackRSLEnum : return "StackRSL";
+		case StackUEnum : return "StackU";
+		case StackNEnum : return "StackN";
+		case StackEEnum : return "StackE";
 		case StrainRateeffectiveEnum : return "StrainRateeffective";
 		case StrainRateparallelEnum : return "StrainRateparallel";
Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26271)
+++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp	(revision 26272)
@@ -372,4 +372,5 @@
 	      else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum;
 	      else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
+	      else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
 	      else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
 	      else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
@@ -381,13 +382,16 @@
 	      else if (strcmp(name,"LoadLoveK")==0) return LoadLoveKEnum;
 	      else if (strcmp(name,"LoadLoveL")==0) return LoadLoveLEnum;
-	      else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
-	      else if (strcmp(name,"SealevelchangeGElastic")==0) return SealevelchangeGElasticEnum;
+	      else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum;
          else stage=4;
    }
    if(stage==4){
-	      if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
+	      if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum;
+	      else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
+	      else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
+	      else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
 	      else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
 	      else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
-	      else if (strcmp(name,"SealevelchangeHElastic")==0) return SealevelchangeHElasticEnum;
+	      else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
+	      else if (strcmp(name,"SealevelchangeHViscoElastic")==0) return SealevelchangeHViscoElasticEnum;
 	      else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
 	      else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
@@ -400,5 +404,5 @@
 	      else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
 	      else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
-	      else if (strcmp(name,"SealevelchangeUElastic")==0) return SealevelchangeUElasticEnum;
+	      else if (strcmp(name,"SealevelchangeUViscoElastic")==0) return SealevelchangeUViscoElasticEnum;
 	      else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
 	      else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum;
@@ -409,4 +413,7 @@
 	      else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum;
 	      else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
+	      else if (strcmp(name,"StackNumSteps")==0) return StackNumStepsEnum;
+	      else if (strcmp(name,"StackTimes")==0) return StackTimesEnum;
+	      else if (strcmp(name,"StackIndex")==0) return StackIndexEnum;
 	      else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
 	      else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
@@ -499,5 +506,8 @@
 	      else if (strcmp(name,"TimesteppingFinalTime")==0) return TimesteppingFinalTimeEnum;
 	      else if (strcmp(name,"TimesteppingInterpForcing")==0) return TimesteppingInterpForcingEnum;
-	      else if (strcmp(name,"TimesteppingCycleForcing")==0) return TimesteppingCycleForcingEnum;
+         else stage=5;
+   }
+   if(stage==5){
+	      if (strcmp(name,"TimesteppingCycleForcing")==0) return TimesteppingCycleForcingEnum;
 	      else if (strcmp(name,"TimesteppingStartTime")==0) return TimesteppingStartTimeEnum;
 	      else if (strcmp(name,"TimesteppingTimeStep")==0) return TimesteppingTimeStepEnum;
@@ -506,8 +516,5 @@
 	      else if (strcmp(name,"TimesteppingType")==0) return TimesteppingTypeEnum;
 	      else if (strcmp(name,"ToMITgcmComm")==0) return ToMITgcmCommEnum;
-         else stage=5;
-   }
-   if(stage==5){
-	      if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
+	      else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
 	      else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
 	      else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
@@ -622,5 +629,8 @@
 	      else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
 	      else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
-	      else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
+         else stage=6;
+   }
+   if(stage==6){
+	      if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
 	      else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
 	      else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
@@ -629,8 +639,5 @@
 	      else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum;
 	      else if (strcmp(name,"DistanceToGroundingline")==0) return DistanceToGroundinglineEnum;
-         else stage=6;
-   }
-   if(stage==6){
-	      if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
+	      else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
 	      else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
 	      else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
@@ -745,5 +752,8 @@
 	      else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
 	      else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
-	      else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
+         else stage=7;
+   }
+   if(stage==7){
+	      if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
 	      else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
 	      else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
@@ -752,8 +762,5 @@
 	      else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
 	      else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
-         else stage=7;
-   }
-   if(stage==7){
-	      if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
+	      else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
 	      else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
 	      else if (strcmp(name,"Node")==0) return NodeEnum;
@@ -868,5 +875,8 @@
 	      else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
 	      else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
-	      else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
+         else stage=8;
+   }
+   if(stage==8){
+	      if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
 	      else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
 	      else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
@@ -875,8 +885,5 @@
 	      else if (strcmp(name,"SmbC")==0) return SmbCEnum;
 	      else if (strcmp(name,"SmbCcsnowValue")==0) return SmbCcsnowValueEnum;
-         else stage=8;
-   }
-   if(stage==8){
-	      if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
+	      else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
 	      else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
 	      else if (strcmp(name,"SmbD")==0) return SmbDEnum;
@@ -972,4 +979,8 @@
 	      else if (strcmp(name,"SolidearthExternalGeoidRate")==0) return SolidearthExternalGeoidRateEnum;
 	      else if (strcmp(name,"SolidearthExternalBarystaticSeaLevelRate")==0) return SolidearthExternalBarystaticSeaLevelRateEnum;
+	      else if (strcmp(name,"StackRSL")==0) return StackRSLEnum;
+	      else if (strcmp(name,"StackU")==0) return StackUEnum;
+	      else if (strcmp(name,"StackN")==0) return StackNEnum;
+	      else if (strcmp(name,"StackE")==0) return StackEEnum;
 	      else if (strcmp(name,"StrainRateeffective")==0) return StrainRateeffectiveEnum;
 	      else if (strcmp(name,"StrainRateparallel")==0) return StrainRateparallelEnum;
@@ -987,5 +998,8 @@
 	      else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
 	      else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
-	      else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+         else stage=9;
+   }
+   if(stage==9){
+	      if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
 	      else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
 	      else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
@@ -998,8 +1012,5 @@
 	      else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
 	      else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
-         else stage=9;
-   }
-   if(stage==9){
-	      if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
+	      else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
 	      else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
 	      else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
@@ -1110,5 +1121,8 @@
 	      else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
 	      else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
-	      else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
+         else stage=10;
+   }
+   if(stage==10){
+	      if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum;
 	      else if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum;
 	      else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum;
@@ -1121,8 +1135,5 @@
 	      else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
 	      else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
-         else stage=10;
-   }
-   if(stage==10){
-	      if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
+	      else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
 	      else if (strcmp(name,"Outputdefinition75")==0) return Outputdefinition75Enum;
 	      else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
@@ -1233,5 +1244,8 @@
 	      else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
 	      else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum;
-	      else if (strcmp(name,"Element")==0) return ElementEnum;
+         else stage=11;
+   }
+   if(stage==11){
+	      if (strcmp(name,"Element")==0) return ElementEnum;
 	      else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
 	      else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
@@ -1244,8 +1258,5 @@
 	      else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
 	      else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
-         else stage=11;
-   }
-   if(stage==11){
-	      if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
+	      else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
 	      else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
 	      else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
@@ -1356,5 +1367,8 @@
 	      else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
 	      else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
-	      else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
+         else stage=12;
+   }
+   if(stage==12){
+	      if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
 	      else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
 	      else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
@@ -1367,8 +1381,5 @@
 	      else if (strcmp(name,"Melange")==0) return MelangeEnum;
 	      else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
-         else stage=12;
-   }
-   if(stage==12){
-	      if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
+	      else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
 	      else if (strcmp(name,"MeshX")==0) return MeshXEnum;
 	      else if (strcmp(name,"MeshY")==0) return MeshYEnum;
@@ -1479,5 +1490,8 @@
 	      else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
 	      else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
-	      else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
+         else stage=13;
+   }
+   if(stage==13){
+	      if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
 	      else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
 	      else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
@@ -1490,8 +1504,5 @@
 	      else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
 	      else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
-         else stage=13;
-   }
-   if(stage==13){
-	      if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
+	      else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
 	      else if (strcmp(name,"Tetra")==0) return TetraEnum;
 	      else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;
Index: /issm/trunk-jpl/src/m/classes/lovenumbers.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/lovenumbers.m	(revision 26271)
+++ /issm/trunk-jpl/src/m/classes/lovenumbers.m	(revision 26272)
@@ -9,5 +9,5 @@
 classdef lovenumbers
 	properties (SetAccess=public) 
-		
+
 		%regular love numbers:
 		h           = []; %provided by PREM model
@@ -20,4 +20,8 @@
 		tl          = []; 
 		tk2secular  = 0;  %deg 2 secular number.
+
+		%time/frequency for visco-elastic love numbers
+		timefreq    = [];
+		istime      = 1; 
 
 	end
@@ -41,5 +45,8 @@
 			%secular fluid love number: 
 			self.tk2secular=0.942; 
-
+			
+			%time: 
+			self.istime=1; %temporal love numbers by default
+			self.timefreq=0; %elastic case by default.
 
 		end % }}}
@@ -59,4 +66,6 @@
 			md = checkfield(md,'fieldname','solidearth.lovenumbers.tl','NaN',1,'Inf',1);
 			md = checkfield(md,'fieldname','solidearth.lovenumbers.tk2secular','NaN',1,'Inf',1);
+			md = checkfield(md,'fieldname','solidearth.lovenumbers.timefreq','NaN',1,'Inf',1);
+			md = checkfield(md,'fieldname','solidearth.lovenumbers.istime','NaN',1,'Inf',1,'values',[0 1]);
 
 			%check that love numbers are provided at the same level of accuracy: 
@@ -65,4 +74,8 @@
 			end
 			
+			ntf=length(self.timefreq);
+			if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
+				error('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector');
+			end
 
 		end % }}}
@@ -81,4 +94,7 @@
 			fielddisplay(self,'tl','tidal load Love number (deg 2)');
 			fielddisplay(self,'tk2secular','secular fluid Love number');
+			
+			fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)');
+			fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
 
 		end % }}}
@@ -94,4 +110,12 @@
 			WriteData(fid,prefix,'object',self,'data',self.tk2secular,'fieldname','lovenumbers.tk2secular','format','Double');
 
+			if self.istime,
+				scale=md.constants.yts; 
+			else
+				scale=1.0/md.constants.yts;
+			end
+			WriteData(fid,prefix,'object',self,'fieldname','istime','name','md.solidearth.lovenumbers.istime','format','Boolean');
+			WriteData(fid,prefix,'object',self,'fieldname','timefreq','name','md.solidearth.lovenumbers.timefreq','format','DoubleMat','mattype',1,'scale',scale);
+
 		end % }}}
 		function savemodeljs(self,fid,modelname) % {{{
@@ -99,4 +123,6 @@
 			writejs1Darray(fid,[modelname '.lovenumbers.k'],self.k);
 			writejs1Darray(fid,[modelname '.lovenumbers.l'],self.l);
+			writejs1Darray(fid,[modelname '.lovenumbers.istime'],self.istime);
+			writejs1Darray(fid,[modelname '.lovenumbers.time'],self.time);
 		end % }}}
 		function self = extrude(self,md) % {{{
Index: /issm/trunk-jpl/src/m/classes/solidearthsettings.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 26271)
+++ /issm/trunk-jpl/src/m/classes/solidearthsettings.m	(revision 26272)
@@ -11,4 +11,5 @@
 		rigid                  = 0;
 		elastic                = 0;
+		viscous                = 0;
 		rotation               = 0;
 		ocean_area_scaling     = 0;
@@ -17,5 +18,6 @@
 		isgrd                  = 0; %will GRD patterns be computed? 
 		compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed? 
-		degacc                 = 0; %degree increment for resolution of Green tables
+		degacc                 = 0; %degree increment for resolution of Green tables.
+		timeacc                = 0; %time step accurary required to compute Green tables
 		horiz                  = 0; %compute horizontal deformation
 		grdmodel               = 0; %grd model (0 by default, 1 for elastic, 2 for Ivins)
@@ -43,4 +45,5 @@
 		self.rigid=1;
 		self.elastic=1;
+		self.viscous=1;
 		self.rotation=1;
 		self.ocean_area_scaling=0;
@@ -51,4 +54,5 @@
 		%numerical discretization accuracy
 		self.degacc=.01;
+		self.timeacc=1; 
 
 		%how many time steps we skip before we run solidearthsettings solver during transient
@@ -75,4 +79,5 @@
 			md = checkfield(md,'fieldname','solidearth.settings.runfrequency','size',[1 1],'>=',1);
 			md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10);
+			md = checkfield(md,'fieldname','solidearth.settings.timeacc','size',[1 1],'>',0);
 			md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]);
 			md = checkfield(md,'fieldname','solidearth.settings.grdmodel','>=',0,'<=',2);
@@ -82,4 +87,7 @@
 			if self.elastic==1 & self.rigid==0,
 				error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
+			end
+			if self.viscous==1 & self.elastic==0,
+				error('solidearthsettings checkconsistency error message: need elastic on if viscous flag is set');
 			end
 
@@ -101,5 +109,4 @@
 			end
 
-
 		end % }}}
 		function disp(self) % {{{
@@ -116,6 +123,8 @@
 			fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
 			fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
+			fielddisplay(self,'viscous','viscous earth graviational potential perturbation');
 			fielddisplay(self,'rotation','earth rotational potential perturbation');
 			fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
+			fielddisplay(self,'timeacc','time accuracy (default 1 yr) for numerical discretization of the Green''s functions');
 			fielddisplay(self,'grdmodel','type of deformation model, 1 for elastic, 2 for visco-elastic from Ivins');
 			fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
@@ -127,8 +136,10 @@
 			WriteData(fid,prefix,'object',self,'fieldname','rigid','name','md.solidearth.settings.rigid','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','elastic','name','md.solidearth.settings.elastic','format','Boolean');
+			WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean');
 			WriteData(fid,prefix,'object',self,'fieldname','runfrequency','name','md.solidearth.settings.runfrequency','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','degacc','name','md.solidearth.settings.degacc','format','Double');
+			WriteData(fid,prefix,'object',self,'fieldname','timeacc','name','md.solidearth.settings.timeacc','format','Double','scale',md.constants.yts);
 			WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer');
 			WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer');
@@ -145,8 +156,10 @@
 			writejsdouble(fid,[modelname '.solidearth.settings.rigid'],self.rigid);
 			writejsdouble(fid,[modelname '.solidearth.settings.elastic'],self.elastic);
+			writejsdouble(fid,[modelname '.solidearth.settings.viscous'],self.viscous);
 			writejsdouble(fid,[modelname '.solidearth.settings.rotation'],self.rotation);
 			writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling);
 			writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);
 			writejsdouble(fid,[modelname '.solidearth.settings.degacc'],self.degacc);
+			writejsdouble(fid,[modelname '.solidearth.settings.timeacc'],self.timeacc);
 			writejsdouble(fid,[modelname '.solidearth.settings.cross_section_shape'],self.cross_section_shape);
 		end % }}}
Index: /issm/trunk-jpl/test/NightlyRun/test2090.m
===================================================================
--- /issm/trunk-jpl/test/NightlyRun/test2090.m	(revision 26272)
+++ /issm/trunk-jpl/test/NightlyRun/test2090.m	(revision 26272)
@@ -0,0 +1,114 @@
+%Test Name: ViscoElasticEarthSlc
+
+%mesh earth:
+md=model;
+md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',500.); %700 km resolution mesh
+
+%Geometry for the bed, arbitrary thickness of 1000: 
+md.geometry.bed=-ones(md.mesh.numberofvertices,1);
+md.geometry.base=md.geometry.bed;
+md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1);
+md.geometry.surface=md.geometry.bed+md.geometry.thickness;
+
+
+%parameterize solidearth solution:
+%solidearth loading:  {{{
+md.masstransport.spcthickness=[md.geometry.thickness;0];
+md.dsl.global_average_thermosteric_sea_level=[0;0];
+md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1);
+md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1);
+md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
+%antarctica
+xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
+ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
+ze=md.mesh.z(md.mesh.elements)*[1;1;1]/3;
+re=sqrt(xe.^2+ye.^2+ze.^2);
+
+late=asind(ze./re);
+longe=atan2d(ye,xe);
+pos=find(late < -80);
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posant=pos;
+%greenland
+pos=find(late>60 & late<90 & longe>-75 & longe<-15);
+md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
+posgre=pos;
+
+%visco-elastic loading from love numbers that we load ourselves. We 
+%still use the lovenumbers constructor to initialize the fields: 
+load ../Data/lnb_temporal.mat
+md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
+md.solidearth.lovenumbers.h=[zeros(1,100); h1];
+md.solidearth.lovenumbers.k=[zeros(1,100); k1];
+md.solidearth.lovenumbers.l=[zeros(1,100); l1];
+md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,100);
+md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,100);
+md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,100);
+md.solidearth.lovenumbers.timefreq=time;
+
+%}}}
+%mask:  {{{
+mask=gmtmask(md.mesh.lat,md.mesh.long);
+oceanmask=-ones(md.mesh.numberofvertices,1);
+pos=find(mask==0); oceanmask(pos)=1;
+
+icemask=ones(md.mesh.numberofvertices,1);
+icemask(md.mesh.elements(posant))=-1;
+icemask(md.mesh.elements(posgre))=-1;
+
+md.mask.ice_levelset=icemask;
+md.mask.ocean_levelset=oceanmask;
+% }}}
+
+%time stepping: 
+md.timestepping.start_time=0;
+md.timestepping.time_step=1;
+md.timestepping.final_time=10;
+
+md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1);
+md.initialization.vx=zeros(md.mesh.numberofvertices,1);
+md.initialization.vy=zeros(md.mesh.numberofvertices,1);
+md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
+md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1);
+md.initialization.dsl=zeros(md.mesh.numberofvertices,1);
+md.initialization.str=0;
+
+%Materials: 
+md.materials=materials('hydro');
+
+%Miscellaneous
+md.miscellaneous.name='test2090';
+
+%Solution parameters
+md.cluster.np=3;
+md.solidearth.settings.reltol=NaN;
+md.solidearth.settings.abstol=1e-3;
+md.solidearth.settings.computesealevelchange=1;
+md.solidearth.settings.isgrd=1;
+md.solidearth.settings.ocean_area_scaling=0;
+md.solidearth.settings.grdmodel=1;
+
+%Physics: 
+md.transient.issmb=0; 
+md.transient.isstressbalance=0;
+md.transient.isthermal=0;
+md.transient.ismasstransport=1;
+md.transient.isslc=1;
+md.solidearth.requested_outputs={'Sealevel','Bed','BedGRD'};
+
+
+% max number of iteration reverted back to 10 (i.e., the original default value)
+md.solidearth.settings.maxiter=10;
+
+%eustatic + rigid + elastic + rotation run:
+md.solidearth.settings.rigid=1;
+md.solidearth.settings.elastic=1;
+md.solidearth.settings.rotation=0;
+md=solve(md,'tr');
+S=md.results.TransientSolution.Sealevel;
+
+%Fields and tolerances to track changes
+field_names={'Sealevel'};
+field_tolerances={1e-13};
+field_values={S};