Changeset 26242

Timestamp:

05/05/21 15:40:27 (4 years ago)

Author:

Eric.Larour

Message:

CHG: implemented templated version of the love_core routine, to allow
for IssmComplex and (IssmQuad tuypes in the future). Created new IssmComplex
type.

Location:

issm/trunk-jpl

Files:

• src/c/analyses/LoveAnalysis.cpp (modified) (1 diff)
• src/c/analyses/SealevelchangeAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Tria.cpp (modified) (3 diffs)
• src/c/classes/ExternalResults/GenericExternalResult.h (modified) (1 diff)
• src/c/cores/love_core.cpp (modified) (37 diffs)
• src/c/shared/Enum/Enum.vim (modified) (3 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (22 diffs)
• src/c/shared/Numerics/types.h (modified) (2 diffs)
• src/m/classes/fourierlove.m (modified) (4 diffs)
• src/m/solve/loadresultfromdisk.m (modified) (1 diff)
• test/NightlyRun/test2084.m (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/LoveAnalysis.cpp

@@ -42 +42 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.love.inner_core_boundary",LoveInnerCoreBoundaryEnum));
         parameters->AddObject(iomodel->CopyConstantObject("md.love.core_mantle_boundary",LoveCoreMantleBoundaryEnum));
+        parameters->AddObject(iomodel->CopyConstantObject("md.love.complex_computation",LoveComplexComputationEnum));
 }/*}}}*/
 void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp

@@ -66 +66 @@
 
         /*Initialize solid earth motion and sea level: */
-        InputUpdateFromConstantx(inputs,elements,0.,BedEastEnum);
-        InputUpdateFromConstantx(inputs,elements,0.,BedNorthEnum);
+        iomodel->ConstantToInput(inputs,elements,0.,BedEastEnum,P1Enum);
+        iomodel->ConstantToInput(inputs,elements,0.,BedNorthEnum,P1Enum);
     iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum);
 

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -7169 +7169 @@
 
         IssmDouble sealevel[3]={0};
+        int i,e,l,nbar;
 
         if(loads->sealevelloads){
-                for(int i=0;i<NUMVERTICES;i++) {
+                for(i=0;i<NUMVERTICES;i++) {
                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                        for (int e=0;e<nel;e++){
+                        for (e=0;e<nel;e++){
                                 sealevel[i]+=G[i*nel+e]*(loads->sealevelloads[e]+loads->loads[e]);
                         }
-                        for(int l=0;l<SLGEOM_NUMLOADS;l++){
-                                int nbar=slgeom->nbar[l];
-                                for (int e=0;e<nbar;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]);
                                 }
                                 if(l==SLGEOM_OCEAN){
-                                        for (int e=0;e<nbar;e++){
+                                        for (e=0;e<nbar;e++){
                                                 sealevel[i]+=Gsub[l][i*nbar+e]*(loads->subsealevelloads[e]);
                                         }
@@ -7190 +7191 @@
         }
         else{  //this is the initial convolution where only loads are provided
-                for(int i=0;i<NUMVERTICES;i++) {
+                for(i=0;i<NUMVERTICES;i++) {
                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
-                        for (int e=0;e<nel;e++){
+                        for (e=0;e<nel;e++){
                                 sealevel[i]+=G[i*nel+e]*(loads->loads[e]);
                         }
-                        for(int l=0;l<SLGEOM_NUMLOADS;l++){
-                                int nbar=slgeom->nbar[l];
-                                for (int e=0;e<nbar;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]);
                                 }
@@ -7206 +7207 @@
         /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
         if(percpu){
-                for(int i=0;i<NUMVERTICES;i++){
+                for(i=0;i<NUMVERTICES;i++){
                         if(slgeom->lids[this->vertices[i]->lid]==this->lid)sealevelout[this->vertices[i]->lid]=sealevel[i];
                 }
         }
         else{
-                for(int i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i];
+                for(i=0;i<NUMVERTICES;i++) sealevelout[i]=sealevel[i];
         }
 

issm/trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h

@@ -679 +679 @@
 #endif
 
+        /*Specific instantiations for IssmComplex*: */
+template <> inline GenericExternalResult<IssmComplex*>::GenericExternalResult(int in_id, int in_enum_type,IssmComplex* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/
+
+        id = in_id;
+        M  = in_M;
+        N  = in_N;
+
+        EnumToStringx(&this->result_name,in_enum_type);
+
+        step = in_step;
+        time = in_time;
+
+        /*Copy result in values*/
+        if(M*N){
+                value=xNew<IssmComplex>(M*N);
+                xMemCpy<IssmComplex>(value,in_values,M*N);
+        }
+        else value=NULL;
+}
+/*}}}*/
+template <> inline void GenericExternalResult<IssmComplex*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/
+
+        int     my_rank;
+        int     type;
+        int     rows,cols;
+        char   *name    = NULL;
+        IssmPDouble passiveDouble;
+        IssmDouble* reals=NULL;
+        IssmDouble* imags=NULL;
+
+        /*recover my_rank:*/
+        my_rank=IssmComm::GetRank();
+
+        if(io_gather){
+                /*we are gathering the data on cpu 0, don't write on other cpus: */
+                if(my_rank) return;
+        }
+
+        /*First write enum: */
+        int length=(strlen(this->result_name)+1)*sizeof(char);
+        fwrite(&length,sizeof(int),1,fid);
+        fwrite(this->result_name,length,1,fid);
+
+        /*Now write time and step: */
+        passiveDouble=reCast<IssmPDouble>(time);
+        fwrite(&passiveDouble,sizeof(IssmPDouble),1,fid);
+        fwrite(&step,sizeof(int),1,fid);
+
+        /*writing a IssmComplex array, type is 3:*/
+        type=5;
+        fwrite(&type,sizeof(int),1,fid);
+        rows=this->M;
+        fwrite(&rows,sizeof(int),1,fid);
+        cols=this->N;
+        fwrite(&cols,sizeof(int),1,fid);
+
+        /*write complex array into two real arrays:*/
+        reals=xNew<IssmDouble>(cols*rows);
+        imags=xNew<IssmDouble>(cols*rows);
+        for(int i=0;i<cols*rows;i++){
+                reals[i]=value[i].real();
+                imags[i]=value[i].imag();
+        }
+        fwrite(reals,cols*rows*sizeof(IssmComplex),1,fid);
+        fwrite(imags,cols*rows*sizeof(IssmComplex),1,fid);
+
+}
+/*}}}*/
+
         /*Specifics instantiations for Vector*/
         template <> inline GenericExternalResult<Vector<IssmPDouble>*>::GenericExternalResult(int in_id, int in_enum_type,Vector<IssmPDouble>* in_values, int in_M,int in_N,int in_step,IssmDouble in_time){/*{{{*/

issm/trunk-jpl/src/c/cores/love_core.cpp

@@ -9 +9 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
-#include <petscblaslapack.h>
-
-void love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/
-        // initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system
-        bool        verbosemod = (int)VerboseModule();
-
-        int numlayers = matlitho->numlayers;
-        IssmDouble* r=matlitho->radius;
-        IssmDouble* EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
-        IssmDouble r1,r2,ro, GG;
-        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
-
-        for (int i=0;i<numlayers;i++){
-                  r2 = r[i+1];
-                  ro = matlitho->density[i];
-                if (i==0){
-                        EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0;
-                }else{
-                        r1=r[i];
-                        EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;;
-                }
-        }
-
-        IssmDouble g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0);
-        femmodel->parameters->SetParam(g0,LoveG0Enum);
-        femmodel->parameters->SetParam(EarthMass,numlayers,LoveEarthMassEnum);
-        xDelete<IssmDouble>(EarthMass);
-
-        int nyi=6*(numlayers+1);
-        int starting_layer=0;
-        femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
-        femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
-}/*}}}*/
-void GetEarthRheology(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+#include "petscblaslapack.h"
+
+/*local definitions:*/
+class LoveVariables{  /*{{{*/
+
+        public: 
+                IssmDouble g0; 
+                IssmDouble r0;
+                IssmDouble* EarthMass; 
+                int nyi; 
+                int starting_layer;
+                
+                LoveVariables(){  /*{{{*/
+                        g0=0;
+                        r0=0;
+                        EarthMass=NULL;
+                        nyi=0;
+                        starting_layer=0;
+                } /*}}}*/
+                LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin){  /*{{{*/
+                        EarthMass=EarthMassin;
+                        g0=g0in;
+                        r0=r0in;
+                        nyi=nyiin;
+                        starting_layer=starting_layerin;
+                } /*}}}*/
+                ~LoveVariables(){};
+}; /*}}}*/
+
+/*self contained support routines used by cores below:*/
+template<typename doubletype> doubletype                 angular_frequency(IssmDouble frequency); //pure declaration
+template <> IssmDouble                     angular_frequency<IssmDouble>(IssmDouble frequency){ /*{{{*/
+        return 2.0*PI*frequency;
+} /*}}}*/
+template <> IssmComplex                    angular_frequency<IssmComplex>(IssmDouble frequency){ /*{{{*/
+        IssmComplex value=reCast<IssmComplex>(complex<double>(0,1))*2.0*PI*reCast<IssmComplex>(frequency);
+        return value;
+} /*}}}*/
+template<typename doubletype> void                       allgesv(int* pnyi, int* pnrhs, doubletype* yilocal, int* plda, int* ipiv, doubletype* rhslocal, int* pldb, int* pinfo);
+template <> void                           allgesv<IssmDouble>(int* pnyi, int* pnrhs, IssmDouble* yilocal, int* plda, int* ipiv, IssmDouble* rhslocal, int* pldb, int* pinfo){ /*{{{*/
+        dgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo);
+} /*}}}*/
+template <> void                           allgesv<IssmComplex>(int* pnyi, int* pnrhs, IssmComplex* yilocal, int* plda, int* ipiv, IssmComplex* rhslocal, int* pldb, int* pinfo){ /*{{{*/
+        _error_("zgesv_ not linked correctly yet! ");
+        //zgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo);
+} /*}}}*/
+template<typename doubletype> doubletype   factorial(int n){ /*{{{*/
+        doubletype prod=1;
+        for (int i=2;i<n+1;i++) prod*=i;
+        return prod;
+}/*}}}*/
+template<typename doubletype> doubletype   n_C_r(int n, int r){ /*{{{*/ 
+        //n choose r
+        int primes[169] = 
+        {2,    3,    5,    7,   11,   13,   17,   19,   23,   29,  
+                31,   37,   41,   43,   47,   53,   59,   61,   67,   71,  
+                73,   79,   83,   89,   97,  101,  103,  107,  109,  113,  
+                127,  131,  137,  139,  149,  151,  157,  163,  167,  173,  
+                179,  181,  191,  193,  197,  199,  211,  223,  227,  229,  
+                233,  239,  241,  251,  257,  263,  269,  271,  277,  281,  
+                283,  293,  307,  311,  313,  317,  331,  337,  347,  349,  
+                353,  359,  367,  373,  379,  383,  389,  397,  401,  409,  
+                419,  421,  431,  433,  439,  443,  449,  457,  461,  463,  
+                467,  479,  487,  491,  499,  503,  509,  521,  523,  541,  
+                547,  557,  563,  569,  571,  577,  587,  593,  599,  601,  
+                607,  613,  617,  619,  631,  641,  643,  647,  653,  659,  
+                661,  673,  677,  683,  691,  701,  709,  719,  727,  733,  
+                739,  743,  751,  757,  761,  769,  773,  787,  797,  809,  
+                811,  821,  823,  827,  829,  839,  853,  857,  859,  863,  
+                877,  881,  883,  887,  907,  911,  919,  929,  937,  941,  
+                947,  953,  967,  971,  977,  983,  991,  997, 1009};
+        int num, den;
+        num = 1;
+        den = 1;
+
+        for (int i=0;i<r;i++){
+                num = num*(n-i);
+                den = den*(i+1);
+                if (i>0) {
+                        // Divide out common prime factors
+                        for (int k=0;k<169;k++){ //169 is the length of the prime vector here
+                                if ( i % primes[k] == 0) { // modulo
+                                        num = num/primes[k];
+                                        den = den/primes[k];
+                                }
+                        }
+                }
+        }
+
+        doubletype res;        
+        return res = num/den;
+}/*}}}*/
+template<typename doubletype> doubletype*  postwidder_coef(int NTit){ /*{{{*/
+        //Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform:
+        //The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k))
+        //The method is based on equations (2), (6), (7) in: 
+        //Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004)
+        //Note that the coefficients xi lack the factor s=k*log(2)/t. 
+        //That is because we are computing the heaviside response of the system rather than its impulse response, 
+        //and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M.
+
+        doubletype* xi=xNewZeroInit<doubletype>(2*NTit*NTit);
+        int indxi;
+        for (int M=1;M<NTit+1;M++){
+                for (int k=1;k<2*M+1;k++){
+                        indxi=(M-1)*(2*NTit)+k-1;
+                        for (int j=floor((k+1)/2);j<min(k,M)+1;j++){
+                                xi[indxi]+=pow(j,M+1.0)/factorial<doubletype>(M)*n_C_r<doubletype>(M,j)*n_C_r<doubletype>(2*j,j)*n_C_r<doubletype>(j,k-j);
+                        }
+                        xi[indxi]*=pow(-1.0,k+M)/k;
+                }
+        }
+        return xi;
+}/*}}}*/
+template<typename doubletype> void         postwidder_transform(doubletype* Lovet, doubletype* Lovef,int d, int t, int NTit, doubletype* xi, FemModel* femmodel){ /*{{{*/
+        //Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef
+
+        int nfreq, indxi, indf;
+        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
+        int nt=nfreq/2/NTit;
+
+        indf=d*nfreq+t*2*NTit;
+        doubletype LoveM[NTit];
+
+        for (int M=1;M<NTit+1;M++){
+                LoveM[M-1]=0.0;
+                for (int k=1;k<2*M+1;k++){
+                        indxi=(M-1)*(2*NTit)+k-1;
+                        LoveM[M-1]+=xi[indxi]*Lovef[indf+k-1];
+                }
+
+                //Make sure we are not getting into numerical instability
+                //Diverging once: ok, we'll give that the benefit of the doubt, it can be an inflexion point in the convergence series
+                //Diverging twice in a row: we are definitely propagating numerical error: revert to the last stable value and exit
+                if (M>3){ 
+                        if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) &&
+                                        abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){
+                                Lovet[d*nt+t]=LoveM[M-3];
+                                return;
+                        }
+                }
+        }
+        Lovet[d*nt+t]=LoveM[NTit-1];
+}/*}}}*/
+template <typename doubletype> void        GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega,  Matlitho* matlitho){ /*{{{*/
 
         //returns lame parameters (material rigity) lambda and mu for the right frequency and layer
-        IssmDouble vi=matlitho->viscosity[layer_index];
-        IssmDouble mu=matlitho->lame_mu[layer_index];
-        IssmDouble la=matlitho->lame_lambda[layer_index];
-        int rheo=matlitho->rheologymodel[layer_index];
+        doubletype mu,la;
+
+        IssmDouble vi=matlitho->viscosity[layer_index];
+        IssmDouble mu0=matlitho->lame_mu[layer_index];
+        IssmDouble la0=matlitho->lame_lambda[layer_index];
+        int rheo=matlitho->rheologymodel[layer_index];
 
         if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
-                                   //Otherwise return the elastic value
-                IssmDouble ka=la + 2.0/3.0*mu; //Bulk modulus
+                //Otherwise return the elastic value
+                IssmDouble ka=la0 + 2.0/3.0*mu0; //Bulk modulus
                 if (rheo==2){//EBM
                         IssmDouble alpha=matlitho->ebm_alpha[layer_index];
@@ -58 +171 @@
                         IssmDouble taul=matlitho->ebm_taul[layer_index];
                         IssmDouble tauh=matlitho->ebm_tauh[layer_index];
-                        IssmDouble hf1,hf2,U1,U2;
+                        IssmDouble hf1,hf2;
+                        doubletype U1,U2;
                         IssmDouble* a=xNew<IssmDouble>(2);
                         IssmDouble* b=xNew<IssmDouble>(1);
@@ -79 +193 @@
 
                         mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega);
-                        la=ka-2.0/3.0*mu;
+                        la=ka-2.0/3.0*mu;
 
                         xDelete<IssmDouble>(a);
@@ -88 +202 @@
                         IssmDouble mu2=matlitho->burgers_mu[layer_index];
 
-                        mu=mu*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu/vi)+mu/vi2*omega);
-                        la=ka-2.0/3.0*mu;
-                }
+                        mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega);
+                        la=ka-2.0/3.0*mu0;
+                }
                 else{//Maxwell
-                        la = (la + mu*ka/vi/omega)/(1.0 + mu/vi/omega);
-                        mu = mu/(1.0+mu/vi/omega);
-                }
-         }
-
-         *pla=la;
-         *pmu=mu;
+                        la = (la0 + mu0*ka/vi/omega)/(1.0 + mu0/vi/omega);
+                        mu = mu0/(1.0+mu0/vi/omega);
+                }
+        }
+
+        *pla=la;
+        *pmu=mu;
 
 } /*}}}*/
