Changeset 26800


Ignore:
Timestamp:
01/20/22 17:02:07 (3 years ago)
Author:
caronlam
Message:

BUG: Sealevelchange core: load weights did not add up correctly; BUG: Sealevel change core viscous stacks were offset in time; partial load area phi was double counted, appearing both in Green functions and loads; NEW: love core optimization, paralelization supported, support for viscous rotational feedback; NEW:Sea level change core: support for viscous rotational feedback; CHG: grd loads now contain the average kg.m-2 over the (sub-)elemental area they apply in, not the average*phi, phi is simply accounted for in the loadarea or Green's function; NEW: nightly runs 2070, 2071, 2072, 2091, 2092 test Spada et al (2011) benchmark cases

Location:
issm/trunk-jpl
Files:
162 added
39 edited

Legend:

Unmodified
Added
Removed

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

    r26468 r26800  
    4444}/*}}}*/
    4545void HydrologyTwsAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
     46
     47        /*retrieve some parameters: */
     48        int    hydrology_model;
     49        int    numoutputs;
     50        char** requestedoutputs = NULL;
     51        iomodel->FindConstant(&hydrology_model,"md.hydrology.model");
     52
     53        /*Now, do we really want Tws?*/
     54        if(hydrology_model!=HydrologyTwsEnum) return;
     55
     56        parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model));
     57
     58        /*Requested outputs*/
     59        iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs");
     60        parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs));
     61        if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs));
     62        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs");
    4663
    4764}/*}}}*/

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

    r26242 r26800  
    3333        parameters->AddObject(iomodel->CopyConstantObject("md.love.mu0",LoveMu0Enum));
    3434        parameters->AddObject(iomodel->CopyConstantObject("md.love.Gravitational_Constant",LoveGravitationalConstantEnum));
     35        parameters->AddObject(iomodel->CopyConstantObject("md.love.chandler_wobble",LoveChandlerWobbleEnum));
    3536        parameters->AddObject(iomodel->CopyConstantObject("md.love.allow_layer_deletion",LoveAllowLayerDeletionEnum));
    3637        parameters->AddObject(iomodel->CopyConstantObject("md.love.underflow_tol",LoveUnderflowTolEnum));
     38        parameters->AddObject(iomodel->CopyConstantObject("md.love.pw_threshold",LovePostWidderThresholdEnum));
    3739        parameters->AddObject(iomodel->CopyConstantObject("md.love.integration_steps_per_layer",LoveIntStepsPerLayerEnum));
    3840        parameters->AddObject(iomodel->CopyConstantObject("md.love.istemporal",LoveIsTemporalEnum));
    … …  
    4345        parameters->AddObject(iomodel->CopyConstantObject("md.love.core_mantle_boundary",LoveCoreMantleBoundaryEnum));
    4446        parameters->AddObject(iomodel->CopyConstantObject("md.love.complex_computation",LoveComplexComputationEnum));
     47        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.equatorialmoi",RotationalEquatorialMoiEnum));
     48        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.polarmoi",RotationalPolarMoiEnum));
     49        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.rotational.angularvelocity",RotationalAngularVelocityEnum));
     50        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    4551}/*}}}*/
    4652void LoveAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/

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

    r26468 r26800  
    8585        IssmDouble* love_tk=NULL;
    8686        IssmDouble* love_tl=NULL;
     87        IssmDouble* love_pmtf_colinear=NULL;
     88        IssmDouble* love_pmtf_ortho=NULL;
    8789        IssmDouble* love_timefreq=NULL;
    8890        bool        love_istime=true;
    … …  
    101103        IssmDouble* H_viscoelastic_interpolated= NULL;
    102104        IssmDouble* H_viscoelastic_local = NULL;
     105        IssmDouble* Pmtf_col_interpolated = NULL;
     106        IssmDouble* Pmtf_ortho_interpolated = NULL;
     107        IssmDouble* Pmtf_z_interpolated = NULL;
     108        IssmDouble* Love_th2_interpolated = NULL;
     109        IssmDouble* Love_tk2_interpolated = NULL;
     110        IssmDouble* Love_tl2_interpolated = NULL;
     111
    103112        int         M,m,lower_row,upper_row;
    104113        IssmDouble  degacc=.01;
    … …  
    166175        /*compute planet area and plug into parameters:*/
    167176        iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
    168         planetarea=4*PI*planetradius*planetradius;
     177        planetarea=4*M_PI*planetradius*planetradius;
    169178        parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
    170179
    … …  
    253262                        iomodel->FetchData(&love_k,&ndeg,&precomputednt,"md.solidearth.lovenumbers.k");
    254263                        iomodel->FetchData(&love_l,&ndeg,&precomputednt,"md.solidearth.lovenumbers.l");
    255                         iomodel->FetchData(&love_th,&ndeg,&precomputednt,"md.solidearth.lovenumbers.th");
    256                         iomodel->FetchData(&love_tk,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tk");
    257                         iomodel->FetchData(&love_tl,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tl");
    258264
    259265                        parameters->AddObject(new DoubleParam(SolidearthSettingsTimeAccEnum,timeacc));
    … …  
    261267                        parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
    262268                        parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
    263                         parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
    264                         parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
    265                         parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
     269
     270                        if (rotation){
     271                                iomodel->FetchData(&love_th,&ndeg,&precomputednt,"md.solidearth.lovenumbers.th");
     272                                iomodel->FetchData(&love_tk,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tk");
     273                                iomodel->FetchData(&love_tl,&ndeg,&precomputednt,"md.solidearth.lovenumbers.tl");
     274                                iomodel->FetchData(&love_pmtf_colinear,&dummy,&precomputednt,"md.solidearth.lovenumbers.pmtf_colinear");
     275                                iomodel->FetchData(&love_pmtf_ortho,&dummy,&precomputednt,"md.solidearth.lovenumbers.pmtf_ortho");
     276
     277                                parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,love_pmtf_colinear,1,precomputednt));
     278                                parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,love_pmtf_ortho,1,precomputednt));
     279                                parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
     280                                parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
     281                                parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
     282                        }
     283
    266284                        parameters->AddObject(new DoubleMatParam(LoveTimeFreqEnum,love_timefreq,precomputednt,1));
    267285                        parameters->AddObject(new BoolParam(LoveIsTimeEnum,love_istime));
    … …  
    283301                                parameters->AddObject(new IntParam(SealevelchangeViscousIndexEnum,0));
    284302                                xDelete<IssmDouble>(viscoustimes);
     303                                if (rotation){
     304                                        IssmDouble* viscouspolarmotion=NULL;
     305                                        viscouspolarmotion=xNewZeroInit<IssmDouble>(3*nt);
     306                                        parameters->AddObject(new DoubleMatParam(SealevelchangeViscousPolarMotionEnum,viscouspolarmotion,3,nt));
     307                                        xDelete<IssmDouble>(viscouspolarmotion);
     308                                }
    285309                        }
    286310                        else {
    … …  
    307331                }
    308332
    309                 if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    310 
     333                if(rotation) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    311334                if(selfattraction){
    312335
    313                         /*compute combined legendre + love number (elastic green function:*/
     336                        /*compute combined legendre + love number (elastic green function):*/
    314337                        m=DetermineLocalSize(M,IssmComm::GetComm());
    315338                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    … …  
    441464                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
    442465                                if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
     466                                if(rotation){
     467                                        Pmtf_col_interpolated=xNew<IssmDouble>(nt,"t");
     468                                        Pmtf_ortho_interpolated=xNew<IssmDouble>(nt,"t");
     469                                        Pmtf_z_interpolated=xNew<IssmDouble>(nt,"t");
     470                                        Love_tk2_interpolated=xNew<IssmDouble>(nt,"t");
     471                                        Love_th2_interpolated=xNew<IssmDouble>(nt,"t");
     472                                        if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(nt,"t");
     473                                }
    443474#else
    444475                                G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
    445476                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
    446477                                if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
     478                                if(rotation){
     479                                        Pmtf_col_interpolated=xNew<IssmDouble>(nt);
     480                                        Pmtf_ortho_interpolated=xNew<IssmDouble>(nt);
     481                                        Pmtf_z_interpolated=xNew<IssmDouble>(nt);
     482                                        Love_tk2_interpolated=xNew<IssmDouble>(nt);
     483                                        Love_th2_interpolated=xNew<IssmDouble>(nt);
     484                                        if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(nt);
     485                                }
    447486#endif
    448487
    … …  
    477516                                                if(horiz)H_viscoelastic_interpolated[timeindex]=(1-lincoeff)*H_viscoelastic[timepreindex]+lincoeff*H_viscoelastic[timepreindex+1];
    478517                                        }
     518
     519                                        if(rotation){
     520                                                int timepreindex= 2*precomputednt+timeindex2;
     521                                                Pmtf_col_interpolated[t]=(1.0-lincoeff)*love_pmtf_colinear[timeindex2]+lincoeff*love_pmtf_colinear[timeindex2+1];
     522                                                Pmtf_ortho_interpolated[t]=(1.0-lincoeff)*love_pmtf_ortho[timeindex2]+lincoeff*love_pmtf_ortho[timeindex2+1];
     523                                                Pmtf_z_interpolated[t]=1.0+(1.0-lincoeff)*love_k[timepreindex]+lincoeff*love_k[timepreindex+1];
     524                                                Love_tk2_interpolated[t]=(1.0-lincoeff)*love_tk[timepreindex]+lincoeff*love_tk[timepreindex+1];
     525                                                Love_th2_interpolated[t]=(1.0-lincoeff)*love_th[timepreindex]+lincoeff*love_th[timepreindex+1];
     526                                                if (horiz) Love_tl2_interpolated[t]=(1.0-lincoeff)*love_tl[timepreindex]+lincoeff*love_tl[timepreindex+1];
     527                                        }
    479528                                }
    480529
    … …  
    491540                                if(horiz)H_viscoelastic_interpolated=H_viscoelastic;
    492541#endif
     542
     543                                if(rotation){ //if this cpu handles degree 2
     544#ifdef _HAVE_AD_
     545                                        Pmtf_col_interpolated=xNew<IssmDouble>(1,"t");
     546                                        Pmtf_ortho_interpolated=xNew<IssmDouble>(1,"t");
     547                                        Pmtf_z_interpolated=xNew<IssmDouble>(1,"t");
     548                                        Love_tk2_interpolated=xNew<IssmDouble>(1,"t");
     549                                        Love_th2_interpolated=xNew<IssmDouble>(1,"t");
     550                                        if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1,"t");
     551#else
     552                                        Pmtf_col_interpolated=xNew<IssmDouble>(1);
     553                                        Pmtf_ortho_interpolated=xNew<IssmDouble>(1);
     554                                        Pmtf_z_interpolated=xNew<IssmDouble>(1);
     555                                        Love_tk2_interpolated=xNew<IssmDouble>(1);
     556                                        Love_th2_interpolated=xNew<IssmDouble>(1);
     557                                        if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1);
     558#endif
     559
     560                                        Pmtf_col_interpolated=love_pmtf_colinear;
     561                                        Pmtf_ortho_interpolated=love_pmtf_ortho;
     562                                        Pmtf_z_interpolated[0]=1.0+love_k[2];
     563                                        Love_tk2_interpolated[0]=love_tk[2];
     564                                        Love_th2_interpolated[0]=love_th[2];
     565                                        if (horiz) Love_tl2_interpolated[0]=love_tl[2];
     566                                }
    493567                        }       
    494568
    … …  
    499573                                parameters->AddObject(new DoubleVecParam(SealevelchangeUViscoElasticEnum,U_viscoelastic_interpolated,M*nt));
    500574                                if(horiz)parameters->AddObject(new DoubleVecParam(SealevelchangeHViscoElasticEnum,H_viscoelastic_interpolated,M*nt));
     575                                if(rotation){
     576                                        parameters->AddObject(new DoubleVecParam(SealevelchangePolarMotionTransferFunctionColinearEnum,Pmtf_col_interpolated,nt));
     577                                        parameters->AddObject(new DoubleVecParam(SealevelchangePolarMotionTransferFunctionOrthogonalEnum,Pmtf_ortho_interpolated,nt));
     578                                        parameters->AddObject(new DoubleVecParam(SealevelchangePolarMotionTransferFunctionZEnum,Pmtf_z_interpolated,nt));
     579                                        parameters->AddObject(new DoubleVecParam(SealevelchangeTidalH2Enum,Love_th2_interpolated,nt));
     580                                        parameters->AddObject(new DoubleVecParam(SealevelchangeTidalK2Enum,Love_tk2_interpolated,nt));
     581                                        if (horiz) parameters->AddObject(new DoubleVecParam(SealevelchangeTidalL2Enum,Love_tl2_interpolated,nt));
     582                                }
    501583                        }
    502584
    … …  
    519601                                        xDelete<IssmDouble>(H_viscoelastic_local);
    520602                                }
     603                                if(rotation){
     604                                        xDelete<IssmDouble>(love_pmtf_colinear);
     605                                        xDelete<IssmDouble>(love_pmtf_ortho);
     606
     607                                }
    521608                        }
    522609                } /*}}}*/
    … …  
    545632        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.solidearth.requested_outputs");
    546633        /*}}}*/
    547 
    548634}/*}}}*/
    549635

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

    r26784 r26800  
    396396
    397397                virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
    398                 virtual void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
    399398                virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
    400399                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    … …  
    405404                virtual void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom)=0;
    406405                virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom)=0;
    407                 virtual void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom)=0;
    408                 virtual void       SealevelchangeUpdateViscousFields()=0;
     406                virtual void       SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset)=0;
    409407                #endif
    410408

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

    r26626 r26800  
    225225                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    226226                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    227                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
    228227                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    229228                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    … …  
    234233                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    235234                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    236                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    237                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
     235                void       SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){_error_("not implemented yet");};
    238236                #endif
    239237

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

    r26296 r26800  
    179179                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    180180                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    181                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
    182181                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    183182                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    … …  
    188187                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    189188                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    190                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    191                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
     189                void       SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){_error_("not implemented yet");};
    192190#endif
    193191

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

    r26296 r26800  
    186186                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");};
    187187                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    188                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom){_error_("not implemented yet");};
    189188                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    190189                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    … …  
    195194                void       SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){_error_("not implemented yet");};
    196195                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads,  IssmDouble* rotationvector,SealevelGeometry* slgeom){_error_("not implemented yet");};
    197                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){_error_("not implemented yet");};
    198                 void       SealevelchangeUpdateViscousFields(){_error_("not implemented yet");};
     196                void       SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){_error_("not implemented yet");};
    199197#endif
    200198

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

    r26794 r26800  
    19701970        //compute sea level load weights
    19711971        this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset);
    1972 
     1972        for (int i=0;i<NUMVERTICES;i++) loadweights[i]/=phi;
    19731973        this->GetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius);
    19741974
    … …  
    20612061                this->GetFractionGeometry(loadweights,&phi2,&point2,&f,&g,&istrapeze2,levelset2);
    20622062                this->GetBarycenterFromLevelset(platbar,plongbar, phi2, f, g, late, longe, point2,istrapeze2,planetradius);
     2063                for (int i=0;i<NUMVERTICES;i++) loadweights[i]/=phi2;
    20632064                *ploadarea=area*phi2;
    20642065                return;
    … …  
    20672068                this->GetFractionGeometry(loadweights,&phi1,&point1,&d,&e,&istrapeze1,levelset1);
    20682069                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius);
     2070                for (int i=0;i<NUMVERTICES;i++) loadweights[i]/=phi1;
    20692071                *ploadarea=area*phi1;
    20702072                return;
    … …  
    20792081                this->GetBarycenterFromLevelset(platbar,plongbar, phi1, d, e, late, longe, point1,istrapeze1,planetradius);
    20802082                *ploadarea=area*phi1;
    2081                 for (int i=0;i<NUMVERTICES;i++) loadweights[i]=weights1[i];
     2083                for (int i=0;i<NUMVERTICES;i++) loadweights[i]=weights1[i]/phi1;
    20822084                return;
    20832085        }
    … …  
    21402142
    21412143        //interpolant weights of each point. Any field F[0,1,2] provided at the original vertices [0,1,2] will be equal on point k to sum_i (F[i] * w[k][i])
    2142         w[0][0]=1; //A
    2143         w[1][1]=1; //B
    2144         w[2][2]=1; //C
     2144        w[i0][i0]=1; //A
     2145        w[i1][i1]=1; //B
     2146        w[i2][i2]=1; //C
    21452147        w[3][i0]=1.0-d; w[3][i1]=d; //D
    21462148        w[4][i0]=1.0-e; w[4][i2]=e; //E
    … …  
    23642366                for (int j=0;j<3;j++){
    23652367                        for (int i=0;i<NUMVERTICES;i++) {
    2366                                 loadweights[j]=w1[i][j]*area1 + w2[i][j]*area2 + w3[i][j]*area3;
     2368                                loadweights[j]+=w1[i][j]*area1 + w2[i][j]*area2 + w3[i][j]*area3;
    23672369                                barycenter[j]+=xyz1[i][j]*area1+xyz2[i][j]*area2+xyz3[i][j]*area3;
    23682370                        }
    2369                         loadweights[j]/=area;
     2371                        loadweights[j]/=areasub*3.0;
    23702372                        barycenter[j]/=areasub *3.0;
    23712373                }
    … …  
    36313633        }
    36323634
    3633        
    36343635        /*Cleanup & return: */
    36353636        xDelete<int>(indices);
    … …  
    62856286        IssmDouble area,planetarea,planetradius;
    62866287        IssmDouble constant,ratioe;
    6287         IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    62886288        IssmDouble rho_earth;
    62896289        IssmDouble NewtonG;
    6290         IssmDouble g;
     6290        IssmDouble g, cent_scaling;
    62916291        IssmDouble lati,longi;
    62926292        IssmDouble latitude[NUMVERTICES];
    … …  
    63396339
    63406340        /*Rotational:*/
     6341        #ifdef _HAVE_RESTRICT_
     6342        IssmDouble* __restrict__ tide_love_h  = NULL;
     6343        IssmDouble* __restrict__ tide_love_k  = NULL;
     6344        IssmDouble* __restrict__ tide_love_l  = NULL;
     6345        IssmDouble* __restrict__ LoveRotRSL   = NULL;
     6346        IssmDouble* __restrict__ LoveRotU     = NULL;
     6347        IssmDouble* __restrict__ LoveRothoriz = NULL;
     6348        IssmDouble* __restrict__  Grot        = NULL;
     6349        IssmDouble* __restrict__ GUrot        = NULL;
     6350        IssmDouble* __restrict__ GNrot        = NULL;
     6351        IssmDouble* __restrict__ GErot        = NULL;
     6352        #else
    63416353        IssmDouble* tide_love_h  = NULL;
    63426354        IssmDouble* tide_love_k  = NULL;
    63436355        IssmDouble* tide_love_l  = NULL;
     6356        IssmDouble* LoveRotRSL   = NULL;
     6357        IssmDouble* LoveRotU     = NULL;
     6358        IssmDouble* LoveRothoriz = NULL;
     6359        IssmDouble*  Grot        = NULL;
     6360        IssmDouble* GUrot        = NULL;
     6361        IssmDouble* GNrot        = NULL;
     6362        IssmDouble* GErot        = NULL;
     6363        #endif
     6364
    63446365        IssmDouble  moi_e, moi_p, omega;
    6345         IssmDouble  Grotm1[3],GUrotm1[3],GNrotm1[3],GErotm1[3];
    6346         IssmDouble  Grotm2[3],GUrotm2[3],GNrotm2[3],GErotm2[3];
    6347         IssmDouble  Grotm3[3],GUrotm3[3],GNrotm3[3],GErotm3[3];
    63486366        IssmDouble  Y21cos     , Y21sin     , Y20;
    63496367        IssmDouble dY21cos_dlat,dY21sin_dlat,dY20_dlat;
    63506368        IssmDouble dY21cos_dlon,dY21sin_dlon;
    6351         IssmDouble LoveRotRSL,LoveRotU,LoveRothoriz;
    63526369        IssmDouble polenudge;
    63536370        /*}}}*/
    … …  
    63726389
    63736390        if(computerotation){
    6374                 parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
    6375                 parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
    6376                 parameters->FindParam(&tide_love_l,NULL,NULL,TidalLoveLEnum);
    63776391                parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
    63786392                parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
    63796393                parameters->FindParam(&omega,RotationalAngularVelocityEnum);
     6394                //parameters->FindParam(&tide_love_h,NULL,NULL,SealevelchangeTidalH2Enum);
     6395                //parameters->FindParam(&tide_love_k,NULL,NULL,SealevelchangeTidalK2Enum);
     6396                //if(horiz) parameters->FindParam(&tide_love_l,NULL,NULL,SealevelchangeTidalL2Enum);
    63806397        }
    63816398        /*}}}*/
    … …  
    63956412                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    63966413                        parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
     6414                }
     6415
     6416                if(computerotation){
     6417                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalH2Enum)); _assert_(parameter);
     6418                        parameter->GetParameterValueByPointer((IssmDouble**)&tide_love_h,NULL);
     6419
     6420                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalK2Enum)); _assert_(parameter);
     6421                        parameter->GetParameterValueByPointer((IssmDouble**)&tide_love_k,NULL);
     6422
     6423                        if (horiz) {
     6424                                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalL2Enum)); _assert_(parameter);
     6425                                parameter->GetParameterValueByPointer((IssmDouble**)&tide_love_l,NULL);
     6426                        }
    63976427                }
    63986428        }
    … …  
    64356465                        IssmDouble alpha;
    64366466                        IssmDouble delPhi,delLambda;
    6437                         /*recover info for this element and vertex:*/
     6467                        /*recovers info for this element and vertex:*/
    64386468                        IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
    64396469                        IssmDouble longe= atan2(yye[e],xxe[e]);
    … …  
    64426472                        longi=longitude[i];
    64436473
    6444                         /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6474                        /*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
    64456475                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    64466476                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    … …  
    65106540        /*Compute rotation kernel:{{{*/
    65116541        if(computerotation){
     6542                //initialization
     6543                LoveRotRSL  = xNewZeroInit<IssmDouble>(nt);
     6544                LoveRotU    = xNewZeroInit<IssmDouble>(nt);
     6545                LoveRothoriz= xNewZeroInit<IssmDouble>(nt);
     6546                Grot        = xNewZeroInit<IssmDouble>(3*3*nt); //3 polar motion components * 3 vertices * number of time steps
     6547                GUrot       = xNewZeroInit<IssmDouble>(3*3*nt);
     6548
     6549                if (horiz){
     6550                        GErot=xNewZeroInit<IssmDouble>(3*3*nt);
     6551                        GNrot=xNewZeroInit<IssmDouble>(3*3*nt);
     6552                }
    65126553
    65136554                /*What is the gravity at planet's surface: */
    65146555                g=4.0/3.0*M_PI*rho_earth*NewtonG*planetradius;
    6515 
    6516                 //Amplitude of the rotational feedback
    6517                 LoveRotRSL=((1.0+tide_love_k[2]-tide_love_h[2])/g)*pow(omega*planetradius,2.0);
    6518                 LoveRotU=(tide_love_h[2]/g)*pow(omega*planetradius,2.0);
    6519                 LoveRothoriz=(tide_love_l[2]/g)*pow(omega*planetradius,2.0);
    6520 
     6556                cent_scaling=pow(omega*planetradius,2.0); //centrifugal potential dimensioning constant
     6557                for(int t=0;t<nt;t++){
     6558                        //Amplitude of the rotational feedback
     6559                        //to speed up calculation we include the dimension constant r^2*Omega^2/g, so all 3 of these are in meters
     6560                        LoveRotRSL[t]=((1.0+tide_love_k[t]-tide_love_h[t])/g)*cent_scaling;
     6561                        LoveRotU[t]=(tide_love_h[t]/g)*cent_scaling;
     6562                        if (horiz) LoveRothoriz[t]=(tide_love_l[t]/g)*cent_scaling;
     6563                }
    65216564                for(int i=0;i<3;i++){
    65226565
    … …  
    65376580                        longi=longitude[i];
    65386581
    6539                         //Spherical harmonic functions of degree 2
     6582                        //Spherical harmonic functions of degree 2 (spatial pattern of rotation)
    65406583                        Y21cos= -0.5*sin(2.*lati)*cos(longi);
    65416584                        Y21sin= -0.5*sin(2.*lati)*sin(longi);
    65426585                        Y20   = -(1.0/6.0 - 0.5*cos(2.0*lati));
    65436586
    6544                         Grotm1[i]= LoveRotRSL*Y21cos;
    6545                         Grotm2[i]= LoveRotRSL*Y21sin;
    6546                         Grotm3[i]= LoveRotRSL*Y20;
    6547 
    6548                         if (computeelastic){
    6549                                 GUrotm1[i]= LoveRotU*Y21cos;
    6550                                 GUrotm2[i]= LoveRotU*Y21sin;
    6551                                 GUrotm3[i]= LoveRotU*Y20;
    6552                                 if (horiz){
     6587                        if (computeelastic && horiz){
    65536588                                //bed_N = Love_l * d(Y)/dlat ;
    65546589                                dY21cos_dlat=-cos(2.*lati)*cos(longi);
    65556590                                dY21sin_dlat=-cos(2.*lati)*sin(longi);
    65566591                                dY20_dlat=   -sin(2.*lati);
    6557                                 GNrotm1[i]= LoveRothoriz*dY21cos_dlat;
    6558                                 GNrotm2[i]= LoveRothoriz*dY21sin_dlat;
    6559                                 GNrotm3[i]= LoveRothoriz*dY20_dlat;
    65606592
    65616593                                //bed_E = Love_l * 1/cos(lat) * d(Y)/dlon ;
    … …  
    65636595                                dY21sin_dlon=Y21cos/cos(lati);
    65646596                                //dY20_dlon=0.;
    6565 
    6566                                 GErotm1[i]= LoveRothoriz*dY21cos_dlon;
    6567                                 GErotm2[i]= LoveRothoriz*dY21sin_dlon;
    6568                                 GErotm3[i]= 0.0;
     6597                        }
     6598
     6599                        for(int t=0;t<nt;t++){
     6600
     6601                                Grot[0*3*nt+i*nt+t]= LoveRotRSL[t]*Y21cos; //x component of polar motion
     6602                                Grot[1*3*nt+i*nt+t]= LoveRotRSL[t]*Y21sin; //y
     6603                                Grot[2*3*nt+i*nt+t]= LoveRotRSL[t]*Y20;    //z
     6604
     6605                                if (computeelastic){
     6606                                        GUrot[0*3*nt+i*nt+t]= LoveRotU[t]*Y21cos;
     6607                                        GUrot[1*3*nt+i*nt+t]= LoveRotU[t]*Y21sin;
     6608                                        GUrot[2*3*nt+i*nt+t]= LoveRotU[t]*Y20;
     6609                                        if (horiz){
     6610                                                //bed_N = Love_l * d(Y)/dlat ;
     6611                                                GNrot[0*3*nt+i*nt+t]= LoveRothoriz[t]*dY21cos_dlat;
     6612                                                GNrot[1*3*nt+i*nt+t]= LoveRothoriz[t]*dY21sin_dlat;
     6613                                                GNrot[2*3*nt+i*nt+t]= LoveRothoriz[t]*dY20_dlat;
     6614
     6615                                                //bed_E = Love_l * 1/cos(lat) * d(Y)/dlon ;
     6616                                                GErot[0*3*nt+i*nt+t]= LoveRothoriz[t]*dY21cos_dlon;
     6617                                                GErot[1*3*nt+i*nt+t]= LoveRothoriz[t]*dY21sin_dlon;
     6618                                                GErot[2*3*nt+i*nt+t]= 0.0;
     6619                                        }
    65696620                                }
    65706621                        }
    65716622                }
    6572                 this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
    6573                 this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
    6574                 this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
     6623                this->inputs->SetArrayInput(SealevelchangeGrotEnum,this->lid,Grot,3*3*nt);
    65756624                if (computeelastic){
    6576                         this->AddInput(SealevelGUrotm1Enum,&GUrotm1[0],P1Enum);
    6577                         this->AddInput(SealevelGUrotm2Enum,&GUrotm2[0],P1Enum);
    6578                         this->AddInput(SealevelGUrotm3Enum,&GUrotm3[0],P1Enum);
     6625                        this->inputs->SetArrayInput(SealevelchangeGUrotEnum,this->lid,GUrot,3*3*nt);
    65796626                        if(horiz){
    6580                                 this->AddInput(SealevelGNrotm1Enum,&GNrotm1[0],P1Enum);
    6581                                 this->AddInput(SealevelGNrotm2Enum,&GNrotm2[0],P1Enum);
    6582                                 this->AddInput(SealevelGNrotm3Enum,&GNrotm3[0],P1Enum);
    6583                                 this->AddInput(SealevelGErotm1Enum,&GErotm1[0],P1Enum);
    6584                                 this->AddInput(SealevelGErotm2Enum,&GErotm2[0],P1Enum);
    6585                                 this->AddInput(SealevelGErotm3Enum,&GErotm3[0],P1Enum);
     6627                                this->inputs->SetArrayInput(SealevelchangeGNrotEnum,this->lid,GNrot,3*3*nt);
     6628                                this->inputs->SetArrayInput(SealevelchangeGErotEnum,this->lid,GErot,3*3*nt);
    65866629                        }
    65876630                }
    … …  
    65996642                        viscousN=xNewZeroInit<IssmDouble>(3*nt);
    66006643                        viscousE=xNewZeroInit<IssmDouble>(3*nt);
    6601                         this->inputs->SetArrayInput(SealevelchangeViscousNEnum,this->lid,viscousRSL,3*nt);
    6602                         this->inputs->SetArrayInput(SealevelchangeViscousEEnum,this->lid,viscousU,3*nt);
     6644                        this->inputs->SetArrayInput(SealevelchangeViscousNEnum,this->lid,viscousN,3*nt);
     6645                        this->inputs->SetArrayInput(SealevelchangeViscousEEnum,this->lid,viscousE,3*nt);
    66036646                }
    66046647        }
    … …  
    66146657                        delete GE;
    66156658                }
     6659                if(computerotation){
     6660                        delete Grot;
     6661                        delete GUrot;
     6662                        if (horiz){
     6663                                delete GNrot;
     6664                                delete GErot;
     6665                        }
     6666                }
    66166667        }
    66176668        #else
    … …  
    66226673                        xDelete(GN);
    66236674                        xDelete(GE);
     6675                }
     6676                if(computerotation){
     6677                        xDelete(Grot);
     6678                        xDelete(GUrot);
     6679                        if (horiz){
     6680                                xDelete(GNrot);
     6681                                xDelete(GErot);
     6682                        }
    66246683                }
    66256684        }
    … …  
    71757234}
    71767235/*}}}*/
    7177 void       Tria::SealevelchangeUpdateViscousFields(){ /*{{{*/
     7236void       Tria::SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){ /*{{{*/
    71787237
    71797238        /*Inputs:*/
    … …  
    71827241        IssmDouble* viscousN=NULL;
    71837242        IssmDouble* viscousE=NULL;
    7184         IssmDouble* viscoustimes=NULL;
    71857243        int         viscousnumsteps;
    7186         int         viscousindex=0;
    7187         int         newindex=0;
    71887244        int         dummy;
    71897245        bool        viscous=false;
    7190         IssmDouble  currenttime;
    7191         IssmDouble  lincoeff=0;
    7192         int horiz;
     7246        int         horiz=0;
    71937247
    71947248        this->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    … …  
    71967250        if(viscous){
    71977251                this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    7198 
    71997252                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    7200                 this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
    7201                 this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
    7202                 this->parameters->FindParam(&currenttime,TimeEnum);
    72037253
    72047254                this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
    … …  
    72097259                }
    72107260
    7211                 bool foundtime=false;
    7212                 int offset=1;
    7213                 lincoeff=0;
    7214                 newindex=viscousnumsteps-2;
    7215 
    7216                 for(int t=viscousindex;t<viscousnumsteps;t++){
    7217                         if (viscoustimes[t]>currenttime){
    7218                                 newindex=t-1;
    7219                                 lincoeff=(currenttime-viscoustimes[newindex])/(viscoustimes[t]-viscoustimes[newindex]);
    7220                                 foundtime=true;
    7221                                 offset=0;
    7222                                 break;
    7223                         }
    7224                 }
    7225 
    7226                 if(!foundtime) lincoeff=1;
    7227                 viscoustimes[newindex]=currenttime;
    72287261                for(int i=0;i<NUMVERTICES;i++){
    72297262                        viscousRSL[i*viscousnumsteps+newindex+offset]=(1-lincoeff)*viscousRSL[i*viscousnumsteps+newindex]+lincoeff*viscousRSL[i*viscousnumsteps+newindex+1];
    … …  
    72347267                        }
    72357268                }
    7236                 viscousindex=newindex+offset;
    7237 
    7238                 this->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum);
    7239                 this->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum);
    7240 
    7241                 /*free allocations:*/
    7242                 xDelete<IssmDouble>(viscoustimes);
     7269
    72437270        }
    72447271
    … …  
    72467273/*}}}*/
    72477274void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
     7275
     7276        int nel;
    72487277
    72497278        /*Inputs:*/
    … …  
    72517280        IssmDouble W[NUMVERTICES];
    72527281        IssmDouble BP[NUMVERTICES];
     7282        IssmDouble* areae=NULL;
    72537283
    72547284        /*output: */
    … …  
    72677297        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    72687298        this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
     7299        this->parameters->FindParam(&areae,&nel,AreaeEnum);
    72697300
    72707301        /*Retrieve inputs:*/
    … …  
    72917322
    72927323        /*Compute barystatic component in kg:*/
     7324        // Note: Iavg, etc, already include partial area factor phi for subelement loading
    72937325        bslcice =   -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg;
    72947326        bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
    … …  
    73157347                BPavg=0;
    73167348        }
    7317         /*Plug remaining values into centroi load vector:*/
     7349        /*Plug remaining values into centroid load vector:*/
    73187350        loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
    73197351
    … …  
    73277359        /*sal green function:*/
    73287360        IssmDouble* G=NULL;
     7361        IssmDouble* Grot=NULL;
    73297362        IssmDouble* Gsub[SLGEOM_NUMLOADS];
    73307363        bool computefuture=false;
    … …  
    73367369        int  size;
    73377370        int  nel,nbar;
    7338         IssmDouble Grotm1[3];
    7339         IssmDouble Grotm2[3];
    7340         IssmDouble Grotm3[3];
     7371
    73417372
    73427373        this->parameters->FindParam(&sal,SolidearthSettingsSelfAttractionEnum);
    … …  
    73507381                this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    73517382                this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
    7352 
    7353                 this->SealevelchangeGxL(sealevelpercpu, G, Gsub, loads, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    7354         }
    7355 
    7356         if(rotation){
    7357                 Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
    7358                 Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    7359                 Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    7360 
    7361                 for(int i=0;i<NUMVERTICES;i++){
    7362                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7363                                 sealevelpercpu[this->vertices[i]->lid]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2];
    7364                         }
    7365                 }
    7366         }
     7383                if (rotation)   this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
     7384
     7385                this->SealevelchangeGxL(sealevelpercpu, G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
     7386        }
     7387
    73677388        return;
    73687389} /*}}}*/
    73697390void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
    73707391
    7371         /*sal green function:*/
    73727392        IssmDouble oceanaverage=0;
    73737393        IssmDouble oceanarea=0;
    … …  
    74177437        IssmDouble* GE=NULL;
    74187438        IssmDouble* GN=NULL;
     7439        IssmDouble* Grot=NULL;
     7440        IssmDouble* GUrot=NULL;
     7441        IssmDouble* GNrot=NULL;
     7442        IssmDouble* GErot=NULL;
    74197443        IssmDouble* Gsub[SLGEOM_NUMLOADS];
    74207444        IssmDouble* GUsub[SLGEOM_NUMLOADS];
    … …  
    74257449        int horiz;
    74267450        int size;
    7427         IssmDouble Grotm1[3];
    7428         IssmDouble Grotm2[3];
    7429         IssmDouble Grotm3[3];
    7430         IssmDouble GUrotm1[3];
    7431         IssmDouble GUrotm2[3];
    7432         IssmDouble GUrotm3[3];
    7433         IssmDouble GNrotm1[3];
    7434         IssmDouble GNrotm2[3];
    7435         IssmDouble GNrotm3[3];
    7436         IssmDouble GErotm1[3];
    7437         IssmDouble GErotm2[3];
    7438         IssmDouble GErotm3[3];
     7451
    74397452        bool rotation= false;
    74407453        bool elastic=false;
    … …  
    74727485                                this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size);
    74737486                        }
    7474                 }
    7475                 this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
     7487                        if (rotation) {
     7488                                this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
     7489                                this->inputs->GetArrayPtr(SealevelchangeGUrotEnum,this->lid,&GUrot,&size);
     7490                                if (horiz){
     7491                                        this->inputs->GetArrayPtr(SealevelchangeGErotEnum,this->lid,&GErot,&size);
     7492                                        this->inputs->GetArrayPtr(SealevelchangeGNrotEnum,this->lid,&GNrot,&size);
     7493                                }
     7494                        }
     7495                }
     7496                this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
    74767497
    74777498                if(elastic){
    7478                         this->SealevelchangeGxL(&UGrd[0],GU, GUsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
    7479                         if(horiz ){
    7480                                 this->SealevelchangeGxL(&NGrd[0],GN, GNsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
    7481                                 this->SealevelchangeGxL(&EGrd[0],GE, GEsub, loads, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    7482                         }
    7483                 }
    7484         }
    7485 
    7486         if(rotation){
    7487                 Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
    7488                 Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
    7489                 Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
    7490 
    7491                 for(int i=0;i<NUMVERTICES;i++){
    7492                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7493                                 RSLGrd[i]+=Grotm1[i]*polarmotionvector[0]+Grotm2[i]*polarmotionvector[1]+Grotm3[i]*polarmotionvector[2];
    7494                         }
    7495                 }
    7496 
    7497                 if (elastic){
    7498                         Element::GetInputListOnVertices(&GUrotm1[0],SealevelGUrotm1Enum);
    7499                         Element::GetInputListOnVertices(&GUrotm2[0],SealevelGUrotm2Enum);
    7500                         Element::GetInputListOnVertices(&GUrotm3[0],SealevelGUrotm3Enum);
    7501 
    7502                         for(int i=0;i<NUMVERTICES;i++){
    7503                                 if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7504                                         UGrd[i]+=GUrotm1[i]*polarmotionvector[0]+GUrotm2[i]*polarmotionvector[1]+GUrotm3[i]*polarmotionvector[2];
    7505                                 }
    7506                         }
    7507                         if (horiz){
    7508                                 Element::GetInputListOnVertices(&GNrotm1[0],SealevelGNrotm1Enum);
    7509                                 Element::GetInputListOnVertices(&GNrotm2[0],SealevelGNrotm2Enum);
    7510                                 Element::GetInputListOnVertices(&GNrotm3[0],SealevelGNrotm3Enum);
    7511                                 Element::GetInputListOnVertices(&GErotm1[0],SealevelGErotm1Enum);
    7512                                 Element::GetInputListOnVertices(&GErotm2[0],SealevelGErotm2Enum);
    7513                                 Element::GetInputListOnVertices(&GErotm3[0],SealevelGErotm3Enum);
    7514 
    7515                                 for(int i=0;i<NUMVERTICES;i++){
    7516                                         if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7517                                                 NGrd[i]+=GNrotm1[i]*polarmotionvector[0]+GNrotm2[i]*polarmotionvector[1]+GNrotm3[i]*polarmotionvector[2];
    7518                                                 EGrd[i]+=GErotm1[i]*polarmotionvector[0]+GErotm2[i]*polarmotionvector[1]+GErotm3[i]*polarmotionvector[2];
    7519                                         }
    7520                                 }
     7499                        this->SealevelchangeGxL(&UGrd[0],GU, GUsub, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
     7500                        if(horiz){
     7501                                this->SealevelchangeGxL(&NGrd[0],GN, GNsub, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
     7502                                this->SealevelchangeGxL(&EGrd[0],GE, GEsub, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    75217503                        }
    75227504                }
    … …  
    75437525
    75447526} /*}}}*/
    7545 void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
     7527void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
     7528
     7529        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    75467530
    75477531        IssmDouble* grdfield=NULL;
    75487532        int i,e,l,t,nbar;
    75497533        bool computeviscous=false;
     7534        bool rotation=false;
    75507535        IssmDouble* viscousfield=NULL;
    7551         int nt=1; //important
     7536        int nt=1; //important, ensures there is a defined value if computeviscous is false
    75527537        int viscousindex=0; //important
    75537538        int viscousnumsteps=1; //important
    75547539
    75557540        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
     7541        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    75567542        if(computeviscous){
    75577543                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    7558                 if(computefuture) nt=viscousnumsteps;
     7544                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     7545                if(computefuture) {
     7546                        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
     7547                        if (nt>viscousnumsteps) nt=viscousnumsteps;
     7548                }
    75597549                else nt=1;
    7560 
    7561                 //allocate
    7562                 grdfield=xNewZeroInit<IssmDouble>(3*nt);
    7563         }
    7564         else grdfield=xNewZeroInit<IssmDouble>(3*nt);
    7565 
    7566         if(loads->sealevelloads){
    7567 
     7550        }
     7551        //allocate
     7552        grdfield=xNewZeroInit<IssmDouble>(3*nt);
     7553
     7554        if(rotation){ //add rotational feedback
     7555                for(int i=0;i<NUMVERTICES;i++){ //vertices
     7556                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     7557                                for (int m=0;m<3;m++){ //polar motion components
     7558                                        for(t=0;t<nt;t++){ //time
     7559                                                int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
     7560                                                grdfield[i*nt+t]+=Grot[index]*polarmotionvector[m];
     7561                                        }
     7562                                }
     7563                        }
     7564                }
     7565        }
     7566
     7567
     7568        if(loads->sealevelloads){ // general case: loads + sealevel loads
    75687569                for(i=0;i<NUMVERTICES;i++) {
    75697570                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    … …  
    75947595        }
    75957596        else{  //this is the initial convolution where only loads are provided
    7596 
    75977597                for(i=0;i<NUMVERTICES;i++) {
    75987598                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    … …  
    76157615        }
    76167616
    7617         if(computeviscous){
     7617        if(computeviscous){ /*{{{*/
     7618                // we need to do up to 3 things (* = only if computefuture)
     7619                // 1*: add new grdfield contribution to the viscous stack for future time steps
     7620                // 2: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time]
     7621                // 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
     7622
    76187623                IssmDouble* grdfieldinterp=NULL;
    76197624                IssmDouble* viscoustimes=NULL;
    … …  
    76227627                IssmDouble  timeacc;
    76237628
    7624                 this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
    76257629                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
    76267630                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    76277631                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
    76287632                this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
     7633                /* Map new grdfield generated by present-day loads onto viscous time vector*/
    76297634                if(computefuture){
    7630                         grdfieldinterp=xNew<IssmDouble>(3*nt);
     7635                        grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps);
     7636                        //viscousindex time and first time step of grdfield coincide, so just copy that value
     7637                        for(int i=0;i<NUMVERTICES;i++){
     7638                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7639                                grdfieldinterp[i*viscousnumsteps+viscousindex]=  grdfield[i*nt+0];
     7640                        }
    76317641                        if(viscoustimes[viscousindex]<final_time){
     7642                                //And interpolate the rest of the points in the future
    76327643                                lincoeff=(viscoustimes[viscousindex+1]-viscoustimes[viscousindex])/timeacc;
    7633                                 for(int t=viscousindex;t<nt;t++){
    7634                                         if(t==viscousindex){
    7635                                                 for(int i=0;i<NUMVERTICES;i++){
    7636                                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7637                                                         grdfieldinterp[i*nt+t]=  grdfield[i*nt+0];
    7638                                                 }
    7639                                         }
    7640                                         else{
    7641                                                 for(int i=0;i<NUMVERTICES;i++){
    7642                                                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7643                                                         grdfieldinterp[i*nt+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
    7644                                                 }
     7644                                for(int t=viscousindex+1;t<viscousnumsteps;t++){
     7645                                        for(int i=0;i<NUMVERTICES;i++){
     7646                                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7647                                                grdfieldinterp[i*viscousnumsteps+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
    76457648                                        }
    76467649                                }
    … …  
    76547657                }
    76557658
    7656                 if(computefuture){ /*update viscous stack with future deformation from present load: */
    7657 
    7658                         for(int t=nt-1;t>=viscousindex;t--){
     7659                /*update viscous stack with future deformation from present load: */
     7660                if(computefuture){
     7661                        for(int t=viscousnumsteps-1;t>=viscousindex;t--){
    76597662                                for(int i=0;i<NUMVERTICES;i++){
    76607663                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    76617664                                        //offset viscousfield to remove all deformation that has already been added
    7662                                         viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*nt+t]-grdfieldinterp[i*nt+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
     7665                                        viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]-grdfieldinterp[i*viscousnumsteps+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
    76637666                                }
    76647667                        }
    7665                         /*Re-add viscous stack now that we updated:*/
     7668                        /*Save viscous stack now that we updated the values:*/
    76667669                        this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps);
    76677670                }
    … …  
    76697672                /*Free allocatoins:*/
    76707673                xDelete<IssmDouble>(viscoustimes);
    7671         }
     7674        }
     7675        /*}}}*/
    76727676
    76737677        /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
    … …  
    76847688
    76857689} /*}}}*/
    7686 void       Tria::SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom){/*{{{*/
    7687 
    7688         IssmDouble S=0;
    7689 
    7690         /*Compute area of element:*/
    7691         IssmDouble area,planetarea,re;
    7692         IssmDouble late,longe;
    7693         this->Element::GetInputValue(&area,AreaEnum);
    7694 
    7695         /*recover parameters: */
    7696         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    7697         this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
    7698         late=slgeom->late[this->lid]/180*M_PI;
    7699         longe=slgeom->longe[this->lid]/180*M_PI;
    7700 
    7701         /*recover total load: */
    7702         if(loads->loads) S+=loads->loads[this->Sid()];
    7703         if(loads->sealevelloads) S+=loads->sealevelloads[this->Sid()];
    7704 
    7705         /* Perturbation terms for moment of inertia (moi_list):
    7706          * computed analytically (see Wu & Peltier, eqs 10 & 32)
    7707          * also consistent with my GMD formulation!
    7708          * ALL in geographic coordinates
    7709          * */
    7710         dI_list[0] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    7711         dI_list[1] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    7712         dI_list[2] = +4*M_PI*(S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    7713         return;
    7714 }/*}}}*/
    7715 void       Tria::SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* grdloads,  SealevelGeometry* slgeom){/*{{{*/
    7716 
    7717         IssmDouble  SA=0;
    7718         IssmDouble* loads=NULL;
    7719         IssmDouble* sealevelloads=NULL;
    7720         IssmDouble  late,longe,re;
    7721         int         intj;
    7722         IssmDouble  area;
    7723         IssmDouble  planetarea;
    7724 
    7725         /*recover parameters: */
    7726         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    7727         this->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
    7728 
    7729         /*Initalize:*/
    7730         for(int i=0;i<3;i++)dI_list[i]=0;
    7731 
    7732         /*Go through our loads:*/
    7733         for(int i=0;i<SLGEOM_NUMLOADS;i++){
    7734                 if(slgeom->issubelement[i][this->lid]){
    7735                         loads=grdloads->subloads[i];
    7736                         if(i==SLGEOM_OCEAN) sealevelloads=grdloads->subsealevelloads;
    7737                         intj=slgeom->subelementmapping[i][this->lid];
    7738                         late=slgeom->latbarycentre[i][intj]/180*M_PI;
    7739                         longe=slgeom->longbarycentre[i][intj]/180*M_PI;
    7740                         area=slgeom->area_subel[i][intj];
    7741 
    7742                         /*recover total load: */
    7743                         if(loads) SA+=loads[intj]*area;
    7744                         if(sealevelloads) SA+=sealevelloads[intj]*area;
    7745                 }
    7746         }
    7747 
    7748         /* Perturbation terms for moment of inertia (moi_list):
    7749          * computed analytically (see Wu & Peltier, eqs 10 & 32)
    7750          * also consistent with my GMD formulation!
    7751          * ALL in geographic coordinates
    7752          * */
    7753         dI_list[0] += -4*M_PI*(SA)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
    7754         dI_list[1] += -4*M_PI*(SA)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
    7755         dI_list[2] += +4*M_PI*(SA)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
    7756 
    7757         return;
    7758 }/*}}}*/
     7690
    77597691void       Tria::SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    77607692

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

    r26784 r26800  
    178178                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelpercpu, GrdLoads* loads, IssmDouble* rotationvector,SealevelGeometry* slgeom);
    179179                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom);
    180                 void       SealevelchangeMomentOfInertiaCentroid(IssmDouble* dI_list, GrdLoads* loads,  SealevelGeometry* slgeom);
    181                 void       SealevelchangeMomentOfInertiaSubElement(IssmDouble* dI_list, GrdLoads* loads, SealevelGeometry* slgeom);
    182                 void       SealevelchangeUpdateViscousFields();
     180                void       SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset);
    183181                #endif
    184182                /*}}}*/
    … …  
    244242                void           UpdateConstraintsExtrudeFromBase(void);
    245243                void           UpdateConstraintsExtrudeFromTop(void);
    246                 void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, GrdLoads* loads, SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
     244                void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
    247245                /*}}}*/
    248246

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

    r26468 r26800  
    6868
    6969} /*}}}*/
     70void GrdLoads::SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom){
     71
     72        IssmDouble area,re, S;
     73        int ylmindex, intj;
     74        IssmDouble deg2coeff_local[5];
     75
     76        femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
     77
     78        for (int c=0;c<5;c++) deg2coeff_local[c]=0;
     79
     80        for(Object* & object : femmodel->elements->objects){
     81                Element* element = xDynamicCast<Element*>(object);
     82                //first, loads on centroids
     83                S=0;
     84                S+=loads[element->Sid()];
     85
     86                if(sealevelloads) S+=sealevelloads[element->Sid()];
     87                if(S!=0){
     88                        element->Element::GetInputValue(&area,AreaEnum);
     89
     90                        for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients: 2,0; 2,1cos; 2,1sin; 2,2cos; 2,2sin
     91                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
     92                                deg2coeff_local[c] += S*area/re/re*slgeom->Ylm[ylmindex];
     93                        }
     94                }
     95                //add loads on subelement barycenters
     96                for (int i=0;i<SLGEOM_NUMLOADS;i++){
     97                        if (slgeom->issubelement[i][element->lid]){
     98                                intj=slgeom->subelementmapping[i][element->lid];
     99                                S=0;
     100                                S+=subloads[i][intj];
     101                                if(i==SLGEOM_OCEAN && sealevelloads) S+=subsealevelloads[intj];
     102                                if(S!=0){
     103                                        area=slgeom->area_subel[i][intj];
     104                                        for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients
     105                                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
     106                                                deg2coeff_local[c] += S*area/re/re*slgeom->Ylm_subel[i][ylmindex];
     107                                        }
     108                                }
     109                        }
     110                }
     111        }
     112
     113        //sum each degree 2 coefficient across all cpus to get global total
     114        ISSM_MPI_Reduce (&deg2coeff_local[0],&deg2coeff[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     115        ISSM_MPI_Bcast(&deg2coeff[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     116        ISSM_MPI_Reduce (&deg2coeff_local[1],&deg2coeff[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     117        ISSM_MPI_Bcast(&deg2coeff[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     118        ISSM_MPI_Reduce (&deg2coeff_local[2],&deg2coeff[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     119        ISSM_MPI_Bcast(&deg2coeff[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     120
     121}

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

    r26227 r26800  
    2929                void BroadcastLoads(void);
    3030                void BroadcastSealevelLoads(void);
     31                void SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom);
    3132};
    3233#endif  /* _SEALEVELGRDLOADS_H_ */

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

    r26468 r26800  
    3535                nsubel[i]=0;
    3636                nbar[i]=0;
     37                Ylm_subel[i]= xNewZeroInit<IssmDouble>(localnel*9);
    3738        }
    3839        late=xNew<IssmDouble>(localnel);
    … …  
    4041        isoceanin=xNew<bool>(localnel);
    4142        lids=xNew<int>(localnodsin);
     43        Ylm=xNewZeroInit<IssmDouble>(localnel*9); // (degmax+1)^2 terms, degmax=2
    4244
    4345}; /*}}}*/
    … …  
    5658                xDelete<IssmDouble>(longbarycentre[i]);
    5759                xDelete<IssmDouble>(area_subel[i]);
    58         }
     60                xDelete<IssmDouble>(Ylm_subel[i]);
     61        }
     62        xDelete<IssmDouble>(Ylm);
    5963        xDelete<IssmDouble>(late);
    6064        xDelete<IssmDouble>(longe);
    … …  
    155159
    156160} /*}}}*/
     161void SealevelGeometry::BuildSphericalHarmonics(){ /*{{{*/
     162        //builds spherical harmonics functions for degrees 0, 1, 2 on centroids/barycenters
     163        //0: used for global average
     164        //1: used for geocenter motion
     165        //2: used for rotational feedback
     166        int intj, count;
     167
     168        IssmDouble YlmNorm[9];
     169
     170        //YlmNormalization: N^2=(2*l+1)/4/pi * factorial(l-m)/factorial(l+m) if m==0
     171        //             : 2*N^2 if m>0
     172        // such that integral(Ylm * Ylm *YlmNorm dS) = 1 on the unit sphere.
     173        YlmNorm[0]=(0.25/M_PI); //Y00
     174        YlmNorm[1]=(0.75/M_PI); //Y10
     175        YlmNorm[2]=(0.75/M_PI); //Y11c
     176        YlmNorm[3]=YlmNorm[2];   //Y11s
     177        YlmNorm[4]=(1.25/M_PI); //Y20
     178        YlmNorm[5]=(1.25/3./M_PI); //Y21c
     179        YlmNorm[6]=YlmNorm[5]; //Y21s
     180        YlmNorm[7]=(1.25/12./M_PI); //Y22c
     181        YlmNorm[8]=YlmNorm[7]; //Y22s
     182
     183        for (int e=0;e<localnel;e++){
     184                IssmDouble lat=late[e]*M_PI/180.;
     185                IssmDouble lon=longe[e]*M_PI/180.;
     186                Ylm[0*localnel+e] = 1.0 *YlmNorm[0]; //Y00
     187
     188                Ylm[1*localnel+e] = sin(lat)*YlmNorm[1]; //Y10
     189                Ylm[2*localnel+e] = cos(lat)*cos(lon)*YlmNorm[2]; //Y11cos
     190                Ylm[3*localnel+e] = cos(lat)*sin(lon)*YlmNorm[3]; //Y11sin
     191
     192                //Ylm[4*localnel+e] = 0.25 - 0.75*cos(2.0*lat) ; //Y20
     193                Ylm[4*localnel+e] = (1.5*pow(sin(lat),2.)-0.5)*YlmNorm[4]; //Y20
     194                Ylm[5*localnel+e] = 1.5*sin(2.*lat)*cos(lon)*YlmNorm[5]; //Y21cos
     195                Ylm[6*localnel+e] = 1.5*sin(2.*lat)*sin(lon)*YlmNorm[6]; //Y21sin
     196                Ylm[7*localnel+e] = 1.5*(1.+cos(2.*lat))*cos(2.*lon)*YlmNorm[7]; //Y22cos
     197                Ylm[8*localnel+e] = 1.5*(1.+cos(2.*lat))*sin(2.*lon)*YlmNorm[8]; //Y22sin
     198        }
     199
     200        for (int i=0;i<SLGEOM_NUMLOADS;i++){
     201                for (int e=0;e<localnel;e++){
     202                        if (issubelement[i][e]){
     203                                intj=subelementmapping[i][e];
     204                                IssmDouble lat=latbarycentre[i][intj]*M_PI/180.;
     205                                IssmDouble lon=longbarycentre[i][intj]*M_PI/180.;
     206                                Ylm_subel[i][0*localnel+e] = 1.0*YlmNorm[0]; //Y00
     207
     208                                Ylm_subel[i][1*localnel+e] = sin(lat)*YlmNorm[1]; //Y10
     209                                Ylm_subel[i][2*localnel+e] = cos(lat)*cos(lon)*YlmNorm[2]; //Y11cos
     210                                Ylm_subel[i][3*localnel+e] = cos(lat)*sin(lon)*YlmNorm[3]; //Y11sin
     211
     212                                Ylm_subel[i][4*localnel+e] = (1.5*pow(sin(lat),2.)-0.5)*YlmNorm[4]; //Y20
     213                                Ylm_subel[i][5*localnel+e] = 1.5*sin(2.*lat)*cos(lon)*YlmNorm[5]; //Y21cos
     214                                Ylm_subel[i][6*localnel+e] = 1.5*sin(2.*lat)*sin(lon)*YlmNorm[6]; //Y21sin
     215                                Ylm_subel[i][7*localnel+e] = 1.5*(1.+cos(2.*lat))*cos(2.*lon)*YlmNorm[7]; //Y22cos
     216                                Ylm_subel[i][8*localnel+e] = 1.5*(1.+cos(2.*lat))*sin(2.*lon)*YlmNorm[8]; //Y22sin
     217                        }
     218                }
     219        }
     220} /*}}}*/

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

    r26469 r26800  
    1414
    1515#include "../toolkits/toolkits.h"
     16#include "../classes/classes.h"
    1617
    1718class SealevelGeometry{
    … …  
    3031                IssmDouble* late;
    3132                IssmDouble* longe;
     33                IssmDouble* Ylm;
     34                IssmDouble* Ylm_subel[SLGEOM_NUMLOADS];
     35                IssmDouble* YlmNorm[9];
    3236                bool* isoceanin;
    3337                bool*       issubelement[SLGEOM_NUMLOADS];
    … …  
    4549                int GNEnum(int l);
    4650                int GEEnum(int l);
     51                void BuildSphericalHarmonics(void);
    4752};
    4853#endif  /* _SEALEVELGEOMETRY_H_ */

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

    r26468 r26800  
    2020                int nyi;
    2121                int starting_layer;
     22                int* deg_layer_delete;
    2223
    2324                LoveVariables(){  /*{{{*/
    … …  
    2728                        nyi=0;
    2829                        starting_layer=0;
     30                        deg_layer_delete=NULL;
    2931                } /*}}}*/
    30                 LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin){  /*{{{*/
     32                LoveVariables(IssmDouble* EarthMassin,IssmDouble g0in,IssmDouble r0in,int nyiin,int starting_layerin, int* deg_layer_deletein){  /*{{{*/
    3133                        EarthMass=EarthMassin;
    3234                        g0=g0in;
    … …  
    3436                        nyi=nyiin;
    3537                        starting_layer=starting_layerin;
     38                        deg_layer_delete=deg_layer_deletein;
    3639                } /*}}}*/
    3740                ~LoveVariables(){};
     41}; /*}}}*/
     42
     43template <class doubletype> class LoveNumbers{  /*{{{*/
     44
     45        public:
     46                doubletype* H;
     47                doubletype* K;
     48                doubletype* L;
     49                doubletype* Kernels;
     50                int sh_nmin, sh_nmax, nfreq, nkernels;
     51
     52                LoveNumbers(){  /*{{{*/
     53                        H=NULL;
     54                        K=NULL;
     55                        L=NULL;
     56                        Kernels=NULL;
     57                        sh_nmin=0;
     58                        sh_nmax=0;
     59                        nfreq=0;
     60                        nkernels=0;
     61                } /*}}}*/
     62                LoveNumbers(int sh_nminin, int sh_nmaxin, int nfreqin, Matlitho* matlitho){  /*{{{*/
     63                        sh_nmax=sh_nmaxin;
     64                        sh_nmin=sh_nminin;
     65                        nfreq=nfreqin;
     66                        nkernels=(sh_nmax+1)*(matlitho->numlayers+1)*6;
     67                        H=xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
     68                        K=xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
     69                        L=xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
     70                        Kernels=xNewZeroInit<doubletype>(nfreq*nkernels);
     71                } /*}}}*/
     72                void DownCastToDouble(LoveNumbers<IssmDouble>* LoveDouble){
     73                        for(int i=0;i<(sh_nmax+1)*nfreq;i++){
     74                                LoveDouble->H[i]=std::real(H[i]);
     75                                LoveDouble->K[i]=std::real(K[i]);
     76                                LoveDouble->L[i]=std::real(L[i]);
     77                        }
     78                }
     79                void LoveMPI_Gather(LoveNumbers<doubletype>* Love_local, int lower_row){
     80                        int* recvcounts=xNew<int>(IssmComm::GetSize());
     81                        int* displs=xNew<int>(IssmComm::GetSize());
     82                        int  rc;
     83                        int  offset;
     84                        int nf_local = Love_local->nfreq;
     85
     86                        /*Deal H, K, L first, as they share the same size*/
     87                        rc=(sh_nmax+1)*nf_local;
     88                        offset=(sh_nmax+1)*lower_row;
     89                        ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     90                        ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     91                        ISSM_MPI_Allgatherv(Love_local->H, rc, ISSM_MPI_DOUBLE, H, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     92                        ISSM_MPI_Allgatherv(Love_local->K, rc, ISSM_MPI_DOUBLE, K, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     93                        ISSM_MPI_Allgatherv(Love_local->L, rc, ISSM_MPI_DOUBLE, L, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     94
     95                        /* deal with love kernels now */
     96                        rc=nf_local*nkernels;
     97                        offset=lower_row*nkernels;
     98                        ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     99                        ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     100                        ISSM_MPI_Allgatherv(Love_local->Kernels, rc, ISSM_MPI_DOUBLE, Kernels, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     101
     102                        xDelete<int>(recvcounts);
     103                        xDelete<int>(displs);
     104                }
     105                void Copy(LoveNumbers<doubletype>* LoveDup){
     106                        for(int i=0;i<(sh_nmax+1)*nfreq;i++){
     107                                H[i]=LoveDup->H[i];
     108                                K[i]=LoveDup->K[i];
     109                                L[i]=LoveDup->L[i];
     110                        }
     111                        for(int i=0;i<nkernels*nfreq;i++){
     112                                Kernels[i]=LoveDup->Kernels[i];
     113                        }
     114                }
     115                ~LoveNumbers(){
     116                        xDelete<doubletype>(H);
     117                        xDelete<doubletype>(K);
     118                        xDelete<doubletype>(L);
     119                        xDelete<doubletype>(Kernels);
     120                };
    38121}; /*}}}*/
    39122
    … …  
    123206        return xi;
    124207}/*}}}*/
    125 template<typename doubletype> void         postwidder_transform(doubletype* Lovet, doubletype* Lovef,int d, int t, int NTit, doubletype* xi, FemModel* femmodel){ /*{{{*/
     208template<typename doubletype> void         postwidder_transform(doubletype* Lovet, doubletype* Lovef,int d, int t, int sh_nmax,int NTit, doubletype* xi, FemModel* femmodel){ /*{{{*/
    126209        //Computes Lovet for time step t and degree d from the PW coefficients xi and the corresponding 2*NTit frequency samples in Lovef
    127210
    128         int nfreq, indxi, indf;
    129         femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
    130         int nt=nfreq/2/NTit;
    131 
    132         indf=d*nfreq+t*2*NTit;
     211        int indxi, indf;
     212        IssmDouble PW_test;
     213        IssmDouble PW_threshold;
     214        femmodel->parameters->FindParam(&PW_threshold,LovePostWidderThresholdEnum);
     215
     216        indf=(t*2*NTit)*(sh_nmax+1)+d;
    133217        doubletype* LoveM = xNew<doubletype>(NTit);
     218
     219
     220        // test variation across frequencies tested, something with little frequency dependence is not worth going through PW tranform
     221        // that transform would also be numerically unstable
     222        PW_test = abs((Lovef[indf+(2*NTit-1)*(sh_nmax+1)]-Lovef[indf])/Lovef[indf]);
     223
     224        //if (PW_test < PW_threshold){ //elastic or fluid response: Love(t) = Love(s), we can do an early return
     225        //      Lovet[t*(sh_nmax+1)+d]=Lovef[indf];
     226        //      return;
     227        //}
     228
     229        if (PW_test==0){ //elastic or fluid response: Love(t) = Love(s), we can do an early return
     230                Lovet[t*(sh_nmax+1)+d]=Lovef[indf];
     231                return;
     232        }
    134233
    135234        for (int M=1;M<NTit+1;M++){
    … …  
    137236                for (int k=1;k<2*M+1;k++){
    138237                        indxi=(M-1)*(2*NTit)+k-1;
    139                         LoveM[M-1]+=xi[indxi]*Lovef[indf+k-1];
     238                        LoveM[M-1]+=xi[indxi]*Lovef[indf+(k-1)*(sh_nmax+1)];
    140239                }
    141240
    … …  
    146245                        if ( abs(LoveM[M-1]-LoveM[M-2]) > abs(LoveM[M-2]-LoveM[M-3]) &&
    147246                                        abs(LoveM[M-2]-LoveM[M-3]) > abs(LoveM[M-3]-LoveM[M-4]) ){
    148                                 Lovet[d*nt+t]=LoveM[M-3];
     247                                Lovet[t*(sh_nmax+1)+d]=LoveM[M-3];
    149248                                return;
    150249                        }
    151250                }
    152251        }
    153         Lovet[d*nt+t]=LoveM[NTit-1];
     252        Lovet[t*(sh_nmax+1)+d]=LoveM[NTit-1];
    154253}/*}}}*/
    155254template <typename doubletype> void        GetEarthRheology(doubletype* pla, doubletype* pmu, int layer_index, doubletype omega,  Matlitho* matlitho){ /*{{{*/
    … …  
    726825
    727826}/*}}}*/
    728 template <typename doubletype> void        yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars){ /*{{{*/
     827template <typename doubletype> void        yi_boundary_conditions(doubletype* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type){ /*{{{*/
    729828
    730829        IssmDouble  g0,r0,mu0,ra,rb,rc;
    731         int nyi, forcing_type,icb,cmb,starting_layer;
     830        int nyi,icb,cmb,starting_layer;
    732831        IssmDouble* EarthMass;
    733832
    … …  
    739838
    740839        femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
    741         femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
    742840        femmodel->parameters->FindParam(&icb,LoveInnerCoreBoundaryEnum);
    743841        femmodel->parameters->FindParam(&cmb,LoveCoreMantleBoundaryEnum);
    … …  
    810908        }
    811909}/*}}}*/
    812 template <typename doubletype> void        solve_yi_system(doubletype* loveh, doubletype* lovel, doubletype* lovek, int deg, doubletype omega, doubletype* yi, doubletype* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars){ /*{{{*/
     910template <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){ /*{{{*/
    813911
    814912        IssmDouble  g0,r0,mu0,loveratio,underflow_tol;
    815         IssmDouble* frequencies;
     913        //IssmDouble* frequencies;
    816914        int nyi,starting_layer, dummy;
    817915        bool allow_layer_deletion;
    … …  
    827925        femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
    828926        femmodel->parameters->FindParam(&underflow_tol,LoveUnderflowTolEnum);
    829         femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
     927        //femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum);
    830928        IssmDouble ra=matlitho->radius[matlitho->numlayers];
    831929        bool exit=false;
    832930        int lda,ldb;
     931
    833932
    834933        for(;!exit;){ //cycles of: attempt to solve the yi system, then delete a layer if necessary
    … …  
    854953                xDelete<int>(ipiv);
    855954
    856                 if(VerboseModule() && info!=0){
    857                         _printf0_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
     955                        /*_printf_("i j yi[i+nyi*j] rhs[i]");
     956                        for (int i=0;i<nyi;i++){
     957                                        _printf_(i<<" "<<rhs[i]<<"\n");
     958                        }
     959
     960                        for (int i=0;i<nyi;i++){
     961                                for (int j=0;j<nyi;j++){
     962                                        _printf_(i<<" "<<j<<" "<<yi[i+nyi*j]<<" "<<rhs[i]<<"\n");
     963                                }
     964                        }
     965                        _error_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");*/
     966
     967                if(VerboseSolution() && info!=0){
    858968                        _printf_("i j yi[i+nyi*j] rhs[i]");
    859969                        for (int i=0;i<nyi;i++){
    … …  
    862972                                }
    863973                        }
     974                        _error_("love core warning in DGESV : LAPACK linear equation solver couldn't resolve the system");
    864975                }
    865976
    … …  
    880991                if (abs(lovek1/lovek1s) < loveratio) loveratio = abs(lovek1/lovek1s);
    881992
    882                 if (!allow_layer_deletion || nyi<=12 || omega!=angular_frequency<doubletype>(frequencies[0])){
     993               
     994
     995                if (!allow_layer_deletion || nyi<=12 || omega!=angular_frequency<doubletype>(frequencies[0]) || deg==0){
    883996                        goto save_results;
    884997                        /*We are not allowed to delete layers, or there is only one layer left. We also don't want to delete
    … …  
    8871000                }
    8881001
    889                 if (loveratio<=underflow_tol || xIsNan(loveratio) || xIsInf(loveratio)){//We need to delete a layer and try again if the ratio between deepest love number to surface love number is too low (risk of underflow) or garbage
    890                         if(VerboseSolution()){
    891                                 _printf_("\n   Degree: " << deg <<", surface/Depth Love number ratio small: "<<loveratio<<"\n");
    892                                 _printf_("    Changing the interface where the integration starts\n");
    893                                 _printf_("    New start interface: r="<< matlitho->radius[starting_layer+1] <<"m\n\n");
    894                         }
     1002                if (omega==0){ // if running elastic love_numbers, record at which degree we must delete layers, this way we synch layer deletion between cpus next time we calculate love numbers
     1003                        //We need to delete a layer and try again if the ratio between deepest love number to surface love number is too low (risk of underflow) or garbage
     1004                        if (loveratio<=underflow_tol || xIsNan(loveratio) || xIsInf(loveratio)) {
     1005                                vars->deg_layer_delete[starting_layer]=deg;
     1006                                if(VerboseSolution() && verbosecpu){
     1007                                        _printf_("\n   Degree: " << deg <<", surface/Depth Love number ratio small: "<<loveratio<<"\n");
     1008                                        _printf_("    Changing the interface where integration starts\n");
     1009                                        _printf_("    New start interface: r="<< matlitho->radius[starting_layer+1] <<"m\n\n");
     1010                                }
     1011                        }
     1012                }
     1013
     1014                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
     1015                        if (omega!=0 && VerboseSolution()  && verbosecpu) _printf_(", deleting layer " << starting_layer << "\n");
    8951016                        nyi-=6;
    8961017                        starting_layer+=1;
    … …  
    9341055                xDelete<doubletype>(rhslocal);
    9351056        }
    936         xDelete<IssmDouble>(frequencies);       
     1057        //xDelete<IssmDouble>(frequencies);     
    9371058
    9381059}/*}}}*/
    939 template <typename doubletype> void        love_freq_to_temporal(doubletype* LoveHt,doubletype* LoveLt,doubletype* LoveKt,doubletype* LoveHf,doubletype* LoveLf,doubletype* LoveKf, FemModel* femmodel){ /*{{{*/
     1060template <typename doubletype> void        love_freq_to_temporal(LoveNumbers<doubletype>* Lovet, LoveNumbers<doubletype>* Tidalt, doubletype* pmtf_colineart, doubletype* pmtf_orthot, LoveNumbers<doubletype>* Lovef,LoveNumbers<doubletype>* Tidalf, IssmDouble* frequencies, FemModel* femmodel, bool verbosecpu){ /*{{{*/
    9401061        //Transforms all frequency-dependent love numbers into time-dependent love numbers
    941         int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit;
    942 
    943         femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
    944         femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
    945         femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
     1062        int nfreq,sh_nmax,sh_nmin,indxi,indf, NTit, forcing_type, nt;
     1063        IssmDouble kf,Omega,moi_e,moi_p,alpha;
     1064        doubletype* pmtf_colinearf=NULL;
     1065        doubletype* pmtf_orthof=NULL;
     1066        bool chandler_wobble=false;
     1067
     1068        nfreq=Lovef->nfreq;
     1069        sh_nmin=Lovef->sh_nmin;
     1070        sh_nmax=Lovef->sh_nmax;
     1071
    9461072        femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
    9471073
    948         int nt=nfreq/2/NTit;
     1074        //Parameters for the rotationnal feedback
     1075        femmodel->parameters->FindParam(&chandler_wobble,LoveChandlerWobbleEnum);
     1076        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
     1077        femmodel->parameters->FindParam(&kf,TidalLoveK2SecularEnum);
     1078        femmodel->parameters->FindParam(&Omega,RotationalAngularVelocityEnum);
     1079        femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
     1080        femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
     1081
     1082        nt=Lovet->nfreq;
    9491083        doubletype* xi=postwidder_coef<doubletype>(NTit);
     1084
     1085        if(VerboseSolution()  && verbosecpu) _printf_("   Inverse Laplace Transform... ");
    9501086
    9511087        for (int d=sh_nmin;d<sh_nmax+1;d++){
    9521088                for (int t=0;t<nt;t++){
    953                         postwidder_transform<doubletype>(LoveHt,LoveHf,d,t,NTit,xi,femmodel);
    954                         postwidder_transform<doubletype>(LoveLt,LoveLf,d,t,NTit,xi,femmodel);
    955                         postwidder_transform<doubletype>(LoveKt,LoveKf,d,t,NTit,xi,femmodel);
    956                 }
    957         }
     1089                        postwidder_transform<doubletype>(Lovet->H,Lovef->H,d,t,sh_nmax,NTit,xi,femmodel);
     1090                        postwidder_transform<doubletype>(Lovet->K,Lovef->K,d,t,sh_nmax,NTit,xi,femmodel);
     1091                        postwidder_transform<doubletype>(Lovet->L,Lovef->L,d,t,sh_nmax,NTit,xi,femmodel);
     1092                }
     1093        }
     1094
     1095        if(VerboseSolution()  && verbosecpu) _printf_("done!\n");
     1096
     1097        if (forcing_type==11){ //Let's retrieve the functions necessary for the rotational_feedback
     1098                if(VerboseSolution()  && verbosecpu) _printf_("     Transforming PMTF and tidal love numbers... ");
     1099                pmtf_colinearf = xNewZeroInit<doubletype>(3*nfreq);
     1100                pmtf_orthof = xNewZeroInit<doubletype>(3*nfreq);
     1101                int d=2;
     1102                doubletype s,R1,R2;
     1103                alpha=(moi_p-moi_e)/moi_e; //Earth flattening
     1104
     1105                if (chandler_wobble){ //Chandler Wobble is untested yet
     1106                        for (int fr=0;fr<nfreq;fr++){           
     1107                                s=angular_frequency<doubletype>(frequencies[fr]);
     1108                                R1=alpha*Omega*(1.0-Tidalf->K[fr*3+d]/kf);
     1109                                R2=1.0+alpha*Tidalf->K[fr*3+d]/kf;
     1110                                pmtf_colinearf[fr*3+d]=alpha*(1.0+Lovef->K[fr*(sh_nmax+1)+d])*(Omega*R1-pow(s,2)*R2)/(pow(R1,2.0)+pow(s*R2,2.0));
     1111                                pmtf_orthof[fr*3+d]=alpha*(1.0+Lovef->K[fr*(sh_nmax+1)+d])*s*Omega*(1.0+alpha)/(pow(R1,2.0)+pow(s*R2,2.0));
     1112                        }
     1113                }
     1114                else {
     1115                        for (int fr=0;fr<nfreq;fr++){           
     1116                                pmtf_colinearf[fr*3+d]=(1.0+Lovef->K[fr*(sh_nmax+1)+d])/(1.0-Tidalf->K[fr*3+d]/kf);
     1117                                pmtf_orthof[fr*3+d]=0.0;
     1118                        }
     1119                }
     1120                for (int t=0;t<nt;t++){
     1121                        postwidder_transform<doubletype>(Tidalt->H,Tidalf->H,2,t,2,NTit,xi,femmodel);
     1122                        postwidder_transform<doubletype>(Tidalt->K,Tidalf->K,2,t,2,NTit,xi,femmodel);
     1123                        postwidder_transform<doubletype>(Tidalt->L,Tidalf->L,2,t,2,NTit,xi,femmodel);
     1124                        postwidder_transform<doubletype>(pmtf_colineart,pmtf_colinearf,2,t,2,NTit,xi,femmodel);
     1125                        postwidder_transform<doubletype>(pmtf_orthot,pmtf_orthof,2,t,2,NTit,xi,femmodel);
     1126                }
     1127                if(VerboseSolution() && verbosecpu) _printf_("done!\n");
     1128        }
     1129
    9581130        xDelete<doubletype>(xi);
     1131        xDelete<doubletype>(pmtf_colinearf);
     1132        xDelete<doubletype>(pmtf_orthof);
     1133}/*}}}*/
     1134
     1135template <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){
     1136
     1137        int nstep, kernel_index,kernel_indexe,deleted_layer_offset, deg, sh_nmin, sh_nmax, nfreq;
     1138        doubletype  lovek, loveh, lovel, loveratio;
     1139        doubletype  omega;
     1140        doubletype* yi_prefactor=NULL;
     1141        doubletype* yi_righthandside=NULL;
     1142        doubletype* yi=NULL;
     1143        IssmDouble  underflow_tol;
     1144        bool freq_skip, istemporal;
     1145
     1146        femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
     1147        femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
     1148
     1149        nfreq=Lovef->nfreq;
     1150        sh_nmin=Lovef->sh_nmin;
     1151        sh_nmax=Lovef->sh_nmax;
     1152        if (Elastic==NULL) sh_cutoff=sh_nmax;
     1153
     1154        // reset deleted layers in case we have called this function before;
     1155        vars->starting_layer=0;
     1156        vars->nyi=6*(matlitho->numlayers+1);
     1157
     1158        yi_prefactor=xNewZeroInit<doubletype>(6*6*nstep*matlitho->numlayers);
     1159        yi_righthandside=xNewZeroInit<doubletype>(vars->nyi);
     1160        yi=xNewZeroInit<doubletype>(vars->nyi*vars->nyi);
     1161
     1162        //precalculate yi coefficients that do not depend on degree or frequency
     1163        fill_yi_prefactor<doubletype>(yi_prefactor, NULL, NULL,femmodel, matlitho,vars);
     1164
     1165        if (VerboseSolution() && Elastic  && verbosecpu) _printf_("\n");
     1166
     1167        for(int deg=0;deg<2;deg++){ // calculation is in the center of mass reference frame, neutralize degree 0 and 1 mass changes
     1168                for (int fr=0;fr<nfreq;fr++){
     1169                        Lovef->K[fr*(sh_nmax+1)+deg]=-1.0;
     1170                }
     1171        }
     1172
     1173        for(int deg=sh_nmin;deg<sh_cutoff+1;deg++){
     1174                if (VerboseSolution() && Elastic && verbosecpu) {
     1175                        _printf_("\r   Degree: " << deg << "/" << sh_nmax << "    ");
     1176                }
     1177
     1178                //precalculate yi coefficients that depend on degree but not frequency
     1179                fill_yi_prefactor<doubletype>(yi_prefactor, &deg, NULL,femmodel, matlitho,vars);
     1180
     1181                for (int fr=0;fr<nfreq;fr++){
     1182                        omega=angular_frequency<doubletype>(frequencies[fr]);
     1183                       
     1184                        //precalculate yi coefficients that depend on degree and frequency
     1185                        fill_yi_prefactor<doubletype>(yi_prefactor, &deg,&omega,femmodel, matlitho,vars);
     1186
     1187                        //solve the system
     1188                        yi_boundary_conditions<doubletype>(yi_righthandside,deg,femmodel,matlitho,vars,forcing_type);
     1189                        build_yi_system<doubletype>(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars);
     1190                        solve_yi_system<doubletype>(&loveh,&lovel,&lovek, deg, omega, frequencies, yi, yi_righthandside,femmodel, matlitho,vars,verbosecpu);
     1191                        Lovef->H[fr*(sh_nmax+1)+deg]=loveh;
     1192                        Lovef->K[fr*(sh_nmax+1)+deg]=lovek-1.0;
     1193                        Lovef->L[fr*(sh_nmax+1)+deg]=lovel;
     1194                        deleted_layer_offset=(matlitho->numlayers+1)*6-vars->nyi;// =6 per deleted layer
     1195                        kernel_index=fr*(sh_nmax+1)*(matlitho->numlayers+1)*6 + deg*(matlitho->numlayers+1)*6 + deleted_layer_offset;
     1196                        for (int i=0;i<vars->nyi;i++){
     1197                                Lovef->Kernels[kernel_index+i]=yi_righthandside[i];
     1198                        }
     1199                }
     1200        }
     1201
     1202        if (Elastic) { // if elastic values were provided, we copy elastic love numbers above the cutoff degree instead of calculating them
     1203                for(int deg=sh_cutoff+1;deg<sh_nmax+1;deg++){
     1204                        if (VerboseSolution() && Elastic  && verbosecpu) {
     1205                                if (deg==sh_nmax || deg%100==0) _printf_("\r   Degree: " << deg << "/" << Lovef->sh_nmax << "    ");
     1206                        }
     1207                        for (int fr=0;fr<nfreq;fr++){
     1208                                // just copy the elastic values
     1209                                Lovef->H[fr*(sh_nmax+1)+deg]=Elastic->H[deg];
     1210                                Lovef->K[fr*(sh_nmax+1)+deg]=Elastic->K[deg];
     1211                                Lovef->L[fr*(sh_nmax+1)+deg]=Elastic->L[deg];
     1212                                deleted_layer_offset=(matlitho->numlayers+1)*6-vars->nyi;// =6 per deleted layer
     1213                                kernel_index=fr*(sh_nmax+1)*(matlitho->numlayers+1)*6 + deg*(matlitho->numlayers+1)*6 + deleted_layer_offset;
     1214                                kernel_indexe=deg*(matlitho->numlayers+1)*6 + deleted_layer_offset;
     1215                                for (int i=0;i<vars->nyi;i++){
     1216                                        Lovef->Kernels[kernel_index+i]=Elastic->Kernels[kernel_indexe+i];
     1217                                }
     1218                        }
     1219                }
     1220        }
     1221
     1222        if (VerboseSolution() && Elastic  && verbosecpu) _printf_("\n");
     1223        xDelete<doubletype>(yi);
     1224        xDelete<doubletype>(yi_righthandside);
     1225        xDelete<doubletype>(yi_prefactor);
    9591226}/*}}}*/
    9601227
    … …  
    9731240        IssmDouble  g0,r0;
    9741241        int         nyi,starting_layer;
     1242        int*        deg_layer_delete;
    9751243
    9761244        femmodel->parameters->FindParam(&GG,LoveGravitationalConstantEnum);
    9771245        EarthMass=xNewZeroInit<IssmDouble>(numlayers+1);
     1246        deg_layer_delete=xNewZeroInit<int>(numlayers);
    9781247
    9791248        for (int i=0;i<numlayers;i++){
    … …  
    9931262        starting_layer=0;
    9941263
    995         return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer);
     1264        if(VerboseSolution()){
     1265                _printf_("     Surface gravity: " << g0 << " m.s^-2\n");
     1266                _printf_("     Mean density: " << EarthMass[numlayers-1]/(4.0/3.0*PI*pow(r0,3.0)) << " kg.m^-3\n");
     1267        }
     1268
     1269        return new LoveVariables(EarthMass,g0,r0,nyi,starting_layer,deg_layer_delete);
    9961270
    9971271} /*}}}*/
    … …  
    9991273
    10001274        Matlitho*   matlitho=NULL;
    1001         int         nfreq,nyi,starting_layer,nstep,forcing_type,dummy;
     1275        int         nfreq, NTit,nt, forcing_type,dummy, sh_cutoff;
    10021276        int         sh_nmin,sh_nmax,kernel_index,deleted_layer_offset;
    1003         bool        allow_layer_deletion,love_kernels, istemporal;
     1277        bool        allow_layer_deletion,love_kernels, istemporal, freq_skip;
    10041278        bool        verbosemod = (int)VerboseModule();
    1005         IssmDouble  mu0, underflow_tol;
    10061279        IssmDouble *frequencies = NULL;
     1280        IssmDouble *frequencies_local=NULL;
    10071281        bool        save_results;
    10081282        bool        complex_computation;
     1283        bool        verbosecpu=false;
    10091284
    10101285        doubletype  omega;
    10111286        doubletype  lovek, loveh, lovel, loveratio;
    1012 
    1013         doubletype*  LoveKf = NULL;
    1014         doubletype*  LoveHf = NULL;
    1015         doubletype*  LoveLf = NULL;
    1016 
    1017         doubletype*  LoveKernels= NULL;
     1287        IssmDouble pw_threshold, pw_test_h, pw_test_l,pw_test_k;
     1288
     1289        LoveNumbers<doubletype>* Lovef=NULL;
     1290        LoveNumbers<doubletype>* Tidalf=NULL;
     1291
     1292        /* parallel computing */
     1293        LoveNumbers<doubletype>* Lovef_local=NULL;
     1294        LoveNumbers<doubletype>* Tidalf_local=NULL;
     1295
     1296        /*elastic & fluid love numbers*/
     1297        LoveNumbers<doubletype>* Elastic=NULL;
     1298        IssmDouble* frequencies_elastic=NULL;
     1299        LoveNumbers<doubletype>* Fluid=NULL;
     1300        IssmDouble* frequencies_fluid=NULL;
     1301
    10181302        LoveVariables* vars=NULL;
    1019         doubletype* yi_prefactor=NULL;
    1020         doubletype* yi_righthandside=NULL;
    1021         doubletype* yi=NULL;
    1022 
    1023         if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");
    10241303
    10251304        /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
    … …  
    10391318        femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
    10401319        femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
    1041         femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
    10421320        femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
    10431321        femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
    10441322        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
    10451323        femmodel->parameters->FindParam(&istemporal,LoveIsTemporalEnum);
    1046         femmodel->parameters->FindParam(&nstep,LoveIntStepsPerLayerEnum);
    10471324        femmodel->parameters->FindParam(&complex_computation,LoveComplexComputationEnum);
    1048 
    1049         /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
    1050         LoveKf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
    1051         LoveHf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
    1052         LoveLf = xNewZeroInit<doubletype>(nfreq*(sh_nmax+1));
     1325        femmodel->parameters->FindParam(&pw_threshold,LovePostWidderThresholdEnum);
     1326        if (istemporal) femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
     1327
     1328        /*Initialize three love matrices: geoid, vertical and horizontal displacement*/
     1329        Lovef= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nfreq, matlitho);
     1330        Tidalf= new LoveNumbers<doubletype>(2,2,nfreq, matlitho);
     1331        Elastic= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,matlitho);
     1332        Fluid= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,1,matlitho);
    10531333
    10541334        /*Initialize love kernels (real and imaginary parts): */
    1055         LoveKernels= xNewZeroInit<doubletype>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
    1056 
    10571335        vars=love_init(femmodel,matlitho);
    10581336
    1059         yi_prefactor=xNewZeroInit<doubletype>(6*6*nstep*matlitho->numlayers);
    1060         yi_righthandside=xNewZeroInit<doubletype>(vars->nyi);
    1061         yi=xNewZeroInit<doubletype>(vars->nyi*vars->nyi);
    1062 
    1063         //precalculate yi coefficients that do not depend on degree or frequency
    1064         fill_yi_prefactor<doubletype>(yi_prefactor, NULL, NULL,femmodel, matlitho,vars);
    1065         for(int deg=sh_nmin;deg<sh_nmax+1;deg++){
    1066 
    1067                 //precalculate yi coefficients that depend on degree but not frequency
    1068                 fill_yi_prefactor<doubletype>(yi_prefactor, &deg, NULL,femmodel, matlitho,vars);
    1069 
    1070                 for (int fr=0;fr<nfreq;fr++){
    1071                         omega=angular_frequency<doubletype>(frequencies[fr]);
    1072 
    1073                         //precalculate yi coefficients that depend on degree and frequency
    1074                         fill_yi_prefactor<doubletype>(yi_prefactor, &deg,&omega,femmodel, matlitho,vars);
    1075 
    1076                         yi_boundary_conditions<doubletype>(yi_righthandside,deg,femmodel,matlitho,vars);
    1077 
    1078                         build_yi_system<doubletype>(yi,deg,omega,yi_prefactor,femmodel,matlitho,vars);
    1079 
    1080                         solve_yi_system<doubletype>(&loveh,&lovel,&lovek, deg, omega, yi, yi_righthandside,femmodel, matlitho,vars);
    1081 
    1082                         LoveHf[deg*nfreq+fr]=loveh;
    1083                         LoveKf[deg*nfreq+fr]=lovek-1.0;
    1084                         LoveLf[deg*nfreq+fr]=lovel;
    1085 
    1086                         deleted_layer_offset=(matlitho->numlayers+1)*6-vars->nyi;// =6 per deleted layer
    1087                         kernel_index=deg*nfreq*(matlitho->numlayers+1)*6 + fr*(matlitho->numlayers+1)*6 + deleted_layer_offset;
    1088                         for (int i=0;i<vars->nyi;i++){
    1089                                 LoveKernels[kernel_index+i]=yi_righthandside[i];
    1090                         }
    1091                 }
    1092         }
    1093 
    1094         xDelete<doubletype>(yi);
    1095         xDelete<doubletype>(yi_righthandside);
     1337        //distribute frequencies for parallel computation /*{{{*/
     1338        int nt_local, nf_local, lower_row, upper_row;
     1339        if (istemporal){
     1340                //temporal love numbers are obtained via blocks of 2*NTit frequencies samples
     1341                //here we are making sure no block is split between different cpus, which would make the inverse laplace transform impossible
     1342                nt=nfreq/2/NTit;
     1343                nt_local=DetermineLocalSize(nt,IssmComm::GetComm());
     1344                nf_local=nt_local*2*NTit; // number of local frequencies
     1345                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,nt_local,IssmComm::GetComm());
     1346                lower_row*=2*NTit;
     1347                upper_row*=2*NTit;
     1348        }
     1349        else{
     1350                //in this case frequency samples are completely independent so we can split them evenly across cpus
     1351                nf_local=DetermineLocalSize(nfreq,IssmComm::GetComm());
     1352                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,nf_local,IssmComm::GetComm());
     1353        }
     1354
     1355        if (lower_row==0) verbosecpu=true; //let only cpu1 be verbose
     1356        if(VerboseSolution() && verbosecpu) _printf0_("   computing LOVE numbers\n");
     1357
     1358        frequencies_local=xNewZeroInit<IssmDouble>(nf_local);
     1359        for (int fr=0;fr<nf_local;fr++) frequencies_local[fr]=frequencies[lower_row+fr];
     1360
     1361        Lovef_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nf_local, matlitho);
     1362        Tidalf_local= new LoveNumbers<doubletype>(2,2,nf_local, matlitho);
     1363
     1364        /*}}}*/
     1365        frequencies_elastic=xNewZeroInit<IssmDouble>(1);
     1366        frequencies_fluid=xNewZeroInit<IssmDouble>(1);
     1367        for (int fr=0;fr<nfreq;fr++){ // find the lowest non-zero frequency requested
     1368                if (frequencies_fluid[0]==0) frequencies_fluid[0]=frequencies[fr];
     1369                else if(frequencies[fr]!=0 && frequencies_fluid[0]>frequencies[fr]) frequencies_fluid[0]=frequencies[fr];
     1370        }
     1371
     1372        // run elastic and fluid love numbers
     1373        if(VerboseSolution() && verbosecpu) _printf_("     elastic\n");
     1374        compute_love_numbers<doubletype>(Elastic, NULL, forcing_type, sh_nmax,frequencies_elastic, femmodel, matlitho, vars,verbosecpu);
     1375
     1376        if (nfreq>1){
     1377                compute_love_numbers<doubletype>(Fluid, NULL, forcing_type, sh_nmax,frequencies_fluid, femmodel, matlitho, vars,verbosecpu);
     1378                sh_cutoff=sh_nmax;
     1379                for (int deg=100;deg<sh_nmax+1;deg++){
     1380                        pw_test_h=abs((Fluid->H[deg]-Elastic->H[deg])/Elastic->H[deg]);
     1381                        pw_test_k=abs((Fluid->K[deg]-Elastic->K[deg])/Elastic->K[deg]);
     1382                        pw_test_l=abs((Fluid->L[deg]-Elastic->L[deg])/Elastic->L[deg]);
     1383                        if (pw_test_h<pw_threshold && pw_test_k<pw_threshold && pw_test_l<pw_threshold){
     1384                                sh_cutoff=deg;
     1385                                if(VerboseSolution() && verbosecpu){
     1386                                        _printf_("   Degree: " << deg << "/" << sh_nmax << "    ");
     1387                                        _printf_("      found negligible variation across frequencies, will copy elastic values after this degree\n");
     1388                                        _printf_("      Delta_h/h=" << pw_test_h << "; Delta_k/k="<< pw_test_k << "; Delta_l/l=" << pw_test_l << "; threshold set to " << pw_threshold << "\n");
     1389                                }
     1390                                break;
     1391                        }
     1392                }
     1393        }
     1394        else sh_cutoff=sh_nmax;
     1395               
     1396
     1397
     1398        //Take care of rotationnal feedback love numbers first, if relevant /*{{{*/
     1399        if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){ // if forcing is surface loading and we have degree 2
     1400                if(VerboseSolution() && verbosecpu) _printf_("     tidal\n");
     1401                int tidal_forcing_type=9;
     1402                compute_love_numbers<doubletype>(Tidalf_local, NULL,tidal_forcing_type=9, 2,frequencies_local, femmodel, matlitho, vars,verbosecpu);
     1403        }
     1404        /*}}}*/
     1405
     1406        //Resume requested forcing_type
     1407        if (nfreq>1){ // if we are not running just elastic love numbers
     1408                if(VerboseSolution() && verbosecpu){
     1409                        if (forcing_type==11) _printf_("     loading\n");
     1410                        else if(forcing_type==9) _printf_("     tidal\n");
     1411                        else _printf_("     love\n");
     1412                }
     1413                compute_love_numbers<doubletype>(Lovef_local, Elastic, forcing_type, sh_cutoff, frequencies_local, femmodel, matlitho, vars,verbosecpu);
     1414        }
     1415        else{
     1416                Lovef_local->Copy(Elastic);
     1417        }
     1418        /*}}}*/
     1419
    10961420
    10971421        //Temporal love numbers
    10981422        if (istemporal && !complex_computation){
    1099 
    1100                 IssmDouble*  LoveHtDouble=NULL;
    1101                 IssmDouble*  LoveKtDouble=NULL;
    1102                 IssmDouble*  LoveLtDouble=NULL;
    1103                 doubletype*  LoveHt=NULL;
    1104                 doubletype*  LoveLt=NULL;
    1105                 doubletype*  LoveKt=NULL;
    1106 
    1107                 int NTit;
    1108                 femmodel->parameters->FindParam(&NTit,LoveNTemporalIterationsEnum);
    1109                 int nt = nfreq/2/NTit;
    1110 
    1111                 LoveHt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    1112                 LoveLt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    1113                 LoveKt=xNewZeroInit<doubletype>((sh_nmax+1)*nt);
    1114 
    1115                 love_freq_to_temporal<doubletype>(LoveHt,LoveLt,LoveKt,LoveHf,LoveLf,LoveKf,femmodel);
    1116 
    1117                 /*Downcast and add into parameters:*/
    1118                 LoveHtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
    1119                 LoveLtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
    1120                 LoveKtDouble=xNew<IssmDouble>((sh_nmax+1)*nt);
    1121                 for(int i=0;i<(sh_nmax+1)*nt;i++){
    1122                         LoveHtDouble[i]=std::real(LoveHt[i]);
    1123                         LoveLtDouble[i]=std::real(LoveLt[i]);
    1124                         LoveKtDouble[i]=std::real(LoveKt[i]);
    1125                 }
     1423                /*Initialize*/
     1424                doubletype*  pmtf_colineart=NULL;
     1425                doubletype*  pmtf_orthot=NULL;
     1426                LoveNumbers<doubletype>* Lovet=NULL;
     1427                LoveNumbers<doubletype>* Tidalt=NULL;
     1428                /*Downcast arrays to be exported in parameters*/
     1429                LoveNumbers<IssmDouble>* LovetDouble=NULL;
     1430                LoveNumbers<IssmDouble>* TidaltDouble=NULL;
     1431                IssmDouble*  pmtf_colineartDouble=NULL;
     1432                IssmDouble*  pmtf_orthotDouble=NULL;
     1433                /* parallel computing */
     1434                doubletype*  pmtf_colineart_local=NULL;
     1435                doubletype*  pmtf_orthot_local=NULL;   
     1436                LoveNumbers<doubletype>* Lovet_local=NULL;
     1437                LoveNumbers<doubletype>* Tidalt_local=NULL;     
     1438
     1439                Lovet= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nt,matlitho);
     1440                Tidalt= new LoveNumbers<doubletype>(2,2,nt,matlitho);
     1441                pmtf_colineart=xNewZeroInit<doubletype>(3*nt);
     1442                pmtf_orthot=xNewZeroInit<doubletype>(3*nt);     
     1443
     1444                Lovet_local= new LoveNumbers<doubletype>(sh_nmin,sh_nmax,nt_local,matlitho);
     1445                Tidalt_local= new LoveNumbers<doubletype>(2,2,nt_local,matlitho);       
     1446                pmtf_colineart_local=xNewZeroInit<doubletype>(3*nt_local);
     1447                pmtf_orthot_local=xNewZeroInit<doubletype>(3*nt_local);
     1448
     1449                love_freq_to_temporal<doubletype>(Lovet_local,Tidalt_local,pmtf_colineart_local,pmtf_orthot_local,Lovef_local,Tidalf_local,frequencies_local,femmodel,verbosecpu);
     1450
     1451                /* MPI Gather */ /*{{{*/
     1452                Lovef->LoveMPI_Gather(Lovef_local, lower_row);
     1453                if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){
     1454                        Tidalf->LoveMPI_Gather(Tidalf_local, lower_row);               
     1455                }
     1456                Lovet->LoveMPI_Gather(Lovet_local, lower_row/2/NTit);
     1457                Tidalt->LoveMPI_Gather(Tidalt_local, lower_row/2/NTit);
     1458                //pmtf:
     1459                int* recvcounts=xNew<int>(IssmComm::GetSize());
     1460                int* displs=xNew<int>(IssmComm::GetSize());
     1461                int  rc;
     1462                int  offset;
     1463                rc=3*nt_local;
     1464                offset=3*lower_row/2/NTit;
     1465                ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     1466                ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     1467                ISSM_MPI_Allgatherv(pmtf_colineart_local, rc, ISSM_MPI_DOUBLE, pmtf_colineart, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     1468                ISSM_MPI_Allgatherv(pmtf_orthot_local, rc, ISSM_MPI_DOUBLE, pmtf_orthot, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
     1469                xDelete<int>(recvcounts);
     1470                xDelete<int>(displs);
     1471                /*}}}*/
     1472
     1473                /*Downcast and add into parameters:*/ /*{{{*/
     1474                LovetDouble= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nt,matlitho);
     1475                TidaltDouble= new LoveNumbers<IssmDouble>(2,2,nt,matlitho);     
     1476
     1477                pmtf_colineartDouble=xNew<IssmDouble>(nt);
     1478                pmtf_orthotDouble=xNew<IssmDouble>(nt);
     1479
     1480                Lovet->DownCastToDouble(LovetDouble);
     1481                Tidalt->DownCastToDouble(TidaltDouble);
     1482                for(int i=0;i<nt;i++){
     1483                        pmtf_colineartDouble[i]=std::real(pmtf_colineart[2*nt+i]);
     1484                        pmtf_orthotDouble[i]=std::real(pmtf_orthot[2*nt+i]);
     1485                }
     1486
    11261487                if(forcing_type==9){ //tidal loading
    1127                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
    1128                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
    1129                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
     1488                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1));
     1489                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1));
     1490                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1));
    11301491                }
    11311492                else if(forcing_type==11){ //surface loading
    1132                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHtDouble,(sh_nmax+1)*nt,1));
    1133                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKtDouble,(sh_nmax+1)*nt,1));
    1134                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLtDouble,(sh_nmax+1)*nt,1));
    1135                 }
    1136                 xDelete<IssmDouble>(LoveHtDouble);
    1137                 xDelete<IssmDouble>(LoveKtDouble);
    1138                 xDelete<IssmDouble>(LoveLtDouble);
    1139 
     1493                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovetDouble->H,(sh_nmax+1)*nt,1));
     1494                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovetDouble->K,(sh_nmax+1)*nt,1));
     1495                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovetDouble->L,(sh_nmax+1)*nt,1));
     1496                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,TidaltDouble->H,3*nt,1));
     1497                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,TidaltDouble->K,3*nt,1));
     1498                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,TidaltDouble->L,3*nt,1));
     1499                        femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,pmtf_colineartDouble,nt,1));
     1500                        femmodel->parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,pmtf_orthotDouble,nt,1));
     1501                }
     1502
     1503                xDelete<IssmDouble>(pmtf_colineartDouble);
     1504                xDelete<IssmDouble>(pmtf_orthotDouble);
     1505                /*}}}*/
     1506       
    11401507                /*Add into external results, no need to downcast, we can handle complexes/quad precision: */
    1141                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKrEnum,LoveKt,sh_nmax+1,nt,0,0));
    1142                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHrEnum,LoveHt,sh_nmax+1,nt,0,0));
    1143                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLrEnum,LoveLt,sh_nmax+1,nt,0,0));
    1144 
    1145                 xDelete<doubletype>(LoveHt);
    1146                 xDelete<doubletype>(LoveLt);
    1147                 xDelete<doubletype>(LoveKt);
     1508                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKtEnum,Lovet->K,nt,sh_nmax+1,0,0));
     1509                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHtEnum,Lovet->H,nt,sh_nmax+1,0,0));
     1510                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLtEnum,Lovet->L,nt,sh_nmax+1,0,0));
     1511
     1512                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalKtEnum,Tidalt->K,nt,3,0,0));
     1513                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalHtEnum,Tidalt->H,nt,3,0,0));
     1514                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveTidalLtEnum,Tidalt->L,nt,3,0,0));
     1515                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LovePMTF1tEnum,pmtf_colineart,nt,3,0,0));
     1516                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LovePMTF2tEnum,pmtf_orthot,nt,3,0,0));
     1517
     1518                delete Lovet;
     1519                delete Tidalt;
     1520                delete Lovet_local;
     1521                delete Tidalt_local;
     1522                delete LovetDouble;
     1523                delete TidaltDouble;
     1524
     1525                xDelete<doubletype>(pmtf_colineart);
     1526                xDelete<doubletype>(pmtf_orthot);
    11481527        }
    11491528        else{
    1150 
    1151                 IssmDouble*  LoveHfDouble=NULL;
    1152                 IssmDouble*  LoveKfDouble=NULL;
    1153                 IssmDouble*  LoveLfDouble=NULL;
     1529                LoveNumbers<IssmDouble>* LovefDouble=NULL;
     1530                LovefDouble= new LoveNumbers<IssmDouble>(sh_nmin,sh_nmax,nfreq,matlitho);
     1531
     1532                /*MPI_Gather*/
     1533                if (nfreq>1){
     1534                        Lovef->LoveMPI_Gather(Lovef_local, lower_row);
     1535                        if (forcing_type==11 && sh_nmin<=2 && sh_nmax>=2){
     1536                                Tidalf->LoveMPI_Gather(Tidalf_local, lower_row);               
     1537                        }
     1538                }
     1539                else{
     1540                        Lovef->Copy(Elastic);
     1541                        Tidalf->Copy(Tidalf_local);
     1542                }
    11541543
    11551544                /*Add into parameters:*/
    1156                 LoveHfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
    1157                 LoveLfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
    1158                 LoveKfDouble=xNew<IssmDouble>((sh_nmax+1)*nfreq);
    1159                 for(int i=0;i<(sh_nmax+1)*nfreq;i++){
    1160                         LoveHfDouble[i]=std::real(LoveHf[i]);
    1161                         LoveLfDouble[i]=std::real(LoveLf[i]);
    1162                         LoveKfDouble[i]=std::real(LoveKf[i]);
    1163                 }
     1545                Lovef->DownCastToDouble(LovefDouble);
    11641546                if(forcing_type==9){ //tidal loading
    1165                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
    1166                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
    1167                         femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
     1547                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,LovefDouble->H,(sh_nmax+1)*nfreq,1));
     1548                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,LovefDouble->K,(sh_nmax+1)*nfreq,1));
     1549                        femmodel->parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,LovefDouble->L,(sh_nmax+1)*nfreq,1));
    11681550                }
    11691551                else if(forcing_type==11){ //surface loading
    1170                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LoveHfDouble,(sh_nmax+1)*nfreq,1));
    1171                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LoveKfDouble,(sh_nmax+1)*nfreq,1));
    1172                         femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LoveLfDouble,(sh_nmax+1)*nfreq,1));
    1173                 }
    1174                 xDelete<IssmDouble>(LoveHfDouble);
    1175                 xDelete<IssmDouble>(LoveKfDouble);
    1176                 xDelete<IssmDouble>(LoveLfDouble);
     1552                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,LovefDouble->H,(sh_nmax+1)*nfreq,1));
     1553                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,LovefDouble->K,(sh_nmax+1)*nfreq,1));
     1554                        femmodel->parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,LovefDouble->L,(sh_nmax+1)*nfreq,1));
     1555                }
     1556                delete LovefDouble;
    11771557        }
    11781558
    11791559        /*Add into external results:*/
    1180         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKiEnum,LoveKf,sh_nmax+1,nfreq,0,0));
    1181         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHiEnum,LoveHf,sh_nmax+1,nfreq,0,0));
    1182         femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLiEnum,LoveLf,sh_nmax+1,nfreq,0,0));
    1183 
     1560        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKfEnum,Lovef->K,nfreq,sh_nmax+1,0,0));
     1561        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveHfEnum,Lovef->H,nfreq,sh_nmax+1,0,0));
     1562        femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveLfEnum,Lovef->L,nfreq,sh_nmax+1,0,0));
    11841563        /*Only when love_kernels is on*/
    11851564        if (love_kernels==1) {
    1186                 femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,LoveKernels,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
    1187         }
    1188 
     1565                femmodel->results->AddObject(new GenericExternalResult<doubletype*>(femmodel->results->Size()+1,LoveKernelsEnum,Lovef->Kernels,nfreq,(sh_nmax+1)*(matlitho->numlayers+1)*6,0,0));
     1566        }
    11891567        /*Free resources:*/
    11901568        xDelete<IssmDouble>(frequencies);
    1191         xDelete<doubletype>(LoveKf);
    1192         xDelete<doubletype>(LoveHf);
    1193         xDelete<doubletype>(LoveLf);
    1194         xDelete<doubletype>(LoveKernels);
    1195 
     1569        xDelete<IssmDouble>(frequencies_local);
     1570        xDelete<IssmDouble>(frequencies_elastic);
     1571        delete Lovef;
     1572        delete Lovef_local;
     1573        delete Tidalf;
     1574        delete Tidalf_local;
     1575        delete Elastic;
    11961576        /* Legacy for fortran core, to be removed after complete validation */
    11971577
    … …  
    12251605
    12261606        /*Add love matrices to results:*/
    1227         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
    1228         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
    1229         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
    1230         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
    1231         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
    1232         //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
     1607        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKtEnum,LoveKr,sh_nmax+1,nfreq,0,0));
     1608        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHtEnum,LoveHr,sh_nmax+1,nfreq,0,0));
     1609        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLtEnum,LoveLr,sh_nmax+1,nfreq,0,0));
     1610        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKfEnum,LoveKi,sh_nmax+1,nfreq,0,0));
     1611        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHfEnum,LoveHi,sh_nmax+1,nfreq,0,0));
     1612        //femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLfEnum,LoveLi,sh_nmax+1,nfreq,0,0));
    12331613
    12341614        //xDelete<IssmDouble>(LoveKr);
    … …  
    12511631template void        GetEarthRheology<IssmDouble>(IssmDouble* pla, IssmDouble* pmu, int layer_index, IssmDouble omega,  Matlitho* matlitho);
    12521632template void        GetEarthRheology<IssmComplex>(IssmComplex* pla, IssmComplex* pmu, int layer_index, IssmComplex omega,  Matlitho* matlitho);
    1253 template void        yi_boundary_conditions<IssmDouble>(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
    1254 template void        yi_boundary_conditions<IssmComplex>(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
     1633template void        yi_boundary_conditions<IssmDouble>(IssmDouble* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type);
     1634template void        yi_boundary_conditions<IssmComplex>(IssmComplex* yi_righthandside, int deg, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars, int forcing_type);
    12551635template void        yi_derivatives<IssmDouble>(IssmDouble* dydx, IssmDouble* y, int layer_index, int n, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    12561636template void        yi_derivatives<IssmComplex>(IssmComplex* dydx, IssmComplex* y, int layer_index, int n, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho);
    … …  
    12631643template void        build_yi_system<IssmDouble>(IssmDouble* yi, int deg, IssmDouble omega, IssmDouble* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
    12641644template void        build_yi_system<IssmComplex>(IssmComplex* yi, int deg, IssmComplex omega, IssmComplex* yi_prefactor, FemModel* femmodel, Matlitho* matlitho,LoveVariables* vars);
    1265 template void        solve_yi_system<IssmDouble>(IssmDouble* loveh, IssmDouble* lovel, IssmDouble* lovek, int deg, IssmDouble omega, IssmDouble* yi, IssmDouble* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
    1266 template void        solve_yi_system<IssmComplex>(IssmComplex* loveh, IssmComplex* lovel, IssmComplex* lovek, int deg, IssmComplex omega, IssmComplex* yi, IssmComplex* rhs, FemModel* femmodel, Matlitho* matlitho, LoveVariables* vars);
     1645template 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);
     1646template 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);
     1647template 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);
     1648template 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);
    12671649template IssmDouble  factorial<IssmDouble>(int n);
    12681650template IssmDouble* postwidder_coef<IssmDouble>(int NTit);
    12691651template IssmDouble  n_C_r<IssmDouble>(int n, int r);
    1270 template void         postwidder_transform<IssmDouble>(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int NTit, IssmDouble* xi, FemModel* femmodel);
     1652template void         postwidder_transform<IssmDouble>(IssmDouble* Lovet, IssmDouble* Lovef,int d, int t, int sh_nmax,int NTit, IssmDouble* xi, FemModel* femmodel);
    12711653
    12721654/*}}}*/

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

    r26468 r26800  
    2424bool slcconvergence(Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    2525IssmDouble  SealevelloadsOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble totaloceanarea);
    26 void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom);
     26void PolarMotion(IssmDouble* m, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom, bool computefuture);
     27void SealevelchangeUpdateViscousTimeSeries(FemModel* femmodel);
    2728void ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom);
    2829void ivins_deformation_core(FemModel* femmodel);
    … …  
    247248        bool rotation=false;
    248249        bool planethasocean=false;
     250        bool computefuture=false;
    249251        IssmDouble*           sealevelpercpu=NULL;
    250252
    … …  
    304306        if(VerboseSolution()) _printf0_("         starting  GRD convolutions\n");
    305307
    306         /*update viscous RSL:*/
    307         for(Object* & object : femmodel->elements->objects){
    308                 Element* element = xDynamicCast<Element*>(object);
    309                 element->SealevelchangeUpdateViscousFields();
    310         }
     308        /*update viscous time series to keep up with time stepping:*/
     309        SealevelchangeUpdateViscousTimeSeries(femmodel);
    311310
    312311        /*buildup loads: */
    … …  
    319318        loads->BroadcastLoads();
    320319
    321         //compute polar motion:
    322         PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom);
    323 
    324320        /*skip computation of sea level equation if requested, which means sea level loads should be zeroed */
    325321        if(!sealevelloading){
    326322                loads->sealevelloads=xNewZeroInit<IssmDouble>(nel);
    327323                loads->subsealevelloads=xNewZeroInit<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);
    328 
     324                PolarMotion(&polarmotionvector[0],femmodel,loads, slgeom, computefuture=true);
    329325                goto deformation;
    330326        }
    … …  
    332328        if(VerboseSolution()) _printf0_("         converging GRD convolutions\n");
    333329        for(;;){
     330
     331                //compute polar motion:
     332                PolarMotion(&polarmotionvector[0],femmodel,loads,slgeom,computefuture=false);
    334333
    335334                oldsealevelloads=loads->vsealevelloads->Duplicate(); loads->vsealevelloads->Copy(oldsealevelloads);
    … …  
    338337                for(Object* & object : femmodel->elements->objects){
    339338                        Element* element = xDynamicCast<Element*>(object);
    340                         element->SealevelchangeConvolution(sealevelpercpu, loads , polarmotionvector,slgeom);
     339                        element->SealevelchangeConvolution(sealevelpercpu, loads, polarmotionvector,slgeom);
    341340                }
    342341
    … …  
    364363                loads->BroadcastSealevelLoads();
    365364
    366                 //compute polar motion:
    367                 PolarMotion(&polarmotionvector[0],femmodel,loads, slgeom);
    368 
    369365                //convergence?
    370366                if(slcconvergence(loads->vsealevelloads,oldsealevelloads,eps_rel,eps_abs))break;
    … …  
    374370                iterations++;
    375371        }
     372
     373        //recompute polar motion one final time, this time updating viscous stacks for future time steps
     374        if (viscous)    PolarMotion(&polarmotionvector[0],femmodel,loads, slgeom, computefuture=true);
    376375
    377376        deformation:
    … …  
    418417
    419418        if (rotation) {
    420                 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,polarmotionvector[0],step,time));
    421                 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,polarmotionvector[1],step,time));
    422                 femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,polarmotionvector[2],step,time));
     419                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelchangePolarMotionXEnum,polarmotionvector[0],step,time));
     420                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelchangePolarMotionYEnum,polarmotionvector[1],step,time));
     421                femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelchangePolarMotionZEnum,polarmotionvector[2],step,time));
    423422        }
    424423
    … …  
    666665        }
    667666
     667        /*Compute spherical harmonic functions for spatial integrations of the loads*/
     668        slgeom->BuildSphericalHarmonics();
     669
    668670        femmodel->parameters->AddObject(new DoubleVecParam(XxeEnum,xxe,nel));
    669671        femmodel->parameters->AddObject(new DoubleVecParam(YyeEnum,yye,nel));
    … …  
    731733
    732734        vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
    733         vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas);
     735        vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas); 
    734736
    735737        vsealevelloadsvolume->Sum(&sealevelloadsaverage);
    … …  
    740742        return (sealevelloadsaverage+subsealevelloadsaverage)/totaloceanarea;
    741743} /*}}}*/
    742 void PolarMotion(IssmDouble* polarmotionvector, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom){ /*{{{*/
    743 
    744         IssmDouble  moi_list[3]={0,0,0};
    745         IssmDouble  moi_list_sub[3]={0,0,0};
    746         IssmDouble  moi_list_cpu[3]={0,0,0};
    747         IssmDouble*     tide_love_h  = NULL;
    748         IssmDouble*     tide_love_k  = NULL;
    749         IssmDouble*     load_love_k  = NULL;
    750         IssmDouble  tide_love_k2secular;
    751         IssmDouble  moi_e, moi_p;
    752         IssmDouble      m1, m2, m3;
     744void PolarMotion(IssmDouble* polarmotionvector, FemModel* femmodel,GrdLoads* loads, SealevelGeometry* slgeom, bool computefuture){ /*{{{*/
     745        //The purpose of this routine is to get the polar motion vector m=(m1, m2, m3) induced by the GrdLoads
     746        IssmDouble  S2coef[3];
     747        IssmDouble*     pmtf_col= NULL;
     748        IssmDouble*     pmtf_ortho   = NULL;
     749        IssmDouble*     pmtf_z   = NULL;
     750        IssmDouble* m1=NULL;
     751        IssmDouble* m2=NULL;
     752        IssmDouble* m3=NULL;
     753        IssmDouble* m1interp=NULL;
     754        IssmDouble* m2interp=NULL;
     755        IssmDouble* m3interp=NULL;
     756
     757        IssmDouble  moi_e, moi_p, re;
     758        IssmDouble mhprefactor, mzprefactor;
    753759        bool rotation=false;
     760        bool viscous=false;
     761        int nt=1;
     762
     763        IssmDouble* viscoustimes=NULL;
     764        IssmDouble* viscouspolarmotion=NULL;
     765        int         viscousnumsteps;
     766        int         viscousindex=0;
     767        int         dummy;
     768        IssmDouble  currenttime, final_time, lincoeff, timeacc;
    754769
    755770        /*early return?:*/
    … …  
    758773
    759774        /*retrieve parameters: */
    760         femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
    761         femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
    762         femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
    763         femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
     775        femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
    764776        femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
    765777        femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
    766 
    767         for(Object* & object : femmodel->elements->objects){
    768                 Element* element = xDynamicCast<Element*>(object);
    769 
    770                 element->SealevelchangeMomentOfInertiaCentroid(&moi_list[0],loads,slgeom);
    771                 element->SealevelchangeMomentOfInertiaSubElement(&moi_list_sub[0],loads, slgeom);
    772 
    773                 moi_list_cpu[0] += moi_list[0]+moi_list_sub[0];
    774                 moi_list_cpu[1] += moi_list[1]+moi_list_sub[1];
    775                 moi_list_cpu[2] += moi_list[2]+moi_list_sub[2];
    776         }
    777 
    778         for (int i=0;i<3;i++) moi_list[i]=0;
    779 
    780         ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    781         ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    782 
    783         ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    784         ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    785 
    786         ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    787         ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    788 
    789         /*compute perturbation terms for angular velocity vector: */
    790         m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
    791         m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
    792         m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
     778        femmodel->parameters->FindParam(&pmtf_col,NULL,SealevelchangePolarMotionTransferFunctionColinearEnum);
     779        femmodel->parameters->FindParam(&pmtf_ortho,NULL,SealevelchangePolarMotionTransferFunctionOrthogonalEnum);
     780        femmodel->parameters->FindParam(&pmtf_z,NULL,SealevelchangePolarMotionTransferFunctionZEnum);
     781        femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
     782
     783        if (viscous){
     784                femmodel->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     785                femmodel->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     786                femmodel->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
     787                femmodel->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
     788                femmodel->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     789                femmodel->parameters->FindParam(&currenttime,TimeEnum);
     790                femmodel->parameters->FindParam(&viscouspolarmotion,NULL,NULL,SealevelchangeViscousPolarMotionEnum);
     791
     792                if (computefuture){
     793                        nt = viscousnumsteps;
     794                        m1interp=xNewZeroInit<IssmDouble>(nt);
     795                        m2interp=xNewZeroInit<IssmDouble>(nt);
     796                        m3interp=xNewZeroInit<IssmDouble>(nt);
     797                }
     798        }
     799
     800        m1=xNewZeroInit<IssmDouble>(nt);
     801        m2=xNewZeroInit<IssmDouble>(nt);
     802        m3=xNewZeroInit<IssmDouble>(nt);
     803
     804        //Isolate degree 2 load coefficients
     805        for (int i=0;i<3;i++) S2coef[i]=0;
     806        loads->SHDegree2Coefficients(&S2coef[0],femmodel,slgeom);
     807
     808        //compute present (& future) polar motion from present loads
     809        mhprefactor=-4.*M_PI/5.*pow(re,4.)/(moi_p-moi_e);
     810        mzprefactor=8.*M_PI/15.*pow(re,4.)/moi_p;
     811        for (int tprime=0;tprime<nt;tprime++){
     812                m1[tprime] = mhprefactor * (pmtf_col[tprime] * S2coef[1] + pmtf_ortho[tprime] * S2coef[2]); //x-component
     813                m2[tprime] = mhprefactor * (pmtf_col[tprime] * S2coef[2] - pmtf_ortho[tprime] * S2coef[1]); //y-component
     814                m3[tprime] = mzprefactor *  pmtf_z[tprime]   * S2coef[0];                                   //z-component
     815        }
     816
     817        if(viscous){
     818                // we need to do up to 3 things (* = only if computefuture)
     819                // 1*: add new PM contribution to the viscous stack for future time steps
     820                // 2: collect viscous PM from past loads due at present-day and add it to PM[current_time]
     821                // 3*: subtract from viscous stack PM that has already been accounted for so we don't add it again at the next time step
     822                if(computefuture){
     823                        if(viscoustimes[viscousindex]<final_time){
     824                                lincoeff=(viscoustimes[viscousindex+1]-viscoustimes[viscousindex])/timeacc;
     825                                for(int t=viscousindex;t<nt;t++){ //we resynchronize m from the relative time above to the absolute time where t=0 <=> beginning of the simulation
     826                                        if(t==viscousindex){
     827                                                m1interp[t]=  m1[0];
     828                                                m2interp[t]=  m2[0];
     829                                                m3interp[t]=  m3[0];
     830                                        }
     831                                        else{ //we reinterpolate PM on viscoustimes, so we can handle the case where we are running with adaptative/uneven time steps
     832                                                int tprime=t-viscousindex-1;
     833                                                m1interp[t]=  (1.0-lincoeff)*m1[tprime]+lincoeff*m1[tprime+1];
     834                                                m2interp[t]=  (1.0-lincoeff)*m2[tprime]+lincoeff*m2[tprime+1];
     835                                                m3interp[t]=  (1.0-lincoeff)*m3[tprime]+lincoeff*m3[tprime+1];
     836                                        }
     837                                }
     838                        }
     839                }
     840                /*update PM at present time using viscous stack at present time: */
     841                m1[0]+=viscouspolarmotion[0*nt+viscousindex];
     842                m2[0]+=viscouspolarmotion[1*nt+viscousindex];
     843                m3[0]+=viscouspolarmotion[2*nt+viscousindex];
     844                if(computefuture){ /*update viscous stack with future deformation from present load: */
     845                        for(int t=nt-1;t>=viscousindex;t--){
     846                                //offset viscous PM to remove all deformation that has already been added
     847                                viscouspolarmotion[0*nt+t]+=m1interp[t]-m1interp[viscousindex]-viscouspolarmotion[0*nt+viscousindex];
     848                                viscouspolarmotion[1*nt+t]+=m2interp[t]-m2interp[viscousindex]-viscouspolarmotion[1*nt+viscousindex];
     849                                viscouspolarmotion[2*nt+t]+=m3interp[t]-m3interp[viscousindex]-viscouspolarmotion[2*nt+viscousindex];
     850                        }
     851                        // save updated viscous PM
     852                        femmodel->parameters->SetParam(viscouspolarmotion,viscousnumsteps,3,SealevelchangeViscousPolarMotionEnum);
     853                }
     854        }
     855       
    793856
    794857        /*Assign output pointers:*/
    795         polarmotionvector[0]=m1;
    796         polarmotionvector[1]=m2;
    797         polarmotionvector[2]=m3;
     858        polarmotionvector[0]=m1[0];
     859        polarmotionvector[1]=m2[0];
     860        polarmotionvector[2]=m3[0];
     861
     862        /*Free allocations:*/
     863        xDelete<IssmDouble>(m1);
     864        xDelete<IssmDouble>(m2);
     865        xDelete<IssmDouble>(m3);
     866        if (viscous){
     867                if (computefuture){
     868                        xDelete<IssmDouble>(m1interp);
     869                        xDelete<IssmDouble>(m2interp);
     870                        xDelete<IssmDouble>(m3interp);
     871                }
     872        }
     873
    798874} /*}}}*/
     875void       SealevelchangeUpdateViscousTimeSeries(FemModel* femmodel){ /*{{{*/
     876       
     877        IssmDouble* viscouspolarmotion=NULL;
     878        IssmDouble* viscoustimes=NULL;
     879        int         viscousnumsteps;
     880        int         viscousindex=0;
     881        int         newindex=0;
     882        int         dummy;
     883        bool        viscous=false;
     884        bool        rotation=false;
     885        IssmDouble  currenttime;
     886        IssmDouble  lincoeff=0;
     887               
     888        femmodel->parameters->FindParam(&viscous,SolidearthSettingsViscousEnum);
     889        femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     890       
     891        if(viscous){
     892                femmodel->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
     893                femmodel->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
     894                femmodel->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     895                femmodel->parameters->FindParam(&currenttime,TimeEnum);
     896                if (rotation)   femmodel->parameters->FindParam(&viscouspolarmotion,NULL,NULL,SealevelchangeViscousPolarMotionEnum);
     897
     898                bool foundtime=false;
     899                int offset=1; //handles the egde case where time found = max time in viscoustimes
     900                lincoeff=0;
     901                newindex=viscousnumsteps-2;
     902
     903                for(int t=viscousindex;t<viscousnumsteps;t++){
     904                        if (viscoustimes[t]>=currenttime){
     905                                newindex=t-1;
     906                                foundtime=true;
     907                                lincoeff=(currenttime-viscoustimes[newindex])/(viscoustimes[t]-viscoustimes[newindex]);
     908                                offset=0;
     909                                break;
     910                        }
     911                }
     912
     913                if(rotation){
     914                        int index=0;
     915                        for (int i=0;i<3;i++){
     916                                index=i*viscousnumsteps+newindex;
     917                                viscouspolarmotion[index+offset]=(1-lincoeff)*viscouspolarmotion[index]+lincoeff*viscouspolarmotion[index+1];
     918                        }
     919                        femmodel->parameters->SetParam(viscouspolarmotion,viscousnumsteps,3,SealevelchangeViscousPolarMotionEnum);
     920                }
     921
     922
     923                /*update viscous inputs:*/
     924                for(Object* & object : femmodel->elements->objects){
     925                        Element* element = xDynamicCast<Element*>(object);
     926                        element->SealevelchangeUpdateViscousFields(lincoeff,newindex,offset);
     927                }
     928
     929                viscoustimes[newindex]=currenttime;
     930                viscousindex=newindex+offset;
     931
     932                femmodel->parameters->SetParam(viscousindex,SealevelchangeViscousIndexEnum);
     933                femmodel->parameters->SetParam(viscoustimes,viscousnumsteps,SealevelchangeViscousTimesEnum);
     934
     935                /*free allocations:*/
     936                xDelete<IssmDouble>(viscoustimes);
     937        }
     938
     939
     940}
    799941void        ConserveOceanMass(FemModel* femmodel,GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    800942
    … …  
    809951IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel){ /*{{{*/
    810952
     953        //merges loads on centroids and subelements onto a single variable loadcopy
    811954        int* indices=xNew<int>(nel);
    812955        for(int i=0;i<nel;i++)indices[i]=i;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r26676 r26800  
    488488                /*Nothing to add*/
    489489        }
     490        else if(hydrology_model==HydrologyTwsEnum){
     491                /*Nothing to add*/
     492        }
    490493        else{
    491494                _error_("Hydrology model "<<EnumToStringx(hydrology_model)<<" not supported yet");

  • issm/trunk-jpl/src/c/modules/modules.h

    r26526 r26800  
    3434#include "./GiaDeflectionCorex/GiaDeflectionCorex.h"
    3535#include "./FourierLoveCorex/FourierLoveCorex.h"
     36#include "./HypergeometricFunctionx/HypergeometricFunctionx.h"
    3637#include "./SetControlInputsFromVectorx/SetControlInputsFromVectorx.h"
    3738#include "./SetActiveNodesLSMx/SetActiveNodesLSMx.h"

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

    r26789 r26800  
    264264syn keyword cConstant LockFileNameEnum
    265265syn keyword cConstant LoveAllowLayerDeletionEnum
     266syn keyword cConstant LoveChandlerWobbleEnum
    266267syn keyword cConstant LoveCoreMantleBoundaryEnum
    267268syn keyword cConstant LoveEarthMassEnum
    … …  
    284285syn keyword cConstant LoveStartingLayerEnum
    285286syn keyword cConstant LoveUnderflowTolEnum
     287syn keyword cConstant LovePostWidderThresholdEnum
    286288syn keyword cConstant MassFluxSegmentsEnum
    287289syn keyword cConstant MassFluxSegmentsPresentEnum
    … …  
    368370syn keyword cConstant SolidearthSettingsAbstolEnum
    369371syn keyword cConstant SolidearthSettingsCrossSectionShapeEnum
    370 syn keyword cConstant RotationalAngularVelocityEnum
    371372syn keyword cConstant SolidearthSettingsElasticEnum
    372373syn keyword cConstant SolidearthSettingsViscousEnum
    … …  
    375376syn keyword cConstant SealevelchangeViscousTimesEnum
    376377syn keyword cConstant SealevelchangeViscousIndexEnum
     378syn keyword cConstant SealevelchangeViscousPolarMotionEnum
     379syn keyword cConstant SealevelchangeRunCountEnum
     380syn keyword cConstant SealevelchangeTransitionsEnum
     381syn keyword cConstant SealevelchangeRequestedOutputsEnum
     382syn keyword cConstant RotationalAngularVelocityEnum
    377383syn keyword cConstant RotationalEquatorialMoiEnum
     384syn keyword cConstant RotationalPolarMoiEnum
     385syn keyword cConstant LovePolarMotionTransferFunctionColinearEnum
     386syn keyword cConstant LovePolarMotionTransferFunctionOrthogonalEnum
    378387syn keyword cConstant TidalLoveHEnum
    379388syn keyword cConstant TidalLoveKEnum
    … …  
    387396syn keyword cConstant SealevelchangeGSelfAttractionEnum
    388397syn keyword cConstant SealevelchangeGViscoElasticEnum
     398syn keyword cConstant SealevelchangeUViscoElasticEnum
     399syn keyword cConstant SealevelchangeHViscoElasticEnum
     400syn keyword cConstant SealevelchangePolarMotionTransferFunctionColinearEnum
     401syn keyword cConstant SealevelchangePolarMotionTransferFunctionOrthogonalEnum
     402syn keyword cConstant SealevelchangePolarMotionTransferFunctionZEnum
     403syn keyword cConstant SealevelchangeTidalK2Enum
     404syn keyword cConstant SealevelchangeTidalH2Enum
     405syn keyword cConstant SealevelchangeTidalL2Enum
    389406syn keyword cConstant SolidearthSettingsSealevelLoadingEnum
    390407syn keyword cConstant SolidearthSettingsGRDEnum
    391408syn keyword cConstant SolidearthSettingsRunFrequencyEnum
    392409syn keyword cConstant SolidearthSettingsTimeAccEnum
    393 syn keyword cConstant SealevelchangeHViscoElasticEnum
    394410syn keyword cConstant SolidearthSettingsHorizEnum
    395411syn keyword cConstant SolidearthSettingsMaxiterEnum
    … …  
    407423syn keyword cConstant RotationalPolarMoiEnum
    408424syn keyword cConstant SolidearthSettingsReltolEnum
    409 syn keyword cConstant SealevelchangeRequestedOutputsEnum
    410425syn keyword cConstant SolidearthSettingsSelfAttractionEnum
    411426syn keyword cConstant SolidearthSettingsRotationEnum
    412427syn keyword cConstant SolidearthSettingsMaxSHCoeffEnum
    413 syn keyword cConstant SealevelchangeRunCountEnum
    414 syn keyword cConstant SealevelchangeTransitionsEnum
    415 syn keyword cConstant SealevelchangeUViscoElasticEnum
    416428syn keyword cConstant SettingsIoGatherEnum
    417429syn keyword cConstant SettingsNumResultsOnNodesEnum
    … …  
    847859syn keyword cConstant BslcRateEnum
    848860syn keyword cConstant GmtslcEnum
    849 syn keyword cConstant SealevelGrotm1Enum
    850 syn keyword cConstant SealevelGrotm2Enum
    851 syn keyword cConstant SealevelGrotm3Enum
    852 syn keyword cConstant SealevelGUrotm1Enum
    853 syn keyword cConstant SealevelGUrotm2Enum
    854 syn keyword cConstant SealevelGUrotm3Enum
    855 syn keyword cConstant SealevelGNrotm1Enum
    856 syn keyword cConstant SealevelGNrotm2Enum
    857 syn keyword cConstant SealevelGNrotm3Enum
    858 syn keyword cConstant SealevelGErotm1Enum
    859 syn keyword cConstant SealevelGErotm2Enum
    860 syn keyword cConstant SealevelGErotm3Enum
    861861syn keyword cConstant SealevelRSLBarystaticEnum
    862862syn keyword cConstant SealevelRSLRateEnum
    … …  
    870870syn keyword cConstant SealevelchangeGEEnum
    871871syn keyword cConstant SealevelchangeGNEnum
     872syn keyword cConstant SealevelchangeGrotEnum
     873syn keyword cConstant SealevelchangeGUrotEnum
     874syn keyword cConstant SealevelchangeGNrotEnum
     875syn keyword cConstant SealevelchangeGErotEnum
    872876syn keyword cConstant SealevelchangeGsubelOceanEnum
    873877syn keyword cConstant SealevelchangeGUsubelOceanEnum
    … …  
    13671371syn keyword cConstant LoadsEnum
    13681372syn keyword cConstant LoveAnalysisEnum
    1369 syn keyword cConstant LoveHiEnum
    1370 syn keyword cConstant LoveHrEnum
     1373syn keyword cConstant LoveHfEnum
     1374syn keyword cConstant LoveHtEnum
    13711375syn keyword cConstant LoveKernelsImagEnum
    13721376syn keyword cConstant LoveKernelsRealEnum
    1373 syn keyword cConstant LoveKiEnum
    1374 syn keyword cConstant LoveKrEnum
    1375 syn keyword cConstant LoveLiEnum
    1376 syn keyword cConstant LoveLrEnum
     1377syn keyword cConstant LoveKfEnum
     1378syn keyword cConstant LoveKtEnum
     1379syn keyword cConstant LoveLfEnum
     1380syn keyword cConstant LoveLtEnum
     1381syn keyword cConstant LoveTidalHtEnum
     1382syn keyword cConstant LoveTidalKtEnum
     1383syn keyword cConstant LoveTidalLtEnum
     1384syn keyword cConstant LovePMTF1tEnum
     1385syn keyword cConstant LovePMTF2tEnum
    13771386syn keyword cConstant LoveSolutionEnum
    13781387syn keyword cConstant MINIEnum
    … …  
    14861495syn keyword cConstant SealevelAbsoluteEnum
    14871496syn keyword cConstant SealevelEmotionEnum
    1488 syn keyword cConstant SealevelInertiaTensorXZEnum
    1489 syn keyword cConstant SealevelInertiaTensorYZEnum
    1490 syn keyword cConstant SealevelInertiaTensorZZEnum
     1497syn keyword cConstant SealevelchangePolarMotionXEnum
     1498syn keyword cConstant SealevelchangePolarMotionYEnum
     1499syn keyword cConstant SealevelchangePolarMotionZEnum
    14911500syn keyword cConstant SealevelchangePolarMotionEnum
    14921501syn keyword cConstant SealevelNmotionEnum

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

    r26789 r26800  
    258258        LockFileNameEnum,
    259259        LoveAllowLayerDeletionEnum,
     260        LoveChandlerWobbleEnum,
    260261        LoveCoreMantleBoundaryEnum,
    261262        LoveEarthMassEnum,
    … …  
    278279        LoveStartingLayerEnum,
    279280        LoveUnderflowTolEnum,
     281        LovePostWidderThresholdEnum,
    280282        MassFluxSegmentsEnum,
    281283        MassFluxSegmentsPresentEnum,
    … …  
    362364        SolidearthSettingsAbstolEnum,
    363365        SolidearthSettingsCrossSectionShapeEnum,
    364         RotationalAngularVelocityEnum,
    365366        SolidearthSettingsElasticEnum,
    366367        SolidearthSettingsViscousEnum,
    … …  
    369370        SealevelchangeViscousTimesEnum,
    370371        SealevelchangeViscousIndexEnum,
     372        SealevelchangeViscousPolarMotionEnum,
     373        SealevelchangeRunCountEnum,
     374        SealevelchangeTransitionsEnum,
     375        SealevelchangeRequestedOutputsEnum,
     376        RotationalAngularVelocityEnum,
    371377        RotationalEquatorialMoiEnum,
     378        RotationalPolarMoiEnum,
     379        LovePolarMotionTransferFunctionColinearEnum,
     380        LovePolarMotionTransferFunctionOrthogonalEnum,
    372381        TidalLoveHEnum,
    373382        TidalLoveKEnum,
    … …  
    381390        SealevelchangeGSelfAttractionEnum,
    382391        SealevelchangeGViscoElasticEnum,
     392        SealevelchangeUViscoElasticEnum,
     393        SealevelchangeHViscoElasticEnum,
     394        SealevelchangePolarMotionTransferFunctionColinearEnum,
     395        SealevelchangePolarMotionTransferFunctionOrthogonalEnum,
     396        SealevelchangePolarMotionTransferFunctionZEnum,
     397        SealevelchangeTidalK2Enum,
     398        SealevelchangeTidalH2Enum,
     399        SealevelchangeTidalL2Enum,
    383400        SolidearthSettingsSealevelLoadingEnum,
    384401        SolidearthSettingsGRDEnum,
    385402        SolidearthSettingsRunFrequencyEnum,
    386403        SolidearthSettingsTimeAccEnum,
    387         SealevelchangeHViscoElasticEnum,
    388404        SolidearthSettingsHorizEnum,
    389405        SolidearthSettingsMaxiterEnum,
    … …  
    399415        StochasticForcingNumFieldsEnum,
    400416        StochasticForcingRandomflagEnum,
    401         RotationalPolarMoiEnum,
    402417        SolidearthSettingsReltolEnum,
    403         SealevelchangeRequestedOutputsEnum,
    404418        SolidearthSettingsSelfAttractionEnum,
    405419        SolidearthSettingsRotationEnum,
    406420        SolidearthSettingsMaxSHCoeffEnum,
    407         SealevelchangeRunCountEnum,
    408         SealevelchangeTransitionsEnum,
    409         SealevelchangeUViscoElasticEnum,
    410421        SettingsIoGatherEnum,
    411422        SettingsNumResultsOnNodesEnum,
    … …  
    843854        BslcRateEnum,
    844855        GmtslcEnum,
    845         SealevelGrotm1Enum,
    846         SealevelGrotm2Enum,
    847         SealevelGrotm3Enum,
    848         SealevelGUrotm1Enum,
    849         SealevelGUrotm2Enum,
    850         SealevelGUrotm3Enum,
    851         SealevelGNrotm1Enum,
    852         SealevelGNrotm2Enum,
    853         SealevelGNrotm3Enum,
    854         SealevelGErotm1Enum,
    855         SealevelGErotm2Enum,
    856         SealevelGErotm3Enum,
    857856        SealevelRSLBarystaticEnum,
    858857        SealevelRSLRateEnum,
    … …  
    866865        SealevelchangeGEEnum,
    867866        SealevelchangeGNEnum,
     867        SealevelchangeGrotEnum,
     868        SealevelchangeGUrotEnum,
     869        SealevelchangeGNrotEnum,
     870        SealevelchangeGErotEnum,
    868871        SealevelchangeGsubelOceanEnum,
    869872        SealevelchangeGUsubelOceanEnum,
    … …  
    13661369        LoadsEnum,
    13671370        LoveAnalysisEnum,
    1368         LoveHiEnum,
    1369         LoveHrEnum,
     1371        LoveHfEnum,
     1372        LoveHtEnum,
    13701373        LoveKernelsImagEnum,
    13711374        LoveKernelsRealEnum,
    1372         LoveKiEnum,
    1373         LoveKrEnum,
    1374         LoveLiEnum,
    1375         LoveLrEnum,
     1375        LoveKfEnum,
     1376        LoveKtEnum,
     1377        LoveLfEnum,
     1378        LoveLtEnum,
     1379        LoveTidalHtEnum,
     1380        LoveTidalKtEnum,
     1381        LoveTidalLtEnum,
     1382        LovePMTF1tEnum,
     1383        LovePMTF2tEnum,
    13761384        LoveSolutionEnum,
    13771385        MINIEnum,
    … …  
    14851493        SealevelAbsoluteEnum,
    14861494        SealevelEmotionEnum,
    1487         SealevelInertiaTensorXZEnum,
    1488         SealevelInertiaTensorYZEnum,
    1489         SealevelInertiaTensorZZEnum,
     1495        SealevelchangePolarMotionXEnum,
     1496        SealevelchangePolarMotionYEnum,
     1497        SealevelchangePolarMotionZEnum,
    14901498        SealevelchangePolarMotionEnum,
    14911499        SealevelNmotionEnum,

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

    r26791 r26800  
    266266                case LockFileNameEnum : return "LockFileName";
    267267                case LoveAllowLayerDeletionEnum : return "LoveAllowLayerDeletion";
     268                case LoveChandlerWobbleEnum : return "LoveChandlerWobble";
    268269                case LoveCoreMantleBoundaryEnum : return "LoveCoreMantleBoundary";
    269270                case LoveEarthMassEnum : return "LoveEarthMass";
    … …  
    286287                case LoveStartingLayerEnum : return "LoveStartingLayer";
    287288                case LoveUnderflowTolEnum : return "LoveUnderflowTol";
     289                case LovePostWidderThresholdEnum : return "LovePostWidderThreshold";
    288290                case MassFluxSegmentsEnum : return "MassFluxSegments";
    289291                case MassFluxSegmentsPresentEnum : return "MassFluxSegmentsPresent";
    … …  
    370372                case SolidearthSettingsAbstolEnum : return "SolidearthSettingsAbstol";
    371373                case SolidearthSettingsCrossSectionShapeEnum : return "SolidearthSettingsCrossSectionShape";
    372                 case RotationalAngularVelocityEnum : return "RotationalAngularVelocity";
    373374                case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
    374375                case SolidearthSettingsViscousEnum : return "SolidearthSettingsViscous";
    … …  
    377378                case SealevelchangeViscousTimesEnum : return "SealevelchangeViscousTimes";
    378379                case SealevelchangeViscousIndexEnum : return "SealevelchangeViscousIndex";
     380                case SealevelchangeViscousPolarMotionEnum : return "SealevelchangeViscousPolarMotion";
     381                case SealevelchangeRunCountEnum : return "SealevelchangeRunCount";
     382                case SealevelchangeTransitionsEnum : return "SealevelchangeTransitions";
     383                case SealevelchangeRequestedOutputsEnum : return "SealevelchangeRequestedOutputs";
     384                case RotationalAngularVelocityEnum : return "RotationalAngularVelocity";
    379385                case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi";
     386                case RotationalPolarMoiEnum : return "RotationalPolarMoi";
     387                case LovePolarMotionTransferFunctionColinearEnum : return "LovePolarMotionTransferFunctionColinear";
     388                case LovePolarMotionTransferFunctionOrthogonalEnum : return "LovePolarMotionTransferFunctionOrthogonal";
    380389                case TidalLoveHEnum : return "TidalLoveH";
    381390                case TidalLoveKEnum : return "TidalLoveK";
    … …  
    389398                case SealevelchangeGSelfAttractionEnum : return "SealevelchangeGSelfAttraction";
    390399                case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
     400                case SealevelchangeUViscoElasticEnum : return "SealevelchangeUViscoElastic";
     401                case SealevelchangeHViscoElasticEnum : return "SealevelchangeHViscoElastic";
     402                case SealevelchangePolarMotionTransferFunctionColinearEnum : return "SealevelchangePolarMotionTransferFunctionColinear";
     403                case SealevelchangePolarMotionTransferFunctionOrthogonalEnum : return "SealevelchangePolarMotionTransferFunctionOrthogonal";
     404                case SealevelchangePolarMotionTransferFunctionZEnum : return "SealevelchangePolarMotionTransferFunctionZ";
     405                case SealevelchangeTidalK2Enum : return "SealevelchangeTidalK2";
     406                case SealevelchangeTidalH2Enum : return "SealevelchangeTidalH2";
     407                case SealevelchangeTidalL2Enum : return "SealevelchangeTidalL2";
    391408                case SolidearthSettingsSealevelLoadingEnum : return "SolidearthSettingsSealevelLoading";
    392409                case SolidearthSettingsGRDEnum : return "SolidearthSettingsGRD";
    393410                case SolidearthSettingsRunFrequencyEnum : return "SolidearthSettingsRunFrequency";
    394411                case SolidearthSettingsTimeAccEnum : return "SolidearthSettingsTimeAcc";
    395                 case SealevelchangeHViscoElasticEnum : return "SealevelchangeHViscoElastic";
    396412                case SolidearthSettingsHorizEnum : return "SolidearthSettingsHoriz";
    397413                case SolidearthSettingsMaxiterEnum : return "SolidearthSettingsMaxiter";
    … …  
    407423                case StochasticForcingNumFieldsEnum : return "StochasticForcingNumFields";
    408424                case StochasticForcingRandomflagEnum : return "StochasticForcingRandomflag";
    409                 case RotationalPolarMoiEnum : return "RotationalPolarMoi";
    410425                case SolidearthSettingsReltolEnum : return "SolidearthSettingsReltol";
    411                 case SealevelchangeRequestedOutputsEnum : return "SealevelchangeRequestedOutputs";
    412426                case SolidearthSettingsSelfAttractionEnum : return "SolidearthSettingsSelfAttraction";
    413427                case SolidearthSettingsRotationEnum : return "SolidearthSettingsRotation";
    414428                case SolidearthSettingsMaxSHCoeffEnum : return "SolidearthSettingsMaxSHCoeff";
    415                 case SealevelchangeRunCountEnum : return "SealevelchangeRunCount";
    416                 case SealevelchangeTransitionsEnum : return "SealevelchangeTransitions";
    417                 case SealevelchangeUViscoElasticEnum : return "SealevelchangeUViscoElastic";
    418429                case SettingsIoGatherEnum : return "SettingsIoGather";
    419430                case SettingsNumResultsOnNodesEnum : return "SettingsNumResultsOnNodes";
    … …  
    849860                case BslcRateEnum : return "BslcRate";
    850861                case GmtslcEnum : return "Gmtslc";
    851                 case SealevelGrotm1Enum : return "SealevelGrotm1";
    852                 case SealevelGrotm2Enum : return "SealevelGrotm2";
    853                 case SealevelGrotm3Enum : return "SealevelGrotm3";
    854                 case SealevelGUrotm1Enum : return "SealevelGUrotm1";
    855                 case SealevelGUrotm2Enum : return "SealevelGUrotm2";
    856                 case SealevelGUrotm3Enum : return "SealevelGUrotm3";
    857                 case SealevelGNrotm1Enum : return "SealevelGNrotm1";
    858                 case SealevelGNrotm2Enum : return "SealevelGNrotm2";
    859                 case SealevelGNrotm3Enum : return "SealevelGNrotm3";
    860                 case SealevelGErotm1Enum : return "SealevelGErotm1";
    861                 case SealevelGErotm2Enum : return "SealevelGErotm2";
    862                 case SealevelGErotm3Enum : return "SealevelGErotm3";
    863862                case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic";
    864863                case SealevelRSLRateEnum : return "SealevelRSLRate";
    … …  
    872871                case SealevelchangeGEEnum : return "SealevelchangeGE";
    873872                case SealevelchangeGNEnum : return "SealevelchangeGN";
     873                case SealevelchangeGrotEnum : return "SealevelchangeGrot";
     874                case SealevelchangeGUrotEnum : return "SealevelchangeGUrot";
     875                case SealevelchangeGNrotEnum : return "SealevelchangeGNrot";
     876                case SealevelchangeGErotEnum : return "SealevelchangeGErot";
    874877                case SealevelchangeGsubelOceanEnum : return "SealevelchangeGsubelOcean";
    875878                case SealevelchangeGUsubelOceanEnum : return "SealevelchangeGUsubelOcean";
    … …  
    13691372                case LoadsEnum : return "Loads";
    13701373                case LoveAnalysisEnum : return "LoveAnalysis";
    1371                 case LoveHiEnum : return "LoveHi";
    1372                 case LoveHrEnum : return "LoveHr";
     1374                case LoveHfEnum : return "LoveHf";
     1375                case LoveHtEnum : return "LoveHt";
    13731376                case LoveKernelsImagEnum : return "LoveKernelsImag";
    13741377                case LoveKernelsRealEnum : return "LoveKernelsReal";
    1375                 case LoveKiEnum : return "LoveKi";
    1376                 case LoveKrEnum : return "LoveKr";
    1377                 case LoveLiEnum : return "LoveLi";
    1378                 case LoveLrEnum : return "LoveLr";
     1378                case LoveKfEnum : return "LoveKf";
     1379                case LoveKtEnum : return "LoveKt";
     1380                case LoveLfEnum : return "LoveLf";
     1381                case LoveLtEnum : return "LoveLt";
     1382                case LoveTidalHtEnum : return "LoveTidalHt";
     1383                case LoveTidalKtEnum : return "LoveTidalKt";
     1384                case LoveTidalLtEnum : return "LoveTidalLt";
     1385                case LovePMTF1tEnum : return "LovePMTF1t";
     1386                case LovePMTF2tEnum : return "LovePMTF2t";
    13791387                case LoveSolutionEnum : return "LoveSolution";
    13801388                case MINIEnum : return "MINI";
    … …  
    14881496                case SealevelAbsoluteEnum : return "SealevelAbsolute";
    14891497                case SealevelEmotionEnum : return "SealevelEmotion";
    1490                 case SealevelInertiaTensorXZEnum : return "SealevelInertiaTensorXZ";
    1491                 case SealevelInertiaTensorYZEnum : return "SealevelInertiaTensorYZ";
    1492                 case SealevelInertiaTensorZZEnum : return "SealevelInertiaTensorZZ";
     1498                case SealevelchangePolarMotionXEnum : return "SealevelchangePolarMotionX";
     1499                case SealevelchangePolarMotionYEnum : return "SealevelchangePolarMotionY";
     1500                case SealevelchangePolarMotionZEnum : return "SealevelchangePolarMotionZ";
    14931501                case SealevelchangePolarMotionEnum : return "SealevelchangePolarMotion";
    14941502                case SealevelNmotionEnum : return "SealevelNmotion";

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

    r26789 r26800  
    257257syn keyword juliaConstC LockFileNameEnum
    258258syn keyword juliaConstC LoveAllowLayerDeletionEnum
     259syn keyword juliaConstC LoveChandlerWobbleEnum
    259260syn keyword juliaConstC LoveCoreMantleBoundaryEnum
    260261syn keyword juliaConstC LoveEarthMassEnum
    … …  
    277278syn keyword juliaConstC LoveStartingLayerEnum
    278279syn keyword juliaConstC LoveUnderflowTolEnum
     280syn keyword juliaConstC LovePostWidderThresholdEnum
    279281syn keyword juliaConstC MassFluxSegmentsEnum
    280282syn keyword juliaConstC MassFluxSegmentsPresentEnum
    … …  
    361363syn keyword juliaConstC SolidearthSettingsAbstolEnum
    362364syn keyword juliaConstC SolidearthSettingsCrossSectionShapeEnum
    363 syn keyword juliaConstC RotationalAngularVelocityEnum
    364365syn keyword juliaConstC SolidearthSettingsElasticEnum
    365366syn keyword juliaConstC SolidearthSettingsViscousEnum
    … …  
    368369syn keyword juliaConstC SealevelchangeViscousTimesEnum
    369370syn keyword juliaConstC SealevelchangeViscousIndexEnum
     371syn keyword juliaConstC SealevelchangeViscousPolarMotionEnum
     372syn keyword juliaConstC SealevelchangeRunCountEnum
     373syn keyword juliaConstC SealevelchangeTransitionsEnum
     374syn keyword juliaConstC SealevelchangeRequestedOutputsEnum
     375syn keyword juliaConstC RotationalAngularVelocityEnum
    370376syn keyword juliaConstC RotationalEquatorialMoiEnum
     377syn keyword juliaConstC RotationalPolarMoiEnum
     378syn keyword juliaConstC LovePolarMotionTransferFunctionColinearEnum
     379syn keyword juliaConstC LovePolarMotionTransferFunctionOrthogonalEnum
    371380syn keyword juliaConstC TidalLoveHEnum
    372381syn keyword juliaConstC TidalLoveKEnum
    … …  
    380389syn keyword juliaConstC SealevelchangeGSelfAttractionEnum
    381390syn keyword juliaConstC SealevelchangeGViscoElasticEnum
     391syn keyword juliaConstC SealevelchangeUViscoElasticEnum
     392syn keyword juliaConstC SealevelchangeHViscoElasticEnum
     393syn keyword juliaConstC SealevelchangePolarMotionTransferFunctionColinearEnum
     394syn keyword juliaConstC SealevelchangePolarMotionTransferFunctionOrthogonalEnum
     395syn keyword juliaConstC SealevelchangePolarMotionTransferFunctionZEnum
     396syn keyword juliaConstC SealevelchangeTidalK2Enum
     397syn keyword juliaConstC SealevelchangeTidalH2Enum
     398syn keyword juliaConstC SealevelchangeTidalL2Enum
    382399syn keyword juliaConstC SolidearthSettingsSealevelLoadingEnum
    383400syn keyword juliaConstC SolidearthSettingsGRDEnum
    384401syn keyword juliaConstC SolidearthSettingsRunFrequencyEnum
    385402syn keyword juliaConstC SolidearthSettingsTimeAccEnum
    386 syn keyword juliaConstC SealevelchangeHViscoElasticEnum
    387403syn keyword juliaConstC SolidearthSettingsHorizEnum
    388404syn keyword juliaConstC SolidearthSettingsMaxiterEnum
    … …  
    398414syn keyword juliaConstC StochasticForcingNumFieldsEnum
    399415syn keyword juliaConstC StochasticForcingRandomflagEnum
    400 syn keyword juliaConstC RotationalPolarMoiEnum
    401416syn keyword juliaConstC SolidearthSettingsReltolEnum
    402 syn keyword juliaConstC SealevelchangeRequestedOutputsEnum
    403417syn keyword juliaConstC SolidearthSettingsSelfAttractionEnum
    404418syn keyword juliaConstC SolidearthSettingsRotationEnum
    405419syn keyword juliaConstC SolidearthSettingsMaxSHCoeffEnum
    406 syn keyword juliaConstC SealevelchangeRunCountEnum
    407 syn keyword juliaConstC SealevelchangeTransitionsEnum
    408 syn keyword juliaConstC SealevelchangeUViscoElasticEnum
    409420syn keyword juliaConstC SettingsIoGatherEnum
    410421syn keyword juliaConstC SettingsNumResultsOnNodesEnum
    … …  
    840851syn keyword juliaConstC BslcRateEnum
    841852syn keyword juliaConstC GmtslcEnum
    842 syn keyword juliaConstC SealevelGrotm1Enum
    843 syn keyword juliaConstC SealevelGrotm2Enum
    844 syn keyword juliaConstC SealevelGrotm3Enum
    845 syn keyword juliaConstC SealevelGUrotm1Enum
    846 syn keyword juliaConstC SealevelGUrotm2Enum
    847 syn keyword juliaConstC SealevelGUrotm3Enum
    848 syn keyword juliaConstC SealevelGNrotm1Enum
    849 syn keyword juliaConstC SealevelGNrotm2Enum
    850 syn keyword juliaConstC SealevelGNrotm3Enum
    851 syn keyword juliaConstC SealevelGErotm1Enum
    852 syn keyword juliaConstC SealevelGErotm2Enum
    853 syn keyword juliaConstC SealevelGErotm3Enum
    854853syn keyword juliaConstC SealevelRSLBarystaticEnum
    855854syn keyword juliaConstC SealevelRSLRateEnum
    … …  
    863862syn keyword juliaConstC SealevelchangeGEEnum
    864863syn keyword juliaConstC SealevelchangeGNEnum
     864syn keyword juliaConstC SealevelchangeGrotEnum
     865syn keyword juliaConstC SealevelchangeGUrotEnum
     866syn keyword juliaConstC SealevelchangeGNrotEnum
     867syn keyword juliaConstC SealevelchangeGErotEnum
    865868syn keyword juliaConstC SealevelchangeGsubelOceanEnum
    866869syn keyword juliaConstC SealevelchangeGUsubelOceanEnum
    … …  
    13601363syn keyword juliaConstC LoadsEnum
    13611364syn keyword juliaConstC LoveAnalysisEnum
    1362 syn keyword juliaConstC LoveHiEnum
    1363 syn keyword juliaConstC LoveHrEnum
     1365syn keyword juliaConstC LoveHfEnum
     1366syn keyword juliaConstC LoveHtEnum
    13641367syn keyword juliaConstC LoveKernelsImagEnum
    13651368syn keyword juliaConstC LoveKernelsRealEnum
    1366 syn keyword juliaConstC LoveKiEnum
    1367 syn keyword juliaConstC LoveKrEnum
    1368 syn keyword juliaConstC LoveLiEnum
    1369 syn keyword juliaConstC LoveLrEnum
     1369syn keyword juliaConstC LoveKfEnum
     1370syn keyword juliaConstC LoveKtEnum
     1371syn keyword juliaConstC LoveLfEnum
     1372syn keyword juliaConstC LoveLtEnum
     1373syn keyword juliaConstC LoveTidalHtEnum
     1374syn keyword juliaConstC LoveTidalKtEnum
     1375syn keyword juliaConstC LoveTidalLtEnum
     1376syn keyword juliaConstC LovePMTF1tEnum
     1377syn keyword juliaConstC LovePMTF2tEnum
    13701378syn keyword juliaConstC LoveSolutionEnum
    13711379syn keyword juliaConstC MINIEnum
    … …  
    14791487syn keyword juliaConstC SealevelAbsoluteEnum
    14801488syn keyword juliaConstC SealevelEmotionEnum
    1481 syn keyword juliaConstC SealevelInertiaTensorXZEnum
    1482 syn keyword juliaConstC SealevelInertiaTensorYZEnum
    1483 syn keyword juliaConstC SealevelInertiaTensorZZEnum
     1489syn keyword juliaConstC SealevelchangePolarMotionXEnum
     1490syn keyword juliaConstC SealevelchangePolarMotionYEnum
     1491syn keyword juliaConstC SealevelchangePolarMotionZEnum
    14841492syn keyword juliaConstC SealevelchangePolarMotionEnum
    14851493syn keyword juliaConstC SealevelNmotionEnum

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

    r26791 r26800  
    272272              else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
    273273              else if (strcmp(name,"LoveAllowLayerDeletion")==0) return LoveAllowLayerDeletionEnum;
     274              else if (strcmp(name,"LoveChandlerWobble")==0) return LoveChandlerWobbleEnum;
    274275              else if (strcmp(name,"LoveCoreMantleBoundary")==0) return LoveCoreMantleBoundaryEnum;
    275276              else if (strcmp(name,"LoveEarthMass")==0) return LoveEarthMassEnum;
    … …  
    292293              else if (strcmp(name,"LoveStartingLayer")==0) return LoveStartingLayerEnum;
    293294              else if (strcmp(name,"LoveUnderflowTol")==0) return LoveUnderflowTolEnum;
     295              else if (strcmp(name,"LovePostWidderThreshold")==0) return LovePostWidderThresholdEnum;
    294296              else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
    295297              else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
    … …  
    376378              else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum;
    377379              else if (strcmp(name,"SolidearthSettingsCrossSectionShape")==0) return SolidearthSettingsCrossSectionShapeEnum;
    378               else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum;
    379380              else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
    380381              else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
    … …  
    382383              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
    383384              else if (strcmp(name,"SealevelchangeViscousTimes")==0) return SealevelchangeViscousTimesEnum;
    384               else if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
     388              if (strcmp(name,"SealevelchangeViscousIndex")==0) return SealevelchangeViscousIndexEnum;
     389              else if (strcmp(name,"SealevelchangeViscousPolarMotion")==0) return SealevelchangeViscousPolarMotionEnum;
     390              else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
     391              else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
     392              else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
     393              else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum;
     394              else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
     395              else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum;
     396              else if (strcmp(name,"LovePolarMotionTransferFunctionColinear")==0) return LovePolarMotionTransferFunctionColinearEnum;
     397              else if (strcmp(name,"LovePolarMotionTransferFunctionOrthogonal")==0) return LovePolarMotionTransferFunctionOrthogonalEnum;
    389398              else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum;
    390399              else if (strcmp(name,"TidalLoveK")==0) return TidalLoveKEnum;
    … …  
    398407              else if (strcmp(name,"SealevelchangeGSelfAttraction")==0) return SealevelchangeGSelfAttractionEnum;
    399408              else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
     409              else if (strcmp(name,"SealevelchangeUViscoElastic")==0) return SealevelchangeUViscoElasticEnum;
     410              else if (strcmp(name,"SealevelchangeHViscoElastic")==0) return SealevelchangeHViscoElasticEnum;
     411              else if (strcmp(name,"SealevelchangePolarMotionTransferFunctionColinear")==0) return SealevelchangePolarMotionTransferFunctionColinearEnum;
     412              else if (strcmp(name,"SealevelchangePolarMotionTransferFunctionOrthogonal")==0) return SealevelchangePolarMotionTransferFunctionOrthogonalEnum;
     413              else if (strcmp(name,"SealevelchangePolarMotionTransferFunctionZ")==0) return SealevelchangePolarMotionTransferFunctionZEnum;
     414              else if (strcmp(name,"SealevelchangeTidalK2")==0) return SealevelchangeTidalK2Enum;
     415              else if (strcmp(name,"SealevelchangeTidalH2")==0) return SealevelchangeTidalH2Enum;
     416              else if (strcmp(name,"SealevelchangeTidalL2")==0) return SealevelchangeTidalL2Enum;
    400417              else if (strcmp(name,"SolidearthSettingsSealevelLoading")==0) return SolidearthSettingsSealevelLoadingEnum;
    401418              else if (strcmp(name,"SolidearthSettingsGRD")==0) return SolidearthSettingsGRDEnum;
    402419              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    403420              else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
    404               else if (strcmp(name,"SealevelchangeHViscoElastic")==0) return SealevelchangeHViscoElasticEnum;
    405421              else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
    406422              else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
    … …  
    416432              else if (strcmp(name,"StochasticForcingNumFields")==0) return StochasticForcingNumFieldsEnum;
    417433              else if (strcmp(name,"StochasticForcingRandomflag")==0) return StochasticForcingRandomflagEnum;
    418               else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum;
    419434              else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    420               else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
    421435              else if (strcmp(name,"SolidearthSettingsSelfAttraction")==0) return SolidearthSettingsSelfAttractionEnum;
    422436              else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
    423437              else if (strcmp(name,"SolidearthSettingsMaxSHCoeff")==0) return SolidearthSettingsMaxSHCoeffEnum;
    424               else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
    425               else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
    426               else if (strcmp(name,"SealevelchangeUViscoElastic")==0) return SealevelchangeUViscoElasticEnum;
    427438              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
    428439              else if (strcmp(name,"SettingsNumResultsOnNodes")==0) return SettingsNumResultsOnNodesEnum;
    … …  
    495506              else if (strcmp(name,"SmbTdiff")==0) return SmbTdiffEnum;
    496507              else if (strcmp(name,"SmbThermoDeltaTScaling")==0) return SmbThermoDeltaTScalingEnum;
    497               else if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum;
    498512              else if (strcmp(name,"SmbPrecipitationsReconstructedYears")==0) return SmbPrecipitationsReconstructedYearsEnum;
    499513              else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
    … …  
    506520              else if (strcmp(name,"Steps")==0) return StepsEnum;
    507521              else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
     522              else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum;
    512523              else if (strcmp(name,"StressbalanceIsnewton")==0) return StressbalanceIsnewtonEnum;
    513524              else if (strcmp(name,"StressbalanceMaxiter")==0) return StressbalanceMaxiterEnum;
    … …  
    618629              else if (strcmp(name,"Base")==0) return BaseEnum;
    619630              else if (strcmp(name,"BaseOld")==0) return BaseOldEnum;
    620               else if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"BaseSlopeX")==0) return BaseSlopeXEnum;
    621635              else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum;
    622636              else if (strcmp(name,"BaselineBasalforcingsFloatingiceMeltingRate")==0) return BaselineBasalforcingsFloatingiceMeltingRateEnum;
    … …  
    629643              else if (strcmp(name,"BedEastGRD")==0) return BedEastGRDEnum;
    630644              else if (strcmp(name,"BedNorth")==0) return BedNorthEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum;
     645              else if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum;
    635646              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    636647              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    … …  
    741752              else if (strcmp(name,"HydrologydcEplThicknessSubstep")==0) return HydrologydcEplThicknessSubstepEnum;
    742753              else if (strcmp(name,"HydrologydcEplThicknessTransient")==0) return HydrologydcEplThicknessTransientEnum;
    743               else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum;
    744758              else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
    745759              else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
    … …  
    752766              else if (strcmp(name,"HydrologyGapHeightXX")==0) return HydrologyGapHeightXXEnum;
    753767              else if (strcmp(name,"HydrologyGapHeightY")==0) return HydrologyGapHeightYEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
     768              else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
    758769              else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
    759770              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    … …  
    864875              else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
    865876              else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
    866               else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
    867               else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
    868               else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
    869               else if (strcmp(name,"SealevelGrotm1")==0) return SealevelGrotm1Enum;
    870               else if (strcmp(name,"SealevelGrotm2")==0) return SealevelGrotm2Enum;
    871               else if (strcmp(name,"SealevelGrotm3")==0) return SealevelGrotm3Enum;
    872               else if (strcmp(name,"SealevelGUrotm1")==0) return SealevelGUrotm1Enum;
    873               else if (strcmp(name,"SealevelGUrotm2")==0) return SealevelGUrotm2Enum;
    874               else if (strcmp(name,"SealevelGUrotm3")==0) return SealevelGUrotm3Enum;
    875               else if (strcmp(name,"SealevelGNrotm1")==0) return SealevelGNrotm1Enum;
    876               else if (strcmp(name,"SealevelGNrotm2")==0) return SealevelGNrotm2Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"SealevelGNrotm3")==0) return SealevelGNrotm3Enum;
    881               else if (strcmp(name,"SealevelGErotm1")==0) return SealevelGErotm1Enum;
    882               else if (strcmp(name,"SealevelGErotm2")==0) return SealevelGErotm2Enum;
    883               else if (strcmp(name,"SealevelGErotm3")==0) return SealevelGErotm3Enum;
     880              if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
     881              else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
     882              else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
    884883              else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
    885884              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    … …  
    893892              else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    894893              else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
     894              else if (strcmp(name,"SealevelchangeGrot")==0) return SealevelchangeGrotEnum;
     895              else if (strcmp(name,"SealevelchangeGUrot")==0) return SealevelchangeGUrotEnum;
     896              else if (strcmp(name,"SealevelchangeGNrot")==0) return SealevelchangeGNrotEnum;
     897              else if (strcmp(name,"SealevelchangeGErot")==0) return SealevelchangeGErotEnum;
    895898              else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
    896899              else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
    … …  
    995998              else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum;
    996999              else if (strcmp(name,"SmbRain")==0) return SmbRainEnum;
    997               else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
    998               else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
    999               else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
     1003              if (strcmp(name,"SmbRe")==0) return SmbReEnum;
     1004              else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
     1005              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
     1006              else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
    10041007              else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    10051008              else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
    … …  
    11181121              else if (strcmp(name,"OldAccumulatedDeltaTws")==0) return OldAccumulatedDeltaTwsEnum;
    11191122              else if (strcmp(name,"Outputdefinition1")==0) return Outputdefinition1Enum;
    1120               else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum;
    1121               else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum;
    1122               else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
     1126              if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum;
     1127              else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum;
     1128              else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
     1129              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    11271130              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    11281131              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    … …  
    12411244              else if (strcmp(name,"BalancevelocityAnalysis")==0) return BalancevelocityAnalysisEnum;
    12421245              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    1243               else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
    1244               else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
    1245               else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     1249              if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
     1250              else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
     1251              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
     1252              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    12501253              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    12511254              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    … …  
    13641367              else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum;
    13651368              else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum;
    1366               else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum;
    1367               else if (strcmp(name,"IceMass")==0) return IceMassEnum;
    1368               else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
     1372              if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum;
     1373              else if (strcmp(name,"IceMass")==0) return IceMassEnum;
     1374              else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
     1375              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    13731376              else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
    13741377              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    … …  
    14021405              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    14031406              else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
    1404               else if (strcmp(name,"LoveHi")==0) return LoveHiEnum;
    1405               else if (strcmp(name,"LoveHr")==0) return LoveHrEnum;
     1407              else if (strcmp(name,"LoveHf")==0) return LoveHfEnum;
     1408              else if (strcmp(name,"LoveHt")==0) return LoveHtEnum;
    14061409              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    14071410              else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    1408               else if (strcmp(name,"LoveKi")==0) return LoveKiEnum;
    1409               else if (strcmp(name,"LoveKr")==0) return LoveKrEnum;
    1410               else if (strcmp(name,"LoveLi")==0) return LoveLiEnum;
    1411               else if (strcmp(name,"LoveLr")==0) return LoveLrEnum;
     1411              else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
     1412              else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
     1413              else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
     1414              else if (strcmp(name,"LoveLt")==0) return LoveLtEnum;
     1415              else if (strcmp(name,"LoveTidalHt")==0) return LoveTidalHtEnum;
     1416              else if (strcmp(name,"LoveTidalKt")==0) return LoveTidalKtEnum;
     1417              else if (strcmp(name,"LoveTidalLt")==0) return LoveTidalLtEnum;
     1418              else if (strcmp(name,"LovePMTF1t")==0) return LovePMTF1tEnum;
     1419              else if (strcmp(name,"LovePMTF2t")==0) return LovePMTF2tEnum;
    14121420              else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    14131421              else if (strcmp(name,"MINI")==0) return MINIEnum;
    … …  
    14821490              else if (strcmp(name,"P2")==0) return P2Enum;
    14831491              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    1484               else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    14851496              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
    14861497              else if (strcmp(name,"P2xP4")==0) return P2xP4Enum;
    … …  
    14901501              else if (strcmp(name,"Penta")==0) return PentaEnum;
    14911502              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"Profiler")==0) return ProfilerEnum;
     1503              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    14961504              else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
    14971505              else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;
    … …  
    15241532              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    15251533              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    1526               else if (strcmp(name,"SealevelInertiaTensorXZ")==0) return SealevelInertiaTensorXZEnum;
    1527               else if (strcmp(name,"SealevelInertiaTensorYZ")==0) return SealevelInertiaTensorYZEnum;
    1528               else if (strcmp(name,"SealevelInertiaTensorZZ")==0) return SealevelInertiaTensorZZEnum;
     1534              else if (strcmp(name,"SealevelchangePolarMotionX")==0) return SealevelchangePolarMotionXEnum;
     1535              else if (strcmp(name,"SealevelchangePolarMotionY")==0) return SealevelchangePolarMotionYEnum;
     1536              else if (strcmp(name,"SealevelchangePolarMotionZ")==0) return SealevelchangePolarMotionZEnum;
    15291537              else if (strcmp(name,"SealevelchangePolarMotion")==0) return SealevelchangePolarMotionEnum;
    15301538              else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
    … …  
    16051613              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    16061614              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
    1607               else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
    16081619              else if (strcmp(name,"StrainRate")==0) return StrainRateEnum;
    16091620              else if (strcmp(name,"StressTensor")==0) return StressTensorEnum;
    … …  
    16131624              else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum;
    16141625              else if (strcmp(name,"MaximumNumberOfDefinitions")==0) return MaximumNumberOfDefinitionsEnum;
    1615          else stage=14;
     1626         else stage=15;
    16161627   }
    16171628        /*If we reach this point, the string provided has not been found*/

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

    r26301 r26800  
    1414                mu0                         = 0;
    1515                Gravitational_Constant      = 0;
     16                chandler_wobble             = 0;
    1617                allow_layer_deletion        = 0;
    1718                underflow_tol               = 0;
     19                pw_threshold                = 0;
    1820                integration_steps_per_layer = 0;
    1921                istemporal                  = 0;
    … …  
    5355                        self.mu0=1e11; % Pa
    5456                        self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2
     57                        self.chandler_wobble=0;
    5558                        self.allow_layer_deletion=1;
    5659                        self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer
     60                        self.pw_threshold=1e-3; %if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed
    5761                        self.integration_steps_per_layer=100;
    5862                        self.istemporal=0;
    … …  
    7478                        fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
    7579                        fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
    76                         fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');
     80                        fielddisplay(self,'chandler_wobble','includes the inertial terms for the chandler wobble in the rotational feedback love numbers, only for forcing_type=11 (default: 0) (/!\ 1 is untested yet)');
     81                        fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');                   
    7782                        fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
     83                        fielddisplay(self,'pw_threshold','if relative variation across frequencies is smaller than this ratio, the post-widder transform for time-dependent love numbers is bypassed (default (1e-3)');
    7884                        fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)');
    7985                        fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
    … …  
    98104                        md = checkfield(md,'fieldname','love.mu0','NaN',1,'Inf',1,'numel',1,'>',0);
    99105                        md = checkfield(md,'fieldname','love.Gravitational_Constant','NaN',1,'Inf',1,'numel',1,'>',0);
     106                        md = checkfield(md,'fieldname','love.chandler_wobble','values',[0 1]);
    100107                        md = checkfield(md,'fieldname','love.allow_layer_deletion','values',[0 1]);
    101108                        md = checkfield(md,'fieldname','love.underflow_tol','NaN',1,'Inf',1,'numel',1,'>',0);
     109                        md = checkfield(md,'fieldname','love.pw_threshold','NaN',1,'Inf',1,'numel',1,'>',0);
    102110                        md = checkfield(md,'fieldname','love.integration_steps_per_layer','NaN',1,'Inf',1,'numel',1,'>',0);
    103111                        md = checkfield(md,'fieldname','love.love_kernels','values',[0 1]);
    … …  
    112120                        if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9)
    113121                                error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.'])
     122                        end
     123
     124                        if md.love.chandler_wobble==1
     125                                disp('Warning, Chandler Wobble in Love number calculator has not been validated yet');
    114126                        end
    115127
    … …  
    137149                        WriteData(fid,prefix,'object',self,'fieldname','mu0','format','Double');
    138150                        WriteData(fid,prefix,'object',self,'fieldname','Gravitational_Constant','format','Double');
     151                        WriteData(fid,prefix,'object',self,'fieldname','chandler_wobble','format','Boolean');
    139152                        WriteData(fid,prefix,'object',self,'fieldname','allow_layer_deletion','format','Boolean');
    140153                        WriteData(fid,prefix,'object',self,'fieldname','underflow_tol','format','Double');
     154                        WriteData(fid,prefix,'object',self,'fieldname','pw_threshold','format','Double');
    141155                        WriteData(fid,prefix,'object',self,'fieldname','integration_steps_per_layer','format','Integer');
    142156                        WriteData(fid,prefix,'object',self,'fieldname','istemporal','format','Boolean');
    … …  
    164178                        for i=1:length(self.time)
    165179                                for j=1:2*self.n_temporal_iterations
    166                                         self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi;
     180                                        if self.time(i)==0
     181                                                self.frequencies((i-1)*2*self.n_temporal_iterations +j) =0; % convention to avoid marshalling infinite numbers
     182                                        else                                           
     183                                                self.frequencies((i-1)*2*self.n_temporal_iterations +j) = j*log(2)/self.time(i)/2/pi;
     184                                        end
    167185                                end
    168186                        end

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

    r26310 r26800  
    3535                end % }}}
    3636                function md = checkconsistency(self,md,solution,analyses) % {{{
    37 
    3837                        %Early return
    3938                        if ~ismember('HydrologyShreveAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.ishydrology==0), return; end

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

    r26059 r26800  
    4444                end % }}}
    4545                function marshall(self,prefix,md,fid) % {{{
    46                         WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
     46                        WriteData(fid,prefix,'name','md.hydrology.model','data',6,'format','Integer');
    4747                        WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    4848                        outputs = self.requested_outputs;

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

    r26301 r26800  
    1010        properties (SetAccess=public)
    1111
    12                 %regular love numbers:
    13                 h           = []; %provided by PREM model
    14                 k           = []; %idem
    15                 l           = []; %idem
     12                %loading love numbers:
     13                h               = []; %provided by PREM model
     14                k               = []; %idem
     15                l               = []; %idem
    1616               
    1717                %tidal love numbers for computing rotational feedback:
    18                 th          = [];
    19                 tk          = [];
    20                 tl          = [];
    21                 tk2secular  = 0; %deg 2 secular number.
     18                th              = [];
     19                tk              = [];
     20                tl              = [];
     21                tk2secular      = 0; %deg 2 secular number.
     22                pmtf_colinear   = [];
     23                pmtf_ortho      = [];
    2224
    2325                %time/frequency for visco-elastic love numbers
    … …  
    4446                        fielddisplay(self,'tl','tidal load Love number (deg 2)');
    4547                        fielddisplay(self,'tk2secular','secular fluid Love number');
     48                        fielddisplay(self,'pmtf_colinear','Colinear component of the Polar Motion Transfer Function (e.g. x-motion due to x-component perturbation of the inertia tensor)');
     49                        fielddisplay(self,'pmtf_ortho','Orthogonal component of the Polar Motion Transfer Function (couples x and y components, only used for Chandler Wobble)');
     50
    4651
    4752                        fielddisplay(self,'istime','time (default: 1) or frequency love numbers (0)');
    … …  
    6267                        self.tk2secular=0.942;
    6368
     69                        self.pmtf_colinear=0.0;
     70                        self.pmtf_ortho=0.0;
     71                        if maxdeg>=2
     72                                self.pmtf_colinear= (1.0+self.k(3,:))/(1.0-self.tk(3,:)/self.tk2secular); %valid only for elastic regime, not viscous. Also neglects chandler wobble
     73                                self.pmtf_ortho= 0.0;
     74                        end
    6475                        %time:
    6576                        self.istime=1; %temporal love numbers by default
    … …  
    8293                        md = checkfield(md,'fieldname','solidearth.lovenumbers.tl','NaN',1,'Inf',1);
    8394                        md = checkfield(md,'fieldname','solidearth.lovenumbers.tk2secular','NaN',1,'Inf',1);
     95                        md = checkfield(md,'fieldname','solidearth.lovenumbers.pmtf_colinear','NaN',1,'Inf',1);
     96                        md = checkfield(md,'fieldname','solidearth.lovenumbers.pmtf_ortho','NaN',1,'Inf',1);
    8497                        md = checkfield(md,'fieldname','solidearth.lovenumbers.timefreq','NaN',1,'Inf',1);
    8598                        md = checkfield(md,'fieldname','solidearth.lovenumbers.istime','NaN',1,'Inf',1,'values',[0 1]);
    … …  
    91104
    92105                        ntf=length(self.timefreq);
    93                         if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
     106                        if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf | size(self.pmtf_colinear,2) ~= ntf | size(self.pmtf_ortho,2) ~= ntf),
    94107                                error('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector');
     108                        end
     109
     110                        if self.istime && self.timefreq(1)~=0
     111                                error('temporal love numbers must start with elastic response, i.e timefreq(1)=0');
    95112                        end
    96113
    … …  
    108125                        WriteData(fid,prefix,'object',self,'fieldname','tk','name','md.solidearth.lovenumbers.tk','format','DoubleMat','mattype',1);
    109126                        WriteData(fid,prefix,'object',self,'fieldname','tl','name','md.solidearth.lovenumbers.tl','format','DoubleMat','mattype',1);
     127                        WriteData(fid,prefix,'object',self,'fieldname','pmtf_colinear','name','md.solidearth.lovenumbers.pmtf_colinear','format','DoubleMat','mattype',1);
     128                        WriteData(fid,prefix,'object',self,'fieldname','pmtf_ortho','name','md.solidearth.lovenumbers.pmtf_ortho','format','DoubleMat','mattype',1);
    110129                        WriteData(fid,prefix,'object',self,'data',self.tk2secular,'fieldname','lovenumbers.tk2secular','format','Double');
    111130

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

    r26358 r26800  
    264264                                        end
    265265                                        ind=find(md.materials.issolid==0);
    266                                         if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
    267                                                 error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
    268                                         end
     266                                        %if sum(ismember(diff(ind),1)>=1) %if there are at least two consecutive indices that contain issolid = 0
     267                                        %               error(['Fluid layers detected at layers #', num2str(ind'), ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'])
     268                                        %end
    269269
    270270                                case 'hydro'

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

    r26358 r26800  
    8383md.transient.ismasstransport=1;
    8484md.transient.isslc=1;
    85 md.solidearth.requested_outputs={'Sealevel','Bed'};
     85md.solidearth.requested_outputs={'Sealevel','Bed','SealevelBarystaticIceLoad', 'SealevelBarystaticIceArea', 'SealevelBarystaticIceWeights'};
    8686
    8787% max number of iteration reverted back to 10 (i.e., the original default value)

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

    r26358 r26800  
    11%Test Name: EarthSlc_rotationalFeedback
     2
    23
    34%mesh earth:

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

    r26358 r26800  
    362362if plotting,
    363363        flags=ones(sl.earth.mesh.numberofelements,1);
    364         for i=1:length(sl.eltransitions),
     364        for i=1:length(sl.eltransitions)
    365365                flags(sl.eltransitions{i})=i;
    366366        end
    … …  
    422422md.solidearth.settings.viscous=0;
    423423md.solidearth.requested_outputs= {'default',...
    424         'DeltaIceThickness','Sealevel','SealevelUGrd',...
    425         'SealevelchangeBarystaticMask','SealevelchangeBarystaticOceanMask'};
     424        'DeltaIceThickness','Sealevel','Bed',...
     425        'SealevelBarystaticIceMask','SealevelBarystaticOceanMask'};
    426426md=solve(md,'Transient');
    427427Seustatic=md.results.TransientSolution.Sealevel;

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

    r26358 r26800  
    1515longe=md.mesh.long;
    1616time=0:0.5:5;
     17%The offline solution pattern is a degree (2,1) spherical harmonic
    1718md.solidearth.external=offlinesolidearthsolution;
    1819md.solidearth.external.displacementup=.5*sind(late).*cosd(late).*cosd(longe) .*time;

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

    r26358 r26800  
    9999longe=md.mesh.long;
    100100time=0:1;
    101 Y22=1.5*(1.+cosd(2.*late)).*cosd(2.*longe);
     101Y22=1.5*(1.+cosd(2.*late)).*cosd(2.*longe);%The additional solidearth signal is a degree (2,2) spherical harmonic
    102102md.solidearth.external=additionalsolidearthsolution;
    103103md.solidearth.external.displacementup=0.5*Y22 .*time;

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

    r26358 r26800  
    9898eus=md.results.TransientSolution.Bslc;
    9999slc=md.results.TransientSolution.Sealevel;
    100 moixz=md.results.TransientSolution.SealevelInertiaTensorXZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
    101 moiyz=md.results.TransientSolution.SealevelInertiaTensorYZ / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
    102 moizz=md.results.TransientSolution.SealevelInertiaTensorZZ / ( -(1+load_love_k2)/moi_p);
     100moixz=md.results.TransientSolution.SealevelchangePolarMotionX / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
     101moiyz=md.results.TransientSolution.SealevelchangePolarMotionY / (1/(1-tide_love_k2/tide_love_k2secular) * (1+load_love_k2)/(moi_p-moi_e) );
     102moizz=md.results.TransientSolution.SealevelchangePolarMotionZ / ( -(1+load_love_k2)/moi_p);
    103103
    104104areaice=md.results.TransientSolution.SealevelBarystaticIceArea;
     105areaice(isnan(areaice))=0;
    105106loadice=md.results.TransientSolution.SealevelBarystaticIceLoad;
    106 
    107 % analytical moi => just checking FOR ICE only!!! {{{
    108 % ...have to mute ** slc induced MOI in Tria.cpp ** prior to the comparison
    109107rad_e = md.solidearth.planetradius;
    110108
    … …  
    113111moi_xz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*cos(lon));
    114112moi_yz = sum(-loadice.*areaice.*rad_e^2.*sin(lat).*cos(lat).*sin(lon));
    115 moi_zz = sum(-loadice.*areaice.*rad_e^2.*(1.0-sin(lat).^2));
    116 theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moi_zz]
     113moi_zz = sum(-loadice.*areaice.*rad_e^2.*(-1.0/3.0+sin(lat).^2));
     114theoretical_value_check=[moixz/moi_xz moiyz/moi_yz moizz/moi_zz] % should yield [1.0 1.0 1.0]
    117115% }}}
    118116

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

    r25300 r26800  
    55md=parameterize(md,'../Par/GiaIvinsBenchmarksAB.par');
    66
     7%GIA Ivins, 2 layer model.
     8md.solidearth.settings.grdmodel=2;
     9md.solidearth.settings.isgrd=1;
     10md.initialization.sealevel=zeros(md.mesh.numberofvertices,1);
     11
    712%% indicate what you want to compute
    8 md.gia.cross_section_shape=1;    % for square-edged x-section
     13md.solidearth.settings.cross_section_shape=1;    % for square-edged x-section
     14md.solidearth.settings.grdocean=0;  %do not compute sea level, only deformation
     15md.solidearth.settings.sealevelloading=0;  %do not compute sea level, only deformation
    916
    1017% evaluation time (termed start_time)
    11 md.timestepping.start_time=2002100; % after 2 kyr of deglaciation
     18
     19md.timestepping.time_step=2002100; % after 2 kyr of deglaciation
     20md.timestepping.start_time=-md.timestepping.time_step; % need one time step before t=0 to generate a thickness change at t=0
    1221md.timestepping.final_time=2500000; % 2,500 kyr
    1322
    1423%% define loading history {{{
    15 md.geometry.thickness=[...
     24md.masstransport.spcthickness=[...
    1625        [md.geometry.thickness*0.0; 0.0],...
    1726        [md.geometry.thickness; 1000],...
    1827        [md.geometry.thickness; 2000000],...
    1928        [md.geometry.thickness*0.0; 2000100],...
    20         [md.geometry.thickness*0.0; md.timestepping.start_time],...
     29        [md.geometry.thickness*0.0; md.timestepping.start_time+2*md.timestepping.time_step],...
    2130        ];
    2231% }}}
     32
     33md.geometry.bed=zeros(md.mesh.numberofvertices,1);
    2334
    2435% find out elements that have zero loads throughout the loading history.
    … …  
    2637md.mask.ice_levelset(pos)=1; % no ice
    2738
     39%Physics:
     40md.transient.issmb=0;
     41md.transient.isstressbalance=0;
     42md.transient.isthermal=0;
     43md.transient.ismasstransport=1;
     44md.transient.isslc=1;
     45
    2846md.cluster=generic('name',oshostname(),'np',3);
    2947md.verbose=verbose('1111111');
     48md.solidearth.requested_outputs={'Sealevel','BedGRD'};
    3049
    3150%% solve for GIA deflection
    32 md=solve(md,'Gia');
     51md=solve(md,'Transient');
    3352
    3453%Test Name: GiaIvinsBenchmarksAB2dA1
    35 U_AB2dA1 = md.results.GiaSolution.UGia;
    36 URate_AB2dA1 = md.results.GiaSolution.UGiaRate;
     54U_AB2dA1 = md.results.TransientSolution.BedGRD;
     55%URate_AB2dA1 = md.results.TransientSolution.UGiaRate;
    3756
    3857%Test Name: GiaIvinsBenchmarksAB2dA2
    3958%% different evaluation time. {{{
    40 md.timestepping.start_time=2005100; % after 5 kyr of deglaciation
    41 md.geometry.thickness(end,end) = md.timestepping.start_time;
     59md.timestepping.time_step=2005100;% after 5 kyr of deglaciation
     60md.timestepping.start_time=-md.timestepping.time_step; % need one time step before t=0 to generate a thickness change at t=0
     61md.masstransport.spcthickness(end,end) = md.timestepping.start_time+2*md.timestepping.time_step;
    4262
    43 md=solve(md,'Gia');
     63md=solve(md,'Transient');
    4464
    45 U_AB2dA2 = md.results.GiaSolution.UGia;
    46 URate_AB2dA2 = md.results.GiaSolution.UGiaRate;
     65U_AB2dA2 = md.results.TransientSolution.BedGRD;
     66%URate_AB2dA2 = md.results.TransientSolution.BedGRDRate;
    4767% }}}
    4868
    4969%Test Name: GiaIvinsBenchmarksAB2dA3
    5070%% different evaluation time. {{{
    51 md.timestepping.start_time=2010100; % after 10 kyr of deglaciation
    52 md.geometry.thickness(end,end) = md.timestepping.start_time;
     71md.timestepping.time_step=2010100;% after 10 kyr of deglaciation
     72md.timestepping.start_time=-md.timestepping.time_step; % need one time step before t=0 to generate a thickness change at t=0
     73md.masstransport.spcthickness(end,end) = md.timestepping.start_time+2*md.timestepping.time_step;
    5374
    54 md=solve(md,'Gia');
     75md=solve(md,'Transient');
    5576
    56 U_AB2dA3 = md.results.GiaSolution.UGia;
    57 URate_AB2dA3 = md.results.GiaSolution.UGiaRate;
     77U_AB2dA3 = md.results.TransientSolution.BedGRD;
     78%URate_AB2dA3 = md.results.TransientSolution.UGiaRate;
    5879% }}}
    5980

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

    r26358 r26800  
    77
    88md=model();
    9 md.cluster=generic('name',oshostname(),'np',1);
     9md.cluster=generic('name',oshostname(),'np',8);
    1010
    11 % set validation=1 for comparing against the Spada benchark.
     11% set validation=1 for comparing against the Spada benchmark.
    1212validation=0;
    1313
    … …  
    1717
    1818md.verbose=verbose('all');
     19md.verbose=verbose('1111111111111111');
    1920cst=365.25*24*3600*1000;
    2021
    … …  
    3637
    3738md.love.allow_layer_deletion=1;
    38 md.love.frequencies=([0]*2*pi)'/cst;
     39md.love.frequencies=[0; (logspace(-6,3, 1000))'/cst];
    3940md.love.nfreq=length(md.love.frequencies);
    4041md.love.sh_nmin=1;
    41 md.love.sh_nmax=256;
     42md.love.sh_nmax=1000;
    4243md.love.underflow_tol=1e-20;
     44md.love.pw_threshold=1e-3;
    4345md.love.Gravitational_Constant=6.6732e-11;
    44 md.love.integration_steps_per_layer=200;
     46md.love.integration_steps_per_layer=100;
     47md.love.allow_layer_deletion=1;
     48md.love.forcing_type=11;
     49md.love.chandler_wobble=0;
     50md.love.complex_computation=0;
    4551
    4652md.love.istemporal=1;
    4753md.love.n_temporal_iterations=8;
    48 %md.love.time=(logspace(-4,5, 2))'*cst;
    49 md.love.time=(logspace(-1,2, 50))'*cst;
     54%md.love.time=(logspace(-6,5, 2))'*cst;
     55md.love.time=[0; (logspace(-3,5, 99))'*cst];
     56
     57%md.love.time=(linspace(1/12,10, 10*12))'*cst/1e3;
    5058md.love.love_kernels=1;
    5159if md.love.istemporal
    … …  
    5563md=solve(md,'lv');
    5664
    57 ht2=md.results.LoveSolution.LoveHr;
    58 lt2=md.results.LoveSolution.LoveLr;
    59 kt2=md.results.LoveSolution.LoveKr;
     65ht=md.results.LoveSolution.LoveHt';
     66lt=md.results.LoveSolution.LoveLt';
     67kt=md.results.LoveSolution.LoveKt';
    6068t=md.love.time/cst*1e3;
    6169
    … …  
    6977field_tolerances={2.0e-8,2.0e-8,2.0e-8};
    7078field_values={...
    71         (md.results.LoveSolution.LoveHr(:,1)),...
    72         (md.results.LoveSolution.LoveKr(:,1)),...
    73         (md.results.LoveSolution.LoveLr(:,1)),...
     79        (md.results.LoveSolution.LoveHt(:,1)),...
     80        (md.results.LoveSolution.LoveKt(:,1)),...
     81        (md.results.LoveSolution.LoveLt(:,1)),...
    7482        };
     83
     84return
     85
     86if false
     87addpath('../../../../invGIA/Spada_benchmark/')
     88load spada.mat
     89s_weak=[9 12 15];
     90
     91hspada(:,s_weak)=[];
     92kspada(:,s_weak)=[];
     93lspada(:,s_weak)=[];
     94
     95hts=zeros(length(t),md.love.sh_nmax+1);
     96lts=zeros(length(t),md.love.sh_nmax+1);
     97kts=zeros(length(t),md.love.sh_nmax+1);
     98for d=max(2,md.love.sh_nmin):md.love.sh_nmax
     99hts(:,d+1)=hspada(d-1,2);
     100lts(:,d+1)=lspada(d-1,2);
     101kts(:,d+1)=kspada(d-1,2);
     102for mo=1:9
     103hts(:,d+1)=hts(:,d+1)-hspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     104lts(:,d+1)=lts(:,d+1)-lspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     105kts(:,d+1)=kts(:,d+1)-kspada(d-1,3+mo)./sspada(d-1,1+mo).*(1-exp(t/1e3*sspada(d-1,1+mo)));
     106end
     107end
     108close all
     109
     110if md.love.sh_nmin==md.love.sh_nmax
     111subplot(2,3,1)
     112plot(t/cst*1e3,ht(:,md.love.sh_nmin+1),'x-',t/cst*1e3,hts);
     113set(gca, 'xscale', 'log'); shading flat;
     114title('h')
     115
     116subplot(2,3,2)
     117plot(t/cst*1e3,lt(:,md.love.sh_nmin+1),'x-',t/cst*1e3,lts);
     118set(gca, 'xscale', 'log'); shading flat;
     119title('l')
     120
     121subplot(2,3,3)
     122plot(t/cst*1e3,kt(:,md.love.sh_nmin+1),'x-',t/cst*1e3,kts);
     123set(gca, 'xscale', 'log'); shading flat;
     124title('k')
     125
     126subplot(2,3,4)
     127plot(t/cst*1e3,log10(abs((ht-hts)./hts))');
     128set(gca, 'xscale', 'log'); shading flat;
     129
     130subplot(2,3,5)
     131plot(t/cst*1e3,log10(abs((lt-lts)./lts))');
     132set(gca, 'xscale', 'log'); shading flat;
     133title('l')
     134
     135subplot(2,3,6)
     136plot(t/cst*1e3,log10(abs((kt-kts)./kts))');
     137set(gca, 'xscale', 'log'); shading flat;
     138title('k')
     139else
     140
     141[T,N]=meshgrid(t/cst*1e3,0:md.love.sh_nmax);
     142subplot(1,3,1)
     143pcolor(T,N,log10(abs((ht-hts)./hts))');
     144set(gca, 'xscale', 'log'); shading flat;colorbar;
     145title('h')
     146
     147subplot(1,3,2)
     148pcolor(T,N,log10(abs((lt-lts)./lts))');
     149set(gca, 'xscale', 'log'); shading flat;colorbar;
     150title('l')
     151
     152subplot(1,3,3)
     153pcolor(T,N,log10(abs((kt-kts)./kts))');
     154set(gca, 'xscale', 'log'); shading flat;colorbar;
     155title('k')
     156end
     157
    75158
    76159% validate elastic loading solutions against the Spada benchmark. {{{
    … …  
    100183        %
    101184end
     185
     186end
    102187% }}}
    103 
    104 md.love.frequencies=([1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
    105 md.love.nfreq=length(md.love.frequencies);
    106 md.love.sh_nmax=256;
    107 md.materials.burgers_mu=md.materials.lame_mu;
    108 md.materials.burgers_viscosity=md.materials.viscosity;
    109 
    110 md=solve(md,'lv');
    111 
    112 %Fields and tolerances to track changes
    113 field_names     ={field_names{:},'LoveH_loading_realpart','LoveK_loading_realpart','LoveL_loading_realpart','LoveH_loading_imagpart','LoveK_loading_imagpart','LoveL_loading_imagpart'};
    114 field_tolerances={field_tolerances{:},5e-7,5e-7,5e-7,5e-7,5e-7,5e-7};
    115 field_values={field_values{:},...
    116         (md.results.LoveSolution.LoveHr(:,:)),...
    117         (md.results.LoveSolution.LoveKr(:,:)),...
    118         (md.results.LoveSolution.LoveLr(:,:)),...
    119         (md.results.LoveSolution.LoveHi(:,:)),...
    120         (md.results.LoveSolution.LoveKi(:,:)),...
    121         (md.results.LoveSolution.LoveLi(:,:)),...
    122         };
    123 
    124 md.love.forcing_type=9;
    125 md.love.sh_nmin=2;
    126 md.love.frequencies=([0 1e-3 1e-2 1e-1 1 -1e-3 -1e-2 -1e-1 -1]*2*pi)'/cst;
    127 md.love.nfreq=length(md.love.frequencies);
    128 md=solve(md,'lv');
    129 
    130 % validate elastic tidal solutions against the Spada benchmark. {{{
    131 if validation
    132         spada_solutions = load('spada_elastic_tidal_deg_h_l_k');
    133         spada_d = spada_solutions(:,1);
    134         spada_h = spada_solutions(:,2);
    135         spada_l = spada_solutions(:,3);
    136         spada_k = spada_solutions(:,4);
    137 
    138         %rename ISSM solutions. 
    139         issm_d = [2:md.love.sh_nmax]; 
    140         issm_h = md.results.LoveSolution.LoveHr(3:end,1);
    141         issm_l = md.results.LoveSolution.LoveLr(3:end,1);
    142         issm_k = md.results.LoveSolution.LoveKr(3:end,1);
    143 
    144         % relative difference for each degree, except for zero and one.
    145         diff_h = 1 - issm_h./spada_h; 
    146         diff_l = 1 - issm_l./spada_l; 
    147         diff_k = 1 - issm_k./spada_k; 
    148 
    149         % k seems to have larger relative difference at higher degrees, although values approach zero.
    150         figure
    151         plot(spada_d,issm_k./spada_k);
    152 
    153         figure
    154         plot(spada_d,[diff_h diff_l diff_k]); grid on;
    155         legend('h','l','k'); title('tidal love');
    156 else
    157         %
    158 end
    159 % }}}
    160 
    161 %tidal love numbers
    162 field_names     ={field_names{:},'LoveH_tidal_elastic','LoveK_tidal_elastic','LoveL_tidal_elastic','LoveH_tidal_realpart','LoveK_tidal_realpart','LoveL_tidal_realpart','LoveH_tidal_imagpart','LoveK_tidal_imagpart','LoveL_tidal_imagpart'};
    163 field_tolerances={field_tolerances{:},8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6,8e-6};
    164 field_values={field_values{:},...
    165         (md.results.LoveSolution.LoveHr(:,1)),...
    166         (md.results.LoveSolution.LoveKr(:,1)),...
    167         (md.results.LoveSolution.LoveLr(:,1)),...
    168         (md.results.LoveSolution.LoveHr(:,2:end)),...
    169         (md.results.LoveSolution.LoveKr(:,2:end)),...
    170         (md.results.LoveSolution.LoveLr(:,2:end)),...
    171         (md.results.LoveSolution.LoveHi(:,2:end)),...
    172         (md.results.LoveSolution.LoveKi(:,2:end)),...
    173         (md.results.LoveSolution.LoveLi(:,2:end)),...
    174         };
    175 
    176 %Many layers PREM-based model
    177 %data=load('../Data/PREM_500layers');
    178 %md.love.sh_nmin=1;
    179 %md.materials.radius=data(2:end-2,1);
    180 %md.materials.density=data(3:end-2,2);
    181 %md.materials.lame_lambda=data(3:end-2,3);
    182 %md.materials.lame_mu=data(3:end-2,4);
    183 %md.materials.issolid=data(3:end-2,4)>0;
    184 %ind=find(md.materials.issolid==0);
    185 %md.materials.density(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.density(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
    186 %md.materials.lame_lambda(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_lambda(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)+1).^3);
    187 %md.materials.lame_mu(ind(1))=sum((md.materials.radius(ind+1).^3-md.materials.radius(ind).^3).*md.materials.lame_mu(ind))/(md.materials.radius(ind(end)+1).^3-md.materials.radius(ind(1)).^3);
    188 %md.materials.radius(ind(2:end)+1)=[];
    189 %md.materials.density(ind(2:end))=[];
    190 %md.materials.lame_lambda(ind(2:end))=[];
    191 %md.materials.lame_mu(ind(2:end))=[];
    192 %md.materials.issolid(ind(2:end))=[];
    193 %md.materials.viscosity=10.^interp1([0 3479e3 3480e3 3680e3 5720e3 5800e3 6270e3 6371e3], log10([1e8 1e8 5e21 1e23 1e22 1e20 1e21 1e40]), md.materials.radius(2:end),'PCHIP');
    194 %md.materials.viscosity=md.materials.viscosity.*md.materials.issolid;
    195 %md.materials.burgers_mu=md.materials.lame_mu;
    196 %md.materials.burgers_viscosity=md.materials.viscosity;
    197 %md.materials.rheologymodel=md.materials.issolid*0;
    198 %md.love.forcing_type=11;
    199 %md.materials.numlayers=length(md.materials.viscosity);
    200 %md=solve(md,'lv');
    201 %
    202 %field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
    203 %field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
    204 %field_values={field_values{:},...
    205 %       (md.results.LoveSolution.LoveHr(:,1)),...
    206 %       (md.results.LoveSolution.LoveKr(:,1)),...
    207 %       (md.results.LoveSolution.LoveLr(:,1)),...
    208 %       (md.results.LoveSolution.LoveHr(:,2:end)),...
    209 %       (md.results.LoveSolution.LoveKr(:,2:end)),...
    210 %       (md.results.LoveSolution.LoveLr(:,2:end)),...
    211 %       (md.results.LoveSolution.LoveHi(:,2:end)),...
    212 %       (md.results.LoveSolution.LoveKi(:,2:end)),...
    213 %       (md.results.LoveSolution.LoveLi(:,2:end)),...
    214 %       };
    215 
    216 %Model VSS96 from Vermeersen, L.L.A., Sabadini, R. & Spada, G., 1996a. Analytical visco-elastic relaxation models, Geophys. Res. Lett., 23, 697–700.
    217 md.materials.radius=[10, 1222.5, 3480., 3600., 3630.5, 3700., 3900., 4000., 4200., 4300., 4500., 4600., 4800., 4900., 5100., 5200., 5400., 5500., 5600.5, 5650., 5701., 5736., 5771.5, 5821., 5951., 5970.5, 6016., 6061., 6150.5, 6151.5, 6251., 6371.]'*1e3;
    218 md.materials.lame_mu=[1e-5, 0., 2.933, 2.8990002, 2.8550003, 2.7340002, 2.675, 2.559, 2.502, 2.388, 2.331, 2.215, 2.157, 2.039, 1.979, 1.8560001, 1.794, 1.73, 1.639, 1.2390001, 1.224, 1.21, 1.128, 0.97700006, 0.906, 0.79, 0.773, 0.741, 0.656, 0.665, 0.602]'*1e11;
    219 md.materials.density=[10925., 10925., 5506.42, 5491.45, 5456.57, 5357.06, 5307.24, 5207.13, 5156.69, 5054.69, 5002.99, 4897.83, 4844.22, 4734.6, 4678.44, 4563.07, 4503.72, 4443.16, 4412.41, 3992.14, 3983.99, 3975.84, 3912.82, 3786.78, 3723.78, 3516.39, 3489.51, 3435.78, 3359.5, 3367.1, 3184.3]';
    220 md.materials.viscosity=[0., 0., 7.999999999999999E+21, 8.5E+21, 8.999999999999999E+21, 3.E+22, 4.E+22, 5.0000000000000004E+22, 6.E+22, 5.0000000000000004E+22, 4.5E+22, 3.E+22, 2.5000000000000002E+22, 1.7999999999999998E+22, 1.3E+22, 7.999999999999999E+21, 6.999999999999999E+21, 6.5E+21, 6.E+21, 5.5E+21, 5.E+21, 4.4999999999999995E+21, 3.9999999999999995E+21, 2.5E+21, 1.9999999999999997E+21, 1.5E+21, 9.999999999999999E+20, 6.E+20, 5.5000000000000007E+20, 2.E+20, 1.E40]';
    221 md.materials.lame_lambda=md.materials.lame_mu*0+5e14;
    222 md.materials.issolid=md.materials.lame_mu>0;
    223 md.materials.numlayers=length(md.materials.lame_mu);
    224 md.materials.burgers_mu=md.materials.lame_mu;
    225 md.materials.burgers_viscosity=md.materials.viscosity;
    226 md.materials.rheologymodel=md.materials.issolid*0;
    227 md.love.forcing_type=11;
    228 md.love.sh_nmin=1;
    229 md.love.sh_nmax=100;
    230 md=solve(md,'lv');
    231 md.love.frequencies=([0 1e-3 1e-2 1 -1e-3 -1e-2 -1]*2*pi)'/cst;
    232 md.love.nfreq=length(md.love.frequencies);
    233 
    234 field_names     ={field_names{:},'LoveH_loadingVSS96_elastic','LoveK_loadingVSS96_elastic','LoveL_loadingVSS96_elastic','LoveH_loadingVSS96_realpart','LoveK_loadingVSS96_realpart','LoveL_loadingVSS96_realpart','LoveH_loadingVSS96_imagpart','LoveK_loadingVSS96_imagpart','LoveL_loadingVSS96_imagpart'};
    235 field_tolerances={field_tolerances{:},2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6,2e-6};
    236 field_values={field_values{:},...
    237         (md.results.LoveSolution.LoveHr(:,1)),...
    238         (md.results.LoveSolution.LoveKr(:,1)),...
    239         (md.results.LoveSolution.LoveLr(:,1)),...
    240         (md.results.LoveSolution.LoveHr(:,2:end)),...
    241         (md.results.LoveSolution.LoveKr(:,2:end)),...
    242         (md.results.LoveSolution.LoveLr(:,2:end)),...
    243         (md.results.LoveSolution.LoveHi(:,2:end)),...
    244         (md.results.LoveSolution.LoveKi(:,2:end)),...
    245         (md.results.LoveSolution.LoveLi(:,2:end)),...
    246         };

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

    r25356 r26800  
    99% for volumetric potential
    1010        md=model();
     11        md.groundingline.migration='None';
    1112       
    1213        md.materials=materials('litho');
     14        cst=365.25*24*3600*1000;
    1315
    1416        md.materials.numlayers = 40;
    … …  
    1719        md.materials.density=zeros(md.materials.numlayers,1)+5511;
    1820        md.materials.lame_mu=zeros(md.materials.numlayers,1)+0.75e11;
    19         md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
    2021        md.materials.lame_lambda=zeros(md.materials.numlayers,1)+5e17;
    2122        md.materials.issolid=ones(md.materials.numlayers,1);
    22         md.materials.isburgers=zeros(md.materials.numlayers,1);
     23        md.materials.rheologymodel=zeros(md.materials.numlayers,1);
     24
     25        %the following isn't used here but needs to have arrays of consistent size with the rest of the materials
     26        md.materials.viscosity=zeros(md.materials.numlayers,1)+1e21;
    2327        md.materials.burgers_mu=md.materials.lame_mu/3;
    2428        md.materials.burgers_viscosity=md.materials.viscosity/10;
     29        md.materials.ebm_alpha= ones(md.materials.numlayers,1)*.9;
     30        md.materials.ebm_delta= ones(md.materials.numlayers,1)*0.2;
     31        md.materials.ebm_taul= ones(md.materials.numlayers,1)*54*60; %54min
     32        md.materials.ebm_tauh= ones(md.materials.numlayers,1)*18.6*cst/1e3; %18.6yr
     33
    2534
    2635        md.materials.radius =  linspace(10e3,6371e3,md.materials.numlayers+1)';
    … …  
    2938        md.love.allow_layer_deletion=1;
    3039        md.love.frequencies=0;
    31         md.love.nfreq=length(md.love.frequencies); % TODO: Why are we grabbing length of a scalar here?
     40        md.love.nfreq=1;
     41        md.love.istemporal=0;
    3242
    3343        md.love.sh_nmin = 2;
    … …  
    3646
    3747        md.miscellaneous.name='kernels';
    38         md.cluster=generic('name',oshostname(),'np',1);
     48        md.cluster=generic('name',oshostname(),'np',3);
    3949        md.verbose=verbose('111111101');
    4050       
    … …  
    4353        % save yi's for all layers except for the inner-most one, at select degrees.
    4454        degrees = [2 20 200];  % we archive solutions for degrees 2, 20, 200
    45 
     55        kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
    4656        % extract love kernels {{{
    4757        % degree 2.
    48         y1_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,1)));
    49         y2_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,2)));
    50         y3_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,3)));
    51         y4_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,4)));
    52         y5_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,5)));
    53         y6_tidal_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,6)));
     58        y1_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
     59        y2_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
     60        y3_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
     61        y4_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
     62        y5_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
     63        y6_tidal_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
    5464
    5565        % degree 20.
    56         y1_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,1)));
    57         y2_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,2)));
    58         y3_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,3)));
    59         y4_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,4)));
    60         y5_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,5)));
    61         y6_tidal_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,6)));
     66        y1_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
     67        y2_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
     68        y3_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
     69        y4_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
     70        y5_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
     71        y6_tidal_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
    6272
    6373        % degree 200.
    64         y1_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,1)));
    65         y2_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,2)));
    66         y3_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,3)));
    67         y4_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,4)));
    68         y5_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,5)));
    69         y6_tidal_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,6)));
     74        y1_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
     75        y2_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
     76        y3_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
     77        y4_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
     78        y5_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
     79        y6_tidal_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
    7080        % }}}
    7181
    … …  
    100110        depth = (max(param.radius)-param.radius)/1000; % km.
    101111
    102         y1 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
    103         y2 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
    104         y3 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
    105         y4 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
    106         y5 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,5)));
    107         y6 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,6)));
     112        kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
     113        y1 = squeeze(kernels(:,1,:,1));
     114        y2 = squeeze(kernels(:,1,:,2));
     115        y3 = squeeze(kernels(:,1,:,3));
     116        y4 = squeeze(kernels(:,1,:,4));
     117        y5 = squeeze(kernels(:,1,:,5));
     118        y6 = squeeze(kernels(:,1,:,6));
    108119       
    109120        set(0,'DefaultAxesFontSize',16,'DefaultTextFontSize',15,'DefaultAxesLineWidth',1,...
    … …  
    175186        md.love.forcing_type = 11;
    176187        md=solve(md,'lv');
     188        kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
    177189
    178190        % extract love kernels {{{
    179191        % degree 2.
    180         y1_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,1)));
    181         y2_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,2)));
    182         y3_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,3)));
    183         y4_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,4)));
    184         y5_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,5)));
    185         y6_loading_degree002 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(1)+1,1,2:end,6)));
     192        y1_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,1));
     193        y2_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,2));
     194        y3_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,3));
     195        y4_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,4));
     196        y5_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,5));
     197        y6_loading_degree002 = squeeze(kernels(degrees(1)+1,1,2:end,6));
    186198
    187199        % degree 20.
    188         y1_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,1)));
    189         y2_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,2)));
    190         y3_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,3)));
    191         y4_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,4)));
    192         y5_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,5)));
    193         y6_loading_degree020 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(2)+1,1,2:end,6)));
     200        y1_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,1));
     201        y2_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,2));
     202        y3_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,3));
     203        y4_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,4));
     204        y5_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,5));
     205        y6_loading_degree020 = squeeze(kernels(degrees(2)+1,1,2:end,6));
    194206
    195207        % degree 200.
    196         y1_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,1)));
    197         y2_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,2)));
    198         y3_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,3)));
    199         y4_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,4)));
    200         y5_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,5)));
    201         y6_loading_degree200 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(degrees(3)+1,1,2:end,6)));
     208        y1_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,1));
     209        y2_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,2));
     210        y3_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,3));
     211        y4_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,4));
     212        y5_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,5));
     213        y6_loading_degree200 = squeeze(kernels(degrees(3)+1,1,2:end,6));
    202214        % }}}
    203215
    … …  
    225237        depth = (max(param.radius)-param.radius)/1000; % km.
    226238
    227         y1 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,1)));
    228         y2 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,2)));
    229         y3 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,3)));
    230         y4 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,4)));
    231         y5 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,5)));
    232         y6 = squeeze(cell2mat(md.results.LoveSolution.LoveKernelsReal(:,1,:,6)));
     239        kernels=reshape(md.results.LoveSolution.LoveKernels, [md.love.sh_nmax+1 md.love.nfreq md.materials.numlayers+1 6]);
     240        y1 = squeeze(kernels(:,1,:,1));
     241        y2 = squeeze(kernels(:,1,:,2));
     242        y3 = squeeze(kernels(:,1,:,3));
     243        y4 = squeeze(kernels(:,1,:,4));
     244        y5 = squeeze(kernels(:,1,:,5));
     245        y6 = squeeze(kernels(:,1,:,6));
    233246       
    234247        set(0,'DefaultAxesFontSize',16,'DefaultTextFontSize',15,'DefaultAxesLineWidth',1,...

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

    r26358 r26800  
    3737%still use the lovenumbers constructor to initialize the fields:
    3838load ../Data/lnb_temporal.mat
    39 md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
    40 md.solidearth.lovenumbers.timefreq=[0];
     39%md.solidearth.lovenumbers=lovenumbers('maxdeg',1000);
     40%md.solidearth.lovenumbers.timefreq=[0];
    4141
    42 md.solidearth.lovenumbers.h=[h1'];
    43 md.solidearth.lovenumbers.k=[k1'];
    44 md.solidearth.lovenumbers.l=[l1'];
     42md.solidearth.lovenumbers.h=[ht];
     43md.solidearth.lovenumbers.k=[kt];
     44md.solidearth.lovenumbers.l=[lt];
    4545%md.solidearth.lovenumbers.h=repmat(md.solidearth.lovenumbers.h,1,100);
    4646%md.solidearth.lovenumbers.k=repmat(md.solidearth.lovenumbers.k,1,100);
    4747%md.solidearth.lovenumbers.l=repmat(md.solidearth.lovenumbers.l,1,100);
    48 md.solidearth.lovenumbers.th=repmat(md.solidearth.lovenumbers.th,1,101);
    49 md.solidearth.lovenumbers.tk=repmat(md.solidearth.lovenumbers.tk,1,101);
    50 md.solidearth.lovenumbers.tl=repmat(md.solidearth.lovenumbers.tl,1,101);
     48md.solidearth.lovenumbers.th=tht;
     49md.solidearth.lovenumbers.tk=tkt;
     50md.solidearth.lovenumbers.tl=tlt;
     51md.solidearth.lovenumbers.pmtf_colinear=pmtf1;
     52md.solidearth.lovenumbers.pmtf_ortho=pmtf2;
    5153md.solidearth.lovenumbers.timefreq=time;
    5254
    … …  
    9294
    9395%Solution parameters
    94 md.cluster.np=3;
     96md.cluster.np=10;
    9597md.solidearth.settings.reltol=NaN;
    9698md.solidearth.settings.abstol=1e-3;
    … …  
    117119md.solidearth.settings.elastic=1;
    118120md.solidearth.settings.viscous=1;
    119 md.solidearth.settings.rotation=0;
     121md.solidearth.settings.rotation=1;
    120122md=solve(md,'tr');
    121123

  • issm/trunk-jpl/test/Par/GiaIvinsBenchmarksAB.par

    r25290 r26800  
    1414md.geometry.surface=md.geometry.thickness+md.geometry.base;
    1515
     16
     17%GIA parameters specific to Experiments A  and B
     18
     19md.materials=materials('litho','ice');
     20md.materials.numlayers=2;
     21md.materials.radius =  [10 6271e3 6371e3]; %100km litosphere, the rest is mantle material
     22md.materials.density=  [3.38e3 3.36e3];
     23md.materials.lame_mu=  [1.45e11         6.7e10];
     24md.materials.viscosity=[1e21            1e40];
     25
    1626%Ice density used for benchmarking, not 917 kg/m^3
    1727md.materials.rho_ice=1000; %kg m^3
    18 
    19 %GIA parameters specific to Experiments A  and B
    20 md.gia=giaivins();
    21 md.gia.mantle_viscosity=10^21*ones(md.mesh.numberofvertices,1);    %in Pa.s
    22 md.gia.lithosphere_thickness=100*ones(md.mesh.numberofvertices,1); %in km
    23 md.materials.lithosphere_shear_modulus=6.7*10^10;                  %in Pa
    24 md.materials.lithosphere_density=3.36;                             %in g/cm^3
    25 md.materials.mantle_shear_modulus=1.45*10^11;                      %in Pa
    26 md.materials.mantle_density=3.38;                                  %in g/cm^3
    2728
    2829
Note: See TracChangeset for help on using the changeset viewer.