Changeset 27308


Ignore:
Timestamp:
10/13/22 15:01:31 (2 years ago)
Author:
caronlam
Message:

CHG: sea-level change core optimization (runtime and memory use). NEW: love number core, added support for parallelization, quad precision and EBM rheology. CHG: class fourierlove renamed to just love

Location:
issm/trunk-jpl
Files:
2 added
28 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/HydrologyTwsAnalysis.cpp

    r26800 r27308  
    3535
    3636        /*Plug inputs into element:*/
    37         iomodel->FetchDataToInput(inputs,elements,"md.hydrology.spcwatercolum", HydrologyTwsSpcEnum);
     37        iomodel->FetchDataToInput(inputs,elements,"md.hydrology.spcwatercolumn", HydrologyTwsSpcEnum);
    3838
    3939        /*Initialize sea level cumulated sea level loads :*/
    … …  
    5353        /*Now, do we really want Tws?*/
    5454        if(hydrology_model!=HydrologyTwsEnum) return;
    55 
    5655        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
    5756

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

    r26800 r27308  
    2020void LoveAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    2121        IssmDouble* frequencies = NULL;
    22         int         nfreq;
     22        IssmDouble* hypergeom_z = NULL;
     23        IssmDouble* hypergeom_table1 = NULL;
     24        IssmDouble* hypergeom_table2 = NULL;
     25        int         nfreq,nz, nalpha;
    2326        iomodel->FetchData(&nfreq,"md.love.nfreq");
     27        iomodel->FetchData(&nz,"md.love.hypergeom_nz");
     28        iomodel->FetchData(&nalpha,"md.love.hypergeom_nalpha");
    2429        iomodel->FetchData(&frequencies,NULL,NULL,"md.love.frequencies");
     30        iomodel->FetchData(&hypergeom_z,NULL,NULL,"md.love.hypergeom_z");
     31        iomodel->FetchData(&hypergeom_table1,NULL,NULL,"md.love.hypergeom_table1");
     32        iomodel->FetchData(&hypergeom_table2,NULL,NULL,"md.love.hypergeom_table2");
     33
    2534        parameters->AddObject(new DoubleVecParam(LoveFrequenciesEnum,frequencies,nfreq));
     35        parameters->AddObject(new DoubleVecParam(LoveHypergeomZEnum,hypergeom_z,nz));
     36        parameters->AddObject(new DoubleMatParam(LoveHypergeomTable1Enum,hypergeom_table1,nz,nalpha));
     37        parameters->AddObject(new DoubleMatParam(LoveHypergeomTable2Enum,hypergeom_table2,nz,nalpha));
     38
    2639        xDelete<IssmDouble>(frequencies);
     40        xDelete<IssmDouble>(hypergeom_z);
     41        xDelete<IssmDouble>(hypergeom_table1);
     42        xDelete<IssmDouble>(hypergeom_table2);
    2743
    2844        parameters->AddObject(iomodel->CopyConstantObject("md.love.nfreq",LoveNfreqEnum));
    … …  
    3753        parameters->AddObject(iomodel->CopyConstantObject("md.love.underflow_tol",LoveUnderflowTolEnum));
    3854        parameters->AddObject(iomodel->CopyConstantObject("md.love.pw_threshold",LovePostWidderThresholdEnum));
    39         parameters->AddObject(iomodel->CopyConstantObject("md.love.integration_steps_per_layer",LoveIntStepsPerLayerEnum));
     55        parameters->AddObject(iomodel->CopyConstantObject("md.love.min_integration_steps",LoveMinIntegrationStepsEnum));
     56        parameters->AddObject(iomodel->CopyConstantObject("md.love.max_integration_dr",LoveMaxIntegrationdrEnum));
     57        parameters->AddObject(iomodel->CopyConstantObject("md.love.integration_scheme",LoveIntegrationSchemeEnum));
    4058        parameters->AddObject(iomodel->CopyConstantObject("md.love.istemporal",LoveIsTemporalEnum));
    4159        parameters->AddObject(iomodel->CopyConstantObject("md.love.n_temporal_iterations",LoveNTemporalIterationsEnum));
    … …  
    4563        parameters->AddObject(iomodel->CopyConstantObject("md.love.core_mantle_boundary",LoveCoreMantleBoundaryEnum));
    4664        parameters->AddObject(iomodel->CopyConstantObject("md.love.complex_computation",LoveComplexComputationEnum));
     65        parameters->AddObject(iomodel->CopyConstantObject("md.love.quad_precision",LoveQuadPrecisionEnum));
     66        parameters->AddObject(iomodel->CopyConstantObject("md.love.debug",LoveDebugEnum));
    4767        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
    4868        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
    4969        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
    5070        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     71        parameters->AddObject(iomodel->CopyConstantObject("md.love.hypergeom_nz",LoveHypergeomNZEnum));
     72        parameters->AddObject(iomodel->CopyConstantObject("md.love.hypergeom_nalpha",LoveHypergeomNAlphaEnum));
    5173}/*}}}*/
    5274void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

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

    r27131 r27308  
    4949        iomodel->FetchDataToInput(inputs,elements,"md.geometry.bed",BedEnum);
    5050
     51        iomodel->FetchDataToInput(inputs,elements,"md.solidearth.transfercount",CouplingTransferCountEnum);
     52
    5153        /*external solidearthsolution: solid-Earth model*/
    5254        iomodel->FetchData(&isexternal,"md.solidearth.isexternal");
    … …  
    5961
    6062                /*Resolve Mmes using the modelid, if necessary:*/
    61                 if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementEastRateEnum)==DatasetInputEnum){
     63                if (inputs->GetInputObjectEnum(SolidearthExternalDisplacementUpRateEnum)==DatasetInputEnum){
    6264                        int modelid;
    6365

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r27287 r27308  
    406406                virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
    407407                virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
    408                 virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0;
     408                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount)=0;
    409409                virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    410410                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r27131 r27308  
    227227                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    228228                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    229                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
     229                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");};
    230230                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    231231                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r27131 r27308  
    180180                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    181181                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    182                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
     182                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");};
    183183                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    184184                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.h

    r27131 r27308  
    187187                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    188188                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    189                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
     189                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount){_error_("not implemented yet");};
    190190                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    191191                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

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

    r27296 r27308  
    21982198        //compute sea level load weights
    21992199        this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset);
     2200
     2201        //failsafe for phi so small that GetFractionGeometry returns 0 
     2202        if (phi==0) phi=1e-16;
     2203
    22002204        for (int i=0;i<NUMVERTICES;i++) loadweights[i]/=phi;
    22012205        this->GetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius);
    … …  
    64056409}
    64066410/*}}}*/
    6407 void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/
     6411void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* n_activevertices){ /*{{{*/
    64086412
    64096413        /*Declarations:{{{*/
    … …  
    65556559        }
    65566560
    6557         AlphaIndex=xNewZeroInit<int>(3*nel);
    6558         if (horiz) AzimuthIndex=xNewZeroInit<int>(3*nel);
     6561        AlphaIndex=xNewZeroInit<int>(n_activevertices[this->lid]*nel);
     6562        if (horiz) AzimuthIndex=xNewZeroInit<int>(n_activevertices[this->lid]*nel);
    65596563        int intmax=pow(2,16)-1;
    65606564
     6565        int* activevertices=xNew<int>(n_activevertices[this->lid]);
     6566       
     6567        int av=0;
    65616568
    65626569        for (int i=0;i<3;i++){
    65636570                if(lids[this->vertices[i]->lid]==this->lid){
     6571                        activevertices[av]=i;
    65646572                        for(int e=0;e<nel;e++){
    65656573                                IssmDouble alpha;
    … …  
    65856593                                        dy=sin(longe-longi)*cos(late);
    65866594                                        //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
    6587                                         AzimuthIndex[i*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6595                                        AzimuthIndex[av*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    65886596                                }
    6589                                 AlphaIndex[i*nel+e]=index;
    6590                         }                       
     6597                                AlphaIndex[av*nel+e]=index;
     6598                        }
     6599                        av+=1;                 
    65916600                } //for (int i=0;i<3;i++)
    65926601        } //for(int e=0;e<nel;e++)
    65936602
    65946603        /*Add in inputs:*/
    6595         this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3);
    6596         if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3);
     6604        this->inputs->SetIntArrayInput(SealevelchangeConvolutionVerticesEnum,this->lid,activevertices,n_activevertices[this->lid]);
     6605        this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*n_activevertices[this->lid]);
     6606        if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*n_activevertices[this->lid]);
    65976607
    65986608        /*}}}*/
    … …  
    67136723        /*Free allocations:{{{*/
    67146724        #ifdef _HAVE_RESTRICT_
     6725        delete activevertices;
    67156726        delete AlphaIndex;
    67166727        if(horiz) AzimuthIndex;
    … …  
    67266737
    67276738        #else
     6739        xDelete<int>(activevertices);
    67286740        xDelete<int>(AlphaIndex);
    67296741        if(horiz){
    … …  
    67576769        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    67586770        IssmDouble xyz_list[NUMVERTICES][3];
     6771        int* activevertices = NULL;
     6772        int n_activevertices, av;
    67596773
    67606774        #ifdef _HAVE_RESTRICT_
    … …  
    68316845        if(horiz) AzimIndex=xNew<int*>(SLGEOM_NUMLOADS);
    68326846
     6847        this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
     6848        // 0<=n_activevertices<=3 is the number of vertices this element is in charge of computing fields in during the sea level convolutions
     6849        // activevertices contains the vertex indices (1,2 and/or 3) in case debugging is required, they are supposed to appear in the same order as slgeom->lids
    68336850
    68346851        //Allocate:
    68356852        for(int l=0;l<SLGEOM_NUMLOADS;l++){
    68366853                int nbar=slgeom->nbar[l];
    6837                 AlphaIndex[l]=xNewZeroInit<int>(3*nbar);
    6838                 if(horiz) AzimIndex[l]=xNewZeroInit<int>(3*nbar);
    6839 
    6840 
    6841                 for (int i=0;i<3;i++){
    6842                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    6843                                 for(int e=0;e<nbar;e++){
    6844                                         IssmDouble alpha;
    6845                                         IssmDouble delPhi,delLambda;
    6846                                         /*recover info for this element and vertex:*/
    6847                                         IssmDouble late= slgeom->latbarycentre[l][e];
    6848                                         IssmDouble longe= slgeom->longbarycentre[l][e];
    6849                                         late=late/180*M_PI;
    6850                                         longe=longe/180*M_PI;
    6851 
    6852                                         lati=latitude[i];
    6853                                         longi=longitude[i];
    6854 
     6854                AlphaIndex[l]=xNewZeroInit<int>(n_activevertices*nbar);
     6855                if(horiz) AzimIndex[l]=xNewZeroInit<int>(n_activevertices*nbar);
     6856
     6857                //av=0;
     6858                //for (int i=0;i<3;i++){
     6859                for (int av=0;av<n_activevertices;av++){
     6860                        //if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     6861                        int i=activevertices[av];
     6862                        for(int e=0;e<nbar;e++){
     6863                                IssmDouble alpha;
     6864                                IssmDouble delPhi,delLambda;
     6865                                /*recover info for this element and vertex:*/
     6866                                IssmDouble late= slgeom->latbarycentre[l][e];
     6867                                IssmDouble longe= slgeom->longbarycentre[l][e];
     6868                                late=late/180*M_PI;
     6869                                longe=longe/180*M_PI;
     6870                                lati=latitude[i];
     6871                                longi=longitude[i];
    68556872                                        if(horiz){
    6856                                                 /*Compute azimuths*/
     6873                                        /*Compute azimuths*/
    68576874                                                dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
    68586875                                                dy=sin(longe-longi)*cos(late);
    68596876                                                //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
    6860                                                 AzimIndex[l][i*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6877                                                AzimIndex[l][av*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    68616878                                        }
    68626879
    6863                                         /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    6864                                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    6865                                         alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
    6866                                         doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    6867                                         index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    6868 
    6869                                         if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
    6870                                         if (index==M-1){ //avoids out of bound case
    6871                                                 index-=1;
    6872                                                 lincoef=1;
    6873                                         }
    6874                                         AlphaIndex[l][i*nbar+e]=index;
     6880                                /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6881                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6882                                alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
     6883                                doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6884                                index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6885
     6886                                if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6887                                if (index==M-1){ //avoids out of bound case
     6888                                        index-=1;
     6889                                        lincoef=1;
    68756890                                }
    6876                         }
     6891                                AlphaIndex[l][av*nbar+e]=index;
     6892                        //}
     6893                        //av+=1;
     6894                        }
     6895
    68776896                }
    68786897        }
    … …  
    68806899        /*Save all these arrayins for each element:*/
    68816900        for (int l=0;l<SLGEOM_NUMLOADS;l++){
    6882                 this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3);
    6883                 if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3);
     6901                this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*n_activevertices);
     6902                if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*n_activevertices);
    68846903        }
    68856904        /*}}}*/
    … …  
    72277246
    72287247                for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i];
     7248
    72297249                #ifdef _ISSM_DEBUG_ /*{{{*/
    72307250                /*Inform mask: */
    … …  
    73377357                oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    73387358        }
     7359
     7360        oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     7361        oceanaverage*=rho_water*oceanarea;
     7362
     7363        /*add ocean average in the global sealevelloads vector:*/
     7364        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7365                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7366                loads->vsubsealevelloads->SetValue(intj,oceanaverage,INS_VAL);
     7367                loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
     7368        }
     7369        else loads->vsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
     7370
    73397371        #ifdef _ISSM_DEBUG_
    73407372        this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    73417373        #endif
    7342         oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
    7343 
    7344         /*add ocean average in the global sealevelloads vector:*/
    7345         if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    7346                 int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7347                 loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL);
    7348                 loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
    7349         }
    7350         else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL);
    73517374
    73527375        /*add ocean area into a global oceanareas vector:*/
    … …  
    73677390        IssmDouble* G=NULL;
    73687391        IssmDouble* Grot=NULL;
     7392        IssmDouble* rslfield=NULL;
    73697393        DoubleVecParam* parameter;
    73707394        bool computefuture=false;
    7371         int spatial_component=0;
    73727395
    73737396        bool sal = false;
    … …  
    73887411                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    73897412
    7390                 this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
    7391                 for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    73927413                if (rotation)   this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
    73937414
    7394                 this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
     7415                rslfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=false);
     7416                this->SealevelchangeCollectGrdfield(sealevelpercpu,rslfield,slgeom,nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
     7417
    73957418        }
    73967419
    … …  
    74067429        int nel,nbar;
    74077430        bool sal = false;
    7408         int* AlphaIndex=NULL;
    7409         int* AzimIndex=NULL;
    7410         int* AlphaIndexsub[SLGEOM_NUMLOADS];
    7411         int* AzimIndexsub[SLGEOM_NUMLOADS];
    74127431        int spatial_component=0;
    74137432        IssmDouble* G=NULL;
    … …  
    74187437        IssmDouble* GNrot=NULL;
    74197438        IssmDouble* GErot=NULL;
     7439        IssmDouble* grdfield=NULL;
    74207440
    74217441        DoubleVecParam* parameter;
    … …  
    74387458
    74397459        if(sal){
    7440                 this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
    7441                 for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    7442 
    74437460                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    74447461                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    … …  
    74497466
    74507467                        if(horiz){
    7451                                 this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
    7452                                 for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
    7453 
    74547468                                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    74557469                                parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL);
    … …  
    74647478                        }
    74657479                }
    7466 
    7467                 this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
     7480                //Relative sea level convolution
     7481                grdfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=true);
     7482                this->SealevelchangeCollectGrdfield(&RSLGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
    74687483
    74697484                if(elastic){
    7470                         this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
     7485                        //Bedrock Uplift
     7486                        grdfield=this->SealevelchangeGxL(GU,GUrot,loads,polarmotionvector,slgeom,nel,computefuture=true);
     7487                        this->SealevelchangeCollectGrdfield(&UGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
     7488
    74717489                        if(horiz){
    7472                                 this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
    7473                                 this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
     7490                                //Bedrock North displacement
     7491                                grdfield=this->SealevelchangeHorizGxL(spatial_component=1,GH,GNrot,loads,polarmotionvector,slgeom,nel,computefuture=true);
     7492                                this->SealevelchangeCollectGrdfield(&NGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
     7493
     7494                                //Bedrock East displacement
     7495                                grdfield=this->SealevelchangeHorizGxL(spatial_component=2,GH,GErot,loads,polarmotionvector,slgeom,nel,computefuture=true);
     7496                                this->SealevelchangeCollectGrdfield(&EGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    74747497                        }
    74757498                }
    … …  
    74967519
    74977520} /*}}}*/
    7498 void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
     7521IssmDouble*       Tria::SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/
    74997522
    75007523        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    7501 
     7524        int* AlphaIndex=NULL;
     7525        int* AlphaIndexsub[SLGEOM_NUMLOADS];
     7526        int* activevertices=NULL;
    75027527        IssmDouble* grdfield=NULL;
    7503         int i,e,l,t,a, index, nbar;
     7528        int i,e,l,t,a, index, nbar, size, av,ae,b,c;
    75047529        bool rotation=false;
    7505         IssmDouble* Centroid_loads=NULL;
    7506         IssmDouble* Centroid_loads_copy=NULL;
    7507         IssmDouble* Subelement_loads[SLGEOM_NUMLOADS];
    7508         IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS];
    7509         IssmDouble* horiz_projection=NULL;
    7510         IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
    75117530        int nt=1; //important, ensures there is a defined value if computeviscous is false
     7531        int n_activevertices=0;
    75127532
    75137533        //viscous
    … …  
    75157535        int viscousindex=0; //important
    75167536        int viscousnumsteps=1; //important
     7537
     7538        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
     7539        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     7540
     7541        //Get green functions indexing & geometry
     7542        this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); //the order in which the vertices appear here should be the same as in slgeom->lids
     7543        this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7544        for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
     7545
     7546        if(computeviscous){
     7547                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
     7548                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     7549                if(computefuture) {
     7550                        nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity
     7551                        if (nt>viscousnumsteps) nt=viscousnumsteps;
     7552                }
     7553                else nt=1;
     7554        }
     7555        //allocate
     7556        grdfield=xNewZeroInit<IssmDouble>(3*nt);
     7557
     7558        //early return
     7559        if (n_activevertices==0) return grdfield;
     7560
     7561        if(rotation){ //add rotational feedback
     7562                for(av=0;av<n_activevertices;av++) { //vertices
     7563                        i=activevertices[av];
     7564                        //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7565                        b=i*nt;
     7566                        for (int m=0;m<3;m++){ //polar motion components
     7567                                for(t=0;t<nt;t++){ //time
     7568                                        int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
     7569                                        grdfield[b+t]+=Grot[index]*polarmotionvector[m];
     7570                                }
     7571                        }
     7572                }
     7573        }
     7574
     7575        //Convolution
     7576        for(av=0;av<n_activevertices;av++) { /*{{{*/
     7577                i=activevertices[av];
     7578                //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7579                b=i*nt;
     7580                c=av*nel;
     7581                for (int ae=0;ae<loads->nactiveloads;ae++){
     7582                        e=loads->combined_loads_index[ae];
     7583                        a=AlphaIndex[c+e]*viscousnumsteps;
     7584                        for(t=0;t<nt;t++){
     7585                                grdfield[b+t]+=G[a+t]*loads->combined_loads[ae];
     7586                        }
     7587                }
     7588                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7589                        nbar=slgeom->nbar[l];
     7590                        c=av*nbar;
     7591                        for (ae=0;ae<loads->nactivesubloads[l];ae++){
     7592                                e=loads->combined_subloads_index[l][ae];
     7593                                a=AlphaIndexsub[l][c+e]*viscousnumsteps;
     7594                                for(t=0;t<nt;t++){
     7595                                        grdfield[b+t]+=G[a+t]*loads->combined_subloads[l][ae];
     7596                                }
     7597                        }
     7598                }
     7599                //av+=1;
     7600        } /*}}}*/
     7601
     7602        return grdfield;
     7603
     7604} /*}}}*/
     7605IssmDouble*       Tria::SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/
     7606
     7607        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
     7608        int* AlphaIndex=NULL;
     7609        int* AzimIndex=NULL;
     7610        int* AlphaIndexsub[SLGEOM_NUMLOADS];
     7611        int* AzimIndexsub[SLGEOM_NUMLOADS];
     7612        int* activevertices = NULL;
     7613        IssmDouble* grdfield=NULL;
     7614        int i,e,l,t,a,b,c, index, nbar, av, ae,n_activevertices, size;
     7615        bool rotation=false;
     7616        IssmDouble* projected_loads=NULL;
     7617        IssmDouble* projected_subloads[SLGEOM_NUMLOADS];
     7618        IssmDouble* horiz_projection=NULL;
     7619        IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
     7620        int nt=1; //important, ensures there is a defined value if computeviscous is false
     7621
     7622        //viscous
     7623        bool computeviscous=false;
     7624        int viscousindex=0; //important
     7625        int viscousnumsteps=1; //important
     7626
     7627        //Get green functions indexing & geometry
     7628        this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
     7629        this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7630        for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
     7631        this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
     7632        for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
     7633
     7634        //First, figure out how many time steps to compute grdfield for
     7635        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
     7636        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     7637        if(computeviscous){
     7638                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
     7639                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     7640                if(computefuture) {
     7641                        nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity
     7642                        if (nt>viscousnumsteps) nt=viscousnumsteps;
     7643                }
     7644                else nt=1;
     7645        }
     7646        //allocate
     7647        grdfield=xNewZeroInit<IssmDouble>(3*nt);
     7648        if (n_activevertices==0) return grdfield;
     7649
     7650        if(rotation){ //add rotational feedback
     7651                for(av=0;av<n_activevertices;av++) { //vertices
     7652                        i=activevertices[av];
     7653                        //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7654                        for (int m=0;m<3;m++){ //polar motion components
     7655                                for(t=0;t<nt;t++){ //time
     7656                                        int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
     7657                                        grdfield[i*nt+t]+=Grot[index]*polarmotionvector[m];
     7658                                }
     7659                        }
     7660                        //}
     7661                }
     7662        }
     7663
     7664        //Initialize projection vectors
     7665        horiz_projection=xNewZeroInit<IssmDouble>(loads->nactiveloads);
     7666        projected_loads=xNewZeroInit<IssmDouble>(loads->nactiveloads);
     7667        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7668                //nbar=slgeom->nbar[l];
     7669                projected_subloads[l]=xNewZeroInit<IssmDouble>(loads->nactivesubloads[l]);
     7670                horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(loads->nactivesubloads[l]);
     7671        }
     7672
     7673
     7674        //Convolution
     7675        //av=0;
     7676        for(av=0;av<n_activevertices;av++) { //vertices
     7677                i=activevertices[av];
     7678                //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7679                b=i*nt;
     7680
     7681                //GxL needs to be projected on the right axis before summation into the grd field
     7682                //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
     7683
     7684                //get projection
     7685                if (spatial_component==1){ //north
     7686                        for (int ae=0;ae<loads->nactiveloads;ae++){
     7687                                e=loads->combined_loads_index[ae];
     7688                                horiz_projection[ae]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[av*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
     7689                        }
     7690                        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7691                                nbar=slgeom->nbar[l];
     7692                                for (ae=0;ae<loads->nactivesubloads[l];ae++){
     7693                                        e=loads->combined_subloads_index[l][ae];
     7694                                        horiz_projectionsub[l][ae]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][av*nbar+e])/65535.0);
     7695                                }
     7696                        }
     7697                }
     7698                else if (spatial_component==2){ //east
     7699                        for (int ae=0;ae<loads->nactiveloads;ae++){
     7700                                e=loads->combined_loads_index[ae];
     7701                                horiz_projection[ae]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[av*nel+e])/65535.0);
     7702                        }
     7703                        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7704                                nbar=slgeom->nbar[l];
     7705                                for (ae=0;ae<loads->nactivesubloads[l];ae++){
     7706                                        e=loads->combined_subloads_index[l][ae];
     7707                                        horiz_projectionsub[l][ae]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][av*nbar+e])/65535.0);
     7708                                }
     7709                        }
     7710                }
     7711
     7712                //project load in the right direction
     7713                for (int ae=0;ae<loads->nactiveloads;ae++){
     7714                        projected_loads[ae]=loads->combined_loads[ae]*horiz_projection[ae];
     7715                }
     7716                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7717                        nbar=slgeom->nbar[l];
     7718                        for (ae=0;ae<loads->nactivesubloads[l];ae++){
     7719                                projected_subloads[l][ae]=loads->combined_subloads[l][ae]*horiz_projectionsub[l][ae];
     7720                        }
     7721                }
     7722
     7723                //do the convolution
     7724                c=av*nel;
     7725                for (int ae=0;ae<loads->nactiveloads;ae++){
     7726                        e=loads->combined_loads_index[ae];
     7727                        a=AlphaIndex[c+e]*viscousnumsteps;
     7728                        for(t=0;t<nt;t++){
     7729                                grdfield[b+t]+=G[a+t]*projected_loads[ae];
     7730                        }
     7731                }
     7732                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7733                        nbar=slgeom->nbar[l];
     7734                        c=av*nbar;
     7735                        for (ae=0;ae<loads->nactivesubloads[l];ae++){
     7736                                e=loads->combined_subloads_index[l][ae];
     7737                                a=AlphaIndexsub[l][c+e]*viscousnumsteps;
     7738                                for(t=0;t<nt;t++){
     7739                                        grdfield[b+t]+=G[a+t]*projected_subloads[l][ae];
     7740                                }
     7741                        }
     7742                }
     7743                //av+=1;
     7744        } /*}}}*/
     7745
     7746
     7747        //free resources
     7748
     7749        xDelete<IssmDouble>(horiz_projection);
     7750        xDelete<IssmDouble>(projected_loads);
     7751        for(l=0;l<SLGEOM_NUMLOADS;l++) {
     7752                xDelete<IssmDouble>(projected_subloads[l]);
     7753                xDelete<IssmDouble>(horiz_projectionsub[l]);
     7754        }
     7755        return grdfield;
     7756
     7757} /*}}}*/
     7758
     7759
     7760void       Tria::SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
     7761
     7762        //This function aligns grdfield with the requested output format: in a size 3 vector or in a size numberofvertices vector
     7763        // if compute viscous is on, we also interpolate the field timewise given the current timestepping as well as collect viscous deformation and update the viscous deformation time series for future time steps
     7764        int i,e,l,t,a, index, nbar, av, n_activevertices;
     7765        int nt=1;
     7766
     7767
     7768        //viscous
     7769        bool computeviscous=false;
     7770        int viscousindex=0; //important
     7771        int viscousnumsteps=1; //important
     7772        int* activevertices = NULL;
    75177773        IssmDouble* viscousfield=NULL;
    75187774        IssmDouble* grdfieldinterp=NULL;
    … …  
    75227778        IssmDouble  timeacc;
    75237779
     7780        //parameters & initialization
    75247781        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    7525         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     7782        this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
     7783
    75267784        if(computeviscous){
    75277785                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    … …  
    75387796                else nt=1;
    75397797        }
    7540         //allocate
    7541         grdfield=xNewZeroInit<IssmDouble>(3*nt);
    7542 
    7543         if(rotation){ //add rotational feedback
    7544                 for(int i=0;i<NUMVERTICES;i++){ //vertices
    7545                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7546                                 for (int m=0;m<3;m++){ //polar motion components
    7547                                         for(t=0;t<nt;t++){ //time
    7548                                                 int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
    7549                                                 grdfield[i*nt+t]+=Grot[index]*polarmotionvector[m];
    7550                                         }
    7551                                 }
    7552                         }
    7553                 }
    7554         }
    7555 
    7556         //Determine loads /*{{{*/
    7557         Centroid_loads=xNewZeroInit<IssmDouble>(nel);
    7558         for (e=0;e<nel;e++){
    7559                 Centroid_loads[e]=loads->loads[e];
    7560         }
    7561         for(l=0;l<SLGEOM_NUMLOADS;l++){
    7562                 nbar=slgeom->nbar[l];
    7563                 Subelement_loads[l]=xNewZeroInit<IssmDouble>(nbar);
    7564                 for (e=0;e<nbar;e++){
    7565                         Subelement_loads[l][e]=(loads->subloads[l][e]);
    7566                 }
    7567         }
    7568         if(loads->sealevelloads){
    7569                 for (e=0;e<nel;e++){
    7570                         Centroid_loads[e]+=(loads->sealevelloads[e]);
    7571                 }
    7572                 nbar=slgeom->nbar[SLGEOM_OCEAN];
    7573                 for (e=0;e<nbar;e++){
    7574                         Subelement_loads[SLGEOM_OCEAN][e]+=(loads->subsealevelloads[e]);
    7575                 }
    7576         }
    7577 
    7578         //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex
    7579         if (spatial_component!=0){
    7580                 horiz_projection=xNewZeroInit<IssmDouble>(3*nel);
    7581                 Centroid_loads_copy=xNewZeroInit<IssmDouble>(nel);
    7582                 for (e=0;e<nel;e++){
    7583                         Centroid_loads_copy[e]=Centroid_loads[e];
    7584                 }
    7585 
    7586                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    7587                         nbar=slgeom->nbar[l];
    7588                         Subelement_loads_copy[l]=xNewZeroInit<IssmDouble>(nbar);
    7589                         horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(3*nbar);
    7590                         for (e=0;e<nbar;e++){
    7591                                 Subelement_loads_copy[l][e]=Subelement_loads[l][e];
    7592                         }
    7593                 }
    7594         }
    7595         /*}}}*/
    7596 
    7597         //Convolution
    7598         for(i=0;i<NUMVERTICES;i++) { /*{{{*/
    7599                 if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7600 
    7601                 if (spatial_component!=0){ //horizontals /*{{{*/
    7602                         //GxL needs to be projected on the right axis before summation into the grd field
    7603                         //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
    7604                         if (spatial_component==1){ //north
    7605                                 for (e=0;e<nel;e++){
    7606                                         horiz_projection[i*nel+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
    7607                                 }
    7608                                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    7609                                         nbar=slgeom->nbar[l];
    7610                                         for (e=0;e<nbar;e++){
    7611                                                 horiz_projectionsub[l][i*nbar+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
    7612                                         }
    7613                                 }
    7614                         }
    7615                         else if (spatial_component==2){ //east
    7616                                 for (e=0;e<nel;e++){
    7617                                         horiz_projection[i*nel+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0);
    7618                                 }
    7619                                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    7620                                         nbar=slgeom->nbar[l];
    7621                                         for (e=0;e<nbar;e++){
    7622                                                 horiz_projectionsub[l][i*nbar+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
    7623                                         }
    7624                                 }
    7625                         }
    7626                         for (e=0;e<nel;e++) Centroid_loads[e]=Centroid_loads_copy[e]*horiz_projection[i*nel+e];
    7627                         for(l=0;l<SLGEOM_NUMLOADS;l++){
    7628                                 nbar=slgeom->nbar[l];
    7629                                 for (e=0;e<nbar;e++){
    7630                                         Subelement_loads[l][e]=Subelement_loads_copy[l][e]*horiz_projectionsub[l][i*nbar+e];
    7631                                 }
    7632                         }
    7633                 } /*}}}*/
    7634 
    7635                 for (e=0;e<nel;e++){
    7636                         for(t=0;t<nt;t++){
    7637                                 a=AlphaIndex[i*nel+e];
    7638                                 grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Centroid_loads[e];
    7639                         }
    7640                 }
    7641                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    7642                         nbar=slgeom->nbar[l];
    7643                         for (e=0;e<nbar;e++){
    7644                                 for(t=0;t<nt;t++){
    7645                                         a=AlphaIndexsub[l][i*nbar+e];
    7646                                         grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Subelement_loads[l][e];
    7647                                 }
    7648                         }
    7649                 }
    7650         } /*}}}*/
    7651 
    7652 
    7653 
    7654         if(computeviscous){ /*{{{*/
     7798       
     7799
     7800        if(!computeviscous){ /*{{{*/
     7801                /*elastic or self attraction only case
     7802                /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
     7803                if(percpu){
     7804                        for(av=0;av<n_activevertices;av++) { //vertices
     7805                                i=activevertices[av];
     7806                                //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7807                                grdfieldout[this->vertices[i]->lid]=grdfield[i];
     7808                                //}
     7809                        }
     7810                }
     7811                else{
     7812                        for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i];
     7813                }
     7814                //free resources
     7815                xDelete<IssmDouble>(grdfield);
     7816                return;
     7817        }
     7818        else { //viscous case
    76557819                // we need to do up to 3 things (* = only if computefuture)
    7656                 // 1*: add new grdfield contribution to the viscous stack for future time steps
    7657                 // 2: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time]
     7820                // 1: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time]
     7821                // 2*: add new grdfield contribution to the viscous stack for future time steps
    76587822                // 3*: subtract from viscous stack the grdfield that has already been accounted for so we don't add it again at the next time step
     7823
     7824                /*update grdfield at present time using viscous stack at present time: */
     7825                for(av=0;av<n_activevertices;av++) { //vertices
     7826                        i=activevertices[av];
     7827                        //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7828                        grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex];
     7829                }
     7830
    76597831
    76607832                /* Map new grdfield generated by present-day loads onto viscous time vector*/
    76617833                if(computefuture){
    76627834                        //viscousindex time and first time step of grdfield coincide, so just copy that value
    7663                         for(int i=0;i<NUMVERTICES;i++){
    7664                                 if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7665                                 grdfieldinterp[i*viscousnumsteps+viscousindex]=  grdfield[i*nt+0];
     7835                        for(av=0;av<n_activevertices;av++) { //vertices
     7836                                i=activevertices[av];
     7837                                //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7838                                grdfieldinterp[i*viscousnumsteps+viscousindex]=grdfield[i*nt+0];
    76667839                        }
    76677840                        if(viscoustimes[viscousindex]<final_time){
    76687841                                //And interpolate the rest of the points in the future
    76697842                                lincoeff=(viscoustimes[viscousindex+1]-viscoustimes[viscousindex])/timeacc;
    7670                                 for(int t=viscousindex+1;t<viscousnumsteps;t++){
    7671                                         for(int i=0;i<NUMVERTICES;i++){
    7672                                                 if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7673                                                 grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]
    7674                                                                                          +lincoeff*grdfield[i*nt+(t-viscousindex)];
     7843                                for(av=0;av<n_activevertices;av++) { //vertices
     7844                                        i=activevertices[av];
     7845                                        //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7846                                        int i_time1= i*nt-viscousindex;
     7847                                        int i_time2= i*viscousnumsteps;
     7848                                        for(int t=viscousindex+1;t<viscousnumsteps;t++){
     7849                                                grdfieldinterp[i_time2+t] = (1-lincoeff)*grdfield[i_time1+t-1]
     7850                                                                          +    lincoeff *grdfield[i_time1+t]
     7851                                                                          +          viscousfield[i_time2+t];
     7852                                                /*update viscous stack with future deformation from present load: */
     7853                                                viscousfield[i_time2+t]=grdfieldinterp[i_time2+t]
     7854                                                                       -grdfieldinterp[i_time2+viscousindex];
    76757855                                        }
    76767856                                }
    76777857                        }
    7678                 }
    7679 
    7680                 /*update grdfield at present time using viscous stack at present time: */
    7681                 for(int i=0;i<NUMVERTICES;i++){
    7682                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7683                         grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex];
    7684                 }
    7685 
    7686                 /*update viscous stack with future deformation from present load: */
    7687                 if(computefuture){
    7688                         for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end
    7689                                 for(int i=0;i<NUMVERTICES;i++){
    7690                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7691                                         //offset viscousfield to remove all deformation that has already been added
    7692                                         viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]
    7693                                                                           -grdfieldinterp[i*viscousnumsteps+viscousindex]
    7694                                                                           -viscousfield[i*viscousnumsteps+viscousindex];
    7695                                 }
    7696                         }
    76977858                        /*Save viscous stack now that we updated the values:*/
    76987859                        this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps);
    7699 
    7700                         /*Free resources:*/
    7701                         xDelete<IssmDouble>(grdfieldinterp);
    7702                 }
    7703         }
    7704         /*}}}*/
    7705 
    7706         /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
    7707         if(percpu){
    7708                 for(i=0;i<NUMVERTICES;i++){
    7709                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     7860                }
     7861
     7862                /*store values computed on vertices*/
     7863                if(percpu){
     7864                        for(av=0;av<n_activevertices;av++) { //vertices
     7865                                i=activevertices[av];
     7866                                //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    77107867                                grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0];
    7711                         }
    7712                 }
    7713         }
    7714         else{
    7715                 for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
    7716         }
    7717         //free resources
    7718         xDelete<IssmDouble>(grdfield);
    7719         xDelete<IssmDouble>(Centroid_loads);
    7720         for(l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Subelement_loads[l]);
    7721         if (spatial_component!=0){
    7722                 xDelete<IssmDouble>(horiz_projection);
    7723                 xDelete<IssmDouble>(Centroid_loads_copy);
    7724                 for(l=0;l<SLGEOM_NUMLOADS;l++) {
    7725                         xDelete<IssmDouble>(Subelement_loads_copy[l]);
    7726                         xDelete<IssmDouble>(horiz_projectionsub[l]);
    7727                 }
    7728         }
    7729         if (computeviscous){
     7868                                //}
     7869                        }
     7870                }
     7871                else{
     7872                        for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
     7873                }
     7874                //free resources
     7875                xDelete<IssmDouble>(grdfield);
    77307876                xDelete<IssmDouble>(viscoustimes);
    77317877                if (computefuture){
    77327878                        xDelete<IssmDouble>(grdfieldinterp);
    77337879                }
    7734         }
    7735 
     7880                /*}}}*/
     7881        }
    77367882} /*}}}*/
    77377883

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r27287 r27308  
    172172                #ifdef _HAVE_SEALEVELCHANGE_
    173173                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    174                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids);
     174                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount);
    175175                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom);
    176176                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
    … …  
    245245                void           UpdateConstraintsExtrudeFromBase(void);
    246246                void           UpdateConstraintsExtrudeFromTop(void);
    247                 void           SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
     247                IssmDouble*    SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture);
     248                IssmDouble*    SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture);
     249                void           SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
    248250                /*}}}*/
    249251

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r27306 r27308  
    784784                        if(hydrology_model==HydrologypismEnum){
    785785                                analyses_temp[numanalyses++]=HydrologyPismAnalysisEnum;
     786                        }
     787                        if(hydrology_model==HydrologyTwsEnum){
     788                                analyses_temp[numanalyses++]=HydrologyTwsAnalysisEnum;
    786789                        }
    787790                        if(hydrology_model==HydrologydcEnum){

  • issm/trunk-jpl/src/c/classes/GrdLoads.cpp

    r27131 r27308  
    2020
    2121        /*allocate:*/
     22        nactiveloads=0;
     23
    2224        vloads=new Vector<IssmDouble>(nel);
    23         for (int i=0;i<SLGEOM_NUMLOADS;i++) vsubloads[i]=new Vector<IssmDouble>(slgeom->nbar[i]);
     25        for (int i=0;i<SLGEOM_NUMLOADS;i++) {
     26                vsubloads[i]=new Vector<IssmDouble>(slgeom->nbar[i]);
     27                combined_subloads[i]=xNewZeroInit<IssmDouble>(slgeom->nbar[i]);
     28        }
    2429
    2530        vsealevelloads=new Vector<IssmDouble>(nel);
    … …  
    2732
    2833        vsubsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
     34
     35        combined_loads=NULL;
     36        combined_loads_index=NULL;
     37        for (int i=0;i<SLGEOM_NUMLOADS;i++) {
     38                nactivesubloads[i]=0;
     39                combined_subloads[i]=NULL;
     40                combined_subloads_index[i]=NULL;
     41        }
    2942
    3043        /*make sure all pointers that are not allocated are NULL:*/
    … …  
    4053        delete vloads;
    4154        xDelete<IssmDouble>(loads);
    42         for(int i=0;i<SLGEOM_NUMLOADS;i++){
    43                 delete vsubloads[i];
    44                 xDelete<IssmDouble>(subloads[i]);
    45         }
    4655        delete vsealevelloads;
    4756        xDelete<IssmDouble>(sealevelloads);
    4857        delete vsubsealevelloads;
    4958        xDelete<IssmDouble>(subsealevelloads);
     59        if (combined_loads) xDelete<IssmDouble>(combined_loads);
     60        if (combined_loads_index) xDelete<int>(combined_loads_index);
     61        for(int i=0;i<SLGEOM_NUMLOADS;i++){
     62                delete vsubloads[i];
     63                xDelete<IssmDouble>(subloads[i]);
     64                if (combined_subloads[i]) xDelete<IssmDouble>(combined_subloads[i]);
     65                if (combined_subloads_index[i]) xDelete<int>(combined_subloads_index[i]);
     66        }
     67
     68
    5069}; /*}}}*/
    5170
    … …  
    139158        ISSM_MPI_Bcast(&deg2coeff[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    140159
    141 } /*}}}*/
     160}; /*}}}*/
     161void GrdLoads::Combineloads(int nel,SealevelGeometry* slgeom){ /*{{{*/
     162
     163        int e,l, nbar, ae;
     164        //Determine loads /*{{{*/
     165        nactiveloads=0;
     166
     167        if (combined_loads) xDelete<IssmDouble>(combined_loads);
     168        if (combined_loads_index) xDelete<int>(combined_loads_index);
     169
     170        //find non zero centroid loads, combine with sealevelloads
     171        if(sealevelloads){
     172                for (e=0;e<nel;e++){
     173                        if (loads[e]+sealevelloads[e]!=0) nactiveloads++;
     174                }
     175        }
     176        else {
     177                for (e=0;e<nel;e++){
     178                        if (loads[e]!=0) nactiveloads++;
     179                }
     180        }
     181
     182        combined_loads=xNewZeroInit<IssmDouble>(nactiveloads);
     183        combined_loads_index=xNewZeroInit<int>(nactiveloads);
     184
     185        ae=0;
     186        if(sealevelloads){
     187                for (e=0;e<nel;e++){
     188                        if (loads[e]+sealevelloads[e]!=0){
     189                                combined_loads[ae]=loads[e]+sealevelloads[e];
     190                                combined_loads_index[ae]=e;
     191                                ae++;
     192                        }
     193                }
     194        }
     195        else {
     196                for (e=0;e<nel;e++){
     197                        if (loads[e]!=0){
     198                                combined_loads[ae]=loads[e];
     199                                combined_loads_index[ae]=e;
     200                                ae++;                   
     201                        }
     202                }
     203        }
     204
     205
     206        //subloads
     207        for(l=0;l<SLGEOM_NUMLOADS;l++){
     208                nactivesubloads[l]=0;
     209                nbar=slgeom->nbar[l];
     210                if (combined_subloads[l]) xDelete<IssmDouble>(combined_subloads[l]);
     211                if (combined_subloads_index[l]) xDelete<int>(combined_subloads_index[l]);
     212
     213                //find non zero subelement loads, combine with subsealevelloads
     214                if(subsealevelloads && l==SLGEOM_OCEAN){
     215                        for (e=0;e<nbar;e++){
     216                                if (subloads[l][e]+subsealevelloads[e]!=0) nactivesubloads[l]++;
     217                        }
     218                }
     219                else {
     220                        for (e=0;e<nbar;e++){
     221                                if (subloads[l][e]!=0) nactivesubloads[l]++;;
     222                        }
     223                }
     224
     225                combined_subloads[l]=xNewZeroInit<IssmDouble>(nactivesubloads[l]);
     226                combined_subloads_index[l]=xNewZeroInit<int>(nactivesubloads[l]);
     227
     228                ae=0;
     229                if(subsealevelloads && l==SLGEOM_OCEAN){
     230                        for (e=0;e<nbar;e++){
     231                                if (subloads[l][e]+sealevelloads[e]!=0){
     232                                        combined_subloads[l][ae]=subloads[l][e]+sealevelloads[e];
     233                                        combined_subloads_index[l][ae]=e;
     234                                        ae++;
     235                                }
     236                        }
     237                }
     238                else {
     239                        for (e=0;e<nbar;e++){
     240                                if (subloads[l][e]!=0){
     241                                        combined_subloads[l][ae]=subloads[l][e];
     242                                        combined_subloads_index[l][ae]=e;
     243                                        ae++;                   
     244                                }
     245                        }
     246                }
     247
     248
     249                /*for(l=0;l<SLGEOM_NUMLOADS;l++){
     250                        nbar=slgeom->nbar[l];
     251                        for (e=0;e<nbar;e++){
     252                                combined_subloads[l][e]=(subloads[l][e]);
     253                        }
     254                }
     255                if(sealevelloads){
     256                        nbar=slgeom->nbar[SLGEOM_OCEAN];
     257                        for (e=0;e<nbar;e++){
     258                                combined_subloads[SLGEOM_OCEAN][e]+=(subsealevelloads[e]);
     259                        }
     260                }*/
     261        }
     262}; /*}}}*/

  • issm/trunk-jpl/src/c/classes/GrdLoads.h

    r26800 r27308  
    1414        public:
    1515
     16                int nactiveloads=0;
     17                int nactivesubloads[SLGEOM_NUMLOADS];
    1618                Vector<IssmDouble>* vloads=NULL;
    1719                IssmDouble*         loads=NULL;
    … …  
    2224                Vector<IssmDouble>* vsubsealevelloads=NULL;
    2325                IssmDouble*         subsealevelloads=NULL;
     26                int*                combined_loads_index=NULL;
     27                int*                combined_subloads_index[SLGEOM_NUMLOADS];
     28                IssmDouble*         combined_loads=NULL;
     29                IssmDouble*         combined_subloads[SLGEOM_NUMLOADS];
    2430
    2531                GrdLoads(int nel, SealevelGeometry* slgeom);
    … …  
    3036                void BroadcastSealevelLoads(void);
    3137                void SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom);
     38                void Combineloads(int nel, SealevelGeometry* slgeom);
    3239};
    3340#endif  /* _SEALEVELGRDLOADS_H_ */

  • issm/trunk-jpl/src/c/classes/Inputs/ArrayInput.cpp

    r27113 r27308  
    4444        ArrayInput* output = new ArrayInput(this->numberofelements_local);
    4545
    46         output->N = xNew<int>(this->numberofelements_local);
    4746        xMemCpy<int>(output->N,this->N,this->numberofelements_local);
    4847
    49         output->values = xNew<IssmDouble*>(this->numberofelements_local);
    5048        for(int i=0;i<this->numberofelements_local;i++){
    5149                if(this->values[i]){

  • issm/trunk-jpl/src/c/classes/Inputs/IntArrayInput.cpp

    r27128 r27308  
    4444        IntArrayInput* output = new IntArrayInput(this->numberofelements_local);
    4545
    46         output->N = xNew<int>(this->numberofelements_local);
    4746        xMemCpy<int>(output->N,this->N,this->numberofelements_local);
    4847
    49         output->values = xNew<int*>(this->numberofelements_local);
    5048        for(int i=0;i<this->numberofelements_local;i++){
    5149                if(this->values[i]){

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

    r26885 r27308  
    1212#ifdef _HAVE_MPLAPACK_
    1313#include <quadmath.h>
     14#include <iostream>
    1415#include "mpblas__Float128.h"
    1516#include "mplapack__Float128.h"
     17#include "mplapack_utils__Float128.h"
    1618#endif
    1719
     20template<typename doubletype> IssmDouble DownCastVarToDouble(doubletype var); // pure declaration
     21
     22template <> IssmDouble DownCastVarToDouble<IssmDouble>(IssmDouble var){
     23        return var;
     24}
     25template <> IssmDouble DownCastVarToDouble<IssmComplex>(IssmComplex var){
     26        return std::real(var);
     27}
    1828#ifdef _HAVE_MPLAPACK_
    19 _Float128 a = 0.2345234534512079875620048770134538q;
     29template <> IssmDouble DownCastVarToDouble<__float128>(__float128 var){
     30        return static_cast<IssmDouble>(var);
     31}
     32__float128 pow(__float128 x, int y){
     33        return powq(x,y);
     34}
     35__float128 pow(__float128 x, double y){
     36        return powq(x,y);
     37}
     38__float128 pow(double x, __float128 y){
     39        return powq(x,y);
     40}
     41
     42ostream& operator<<(ostream& os, __float128 x){
     43        char buf[128];
     44        quadmath_snprintf(buf, sizeof(buf), "%.34Qf", x);
     45
     46        os << buf;
     47        return os;
     48}
    2049#endif
    2150
    2251/*local definitions:*/
    23 class LoveVariables{  /*{{{*/
     52template <class doubletype> class LoveVariables{  /*{{{*/
    2453
    2554        public:
    26                 IssmDouble g0;
    27                 IssmDouble r0;
    28                 IssmDouble* EarthMass;
    29                 int nyi;
     55                doubletype g0;
     56                doubletype r0;
     57                doubletype* EarthMass;
     58                int nyi, ifreq, nfreq;
    3059                int starting_layer;
    3160                int* deg_layer_delete;
     61                int* nstep;
     62                doubletype* mu;
     63                doubletype* la;
    3264
    3365                LoveVariables(){  /*{{{*/
    … …  
    3567                        r0=0;
    3668                        EarthMass=NULL;
     69                        mu=NULL;
     70                        la=NULL;
    3771                        nyi=0;
     72                        nfreq=0;
     73                        ifreq=0;
    3874                        starting_layer=0;
     75                        nstep=NULL;
    3976                        deg_layer_delete=NULL;
    4077                } /*}}}*/
    41                 LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin, int* deg_layer_deletein){  /*{{{*/
     78                LoveVariables(doubletype* EarthMassin,doubletype g0in,doubletype r0in,int nyiin,int starting_layerin, int* deg_layer_deletein, int* nstepin){  /*{{{*/
    4279                        EarthMass=EarthMassin;
    4380                        g0=g0in;
    … …  
    4683                        starting_layer=starting_layerin;
    4784                        deg_layer_delete=deg_layer_deletein;
     85                        nstep=nstepin;
     86                        mu=NULL;
     87                        la=NULL;
     88                        nfreq=0;
     89                        ifreq=0;
    4890                } /*}}}*/
    49                 ~LoveVariables(){};
     91                ~LoveVariables(){
     92                        xDelete<int>(deg_layer_delete);
     93                        xDelete<int>(nstep);
     94                        if(mu)  xDelete<doubletype>(mu);
     95                        if(la)  xDelete<doubletype>(la);
     96                };
    5097}; /*}}}*/
    5198
    … …  
    57104                doubletype* L;
    58105                doubletype* Kernels;
    59                 int sh_nmin, sh_nmax, nfreq, nkernels;
     106                int sh_nmin, sh_nmax, nfreq, nkernels, lower_row, nfreqtotal;
    60107
    61108                LoveNumbers(){  /*{{{*/
    … …  
    69116                        nkernels=0;
    70117                } /*}}}*/
    71                 LoveNumbers(int sh_nminin, int sh_nmaxin, int nfreqin, Matlitho* matlitho){  /*{{{*/
     118                LoveNumbers(int sh_nminin, int sh_nmaxin, int nfreqin, int lower_rowin,int nfreqtotalin,Matlitho* matlitho){  /*{{{*/
    72119                        sh_nmax=sh_nmaxin;
    73120                        sh_nmin=sh_nminin;
    74121                        nfreq=nfreqin;
     122                        lower_row=lower_rowin;
     123                        nfreqtotal=nfreqtotalin;
    75124                        nkernels=(sh_nmax+1)*(matlitho->numlayers+1)*6;
    76125                        H=xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
    … …  
    81130                void DownCastToDouble(LoveNumbers<IssmDouble>* LoveDouble){
    82131                        for(int i=0;i<(sh_nmax+1)*nfreq;i++){
    83                                 LoveDouble->H[i]=std::real(H[i]);
    84                                 LoveDouble->K[i]=std::real(K[i]);
    85                                 LoveDouble->L[i]=std::real(L[i]);
    86                         }
     132                                LoveDouble->H[i]=DownCastVarToDouble<doubletype>(H[i]);
     133                                LoveDouble->K[i]=DownCastVarToDouble<doubletype>(K[i]);
     134                                LoveDouble->L[i]=DownCastVarToDouble<doubletype>(L[i]);
     135                        }
     136                        for(int i=0;i<nkernels*nfreq;i++){
     137                                LoveDouble->Kernels[i]=DownCastVarToDouble<doubletype>(Kernels[i]);
     138                        }
     139
    87140                }
    88141                void LoveMPI_Gather(LoveNumbers<doubletype>* Love_local, int lower_row){
    … …  
    112165                        xDelete<int>(displs);
    113166                }
     167                void Broadcast(void){
     168                        //Intended only for IssmDouble type
     169                        Vector<IssmDouble>* vH;
     170                        Vector<IssmDouble>* vK;
     171                        Vector<IssmDouble>* vL;
     172
     173                        vH= new Vector<IssmDouble>((sh_nmax+1)*nfreqtotal);
     174                        for(int i=0;i<(sh_nmax+1)*nfreq;i++) vH->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble<doubletype>(H[i]),INS_VAL);
     175                        xDelete(H);
     176                        vH->Assemble();
     177                        H=vH->ToMPISerial();
     178                        delete vH;
     179
     180                        vK= new Vector<IssmDouble>((sh_nmax+1)*nfreqtotal);
     181                        for(int i=0;i<(sh_nmax+1)*nfreq;i++) vK->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble<doubletype>(K[i]),INS_VAL);
     182                        xDelete(K);
     183                        vK->Assemble();
     184                        K=vK->ToMPISerial();
     185                        delete vK;
     186
     187                        vL= new Vector<IssmDouble>((sh_nmax+1)*nfreqtotal);
     188                        for(int i=0;i<(sh_nmax+1)*nfreq;i++) vL->SetValue(lower_row*(sh_nmax+1)+i,DownCastVarToDouble<doubletype>(L[i]),INS_VAL);
     189                        xDelete(L);
     190                        vL->Assemble();
     191                        L=vL->ToMPISerial();
     192                        delete vL;
     193                }
     194                void KernelBroadcast(void){
     195                        Vector<IssmDouble>* vKernels;
     196                        vKernels= new Vector<IssmDouble>(nkernels*nfreqtotal);
     197                        for(int i=0;i<nkernels*nfreq;i++){
     198                                vKernels->SetValue(lower_row*nkernels+i,DownCastVarToDouble<doubletype>(Kernels[i]),INS_VAL);
     199                        }
     200                        xDelete(Kernels);
     201                        vKernels->Assemble();
     202                        Kernels=vKernels->ToMPISerial();
     203                        delete vKernels;
     204                }
    114205                void Copy(LoveNumbers<doubletype>* LoveDup){
    115206                        for(int i=0;i<(sh_nmax+1)*nfreq;i++){
    … …  
    128219                        xDelete<doubletype>(Kernels);
    129220                };
    130 }; /*}}}*/
     221};
     222
     223 /*}}}*/
    131224
    132225/*self contained support routines used by cores below:*/
    … …  
    139232        return value;
    140233} /*}}}*/
     234#ifdef _HAVE_MPLAPACK_
     235template <> __float128                     angular_frequency<__float128>(IssmDouble frequency){ /*{{{*/
     236        return 2.0*PI*frequency;
     237} /*}}}*/
     238#endif
    141239template<typename doubletype> void                       allgesv(int* pnyi, int* pnrhs, doubletype* yilocal, int* plda, int* ipiv, doubletype* rhslocal, int* pldb, int* pinfo);
    142240template <> void                           allgesv<IssmDouble>(int* pnyi, int* pnrhs, IssmDouble* yilocal, int* plda, int* ipiv, IssmDouble* rhslocal, int* pldb, int* pinfo){ /*{{{*/
    … …  
    147245        //zgesv_(pnyi, pnrhs, yilocal, plda, ipiv, rhslocal, pldb, pinfo);
    148246} /*}}}*/
     247#ifdef _HAVE_MPLAPACK_
     248template <> void                           allgesv<__float128>(int* pnyi, int* pnrhs, __float128* yilocal, int* plda, int* ipiv, __float128* rhslocal, int* pldb, int* pinfo){ /*{{{*/
     249        mplapackint nyi=*pnyi;
     250        mplapackint nrhs=*pnrhs;
     251        mplapackint lda=*plda;
     252        mplapackint* qipiv=NULL;
     253        mplapackint ldb=*pldb;
     254        mplapackint info=0;
     255        qipiv=xNewZeroInit<mplapackint>(*pnyi);
     256       
     257        Rgesv(nyi, nrhs, yilocal, lda, qipiv, rhslocal, ldb, info);
     258
     259        for (int i;i=0;i<*pnyi) ipiv[i]=qipiv[i];
     260        *pinfo=info;
     261        xDelete<mplapackint>(qipiv);
     262} /*}}}*/
     263#endif
     264
    149265template<typename doubletype> doubletype   factorial(int n){ /*{{{*/
    150266        doubletype prod=1;
    … …  
    219335
    220336        int indxi, indf;
    221         IssmDouble PW_test;
     337        doubletype PW_test;
    222338        IssmDouble PW_threshold;
    223339        femmodel->parameters->FindParam(&PW_threshold,LovePostWidderThresholdEnum);
    … …  
    261377        Lovet[t*(sh_nmax+1)+d]=LoveM[NTit-1];
    262378}/*}}}*/
    263 template <typename doubletype> void        GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega,  Matlitho* matlitho){ /*{{{*/
     379
     380template <typename doubletype> doubletype HypergeomTableLookup(doubletype z1, doubletype alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha){/*{{{*/
     381        int iz1, iz2, ialpha;   
     382        doubletype lincoef;
     383        doubletype hf,h00,h10, h01, h11, za, zd, ha, hb,hc,hd, m0,m1,t;
     384        doubletype dalpha=1.0/(nalpha-1); // alpha table resolution given 0 <= alpha <= 1
     385        ialpha= static_cast<int>(DownCastVarToDouble(alpha/dalpha));
     386        lincoef=alpha/dalpha-ialpha;//linear fraction in [0;1] for alpha interpolation
     387        iz1=nz;
     388        for (int i=0;i<nz;i++){
     389                if (abs(z[i])>abs(z1)) {
     390                        iz1=i-1;
     391                        break;
     392                }
     393        }
     394
     395        if (iz1<0){
     396                //1-hf for very small abs(z) tends to 0, and is very log-linear with respect to log(z), so we can simply extrapolate the value of hf via the loglog slope
     397                hf=(1.0-lincoef)*h1[ialpha*nz+0]+lincoef*h1[(ialpha+1)*nz+0];
     398                hf=1.0- (1.0-hf)*pow(10.0,(log10(abs(z1))-log10(abs(z[0]))));
     399                //hf[0]=1.0;
     400        }
     401        else if (iz1==nz){
     402                //hf for very large abs(z) tends to 0, and is very log-linear with respect to log(z), so we can simply extrapolate the value of hf via the loglog slope
     403                hf=(1.0-lincoef)*h1[ialpha*nz+nz-1]+lincoef*h1[(ialpha+1)*nz+nz-1];
     404                hf=hf *pow(10.0,-(log10(abs(z1))-log10(abs(z[nz-1]))));
     405                //hf[0]=0;
     406        }
     407        else{ //cubic spline interpolation
     408                //edge cases: extrapolate 1 point
     409                if (iz1==0){
     410                        za=2.0*z[0]-z[1];
     411                        ha=(1.0-lincoef)*h1[ialpha*nz+0]+lincoef*h1[(ialpha+1)*nz+0];
     412                        ha=1.0- (1.0-ha) *pow(10.0,log10(abs(za))-log10(abs(z[0])));
     413                }
     414                else {
     415                        za=z[iz1-1];
     416                        ha=(1.0-lincoef)*h1[ialpha*nz+iz1-1] + lincoef*h1[(ialpha+1)*nz+iz1-1];
     417                }
     418
     419                if (iz1==nz-2){
     420                        zd=2.0*z[nz-1]-z[nz-2];
     421                        hd=(1.0-lincoef)*h1[ialpha*nz+nz-1]+lincoef*h1[(ialpha+1)*nz+nz-1];
     422                        hd=hd *pow(10.0,-(log10(abs(zd))-log10(abs(z[nz-1]))));
     423                }
     424                else {
     425                        zd=z[iz1+2];
     426                        hd=(1.0-lincoef)*h1[ialpha*nz+iz1+2]+lincoef*h1[(ialpha+1)*nz+iz1+2];
     427                }
     428
     429                hb=(1.0-lincoef)*h1[ialpha*nz+iz1] +lincoef*h1[(ialpha+1)*nz+iz1];
     430                hc=(1.0-lincoef)*h1[ialpha*nz+iz1+1] +lincoef*h1[(ialpha+1)*nz+iz1+1];
     431
     432                //left derivative
     433                m0= 0.5*(z[iz1+1]-z[iz1])*((hc-hb)/(z[iz1+1]-z[iz1]) + (hb-ha)/(z[iz1]-za));
     434                //right derivative
     435                m1= 0.5*(z[iz1+1]-z[iz1])*((hd-hc)/(zd-z[iz1+1]) + (hc-hb)/(z[iz1+1]-z[iz1]));
     436
     437                //interpolation abscissa
     438                t=(z1-z[iz1])/(z[iz1+1]-z[iz1]);
     439               
     440                //cubic spline functions
     441                h00=2*pow(t,3)-3*pow(t,2)+1;
     442                h10=pow(t,3)-2*pow(t,2)+t;
     443                h01=-2*pow(t,3)+3*pow(t,2);
     444                h11=pow(t,3)-pow(t,2);
     445
     446                hf=h00*hb + h10*m0 + h01*hc + h11*m1;
     447        }
     448        return hf;
     449
     450}/*}}}*/
     451
     452template <typename doubletype> doubletype muEBM(int layer_index, doubletype omega, Matlitho* matlitho, FemModel* femmodel); //pure declaration
     453template <> IssmComplex muEBM<IssmComplex>(int layer_index, IssmComplex omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/
     454        // Initialization
     455        int nz, nalpha, dummy1, dummy2;
     456        IssmComplex mu;
     457        IssmDouble* z=NULL;
     458        IssmDouble* h1=NULL;
     459        IssmDouble* h2=NULL;
     460        IssmComplex hf11, hf12, hf21, hf22;
     461        IssmDouble  factor=0;
     462        IssmComplex z1, z2;
     463        IssmComplex U1, U2;
     464        IssmComplex j=1i;
     465        //Matlitho parameters
     466        IssmDouble alpha=matlitho->ebm_alpha[layer_index];
     467        IssmDouble delta=matlitho->ebm_delta[layer_index];
     468        IssmDouble taul=matlitho->ebm_taul[layer_index];
     469        IssmDouble tauh=matlitho->ebm_tauh[layer_index];
     470        IssmDouble vi=matlitho->viscosity[layer_index];
     471        IssmDouble mu0=matlitho->lame_mu[layer_index];
     472        //fetch hypergeometric function tables and parameters
     473        femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum);
     474        femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum);
     475        femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum);
     476        femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum);
     477        femmodel->parameters->FindParam(&h2,&dummy1,&dummy2,LoveHypergeomTable2Enum);
     478        omega=omega/j;
     479
     480        z1= -pow(omega*tauh,2.0);
     481        z2= -pow(omega*taul,2.0);
     482        //Table1 h1 should be 2F1([1 1+alpha], [2+alpha/2], z)
     483        //Table2 h2 should be 2F1([1 0.5+alpha], [1.5+alpha/2], z)
     484        hf11=HypergeomTableLookup<IssmComplex>(z1, alpha, h1, z, nz, nalpha);
     485        hf21=HypergeomTableLookup<IssmComplex>(z1, alpha, h2, z, nz, nalpha);
     486        hf12=HypergeomTableLookup<IssmComplex>(z2, alpha, h1, z, nz, nalpha);
     487        hf22=HypergeomTableLookup<IssmComplex>(z2, alpha, h2, z, nz, nalpha);
     488
     489
     490        //Ivins et al (2020) p11
     491        U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf11-pow(taul,2.0+alpha)*hf12);
     492        U2=(pow(tauh,1.0+alpha)*hf21-pow(taul,1.0+alpha)*hf22)/(1.0+alpha);
     493
     494        factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha));
     495        U1=(1.0+factor) *U1;
     496        U2=factor*omega*U2 +mu0/vi/omega;
     497        mu=mu0*(U1+j*U2)/(pow(U1,2.0)+pow(U2,2.0));
     498        omega=omega*j;
     499
     500        xDelete<IssmDouble>(z);
     501        xDelete<IssmDouble>(h1);
     502        xDelete<IssmDouble>(h2);
     503        return mu;
     504}/*}}}*/
     505
     506template <> IssmDouble muEBM<IssmDouble>(int layer_index, IssmDouble omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/
     507        // Initialization
     508        int nz, nalpha, dummy1, dummy2;
     509        IssmDouble mu;
     510        IssmDouble* z=NULL;
     511        IssmDouble* h1=NULL;
     512        IssmDouble hf11, hf12;
     513        IssmDouble  factor, B, D, z1, z2;
     514        //Matlitho parameters
     515        IssmDouble alpha=matlitho->ebm_alpha[layer_index];
     516        IssmDouble delta=matlitho->ebm_delta[layer_index];
     517        IssmDouble taul=matlitho->ebm_taul[layer_index];
     518        IssmDouble tauh=matlitho->ebm_tauh[layer_index];
     519        IssmDouble vi=matlitho->viscosity[layer_index];
     520        IssmDouble mu0=matlitho->lame_mu[layer_index];
     521        //fetch hypergeometric function tables and parameters
     522        femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum);
     523        femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum);
     524        femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum);
     525        femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum);
     526
     527        z1=-omega*tauh;
     528        z2=-omega*taul;
     529        //Table1 h1 should be 2F1([1 1+alpha], [2+alpha], z)
     530        hf11=HypergeomTableLookup<IssmDouble>(z1, alpha, h1, z, nz, nalpha);
     531        hf12=HypergeomTableLookup<IssmDouble>(z2, alpha, h1, z, nz, nalpha);
     532
     533        //Ivins et al. (2022) p1979
     534        factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha));
     535        B= factor/(1.0+alpha) *mu0/vi * (pow(tauh,1.0+alpha)*hf11 - pow(taul,1.0+alpha)*hf12);
     536        D= omega*vi/mu0* 1.0/(1.0+omega*vi/mu0*(1.0+delta) -pow(omega*vi/mu0,2.0)*B);
     537
     538        xDelete<IssmDouble>(z);
     539        xDelete<IssmDouble>(h1);
     540        return mu=mu0*D;
     541}/*}}}*/
     542#ifdef _HAVE_MPLAPACK_
     543template <> __float128 muEBM<__float128>(int layer_index, __float128 omega, Matlitho* matlitho, FemModel* femmodel){/*{{{*/
     544        // Initialization
     545        int nz, nalpha, dummy1, dummy2;
     546        IssmDouble* z=NULL;
     547        IssmDouble* h1=NULL;
     548        __float128 mu;
     549        __float128 hf11, hf12;
     550        __float128  factor, B, D, z1, z2;
     551        //Matlitho parameters
     552        __float128 alpha=matlitho->ebm_alpha[layer_index];
     553        __float128 delta=matlitho->ebm_delta[layer_index];
     554        __float128 taul=matlitho->ebm_taul[layer_index];
     555        __float128 tauh=matlitho->ebm_tauh[layer_index];
     556        __float128 vi=matlitho->viscosity[layer_index];
     557        __float128 mu0=matlitho->lame_mu[layer_index];
     558        //fetch hypergeometric function tables and parameters
     559        femmodel->parameters->FindParam(&nz,LoveHypergeomNZEnum);
     560        femmodel->parameters->FindParam(&nalpha,LoveHypergeomNAlphaEnum);
     561        femmodel->parameters->FindParam(&z,&dummy1,LoveHypergeomZEnum);
     562        femmodel->parameters->FindParam(&h1,&dummy1,&dummy2,LoveHypergeomTable1Enum);
     563
     564        z1=-(omega*tauh);
     565        z2=-(omega*taul);
     566        //Table1 h1 should be 2F1([1 1+alpha], [2+alpha], z)
     567        hf11=HypergeomTableLookup<__float128>(z1, alpha, h1, z, nz, nalpha);
     568        hf12=HypergeomTableLookup<__float128>(z2, alpha, h1, z, nz, nalpha);
     569
     570        //Ivins et al. (2022) p1979
     571        //Note: therein, mu(s') = s'*mu~(s'); s'=omega*tauM=omega*vi/mu0
     572        factor= alpha*delta/(pow(tauh,alpha)-pow(taul,alpha));
     573        B= factor/(1.0q+alpha) *mu0/vi * (pow(tauh,1.0q+alpha)*hf11 - pow(taul,1.0q+alpha)*hf12);
     574        D= omega*vi/mu0* 1.0q/(1.0q+omega*vi/mu0*(1.0q+delta) -pow(omega*vi/mu0,2.0q)*B);
     575
     576        xDelete<IssmDouble>(z);
     577        xDelete<IssmDouble>(h1);
     578        return mu=mu0*D;
     579}/*}}}*/
     580#endif
     581template <typename doubletype> void        GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega,  Matlitho* matlitho, FemModel* femmodel){ /*{{{*/
    264582
    265583        //returns lame parameters (material rigity) lambda and mu for the right frequency and layer
    266584        doubletype mu,la;
    267585
    268         IssmDouble vi=matlitho->viscosity[layer_index];
    269         IssmDouble mu0=matlitho->lame_mu[layer_index];
    270         IssmDouble la0=matlitho->lame_lambda[layer_index];
     586        doubletype vi=matlitho->viscosity[layer_index];
     587        doubletype mu0=matlitho->lame_mu[layer_index];
     588        doubletype la0=matlitho->lame_lambda[layer_index];
    271589        int rheo=matlitho->rheologymodel[layer_index];
    272590
    273591        if (vi!=0 && omega!=0.0){ //take into account viscosity in the rigidity if the material isn't a perfect fluid
    274                 IssmDouble ka=la0 + 2.0/3.0*mu0; //Bulk modulus
     592                doubletype ka=la0 + 2.0/3.0*mu0; //Bulk modulus
    275593                if (rheo==2){//EBM
    276                         IssmDouble alpha=matlitho->ebm_alpha[layer_index];
    277                         IssmDouble delta=matlitho->ebm_delta[layer_index];
    278                         IssmDouble taul=matlitho->ebm_taul[layer_index];
    279                         IssmDouble tauh=matlitho->ebm_tauh[layer_index];
    280                         IssmDouble hf1,hf2;
    281                         doubletype U1,U2;
    282                         IssmDouble* a=xNew<IssmDouble>(2);
    283                         IssmDouble* b=xNew<IssmDouble>(1);
    284 
    285                         a[0]=1.0;a[1]=1.0+alpha/2.0;
    286                         b[0]=2.0+alpha/2.0;
    287                         //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
    288                         //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
    289                         //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
    290                         //hf2=hypergeometric_pFq_({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
    291                         U1=(pow(tauh,alpha)-pow(taul,alpha))/alpha-pow(omega,2.0)/(2.0+alpha)*(pow(tauh,2.0+alpha)*hf1-pow(taul,2.0+alpha)*hf2);
    292 
    293                         a[0]=1.0;a[1]=.5+alpha/2.0;
    294                         b[0]=1.5+alpha/2.0;
    295                         //hf1=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*tauh,2.0), 0, 0, 15);
    296                         //hf2=HypergeometricFunctionx(a, b, 2, 1, -pow(omega*taul,2.0), 0, 0, 15);
    297                         //hf1=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*tauh,2.0));
    298                         //hf2=hypergeometric_pFq({a[0], a[1]}, {b[0]},-pow(omega*taul,2.0));
    299                         U2=(pow(tauh,1.0+alpha)*hf1-pow(taul,1.0+alpha)*hf2)/(1.0+alpha);
    300 
    301                         mu=mu/(1.0+ alpha*delta/(pow(tauh,alpha)-pow(taul,alpha))*(U1 + omega*U2) -mu/vi/omega);
     594                        mu=muEBM<doubletype>(layer_index, omega, matlitho, femmodel);
    302595                        la=ka-2.0/3.0*mu;
    303 
    304                         xDelete<IssmDouble>(a);
    305                         xDelete<IssmDouble>(b);
    306596                }
    307597                else if (rheo==1){//Burgers
    308                         IssmDouble vi2=matlitho->burgers_viscosity[layer_index];
    309                         IssmDouble mu2=matlitho->burgers_mu[layer_index];
     598                        doubletype vi2=matlitho->burgers_viscosity[layer_index];
     599                        doubletype mu2=matlitho->burgers_mu[layer_index];
    310600
    311601                        mu=mu0*omega*(omega+mu2/vi2)/((omega+mu2/vi2)*(omega+mu0/vi)+mu0/vi2*omega);
    … …  
    321611        mu=mu0;
    322612        }
     613
    323614        *pla=la;
    324615        *pmu=mu;
    325616
    326617} /*}}}*/
    327 IssmDouble                                 GetGravity(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/
     618
     619template <typename doubletype> void        EarthRheology(LoveVariables<doubletype>* vars, IssmDouble* frequencies, int nfreq,  Matlitho* matlitho, FemModel* femmodel){/*{{{*/
     620        doubletype omega;
     621        //reset pointers to NULL if this function was previously called
     622        if(vars->mu)    xDelete<doubletype>(vars->mu);
     623        if(vars->la)    xDelete<doubletype>(vars->la);
     624        //precompute rheology at the requested frequencies
     625        vars->mu=xNewZeroInit<doubletype>(nfreq*matlitho->numlayers);
     626        vars->la=xNewZeroInit<doubletype>(nfreq*matlitho->numlayers);
     627        vars->nfreq=nfreq;
     628        for (int i=0;i<matlitho->numlayers;i++){
     629                for (int fr=0;fr<nfreq;fr++){
     630                        omega=angular_frequency<doubletype>(frequencies[fr]);
     631                        GetEarthRheology<doubletype>(&vars->la[i*nfreq+fr], &vars->mu[i*nfreq+fr], i,omega,matlitho, femmodel);
     632                        //cout << i << " " << fr << " " << vars->mu[i*nfreq+fr] << "\n";
     633                }
     634        }
     635}/*}}}*/
     636
     637template <typename doubletype> doubletype       GetGravity(doubletype r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables<doubletype>* vars){ /*{{{*/
    328638        //computes gravity at radius r2
    329         IssmDouble* EarthMass;
    330         IssmDouble g, GG;
     639        doubletype* EarthMass;
     640        doubletype g, GG;
     641        IssmDouble GGp;
    331642
    332643        EarthMass=vars->EarthMass;
    333         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    334         IssmDouble ro=matlitho->density[layer_index];
    335         IssmDouble M=0;
    336         IssmDouble r1=0;
     644        femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum);
     645        GG=GGp;
     646        doubletype ro=matlitho->density[layer_index];
     647        doubletype M=0;
     648        doubletype r1=0;
    337649        if (layer_index==0){
    338650                M=4.0/3.0*PI*ro*pow(r2,3.0);
    … …  
    344656        return  g= GG*M/pow(r2,2.0);
    345657}/*}}}*/
    346 template <typename doubletype> void        fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
     658template <typename doubletype> void        fill_yi_prefactor(doubletype* yi_prefactor, int* pdeg, doubletype* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
    347659        //precalculates partial derivative factors for function yi_derivatives
    348         IssmDouble ra=matlitho->radius[matlitho->numlayers];
    349         IssmDouble  g0,r0,mu0, GG;
    350         int nstep,nindex, starting_layer;
    351 
    352         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     660        doubletype ra=matlitho->radius[matlitho->numlayers];
     661        doubletype  g0,r0,mu0;
     662        IssmDouble mu0p, GG;
     663        int nstep,nsteps,nindex, starting_layer;
     664
     665        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
    353666        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    354         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     667        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    355668
    356669        g0=vars->g0;
    357670        r0=vars->r0;
     671        mu0=mu0p;
    358672        starting_layer=vars->starting_layer;
    359673
    360674        doubletype frh,frhg0,fgr0,fgr,fn,rm0,rlm,flm;
    361         IssmDouble xmin,xmax,x,dr;
    362         IssmDouble g,ro, issolid;
     675        doubletype xmin,xmax,x,dr;
     676        doubletype g,ro, issolid;
    363677
    364678        if (pomega) { //frequency and degree dependent terms /*{{{*/
    … …  
    370684
    371685                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     686                        nstep=vars->nstep[layer_index];
     687                        nsteps=0;
     688                        for (int i=0;i<layer_index;i++) nsteps+=vars->nstep[i];
    372689
    373690                        ro=matlitho->density[layer_index];
    374691                        issolid=matlitho->issolid[layer_index];
    375692                        if(issolid==1){
    376                                 GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
     693                                //GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
     694                                mu=vars->mu[layer_index*vars->nfreq+vars->ifreq];
     695                                la=vars->la[layer_index*vars->nfreq+vars->ifreq];
    377696
    378697                                /*_______Expressions*/
    … …  
    386705                                f[2]=(la*fn/flm);
    387706                                f[3]=rm0*rlm;
    388                                 f[4]=ro*pow(omega,2.0)*ra/mu0;
     707                                f[4]=-ro*pow(omega,2.0)*ra*ra/mu0;
    389708                                f[5]=(-4.0*mu/flm);
    390709                                f[6]=fn*frh;
    … …  
    399718                                dr = (xmax -xmin)/nstep;
    400719                                x=xmin;
     720
     721                                //fixme
     722                                g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
     723
    401724                                for (int n=0;n<nstep;n++){
    402                                         g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
    403 
    404                                         nindex=layer_index*nstep*36+n*36;
     725
     726                                        g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
     727                                        nindex=nsteps*36+n*36;
    405728                                        yi_prefactor[nindex+ 0*6+0]= f[0]/x;                      // in dy[0*6+0]
    406729                                        yi_prefactor[nindex+ 0*6+1]= f[1];                        // in dy[0*6+1]
    … …  
    423746
    424747                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     748                        nstep=vars->nstep[layer_index];
     749                        nsteps=0;
     750                        for (int i=0;i<layer_index;i++) nsteps+=vars->nstep[i];
     751
    425752                        ro=matlitho->density[layer_index];
    426753                        issolid=matlitho->issolid[layer_index];
    … …  
    433760                        dr = (xmax -xmin)/nstep;
    434761                        x=xmin;
     762
     763
     764                                //fixme
     765                                g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
     766
    435767                        for (int n=0;n<nstep;n++){
    436                                 nindex=layer_index*nstep*36+n*36;
    437                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
     768                                nindex=nsteps*36+n*36;
     769                                g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
    438770
    439771                                if(issolid==1){
    … …  
    450782        } else { // static terms /*{{{*/
    451783                for (int layer_index=starting_layer;layer_index<matlitho->numlayers;layer_index++){
     784                        nstep=vars->nstep[layer_index];
     785                        nsteps=0;
     786                        for (int i=0;i<layer_index;i++) nsteps+=vars->nstep[i];
     787
    452788                        ro=matlitho->density[layer_index];
    453789                        issolid=matlitho->issolid[layer_index];
    … …  
    461797                        dr = (xmax -xmin)/nstep;
    462798                        x=xmin;
     799                                //fixme
     800                                g=GetGravity<doubletype>((xmin+xmax)/2*ra,layer_index,femmodel,matlitho,vars);
    463801                        for (int n=0;n<nstep;n++){
    464                                 g=GetGravity(x*ra,layer_index,femmodel,matlitho,vars);
    465                                 nindex=layer_index*nstep*36+n*36;
     802                                g=GetGravity<doubletype>(x*ra,layer_index,femmodel,matlitho,vars);
     803                                nindex=nsteps*36+n*36;
    466804                                if(issolid==1){
    467805                                        yi_prefactor[nindex+ 1*6+5]= -frhg0;       // in dy[1*6+5]
    … …  
    484822        }
    485823}/*}}}*/
    486 template <typename doubletype> void        yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
     824template <typename doubletype> void        yi_derivatives(doubletype* dydx, doubletype* y, int layer_index, int n, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
    487825        //computes yi derivatives at r=radius[layer_index]+ n/nstep*(radius[layer_index+1]-radius[layer_index])
    488826
    489         IssmDouble issolid=matlitho->issolid[layer_index];
    490         int iy,id,ny, nindex, nstep;
    491         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     827        int issolid=matlitho->issolid[layer_index];
     828        int iy,id,ny, nindex, nstep, nsteps;
     829        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     830
     831        nstep=vars->nstep[layer_index];
     832        nsteps=0;
     833        for (int i=0;i<layer_index;i++) nsteps+=vars->nstep[i];
    492834
    493835        /*{{{*/ /* For reference:
    … …  
    552894                           dy[1*6+0]= (-4.0*(fgr/rg0)+fn/x)/x;
    553895                           dy[1*6+1]= -2.0/x-fgr/rg0;
     896
    554897                           }
    555898        */ /*}}}*/
     899        nindex=nsteps*36+n*36;
    556900
    557901        if(issolid==1){
    … …  
    564908                dydx[id]=0.0;
    565909                for (iy=0;iy<ny;iy++){
    566                         nindex=layer_index*nstep*36+n*36;
    567910                        dydx[id]+=yi_prefactor[nindex+id*6+iy]*y[iy];
    568911                }
    … …  
    570913        return;
    571914}/*}}}*/
    572 template <typename doubletype> void        propagate_yi_euler(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
    573         //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?
     915template <typename doubletype> void        propagate_yi_euler(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
     916        //computes 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?
    574917        //euler method
    575918        int nstep;
    576         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     919        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     920        nstep=vars->nstep[layer_index];
     921
    577922
    578923        doubletype* dydx=xNewZeroInit<doubletype>(6);
    579         IssmDouble dr = (xmax -xmin)/nstep;
    580         IssmDouble x=xmin;
     924        doubletype dr = (xmax -xmin)/nstep;
     925        doubletype x=xmin;
    581926        for(int i = 0;i<nstep;i++){
    582                 yi_derivatives<doubletype>(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho);
     927                yi_derivatives<doubletype>(dydx,y,layer_index, i,yi_prefactor,femmodel,matlitho,vars);
    583928                for (int j=0;j<6;j++){
    584929                        y[j]+=dydx[j]*dr;
    … …  
    588933        xDelete<doubletype>(dydx);
    589934}/*}}}*/
    590 template <typename doubletype> void        propagate_yi_RK2(doubletype* y, IssmDouble xmin, IssmDouble xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
    591         //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?
     935template <typename doubletype> void        propagate_yi_RK2(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
     936        //computes 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?
    592937        //Implements Runge-Kutta 2nd order (midpoint method)
    593938        int nstep;
    594         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     939        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     940        nstep=vars->nstep[layer_index];
    595941
    596942        doubletype k1[6]={0};
    … …  
    602948        doubletype y3[6]={0};
    603949
    604         IssmDouble dr = (xmax -xmin)/nstep;
    605         IssmDouble x=xmin;
     950        doubletype dr = (xmax -xmin)/nstep;
     951        doubletype x=xmin;
    606952
    607953        for(int i = 0;i<nstep/2;i++){
    608                 yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     954                yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars);
    609955                for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
    610                 yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);           
     956                yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);               
    611957
    612958                for (int j=0;j<6;j++){
    … …  
    616962        }
    617963}/*}}}*/
    618 template <typename doubletype> void        propagate_yi_RK4(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho){ /*{{{*/
    619         //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?
     964        template <typename doubletype> void        propagate_yi_RK4(doubletype* y, doubletype xmin, doubletype xmax, int layer_index, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<doubletype>* vars){ /*{{{*/
     965        //computes 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?
    620966        //Implements Runge-Kutta 4th order
    621967        int nstep;
    622         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    623 
    624         IssmDouble* k1=xNewZeroInit<IssmDouble>(6);
    625         IssmDouble* k2=xNewZeroInit<IssmDouble>(6);
    626         IssmDouble* k3=xNewZeroInit<IssmDouble>(6);
    627         IssmDouble* k4=xNewZeroInit<IssmDouble>(6);
    628         IssmDouble* y1=xNewZeroInit<IssmDouble>(6);
    629         IssmDouble* y2=xNewZeroInit<IssmDouble>(6);
    630         IssmDouble* y3=xNewZeroInit<IssmDouble>(6);
    631 
    632         IssmDouble dr = (xmax -xmin)/nstep;
    633         IssmDouble x=xmin;
     968        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     969        nstep=vars->nstep[layer_index];
     970
     971        doubletype k1[6]={0};
     972        doubletype k2[6]={0};
     973        doubletype k3[6]={0};
     974        doubletype k4[6]={0};
     975        doubletype y1[6]={0};
     976        doubletype y2[6]={0};
     977        doubletype y3[6]={0};
     978
     979        doubletype dr = (xmax -xmin)/nstep;
     980        doubletype x=xmin;
    634981        for(int i = 0;i<nstep/2-1;i++){
    635                 yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     982                yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars);
    636983                for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
    637                 yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     984                yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);
    638985                for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
    639                 yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     986                yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);
    640987                for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
    641                 yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho);           
     988                yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+2,yi_prefactor,femmodel,matlitho,vars);               
    642989
    643990                for (int j=0;j<6;j++){
    … …  
    649996        //Last step: we don't know the derivative at xmax, so we will assume the values at xmax-dr
    650997        int i=nstep/2;
    651         yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho);
     998        yi_derivatives<doubletype>(k1,y,layer_index, 2*i,yi_prefactor,femmodel,matlitho,vars);
    652999        for (int j=0;j<6;j++) {y1[j]=y[j]+k1[j]*dr;}
    653         yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     1000        yi_derivatives<doubletype>(k2,y1,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);
    6541001        for (int j=0;j<6;j++) {y2[j]=y[j]+k2[j]*dr;}
    655         yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);
     1002        yi_derivatives<doubletype>(k3,y2,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);
    6561003        for (int j=0;j<6;j++) {y3[j]=y[j]+k3[j]*2.0*dr;}
    657         yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho);           
     1004        yi_derivatives<doubletype>(k4,y3,layer_index, 2*i+1,yi_prefactor,femmodel,matlitho,vars);               
    6581005
    6591006        for (int j=0;j<6;j++){
    … …  
    6631010        x = x + 2.0*dr;
    6641011
    665         xDelete<IssmDouble>(k1);
    666         xDelete<IssmDouble>(k2);
    667         xDelete<IssmDouble>(k3);
    668         xDelete<IssmDouble>(k4);
    669         xDelete<IssmDouble>(y1);
    670         xDelete<IssmDouble>(y2);
    671         xDelete<IssmDouble>(y3);
    6721012}/*}}}*/
    673 template <typename doubletype> void        Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
     1013template <typename doubletype> void        Innersphere_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
    6741014        //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5]
    6751015
    676         IssmDouble r = matlitho->radius[layer_index];
    677         IssmDouble ra=matlitho->radius[matlitho->numlayers];
    678         IssmDouble  g0,r0,mu0, GG;
    6791016        int nyi;
    680 
    681         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
    682         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
     1017        doubletype r = matlitho->radius[layer_index];
     1018        doubletype ra=matlitho->radius[matlitho->numlayers];
     1019        doubletype  g0,r0,mu0, GG;
     1020        IssmDouble mu0p, GGp;
     1021
     1022
     1023        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
     1024        femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum);
    6831025
    6841026        g0=vars->g0;
    6851027        r0=vars->r0;
     1028        mu0=mu0p;
     1029        GG=GGp;
    6861030        nyi=vars->nyi;
    6871031
    688         IssmDouble ro=matlitho->density[layer_index];
    689         IssmDouble issolid=matlitho->issolid[layer_index];
    690         IssmDouble g=GetGravity(r,layer_index,femmodel,matlitho,vars);
    691         doubletype la,mu;
    692 
    693         if (layer_index==0){
    694                 // 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
    695                 GetEarthRheology<doubletype>(&la, &mu,layer_index,omega,matlitho);   
    696         } else {
    697                 GetEarthRheology<doubletype>(&la, &mu,layer_index-1,omega,matlitho);   
    698         }   
    699 
    700         IssmDouble cst = 4.0*PI*GG*ro;
    701         IssmDouble r2=pow(r,2.0);
     1032
     1033        doubletype g=GetGravity<doubletype>(r,layer_index,femmodel,matlitho,vars);
     1034        doubletype la,mu,ro;
     1035       
     1036        int i=layer_index-1;
     1037        if (layer_index==0) i=layer_index;
     1038
     1039        ro=matlitho->density[i];
     1040
     1041        //elastic values
     1042        la=matlitho->lame_lambda[i];
     1043        mu=matlitho->lame_mu[i];
     1044        doubletype Kappa=(la+2.0/3.0*mu);
     1045
     1046        //update to viscoelastic values
     1047        mu=vars->mu[i*vars->nfreq+vars->ifreq];
     1048        la = Kappa-2.0/3.0*mu;
     1049
     1050        doubletype cst = 4.0*PI*GG*ro;
     1051        doubletype r2=pow(r,2.0);
     1052
     1053        //Greff-Lefftz and Legros 1997, p701, analytical solution for incompressible elastic layer for y3, y4, y5 ensuring they =0 at r=0
     1054        //These equations are then divided by r^n for numerical stability at higher degrees
     1055
     1056        //all terms in y1 y2 y6 are 0 in that layer because they are of the type r^l with l<0 and would diverge at the origin, that's why we only need 3 equations
    7021057
    7031058        yi[0+nyi*0]=1.0*r/ra;
    … …  
    7251080        yi[5+nyi*2]=deg/(r*g0);
    7261081
     1082
     1083        /*doubletype vp2 = (la + 2.0*mu)/ro;
     1084        yi[0+nyi*0]=1.0*r/ra;
     1085        yi[0+nyi*1]=1.0/(r*ra);
     1086        yi[0+nyi*2]=0.0;
     1087
     1088        yi[1+nyi*0]=(2.0*mu*(deg-1.0-3.0/deg) + cst/3.0*ro*r2)/mu0;
     1089        yi[1+nyi*1]=(2.0*mu*(deg-1.0)/r2 + cst/3.0*ro)/mu0;
     1090        yi[1+nyi*2]=-ro/mu0;
     1091
     1092        yi[2+nyi*0]=(deg+3.0)/(deg*(deg+1.0))*r/ra;
     1093        yi[2+nyi*1]=1.0/(deg*r*ra);
     1094        yi[2+nyi*2]=0.0;
     1095
     1096        yi[3+nyi*0]=2.0*mu*(deg+2.0)/((deg+1.0)*mu0);
     1097        yi[3+nyi*1]=2.0*mu*(deg-1.0)/(deg*r2*mu0);
     1098        yi[3+nyi*2]=0.0;
     1099
     1100        yi[4+nyi*0]=0.0;
     1101        yi[4+nyi*1]=0.0;
     1102        yi[4+nyi*2]=1.0/(g0*ra);
     1103
     1104        yi[5+nyi*0]=-cst*r/g0;
     1105        yi[5+nyi*1]=-cst/(r*g0);
     1106        yi[5+nyi*2]=deg/(r*g0);*/
     1107
     1108
     1109
    7271110}/*}}}*/
    728 template <typename doubletype> void        build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars) { /*{{{*/
    729 
    730         IssmDouble  g0,r0,mu0,x,ro1, GG;
    731         int nyi,starting_layer, nstep;
     1111template <typename doubletype> void        Coremantle_boundaryconditions(doubletype* yi, int layer_index, int deg, doubletype omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars){ /*{{{*/
     1112        //fills the boundary conditions at the bottom of layer[layer_index] in yi[0:2][0:5]
     1113
     1114        int nyi;
     1115        doubletype r = matlitho->radius[layer_index];
     1116        doubletype ra=matlitho->radius[matlitho->numlayers];
     1117        doubletype  g0,r0,mu0, GG;
     1118        IssmDouble mu0p, GGp;
     1119
     1120
     1121        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
     1122        femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum);
     1123
     1124        g0=vars->g0;
     1125        r0=vars->r0;
     1126        mu0=mu0p;
     1127        GG=GGp;
     1128        nyi=vars->nyi;
     1129
     1130        doubletype ro=matlitho->density[layer_index-1];
     1131
     1132       
     1133        if (!matlitho->issolid[layer_index]) _error_("Love core error: CMB conditions requested but layer " << layer_index << " (mantle side) is not solid");
     1134        if (matlitho->issolid[layer_index-1]) _error_("Love core error: CMB conditions requested but layer " << layer_index-1 << " (outer core) is solid");
     1135
     1136        doubletype cst = 4.0/3.0*PI*GG*ro;
     1137        doubletype rl1=pow(r,deg-1);
     1138
     1139        yi[0+nyi*0]=-rl1/cst/ra;
     1140        yi[0+nyi*1]=1.0/ra;
     1141        yi[0+nyi*2]=0.0;
     1142
     1143        yi[1+nyi*0]=0.0;
     1144        yi[1+nyi*1]=ro*cst*r/mu0;
     1145        yi[1+nyi*2]=0.0;
     1146
     1147        yi[2+nyi*0]=0.0;
     1148        yi[2+nyi*1]=0.0;
     1149        yi[2+nyi*2]=1.0/ra;
     1150
     1151        yi[3+nyi*0]=0.0;
     1152        yi[3+nyi*1]=0.0;
     1153        yi[3+nyi*2]=0.0;
     1154
     1155        yi[4+nyi*0]=r*rl1/(g0*ra);
     1156        yi[4+nyi*1]=0.0;
     1157        yi[4+nyi*2]=0.0;
     1158
     1159        yi[5+nyi*0]=2.0*(deg-1)*rl1/g0;
     1160        yi[5+nyi*1]=3.0*cst/g0;
     1161        yi[5+nyi*2]=0.0;
     1162
     1163}/*}}}*/
     1164template <typename doubletype> void        build_yi_system(doubletype* yi, int deg, doubletype omega, doubletype* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<doubletype>* vars) { /*{{{*/
     1165
     1166        doubletype      g0,r0,mu0,x,ro1, GG;
     1167        int             nyi,starting_layer, nstep;
     1168        doubletype      xmin,xmax,one,ro,g, ra;
     1169        IssmDouble      mu0p, GGp;
     1170        bool            debug;
     1171        int ny,is,ii,jj;
     1172        int scheme;
     1173        int ici = 0;   // Index of current interface
     1174        int cmb=0;
     1175
     1176        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
     1177        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
     1178        femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum);
     1179        femmodel->parameters->FindParam(&debug,LoveDebugEnum);
     1180        femmodel->parameters->FindParam(&scheme,LoveIntegrationSchemeEnum);
    7321181
    7331182        g0=vars->g0;
    … …  
    7351184        nyi=vars->nyi;
    7361185        starting_layer=vars->starting_layer;
    737 
    738         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
    739         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    740         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    741 
    742         IssmDouble xmin,xmax,one,ro,g;
    743         IssmDouble ra=matlitho->radius[matlitho->numlayers];
     1186        mu0=mu0p;
     1187        GG=GGp;
     1188        ra=matlitho->radius[matlitho->numlayers];
    7441189
    7451190        for (int i=0;i<6*(matlitho->numlayers+1);i++){
    … …  
    7491194        }
    7501195
    751         int ny,is,ii,jj;
    7521196        doubletype ystart[6];
    7531197        for (int k=0;k<6;k++) ystart[k]=0.0;           
    7541198
    755         int ici = 0;   // Index of current interface
     1199
     1200
    7561201        for (int i = starting_layer; i<matlitho->numlayers;i++){
    757 
     1202                ici=i-starting_layer;
    7581203                xmin=matlitho->radius[i]/ra;
    7591204                xmax=(matlitho->radius[i+1])/ra;
    … …  
    7741219
    7751220                        // Numerical Integration
    776                         //propagate_yi_euler(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
    777                         propagate_yi_RK2<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
    778                         //propagate_yi_RK4(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho);
     1221                        if (debug) propagate_yi_euler<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars);
     1222                        else {
     1223                                if (scheme==0) propagate_yi_euler<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars);
     1224                                else if (scheme==1) propagate_yi_RK2<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars);
     1225                                else if (scheme==2) propagate_yi_RK4<doubletype>(&ystart[0], xmin, xmax, i, yi_prefactor,femmodel, matlitho, vars);
     1226                                else _error_("Love core error: integration scheme not found");
     1227                        }
    7791228                        // Boundary Condition matrix - propagation part
    7801229                        ii = 6*(ici+1)+is;
    … …  
    7941243                } else { // Boundary Condition matrix - liquid regions
    7951244                        ro1=matlitho->density[i];
    796                         g=GetGravity(matlitho->radius[i], i, femmodel,matlitho,vars);
     1245                        g=GetGravity<doubletype>(matlitho->radius[i], i, femmodel,matlitho,vars);
    7971246                        ii = 6*ici;
    798                         yi[ii+nyi*(ii+3)] = -1.0;
    799                         yi[ii+nyi*(ii+4+3)] = -g0/g;
    800                         yi[(ii+1)+nyi*(ii+3)]=-ro1*g*ra/mu0;
    801                         yi[(ii+2)+nyi*(ii+1+3)]=-1.0;
    802                         yi[(ii+5)+nyi*(ii+3)]= 4.0*PI*GG*ro1*ra/g0;
    803                         yi[(ii+4)+nyi*(ii+4+3)]=-1.0;
    804                         yi[(ii+5)+nyi*(ii+5+3)]=-1.0;
    805                         g=GetGravity(matlitho->radius[i+1], i,femmodel,matlitho,vars);
     1247                        jj = 6*ici+3;
     1248                        yi[ii+nyi*(jj)] = -1.0;
     1249                        yi[ii+nyi*(jj+4)] = -g0/g;
     1250                        yi[(ii+1)+nyi*(jj)]=-ro1*g*ra/mu0;
     1251                        yi[(ii+2)+nyi*(jj+1)]=-1.0;
     1252                        yi[(ii+5)+nyi*(jj)]= 4.0*PI*GG*ro1*ra/g0;
     1253                        yi[(ii+4)+nyi*(jj+4)]=-1.0;
     1254                        yi[(ii+5)+nyi*(jj+5)]=-1.0;
     1255                        g=GetGravity<doubletype>(matlitho->radius[i+1], i,femmodel,matlitho,vars);
    8061256                        ii = 6*(ici+1);
    807                         yi[ii+nyi*(ii-1)]=-1.0;
    808                         yi[ii+nyi*(ii+1)]=yi[(ii+4)+nyi*(ii+1)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB
    809                         yi[ii+nyi*(ii+2)]=yi[(ii+4)+nyi*(ii+2)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB
     1257
     1258                        yi[ii+nyi*(jj+2)]=-1.0;
     1259                        yi[ii+nyi*(jj+4)]=yi[(ii+4)+nyi*(jj+4)]*g0/g; // yi(17,14) solution integration 1 of z5 CMB
     1260                        yi[ii+nyi*(jj+5)]=yi[(ii+4)+nyi*(jj+5)]*g0/g; // yi(17,15) solution integration 2 of z5 CMB
    8101261                        // yi(13,..) y1 CMB
    811                         yi[(ii+1)+nyi*(ii-1)]=-ro1*g*ra/mu0;
    812                         yi[(ii+2)+nyi*(ii)]=-1.0;
    813                         yi[(ii+5)+nyi*(ii-1)]= 4.0*PI*GG*ro1*ra/g0;
     1262                        yi[(ii+1)+nyi*(jj+2)]=-ro1*g*ra/mu0;
     1263                        yi[(ii+2)+nyi*(jj+3)]=-1.0;
     1264                        yi[(ii+5)+nyi*(jj+2)]= 4.0*PI*GG*ro1*ra/g0;
    8141265                }       
    8151266                ici = ici+1;
    … …  
    8171268
    8181269        //-- Internal sphere: integration starts here rather than r=0 for numerical reasons
     1270        /*if (starting_layer==cmb) Coremantle_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
     1271        else Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);*/
     1272
    8191273        Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
    8201274
    … …  
    8341288
    8351289}/*}}}*/
    836 template <typename doubletype> void        yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type){ /*{{{*/
    837 
    838         IssmDouble  g0,r0,mu0,ra,rb,rc;
     1290template <typename doubletype> void        yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables<doubletype>* vars, int forcing_type){ /*{{{*/
     1291
     1292        doubletype  g0,r0,mu0,ra,rb,rc;
    8391293        int nyi,icb,cmb,starting_layer;
    840         IssmDouble* EarthMass;
     1294        doubletype* EarthMass;
     1295        IssmDouble mu0p;
    8411296
    8421297        g0=vars->g0;
    … …  
    8461301        EarthMass=vars->EarthMass;
    8471302
    848         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     1303        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
    8491304        femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum);
    8501305        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
    8511306
     1307        mu0=mu0p;
    8521308        // In Case of a Inner Core - Outer Core - Mantle planet and Boundary conditions on these 3 interfaces
    8531309        ra=matlitho->radius[matlitho->numlayers];       
    … …  
    8611317        }
    8621318
    863         IssmDouble ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0));
     1319        doubletype ro_mean=EarthMass[matlitho->numlayers-1]/(4.0/3.0*PI*pow(ra,3.0));
    8641320
    8651321        for (int i=0;i<(matlitho->numlayers+1)*6;i++) yi_righthandside[i]=0.0;
    … …  
    9171373        }
    9181374}/*}}}*/
    919 template <typename doubletype> void        solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, IssmDouble* frequencies, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){ /*{{{*/
    920 
    921         IssmDouble  g0,r0,mu0,loveratio,underflow_tol;
     1375template <typename doubletype> void        solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, IssmDouble* frequencies, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars, bool verbosecpu){ /*{{{*/
     1376
     1377        doubletype  g0,r0,mu0;
    9221378        //IssmDouble* frequencies;
    923         int nyi,starting_layer, dummy;
    924         bool allow_layer_deletion;
    925         IssmDouble* EarthMass=NULL;
     1379        int nyi,starting_layer, dummy,cmb;
     1380        bool allow_layer_deletion, debug;
     1381        doubletype* EarthMass=NULL;
     1382        IssmDouble mu0p,loveratio,underflow_tol;
    9261383
    9271384        g0=vars->g0;
    … …  
    9311388        EarthMass=vars->EarthMass;
    9321389
    933         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
     1390        femmodel->parameters->FindParam(&mu0p,LoveMu0Enum);
    9341391        femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
    9351392        femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum);
     1393        femmodel->parameters->FindParam(&debug,LoveDebugEnum);
     1394        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
    9361395        //femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
    937         IssmDouble ra=matlitho->radius[matlitho->numlayers];
     1396        mu0=mu0p;
     1397        doubletype ra=matlitho->radius[matlitho->numlayers];
    9381398        bool exit=false;
    9391399        int lda,ldb;
    … …  
    9531413                        }
    9541414                }
     1415               
     1416                if (debug){
     1417                        IssmDouble*  yidebug=xNew<IssmDouble>(nyi*nyi);
     1418                        IssmDouble*  rhsdebug=xNew<IssmDouble>(nyi);
     1419                        for (int i=0;i<nyi;i++){
     1420                                rhsdebug[i]=DownCastVarToDouble(rhs[i]);
     1421                                for (int j=0;j<nyi;j++){
     1422                                        yidebug[i+j*nyi]=DownCastVarToDouble(yi[i+j*nyi]);
     1423                                }
     1424                        }
     1425                        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveYiEnum,yidebug,nyi,nyi,0,0));
     1426                        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveRhsEnum,rhsdebug,nyi,1,0,0));
     1427                        xDelete<IssmDouble>(yidebug);
     1428                        xDelete<IssmDouble>(rhsdebug);
     1429                }
     1430
    9551431                //-- Resolution
    9561432                int* ipiv=xNewZeroInit<int>(nyi); //pivot index vector
    … …  
    9741450                        _error_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");*/
    9751451
    976                 if(VerboseSolution() && info!=0){
    977                         _printf_("i j yi[i+nyi*j] rhs[i]");
     1452                if(VerboseSolution() && verbosecpu && info!=0){
     1453                        _printf_("i j yi[i+nyi*j] rhs[i]\n");
    9781454                        for (int i=0;i<nyi;i++){
    9791455                                for (int j=0;j<nyi;j++){
    … …  
    9811457                                }
    9821458                        }
    983                         _error_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
    984                 }
     1459                        _error_("love core error in DGESV : LAPACK linear equation solver couldn't resolve the system");
     1460                }
     1461
    9851462
    9861463                *loveh = rhslocal[nyi-3]*ra*g0;
    … …  
    10001477                if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s);
    10011478
    1002                
     1479                if (debug) goto save_results;
    10031480
    10041481                if (!allow_layer_deletion || nyi<=12 || omega!=angular_frequency<doubletype>(frequencies[0]) || deg==0){
    … …  
    10221499
    10231500                if (deg==vars->deg_layer_delete[starting_layer]){ // if we are at the degree where we should delete the current layer, proceed to delete the bottom layer
    1024                         if (omega!=0 && VerboseSolution()  && verbosecpu) _printf_(", deleting layer " << starting_layer << "\n");
     1501                        //if (omega!=0 && VerboseSolution()  && verbosecpu) _printf_(", deleting layer " << starting_layer << "\n");
    10251502                        nyi-=6;
    10261503                        starting_layer+=1;
    … …  
    10351512                        }
    10361513
    1037                         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
     1514        /*if (starting_layer==cmb) Coremantle_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
     1515        else 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
     1516        */
     1517        Innersphere_boundaryconditions<doubletype>(yi, starting_layer, deg, omega, femmodel, matlitho,vars);
    10381518                } else { //we are ready to save the outputs and break the main loop
    10391519
    … …  
    11341614                        postwidder_transform<doubletype>(pmtf_orthot,pmtf_orthof,2,t,2,NTit,xi,femmodel);
    11351615                }
     1616                xDelete<doubletype>(pmtf_colinearf);
     1617                xDelete<doubletype>(pmtf_orthof);
    11361618                if(VerboseSolution() && verbosecpu) _printf_("done!\n");
    11371619        }
    11381620
    11391621        xDelete<doubletype>(xi);
    1140         xDelete<doubletype>(pmtf_colinearf);
    1141         xDelete<doubletype>(pmtf_orthof);
    11421622}/*}}}*/
    11431623
    1144 template <typename doubletype> void        compute_love_numbers(LoveNumbers<doubletype>* Lovef, LoveNumbers<doubletype>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu){
    1145 
    1146         int nstep, kernel_index,kernel_indexe,deleted_layer_offset, deg, sh_nmin, sh_nmax, nfreq;
     1624template <typename doubletype> void        compute_love_numbers(LoveNumbers<doubletype>* Lovef, LoveNumbers<doubletype>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables<doubletype>* vars, bool verbosecpu){
     1625
     1626        int nsteps, kernel_index,kernel_indexe,deleted_layer_offset, deg, sh_nmin, sh_nmax, nfreq;
    11471627        doubletype  lovek, loveh, lovel, loveratio;
    11481628        doubletype  omega;
    … …  
    11501630        doubletype* yi_righthandside=NULL;
    11511631        doubletype* yi=NULL;
    1152         IssmDouble  underflow_tol;
     1632        doubletype  underflow_tol;
     1633        IssmDouble dr;
    11531634        bool freq_skip, istemporal;
    1154 
    1155         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     1635        int cmb=0;
     1636        int nyi_init=0;
     1637
     1638        //femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    11561639        femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
     1640        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
    11571641
    11581642        nfreq=Lovef->nfreq;
    … …  
    11621646
    11631647        // reset deleted layers in case we have called this function before;
    1164         vars->starting_layer=0;
    1165         vars->nyi=6*(matlitho->numlayers+1);
    1166 
    1167         yi_prefactor=xNewZeroInit<doubletype>(6*6*nstep*matlitho->numlayers);
    1168         yi_righthandside=xNewZeroInit<doubletype>(vars->nyi);
    1169         yi=xNewZeroInit<doubletype>(vars->nyi*vars->nyi);
    1170 
    1171         //precalculate yi coefficients that do not depend on degree or frequency
    1172         fill_yi_prefactor<doubletype>(yi_prefactor, NULL, NULL,femmodel, matlitho,vars);
     1648        vars->starting_layer=0;
     1649        vars->nyi=6*(matlitho->numlayers-vars->starting_layer+1);
     1650        nyi_init=6*(matlitho->numlayers+1);
     1651        nsteps=0;
     1652        for (int i=0;i<matlitho->numlayers;i++) nsteps+=vars->nstep[i];
     1653
     1654        //yi_prefactor=xNewZeroInit<doubletype>(6*6*nstep*matlitho->numlayers);
     1655        yi_prefactor=xNewZeroInit<doubletype>(6*6*nsteps);
     1656        yi_righthandside=xNewZeroInit<doubletype>(nyi_init);
     1657        yi=xNewZeroInit<doubletype>(nyi_init*nyi_init);
     1658
     1659        //precompute yi coefficients that do not depend on degree or frequency
     1660        fill_yi_prefactor<doubletype>(yi_prefactor, NULL, NULL,femmodel, matlitho,vars);
    11731661
    11741662        if (VerboseSolution() && Elastic  && verbosecpu) _printf_("\n");
    11751663
    1176         for(int deg=0;deg<2;deg++){ // calculation is in the center of mass reference frame, neutralize degree 0 and 1 mass changes
     1664        for(int deg=0;deg<2;deg++){ // calculation is in the center of mass reference frame, neutralize degree 0 and 1 mass changes, i.e 1+k=0
    11771665                for (int fr=0;fr<nfreq;fr++){
    11781666                        Lovef->K[fr*(sh_nmax+1)+deg]=-1.0;
    … …  
    11851673                }
    11861674
    1187                 //precalculate yi coefficients that depend on degree but not frequency
     1675                //precompute yi coefficients that depend on degree but not frequency
    11881676                fill_yi_prefactor<doubletype>(yi_prefactor, &deg, NULL,femmodel, matlitho,vars);
    11891677
    11901678                for (int fr=0;fr<nfreq;fr++){
    11911679                        omega=angular_frequency<doubletype>(frequencies[fr]);
     1680                        vars->ifreq=fr;
    11921681                       
    1193                         //precalculate yi coefficients that depend on degree and frequency
     1682                        //precompute yi coefficients that depend on degree and frequency
    11941683                        fill_yi_prefactor<doubletype>(yi_prefactor, &deg,&omega,femmodel, matlitho,vars);
    11951684
    11961685                        //solve the system
    1197                         yi_boundary_conditions<doubletype>(yi_righthandside,deg,femmodel,matlitho,vars,forcing_type); 
     1686                        yi_boundary_conditions<doubletype>(yi_righthandside,deg,femmodel,matlitho,vars,forcing_type);
    11981687                        build_yi_system<doubletype>(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars);
    1199                         solve_yi_system<doubletype>(&loveh,&lovel,&lovek, deg, omega, frequencies, yi, yi_righthandside,femmodel, matlitho,vars,verbosecpu);
     1688                        solve_yi_system<doubletype>(&loveh,&lovel,&lovek, deg, omega, frequencies, yi, yi_righthandside,femmodel, matlitho,vars,verbosecpu && !Elastic);
     1689
    12001690                        Lovef->H[fr*(sh_nmax+1)+deg]=loveh;
    12011691                        Lovef->K[fr*(sh_nmax+1)+deg]=lovek-1.0;
    … …  
    12091699        }
    12101700
    1211         if (Elastic) { // if elastic values were provided, we copy elastic love numbers above the cutoff degree instead of calculating them
     1701        if (Elastic) { // if elastic values were provided, we copy elastic love numbers above the cutoff degree instead of computing them
    12121702                for(int deg=sh_cutoff+1;deg<sh_nmax+1;deg++){
    12131703                        if (VerboseSolution() && Elastic  && verbosecpu) {
    … …  
    12291719        }
    12301720
     1721
     1722
    12311723        if (VerboseSolution() && Elastic  && verbosecpu) _printf_("\n");
    12321724        xDelete<doubletype>(yi);
    … …  
    12361728
    12371729/*templated cores:*/
    1238 LoveVariables*                             love_init(FemModel* femmodel, Matlitho* matlitho){/*{{{*/
     1730template <typename doubletype> LoveVariables<doubletype>*       love_init(FemModel* femmodel, Matlitho* matlitho, bool verbosecpu){/*{{{*/
    12391731
    12401732        /*initialize Planet_Mass(r) for efficient computation of gravity, value of surface gravity and inital size of the yi equation system*/
    … …  
    12421734        bool        verbosemod = (int)VerboseModule();
    12431735        int         numlayers  = matlitho->numlayers;
    1244         IssmDouble* r          = matlitho->radius;
    1245         IssmDouble  r1,r2,ro, GG;
     1736        int         minsteps;
     1737        doubletype* r=NULL;
     1738        doubletype  r1,r2,ro, GG;
     1739        IssmDouble GGp;
     1740        IssmDouble dr;
    12461741
    12471742        /*outputs:*/
    1248         IssmDouble* EarthMass=NULL;
    1249         IssmDouble  g0,r0;
    1250         int         nyi,starting_layer;
     1743        doubletype* EarthMass=NULL;
     1744        doubletype  g0,r0;
     1745        int         nyi,starting_layer,cmb;
    12511746        int*        deg_layer_delete;
    1252 
    1253         femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    1254         EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
     1747        int*        nstep;
     1748
     1749       
     1750        femmodel->parameters->FindParam(&GGp,LoveGravitationalConstantEnum);
     1751        femmodel->parameters->FindParam(&minsteps, LoveMinIntegrationStepsEnum);
     1752        femmodel->parameters->FindParam(&dr, LoveMaxIntegrationdrEnum);
     1753        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
     1754        GG=GGp;
     1755        EarthMass=xNewZeroInit<doubletype>(numlayers+1);
    12551756        deg_layer_delete=xNewZeroInit<int>(numlayers);
     1757
     1758        r=xNewZeroInit<doubletype>(numlayers+1);
     1759        nstep=xNewZeroInit<int>(numlayers);
     1760        for (int i=0;i<numlayers+1;i++){
     1761                r[i] = matlitho->radius[i];
     1762                if (i<numlayers) {
     1763                        // nstep[i] is the largest even integer such that (radius[i+1]-radius[i])/nstep[i]<dr
     1764                        nstep[i]=ceil((matlitho->radius[i+1]-matlitho->radius[i])/dr/2)*2;
     1765                        if (nstep[i]<minsteps) nstep[i]=minsteps;
     1766                }
     1767        }
    12561768
    12571769        for (int i=0;i<numlayers;i++){
    … …  
    12651777                }
    12661778        }
    1267 
    12681779        g0=EarthMass[numlayers-1]*GG/pow(r[numlayers],2.0);
    12691780        r0=r[numlayers];
    1270         nyi=6*(numlayers+1);
    12711781        starting_layer=0;
    1272 
    1273         if(VerboseSolution()){
     1782        nyi=6*(numlayers-starting_layer+1);
     1783
     1784
     1785        if(VerboseSolution() && verbosecpu){
    12741786                _printf_("     Surface gravity: " << g0 << " m.s^-2\n");
    12751787                _printf_("     Mean density: " << EarthMass[numlayers-1]/(4.0/3.0*PI*pow(r0,3.0)) << " kg.m^-3\n");
    12761788        }
    12771789
    1278         return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer,deg_layer_delete);
     1790        xDelete<doubletype>(r);
     1791        return new LoveVariables<doubletype>(EarthMass,g0,r0,nyi,starting_layer,deg_layer_delete,nstep);
    12791792
    12801793} /*}}}*/
    … …  
    12901803        bool        save_results;
    12911804        bool        complex_computation;
     1805        bool        quad_precision;
    12921806        bool        verbosecpu=false;
    12931807
    … …  
    12951809        doubletype  lovek, loveh, lovel, loveratio;
    12961810        IssmDouble pw_threshold, pw_test_h, pw_test_l,pw_test_k;
    1297 
    1298         LoveNumbers<doubletype>* Lovef=NULL;
    1299         LoveNumbers<doubletype>* Tidalf=NULL;
    13001811
    13011812        /* parallel computing */
    … …  
    13091820        IssmDouble* frequencies_fluid=NULL;
    13101821
    1311         LoveVariables* vars=NULL;
     1822        LoveVariables<doubletype>* vars=NULL;
    13121823
    13131824        /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
    … …  
    13321843        femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
    13331844        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
     1845        femmodel->parameters->FindParam(&quad_precision,LoveQuadPrecisionEnum);
    13341846        femmodel->parameters->FindParam(&pw_threshold,LovePostWidderThresholdEnum);
    13351847        if (istemporal) femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
    13361848
    1337         /*Initialize three love matrices: geoid, vertical and horizontal displacement*/
    1338         Lovef= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nfreq, matlitho);
    1339         Tidalf= new LoveNumbers<doubletype>(2,2,nfreq, matlitho);
    1340         Elastic= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,matlitho);
    1341         Fluid= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,matlitho);
    1342 
    1343         /*Initialize love kernels (real and imaginary parts): */
    1344         vars=love_init(femmodel,matlitho);
    1345 
     1849        Elastic= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,1,1,matlitho);
     1850        Fluid= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,1,1,matlitho);
    13461851        //distribute frequencies for parallel computation /*{{{*/
    13471852        int nt_local, nf_local, lower_row, upper_row;
    … …  
    13641869        if (lower_row==0) verbosecpu=true; //let only cpu1 be verbose
    13651870        if(VerboseSolution() && verbosecpu) _printf0_("   computing LOVE numbers\n");
     1871        vars=love_init<doubletype>(femmodel,matlitho,verbosecpu);
    13661872
    13671873        frequencies_local=xNewZeroInit<IssmDouble>(nf_local);
    13681874        for (int fr=0;fr<nf_local;fr++) frequencies_local[fr]=frequencies[lower_row+fr];
    13691875
    1370         Lovef_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nf_local, matlitho);
    1371         Tidalf_local= new LoveNumbers<doubletype>(2,2,nf_local, matlitho);
     1876        Lovef_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nf_local,lower_row,nfreq, matlitho);
     1877        Tidalf_local= new LoveNumbers<doubletype>(2,2,nf_local,lower_row,nfreq, matlitho);
    13721878
    13731879        /*}}}*/
    … …  
    13811887        // run elastic and fluid love numbers
    13821888        if(VerboseSolution() && verbosecpu) _printf_("     elastic\n");
     1889        EarthRheology<doubletype>(vars,frequencies_elastic,1,matlitho,femmodel);
    13831890        compute_love_numbers<doubletype>(Elastic, NULL, forcing_type, sh_nmax,frequencies_elastic, femmodel, matlitho, vars,verbosecpu);
    13841891
    13851892        if (nfreq>1){
     1893                EarthRheology<doubletype>(vars,frequencies_fluid,1,matlitho,femmodel);
    13861894                compute_love_numbers<doubletype>(Fluid, NULL, forcing_type, sh_nmax,frequencies_fluid, femmodel, matlitho, vars,verbosecpu);
    13871895                sh_cutoff=sh_nmax;
    … …  
    14021910        }
    14031911        else sh_cutoff=sh_nmax;
    1404                
    1405 
    1406 
    1407         //Take care of rotationnal feedback love numbers first, if relevant /*{{{*/
    1408         if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ // if forcing is surface loading and we have degree 2
    1409                 if(VerboseSolution() && verbosecpu) _printf_("     tidal\n");
    1410                 int tidal_forcing_type=9;
    1411                 compute_love_numbers<doubletype>(Tidalf_local, NULL,tidal_forcing_type=9, 2,frequencies_local, femmodel, matlitho, vars,verbosecpu);
    1412         }
    1413         /*}}}*/
    1414 
    1415         //Resume requested forcing_type
     1912
     1913        delete Fluid;
     1914
     1915        //Requested forcing_type
    14161916        if (nfreq>1){ // if we are not running just elastic love numbers
    14171917                if(VerboseSolution() && verbosecpu){
    … …  
    14201920                        else _printf_("     love\n");
    14211921                }
     1922                EarthRheology<doubletype>(vars,frequencies_local,nf_local,matlitho,femmodel);
    14221923                compute_love_numbers<doubletype>(Lovef_local, Elastic, forcing_type, sh_cutoff, frequencies_local, femmodel, matlitho, vars,verbosecpu);
    14231924        }
    14241925        else{
    14251926                Lovef_local->Copy(Elastic);
     1927        }
     1928        /*}}}*/
     1929
     1930        //Take care of rotationnal feedback love numbers, if relevant /*{{{*/
     1931        if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ // if forcing is surface loading and we have degree 2
     1932                if(VerboseSolution() && verbosecpu) _printf_("     tidal\n");
     1933                int tidal_forcing_type=9;
     1934                //no need to call EarthRheology, we already have the right one
     1935                compute_love_numbers<doubletype>(Tidalf_local, NULL,tidal_forcing_type=9, 2,frequencies_local, femmodel, matlitho, vars,verbosecpu);
    14261936        }
    14271937        /*}}}*/
    … …  
    14311941        if (istemporal && !complex_computation){
    14321942                /*Initialize*/
    1433                 doubletype*  pmtf_colineart=NULL;
    1434                 doubletype*  pmtf_orthot=NULL;
    1435                 LoveNumbers<doubletype>* Lovet=NULL;
    1436                 LoveNumbers<doubletype>* Tidalt=NULL;
    14371943                /*Downcast arrays to be exported in parameters*/
    1438                 LoveNumbers<IssmDouble>* LovetDouble=NULL;
    1439                 LoveNumbers<IssmDouble>* TidaltDouble=NULL;
    14401944                IssmDouble*  pmtf_colineartDouble=NULL;
    14411945                IssmDouble*  pmtf_orthotDouble=NULL;
     1946
    14421947                /* parallel computing */
     1948                LoveNumbers<IssmDouble>* LovefDouble_local=NULL;
     1949                LoveNumbers<IssmDouble>* LovetDouble_local=NULL;
     1950                LoveNumbers<IssmDouble>* TidaltDouble_local=NULL;
     1951                IssmDouble*  pmtf_colineartDouble_local=NULL;
     1952                IssmDouble*  pmtf_orthotDouble_local=NULL;
     1953
    14431954                doubletype*  pmtf_colineart_local=NULL;
    14441955                doubletype*  pmtf_orthot_local=NULL;   
    … …  
    14461957                LoveNumbers<doubletype>* Tidalt_local=NULL;     
    14471958
    1448                 Lovet= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nt,matlitho);
    1449                 Tidalt= new LoveNumbers<doubletype>(2,2,nt,matlitho);
    1450                 pmtf_colineart=xNewZeroInit<doubletype>(3*nt);
    1451                 pmtf_orthot=xNewZeroInit<doubletype>(3*nt);     
    1452 
    1453                 Lovet_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nt_local,matlitho);
    1454                 Tidalt_local= new LoveNumbers<doubletype>(2,2,nt_local,matlitho);       
     1959                Lovet_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nt_local,lower_row/2/NTit,nt,matlitho);
     1960                Tidalt_local= new LoveNumbers<doubletype>(2,2,nt_local,lower_row/2/NTit,nt,matlitho);   
    14551961                pmtf_colineart_local=xNewZeroInit<doubletype>(3*nt_local);
    14561962                pmtf_orthot_local=xNewZeroInit<doubletype>(3*nt_local);
    … …  
    14581964                love_freq_to_temporal<doubletype>(Lovet_local,Tidalt_local,pmtf_colineart_local,pmtf_orthot_local,Lovef_local,Tidalf_local,frequencies_local,femmodel,verbosecpu);
    14591965
    1460                 /* MPI Gather */ /*{{{*/
    1461                 Lovef->LoveMPI_Gather(Lovef_local, lower_row);
    1462                 if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){
    1463                         Tidalf->LoveMPI_Gather(Tidalf_local, lower_row);               
    1464                 }
    1465                 Lovet->LoveMPI_Gather(Lovet_local, lower_row/2/NTit);
    1466                 Tidalt->LoveMPI_Gather(Tidalt_local, lower_row/2/NTit);
     1966
     1967                if(VerboseSolution() && verbosecpu) _printf_("   Assembling parralel vectors...");
     1968
     1969                //delete Lovef_local;
     1970                delete Tidalf_local;
     1971                //Lovet
     1972                LovetDouble_local= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nt_local,lower_row/2/NTit,nt,matlitho);
     1973                Lovet_local->DownCastToDouble(LovetDouble_local);
     1974                delete Lovet_local;
     1975                LovetDouble_local->Broadcast();
     1976
     1977                //Lovef
     1978                LovefDouble_local= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nf_local,lower_row,nfreq,matlitho);
     1979                Lovef_local->DownCastToDouble(LovefDouble_local);
     1980                delete Lovef_local;
     1981                LovefDouble_local->Broadcast();
     1982
     1983                if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){                     
     1984                        TidaltDouble_local= new LoveNumbers<IssmDouble>(2,2,nt_local,lower_row,nfreq,matlitho);
     1985                        Tidalt_local->DownCastToDouble(TidaltDouble_local);
     1986                        delete Tidalt_local;
     1987                        TidaltDouble_local->Broadcast();
     1988                }
     1989
    14671990                //pmtf:
     1991                pmtf_colineartDouble_local=xNew<IssmDouble>(nt_local);
     1992                pmtf_orthotDouble_local=xNew<IssmDouble>(nt_local);
     1993                /*Downcast*/ /*{{{*/
     1994                for(int i=0;i<nt_local;i++){
     1995                        pmtf_colineartDouble_local[i]=DownCastVarToDouble<doubletype>(pmtf_colineart_local[i*3+2]);
     1996                        pmtf_orthotDouble_local[i]=DownCastVarToDouble<doubletype>(pmtf_orthot_local[i*3+2]);
     1997                }
     1998                /*}}}*/
     1999                xDelete<doubletype>(pmtf_colineart_local);
     2000                xDelete<doubletype>(pmtf_orthot_local);
     2001                pmtf_colineartDouble=xNew<IssmDouble>(nt);
     2002                pmtf_orthotDouble=xNew<IssmDouble>(nt);
     2003
    14682004                int* recvcounts=xNew<int>(IssmComm::GetSize());
    14692005                int* displs=xNew<int>(IssmComm::GetSize());
    14702006                int  rc;
    14712007                int  offset;
    1472                 rc=3*nt_local;
    1473                 offset=3*lower_row/2/NTit;
     2008                rc=nt_local;
     2009                offset=lower_row/2/NTit;
    14742010                ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
    14752011                ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
    1476                 ISSM_MPI_Allgatherv(pmtf_colineart_local, rc, ISSM_MPI_DOUBLE, pmtf_colineart, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    1477                 ISSM_MPI_Allgatherv(pmtf_orthot_local, rc, ISSM_MPI_DOUBLE, pmtf_orthot, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     2012                ISSM_MPI_Allgatherv(pmtf_colineartDouble_local, rc, ISSM_MPI_DOUBLE, pmtf_colineartDouble, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     2013                ISSM_MPI_Allgatherv(pmtf_orthotDouble_local, rc, ISSM_MPI_DOUBLE, pmtf_orthotDouble, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    14782014                xDelete<int>(recvcounts);
    14792015                xDelete<int>(displs);
     2016
     2017                xDelete<IssmDouble>(pmtf_colineartDouble_local);
     2018                xDelete<IssmDouble>(pmtf_orthotDouble_local);
    14802019                /*}}}*/
    14812020
    1482                 /*Downcast and add into parameters:*/ /*{{{*/
    1483                 LovetDouble= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nt,matlitho);
    1484                 TidaltDouble= new LoveNumbers<IssmDouble>(2,2,nt,matlitho);     
    1485 
    1486                 pmtf_colineartDouble=xNew<IssmDouble>(nt);
    1487                 pmtf_orthotDouble=xNew<IssmDouble>(nt);
    1488 
    1489                 Lovet->DownCastToDouble(LovetDouble);
    1490                 Tidalt->DownCastToDouble(TidaltDouble);
    1491                 for(int i=0;i<nt;i++){
    1492                         pmtf_colineartDouble[i]=std::real(pmtf_colineart[2*nt+i]);
    1493                         pmtf_orthotDouble[i]=std::real(pmtf_orthot[2*nt+i]);
    1494                 }
    1495 
     2021                if(VerboseSolution() && verbosecpu) _printf_("done\n");
     2022                if(VerboseSolution() && verbosecpu) _printf_("   saving results\n");
     2023
     2024                /* Add to parameters */ /*{{{*/
    14962025                if(forcing_type==9){ //tidal loading
    1497                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1));
    1498                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1));
    1499                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1));
    1500                 }
     2026                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovetDouble_local->H,(sh_nmax+1)*nt,1));
     2027                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovetDouble_local->K,(sh_nmax+1)*nt,1));
     2028                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovetDouble_local->L,(sh_nmax+1)*nt,1));
     2029                                }
    15012030                else if(forcing_type==11){ //surface loading
    1502                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1));
    1503                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1));
    1504                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1));
    1505                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,TidaltDouble->H,3*nt,1));
    1506                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,TidaltDouble->K,3*nt,1));
    1507                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,TidaltDouble->L,3*nt,1));
     2031                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovetDouble_local->H,(sh_nmax+1)*nt,1));
     2032                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovetDouble_local->K,(sh_nmax+1)*nt,1));
     2033                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovetDouble_local->L,(sh_nmax+1)*nt,1));
     2034                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,TidaltDouble_local->H,3*nt,1));
     2035                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,TidaltDouble_local->K,3*nt,1));
     2036                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,TidaltDouble_local->L,3*nt,1));
    15082037                        femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,pmtf_colineartDouble,nt,1));
    15092038                        femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,pmtf_orthotDouble,nt,1));
    15102039                }
     2040                /*}}}*/
     2041       
     2042                /*Add into external results*/
     2043                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKtEnum,LovetDouble_local->K,nt,sh_nmax+1,0,0));
     2044                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHtEnum,LovetDouble_local->H,nt,sh_nmax+1,0,0));
     2045                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLtEnum,LovetDouble_local->L,nt,sh_nmax+1,0,0));
     2046                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKfEnum,LovefDouble_local->K,nfreq,sh_nmax+1,0,0));
     2047                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHfEnum,LovefDouble_local->H,nfreq,sh_nmax+1,0,0));
     2048                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLfEnum,LovefDouble_local->L,nfreq,sh_nmax+1,0,0));
     2049
     2050                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveTidalKtEnum,TidaltDouble_local->K,nt,3,0,0));
     2051                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveTidalHtEnum,TidaltDouble_local->H,nt,3,0,0));
     2052                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveTidalLtEnum,TidaltDouble_local->L,nt,3,0,0));
     2053                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LovePMTF1tEnum,pmtf_colineartDouble,nt,1,0,0));
     2054                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LovePMTF2tEnum,pmtf_orthotDouble,nt,1,0,0));
     2055                /*Only when love_kernels is on*/
     2056                if (love_kernels==1) {
     2057                //      femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsEnum,LovefDouble_local->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0));
     2058                }
    15112059
    15122060                xDelete<IssmDouble>(pmtf_colineartDouble);
    15132061                xDelete<IssmDouble>(pmtf_orthotDouble);
    1514                 /*}}}*/
    1515        
    1516                 /*Add into external results, no need to downcast, we can handle complexes/quad precision: */
    1517                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKtEnum,Lovet->K,nt,sh_nmax+1,0,0));
    1518                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHtEnum,Lovet->H,nt,sh_nmax+1,0,0));
    1519                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLtEnum,Lovet->L,nt,sh_nmax+1,0,0));
    1520 
    1521                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalKtEnum,Tidalt->K,nt,3,0,0));
    1522                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalHtEnum,Tidalt->H,nt,3,0,0));
    1523                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalLtEnum,Tidalt->L,nt,3,0,0));
    1524                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LovePMTF1tEnum,pmtf_colineart,nt,3,0,0));
    1525                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LovePMTF2tEnum,pmtf_orthot,nt,3,0,0));
    1526 
    1527                 delete Lovet;
    1528                 delete Tidalt;
    1529                 delete Lovet_local;
    1530                 delete Tidalt_local;
    1531                 delete LovetDouble;
    1532                 delete TidaltDouble;
    1533 
    1534                 xDelete<doubletype>(pmtf_colineart);
    1535                 xDelete<doubletype>(pmtf_orthot);
     2062                delete LovetDouble_local;
     2063                delete TidaltDouble_local;
     2064
    15362065        }
    15372066        else{
    15382067                LoveNumbers<IssmDouble>* LovefDouble=NULL;
    1539                 LovefDouble= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nfreq,matlitho);
     2068                LoveNumbers<IssmDouble>* LovefDouble_local=NULL;
     2069                LovefDouble= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nfreq,lower_row,nfreq,matlitho);
     2070                LovefDouble_local= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nf_local,lower_row,nfreq,matlitho);
     2071
     2072                LoveNumbers<IssmDouble>* TidalfDouble=NULL;
     2073                LoveNumbers<IssmDouble>* TidalfDouble_local=NULL;
     2074                TidalfDouble= new LoveNumbers<IssmDouble>(2,2,nfreq,lower_row,nfreq,matlitho);
     2075                TidalfDouble_local= new LoveNumbers<IssmDouble>(2,2,nf_local,lower_row,nfreq,matlitho);
     2076
     2077                Lovef_local->DownCastToDouble(LovefDouble_local);
     2078                Tidalf_local->DownCastToDouble(TidalfDouble_local);
    15402079
    15412080                /*MPI_Gather*/
    15422081                if (nfreq>1){
    1543                         Lovef->LoveMPI_Gather(Lovef_local, lower_row);
     2082                        LovefDouble->LoveMPI_Gather(LovefDouble_local, lower_row);
    15442083                        if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){
    1545                                 Tidalf->LoveMPI_Gather(Tidalf_local, lower_row);               
     2084                                TidalfDouble->LoveMPI_Gather(TidalfDouble_local, lower_row);           
    15462085                        }
    15472086                }
    15482087                else{
    1549                         Lovef->Copy(Elastic);
    1550                         Tidalf->Copy(Tidalf_local);
     2088                        Elastic->DownCastToDouble(LovefDouble);
     2089                        Tidalf_local->DownCastToDouble(TidalfDouble);
    15512090                }
    15522091
    15532092                /*Add into parameters:*/
    1554                 Lovef->DownCastToDouble(LovefDouble);
    15552093                if(forcing_type==9){ //tidal loading
    15562094                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovefDouble->H,(sh_nmax+1)*nfreq,1));
    … …  
    15632101                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovefDouble->L,(sh_nmax+1)*nfreq,1));
    15642102                }
     2103
     2104                /*Add into external results:*/
     2105                if (complex_computation){
     2106                        //FIXME: complex external result not supported yet
     2107                        //femmodel->results->AddObject(new GenericExternalResult<IssmComplex*>(femmodel->results->Size()+1,LoveKfEnum,Lovef->K,nfreq,sh_nmax+1,0,0));
     2108                        //femmodel->results->AddObject(new GenericExternalResult<IssmComplex*>(femmodel->results->Size()+1,LoveHfEnum,Lovef->H,nfreq,sh_nmax+1,0,0));
     2109                        //femmodel->results->AddObject(new GenericExternalResult<IssmComplex*>(femmodel->results->Size()+1,LoveLfEnum,Lovef->L,nfreq,sh_nmax+1,0,0));
     2110                        ///*Only when love_kernels is on*/
     2111                        //if (love_kernels==1) {
     2112                        //      femmodel->results->AddObject(new GenericExternalResult<IssmComplex*>(femmodel->results->Size()+1,LoveKernelsEnum,Lovef->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0));
     2113                        //}
     2114                }
     2115                else{
     2116                        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKfEnum,LovefDouble->K,nfreq,sh_nmax+1,0,0));
     2117                        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHfEnum,LovefDouble->H,nfreq,sh_nmax+1,0,0));
     2118                        femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLfEnum,LovefDouble->L,nfreq,sh_nmax+1,0,0));
     2119                        /*Only when love_kernels is on*/
     2120                        if (love_kernels==1) {
     2121                                femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsEnum,LovefDouble->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0));
     2122                        }
     2123                }
     2124
     2125
     2126                delete Lovef_local;
    15652127                delete LovefDouble;
    1566         }
    1567 
    1568         /*Add into external results:*/
    1569         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKfEnum,Lovef->K,nfreq,sh_nmax+1,0,0));
    1570         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHfEnum,Lovef->H,nfreq,sh_nmax+1,0,0));
    1571         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLfEnum,Lovef->L,nfreq,sh_nmax+1,0,0));
    1572         /*Only when love_kernels is on*/
    1573         if (love_kernels==1) {
    1574                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,Lovef->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0));
     2128                delete LovefDouble_local;
     2129                delete Tidalf_local;
     2130                delete TidalfDouble;
     2131                delete TidalfDouble_local;
    15752132        }
    15762133        /*Free resources:*/
    … …  
    15782135        xDelete<IssmDouble>(frequencies_local);
    15792136        xDelete<IssmDouble>(frequencies_elastic);
    1580         delete Lovef;
    1581         delete Lovef_local;
    1582         delete Tidalf;
    1583         delete Tidalf_local;
     2137
    15842138        delete Elastic;
    1585         /* Legacy for fortran core, to be removed after complete validation */
    1586 
    1587         //IssmDouble g0,r0;
    1588         //femmodel->parameters->FindParam(&g0,LoveG0Enum);
    1589         //femmodel->parameters->FindParam(&r0,LoveR0Enum);
    1590         //IssmDouble* LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1591         //IssmDouble* LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1592         //IssmDouble* LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1593         //IssmDouble* LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1594         //IssmDouble* LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1595         //IssmDouble* LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
    1596 
    1597         /*Initialize love kernels (real and imaginary parts): */
    1598         //IssmDouble* LoveKernelsr= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
    1599         //IssmDouble* LoveKernelsi= xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
    1600 
    1601         /*call the main module: */
    1602         //if (false){
    1603         //FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHr,LoveLr,LoveLi,LoveKernelsr,LoveKernelsi,  //output
    1604         //              nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
    1605         //              matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
    1606         //              matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->rheologymodel, matlitho->issolid //matlitho inputs
    1607         //              );
    1608 
    1609         /*Only when love_kernels is on*/
    1610         //if (love_kernels==1) {
    1611                 //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernelsr,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    1612                 //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    1613         //}
    1614 
    1615         /*Add love matrices to results:*/
    1616         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKtEnum,LoveKr,sh_nmax+1,nfreq,0,0));
    1617         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHtEnum,LoveHr,sh_nmax+1,nfreq,0,0));
    1618         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLtEnum,LoveLr,sh_nmax+1,nfreq,0,0));
    1619         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKfEnum,LoveKi,sh_nmax+1,nfreq,0,0));
    1620         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHfEnum,LoveHi,sh_nmax+1,nfreq,0,0));
    1621         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLfEnum,LoveLi,sh_nmax+1,nfreq,0,0));
    1622 
    1623         //xDelete<IssmDouble>(LoveKr);
    1624         //xDelete<IssmDouble>(LoveHr);
    1625         //xDelete<IssmDouble>(LoveLr);
    1626         //xDelete<IssmDouble>(LoveKernelsr);
    1627         //xDelete<IssmDouble>(LoveKi);
    1628         //xDelete<IssmDouble>(LoveHi);
    1629         //xDelete<IssmDouble>(LoveLi);
    1630         //xDelete<IssmDouble>(LoveKernelsi);
    16312139
    16322140} /*}}}*/
    … …  
    16342142/*cores and template instantiations:*/
    16352143/*template instantiations :{{{*/
     2144// IssmDouble
    16362145template void love_core_template<IssmDouble>(FemModel* femmodel);
    1637 template void love_core_template<IssmComplex>(FemModel* femmodel);
    1638 template void        fill_yi_prefactor<IssmDouble>(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
    1639 template void        fill_yi_prefactor<IssmComplex>(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
    1640 template void        GetEarthRheology<IssmDouble>(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,  Matlitho* matlitho);
    1641 template void        GetEarthRheology<IssmComplex>(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega,  Matlitho* matlitho);
    1642 template void        yi_boundary_conditions<IssmDouble>(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type);
    1643 template void        yi_boundary_conditions<IssmComplex>(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type);
    1644 template void        yi_derivatives<IssmDouble>(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1645 template void        yi_derivatives<IssmComplex>(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1646 template void        propagate_yi_RK2<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1647 template void        propagate_yi_RK2<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1648 template void        propagate_yi_euler<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1649 template void        propagate_yi_euler<IssmComplex>(IssmComplex* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    1650 template void        Innersphere_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
    1651 template void        Innersphere_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
    1652 template void        build_yi_system<IssmDouble>(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
    1653 template void        build_yi_system<IssmComplex>(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
    1654 template void        solve_yi_system<IssmDouble>(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* frequencies, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu);
    1655 template void        solve_yi_system<IssmComplex>(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmDouble* frequencies, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars,bool verbosecpu);
    1656 template void        compute_love_numbers<IssmDouble>(LoveNumbers<IssmDouble>* Lovef, LoveNumbers<IssmDouble>* Elastic, int forcing_type, int sh_cutoff,IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu);
    1657 template void        compute_love_numbers<IssmComplex>(LoveNumbers<IssmComplex>* Lovef, LoveNumbers<IssmComplex>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars, bool verbosecpu);
     2146template LoveVariables<IssmDouble>*     love_init<IssmDouble>(FemModel* femmodel, Matlitho* matlitho,bool verbosecpu);
     2147template void        fill_yi_prefactor<IssmDouble>(IssmDouble* yi_prefactor, int* pdeg, IssmDouble* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2148template void        GetEarthRheology<IssmDouble>(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,  Matlitho* matlitho, FemModel* femmodel);
     2149template IssmDouble     GetGravity<IssmDouble>(IssmDouble r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmDouble>* vars);
     2150template void        yi_boundary_conditions<IssmDouble>(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmDouble>* vars, int forcing_type);
     2151template void        yi_derivatives<IssmDouble>(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2152template void        propagate_yi_RK2<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2153template void        propagate_yi_RK4<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2154template void        propagate_yi_euler<IssmDouble>(IssmDouble* y, IssmDouble xmin, IssmDouble xmax, int layer_index, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2155template void        Innersphere_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2156template void        Coremantle_boundaryconditions<IssmDouble>(IssmDouble* yi, int layer_index, int deg, IssmDouble omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars);
     2157template void        build_yi_system<IssmDouble>(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmDouble>* vars);
     2158template void        solve_yi_system<IssmDouble>(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* frequencies, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars,bool verbosecpu);
     2159template void        compute_love_numbers<IssmDouble>(LoveNumbers<IssmDouble>* Lovef, LoveNumbers<IssmDouble>* Elastic, int forcing_type, int sh_cutoff,IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmDouble>* vars, bool verbosecpu);
    16582160template IssmDouble  factorial<IssmDouble>(int n);
    16592161template IssmDouble* postwidder_coef<IssmDouble>(int NTit);
    16602162template IssmDouble  n_C_r<IssmDouble>(int n, int r);
    16612163template void         postwidder_transform<IssmDouble>(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int sh_nmax,int NTit, IssmDouble* xi, FemModel* femmodel);
     2164template void        EarthRheology<IssmDouble>(LoveVariables<IssmDouble>* vars, IssmDouble* frequencies, int nfreq,  Matlitho* matlitho, FemModel* femmodel);
     2165template IssmDouble HypergeomTableLookup(IssmDouble z1, IssmDouble alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha);
     2166
     2167//IssmComplex
     2168template void love_core_template<IssmComplex>(FemModel* femmodel);
     2169template LoveVariables<IssmComplex>*    love_init<IssmComplex>(FemModel* femmodel, Matlitho* matlitho,bool verboscpu);
     2170template void        fill_yi_prefactor<IssmComplex>(IssmComplex* yi_prefactor, int* pdeg, IssmComplex* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2171template void        GetEarthRheology<IssmComplex>(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega,  Matlitho* matlitho, FemModel* femmodel);
     2172template IssmComplex    GetGravity<IssmComplex>(IssmComplex r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmComplex>* vars);
     2173template void        yi_boundary_conditions<IssmComplex>(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmComplex>* vars, int forcing_type);
     2174template void        yi_derivatives<IssmComplex>(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2175template void        propagate_yi_RK2<IssmComplex>(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2176template void        propagate_yi_RK4<IssmComplex>(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2177template void        propagate_yi_euler<IssmComplex>(IssmComplex* y, IssmComplex xmin, IssmComplex xmax, int layer_index, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2178template void        Innersphere_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2179template void        Coremantle_boundaryconditions<IssmComplex>(IssmComplex* yi, int layer_index, int deg, IssmComplex omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars);
     2180template void        build_yi_system<IssmComplex>(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<IssmComplex>* vars);
     2181template void        solve_yi_system<IssmComplex>(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmDouble* frequencies, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars,bool verbosecpu);
     2182template void        compute_love_numbers<IssmComplex>(LoveNumbers<IssmComplex>* Lovef, LoveNumbers<IssmComplex>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables<IssmComplex>* vars, bool verbosecpu);
     2183template void        EarthRheology<IssmComplex>(LoveVariables<IssmComplex>* vars, IssmDouble* frequencies, int nfreq,  Matlitho* matlitho, FemModel* femmodel);
     2184template IssmComplex HypergeomTableLookup(IssmComplex z1, IssmComplex alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha);
     2185
     2186//__float128
     2187#ifdef _HAVE_MPLAPACK_
     2188template void love_core_template<__float128>(FemModel* femmodel);
     2189template LoveVariables<__float128>*     love_init<__float128>(FemModel* femmodel, Matlitho* matlitho, bool verbosecpu);
     2190template void        fill_yi_prefactor<__float128>(__float128* yi_prefactor, int* pdeg, __float128* pomega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2191template void        GetEarthRheology<__float128>(__float128* pla, __float128* pmu, int layer_index, __float128 omega,  Matlitho* matlitho, FemModel* femmodel);
     2192template __float128     GetGravity<__float128>(__float128 r2, int layer_index, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars);
     2193template void        yi_boundary_conditions<__float128>(__float128* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars, int forcing_type);
     2194template void        yi_derivatives<__float128>(__float128* dydx, __float128* y, int layer_index, int n, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2195template void        propagate_yi_RK2<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2196template void        propagate_yi_RK4<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2197template void        propagate_yi_euler<__float128>(__float128* y, __float128 xmin, __float128 xmax, int layer_index, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2198template void        Innersphere_boundaryconditions<__float128>(__float128* yi, int layer_index, int deg, __float128 omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2199template void        Coremantle_boundaryconditions<__float128>(__float128* yi, int layer_index, int deg, __float128 omega, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars);
     2200template void        build_yi_system<__float128>(__float128* yi, int deg, __float128 omega, __float128* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables<__float128>* vars);
     2201template void        solve_yi_system<__float128>(__float128* loveh, __float128* lovel, __float128* lovek, int deg, __float128 omega, IssmDouble* frequencies, __float128* yi, __float128* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars,bool verbosecpu);
     2202template void        compute_love_numbers<__float128>(LoveNumbers<__float128>* Lovef, LoveNumbers<__float128>* Elastic, int forcing_type, int sh_cutoff, IssmDouble* frequencies, FemModel* femmodel, Matlitho* matlitho, LoveVariables<__float128>* vars, bool verbosecpu);
     2203template __float128  factorial<__float128>(int n);
     2204template __float128* postwidder_coef<__float128>(int NTit);
     2205template __float128  n_C_r<__float128>(int n, int r);
     2206template void        postwidder_transform<__float128>(__float128* Lovet, __float128* Lovef,int d, int t, int sh_nmax,int NTit, __float128* xi, FemModel* femmodel);
     2207template void        EarthRheology<__float128>(LoveVariables<__float128>* vars, IssmDouble* frequencies, int nfreq,  Matlitho* matlitho, FemModel* femmodel);
     2208template __float128 HypergeomTableLookup(__float128 z1, __float128 alpha, IssmDouble* h1, IssmDouble* z, int nz, int nalpha);
     2209#endif
    16622210
    16632211/*}}}*/
    16642212void           love_core(FemModel* femmodel){ /*{{{*/
    1665 
    16662213        bool        complex_computation;
     2214        bool        quad_precision;
    16672215
    16682216        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
     2217        femmodel->parameters->FindParam(&quad_precision,LoveQuadPrecisionEnum);
    16692218
    16702219        if(complex_computation) love_core_template<IssmComplex>(femmodel);
     2220#ifdef _HAVE_MPLAPACK_
     2221        else if(quad_precision) love_core_template<__float128>(femmodel);
     2222#endif
    16712223        else                    love_core_template<IssmDouble>(femmodel);
    16722224

  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r27131 r27308  
    2222void TransferForcing(FemModel* femmodel,int forcingenum);
    2323void TransferSealevel(FemModel* femmodel,int forcingenum);
    24 bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
     24bool slcconvergence(IssmDouble* RSLg,IssmDouble* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs, IssmDouble totaloceanarea,FemModel* femmodel);
    2525IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble totaloceanarea);
    2626void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom, bool computefuture);
    … …  
    5252        femmodel->SetCurrentConfiguration(SealevelchangeAnalysisEnum);
    5353
     54        /*Run coupler input transfer:*/
     55        couplerinput_core(femmodel);
     56
    5457        /*run geometry core: */
    5558        slgeom=sealevelchange_geometry(femmodel);
    … …  
    5760        /*any external forcings?:*/
    5861        solidearthexternal_core(femmodel);
    59 
    60         /*Run coupler input transfer:*/
    61         couplerinput_core(femmodel);
    6262
    6363        /*Run geodetic:*/
    … …  
    184184        /*if we are carrying loads but are not yet computing grd core, accumulate them and skip
    185185         * the rest: */
    186         if (count<frequency){
    187                 count++;
    188                 femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum);
     186        if (count!=frequency){
     187                if (count>frequency){
     188                        count=1;
     189                        femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum);
     190                }
    189191                return;
    190192        }
     193       
    191194
    192195        /*Basins are supposed to accumulate loads and hand them over to the Earth
    … …  
    208211                TransferForcing(femmodel,DeltaBottomPressureEnum);
    209212                TransferForcing(femmodel,DeltaTwsEnum);
     213                TransferForcing(femmodel,MaskOceanLevelsetEnum);
     214                TransferForcing(femmodel,MaskIceLevelsetEnum);
    210215
    211216                /*transfer external forcings back to Earth:*/
    … …  
    228233
    229234        GrdLoads*              loads=NULL;
    230         Vector<IssmDouble>*    oldsealevelloads=NULL;
     235        IssmDouble*    oldsealevelloads=NULL;
    231236        Vector<IssmDouble>*    oceanareas=NULL;
    232237        IssmDouble             totaloceanarea;
    … …  
    249254        int  grdmodel;
    250255        int  sealevelloading=0;
     256        bool sal=false;
    251257        bool viscous=false;
    252258        bool rotation=false;
    … …  
    256262
    257263        /*}}}*/
    258 
    259         /*Verbose: */
    260         if(VerboseSolution()) _printf0_("         computing GRD patterns\n");
    261264
    262265        /*retrieve parameters:{{{*/
    … …  
    290293        }
    291294
     295        /*Verbose: */
     296        if(VerboseSolution()) _printf0_("         computing GRD patterns\n");
     297
     298
    292299        /*retrieve parameters: {{{*/
    293300        femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
    … …  
    298305        femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
    299306        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     307        femmodel->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum);
    300308        /*}}}*/
    301309
    … …  
    336344                PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom,computefuture=false);
    337345
    338                 oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads);
     346                oldsealevelloads=xNewZeroInit<IssmDouble>(nel);
     347                if (loads->sealevelloads){
     348                        xMemCpy<IssmDouble>(oldsealevelloads,loads->sealevelloads,nel);
     349                }
     350
    339351
    340352                /*convolve load and sealevel loads on oceans:*/
     353                loads->Combineloads(nel,slgeom); //This combines loads and sealevelloads into a single vector
    341354                for(Object* & object : femmodel->elements->objects){
    342355                        Element* element = xDynamicCast<Element*>(object);
    … …  
    367380                loads->BroadcastSealevelLoads();
    368381
     382                if (!sal) xDelete<IssmDouble>(oldsealevelloads); break;
     383
    369384                //convergence?
    370                 if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs)){
    371                         delete oldsealevelloads; break;
     385                if(slcconvergence(loads->sealevelloads,oldsealevelloads,eps_rel,eps_abs,totaloceanarea,femmodel)){
     386                        xDelete<IssmDouble>(oldsealevelloads); break;
    372387                }
    373388
    374389                //early return?
    375390                if(iterations>=max_nonlinear_iterations){
    376                         delete oldsealevelloads; break;
     391                        xDelete<IssmDouble>(oldsealevelloads); break;
    377392                }
    378393                iterations++;
    379                 delete oldsealevelloads;
     394                xDelete<IssmDouble>(oldsealevelloads);
    380395        }
    381396
    … …  
    388403
    389404        /*convolve loads and sea level loads to get the deformation:*/
     405        loads->Combineloads(nel,slgeom); //This combines loads and sealevelloads into a single vector
    390406        for(Object* & object : femmodel->elements->objects){
    391407                Element* element = xDynamicCast<Element*>(object);
    … …  
    512528        int iscoupling;
    513529        int horiz=0;
     530        int count, frequency;
    514531
    515532        /*retrieve more parameters:*/
    516533        femmodel->parameters->FindParam(&iscoupling,IsSlcCouplingEnum);
    517534        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    518 
     535        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
     536        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
     537
     538        count++;
     539        femmodel->parameters->SetParam(count,SealevelchangeRunCountEnum);
     540       
    519541        if(iscoupling){
    520542                /*transfer sea level back to ice caps:*/
    … …  
    526548                }
    527549        }
     550
    528551}; /*}}}*/
    529552void              ivins_deformation_core(FemModel* femmodel){ /*{{{*/
    … …  
    588611        IssmDouble* areae  = NULL;
    589612        int  nel;
    590         int* lids;
     613        int* lids=NULL;
     614        int* n_activevertices=NULL;
    591615        int  grdmodel=0;
    592616
    … …  
    596620
    597621        /*early return?:*/
    598         if(grdmodel==IvinsEnum) return;
     622        if(grdmodel!=ElasticEnum) return;
    599623
    600624        /*Verbose: */
    … …  
    607631        /*Compute element ids, used to speed up computations in convolution phase:{{{*/
    608632        lids=xNew<int>(femmodel->vertices->Size());
     633        n_activevertices = xNew<int>(nel);
     634        //initialize lids to -1, vertex count to 3
     635        for (int v=0; v<femmodel->vertices->Size();v++) lids[v]=-1;
     636        for (int e=0; e<nel;e++) n_activevertices[e]=3;
    609637
    610638        for(Object* & object : femmodel->elements->objects){
    611639                Element*   element=xDynamicCast<Element*>(object);
    612640                for(int i=0;i<3;i++){
     641                        // if lids where we are looking points to an element id (.i.e. not -1) then we are about to claim that element's vertex
     642                        // and need to lower the number of vertices it is in charge of
     643                        if (lids[element->vertices[i]->lid] !=-1){
     644                                n_activevertices[lids[element->vertices[i]->lid]]-=1;
     645                        }
    613646                        lids[element->vertices[i]->lid]=element->lid;
    614647                }
    … …  
    620653        for(Object* & object : femmodel->elements->objects){
    621654                Element*   element=xDynamicCast<Element*>(object);
    622                 element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids);
     655                element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids,n_activevertices);
    623656        }
    624657
    … …  
    639672        xDelete<IssmDouble>(zze);
    640673        xDelete<IssmDouble>(areae);
    641         xDelete(lids);
     674        xDelete<int>(lids);
     675        xDelete<int>(n_activevertices);
    642676
    643677        return;
    … …  
    658692        int  grdmodel=0;
    659693        int isgrd=0;
     694        int count, frequency;
    660695        SealevelGeometry* slgeom=NULL;
    661696
    … …  
    663698        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    664699        femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
     700        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
     701        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    665702        if(grdmodel!=ElasticEnum || !isgrd) return NULL;
     703        if(count!=frequency)return NULL;
    666704
    667705        /*retrieve parameters:*/
    … …  
    722760        int isgrd=0;
    723761        int horiz=0;
     762        int count, frequency;
    724763
    725764        /*early return?:*/
    … …  
    727766        femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
    728767        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     768        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
     769        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    729770        if(grdmodel!=ElasticEnum || !isgrd) return;
     771        if(count!=frequency)return;
    730772
    731773        for (int l=0;l<SLGEOM_NUMLOADS;l++){
    … …  
    738780
    739781/*subroutines:*/
    740 bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    741 
     782bool slcconvergence(IssmDouble* RSLg,IssmDouble* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs, IssmDouble totaloceanarea, FemModel* femmodel){ /*{{{*/
     783
     784        int nel;
    742785        bool converged=true;
    743         IssmDouble ndS,nS;
    744         Vector<IssmDouble> *dRSLg    = NULL;
     786        IssmDouble ndS,nS, nS_old;
     787        IssmDouble* dRSLg    = NULL;
     788        IssmDouble rho_water =0;
     789
     790        femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     791        femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    745792
    746793        //compute norm(du) and norm(u) if requested
    747         dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
    748         ndS=dRSLg->Norm(NORM_TWO);
     794        dRSLg=xNewZeroInit<IssmDouble>(nel);
     795
     796        ndS=0;
     797        nS=0;
     798        nS_old=0;
     799
     800        for (int e=0;e<nel;e++){
     801                dRSLg[e]=(RSLg[e]-RSLg_old[e])/rho_water/totaloceanarea;
     802                ndS+=pow(dRSLg[e],2.0);
     803                nS+=pow(RSLg[e]/rho_water/totaloceanarea,2.0);
     804                nS_old+=pow(RSLg_old[e]/rho_water/totaloceanarea,2.0);
     805        }
     806       
    749807
    750808        if (xIsNan<IssmDouble>(ndS)){
    751                 _error_("convergence criterion is NaN (RSL_old=" << RSLg_old->Norm(NORM_TWO) << " RSL=" << RSLg->Norm(NORM_TWO) << ")");
     809                _error_("convergence criterion is NaN (RSL_old=" << nS_old << " RSL=" << nS << ")");
    752810        }
    753811
    754812        if(!xIsNan<IssmDouble>(eps_rel)){
    755                 nS=RSLg_old->Norm(NORM_TWO);
    756                 if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN! (check the initial RSL)");
     813                if (xIsNan<IssmDouble>(nS_old)) _error_("convergence criterion is NaN! (check the initial RSL)");
    757814        }
    758815
    759816        //clean up
    760         delete dRSLg;
     817        xDelete<IssmDouble>(dRSLg);
    761818
    762819        //print
    … …  
    10481105        IssmDouble*  forcing=NULL;
    10491106        Vector<IssmDouble>* forcingglobal=NULL;
     1107        IssmDouble* transfercount=NULL;
    10501108        int*         nvs=NULL;
    10511109
    … …  
    10881146                nv=femmodel->vertices->NumberOfVertices();
    10891147                existforcing=reCast<int>(femmodel->inputs->Exist(forcingenum));
    1090                 if(existforcing)GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
     1148                if(existforcing){
     1149                        GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
     1150                        GetVectorFromInputsx(&transfercount,femmodel,CouplingTransferCountEnum,VertexSIdEnum);
     1151                        for (int i=0;i<nv;i++) forcing[i]/=transfercount[i]; //Divide forcing at this vertex by the number of icecaps that share it. This way we average the forcing when adding it into the earth model.
     1152                }
    10911153        }
    10921154
    … …  
    11271189                GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
    11281190
     1191                forcingglobal->Set(0.0);
     1192
    11291193                /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
    11301194                femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelchangeTransitionsEnum);
    … …  
    11461210                                        /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
    11471211                                         * ice cap: */
     1212
    11481213                                        forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
    11491214                                        xDelete<int>(index);

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

    r27299 r27308  
    312312        LoveMinIntegrationStepsEnum,
    313313        LoveMaxIntegrationdrEnum,
     314        LoveIntegrationSchemeEnum,
    314315        LoveKernelsEnum,
    315316        LoveMu0Enum,
    … …  
    935936        SealevelUNorthEsaEnum,
    936937        SealevelchangeIndicesEnum,
     938        SealevelchangeConvolutionVerticesEnum,
    937939        SealevelchangeAlphaIndexEnum,
    938940        SealevelchangeAzimuthIndexEnum,
    … …  
    953955        SealevelchangeViscousNEnum,
    954956        SealevelchangeViscousEEnum,
     957        CouplingTransferCountEnum,
    955958        SedimentHeadEnum,
    956959        SedimentHeadOldEnum,

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

    r27298 r27308  
    320320                case LoveMinIntegrationStepsEnum : return "LoveMinIntegrationSteps";
    321321                case LoveMaxIntegrationdrEnum : return "LoveMaxIntegrationdr";
     322                case LoveIntegrationSchemeEnum : return "LoveIntegrationScheme";
    322323                case LoveKernelsEnum : return "LoveKernels";
    323324                case LoveMu0Enum : return "LoveMu0";
    … …  
    941942                case SealevelUNorthEsaEnum : return "SealevelUNorthEsa";
    942943                case SealevelchangeIndicesEnum : return "SealevelchangeIndices";
     944                case SealevelchangeConvolutionVerticesEnum : return "SealevelchangeConvolutionVertices";
    943945                case SealevelchangeAlphaIndexEnum : return "SealevelchangeAlphaIndex";
    944946                case SealevelchangeAzimuthIndexEnum : return "SealevelchangeAzimuthIndex";
    … …  
    959961                case SealevelchangeViscousNEnum : return "SealevelchangeViscousN";
    960962                case SealevelchangeViscousEEnum : return "SealevelchangeViscousE";
     963                case CouplingTransferCountEnum : return "CouplingTransferCount";
    961964                case SedimentHeadEnum : return "SedimentHead";
    962965                case SedimentHeadOldEnum : return "SedimentHeadOld";

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

    r27298 r27308  
    326326              else if (strcmp(name,"LoveMinIntegrationSteps")==0) return LoveMinIntegrationStepsEnum;
    327327              else if (strcmp(name,"LoveMaxIntegrationdr")==0) return LoveMaxIntegrationdrEnum;
     328              else if (strcmp(name,"LoveIntegrationScheme")==0) return LoveIntegrationSchemeEnum;
    328329              else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum;
    329330              else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
    … …  
    382383              else if (strcmp(name,"OutputBufferSizePointer")==0) return OutputBufferSizePointerEnum;
    383384              else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
    384               else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
     388              if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
     389              else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
    389390              else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
    390391              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
    … …  
    505506              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
    506507              else if (strcmp(name,"SmbARMAInitialTime")==0) return SmbARMAInitialTimeEnum;
    507               else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum;
     511              if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum;
     512              else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum;
    512513              else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum;
    513514              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    … …  
    628629              else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    629630              else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
    630               else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
     634              if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
     635              else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
    635636              else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
    636637              else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum;
    … …  
    751752              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
    752753              else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
    753               else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
     757              if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
     758              else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
    758759              else if (strcmp(name,"DeviatoricStressyz")==0) return DeviatoricStressyzEnum;
    759760              else if (strcmp(name,"DeviatoricStresszz")==0) return DeviatoricStresszzEnum;
    … …  
    874875              else if (strcmp(name,"MasstransportSpcthickness")==0) return MasstransportSpcthicknessEnum;
    875876              else if (strcmp(name,"MaterialsRheologyB")==0) return MaterialsRheologyBEnum;
    876               else if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum;
     880              if (strcmp(name,"MaterialsRheologyBbar")==0) return MaterialsRheologyBbarEnum;
     881              else if (strcmp(name,"MaterialsRheologyE")==0) return MaterialsRheologyEEnum;
    881882              else if (strcmp(name,"MaterialsRheologyEbar")==0) return MaterialsRheologyEbarEnum;
    882883              else if (strcmp(name,"MaterialsRheologyEc")==0) return MaterialsRheologyEcEnum;
    … …  
    962963              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    963964              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
     965              else if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
    964966              else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
    965967              else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
    … …  
    980982              else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
    981983              else if (strcmp(name,"SealevelchangeViscousE")==0) return SealevelchangeViscousEEnum;
     984              else if (strcmp(name,"CouplingTransferCount")==0) return CouplingTransferCountEnum;
    982985              else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    983986              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    … …  
    995998              else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
    996999              else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum;
    997               else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    998               else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum;
    999               else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
     1003              if (strcmp(name,"SmbA")==0) return SmbAEnum;
     1004              else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum;
     1005              else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
     1006              else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum;
    10041007              else if (strcmp(name,"SmbAdiffini")==0) return SmbAdiffiniEnum;
    10051008              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    … …  
    11181121              else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum;
    11191122              else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum;
    1120               else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
    1121               else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
    1122               else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
     1126              if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum;
     1127              else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum;
     1128              else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum;
     1129              else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
    11271130              else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
    11281131              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
    … …  
    12411244              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    12421245              else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
    1243               else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
    1244               else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
    1245               else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
     1249              if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
     1250              else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
     1251              else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
     1252              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
    12501253              else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
    12511254              else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
    … …  
    13641367              else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
    13651368              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    1366               else if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum;
    1367               else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum;
    1368               else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
     1372              if (strcmp(name,"DebrisAnalysis")==0) return DebrisAnalysisEnum;
     1373              else if (strcmp(name,"DebrisSolution")==0) return DebrisSolutionEnum;
     1374              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
     1375              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    13731376              else if (strcmp(name,"Dense")==0) return DenseEnum;
    13741377              else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
    … …  
    14871490              else if (strcmp(name,"LoveHf")==0) return LoveHfEnum;
    14881491              else if (strcmp(name,"LoveHt")==0) return LoveHtEnum;
    1489               else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    1490               else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    1491               else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
     1495              if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
     1496              else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
     1497              else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
     1498              else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
    14961499              else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
    14971500              else if (strcmp(name,"LoveLt")==0) return LoveLtEnum;
    … …  
    16101613              else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
    16111614              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    1612               else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    1613               else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    1614               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
     1618              if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
     1619              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
     1620              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
     1621              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    16191622              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    16201623              else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum;

  • issm/trunk-jpl/src/m/classes/model.m

    r27030 r27308  
    137137                        if isa(md.amr,'double'); md.amr=amr(); end
    138138                        %2017 Aug 29th
    139                         if isa(md.love,'double'); md.love=fourierlove(); end
     139                        if isa(md.love,'double'); md.love=love(); end
    140140                        %2017 Oct 26th
    141141                        if isa(md.calving,'calvingdev')
    … …  
    289289                        md.calving          = calving();
    290290                        md.frontalforcings  = frontalforcings();
    291                         md.love             = fourierlove();
     291                        md.love             = love();
    292292                        md.esa              = esa();
    293293                        md.sampling         = sampling();

  • issm/trunk-jpl/src/m/classes/sealevelmodel.m

    r26358 r27308  
    4848
    4949                        %is the coupler turned on?
    50                         for i=1:length(slm.icecaps),
    51                                 if slm.icecaps{i}.transient.iscoupler==0,
    52                                         warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
    53                                 end
    54                         end
     50                        %for i=1:length(slm.icecaps),
     51                        %       if slm.icecaps{i}.transient.iscoupler==0,
     52                        %               warning(sprintf('sealevelmodel.m::checkconsistency: icecap model %s should have the transient coupler option turned on!',slm.icecaps{i}.miscellaneous.name));
     53                        %       end
     54                        %end
    5555                               
    56                         if slm.earth.transient.iscoupler==0,
    57                                 warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!');
     56                        %if slm.earth.transient.iscoupler==0,
     57                        %       warning('sealevelmodel.m::checkconsistency: earth model should have the transient coupler option turned on!');
     58                        %end
     59
     60                        %check that the transition vectors have the right size:
     61
     62                        if slm.earth.mesh.numberofvertices ~= length(slm.earth.solidearth.transfercount)
     63                                error('sealevelmodel.m::checkconsistency: earth.solidearth.transfercount should be of size earth.mesh.numberofvertices')
    5864                        end
    5965
    … …  
    7581                        for i=1:length(slm.icecaps),
    7682                                md= slm.icecaps{i};
    77                                 if ~isempty(find(md.dsl.steric_rate - slm.earth.dsl.steric_rate(slm.earth.dsl.transitions{i}))),
     83                                if ~isempty(find(md.dsl.sea_surface_height_above_geoid - slm.earth.dsl.sea_surface_height_above_geoid(slm.transitions{i}))),
    7884                                        error(sprintf('sealevelmodel.m::checkconsistency: steric rate on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
    7985                                end
    … …  
    8389                        for i=1:length(slm.icecaps),
    8490                                md= slm.icecaps{i};
    85                                 if md.solidearthsettings.isgrd~=slm.earth.solidearthsettings.isgrd
     91                                if md.solidearth.settings.isgrd~=slm.earth.solidearth.settings.isgrd
    8692                                        error(sprintf('sealevelmodel.m::checkconsistency: isgrd on ice cap %s is not the same as for the earth\n',md.miscellaneous.name));
    8793                                end
    … …  
    228234                        self.transitions={};
    229235                        self.eltransitions={};
     236                        self.earth.solidearth.transfercount=zeros(self.earth.mesh.numberofvertices,1);
    230237
    231238                        %for elements:
    … …  
    248255
    249256                                self.eltransitions{end+1}=meshintersect3d(xe,ye,ze,xei,yei,zei,'force',force);
     257
     258                                self.earth.solidearth.transfercount(self.transitions{i})=self.earth.solidearth.transfercount(self.transitions{i})+1;
     259                        end
     260
     261                        for i=1:length(self.icecaps),
     262                                self.icecaps{i}.solidearth.transfercount=self.earth.solidearth.transfercount(self.transitions{i});
    250263                        end
    251264                end % }}}
    … …  
    255268                                flags(self.transitions{i})=i;
    256269                        end
    257                         plotmodel(self.earth,'data',flags,'coastline','on');
     270                        plotmodel(self.earth,'data',flags,'coastlines','on');
    258271
    259272                end % }}}

  • issm/trunk-jpl/src/m/classes/solidearth.m

    r26358 r27308  
    1414                requested_outputs      = {};
    1515                transitions            = {};
     16                transfercount           = [];
    1617                partitionice           = [];
    1718                partitionhydro         = [];
    … …  
    5152                        fielddisplay(self,'planetradius','planet radius [m]');
    5253                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
     54                        fielddisplay(self,'transfercount','number of icecaps vertices are part of');
    5355                        fielddisplay(self,'requested_outputs','additional outputs requested');
    5456                        fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
    … …  
    7173                        %transitions should be a cell array of vectors:
    7274                        self.transitions={};
     75                        self.transfercount=[0];
    7376
    7477                        %no partitions requested for barystatic contribution:
    … …  
    110113                        WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');
    111114                        WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray');
     115                        %WriteData(fid,prefix,'object',self,'fieldname','transfercount','data',self.transfercount,'format','DoubleMat','name', 'md.solidearth.transfercount','mattype',1);
     116                        %WriteData(fid,prefix,'data',self.transfercount,'format','DoubleMat','name', 'md.solidearth.transfercount');
     117                        WriteData(fid,prefix,'object',self,'fieldname','transfercount','format','DoubleMat','mattype',1);
    112118
    113119                        if ~isempty(self.partitionice),

  • issm/trunk-jpl/src/m/solve/solveslm.m

    r26358 r27308  
    6868nps{end+1}=slm.earth.cluster.np;
    6969
    70 BuildQueueScriptMultipleModels(cluster,slm.private.runtimename,slm.miscellaneous.name,slm.private.solution,valgrind,privateruntimenames,miscellaneousnames,nps);
     70BuildQueueScriptMultipleModels(cluster,slm.private.runtimename,slm.miscellaneous.name,slm.private.solution,privateruntimenames,miscellaneousnames,nps);
    7171
    7272%Upload all required files, given that each individual solution for icecaps and earth model already did:

  • issm/trunk-jpl/test/NightlyRun/test2013.m

    r27133 r27308  
    11%Test Name: EarthSlc_Geometry
    22
    3 step=[1 2];
     3step=[1];
    44if any(step==1)
    55%mesh earth:

  • issm/trunk-jpl/test/NightlyRun/test2070.m

    r26805 r27308  
    4040md.love.pw_threshold=1e-3;
    4141md.love.Gravitational_Constant=6.6732e-11;
    42 md.love.integration_steps_per_layer=100;
     42md.love.min_integration_steps=100;
    4343md.love.allow_layer_deletion=1;
    4444md.love.forcing_type=11;

  • issm/trunk-jpl/test/NightlyRun/test2071.m

    r26884 r27308  
    4040md.love.pw_threshold=1e-3;
    4141md.love.Gravitational_Constant=6.6732e-11;
    42 md.love.integration_steps_per_layer=100;
     42md.love.min_integration_steps=100;
    4343md.love.allow_layer_deletion=1;
    4444md.love.forcing_type=11;

  • issm/trunk-jpl/test/NightlyRun/test2072.m

    r26807 r27308  
    4444md.love.pw_threshold=1e-3;
    4545md.love.Gravitational_Constant=6.6732e-11;
    46 md.love.integration_steps_per_layer=100;
     46md.love.min_integration_steps=100;
    4747md.love.allow_layer_deletion=1;
    4848md.love.forcing_type=11;

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

    r26881 r27308  
    4444md.love.pw_threshold=1e-3;
    4545md.love.Gravitational_Constant=6.6732e-11;
    46 md.love.integration_steps_per_layer=100;
    4746md.love.allow_layer_deletion=1;
    4847md.love.forcing_type=11;
Note: See TracChangeset for help on using the changeset viewer.