-IssmDouble GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+IssmDouble                                 GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/
         //computes gravity at radius r2
         IssmDouble* EarthMass;
         IssmDouble g, GG;
 
-        femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
+        EarthMass=vars->EarthMass;
         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
         IssmDouble ro=matlitho->density[layer_index];
@@ -120 +234 @@
         return  g= GG*M/pow(r2,2.0);
 }/*}}}*/
-void fill_yi_prefactor(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+template <typename doubletype> void        fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
         //precalculates partial derivative factors for function yi_derivatives
- IssmDouble ra=matlitho->radius[matlitho->numlayers];
- IssmDouble  g0,r0,mu0, GG;
- int nstep,nindex, starting_layer;
- femmodel->parameters->FindParam(&g0,LoveG0Enum);
- femmodel->parameters->FindParam(&r0,LoveR0Enum);
- femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
- femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
- femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
- femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
-
- IssmDouble frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
- IssmDouble xmin,xmax,x,dr;
- IssmDouble g,ro, issolid;
+        IssmDouble ra=matlitho->radius[matlitho->numlayers];
+        IssmDouble  g0,r0,mu0, GG;
+        int nstep,nindex, starting_layer;
+
+        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
+        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
+        
+        g0=vars->g0;
+        r0=vars->r0;
+        starting_layer=vars->starting_layer;
+
+        doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
+        IssmDouble xmin,xmax,x,dr;
+        IssmDouble g,ro, issolid;
 
         if (pomega) { //frequency and degree dependent terms /*{{{*/
-                IssmDouble la,mu;
-                IssmDouble* f=xNewZeroInit<IssmDouble>(12);
+                doubletype la,mu;
+                doubletype f[12];
                 int deg=*pdeg;
-                IssmDouble omega=*pomega;       
+                doubletype omega=*pomega;       
                 fn=deg*(deg+1.0);
 
                 for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
 
-                                ro=matlitho->density[layer_index];
-                                issolid=matlitho->issolid[layer_index];
-                                if(issolid==1){
-                                GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
-        
+                        ro=matlitho->density[layer_index];
+                        issolid=matlitho->issolid[layer_index];
+                        if(issolid==1){
+                                GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
+
                                 /*_______Expressions*/
                                 flm=(la+2.0*mu);
@@ -172 +288 @@
                                 xmax=(matlitho->radius[layer_index+1])/ra;
                                 dr = (xmax -xmin)/nstep;
-                                x=xmin;
+                                x=xmin;
                                 for (int n=0;n<nstep;n++){
-                                        g=GetGravity(x*ra,layer_index,femmodel,matlitho);
+                                        g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
 
                                         nindex=layer_index*nstep*36+n*36;
@@ -191 +307 @@
                         }
                 }
-
-                xDelete<IssmDouble>(f);
-        /*}}}*/
+                /*}}}*/
         } else if (pdeg) { // degree dependent terms /*{{{*/
-                int deg=*pdeg;
-                fn=(deg*(deg+1.0));
+                int deg=*pdeg;
+                fn=(deg*(deg+1.0));
 
                 for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
                         ro=matlitho->density[layer_index];
                         issolid=matlitho->issolid[layer_index];
-        
+
                         /*_______Expressions*/
                         fgr=4.0*PI*GG*ro*ra;
-         
+
                         xmin=matlitho->radius[layer_index]/ra;
                         xmax=(matlitho->radius[layer_index+1])/ra;
                         dr = (xmax -xmin)/nstep;
-                        x=xmin;
+                        x=xmin;
                         for (int n=0;n<nstep;n++){
-                                 nindex=layer_index*nstep*36+n*36;
-                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
-        
-                                 if(issolid==1){
-                                 yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
-                                 yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
-                                 yi_prefactor[nindex+ 5*6+4]= fn/(x*x);              // in dy[5*6+4]
-                                 } else {
-                                 yi_prefactor[nindex+ 1*6+0]= (-4.0*(fgr/g)+fn/x)/x; // in dy[1*6+0] liquid layer
-                                 }
+                                nindex=layer_index*nstep*36+n*36;
+                                g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
+
+                                if(issolid==1){
+                                        yi_prefactor[nindex+ 1*6+3]= fn/x;                  // in dy[1*6+3]
+                                        yi_prefactor[nindex+ 5*6+2]= -(fgr/g0*fn)/x;        // in dy[5*6+2]
+                                        yi_prefactor[nindex+ 5*6+4]= fn/(x*x);               // in dy[5*6+4]
+                                } else {
+                                        yi_prefactor[nindex+ 1*6+0]= (-4.0*(fgr/g)+fn/x)/x; // in dy[1*6+0] liquid layer
+                                }
                                 x=x+dr;
                         }
                 }
-        /*}}}*/
+                /*}}}*/
         } else { // static terms /*{{{*/
                 for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
                         ro=matlitho->density[layer_index];
                         issolid=matlitho->issolid[layer_index];
-        
+
                         /*_______Expressions*/
                         frhg0=ro*g0*ra/mu0;
@@ -236 +350 @@
                         xmax=(matlitho->radius[layer_index+1])/ra;
                         dr = (xmax -xmin)/nstep;
-                        x=xmin;
-                        for (int n=0;n<nstep;n++){
-                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho);
-                                 nindex=layer_index*nstep*36+n*36;
-                                 if(issolid==1){
-                                         yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
-                                         yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
-                                         yi_prefactor[nindex+ 2*6+2]= 1.0/x;        // in dy[2*6+2]
-                                         yi_prefactor[nindex+ 3*6+3]= -3.0/x;       // in dy[3*6+3]
-                                         yi_prefactor[nindex+ 3*6+4]= -frhg0/x;     // in dy[3*6+4]
-                                         yi_prefactor[nindex+ 4*6+0]= fgr/g0;       // in dy[4*6+0]
-                                         yi_prefactor[nindex+ 4*6+5]= 1.0;          // in dy[4*6+5]
-                                         yi_prefactor[nindex+ 5*6+5]= -2.0/x;       // in dy[5*6+5]
-                                 } else {
-                                         yi_prefactor[nindex+ 0*6+0]= fgr/g;        // in dy[0*6+0] liquid layer
-                                         yi_prefactor[nindex+ 0*6+1]= 1.0;          // in dy[0*6+1] liquid layer
-                                         yi_prefactor[nindex+ 1*6+1]= -2.0/x-fgr/g; // in dy[1*6+1] liquid layer
-                                 }
+                        x=xmin;
+                        for (int n=0;n<nstep;n++){
+                                g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
+                                nindex=layer_index*nstep*36+n*36;
+                                if(issolid==1){
+                                        yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
+                                        yi_prefactor[nindex+ 2*6+0]= -1.0/x;       // in dy[2*6+0]
+                                        yi_prefactor[nindex+ 2*6+2]= 1.0/x;        // in dy[2*6+2]
+                                        yi_prefactor[nindex+ 3*6+3]= -3.0/x;       // in dy[3*6+3]
+                                        yi_prefactor[nindex+ 3*6+4]= -frhg0/x;     // in dy[3*6+4]
+                                        yi_prefactor[nindex+ 4*6+0]= fgr/g0;       // in dy[4*6+0]
+                                        yi_prefactor[nindex+ 4*6+5]= 1.0;          // in dy[4*6+5]
+                                        yi_prefactor[nindex+ 5*6+5]= -2.0/x;       // in dy[5*6+5]
+                                } else {
+                                        yi_prefactor[nindex+ 0*6+0]= fgr/g;        // in dy[0*6+0] liquid layer
+                                        yi_prefactor[nindex+ 0*6+1]= 1.0;          // in dy[0*6+1] liquid layer
+                                        yi_prefactor[nindex+ 1*6+1]= -2.0/x-fgr/g; // in dy[1*6+1] liquid layer
+                                }
                                 x=x+dr;
                         }
                 }
-        /*}}}*/
-        }
-}/*}}}*/
-IssmDouble* yi_derivatives(IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+                /*}}}*/
+        }
+}/*}}}*/
+template <typename doubletype> void        yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
         //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
 
         IssmDouble issolid=matlitho->issolid[layer_index];
         int iy,id,ny, nindex, nstep;
-        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* dy=xNewZeroInit<IssmDouble>(36);
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
 
         /*{{{*/ /* For reference:
-         flm=(la+2.0*mu);
-         rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
-         rm0=mu/mu0;
-         rg0=g/g0;
-         frh=ro*g*ra/mu0;
-         fgr=4.0*PI*GG*ro*ra/g0;
-         fn=(deg*(deg+1.0));
-         
-
-         if(issolid==1){
-         ny = 6;
-
-         dy[0*6+0]= (-2.0*la/flm)/x;
-         dy[0*6+1]= mu0/flm;
-         dy[0*6+2]= (la*fn/flm)/x;
-         dy[0*6+3]= 0.0;
-         dy[0*6+4]= 0.0;
-         dy[0*6+5]= 0.0;
-
-         dy[1*6+0]=  4.0*(-frh+rm0*rlm/x)/x + ro*pow(omega,2.0)*ra/mu0;
-         dy[1*6+1]=(-4.0*mu/flm)/x;
-         dy[1*6+2]= fn*(frh-2.0*rm0*rlm/x)/x;
-         dy[1*6+3]= fn/x;
-         dy[1*6+4]= 0.0;
-         dy[1*6+5]= -frh/rg0;
-
-         dy[2*6+0]= -1.0/x;
-         dy[2*6+1]= 0.0;
-         dy[2*6+2]= 1.0/x;
-         dy[2*6+3]= 1/rm0;
-         dy[2*6+4]= 0.0;
-         dy[2*6+5]= 0.0;
-
-         dy[3*6+0]= (frh-2.0*rm0*rlm/x)/x;
-         dy[3*6+1]= ( -la/flm)/x;
-         dy[3*6+2]= (2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm)/(x*x) + ro*pow(omega,2.0)*ra/mu0;
-         dy[3*6+3]= -3.0/x;
-         dy[3*6+4]= -(frh/rg0)/x;
-         dy[3*6+5]= 0.0;
-
-         dy[4*6+0]= fgr;
-         dy[4*6+1]= 0.0;
-         dy[4*6+2]= 0.0;
-         dy[4*6+3]= 0.0;
-         dy[4*6+4]= 0.0;
-         dy[4*6+5]= 1.0;
-
-         dy[5*6+0]= 0.0;
-         dy[5*6+1]= 0.0;
-         dy[5*6+2]= -(fgr*fn)/x;
-         dy[5*6+3]= 0.0;
-         dy[5*6+4]= fn/(x*x);
-         dy[5*6+5]= -2.0/x;
-          
-         } else {
-         ny = 2;
-
-         dy[0*6+0]= fgr/rg0;
-         dy[0*6+1]= 1.0;
-         dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x;
-         dy[1*6+1]= -2.0/x-fgr/rg0;
-        }
-        */ /*}}}*/
-
-         if(issolid==1){
-                 ny = 6;
-         } else {
-                 ny = 2;
-         }
-         
-         for (id=0;id<ny;id++){
+                           flm=(la+2.0*mu);
+                           rlm=(3.0*la+2.0*mu)/(la+2.0*mu);
+                           rm0=mu/mu0;
+                           rg0=g/g0;
+                           frh=ro*g*ra/mu0;
+                           fgr=4.0*PI*GG*ro*ra/g0;
+                           fn=(deg*(deg+1.0));
+
+
+                           if(issolid==1){
+                           ny = 6;
+
+                           dy[0*6+0]= (-2.0*la/flm)/x;
+                           dy[0*6+1]= mu0/flm;
+                           dy[0*6+2]= (la*fn/flm)/x;
+                           dy[0*6+3]= 0.0;
+                           dy[0*6+4]= 0.0;
+                           dy[0*6+5]= 0.0;
+
+                           dy[1*6+0]=  4.0*(-frh+rm0*rlm/x)/x + ro*pow(omega,2.0)*ra/mu0;
+                           dy[1*6+1]=(-4.0*mu/flm)/x;
+                           dy[1*6+2]= fn*(frh-2.0*rm0*rlm/x)/x;
+                           dy[1*6+3]= fn/x;
+                           dy[1*6+4]= 0.0;
+                           dy[1*6+5]= -frh/rg0;
+
+                           dy[2*6+0]= -1.0/x;
+                           dy[2*6+1]= 0.0;
+                           dy[2*6+2]= 1.0/x;
+                           dy[2*6+3]= 1/rm0;
+                           dy[2*6+4]= 0.0;
+                           dy[2*6+5]= 0.0;
+
+                           dy[3*6+0]= (frh-2.0*rm0*rlm/x)/x;
+                           dy[3*6+1]= ( -la/flm)/x;
+                           dy[3*6+2]= (2.0*rm0*(la*(2.0*fn-1.0)+2.0*mu*(fn-1.0))/flm)/(x*x) + ro*pow(omega,2.0)*ra/mu0;
+                           dy[3*6+3]= -3.0/x;
+                           dy[3*6+4]= -(frh/rg0)/x;
+                           dy[3*6+5]= 0.0;
+
+                           dy[4*6+0]= fgr;
+                           dy[4*6+1]= 0.0;
+                           dy[4*6+2]= 0.0;
+                           dy[4*6+3]= 0.0;
+                           dy[4*6+4]= 0.0;
+                           dy[4*6+5]= 1.0;
+
+                           dy[5*6+0]= 0.0;
+                           dy[5*6+1]= 0.0;
+                           dy[5*6+2]= -(fgr*fn)/x;
+                           dy[5*6+3]= 0.0;
+                           dy[5*6+4]= fn/(x*x);
+                           dy[5*6+5]= -2.0/x;
+
+                           } else {
+                           ny = 2;
+
+                           dy[0*6+0]= fgr/rg0;
+                           dy[0*6+1]= 1.0;
+                           dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x;
+                           dy[1*6+1]= -2.0/x-fgr/rg0;
+                           }
+                           */ /*}}}*/
+
+        if(issolid==1){
+                ny = 6;
+        } else {
+                ny = 2;
+        }
+
+        for (id=0;id<ny;id++){
                 for (iy=0;iy<ny;iy++){
                         nindex=layer_index*nstep*36+n*36;
                         dydx[id]=dydx[id]+yi_prefactor[nindex+id*6+iy]*y[iy];
-                }
-         }
-        xDelete<IssmDouble>(dy);
-        return dydx;
-}/*}}}*/
-void propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+                }
+        }
+        return;
+}/*}}}*/
+template <typename doubletype> void        propagate_yi_euler(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
         //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
         //euler method
@@ -355 +466 @@
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
 
-        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* dydx=xNewZeroInit<IssmDouble>(6);
         IssmDouble dr = (xmax -xmin)/nstep;
-        IssmDouble x=xmin;
-    for(int i = 0;i<nstep;i++){
-        dydx=yi_derivatives(y,layer_index, i,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++){
-                y[j]+=dydx[j]*dr;
-        }
-        x = x + dr;
-    }
+        IssmDouble x=xmin;
+        for(int i = 0;i<nstep;i++){
+                yi_derivatives<doubletype>(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho);
+                for (int j=0;j<6;j++){
+                        y[j]+=dydx[j]*dr;
+                }
+                x = x + dr;
+        }
         xDelete<IssmDouble>(dydx);
 }/*}}}*/
-void propagate_yi_RK2(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+template <typename doubletype> void        propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
         //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
         //Implements Runge-Kutta 2nd order (midpoint method)
@@ -373 +484 @@
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
 
-        IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
+        doubletype k1[6]={0};
+        doubletype k2[6]={0};
+        doubletype k3[6]={0};
+        doubletype k4[6]={0};
+        doubletype y1[6]={0};
+        doubletype y2[6]={0};
+        doubletype y3[6]={0};
 
         IssmDouble dr = (xmax -xmin)/nstep;
-        IssmDouble x=xmin;
-    for(int i = 0;i<nstep/2;i++){
-        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+        IssmDouble x=xmin;
+
+        for(int i = 0;i<nstep/2;i++){
+                yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+                for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
+                yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);            
+
+                for (int j=0;j<6;j++){
+                        y[j]+=k2[j]*2.0*dr;
+                }
+                x = x + 2.0*dr;
+        }
+}/*}}}*/
+template <typename doubletype> void        propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+        //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
+        //Implements Runge-Kutta 4th order
+        int nstep;
+        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
+
+        IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
+        IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
+
+        IssmDouble dr = (xmax -xmin)/nstep;
+        IssmDouble x=xmin;
+        for(int i = 0;i<nstep/2-1;i++){
+                yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
+                for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
+                yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+                for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
+                yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+                for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
+                yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho);            
+
+                for (int j=0;j<6;j++){
+                        y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
+                }
+                x = x + 2.0*dr;
+        }
+
+        //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr
+        int i=nstep/2;
+        yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
         for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
-        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);                
-        
-        for (int j=0;j<6;j++){
-                  y[j]+=k2[j]*2.0*dr;
-        }
-        x = x + 2.0*dr;
-    }
+        yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
+        yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
+        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
+        yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);            
+
+        for (int j=0;j<6;j++){
+                y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
+        }
+
+        x = x + 2.0*dr;
+
         xDelete<IssmDouble>(k1);
         xDelete<IssmDouble>(k2);
@@ -401 +561 @@
         xDelete<IssmDouble>(y3);
 }/*}}}*/
-void propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
-        //yields this: if we have y[j]=1.0 and y[!j]=0.0 at the bottom of the layer, what is y at the top of the layer?
-        //Implements Runge-Kutta 4th order
-        int nstep;
-        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum); 
-
-        IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
-        IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
-
-        IssmDouble dr = (xmax -xmin)/nstep;
-        IssmDouble x=xmin;
-    for(int i = 0;i<nstep/2-1;i++){
-        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
-        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
-        k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
-        k4=yi_derivatives(y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho);                
-        
-        for (int j=0;j<6;j++){
-                y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
-        }
-        x = x + 2.0*dr;
-    }
-
-        //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr
-        int i=nstep/2;
-        k1=yi_derivatives(y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
-        k2=yi_derivatives(y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
-        k3=yi_derivatives(y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
-        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
-        k4=yi_derivatives(y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);                
-        
-        for (int j=0;j<6;j++){
-                y[j]+=(k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/3.0*dr;         
-        }
-
-        x = x + 2.0*dr;
-
-        xDelete<IssmDouble>(k1);
-        xDelete<IssmDouble>(k2);
-        xDelete<IssmDouble>(k3);
-        xDelete<IssmDouble>(k4);
-        xDelete<IssmDouble>(y1);
-        xDelete<IssmDouble>(y2);
-        xDelete<IssmDouble>(y3);
-}/*}}}*/
-void Innersphere_boundaryconditions(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+template <typename doubletype> void        Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
         //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5]
 
@@ -463 +568 @@
         IssmDouble  g0,r0,mu0, GG;
         int nyi;
-        femmodel->parameters->FindParam(&g0,LoveG0Enum);
-        femmodel->parameters->FindParam(&r0,LoveR0Enum);
-        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
-        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
+
+        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+
+        g0=vars->g0;
+        r0=vars->r0;
+        nyi=vars->nyi;
+
         IssmDouble ro=matlitho->density[layer_index];
         IssmDouble issolid=matlitho->issolid[layer_index];
-        IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho);
-        IssmDouble la,mu;
+        IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho,vars);
+        doubletype la,mu;
 
         if (layer_index==0){
                 // radius[0] cannot be 0 for numerical reasons, but below our first interface at radius[0] would in reality be the same material as in the first layer
-                GetEarthRheology(&la, &mu,layer_index,omega,femmodel,matlitho);   
+                GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
         } else {
-                GetEarthRheology(&la, &mu,layer_index-1,omega,femmodel,matlitho);   
+                GetEarthRheology<doubletype>(&la, &mu,layer_index-1,omega,matlitho);   
         }    
 
@@ -506 +614 @@
         yi[5+nyi*1]=-cst/(r*g0);
         yi[5+nyi*2]=deg/(r*g0);
-}/*}}}*/
-void build_yi_system(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho) { /*{{{*/
+
+}/*}}}*/
+template <typename doubletype> void        build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars) { /*{{{*/
 
         IssmDouble  g0,r0,mu0,x,ro1, GG;
         int nyi,starting_layer, nstep;
-        femmodel->parameters->FindParam(&g0,LoveG0Enum);
-        femmodel->parameters->FindParam(&r0,LoveR0Enum);
+
+        g0=vars->g0;
+        r0=vars->r0;
+        nyi=vars->nyi;
+        starting_layer=vars->starting_layer;
+
         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
-        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
-        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
 
-          IssmDouble xmin,xmax,one,ro,g;
-          IssmDouble ra=matlitho->radius[matlitho->numlayers];
+        IssmDouble xmin,xmax,one,ro,g;
+        IssmDouble ra=matlitho->radius[matlitho->numlayers];
 
         for (int i=0;i<6*(matlitho->numlayers+1);i++){
@@ -529 +640 @@
 
         int ny,is,ii,jj;
-        IssmDouble ystart[6];
+        doubletype ystart[6];
         for (int k=0;k<6;k++) ystart[k]=0.0;            
-  
+
         int ici = 0;   // Index of current interface 
         for (int i = starting_layer; i<matlitho->numlayers;i++){ 
@@ -538 +649 @@
                 xmax=(matlitho->radius[i+1])/ra;
 
-                  
+
                 if (matlitho->issolid[i]){
-                                ny = 6;
-                                is = 0;
-                                one= 1.0;
+                        ny = 6;
+                        is = 0;
+                        one= 1.0;
                 } else {        
-                                ny = 2;
-                                is = 4;
-                                one= -1.0;
-                }
-        
+                        ny = 2;
+                        is = 4;
+                        one= -1.0;
+                }
+
 
                 for (int j = 0;j<ny;j++){
                         for (int k=0;k<6;k++){ystart[k]=0.0;}
                         ystart[j]= 1.0;
-                                        
+
                         // Numerical Integration 
                         //propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
-                        propagate_yi_RK2(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
+                        propagate_yi_RK2<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
                         //propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
                         // Boundary Condition matrix - propagation part 
@@ -562 +673 @@
                         jj = 6*(ici+1)+j+is-3;
                         for (int kk=0;kk<ny;kk++){
-                                yi[(ii+kk)+nyi*jj] = ystart[kk]*one;
-                        }
-                }
-                  
-                  
+                                yi[(ii+kk)+nyi*jj] = ystart[kk]*one;
+                        }
+                }
+
+
                 // Boundary Condition matrix - solid regions
                 if (matlitho->issolid[i]){
@@ -576 +687 @@
                 } else { // Boundary Condition matrix - liquid regions
                         ro1=matlitho->density[i];
-                        g=GetGravity(matlitho->radius[i], i, femmodel,matlitho);
+                        g=GetGravity(matlitho->radius[i], i, femmodel,matlitho,vars);
                         ii = 6*ici;
                         yi[ii+nyi*(ii+3)] = -1.0;
@@ -585 +696 @@
                         yi[(ii+4)+nyi*(ii+4+3)]=-1.0;
                         yi[(ii+5)+nyi*(ii+5+3)]=-1.0;
-                        g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho);
+                        g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho,vars);
                         ii = 6*(ici+1);
                         yi[ii+nyi*(ii-1)]=-1.0;
                         yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB
                         yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB
-                          // yi(13,..) y1 CMB
+                        // yi(13,..) y1 CMB
                         yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0;
                         yi[(ii+2)+nyi*(ii)]=-1.0;
@@ -599 +710 @@
 
 
-        //-- Internal sphere: integration starts here rather than r=0 for numerical reasons
-        Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho);
-
-        //-- Surface conditions
+        //-- Internal sphere: integration starts here rather than r=0 for numerical reasons
+        Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
+
+        //-- Surface conditions
         yi[(nyi-6)+nyi*(nyi-3)]=-1.0;
         yi[(nyi-4)+nyi*(nyi-2)]=-1.0;
@@ -608 +719 @@
         yi[(nyi-1)+nyi*(nyi-1)]=deg+1.0;
 
-        //-- Degree 1 special case
+        //-- Degree 1 special case
         if(deg==1){
                 for (int i=0;i<nyi;i++){
@@ -615 +726 @@
                 yi[(nyi-1)+nyi*(nyi-1)]=1.0;
         }
-}/*}}}*/
-void yi_boundary_conditions(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+
+}/*}}}*/
+template <typename doubletype> void        yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/
 
         IssmDouble  g0,r0,mu0,ra,rb,rc;
-        int nyi, forcing_type,icb,cmb;
+        int nyi, forcing_type,icb,cmb,starting_layer;
         IssmDouble* EarthMass;
-        femmodel->parameters->FindParam(&g0,LoveG0Enum);
-        femmodel->parameters->FindParam(&r0,LoveR0Enum);
+
+        g0=vars->g0;
+        r0=vars->r0;
+        nyi=vars->nyi;
+        starting_layer=vars->starting_layer;
+        EarthMass=vars->EarthMass;
+
         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
-        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
         femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
         femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum);
         femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
-        femmodel->parameters->FindParam(&EarthMass,&matlitho->numlayers,LoveEarthMassEnum);
-          
-        // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
+
+        // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
         ra=matlitho->radius[matlitho->numlayers];       
         rb=0;
@@ -642 +757 @@
 
         IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0));
-        xDelete<IssmDouble>(EarthMass);
 
         for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0;
@@ -648 +762 @@
         switch (forcing_type) {
 
-        //-- forcings at the Inner Core Boundary
-        case 1: //'ICB --Volumetric Potential'
-                yi_righthandside[6*icb+5]=(deg)/(rc*g0);
-                yi_righthandside[6*icb+4]=1.0/(ra*g0);
-                break;
-        case 2: //'ICB --Pressure'
-                yi_righthandside[6*icb+1]=-ro_mean/mu0;
-                break;
-        case 3://'ICB --Loading'
-                yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc;
-                yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0);
-                break;
-        case 4://'ICB --Tangential Traction'
-                yi_righthandside[6*icb+3]= ro_mean/mu0;
-                break;
-
-        //--forcings at the Core Mantle Boundary
-        case 5://'CMB --Volumetric Potential'
-                yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb;
-                yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
-                break;
-        case 6://'CMB --Pressure'
-                yi_righthandside[6*cmb+1]=-ro_mean/mu0;
-                break;
-        case 7://'CMB --Loading'
-                yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb;
-                yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
-                break;
-        case 8://'CMB --Tangential Traction'
-                yi_righthandside[6*cmb+3]=-ro_mean/mu0;
-                break;
-
-        //--forcings at the surface
-        case 9://'SURF--Volumetric Potential'
-                if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0);
-                break;
-        case 10://'SURF--Pressure'
-                yi_righthandside[nyi-5]=-ro_mean/mu0;
-                break;
-        case 11://'SURF--Loading'
-                yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0);
-                if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0);
-                break;
-        case 12://'SURF--Tangential Traction'
-                yi_righthandside[nyi-3]= ro_mean/mu0;
-                break;
-        default:
-                _error_("love core error: forcing_type not supported yet");
-        }
-}/*}}}*/
-void solve_yi_system(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
+                //-- forcings at the Inner Core Boundary
+                case 1: //'ICB --Volumetric Potential'
+                        yi_righthandside[6*icb+5]=(deg)/(rc*g0);
+                        yi_righthandside[6*icb+4]=1.0/(ra*g0);
+                        break;
+                case 2: //'ICB --Pressure'
+                        yi_righthandside[6*icb+1]=-ro_mean/mu0;
+                        break;
+                case 3://'ICB --Loading'
+                        yi_righthandside[6*icb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rc;
+                        yi_righthandside[6*icb+5]= (2.0*deg+1.0)/(rc*g0);
+                        break;
+                case 4://'ICB --Tangential Traction'
+                        yi_righthandside[6*icb+3]= ro_mean/mu0;
+                        break;
+
+                        //--forcings at the Core Mantle Boundary
+                case 5://'CMB --Volumetric Potential'
+                        yi_righthandside[6*cmb+1]=-ro_mean/mu0*ra/rb;
+                        yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
+                        break;
+                case 6://'CMB --Pressure'
+                        yi_righthandside[6*cmb+1]=-ro_mean/mu0;
+                        break;
+                case 7://'CMB --Loading'
+                        yi_righthandside[6*cmb+1]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0)*ra/rb;
+                        yi_righthandside[6*cmb+5]= (2.0*deg+1.0)/(rb*g0);
+                        break;
+                case 8://'CMB --Tangential Traction'
+                        yi_righthandside[6*cmb+3]=-ro_mean/mu0;
+                        break;
+
+                        //--forcings at the surface
+                case 9://'SURF--Volumetric Potential'
+                        if (deg>1) yi_righthandside[nyi-1]=(2.0*deg+1.0)/(ra*g0);
+                        break;
+                case 10://'SURF--Pressure'
+                        yi_righthandside[nyi-5]=-ro_mean/mu0;
+                        break;
+                case 11://'SURF--Loading'
+                        yi_righthandside[nyi-5]=-ro_mean*(2.0*deg+1.0)/(3.0*mu0);
+                        if (deg>1) yi_righthandside[nyi-1]= (2.0*deg+1.0)/(ra*g0);
+                        break;
+                case 12://'SURF--Tangential Traction'
+                        yi_righthandside[nyi-3]= ro_mean/mu0;
+                        break;
+                default:
+                        _error_("love core error: forcing_type not supported yet");
+        }
+}/*}}}*/
+template <typename doubletype> void        solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
 
         IssmDouble  g0,r0,mu0,loveratio,underflow_tol;
@@ -704 +818 @@
         int nyi,starting_layer, dummy;
         bool allow_layer_deletion;
-        femmodel->parameters->FindParam(&g0,LoveG0Enum);
-        femmodel->parameters->FindParam(&r0,LoveR0Enum);
+        IssmDouble* EarthMass=NULL;
+
+        g0=vars->g0;
+        r0=vars->r0;
+        nyi=vars->nyi;
+        starting_layer=vars->starting_layer;
+        EarthMass=vars->EarthMass;
+
         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
-        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
-        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
         femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
         femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum);
         femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
-        IssmDouble ra=matlitho->radius[matlitho->numlayers];
+        IssmDouble ra=matlitho->radius[matlitho->numlayers];
         bool exit=false;
         int lda,ldb;
@@ -719 +837 @@
                 lda=nyi;
                 ldb=nyi;
-                IssmDouble*  yilocal=xNew<IssmDouble>(nyi*nyi); // we will need to redeclare these inside here as nyi changes
-                IssmDouble*  rhslocal=xNew<IssmDouble>(nyi);
+                doubletype*  yilocal=xNew<doubletype>(nyi*nyi); // we will need to redeclare these inside here as nyi changes
+                doubletype*  rhslocal=xNew<doubletype>(nyi);
 
                 //we need to do a local copy of yi,rhs to send them to LAPACK with the appropriate size and to keep the original matrices in case layers are deleted and we need to try again
@@ -734 +852 @@
                 int nrhs=1; // number of right hand size columns
 
-                dgesv_(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info);
+                allgesv<doubletype>(&nyi, &nrhs, yilocal, &lda, ipiv, rhslocal, &ldb, &info);
+                
                 xDelete<int>(ipiv);
-
-                if(VerboseSolution() && info!=0){ 
+                
+                if(VerboseModule() && info!=0){ 
                         _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
                         printf("%s %s %s %s\n", "i","j","yi[i+nyi*j]","rhs[i]");
@@ -751 +870 @@
                 *lovek = rhslocal[nyi-1]*ra*g0;
 
-                IssmDouble loveh1 = rhslocal[3];
-                IssmDouble lovel1 = rhslocal[5];
-                IssmDouble lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra);
-
-                IssmDouble loveh1s = rhslocal[nyi-3];
-                IssmDouble lovel1s = rhslocal[nyi-2];
-                IssmDouble lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra);
+                doubletype loveh1 = rhslocal[3];
+                doubletype lovel1 = rhslocal[5];
+                doubletype lovek1 = rhslocal[7] - pow(matlitho->radius[starting_layer]/ra,deg)/(g0*ra);
+
+                doubletype loveh1s = rhslocal[nyi-3];
+                doubletype lovel1s = rhslocal[nyi-2];
+                doubletype lovek1s = rhslocal[nyi-1] - 1.0/(g0*ra);
 
                 loveratio = abs(loveh1/loveh1s); //ratio of center to surface love numbers, determines if we should remove layers
@@ -763 +882 @@
                 if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s);
 
-                if (!allow_layer_deletion || nyi<=12 || omega!=2.0*PI*frequencies[0]){ //We are not allowed to delete layers, or there is only one layer left. We also don't want to delete layers in the middle of a loop on frequencies, as that can lead to a jump that would compromise the inverse laplace transform.
+                if (!allow_layer_deletion || nyi<=12 || omega!=angular_frequency<doubletype>(frequencies[0])){ 
                         goto save_results;
+                        /*We are not allowed to delete layers, or there is only one layer left. We also don't want to delete 
+                          layers in the middle of a loop on frequencies, as that can lead to a jump that would compromise the 
+                          inverse laplace transform.*/
                 }
 
@@ -775 +897 @@
                         nyi-=6;
                         starting_layer+=1;
-                        femmodel->parameters->SetParam(nyi,LoveNYiEquationsEnum);
-                        femmodel->parameters->SetParam(starting_layer,LoveStartingLayerEnum);
-
+                        vars->nyi=nyi;
+                        vars->starting_layer=starting_layer;
 
                         for (int i=0;i<nyi;i++){//shift everything down by 1 layer
@@ -786 +907 @@
                         }
 
-                        Innersphere_boundaryconditions(yi, starting_layer, deg, omega, femmodel, matlitho); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different
+                        Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars); //we move the first interface to the new starting layer. yi[0:2,0:5] will be different
                 } else { //we are ready to save the outputs and break the main loop
 
-                save_results:
+save_results:
                         for (int i=0;i<nyi;i++){
                                 rhs[i]=rhslocal[i];
@@ -811 +932 @@
 
                         exit = true;
-                }
-                xDelete<IssmDouble>(yilocal);
-                xDelete<IssmDouble>(rhslocal);
-        }
+                }
+                xDelete<doubletype>(yilocal);
+                xDelete<doubletype>(rhslocal);
+        }
         xDelete<IssmDouble>(frequencies);       
-}/*}}}*/
-IssmDouble factorial(int n){ /*{{{*/
-        IssmDouble prod=1;
-        for (int i=2;i<n+1;i++) prod*=i;
-        return prod;
-}/*}}}*/
-IssmDouble n_C_r(int n, int r){ /*{{{*/ 
-        //n choose r
-        int primes[169] = 
-  {2,    3,    5,    7,   11,   13,   17,   19,   23,   29,  
-  31,   37,   41,   43,   47,   53,   59,   61,   67,   71,  
-  73,   79,   83,   89,   97,  101,  103,  107,  109,  113,  
- 127,  131,  137,  139,  149,  151,  157,  163,  167,  173,  
- 179,  181,  191,  193,  197,  199,  211,  223,  227,  229,  
- 233,  239,  241,  251,  257,  263,  269,  271,  277,  281,  
- 283,  293,  307,  311,  313,  317,  331,  337,  347,  349,  
- 353,  359,  367,  373,  379,  383,  389,  397,  401,  409,  
- 419,  421,  431,  433,  439,  443,  449,  457,  461,  463,  
- 467,  479,  487,  491,  499,  503,  509,  521,  523,  541,  
- 547,  557,  563,  569,  571,  577,  587,  593,  599,  601,  
- 607,  613,  617,  619,  631,  641,  643,  647,  653,  659,  
- 661,  673,  677,  683,  691,  701,  709,  719,  727,  733,  
- 739,  743,  751,  757,  761,  769,  773,  787,  797,  809,  
- 811,  821,  823,  827,  829,  839,  853,  857,  859,  863,  
- 877,  881,  883,  887,  907,  911,  919,  929,  937,  941,  
- 947,  953,  967,  971,  977,  983,  991,  997, 1009};
-        int num, den;
-        num = 1;
-        den = 1;
-
-        for (int i=0;i<r;i++){
-                num = num*(n-i);
-                den = den*(i+1);
-                if (i>0) {
-                        // Divide out common prime factors
-                        for (int k=0;k<169;k++){ //169 is the length of the prime vector here
-                                if ( i % primes[k] == 0) { // modulo
-                                        num = num/primes[k];
-                                        den = den/primes[k];
-                                }
-                        }
-                }
-        }
-
-        IssmDouble res;        
-        return res = num/den;
-}/*}}}*/
-IssmDouble* postwidder_coef(int NTit){ /*{{{*/
-        //Coefficients of the Post-Widder method through Saltzer summation for inverse Laplace transform:
-        //The Mth iteration estimate will be: f(t)_M = sum_{k=1:2*M}(xi_[M,k] * f(s_k))
-        //The method is based on equations (2), (6), (7) in: 
-        //Valko PP, Abate J. Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion. Computational Mathematics and Applications. (2004)
-        //Note that the coefficients xi lack the factor s=k*log(2)/t. 
-        //That is because we are computing the heaviside response of the system rather than its impulse response, 
-        //and Laplace_Transform(Heaviside(t).*f(t)) = f(s)/s. So s cancels out in the sum for f(t)_M.
-
-        IssmDouble* xi=xNewZeroInit<IssmDouble>(2*NTit*NTit);
-        int indxi;
-        for (int M=1;M<NTit+1;M++){
-                for (int k=1;k<2*M+1;k++){
-                indxi=(M-1)*(2*NTit)+k-1;
-                    for (int j=floor((k+1)/2);j<min(k,M)+1;j++){
-                            xi[indxi]+=pow(j,M+1.0)/factorial(M)*n_C_r(M,j)*n_C_r(2*j,j)*n_C_r(j,k-j);
-                    }
-                    xi[indxi]*=pow(-1.0,k+M)/k;
-                }
-        }
-        return xi;
-}/*}}}*/
-void postwidder_transform(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel){ /*{{{*/
-        //Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef
-
-        int nfreq, indxi, indf;
-        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
-        int nt=nfreq/2/NTit;
-
-        indf=d*nfreq+t*2*NTit;
-        IssmDouble LoveM[NTit];
-
-        for (int M=1;M<NTit+1;M++){
-                LoveM[M-1]=0.0;
-                for (int k=1;k<2*M+1;k++){
-                        indxi=(M-1)*(2*NTit)+k-1;
-                        LoveM[M-1]+=xi[indxi]*Lovef[indf+k-1];
-                }
-                
-                //Make sure we are not getting into numerical instability
-                //Diverging once: ok, we'll give that the benefit of the doubt, it can be an inflexion point in the convergence series
-                //Diverging twice in a row: we are definitely propagating numerical error: revert to the last stable value and exit
-                if (M>3){ 
-                        if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) &&
-                             abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){
-                                Lovet[d*nt+t]=LoveM[M-3];
-                                return;
-                        }
-                }
-        }
-        Lovet[d*nt+t]=LoveM[NTit-1];
-}/*}}}*/
-void love_freq_to_temporal(IssmDouble* LoveHt,IssmDouble* LoveLt,IssmDouble* LoveKt,IssmDouble* LoveHf,IssmDouble* LoveLf,IssmDouble* LoveKf, FemModel* femmodel){ /*{{{*/
+
+}/*}}}*/
+template <typename doubletype> void        love_freq_to_temporal(doubletype* LoveHt,doubletype* LoveLt,doubletype* LoveKt,doubletype* LoveHf,doubletype* LoveLf,doubletype* LoveKf, FemModel* femmodel){ /*{{{*/
         //Transforms all frequency-dependent love numbers into time-dependent love numbers
         int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit;
-        IssmDouble meanh,meanl,meank,dh,dl,dk;
 
         femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
@@ -924 +947 @@
         femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
         femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
-        
+
         int nt=nfreq/2/NTit;
-        IssmDouble* xi=postwidder_coef(NTit);
+        doubletype* xi=postwidder_coef<doubletype>(NTit);
 
         for (int d=sh_nmin;d<sh_nmax+1;d++){
                 for (int t=0;t<nt;t++){
-                        postwidder_transform(LoveHt,LoveHf,d,t,NTit,xi,femmodel);
-                        postwidder_transform(LoveLt,LoveLf,d,t,NTit,xi,femmodel);
-                        postwidder_transform(LoveKt,LoveKf,d,t,NTit,xi,femmodel);
-                }
-        }
-        xDelete<IssmDouble>(xi);
-}/*}}}*/
-void love_core(FemModel* femmodel){ /*{{{*/
+                        postwidder_transform<doubletype>(LoveHt,LoveHf,d,t,NTit,xi,femmodel);
+                        postwidder_transform<doubletype>(LoveLt,LoveLf,d,t,NTit,xi,femmodel);
+                        postwidder_transform<doubletype>(LoveKt,LoveKf,d,t,NTit,xi,femmodel);
+                }
+        }
+        xDelete<doubletype>(xi);
+}/*}}}*/
+
+/*templated cores:*/
+LoveVariables*                             love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/
+
+        /*initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system*/
+
+        bool        verbosemod = (int)VerboseModule();
+        int         numlayers  = matlitho->numlayers;
+        IssmDouble* r          = matlitho->radius;
+        IssmDouble  r1,r2,ro, GG;
+        
+        /*outputs:*/
+        IssmDouble* EarthMass=NULL;
+        IssmDouble  g0,r0;
+        int         nyi,starting_layer;
+
+        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
+        EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
+
+        for (int i=0;i<numlayers;i++){
+                r2 = r[i+1];
+                ro = matlitho->density[i];
+                if (i==0){
+                        EarthMass[i] = ro*pow(r2,3.0)*4.0*PI/3.0;
+                }else{
+                        r1=r[i];
+                        EarthMass[i] = EarthMass[i-1] + ro*(pow(r2,3.0)-pow(r1,3.0))*4.0*PI/3.0;;
+                }
+        }
+
+        g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0);
+        r0=r[numlayers];
+        nyi=6*(numlayers+1);
+        starting_layer=0;
+        
+        return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer);
+
+} /*}}}*/
+template <typename doubletype> void        love_core_template(FemModel* femmodel){ /*{{{*/
+
+        Matlitho*   matlitho=NULL;
+        int         nfreq,nyi,starting_layer,nstep,forcing_type,dummy;
+        int         sh_nmin,sh_nmax,kernel_index,deleted_layer_offset;
+        bool        allow_layer_deletion,love_kernels, istemporal;
+        bool        verbosemod = (int)VerboseModule();
+        IssmDouble  g0,r0,mu0, underflow_tol;
+        IssmDouble *frequencies = NULL;
+        bool        save_results;
+        bool        complex_computation;
+
+        doubletype  omega;
+        doubletype  lovek, loveh, lovel, loveratio;
+
+        doubletype*  LoveKf = NULL;
+        doubletype*  LoveHf = NULL;
+        doubletype*  LoveLf = NULL;
+        doubletype*  LoveKernels= NULL;
+        LoveVariables* vars=NULL;
+        doubletype* yi_prefactor=NULL;
+        doubletype* yi_righthandside=NULL;
+        doubletype* yi=NULL;
+
+        if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");
 
         /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
-        Matlitho* matlitho=NULL;
         for (Object* & object: femmodel->materials->objects){
                 Material* material=xDynamicCast<Material*>(object);
@@ -950 +1034 @@
         _assert_(matlitho);
 
-        /*love parameters: */
-        int         nfreq,nyi,starting_layer,nstep,forcing_type,dummy;
-        int         sh_nmin,sh_nmax,kernel_index,deleted_layer_offset;
-        bool        allow_layer_deletion,love_kernels, istemporal;
-        bool        verbosemod = (int)VerboseModule();
-        IssmDouble lovek, loveh, lovel, loveratio;
-        IssmDouble  g0,r0,mu0, underflow_tol, omega;
-        IssmDouble *frequencies = NULL;
-
-        /*parameters: */
-        bool save_results;
-
-        if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");
-
-        /*Recover some parameters: */
+        /*recover parameters: */
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
-
-        /*recover love number parameters: */
-
         femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
         femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
@@ -981 +1048 @@
         femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
+        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
+
 
         /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
-        IssmDouble*  LoveKf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
-        IssmDouble*  LoveHf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
-        IssmDouble*  LoveLf = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
+        LoveKf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
+        LoveHf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
+        LoveLf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
 
         /*Initialize love kernels (real and imaginary parts): */
-        IssmDouble*  LoveKernelsReal = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
-        IssmDouble*  LoveKernelsImag = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
-
-        love_init(femmodel,matlitho);
-        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum);
-        femmodel->parameters->FindParam(&starting_layer,LoveStartingLayerEnum);
-        IssmDouble* yi_prefactor=xNewZeroInit<IssmDouble>(6*6*nstep*matlitho->numlayers);
-        IssmDouble* yi_righthandside=xNewZeroInit<IssmDouble>(nyi);
-        IssmDouble* yi=xNewZeroInit<IssmDouble>(nyi*nyi);
+        LoveKernels= xNewZeroInit<doubletype>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
+
+        vars=love_init(femmodel,matlitho);
+
+        yi_prefactor=xNewZeroInit<doubletype>(6*6*nstep*matlitho->numlayers);
+        yi_righthandside=xNewZeroInit<doubletype>(vars->nyi);
+        yi=xNewZeroInit<doubletype>(vars->nyi*vars->nyi);
 
         //precalculate yi coefficients that do not depend on degree or frequency
-        fill_yi_prefactor(yi_prefactor, NULL, NULL,femmodel, matlitho); 
+        fill_yi_prefactor<doubletype>(yi_prefactor, NULL, NULL,femmodel, matlitho,vars); 
         for(int deg=sh_nmin;deg<sh_nmax+1;deg++){
 
                 //precalculate yi coefficients that depend on degree but not frequency
-                fill_yi_prefactor(yi_prefactor, &deg, NULL,femmodel, matlitho); 
+                fill_yi_prefactor<doubletype>(yi_prefactor, &deg, NULL,femmodel, matlitho,vars); 
 
                 for (int fr=0;fr<nfreq;fr++){
-                        omega=2.0*PI*frequencies[fr]; //angular frequency
+                        omega=angular_frequency<doubletype>(frequencies[fr]);
 
                         //precalculate yi coefficients that depend on degree and frequency
-                        fill_yi_prefactor(yi_prefactor, &deg,&omega,femmodel, matlitho); 
-
-                        yi_boundary_conditions(yi_righthandside,deg,femmodel,matlitho); 
-
-                        build_yi_system(yi,deg,omega,yi_prefactor,femmodel,matlitho);
-
-                        solve_yi_system(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho);
+                        fill_yi_prefactor<doubletype>(yi_prefactor, &deg,&omega,femmodel, matlitho,vars); 
+
+                        yi_boundary_conditions<doubletype>(yi_righthandside,deg,femmodel,matlitho,vars); 
+
+                        build_yi_system<doubletype>(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars);
+
+                        solve_yi_system<doubletype>(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho,vars);
 
                         LoveHf[deg*nfreq+fr]=loveh;
@@ -1021 +1088 @@
                         LoveLf[deg*nfreq+fr]=lovel;
 
-                        femmodel->parameters->FindParam(&nyi,LoveNYiEquationsEnum); //synchronize nyi in case we deleted a layer
-                        deleted_layer_offset=(matlitho->numlayers+1)*6-nyi;// =6 per deleted layer
-                        kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 +
-                                            fr*(matlitho->numlayers+1)*6 +
-                                                deleted_layer_offset;
-                        for (int i=0;i<nyi;i++){
-                                LoveKernelsReal[kernel_index+i]=yi_righthandside[i];
-                        }
-                }
-        }
-
-
-        xDelete<IssmDouble>(yi);
-        xDelete<IssmDouble>(yi_righthandside);
+                        deleted_layer_offset=(matlitho->numlayers+1)*6-vars->nyi;// =6 per deleted layer
+                        kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 + fr*(matlitho->numlayers+1)*6 + deleted_layer_offset;
+                        for (int i=0;i<vars->nyi;i++){
+                                LoveKernels[kernel_index+i]=yi_righthandside[i];
+                        }
+                }
+        }
+
+
+        xDelete<doubletype>(yi);
+        xDelete<doubletype>(yi_righthandside);
 
 
         //Temporal love numbers
-        if (istemporal){
+        if (istemporal && !complex_computation){
                 int NTit;
                 femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
                 int nt = nfreq/2/NTit;
 
-                IssmDouble* LoveHt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
-                IssmDouble* LoveLt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
-                IssmDouble* LoveKt=xNewZeroInit<IssmDouble>((sh_nmax+1)*nt);
-
-                love_freq_to_temporal(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
-
-                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
-                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
-                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0));
-
-                xDelete<IssmDouble>(LoveHt);
-                xDelete<IssmDouble>(LoveLt);
-                xDelete<IssmDouble>(LoveKt);
-        }
-
-        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
-        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));
-        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
+                doubletype* LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+                doubletype* LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+                doubletype* LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
+
+                love_freq_to_temporal<doubletype>(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
+
+                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));
+                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0));
+
+                xDelete<doubletype>(LoveHt);
+                xDelete<doubletype>(LoveLt);
+                xDelete<doubletype>(LoveKt);
+        }
+
+        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));
+        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
 
 
@@ -1082 +1146 @@
         /*Only when love_kernels is on*/
         if (love_kernels==1) {
-                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
+                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
                 //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
         }
@@ -1088 +1152 @@
         /*Free resources:*/
         xDelete<IssmDouble>(frequencies);
-        xDelete<IssmDouble>(LoveKf);
-        xDelete<IssmDouble>(LoveHf);
-        xDelete<IssmDouble>(LoveLf);
-        //xDelete<IssmDouble>(LoveKi);
-        //xDelete<IssmDouble>(LoveHi);
+        xDelete<doubletype>(LoveKf);
+        xDelete<doubletype>(LoveHf);
+        xDelete<doubletype>(LoveLf);
         //xDelete<IssmDouble>(LoveLi);
-        xDelete<IssmDouble>(LoveKernelsReal);
-        xDelete<IssmDouble>(LoveKernelsImag);
+        xDelete<doubletype>(LoveKernels);
 } /*}}}*/
+
+/*cores and template instantiations:*/
+/*template instantiations :{{{*/
+template void love_core_template<IssmDouble>(FemModel* femmodel);
+template void love_core_template<IssmComplex>(FemModel* femmodel);
+template void        fill_yi_prefactor<IssmDouble>(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template void        fill_yi_prefactor<IssmComplex>(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template void        GetEarthRheology<IssmDouble>(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,  Matlitho* matlitho);
+template void        GetEarthRheology<IssmComplex>(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega,  Matlitho* matlitho);
+template void        yi_boundary_conditions<IssmDouble>(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
+template void        yi_boundary_conditions<IssmComplex>(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
+template void        yi_derivatives<IssmDouble>(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        yi_derivatives<IssmComplex>(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_RK2<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        propagate_yi_RK2<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
+template void        Innersphere_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template void        Innersphere_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template void        build_yi_system<IssmDouble>(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
+template void        build_yi_system<IssmComplex>(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
+template void        solve_yi_system<IssmDouble>(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template void        solve_yi_system<IssmComplex>(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
+template IssmDouble  factorial<IssmDouble>(int n);
+template IssmDouble* postwidder_coef<IssmDouble>(int NTit);
+template IssmDouble  n_C_r<IssmDouble>(int n, int r);
+template void         postwidder_transform<IssmDouble>(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel);
+
+/*}}}*/
+void           love_core(FemModel* femmodel){ /*{{{*/
+
+        bool        complex_computation;
+
+        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
+
+        if(complex_computation) love_core_template<IssmComplex>(femmodel);
+        else                    love_core_template<IssmDouble>(femmodel);
+
+} /*}}}*/

issm/trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -256 +256 @@
 syn keyword cConstant LockFileNameEnum
 syn keyword cConstant LoveAllowLayerDeletionEnum
+syn keyword cConstant LoveCoreMantleBoundaryEnum
+syn keyword cConstant LoveEarthMassEnum
 syn keyword cConstant LoveForcingTypeEnum
 syn keyword cConstant LoveFrequenciesEnum
+syn keyword cConstant LoveIsTemporalEnum
 syn keyword cConstant LoveG0Enum
+syn keyword cConstant LoveGravitationalConstantEnum
+syn keyword cConstant LoveInnerCoreBoundaryEnum
+syn keyword cConstant LoveComplexComputationEnum
+syn keyword cConstant LoveIntStepsPerLayerEnum
 syn keyword cConstant LoveKernelsEnum
 syn keyword cConstant LoveMu0Enum
 syn keyword cConstant LoveNfreqEnum
+syn keyword cConstant LoveNTemporalIterationsEnum
+syn keyword cConstant LoveNYiEquationsEnum
 syn keyword cConstant LoveR0Enum
 syn keyword cConstant LoveShNmaxEnum
 syn keyword cConstant LoveShNminEnum
+syn keyword cConstant LoveStartingLayerEnum
+syn keyword cConstant LoveUnderflowTolEnum
 syn keyword cConstant MassFluxSegmentsEnum
 syn keyword cConstant MassFluxSegmentsPresentEnum
@@ -783 +794 @@
 syn keyword cConstant SealevelGrotm2Enum
 syn keyword cConstant SealevelGrotm3Enum
+syn keyword cConstant SealevelGUrotm1Enum
+syn keyword cConstant SealevelGUrotm2Enum
+syn keyword cConstant SealevelGUrotm3Enum
+syn keyword cConstant SealevelGNrotm1Enum
+syn keyword cConstant SealevelGNrotm2Enum
+syn keyword cConstant SealevelGNrotm3Enum
+syn keyword cConstant SealevelGErotm1Enum
+syn keyword cConstant SealevelGErotm2Enum
+syn keyword cConstant SealevelGErotm3Enum
 syn keyword cConstant SealevelRSLBarystaticEnum
 syn keyword cConstant SealevelRSLRateEnum
@@ -1528 +1548 @@
 syn keyword cType GenericParam
 syn keyword cType GiaDeflectionCoreArgs
+syn keyword cType GrdLoads
 syn keyword cType Hook
 syn keyword cType Input

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -258 +258 @@
         LoveGravitationalConstantEnum,
         LoveInnerCoreBoundaryEnum,
+        LoveComplexComputationEnum,
         LoveIntStepsPerLayerEnum,
         LoveKernelsEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -258 +258 @@
                 case LockFileNameEnum : return "LockFileName";
                 case LoveAllowLayerDeletionEnum : return "LoveAllowLayerDeletion";
+                case LoveCoreMantleBoundaryEnum : return "LoveCoreMantleBoundary";
+                case LoveEarthMassEnum : return "LoveEarthMass";
                 case LoveForcingTypeEnum : return "LoveForcingType";
                 case LoveFrequenciesEnum : return "LoveFrequencies";
+                case LoveIsTemporalEnum : return "LoveIsTemporal";
                 case LoveG0Enum : return "LoveG0";
+                case LoveGravitationalConstantEnum : return "LoveGravitationalConstant";
+                case LoveInnerCoreBoundaryEnum : return "LoveInnerCoreBoundary";
+                case LoveComplexComputationEnum : return "LoveComplexComputation";
+                case LoveIntStepsPerLayerEnum : return "LoveIntStepsPerLayer";
                 case LoveKernelsEnum : return "LoveKernels";
                 case LoveMu0Enum : return "LoveMu0";
                 case LoveNfreqEnum : return "LoveNfreq";
+                case LoveNTemporalIterationsEnum : return "LoveNTemporalIterations";
+                case LoveNYiEquationsEnum : return "LoveNYiEquations";
                 case LoveR0Enum : return "LoveR0";
                 case LoveShNmaxEnum : return "LoveShNmax";
                 case LoveShNminEnum : return "LoveShNmin";
+                case LoveStartingLayerEnum : return "LoveStartingLayer";
+                case LoveUnderflowTolEnum : return "LoveUnderflowTol";
                 case MassFluxSegmentsEnum : return "MassFluxSegments";
                 case MassFluxSegmentsPresentEnum : return "MassFluxSegmentsPresent";
@@ -785 +796 @@
                 case SealevelGrotm2Enum : return "SealevelGrotm2";
                 case SealevelGrotm3Enum : return "SealevelGrotm3";
+                case SealevelGUrotm1Enum : return "SealevelGUrotm1";
+                case SealevelGUrotm2Enum : return "SealevelGUrotm2";
+                case SealevelGUrotm3Enum : return "SealevelGUrotm3";
+                case SealevelGNrotm1Enum : return "SealevelGNrotm1";
+                case SealevelGNrotm2Enum : return "SealevelGNrotm2";
+                case SealevelGNrotm3Enum : return "SealevelGNrotm3";
+                case SealevelGErotm1Enum : return "SealevelGErotm1";
+                case SealevelGErotm2Enum : return "SealevelGErotm2";
+                case SealevelGErotm3Enum : return "SealevelGErotm3";
                 case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic";
                 case SealevelRSLRateEnum : return "SealevelRSLRate";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -264 +264 @@
    if(stage==3){
               if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
+              else if (strcmp(name,"LoveCoreMantleBoundary")==0) return LoveCoreMantleBoundaryEnum;
+              else if (strcmp(name,"LoveEarthMass")==0) return LoveEarthMassEnum;
               else if (strcmp(name,"LoveForcingType")==0) return LoveForcingTypeEnum;
               else if (strcmp(name,"LoveFrequencies")==0) return LoveFrequenciesEnum;
+              else if (strcmp(name,"LoveIsTemporal")==0) return LoveIsTemporalEnum;
               else if (strcmp(name,"LoveG0")==0) return LoveG0Enum;
+              else if (strcmp(name,"LoveGravitationalConstant")==0) return LoveGravitationalConstantEnum;
+              else if (strcmp(name,"LoveInnerCoreBoundary")==0) return LoveInnerCoreBoundaryEnum;
+              else if (strcmp(name,"LoveComplexComputation")==0) return LoveComplexComputationEnum;
+              else if (strcmp(name,"LoveIntStepsPerLayer")==0) return LoveIntStepsPerLayerEnum;
               else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum;
               else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
               else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum;
+              else if (strcmp(name,"LoveNTemporalIterations")==0) return LoveNTemporalIterationsEnum;
+              else if (strcmp(name,"LoveNYiEquations")==0) return LoveNYiEquationsEnum;
               else if (strcmp(name,"LoveR0")==0) return LoveR0Enum;
               else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum;
               else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum;
+              else if (strcmp(name,"LoveStartingLayer")==0) return LoveStartingLayerEnum;
+              else if (strcmp(name,"LoveUnderflowTol")==0) return LoveUnderflowTolEnum;
               else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
               else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
@@ -372 +383 @@
               else if (strcmp(name,"SealevelchangeGRigid")==0) return SealevelchangeGRigidEnum;
               else if (strcmp(name,"SealevelchangeGElastic")==0) return SealevelchangeGElasticEnum;
-              else if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"SolidearthSettingsComputesealevelchange")==0) return SolidearthSettingsComputesealevelchangeEnum;
               else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
               else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
@@ -383 +397 @@
               else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
               else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
+              else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
               else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
               else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
@@ -495 +506 @@
               else if (strcmp(name,"TimesteppingType")==0) return TimesteppingTypeEnum;
               else if (strcmp(name,"ToMITgcmComm")==0) return ToMITgcmCommEnum;
-              else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
+         else stage=5;
+   }
+   if(stage==5){
+              if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
               else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
               else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
@@ -506 +520 @@
               else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
               else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
-         else stage=5;
-   }
-   if(stage==5){
-              if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
+              else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
               else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
               else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum;
@@ -618 +629 @@
               else if (strcmp(name,"DistanceToCalvingfront")==0) return DistanceToCalvingfrontEnum;
               else if (strcmp(name,"DistanceToGroundingline")==0) return DistanceToGroundinglineEnum;
-              else if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
+         else stage=6;
+   }
+   if(stage==6){
+              if (strcmp(name,"Domain2Dhorizontal")==0) return Domain2DhorizontalEnum;
               else if (strcmp(name,"Domain2Dvertical")==0) return Domain2DverticalEnum;
               else if (strcmp(name,"Domain3D")==0) return Domain3DEnum;
@@ -629 +643 @@
               else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum;
               else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
-         else stage=6;
-   }
-   if(stage==6){
-              if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
+              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
               else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
               else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
@@ -741 +752 @@
               else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
-              else if (strcmp(name,"Node")==0) return NodeEnum;
+         else stage=7;
+   }
+   if(stage==7){
+              if (strcmp(name,"Node")==0) return NodeEnum;
               else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum;
               else if (strcmp(name,"OceantransportSpcbottompressure")==0) return OceantransportSpcbottompressureEnum;
@@ -752 +766 @@
               else if (strcmp(name,"Radar")==0) return RadarEnum;
               else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum;
-         else stage=7;
-   }
-   if(stage==7){
-              if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
+              else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
               else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
               else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
@@ -803 +814 @@
               else if (strcmp(name,"SealevelGrotm2")==0) return SealevelGrotm2Enum;
               else if (strcmp(name,"SealevelGrotm3")==0) return SealevelGrotm3Enum;
+              else if (strcmp(name,"SealevelGUrotm1")==0) return SealevelGUrotm1Enum;
+              else if (strcmp(name,"SealevelGUrotm2")==0) return SealevelGUrotm2Enum;
+              else if (strcmp(name,"SealevelGUrotm3")==0) return SealevelGUrotm3Enum;
+              else if (strcmp(name,"SealevelGNrotm1")==0) return SealevelGNrotm1Enum;
+              else if (strcmp(name,"SealevelGNrotm2")==0) return SealevelGNrotm2Enum;
+              else if (strcmp(name,"SealevelGNrotm3")==0) return SealevelGNrotm3Enum;
+              else if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
+              else if (strcmp(name,"SealevelGErotm2")==0) return SealevelGErotm2Enum;
+              else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum;
               else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
               else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
@@ -855 +875 @@
               else if (strcmp(name,"SmbD")==0) return SmbDEnum;
               else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
-              else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
+         else stage=8;
+   }
+   if(stage==8){
+              if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
               else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
               else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum;
@@ -875 +898 @@
               else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
               else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
-         else stage=8;
-   }
-   if(stage==8){
-              if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
+              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
               else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
               else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
@@ -978 +998 @@
               else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
               else if (strcmp(name,"SurfaceSlopeY")==0) return SurfaceSlopeYEnum;
-              else if (strcmp(name,"Temperature")==0) return TemperatureEnum;
+         else stage=9;
+   }
+   if(stage==9){
+              if (strcmp(name,"Temperature")==0) return TemperatureEnum;
               else if (strcmp(name,"TemperaturePDD")==0) return TemperaturePDDEnum;
               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
@@ -998 +1021 @@
               else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
               else if (strcmp(name,"VxObs")==0) return VxObsEnum;
-         else stage=9;
-   }
-   if(stage==9){
-              if (strcmp(name,"VxShear")==0) return VxShearEnum;
+              else if (strcmp(name,"VxShear")==0) return VxShearEnum;
               else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
@@ -1101 +1121 @@
               else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
               else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
-              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
+         else stage=10;
+   }
+   if(stage==10){
+              if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
               else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum;
               else if (strcmp(name,"Outputdefinition81")==0) return Outputdefinition81Enum;
@@ -1121 +1144 @@
               else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
               else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
-         else stage=10;
-   }
-   if(stage==10){
-              if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
+              else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
               else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
               else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
@@ -1224 +1244 @@
               else if (strcmp(name,"FileParam")==0) return FileParamEnum;
               else if (strcmp(name,"FixedTimestepping")==0) return FixedTimesteppingEnum;
-              else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
+         else stage=11;
+   }
+   if(stage==11){
+              if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
               else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
               else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
@@ -1244 +1267 @@
               else if (strcmp(name,"Gradient1")==0) return Gradient1Enum;
               else if (strcmp(name,"Gradient2")==0) return Gradient2Enum;
-         else stage=11;
-   }
-   if(stage==11){
-              if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
+              else if (strcmp(name,"Gradient3")==0) return Gradient3Enum;
               else if (strcmp(name,"Gradient4")==0) return Gradient4Enum;
               else if (strcmp(name,"GroundedArea")==0) return GroundedAreaEnum;
@@ -1347 +1367 @@
               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
-              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
+         else stage=12;
+   }
+   if(stage==12){
+              if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
               else if (strcmp(name,"Moulin")==0) return MoulinEnum;
               else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
@@ -1367 +1390 @@
               else if (strcmp(name,"OceantransportSolution")==0) return OceantransportSolutionEnum;
               else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
-         else stage=12;
-   }
-   if(stage==12){
-              if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
+              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
               else if (strcmp(name,"Open")==0) return OpenEnum;
               else if (strcmp(name,"Option")==0) return OptionEnum;
@@ -1470 +1490 @@
               else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
               else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
-              else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
+         else stage=13;
+   }
+   if(stage==13){
+              if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
               else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
               else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;
@@ -1490 +1513 @@
               else if (strcmp(name,"VertexSId")==0) return VertexSIdEnum;
               else if (strcmp(name,"Vertices")==0) return VerticesEnum;
-         else stage=13;
-   }
-   if(stage==13){
-              if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
+              else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
               else if (strcmp(name,"Water")==0) return WaterEnum;
               else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;

issm/trunk-jpl/src/c/shared/Numerics/types.h

@@ -13 +13 @@
 
 #include <stdio.h>
+#include <complex>
 
 /*here are our abstracted types: inspired on petsc */
@@ -24 +25 @@
 /*ADOLC typedefs*/
 #include "adolc/adolc.h"
-typedef adouble IssmDouble; /*for active variables*/
-typedef double  IssmPDouble; /*for passive variables*/
+typedef adouble  IssmDouble; /*for active variables*/
+typedef acomplex IssmComplex; /*for active variables*/
+typedef double   IssmPDouble; /*for passive variables*/
+typedef std::complex<double> IssmPComplex; /*for passive variables*/
 
 #elif defined(_HAVE_CODIPACK_) && !defined(_WRAPPERS_)
 /*CoDiPack typedefs*/
 #include <codi.hpp>
-typedef codi::RealReverse IssmDouble;
-typedef double            IssmPDouble;
+typedef codi::RealReverse    IssmDouble;
+//typedef codi::ComplexReverse IssmComplex; //not defined so far
+typedef double               IssmPDouble;
+typedef std::complex<double>  IssmPComplex;
 
 /*Non-AD typedefs*/
 #else 
-typedef double     IssmDouble; 
-typedef IssmDouble IssmPDouble;
+typedef double              IssmDouble; 
+typedef std::complex<double> IssmComplex; 
+typedef IssmDouble          IssmPDouble;
+typedef IssmComplex         IssmPComplex;
 #endif
 

issm/trunk-jpl/src/m/classes/fourierlove.m

@@ -19 +19 @@
                 istemporal                 = 0;
                 n_temporal_iterations      = 0;
-                time                       = 0;
+                time                               = 0;
                 love_kernels               = 0;
                 forcing_type               = 0;
-                inner_core_boundary        = 0;
+                inner_core_boundary            = 0;
                 core_mantle_boundary       = 0;
+                complex_computation        = 0;
+
         end
         methods (Static)
@@ -63 +65 @@
                         self.inner_core_boundary=1;
                         self.core_mantle_boundary=2;
+                        self.complex_computation=0;
                 end % }}}
                 function disp(self) % {{{
@@ -102 +105 @@
                         md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]);
                         md = checkfield(md,'fieldname','love.forcing_type','NaN',1,'Inf',1,'numel',1,'>',0, '<=', 12);
+                        md = checkfield(md,'fieldname','love.complex_computation','NaN',1,'Inf',1,'numel',1,'values',[0 1]);
 
                         md = checkfield(md,'fieldname','love.istemporal','values',[0 1]);
@@ -140 +144 @@
                         WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean');
                         WriteData(fid,prefix,'object',self,'fieldname','n_temporal_iterations','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','complex_computation','format','Boolean');
                         %note: no need to marshall the time vector, we have frequencies
                         WriteData(fid,prefix,'object',self,'fieldname','love_kernels','format','Boolean');

issm/trunk-jpl/src/m/solve/loadresultfromdisk.m

@@ -49 +49 @@
                                         N=fread(fid,1,'int');
                                         field=fread(fid,[N M],'int')';
+                                elseif type==5,
+                                        N=fread(fid,1,'int');
+                                        fieldr=fread(fid,[N M],'double')';
+                                        fieldi=fread(fid,[N M],'double')';
+                                        field=complex(fieldr,fieldi);
                                 else
                                         error(['cannot read data of type ' num2str(type) ]);

issm/trunk-jpl/test/NightlyRun/test2084.m

@@ -26 +26 @@
 md.materials.lame_lambda=md.materials.lame_mu*0+5e17;
 md.materials.issolid=[1 0 1 1 1 1]';
-md.materials.rheologymodel=2*ones(md.materials.numlayers,1);
+md.materials.rheologymodel=0*ones(md.materials.numlayers,1);
 md.materials.burgers_mu=md.materials.lame_mu/3;
 md.materials.burgers_viscosity=md.materials.viscosity/10;