Changeset 27131


Ignore:
Timestamp:
07/01/22 15:15:14 (3 years ago)
Author:
caronlam
Message:

CHG: Memory and runtime optimization of sealevelchange core.

Location:
issm/trunk-jpl/src
Files:
18 edited

Legend:

Unmodified
Added
Removed

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

    r27101 r27131  
    124124        IssmDouble  planetradius=0;
    125125        IssmDouble  planetarea=0;
     126        IssmDouble  constant=0;
     127        IssmDouble  rho_earth;
     128        int             isgrd=0;
    126129        bool            selfattraction=false;
    127130        bool            elastic=false;
    … …  
    137140        int     numoutputs;
    138141        char**  requestedoutputs = NULL;
     142        int* recvcounts = NULL;
     143        int* displs=NULL;
    139144
    140145        /*transition vectors: */
    … …  
    158163        if(isexternal) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.external.nature",SolidearthExternalNatureEnum));
    159164        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.runfrequency",SolidearthSettingsRunFrequencyEnum));
     165        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.degacc",SolidearthSettingsDegreeAccuracyEnum));
    160166        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.reltol",SolidearthSettingsReltolEnum));
    161167        parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.settings.abstol",SolidearthSettingsAbstolEnum));
    … …  
    181187        parameters->AddObject(new DoubleParam(CumBslcHydroEnum,0.0));
    182188        parameters->AddObject(new DoubleParam(CumGmtslcEnum,0.0));
    183 
    184189        /*compute planet area and plug into parameters:*/
    185190        iomodel->FetchData(&planetradius,"md.solidearth.planetradius");
     191        iomodel->FetchData(&rho_earth,"md.materials.earth_density");
    186192        planetarea=4*M_PI*planetradius*planetradius;
    187193        parameters->AddObject(new DoubleParam(SolidearthPlanetAreaEnum,planetarea));
    … …  
    245251
    246252        parameters->FindParam(&grdmodel,GrdModelEnum);
    247         if(grdmodel==ElasticEnum){
     253        parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
     254        if(grdmodel==ElasticEnum && isgrd){
    248255                /*Deal with elasticity {{{*/
    249256                iomodel->FetchData(&selfattraction,"md.solidearth.settings.selfattraction");
    … …  
    259266                }
    260267
     268                //default values
     269                nt=1;
     270                ntimesteps=1;
    261271                /*love numbers: */
    262272                if(viscous || elastic){
    … …  
    273283
    274284                        parameters->AddObject(new DoubleParam(SolidearthSettingsTimeAccEnum,timeacc));
    275                         parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,ndeg,precomputednt));
    276                         parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
    277                         parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
     285                        //parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,ndeg,precomputednt));
     286                        //parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,ndeg,precomputednt));
     287                        //parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,ndeg,precomputednt));
    278288
    279289                        if (rotation){
    … …  
    286296                                parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionColinearEnum,love_pmtf_colinear,1,precomputednt));
    287297                                parameters->AddObject(new DoubleMatParam(LovePolarMotionTransferFunctionOrthogonalEnum,love_pmtf_ortho,1,precomputednt));
    288                                 parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
    289                                 parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
    290                                 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
     298                                //parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,ndeg,precomputednt));
     299                                //parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,ndeg,precomputednt));
     300                                //parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,ndeg,precomputednt));
    291301                        }
    292302
    293303                        parameters->AddObject(new DoubleMatParam(LoveTimeFreqEnum,love_timefreq,precomputednt,1));
    294304                        parameters->AddObject(new BoolParam(LoveIsTimeEnum,love_istime));
    295 
    296                         /*Free allocations:*/
    297                         xDelete<IssmDouble>(love_timefreq);
    298305
    299306                        // AD performance is sensitive to calls to ensurecontiguous.
    … …  
    320327                                }
    321328                        }
    322                         else {
    323                                 ntimesteps=1;
    324                                 nt=1;
    325                         }
    326 
    327329#ifdef _HAVE_AD_
    328                         G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
    329330                        U_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
    330331                        if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
    331332#else
    332                         G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
    333333                        U_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
    334334                        if(horiz)H_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
    … …  
    336336                }
    337337                if(selfattraction){
    338 #ifdef _HAVE_AD_
    339                         G_gravi=xNew<IssmDouble>(M,"t");
    340 #else
    341                         G_gravi=xNew<IssmDouble>(M);
    342 #endif
    343                 }
    344 
    345                 if(rotation) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
    346                 if(selfattraction){
    347 
    348338                        /*compute combined legendre + love number (elastic green function):*/
    349339                        m=DetermineLocalSize(M,IssmComm::GetComm());
    350340                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
     341#ifdef _HAVE_AD_
     342                        G_gravi=xNew<IssmDouble>(M,"t");
     343                        G_gravi_local=xNew<IssmDouble>(m,"t");
     344                        G_viscoelastic=xNew<IssmDouble>(M*ntimesteps,"t");
     345                        G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
     346#else
     347                        G_gravi=xNew<IssmDouble>(M);
     348                        G_gravi_local=xNew<IssmDouble>(m);
     349                        G_viscoelastic=xNew<IssmDouble>(M*ntimesteps);
     350                        G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
     351#endif
    351352                }
    352353                if(viscous | elastic){
    353354#ifdef _HAVE_AD_
    354                         G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
    355355                        U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
    356356                        if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps,"t");
    357357#else
    358                         G_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
    359358                        U_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
    360359                        if(horiz)H_viscoelastic_local=xNew<IssmDouble>(m*ntimesteps);
    361360#endif
    362361                }
    363                 if(selfattraction){
    364 #ifdef _HAVE_AD_
    365                         G_gravi_local=xNew<IssmDouble>(m,"t");
    366 #else
    367                         G_gravi_local=xNew<IssmDouble>(m);
    368 #endif
    369                 }
    370 
     362
     363                if(rotation) parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
     364                constant=3/rho_earth/planetarea;
    371365                if(selfattraction){
    372366                        for(int i=lower_row;i<upper_row;i++){
    373367                                IssmDouble alpha,x;
    374368                                alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
    375                                 G_gravi_local[i-lower_row]= .5/sin(alpha/2.0);
    376                         }
    377                 }
    378                 if(viscous | elastic){
    379                         for(int i=lower_row;i<upper_row;i++){
    380                                 IssmDouble alpha,x;
    381                                 alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
    382 
    383                                 for(int t=0;t<ntimesteps;t++){
    384                                         G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
    385                                         U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
    386                                         if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t]= 0;
    387                                 }
    388 
    389                                 IssmDouble Pn = 0.;
    390                                 IssmDouble Pn1 = 0.;
    391                                 IssmDouble Pn2 = 0.;
    392                                 IssmDouble Pn_p = 0.;
    393                                 IssmDouble Pn_p1 = 0.;
    394                                 IssmDouble Pn_p2 = 0.;
    395 
    396                                 for (int n=0;n<ndeg;n++) {
    397 
    398                                         /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    399                                         if(n==0){
    400                                                 Pn=1;
    401                                                 Pn_p=0;
     369                                G_gravi_local[i-lower_row]= constant*.5/sin(alpha/2.0);
     370                        }
     371                        if(viscous | elastic){
     372                                for(int i=lower_row;i<upper_row;i++){
     373                                        IssmDouble alpha,x;
     374                                        alpha= reCast<IssmDouble>(i)*degacc * M_PI / 180.0;
     375
     376                                        for(int t=0;t<ntimesteps;t++){
     377                                                G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (1.0+love_k[(ndeg-1)*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
     378                                                U_viscoelastic_local[(i-lower_row)*ntimesteps+t]= (love_h[(ndeg-1)*precomputednt+t])*G_gravi_local[i-lower_row];
     379                                                if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t]= 0;
    402380                                        }
    403                                         else if(n==1){
    404                                                 Pn = cos(alpha);
    405                                                 Pn_p = 1;
     381
     382                                        IssmDouble Pn = 0.;
     383                                        IssmDouble Pn1 = 0.;
     384                                        IssmDouble Pn2 = 0.;
     385                                        IssmDouble Pn_p = 0.;
     386                                        IssmDouble Pn_p1 = 0.;
     387                                        IssmDouble Pn_p2 = 0.;
     388
     389                                        for (int n=0;n<ndeg;n++) {
     390
     391                                                /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
     392                                                if(n==0){
     393                                                        Pn=1;
     394                                                        Pn_p=0;
     395                                                }
     396                                                else if(n==1){
     397                                                        Pn = cos(alpha);
     398                                                        Pn_p = 1;
     399                                                }
     400                                                else{
     401                                                        Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
     402                                                        Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     403                                                }
     404                                                Pn2=Pn1; Pn1=Pn;
     405                                                Pn_p2=Pn_p1; Pn_p1=Pn_p;
     406
     407                                                for(int t=0;t<ntimesteps;t++){
     408                                                        IssmDouble deltalove_G;
     409                                                        IssmDouble deltalove_U;
     410
     411                                                        deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
     412                                                        deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
     413
     414                                                        G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*deltalove_G*Pn;                            // gravitational potential
     415                                                        U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*deltalove_U*Pn;                            // vertical (up) displacement
     416                                                        if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += constant*sin(alpha)*love_l[n*precomputednt+t]*Pn_p;                // horizontal displacements
     417                                                }
    406418                                        }
    407                                         else{
    408                                                 Pn = ( (2*n-1)*cos(alpha)*Pn1 - (n-1)*Pn2 ) /n;
    409                                                 Pn_p = ( (2*n-1)*(Pn1+cos(alpha)*Pn_p1) - (n-1)*Pn_p2 ) /n;
     419                                }
     420                        }
     421                        else { //just copy G_gravi into G_viscoelastic
     422                                for(int i=lower_row;i<upper_row;i++){
     423                                        for(int t=0;t<ntimesteps;t++){
     424                                                G_viscoelastic_local[(i-lower_row)*ntimesteps+t]= G_gravi_local[i-lower_row];
    410425                                        }
    411                                         Pn2=Pn1; Pn1=Pn;
    412                                         Pn_p2=Pn_p1; Pn_p1=Pn_p;
    413 
    414                                         for(int t=0;t<ntimesteps;t++){
    415                                                 IssmDouble deltalove_G;
    416                                                 IssmDouble deltalove_U;
    417 
    418                                                 deltalove_G = (love_k[n*precomputednt+t]-love_k[(ndeg-1)*precomputednt+t]-love_h[n*precomputednt+t]+love_h[(ndeg-1)*precomputednt+t]);
    419                                                 deltalove_U = (love_h[n*precomputednt+t]-love_h[(ndeg-1)*precomputednt+t]);
    420 
    421                                                 G_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_G*Pn;                             // gravitational potential
    422                                                 U_viscoelastic_local[(i-lower_row)*ntimesteps+t] += deltalove_U*Pn;                             // vertical (up) displacement
    423                                                 if(horiz)H_viscoelastic_local[(i-lower_row)*ntimesteps+t] += sin(alpha)*love_l[n*precomputednt+t]*Pn_p;         // horizontal displacements
    424                                         }
    425                                 }
    426                         }
    427                 }
    428                 if(selfattraction){
    429 
     426                                }
     427                        }
    430428                        /*merge G_viscoelastic_local into G_viscoelastic; U_viscoelastic_local into U_viscoelastic; H_viscoelastic_local to H_viscoelastic:{{{*/
    431                         int* recvcounts=xNew<int>(IssmComm::GetSize());
    432                         int* displs=xNew<int>(IssmComm::GetSize());
     429                        recvcounts=xNew<int>(IssmComm::GetSize());
     430                        displs=xNew<int>(IssmComm::GetSize());
    433431                        int  rc;
    434432                        int  offset;
    … …  
    436434                        //deal with selfattraction first:
    437435                        ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
    438 
    439436                        /*displs: */
    440437                        ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
    441 
    442438                        /*All gather:*/
    443439                        ISSM_MPI_Allgatherv(G_gravi_local, m, ISSM_MPI_DOUBLE, G_gravi, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    444440
     441                        rc=m*ntimesteps;
     442                        offset=lower_row*ntimesteps;
     443                        ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
     444                        ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
     445                        ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    445446                        if(elastic){
    446                                 rc=m*ntimesteps;
    447                                 offset=lower_row*ntimesteps;
    448                                 ISSM_MPI_Allgather(&rc,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm());
    449                                 ISSM_MPI_Allgather(&offset,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm());
    450 
    451                                 ISSM_MPI_Allgatherv(G_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, G_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    452447                                ISSM_MPI_Allgatherv(U_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, U_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    453448                                if(horiz)ISSM_MPI_Allgatherv(H_viscoelastic_local, m*ntimesteps, ISSM_MPI_DOUBLE, H_viscoelastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
    … …  
    460455                        /*Avoid singularity at 0: */
    461456                        G_gravi[0]=G_gravi[1];
     457                        for(int t=0;t<ntimesteps;t++){
     458                                G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
     459                        }
    462460                        if(elastic){
    463461                                for(int t=0;t<ntimesteps;t++){
    464                                         G_viscoelastic[t]=G_viscoelastic[ntimesteps+t];
    465462                                        U_viscoelastic[t]=U_viscoelastic[ntimesteps+t];
    466463                                        if(horiz)H_viscoelastic[t]=H_viscoelastic[ntimesteps+t];
    … …  
    475472                                G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
    476473                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
    477                                 if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
     474                                if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt,"t");
    478475                                if(rotation){
    479476                                        Pmtf_col_interpolated=xNew<IssmDouble>(nt,"t");
    … …  
    487484                                G_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
    488485                                U_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
    489                                 if(horiz)H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
     486                                if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M*nt);
    490487                                if(rotation){
    491488                                        Pmtf_col_interpolated=xNew<IssmDouble>(nt);
    … …  
    497494                                }
    498495#endif
    499 
    500496                                for(int t=0;t<nt;t++){
    501497                                        IssmDouble lincoeff;
    502498                                        IssmDouble viscoelastic_time=t*timeacc;
    503499                                        int        timeindex2=-1;
    504 
    505500                                        /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
    506501                                        if(t!=0){
    … …  
    521516
    522517                                        for(int index=0;index<M;index++){
    523 
    524518                                                int timeindex=index*nt+t;
    525519                                                int timepreindex= index*ntimesteps+timeindex2;
    … …  
    539533                                        }
    540534                                }
    541 
    542                         }
    543                         else if(elastic){
     535                        }
     536                        else {
     537
    544538                                nt=1; //in elastic, or if we run only selfattraction, we need only one step
    545539#ifdef _HAVE_AD_
    546540                                G_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
    547                                 U_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
    548                                 if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
    549541#else
    550542                                G_viscoelastic_interpolated=xNew<IssmDouble>(M);
    551                                 U_viscoelastic_interpolated=xNew<IssmDouble>(M);
    552                                 if(horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M);
    553543#endif
    554 
    555 
    556544                                xMemCpy<IssmDouble>(G_viscoelastic_interpolated,G_viscoelastic,M);
    557                                 xMemCpy<IssmDouble>(U_viscoelastic_interpolated,U_viscoelastic,M);
    558                                 if(horiz) xMemCpy<IssmDouble>(H_viscoelastic_interpolated,H_viscoelastic,M);
    559 
    560                                 if(rotation){ //if this cpu handles degree 2
     545
     546                                if(elastic){
    561547#ifdef _HAVE_AD_
    562                                         Pmtf_col_interpolated=xNew<IssmDouble>(1,"t");
    563                                         Pmtf_ortho_interpolated=xNew<IssmDouble>(1,"t");
    564                                         Pmtf_z_interpolated=xNew<IssmDouble>(1,"t");
    565                                         Love_tk2_interpolated=xNew<IssmDouble>(1,"t");
    566                                         Love_th2_interpolated=xNew<IssmDouble>(1,"t");
    567                                         if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1,"t");
     548                                        U_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
     549                                        if (horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M,"t");
    568550#else
    569                                         Pmtf_col_interpolated=xNew<IssmDouble>(1);
    570                                         Pmtf_ortho_interpolated=xNew<IssmDouble>(1);
    571                                         Pmtf_z_interpolated=xNew<IssmDouble>(1);
    572                                         Love_tk2_interpolated=xNew<IssmDouble>(1);
    573                                         Love_th2_interpolated=xNew<IssmDouble>(1);
    574                                         if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1);
     551                                        U_viscoelastic_interpolated=xNew<IssmDouble>(M);
     552                                        if (horiz) H_viscoelastic_interpolated=xNew<IssmDouble>(M);
    575553#endif
    576 
    577                                         Pmtf_col_interpolated[0]=love_pmtf_colinear[0];
    578                                         Pmtf_ortho_interpolated[0]=love_pmtf_ortho[0];
    579                                         Pmtf_z_interpolated[0]=1.0+love_k[2];
    580                                         Love_tk2_interpolated[0]=love_tk[2];
    581                                         Love_th2_interpolated[0]=love_th[2];
    582                                         if (horiz) Love_tl2_interpolated[0]=love_tl[2];
     554                                        xMemCpy<IssmDouble>(U_viscoelastic_interpolated,U_viscoelastic,M);
     555                                        if (horiz) xMemCpy<IssmDouble>(H_viscoelastic_interpolated,H_viscoelastic,M);
     556
     557                                        if(rotation){ //if this cpu handles degree 2
     558#ifdef _HAVE_AD_
     559                                                Pmtf_col_interpolated=xNew<IssmDouble>(1,"t");
     560                                                Pmtf_ortho_interpolated=xNew<IssmDouble>(1,"t");
     561                                                Pmtf_z_interpolated=xNew<IssmDouble>(1,"t");
     562                                                Love_tk2_interpolated=xNew<IssmDouble>(1,"t");
     563                                                Love_th2_interpolated=xNew<IssmDouble>(1,"t");
     564                                                if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1,"t");
     565#else
     566                                                Pmtf_col_interpolated=xNew<IssmDouble>(1);
     567                                                Pmtf_ortho_interpolated=xNew<IssmDouble>(1);
     568                                                Pmtf_z_interpolated=xNew<IssmDouble>(1);
     569                                                Love_tk2_interpolated=xNew<IssmDouble>(1);
     570                                                Love_th2_interpolated=xNew<IssmDouble>(1);
     571                                                if (horiz) Love_tl2_interpolated=xNew<IssmDouble>(1);
     572#endif
     573
     574                                                Pmtf_col_interpolated[0]=love_pmtf_colinear[0];
     575                                                Pmtf_ortho_interpolated[0]=love_pmtf_ortho[0];
     576                                                Pmtf_z_interpolated[0]=1.0+love_k[2];
     577                                                Love_tk2_interpolated[0]=love_tk[2];
     578                                                Love_th2_interpolated[0]=love_th[2];
     579                                                if (horiz) Love_tl2_interpolated[0]=love_tl[2];
     580                                        }
     581
    583582                                }
    584583                        }       
    585 
    586584                        /*Save our precomputed tables into parameters*/
    587585                        parameters->AddObject(new DoubleVecParam(SealevelchangeGSelfAttractionEnum,G_gravi,M));
     586                        parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
    588587                        if(viscous || elastic){
    589                                 parameters->AddObject(new DoubleVecParam(SealevelchangeGViscoElasticEnum,G_viscoelastic_interpolated,M*nt));
    590588                                parameters->AddObject(new DoubleVecParam(SealevelchangeUViscoElasticEnum,U_viscoelastic_interpolated,M*nt));
    591589                                if(horiz)parameters->AddObject(new DoubleVecParam(SealevelchangeHViscoElasticEnum,H_viscoelastic_interpolated,M*nt));
    … …  
    599597                                }
    600598                        }
    601 
    602599                        /*free resources: */
    603600                        xDelete<IssmDouble>(G_gravi);
    604601                        xDelete<IssmDouble>(G_gravi_local);
     602                        xDelete<IssmDouble>(G_viscoelastic);
     603                        xDelete<IssmDouble>(G_viscoelastic_local);
     604                        xDelete<IssmDouble>(G_viscoelastic_interpolated);
    605605                        if(elastic){
     606                                xDelete<IssmDouble>(love_timefreq);
    606607                                xDelete<IssmDouble>(love_h);
    607608                                xDelete<IssmDouble>(love_k);
    … …  
    610611                                xDelete<IssmDouble>(love_tk);
    611612                                xDelete<IssmDouble>(love_tl);
    612                                 xDelete<IssmDouble>(G_viscoelastic);
    613                                 xDelete<IssmDouble>(G_viscoelastic_local);
     613
    614614                                xDelete<IssmDouble>(G_viscoelastic_interpolated);
    615615                                xDelete<IssmDouble>(U_viscoelastic);
    616616                                xDelete<IssmDouble>(U_viscoelastic_interpolated);
    617617                                xDelete<IssmDouble>(U_viscoelastic_local);
     618                                xDelete<IssmDouble>(U_viscoelastic_interpolated);
    618619                                if(horiz){
    619620                                        xDelete<IssmDouble>(H_viscoelastic);
    620621                                        xDelete<IssmDouble>(H_viscoelastic_interpolated);
    621622                                        xDelete<IssmDouble>(H_viscoelastic_local);
     623                                        xDelete<IssmDouble>(H_viscoelastic_interpolated);
    622624                                }
    623625                                if(rotation){
    624626                                        xDelete<IssmDouble>(love_pmtf_colinear);
    625627                                        xDelete<IssmDouble>(love_pmtf_ortho);
    626                                        
    627                                         xDelete<IssmDouble>(Love_tk2_interpolated);
    628                                         xDelete<IssmDouble>(Love_th2_interpolated);
    629628                                        xDelete<IssmDouble>(Pmtf_col_interpolated);
    630629                                        xDelete<IssmDouble>(Pmtf_ortho_interpolated);
    631630                                        xDelete<IssmDouble>(Pmtf_z_interpolated);
    632                                         if(horiz){
    633                                                 xDelete<IssmDouble>(Love_tl2_interpolated);
    634                                         }
    635 
     631                                        xDelete<IssmDouble>(Love_tk2_interpolated);
     632                                        xDelete<IssmDouble>(Love_th2_interpolated);
     633                                        if (horiz) xDelete<IssmDouble>(Love_tl2_interpolated);
    636634                                }
    637635                        }

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

    r27102 r27131  
    2222#include "../Inputs/ControlInput.h"
    2323#include "../Inputs/ArrayInput.h"
     24#include "../Inputs/IntArrayInput.h"
    2425/*}}}*/
    2526#define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/

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

    r27031 r27131  
    404404                virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
    405405                virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
    406                 virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
     406                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0;
    407407                virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    408408                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;

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

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

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

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

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

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

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

    r27123 r27131  
    62866286}
    62876287/*}}}*/
    6288 void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
     6288void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/
    62896289
    62906290        /*Declarations:{{{*/
    … …  
    63086308
    63096309        #ifdef _HAVE_RESTRICT_
    6310         IssmDouble* __restrict__ G=NULL;
    6311         IssmDouble* __restrict__ GU=NULL;
    6312         IssmDouble* __restrict__ GN=NULL;
    6313         IssmDouble* __restrict__ GE=NULL;
    6314         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    63156310        IssmDouble* __restrict__ G_gravi_precomputed=NULL;
    6316         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
    6317         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    63186311        #else
    6319         IssmDouble* G=NULL;
    6320         IssmDouble* GU=NULL;
    6321         IssmDouble* GN=NULL;
    6322         IssmDouble* GE=NULL;
    6323         IssmDouble* G_viscoelastic_precomputed=NULL;
    63246312        IssmDouble* G_gravi_precomputed=NULL;
    6325         IssmDouble* U_viscoelastic_precomputed=NULL;
    6326         IssmDouble* H_viscoelastic_precomputed=NULL;
    63276313        #endif
    63286314
    … …  
    63306316        int index;
    63316317        int M;
     6318        IssmDouble degacc;
    63326319        IssmDouble doubleindex,lincoef;
    63336320
    … …  
    63466333        /*Rotational:*/
    63476334        #ifdef _HAVE_RESTRICT_
     6335        IssmDouble* __restrict__ Grot=NULL;
     6336        IssmDouble* __restrict__ GUrot=NULL;
     6337        IssmDouble* __restrict__ GNrot=NULL;
     6338        IssmDouble* __restrict__ GErot=NULL;
    63486339        IssmDouble* __restrict__ tide_love_h  = NULL;
    63496340        IssmDouble* __restrict__ tide_love_k  = NULL;
    … …  
    63526343        IssmDouble* __restrict__ LoveRotU     = NULL;
    63536344        IssmDouble* __restrict__ LoveRothoriz = NULL;
    6354         IssmDouble* __restrict__  Grot        = NULL;
    6355         IssmDouble* __restrict__ GUrot        = NULL;
    6356         IssmDouble* __restrict__ GNrot        = NULL;
    6357         IssmDouble* __restrict__ GErot        = NULL;
     6345        int* __restrict__ AplhaIndex   = NULL;
     6346        int* __restrict__ AzimuthIndex = NULL;
    63586347        #else
     6348        IssmDouble* Grot=NULL;
     6349        IssmDouble* GUrot=NULL;
     6350        IssmDouble* GNrot=NULL;
     6351        IssmDouble* GErot=NULL;
    63596352        IssmDouble* tide_love_h  = NULL;
    63606353        IssmDouble* tide_love_k  = NULL;
    … …  
    63636356        IssmDouble* LoveRotU     = NULL;
    63646357        IssmDouble* LoveRothoriz = NULL;
    6365         IssmDouble*  Grot        = NULL;
    6366         IssmDouble* GUrot        = NULL;
    6367         IssmDouble* GNrot        = NULL;
    6368         IssmDouble* GErot        = NULL;
     6358        int* AlphaIndex   = NULL;
     6359        int* AzimuthIndex = NULL;
    63696360        #endif
    63706361
    … …  
    64056396
    64066397        /*Recover precomputed green function kernels:{{{*/
    6407         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
    6408         parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
    6409 
     6398        parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
     6399        M=reCast<int,IssmDouble>(180.0/degacc+1.);
     6400
     6401        DoubleVecParam* parameter;
    64106402        if(computeelastic){
    6411                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    6412                 parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    6413 
    6414                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
    6415                 parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
    6416 
    6417                 if(horiz){
    6418                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    6419                         parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
    6420                 }
    6421 
    64226403                if(computerotation){
    64236404                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalH2Enum)); _assert_(parameter);
    … …  
    64556436        }
    64566437
    6457         constant=3/rho_earth/planetarea;
    6458 
    6459         G=xNewZeroInit<IssmDouble>(3*nel*nt);
    6460         if(computeelastic){
    6461                 GU=xNewZeroInit<IssmDouble>(3*nel*nt);
    6462                 if(horiz){
    6463                         GN=xNewZeroInit<IssmDouble>(3*nel*nt);
    6464                         GE=xNewZeroInit<IssmDouble>(3*nel*nt);
    6465                 }
    6466         }
    6467 
    6468         for(int e=0;e<nel;e++){
    6469                 ratioe=constant*areae[e];
    6470                 for (int i=0;i<3;i++){
    6471                         IssmDouble alpha;
    6472                         IssmDouble delPhi,delLambda;
    6473                         /*recovers info for this element and vertex:*/
    6474                         IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
    6475                         IssmDouble longe= atan2(yye[e],xxe[e]);
    6476 
    6477                         lati=latitude[i];
    6478                         longi=longitude[i];
    6479 
    6480                         /*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
    6481                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    6482                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    6483                         doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    6484                         index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    6485                         _assert_(index>=0 && index<M);
    6486 
    6487                         lincoef=doubleindex-index; //where between index and index+1 is alpha
    6488                         if (index==M-1){ //avoids out of bound case
    6489                                 index-=1;
    6490                                 lincoef=1;
    6491                         }
    6492 
    6493                         if(horiz){
    6494                                 /*Compute azimuths, both north and east components: */
    6495                                 x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    6496                                 if(lati==M_PI/2){
    6497                                         x=1e-12; y=1e-12;
     6438        AlphaIndex=xNewZeroInit<int>(3*nel);
     6439        if (horiz) AzimuthIndex=xNewZeroInit<int>(3*nel);
     6440        int intmax=pow(2,16)-1;
     6441
     6442
     6443        for (int i=0;i<3;i++){
     6444                if(lids[this->vertices[i]->lid]==this->lid){
     6445                        for(int e=0;e<nel;e++){
     6446                                IssmDouble alpha;
     6447                                IssmDouble delPhi,delLambda;
     6448                                /*recovers info for this element and vertex:*/
     6449                                IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
     6450                                IssmDouble longe= atan2(yye[e],xxe[e]);
     6451
     6452                                lati=latitude[i];
     6453                                longi=longitude[i];
     6454
     6455                                /*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
     6456                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6457                                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     6458                                doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6459                                index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6460                                if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6461                                _assert_(index>=0 && index<M);
     6462
     6463                                if(horiz){
     6464                                        /*Compute azimuths*/
     6465                                        dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
     6466                                        dy=sin(longe-longi)*cos(late);
     6467                                        //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
     6468                                        AzimuthIndex[i*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    64986469                                }
    6499                                 if(lati==-M_PI/2){
    6500                                         x=1e-12; y=1e-12;
    6501                                 }
    6502                                 dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
    6503                                 N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6504                                 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6505                         }
    6506 
    6507                         for(int t=0;t<nt;t++){
    6508                                 int timeindex=i*nel*nt+e*nt+t;
    6509 
    6510                                 /*Rigid earth gravitational perturbation: */
    6511                                 G[timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef)
    6512                                              +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
    6513 
    6514                                 /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
    6515                                 if(computeelastic){
    6516                                         G[timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6517                                                      +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
    6518                                 }
    6519                                 G[timeindex]=G[timeindex]*ratioe;
    6520 
    6521                                 /*Elastic components:*/
    6522                                 if(computeelastic){
    6523                                         GU[timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6524                                                                  +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6525                                         if(horiz){
    6526                                                 GN[timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6527                                                                                +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6528                                                 GE[timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6529                                                                                +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6530                                         }
    6531                                 }
    6532                         } //for(int t=0;t<nt;t++)
     6470                                AlphaIndex[i*nel+e]=index;
     6471                        }                       
    65336472                } //for (int i=0;i<3;i++)
    65346473        } //for(int e=0;e<nel;e++)
    65356474
    65366475        /*Add in inputs:*/
    6537         this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3*nt);
    6538         if(computeelastic){
    6539                 this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3*nt);
    6540                 if(horiz){
    6541                         this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3*nt);
    6542                         this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3*nt);
    6543                 }
    6544         }
     6476        this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3);
     6477        if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3);
     6478
    65456479        /*}}}*/
    65466480        /*Compute rotation kernel:{{{*/
    … …  
    66606594        /*Free allocations:{{{*/
    66616595        #ifdef _HAVE_RESTRICT_
    6662         delete G;
    6663         if(computeelastic){
    6664                 delete GU;
    6665                 if(horiz){
    6666                         delete GN;
    6667                         delete GE;
    6668                 }
    6669                 if(computerotation){
    6670                         delete Grot;
    6671                         delete GUrot;
    6672                         if (horiz){
    6673                                 delete GNrot;
    6674                                 delete GErot;
    6675                         }
    6676                 }
    6677         }
     6596        delete AlphaIndex;
     6597        if(horiz) AzimuthIndex;
     6598
     6599        if(computerotation){
     6600                delete Grot;
     6601                delete GUrot;
     6602                if (horiz){
     6603                        delete GNrot;
     6604                        delete GErot;
     6605                }
     6606        }
     6607
    66786608        #else
    6679         xDelete(G);
    6680         if(computeelastic){
    6681                 xDelete(GU);
    6682                 if(horiz){
    6683                         xDelete(GN);
    6684                         xDelete(GE);
    6685                 }
    6686                 if(computerotation){
    6687                         xDelete(Grot);
    6688                         xDelete(GUrot);
    6689                         if (horiz){
    6690                                 xDelete(GNrot);
    6691                                 xDelete(GErot);
    6692                         }
     6609        xDelete<int>(AlphaIndex);
     6610        if(horiz){
     6611                xDelete<int>(AzimuthIndex);
     6612        }
     6613        if(computerotation){
     6614                xDelete<IssmDouble>(Grot);
     6615                xDelete<IssmDouble>(GUrot);
     6616                if (horiz){
     6617                        xDelete<IssmDouble>(GNrot);
     6618                        xDelete<IssmDouble>(GErot);
    66936619                }
    66946620        }
    66956621        #endif
     6622        /*}}}*/
     6623        return;
     6624
     6625}
     6626/*}}}*/
     6627
     6628void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
     6629
     6630        /*Declarations:{{{*/
     6631        int nel;
     6632        IssmDouble planetarea,planetradius;
     6633        IssmDouble constant,ratioe;
     6634        IssmDouble rho_earth;
     6635        IssmDouble lati,longi;
     6636        IssmDouble latitude[NUMVERTICES];
     6637        IssmDouble longitude[NUMVERTICES];
     6638        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
     6639        IssmDouble xyz_list[NUMVERTICES][3];
     6640
     6641        #ifdef _HAVE_RESTRICT_
     6642        int** __restrict__ AlphaIndex=NULL;
     6643        int** __restrict__ AzimIndex=NULL;
     6644
     6645        #else
     6646        int** AlphaIndex=NULL;
     6647        int** AzimIndex=NULL;
     6648        #endif
     6649
     6650        /*viscoelastic green function:*/
     6651        int index;
     6652        int M;
     6653        IssmDouble doubleindex,lincoef, degacc;
     6654
     6655        /*Computational flags:*/
     6656        bool computeselfattraction = false;
     6657        bool computeelastic = false;
     6658        bool computeviscous = false;
     6659        int  horiz;
     6660        int grd, grdmodel;
     6661
     6662        bool istime=true;
     6663        IssmDouble timeacc=0;
     6664        IssmDouble start_time,final_time;
     6665        int  nt,precomputednt;
     6666        int intmax=pow(2,16)-1;
     6667
     6668        /*}}}*/
     6669        /*Recover parameters:{{{ */
     6670        rho_earth=FindParam(MaterialsEarthDensityEnum);
     6671        this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
     6672        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     6673        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
     6674        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6675        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
     6676        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     6677        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6678        this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
     6679        this->parameters->FindParam(&grdmodel,GrdModelEnum);
     6680        /*}}}*/
     6681
     6682        /*early return:*/
     6683        if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
     6684        if(!computeselfattraction)return;
     6685
     6686        /*Recover precomputed green function kernels:{{{*/
     6687        parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
     6688        M=reCast<int,IssmDouble>(180.0/degacc+1.);
     6689
     6690        /*}}}*/
     6691        /*Compute lat long of all vertices in the element:{{{*/
     6692        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6693        for(int i=0;i<NUMVERTICES;i++){
     6694                latitude[i]= asin(xyz_list[i][2]/planetradius);
     6695                longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
     6696        }
     6697        /*}}}*/
     6698        /*Compute green functions:{{{ */
     6699
     6700        if(computeviscous){
     6701                this->parameters->FindParam(&istime,LoveIsTimeEnum);
     6702                if(!istime)_error_("Frequency love numbers not supported yet!");
     6703                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     6704                this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
     6705                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     6706                nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
     6707        }
     6708        else{
     6709                nt=1; //in elastic, or if we run only selfattraction, we need only one step
     6710        }
     6711        AlphaIndex=xNew<int*>(SLGEOM_NUMLOADS);
     6712        if(horiz) AzimIndex=xNew<int*>(SLGEOM_NUMLOADS);
     6713
     6714
     6715        //Allocate:
     6716        for(int l=0;l<SLGEOM_NUMLOADS;l++){
     6717                int nbar=slgeom->nbar[l];
     6718                AlphaIndex[l]=xNewZeroInit<int>(3*nbar);
     6719                if(horiz) AzimIndex[l]=xNewZeroInit<int>(3*nbar);
     6720
     6721
     6722                for (int i=0;i<3;i++){
     6723                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     6724                                for(int e=0;e<nbar;e++){
     6725                                        IssmDouble alpha;
     6726                                        IssmDouble delPhi,delLambda;
     6727                                        /*recover info for this element and vertex:*/
     6728                                        IssmDouble late= slgeom->latbarycentre[l][e];
     6729                                        IssmDouble longe= slgeom->longbarycentre[l][e];
     6730                                        late=late/180*M_PI;
     6731                                        longe=longe/180*M_PI;
     6732
     6733                                        lati=latitude[i];
     6734                                        longi=longitude[i];
     6735
     6736                                        if(horiz){
     6737                                                /*Compute azimuths*/
     6738                                                dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
     6739                                                dy=sin(longe-longi)*cos(late);
     6740                                                //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
     6741                                                AzimIndex[l][i*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6742                                        }
     6743
     6744                                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6745                                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6746                                        alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
     6747                                        doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6748                                        index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6749
     6750                                        if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6751                                        if (index==M-1){ //avoids out of bound case
     6752                                                index-=1;
     6753                                                lincoef=1;
     6754                                        }
     6755                                        AlphaIndex[l][i*nbar+e]=index;
     6756                                }
     6757                        }
     6758                }
     6759        }
     6760
     6761        /*Save all these arrayins for each element:*/
     6762        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6763                this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3);
     6764                if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3);
     6765        }
     6766        /*}}}*/
     6767        /*Free memory:{{{*/
     6768        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6769                xDelete<int>(AlphaIndex[l]);
     6770                if(horiz) xDelete<int>(AzimIndex[l]);
     6771        }
     6772        xDelete<int*>(AlphaIndex);
     6773        if(horiz) xDelete<int*>(AzimIndex);
     6774       
    66966775        /*}}}*/
    66976776        return;
    … …  
    69066985}
    69076986/*}}}*/
     6987void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
     6988
     6989        int nel;
     6990
     6991        /*Inputs:*/
     6992        IssmDouble I[NUMVERTICES];
     6993        IssmDouble W[NUMVERTICES];
     6994        IssmDouble BP[NUMVERTICES];
     6995        IssmDouble* areae=NULL;
     6996
     6997        /*output: */
     6998        IssmDouble bslcice=0;
     6999        IssmDouble bslchydro=0;
     7000        IssmDouble bslcbp=0;
     7001        IssmDouble BPavg=0;
     7002        IssmDouble Iavg=0;
     7003        IssmDouble Wavg=0;
     7004
     7005        /*ice properties: */
     7006        IssmDouble rho_ice,rho_water,rho_freshwater;
     7007
     7008        /*recover some parameters:*/
     7009        this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     7010        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     7011        this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
     7012        this->parameters->FindParam(&areae,&nel,AreaeEnum);
     7013
     7014        /*Retrieve inputs:*/
     7015        Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
     7016        Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
     7017        Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
     7018
     7019        for(int i=0;i<NUMVERTICES;i++){
     7020                Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]*slgeom->LoadArea[SLGEOM_ICE][this->lid];
     7021                Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]*slgeom->LoadArea[SLGEOM_WATER][this->lid];
     7022                BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]*slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     7023        }
     7024
     7025        /*convert from m^3 to kg:*/
     7026        Iavg*=rho_ice;
     7027        Wavg*=rho_freshwater;
     7028        BPavg*=rho_water;
     7029
     7030        #ifdef _ISSM_DEBUG_
     7031        this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
     7032        this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
     7033        this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
     7034        #endif
     7035
     7036        /*Compute barystatic component in kg:*/
     7037        // Note: Iavg, etc, already include partial area factor phi for subelement loading
     7038        bslcice =   -Iavg;
     7039        bslchydro = -Wavg;
     7040        bslcbp =    -BPavg;
     7041
     7042        _assert_(!xIsNan<IssmDouble>(bslcice));
     7043        _assert_(!xIsNan<IssmDouble>(bslchydro));
     7044        _assert_(!xIsNan<IssmDouble>(bslcbp));
     7045
     7046        /*Plug values into subelement load vector:*/
     7047        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
     7048                int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
     7049                loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
     7050                Iavg=0; //avoid double counting centroid loads and subelement loads!
     7051        }
     7052        if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
     7053                int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
     7054                loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
     7055                Wavg=0;
     7056        }
     7057        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7058                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7059                loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL);
     7060                BPavg=0;
     7061        }
     7062        /*Plug remaining values into centroid load vector:*/
     7063        loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
     7064
     7065        /*Keep track of barystatic contributions:*/
     7066        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
     7067
     7068        /*Free ressources*/
     7069        xDelete<IssmDouble>(areae);
     7070
     7071}/*}}}*/
    69087072void       Tria::SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){ /*{{{*/
    69097073
    … …  
    70067170}
    70077171/*}}}*/
    7008 void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
    7009 
    7010         /*Declarations:{{{*/
    7011         int nel;
    7012         IssmDouble planetarea,planetradius;
    7013         IssmDouble constant,ratioe;
    7014         IssmDouble rho_earth;
    7015         IssmDouble lati,longi;
    7016         IssmDouble latitude[NUMVERTICES];
    7017         IssmDouble longitude[NUMVERTICES];
    7018         IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    7019         IssmDouble xyz_list[NUMVERTICES][3];
    7020 
    7021         #ifdef _HAVE_RESTRICT_
    7022         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    7023         IssmDouble* __restrict__ G_gravi_precomputed=NULL;
    7024         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
    7025         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    7026         IssmDouble** __restrict__ Gsubel=NULL;
    7027         IssmDouble** __restrict__ GUsubel=NULL;
    7028         IssmDouble** __restrict__ GNsubel=NULL;
    7029         IssmDouble** __restrict__ GEsubel=NULL;
    7030 
    7031         #else
    7032         IssmDouble* G_viscoelastic_precomputed=NULL;
    7033         IssmDouble* G_gravi_precomputed=NULL;
    7034         IssmDouble* U_viscoelastic_precomputed=NULL;
    7035         IssmDouble* H_viscoelastic_precomputed=NULL;
    7036         IssmDouble** Gsubel=NULL;
    7037         IssmDouble** GUsubel=NULL;
    7038         IssmDouble** GNsubel=NULL;
    7039         IssmDouble** GEsubel=NULL;
    7040         #endif
    7041 
    7042         /*viscoelastic green function:*/
    7043         int index;
    7044         int M;
    7045         IssmDouble doubleindex,lincoef;
    7046 
    7047         /*Computational flags:*/
    7048         bool computeselfattraction = false;
    7049         bool computeelastic = false;
    7050         bool computeviscous = false;
    7051         int  horiz;
    7052         int grd, grdmodel;
    7053 
    7054         bool istime=true;
    7055         IssmDouble timeacc=0;
    7056         IssmDouble start_time,final_time;
    7057         int  nt,precomputednt;
    7058 
    7059         /*}}}*/
    7060         /*Recover parameters:{{{ */
    7061         rho_earth=FindParam(MaterialsEarthDensityEnum);
    7062         this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
    7063         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    7064         this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    7065         this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    7066         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    7067         this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    7068         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    7069         this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    7070         this->parameters->FindParam(&grdmodel,GrdModelEnum);
    7071         /*}}}*/
    7072 
    7073         /*early return:*/
    7074         if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
    7075         if(!computeselfattraction)return;
    7076 
    7077         /*Recover precomputed green function kernels:{{{*/
    7078         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
    7079         parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
    7080 
    7081         if(computeelastic){
    7082                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    7083                 parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    7084 
    7085                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
    7086                 parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
    7087 
    7088                 if(horiz){
    7089                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    7090                         parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
    7091 
    7092                 }
    7093         }
    7094         /*}}}*/
    7095         /*Compute lat long of all vertices in the element:{{{*/
    7096         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7097         for(int i=0;i<NUMVERTICES;i++){
    7098                 latitude[i]= asin(xyz_list[i][2]/planetradius);
    7099                 longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
    7100         }
    7101         /*}}}*/
    7102         /*Compute green functions:{{{ */
    7103         constant=3/rho_earth/planetarea;
    7104 
    7105         if(computeviscous){
    7106                 this->parameters->FindParam(&istime,LoveIsTimeEnum);
    7107                 if(!istime)_error_("Frequency love numbers not supported yet!");
    7108                 this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
    7109                 this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
    7110                 this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    7111                 nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
    7112         }
    7113         else{
    7114                 nt=1; //in elastic, or if we run only selfattraction, we need only one step
    7115         }
    7116         Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7117         if(computeelastic){
    7118                 GUsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7119                 if(horiz){
    7120                         GNsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7121                         GEsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7122                 }
    7123         }
    7124 
    7125         //Allocate:
    7126         for(int l=0;l<SLGEOM_NUMLOADS;l++){
    7127                 int nbar=slgeom->nbar[l];
    7128                 Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7129                 if(computeelastic){
    7130                         GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7131                         if(horiz){
    7132                                 GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7133                                 GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7134                         }
    7135                 }
    7136 
    7137                 for(int e=0;e<nbar;e++){
    7138                         ratioe=constant*slgeom->area_subel[l][e];
    7139                         for (int i=0;i<3;i++){
    7140                                 IssmDouble alpha;
    7141                                 IssmDouble delPhi,delLambda;
    7142                                 /*recover info for this element and vertex:*/
    7143                                 IssmDouble late= slgeom->latbarycentre[l][e];
    7144                                 IssmDouble longe= slgeom->longbarycentre[l][e];
    7145                                 late=late/180*M_PI;
    7146                                 longe=longe/180*M_PI;
    7147 
    7148                                 lati=latitude[i];
    7149                                 longi=longitude[i];
    7150 
    7151                                 if(horiz){
    7152                                         /*Compute azimuths, both north and east components: */
    7153                                         x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    7154                                         if(lati==90){
    7155                                                 x=1e-12; y=1e-12;
    7156                                         }
    7157                                         if(lati==-90){
    7158                                                 x=1e-12; y=1e-12;
    7159                                         }
    7160                                         IssmDouble xbar=planetradius*cos(late)*cos(longe);
    7161                                         IssmDouble ybar=planetradius*cos(late)*sin(longe);
    7162                                         IssmDouble zbar=planetradius*sin(late);
    7163 
    7164                                         dx = xbar-x; dy = ybar-y; dz = zbar-z;
    7165                                         N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    7166                                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    7167                                 }
    7168 
    7169                                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    7170                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    7171                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
    7172                                 doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    7173                                 index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    7174 
    7175                                 lincoef=doubleindex-index; //where between index and index+1 is alpha
    7176                                 if (index==M-1){ //avoids out of bound case
    7177                                         index-=1;
    7178                                         lincoef=1;
    7179                                 }
    7180 
    7181                                 for(int t=0;t<nt;t++){
    7182                                         int timeindex=i*nbar*nt+e*nt+t;
    7183 
    7184                                         /*Rigid earth gravitational perturbation: */
    7185                                         Gsubel[l][timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef)
    7186                                                              +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
    7187 
    7188                                         if(computeelastic){
    7189                                                 Gsubel[l][timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7190                                                                      +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
    7191                                         }
    7192                                         Gsubel[l][timeindex]*=ratioe;
    7193 
    7194                                         /*Elastic components:*/
    7195                                         if(computeelastic){
    7196                                                 GUsubel[l][timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7197                                                                                  +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7198                                                 if(horiz){
    7199                                                         GNsubel[l][timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7200                                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7201                                                         GEsubel[l][timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7202                                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7203                                                 }
    7204                                         }
    7205                                 }
    7206                         }
    7207                 }
    7208         }
    7209 
    7210         /*Save all these arrayins for each element:*/
    7211         for (int l=0;l<SLGEOM_NUMLOADS;l++){
    7212                 this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3*nt);
    7213                 if(computeelastic){
    7214                         this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3*nt);
    7215                         if(horiz){
    7216                                 this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3*nt);
    7217                                 this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3*nt);
    7218                         }
    7219                 }
    7220         }
    7221         /*}}}*/
    7222         /*Free memory:{{{*/
    7223         for (int l=0;l<SLGEOM_NUMLOADS;l++){
    7224                 xDelete<IssmDouble>(Gsubel[l]);
    7225                 if(computeelastic){
    7226                         xDelete<IssmDouble>(GUsubel[l]);
    7227                         if(horiz){
    7228                                 xDelete<IssmDouble>(GNsubel[l]);
    7229                                 xDelete<IssmDouble>(GEsubel[l]);
    7230                         }
    7231                 }
    7232         }
    7233         xDelete<IssmDouble*>(Gsubel);
    7234         if(computeelastic){
    7235                 xDelete<IssmDouble*>(GUsubel);
    7236                 if(horiz){
    7237                         xDelete<IssmDouble*>(GNsubel);
    7238                         xDelete<IssmDouble*>(GEsubel);
    7239                 }
    7240         }
    7241         /*}}}*/
    7242         return;
    7243 
    7244 }
    7245 /*}}}*/
    72467172void       Tria::SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){ /*{{{*/
    72477173
    … …  
    72527178        IssmDouble* viscousE=NULL;
    72537179        int         viscousnumsteps;
    7254         int         dummy;
     7180        int         size;
    72557181        bool        viscous=false;
    72567182        int         horiz=0;
    … …  
    72627188                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    72637189
    7264                 this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
    7265                 this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy);
     7190                this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&size);
     7191                this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&size);
    72667192                if(horiz){
    7267                         this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&dummy);
    7268                         this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy);
     7193                        this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&size);
     7194                        this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&size);
    72697195                }
    72707196
    … …  
    72777203                        }
    72787204                }
    7279 
    7280         }
    7281 
    7282 }
    7283 /*}}}*/
    7284 void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
    7285 
    7286         int nel;
    7287 
    7288         /*Inputs:*/
    7289         IssmDouble I[NUMVERTICES];
    7290         IssmDouble W[NUMVERTICES];
    7291         IssmDouble BP[NUMVERTICES];
    7292         IssmDouble* areae=NULL;
    7293 
    7294         /*output: */
    7295         IssmDouble bslcice=0;
    7296         IssmDouble bslchydro=0;
    7297         IssmDouble bslcbp=0;
    7298         IssmDouble BPavg=0;
    7299         IssmDouble Iavg=0;
    7300         IssmDouble Wavg=0;
    7301 
    7302         /*ice properties: */
    7303         IssmDouble rho_ice,rho_water,rho_freshwater;
    7304 
    7305         /*recover some parameters:*/
    7306         this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     7205        }
     7206}
     7207/*}}}*/
     7208void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
     7209
     7210        IssmDouble oceanaverage=0;
     7211        IssmDouble oceanarea=0;
     7212        IssmDouble rho_water;
     7213
    73077214        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    7308         this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
    7309         this->parameters->FindParam(&areae,&nel,AreaeEnum);
    7310 
    7311         /*Retrieve inputs:*/
    7312         Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
    7313         Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
    7314         Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
    7315 
     7215
     7216        /*retrieve ocean average and area:*/
    73167217        for(int i=0;i<NUMVERTICES;i++){
    7317                 Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid];
    7318                 Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid];
    7319                 BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    7320         }
    7321 
    7322         /*convert from m to kg/m^2:*/
    7323         Iavg*=rho_ice;
    7324         Wavg*=rho_freshwater;
    7325         BPavg*=rho_water;
    7326 
     7218                oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
     7219        }
    73277220        #ifdef _ISSM_DEBUG_
    7328         this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
    7329         this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
    7330         this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
     7221        this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    73317222        #endif
    7332 
    7333         /*Compute barystatic component in kg:*/
    7334         // Note: Iavg, etc, already include partial area factor phi for subelement loading
    7335         bslcice =   -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg;
    7336         bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
    7337         bslcbp =    -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg;
    7338 
    7339         _assert_(!xIsNan<IssmDouble>(bslcice));
    7340         _assert_(!xIsNan<IssmDouble>(bslchydro));
    7341         _assert_(!xIsNan<IssmDouble>(bslcbp));
    7342 
    7343         /*Plug values into subelement load vector:*/
    7344         if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
    7345                 int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
    7346                 loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
    7347                 Iavg=0; //avoid double counting centroid loads and subelement loads!
    7348         }
    7349         if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
    7350                 int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
    7351                 loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
    7352                 Wavg=0;
    7353         }
     7223        oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     7224
     7225        /*add ocean average in the global sealevelloads vector:*/
    73547226        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    73557227                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7356                 loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL);
    7357                 BPavg=0;
    7358         }
    7359         /*Plug remaining values into centroid load vector:*/
    7360         loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
    7361 
    7362         /*Keep track of barystatic contributions:*/
    7363         barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
    7364 
    7365         /*Free resources:*/
    7366         xDelete<IssmDouble>(areae);
    7367 
     7228                loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL);
     7229                loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
     7230        }
     7231        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL);
     7232
     7233        /*add ocean area into a global oceanareas vector:*/
     7234        if(!loads->sealevelloads){
     7235                oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     7236                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7237                        int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7238                        subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
     7239                }
     7240        }
    73687241}
    73697242/*}}}*/
    … …  
    73717244
    73727245        /*sal green function:*/
     7246        int* AlphaIndex=NULL;
     7247        int* AlphaIndexsub[SLGEOM_NUMLOADS];
    73737248        IssmDouble* G=NULL;
    73747249        IssmDouble* Grot=NULL;
    7375         IssmDouble* Gsub[SLGEOM_NUMLOADS];
     7250        DoubleVecParam* parameter;
    73767251        bool computefuture=false;
     7252        int spatial_component=0;
    73777253
    73787254        bool sal = false;
    … …  
    73907266
    73917267        if(sal){
    7392                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    7393                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
    7394                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    7395                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     7268                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     7269                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
     7270
     7271                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7272                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    73967273                if (rotation)   this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
    73977274
    7398                 this->SealevelchangeGxL(sealevelpercpu, G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    7399         }
    7400 
    7401         return;
    7402 } /*}}}*/
    7403 void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
    7404 
    7405         IssmDouble oceanaverage=0;
    7406         IssmDouble oceanarea=0;
    7407         IssmDouble rho_water;
    7408 
    7409         this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    7410 
    7411         /*retrieve ocean average and area:*/
    7412         for(int i=0;i<NUMVERTICES;i++){
    7413                 oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    7414         }
    7415         #ifdef _ISSM_DEBUG_
    7416         this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    7417         #endif
    7418         oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
    7419 
    7420         /*add ocean average in the global sealevelloads vector:*/
    7421         if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    7422                 int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7423                 loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
    7424                 loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
    7425         }
    7426         else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    7427 
    7428         /*add ocean area into a global oceanareas vector:*/
    7429         if(!loads->sealevelloads){
    7430                 oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    7431                 if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    7432                         int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7433                         subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
    7434                 }
     7275                this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    74357276        }
    74367277
    … …  
    74467287        int nel,nbar;
    74477288        bool sal = false;
     7289        int* AlphaIndex=NULL;
     7290        int* AzimIndex=NULL;
     7291        int* AlphaIndexsub[SLGEOM_NUMLOADS];
     7292        int* AzimIndexsub[SLGEOM_NUMLOADS];
     7293        int spatial_component=0;
    74487294        IssmDouble* G=NULL;
    74497295        IssmDouble* GU=NULL;
    7450         IssmDouble* GE=NULL;
    7451         IssmDouble* GN=NULL;
     7296        IssmDouble* GH=NULL;
    74527297        IssmDouble* Grot=NULL;
    74537298        IssmDouble* GUrot=NULL;
    74547299        IssmDouble* GNrot=NULL;
    74557300        IssmDouble* GErot=NULL;
    7456         IssmDouble* Gsub[SLGEOM_NUMLOADS];
    7457         IssmDouble* GUsub[SLGEOM_NUMLOADS];
    7458         IssmDouble* GNsub[SLGEOM_NUMLOADS];
    7459         IssmDouble* GEsub[SLGEOM_NUMLOADS];
     7301
     7302        DoubleVecParam* parameter;
    74607303        bool computefuture=false;
    74617304
    … …  
    74767319
    74777320        if(sal){
    7478 
    7479                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    7480                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
    7481                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    7482                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     7321                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7322                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
     7323
     7324                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     7325                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    74837326
    74847327                if(elastic){
    7485                         this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
    7486                         this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsub[SLGEOM_ICE],&size);
    7487                         this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsub[SLGEOM_WATER],&size);
    7488                         this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsub[SLGEOM_OCEAN],&size);
     7328                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
     7329                        parameter->GetParameterValueByPointer((IssmDouble**)&GU,NULL);
     7330
    74897331                        if(horiz){
    7490                                 this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
    7491                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsub[SLGEOM_ICE],&size);
    7492                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsub[SLGEOM_WATER],&size);
    7493                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsub[SLGEOM_OCEAN],&size);
    7494 
    7495                                 this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
    7496                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsub[SLGEOM_ICE],&size);
    7497                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsub[SLGEOM_WATER],&size);
    7498                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size);
     7332                                this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
     7333                                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
     7334
     7335                                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
     7336                                parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL);
    74997337                        }
    75007338                        if (rotation) {
    … …  
    75077345                        }
    75087346                }
    7509                 this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
     7347
     7348                this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
    75107349
    75117350                if(elastic){
    7512                         this->SealevelchangeGxL(&UGrd[0],GU, GUsub, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
     7351                        this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
    75137352                        if(horiz){
    7514                                 this->SealevelchangeGxL(&NGrd[0],GN, GNsub, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
    7515                                 this->SealevelchangeGxL(&EGrd[0],GE, GEsub, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
     7353                                this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
     7354                                this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    75167355                        }
    75177356                }
    … …  
    75387377
    75397378} /*}}}*/
    7540 void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
     7379void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
    75417380
    75427381        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    75437382
    75447383        IssmDouble* grdfield=NULL;
    7545         int i,e,l,t,nbar;
     7384        int i,e,l,t,a, index, nbar;
     7385        bool rotation=false;
     7386        IssmDouble* Centroid_loads=NULL;
     7387        IssmDouble* Centroid_loads_copy=NULL;
     7388        IssmDouble* Subelement_loads[SLGEOM_NUMLOADS];
     7389        IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS];
     7390        IssmDouble* horiz_projection=NULL;
     7391        IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
     7392        int nt=1; //important, ensures there is a defined value if computeviscous is false
     7393
     7394        //viscous
    75467395        bool computeviscous=false;
    7547         bool rotation=false;
    7548         IssmDouble* viscousfield=NULL;
    7549         int nt=1; //important, ensures there is a defined value if computeviscous is false
    75507396        int viscousindex=0; //important
    75517397        int viscousnumsteps=1; //important
     7398        IssmDouble* viscousfield=NULL;
     7399        IssmDouble* grdfieldinterp=NULL;
     7400        IssmDouble* viscoustimes=NULL;
     7401        IssmDouble  final_time;
     7402        IssmDouble  lincoeff;
     7403        IssmDouble  timeacc;
    75527404
    75537405        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    … …  
    75567408                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    75577409                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     7410                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
     7411                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     7412                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     7413                this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
    75587414                if(computefuture) {
    75597415                        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
    75607416                        if (nt>viscousnumsteps) nt=viscousnumsteps;
     7417                        grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps);
    75617418                }
    75627419                else nt=1;
    … …  
    75787435        }
    75797436
    7580 
    7581         if(loads->sealevelloads){ // general case: loads + sealevel loads
    7582                 for(i=0;i<NUMVERTICES;i++) {
    7583                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7584                         for (e=0;e<nel;e++){
    7585                                 for(t=0;t<nt;t++){
    7586                                         int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
    7587                                         grdfield[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
     7437        //Determine loads /*{{{*/
     7438        Centroid_loads=xNewZeroInit<IssmDouble>(nel);
     7439        for (e=0;e<nel;e++){
     7440                Centroid_loads[e]=loads->loads[e];
     7441        }
     7442        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7443                nbar=slgeom->nbar[l];
     7444                Subelement_loads[l]=xNewZeroInit<IssmDouble>(nbar);
     7445                for (e=0;e<nbar;e++){
     7446                        Subelement_loads[l][e]=(loads->subloads[l][e]);
     7447                }
     7448        }
     7449        if(loads->sealevelloads){
     7450                for (e=0;e<nel;e++){
     7451                        Centroid_loads[e]+=(loads->sealevelloads[e]);
     7452                }
     7453                nbar=slgeom->nbar[SLGEOM_OCEAN];
     7454                for (e=0;e<nbar;e++){
     7455                        Subelement_loads[SLGEOM_OCEAN][e]+=(loads->subsealevelloads[e]);
     7456                }
     7457        }
     7458
     7459        //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex
     7460        if (spatial_component!=0){
     7461                horiz_projection=xNewZeroInit<IssmDouble>(3*nel);
     7462                Centroid_loads_copy=xNewZeroInit<IssmDouble>(nel);
     7463                for (e=0;e<nel;e++){
     7464                        Centroid_loads_copy[e]=Centroid_loads[e];
     7465                }
     7466
     7467                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7468                        nbar=slgeom->nbar[l];
     7469                        Subelement_loads_copy[l]=xNewZeroInit<IssmDouble>(nbar);
     7470                        horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(3*nbar);
     7471                        for (e=0;e<nbar;e++){
     7472                                Subelement_loads_copy[l][e]=Subelement_loads[l][e];
     7473                        }
     7474                }
     7475        }
     7476        /*}}}*/
     7477
     7478        //Convolution
     7479        for(i=0;i<NUMVERTICES;i++) { /*{{{*/
     7480                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7481
     7482                if (spatial_component!=0){ //horizontals /*{{{*/
     7483                        //GxL needs to be projected on the right axis before summation into the grd field
     7484                        //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
     7485                        if (spatial_component==1){ //north
     7486                                for (e=0;e<nel;e++){
     7487                                        horiz_projection[i*nel+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
    75887488                                }
    7589                         }
     7489                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7490                                        nbar=slgeom->nbar[l];
     7491                                        for (e=0;e<nbar;e++){
     7492                                                horiz_projectionsub[l][i*nbar+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7493                                        }
     7494                                }
     7495                        }
     7496                        else if (spatial_component==2){ //east
     7497                                for (e=0;e<nel;e++){
     7498                                        horiz_projection[i*nel+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0);
     7499                                }
     7500                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7501                                        nbar=slgeom->nbar[l];
     7502                                        for (e=0;e<nbar;e++){
     7503                                                horiz_projectionsub[l][i*nbar+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7504                                        }
     7505                                }
     7506                        }
     7507                        for (e=0;e<nel;e++) Centroid_loads[e]=Centroid_loads_copy[e]*horiz_projection[i*nel+e];
    75907508                        for(l=0;l<SLGEOM_NUMLOADS;l++){
    75917509                                nbar=slgeom->nbar[l];
    75927510                                for (e=0;e<nbar;e++){
    7593                                         for(t=0;t<nt;t++){
    7594                                                 int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7595                                                 grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
    7596                                         }
     7511                                        Subelement_loads[l][e]=Subelement_loads_copy[l][e]*horiz_projectionsub[l][i*nbar+e];
    75977512                                }
    7598                                 if(l==SLGEOM_OCEAN){
    7599                                         for (e=0;e<nbar;e++){
    7600                                                 for(t=0;t<nt;t++){
    7601                                                         int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7602                                                         grdfield[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
    7603                                                 }
    7604                                         }
     7513                        }
     7514                } /*}}}*/
     7515
     7516                for (e=0;e<nel;e++){
     7517                        for(t=0;t<nt;t++){
     7518                                a=AlphaIndex[i*nel+e];
     7519                                grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Centroid_loads[e];
     7520                        }
     7521                }
     7522                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7523                        nbar=slgeom->nbar[l];
     7524                        for (e=0;e<nbar;e++){
     7525                                for(t=0;t<nt;t++){
     7526                                        a=AlphaIndexsub[l][i*nbar+e];
     7527                                        grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Subelement_loads[l][e];
    76057528                                }
    76067529                        }
    76077530                }
    7608         }
    7609         else{  //this is the initial convolution where only loads are provided
    7610                 for(i=0;i<NUMVERTICES;i++) {
    7611                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7612                         for (e=0;e<nel;e++){
    7613                                 for(t=0;t<nt;t++){
    7614                                         int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
    7615                                         grdfield[i*nt+t]+=G[index]*(loads->loads[e]);
    7616                                 }
    7617                         }
    7618                         for(l=0;l<SLGEOM_NUMLOADS;l++){
    7619                                 nbar=slgeom->nbar[l];
    7620                                 for (e=0;e<nbar;e++){
    7621                                         for(t=0;t<nt;t++){
    7622                                                 int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7623                                                 grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
    7624                                         }
    7625                                 }
    7626                         }
    7627                 }
    7628         }
     7531        } /*}}}*/
     7532
     7533
    76297534
    76307535        if(computeviscous){ /*{{{*/
    … …  
    76347539                // 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
    76357540
    7636                 IssmDouble* grdfieldinterp=NULL;
    7637                 IssmDouble* viscoustimes=NULL;
    7638                 IssmDouble  final_time;
    7639                 IssmDouble  lincoeff;
    7640                 IssmDouble  timeacc;
    7641 
    7642                 this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
    7643                 this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    7644                 this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
    7645                 this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
    76467541                /* Map new grdfield generated by present-day loads onto viscous time vector*/
    76477542                if(computefuture){
    7648                         grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps);
    76497543                        //viscousindex time and first time step of grdfield coincide, so just copy that value
    76507544                        for(int i=0;i<NUMVERTICES;i++){
    … …  
    76587552                                        for(int i=0;i<NUMVERTICES;i++){
    76597553                                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7660                                                 grdfieldinterp[i*viscousnumsteps+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
     7554                                                grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]
     7555                                                                                         +lincoeff*grdfield[i*nt+(t-viscousindex)];
    76617556                                        }
    76627557                                }
    … …  
    76727567                /*update viscous stack with future deformation from present load: */
    76737568                if(computefuture){
    7674                         for(int t=viscousnumsteps-1;t>=viscousindex;t--){
     7569                        for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end
    76757570                                for(int i=0;i<NUMVERTICES;i++){
    76767571                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    76777572                                        //offset viscousfield to remove all deformation that has already been added
    7678                                         viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]-grdfieldinterp[i*viscousnumsteps+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
     7573                                        viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]
     7574                                                                          -grdfieldinterp[i*viscousnumsteps+viscousindex]
     7575                                                                          -viscousfield[i*viscousnumsteps+viscousindex];
    76797576                                }
    76807577                        }
    … …  
    76857582                        xDelete<IssmDouble>(grdfieldinterp);
    76867583                }
    7687 
    7688                 /*Free allocatoins:*/
    7689                 xDelete<IssmDouble>(viscoustimes);
    76907584        }
    76917585        /*}}}*/
    … …  
    77027596                for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
    77037597        }
    7704         /*Free resources:*/
     7598        //free resources
    77057599        xDelete<IssmDouble>(grdfield);
     7600        xDelete<IssmDouble>(Centroid_loads);
     7601        for(l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Subelement_loads[l]);
     7602        if (spatial_component!=0){
     7603                xDelete<IssmDouble>(horiz_projection);
     7604                xDelete<IssmDouble>(Centroid_loads_copy);
     7605                for(l=0;l<SLGEOM_NUMLOADS;l++) {
     7606                        xDelete<IssmDouble>(Subelement_loads_copy[l]);
     7607                        xDelete<IssmDouble>(horiz_projectionsub[l]);
     7608                }
     7609        }
     7610        if (computeviscous){
     7611                xDelete<IssmDouble>(viscoustimes);
     7612                if (computefuture){
     7613                        xDelete<IssmDouble>(grdfieldinterp);
     7614                }
     7615        }
    77067616
    77077617} /*}}}*/
    … …  
    77097619void       Tria::SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    77107620
     7621        offset*=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; //kg.m^-2 to kg
    77117622        if(slgeom->isoceanin[this->lid]){
    77127623                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){

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

    r26850 r27131  
    171171                #ifdef _HAVE_SEALEVELCHANGE_
    172172                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    173                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae);
     173                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids);
    174174                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom);
    175175                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
    … …  
    244244                void           UpdateConstraintsExtrudeFromBase(void);
    245245                void           UpdateConstraintsExtrudeFromTop(void);
    246                 void           SealevelchangeGxL(IssmDouble* sealevel, IssmDouble* G, IssmDouble** Gsub, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector,SealevelGeometry* slgeom, int nel, bool percpu,int stackenum,bool computefuture);
     246                void           SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
    247247                /*}}}*/
    248248

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

    r27097 r27131  
    9090void GrdLoads::SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom){ /*{{{*/
    9191
    92         IssmDouble area,re, S;
     92        IssmDouble re, S;
    9393        int ylmindex, intj;
    9494        IssmDouble deg2coeff_local[5];
     95        //IssmDouble area;
    9596
    9697        femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);
    … …  
    106107                if(sealevelloads) S+=sealevelloads[element->Sid()];
    107108                if(S!=0){
    108                         element->Element::GetInputValue(&area,AreaEnum);
    109109
    110110                        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
    111111                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
    112                                 deg2coeff_local[c] += S*area/re/re*slgeom->Ylm[ylmindex];
     112                                deg2coeff_local[c] += S/re/re*slgeom->Ylm[ylmindex];
    113113                        }
    114114                }
    … …  
    121121                                if(i==SLGEOM_OCEAN && sealevelloads) S+=subsealevelloads[intj];
    122122                                if(S!=0){
    123                                         area=slgeom->area_subel[i][intj];
     123                                        //area=slgeom->area_subel[i][intj];
    124124                                        for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients
    125125                                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
    126                                                 deg2coeff_local[c] += S*area/re/re*slgeom->Ylm_subel[i][ylmindex];
     126                                                deg2coeff_local[c] += S/re/re*slgeom->Ylm_subel[i][ylmindex];
    127127                                        }
    128128                                }

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

    r26944 r27131  
    115115
    116116} /*}}}*/
    117 int SealevelGeometry::GEnum(int l){ /*{{{*/
     117int SealevelGeometry::AlphaIndexEnum(int l){ /*{{{*/
    118118
    119119        int output = -1;
    120120        switch(l){
    121                 case SLGEOM_OCEAN: output=SealevelchangeGsubelOceanEnum; break;
    122                 case SLGEOM_ICE:   output=SealevelchangeGsubelIceEnum;   break;
    123                 case SLGEOM_WATER: output=SealevelchangeGsubelHydroEnum; break;
     121                case SLGEOM_OCEAN: output=SealevelchangeAlphaIndexOceanEnum; break;
     122                case SLGEOM_ICE:   output=SealevelchangeAlphaIndexIceEnum;   break;
     123                case SLGEOM_WATER: output=SealevelchangeAlphaIndexHydroEnum; break;
    124124                default: _error_("not supported");
    125125        }
    … …  
    127127
    128128} /*}}}*/
    129 int SealevelGeometry::GUEnum(int l){ /*{{{*/
     129int SealevelGeometry::AzimuthIndexEnum(int l){ /*{{{*/
    130130
    131131        int output = -1;
    132132        switch(l){
    133                 case SLGEOM_OCEAN: output=SealevelchangeGUsubelOceanEnum; break;
    134                 case SLGEOM_ICE:   output=SealevelchangeGUsubelIceEnum;   break;
    135                 case SLGEOM_WATER: output=SealevelchangeGUsubelHydroEnum; break;
    136                 default: _error_("not supported");
    137         }
    138         return output;
    139 
    140 } /*}}}*/
    141 int SealevelGeometry::GNEnum(int l){ /*{{{*/
    142 
    143         int output = -1;
    144         switch(l){
    145                 case SLGEOM_OCEAN: output=SealevelchangeGNsubelOceanEnum; break;
    146                 case SLGEOM_ICE:   output=SealevelchangeGNsubelIceEnum;   break;
    147                 case SLGEOM_WATER: output=SealevelchangeGNsubelHydroEnum; break;
    148                 default: _error_("not supported");
    149         }
    150 
    151         return output;
    152 } /*}}}*/
    153 int SealevelGeometry::GEEnum(int l){ /*{{{*/
    154 
    155         int output = -1;
    156         switch(l){
    157                 case SLGEOM_OCEAN: output=SealevelchangeGEsubelOceanEnum; break;
    158                 case SLGEOM_ICE:   output=SealevelchangeGEsubelIceEnum;   break;
    159                 case SLGEOM_WATER: output=SealevelchangeGEsubelHydroEnum; break;
     133                case SLGEOM_OCEAN: output=SealevelchangeAzimuthIndexOceanEnum; break;
     134                case SLGEOM_ICE:   output=SealevelchangeAzimuthIndexIceEnum;   break;
     135                case SLGEOM_WATER: output=SealevelchangeAzimuthIndexHydroEnum; break;
    160136                default: _error_("not supported");
    161137        }

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

    r26800 r27131  
    4545                void InitializeMappingsAndBarycentres(void);
    4646                void Assemble(void);
    47                 int GEnum(int l);
    48                 int GUEnum(int l);
    49                 int GNEnum(int l);
    50                 int GEEnum(int l);
     47                int AlphaIndexEnum(int l);
     48                int AzimuthIndexEnum(int l);
    5149                void BuildSphericalHarmonics(void);
    5250};

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

    r27105 r27131  
    2929void ivins_deformation_core(FemModel* femmodel);
    3030IssmDouble* CombineLoads(IssmDouble* load,IssmDouble* subload,FemModel* femmodel, SealevelGeometry* slgeom,int loadtype,int nel);
     31void slc_geometry_cleanup(SealevelGeometry* slgeom, FemModel* femmodel);
    3132/*}}}*/
    3233
    … …  
    8384
    8485        /*Free resources:*/
    85         delete slgeom;
     86        slc_geometry_cleanup(slgeom, femmodel);
    8687}
    8788/*}}}*/
    … …  
    361362                //Conserve ocean mass:
    362363                oceanaverage=SealevelloadsOceanAverage(loads, oceanareas,subelementoceanareas, totaloceanarea);
    363                 ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea - oceanaverage,slgeom);
     364                ConserveOceanMass(femmodel,loads,barycontrib->Total()/totaloceanarea -oceanaverage,slgeom);
    364365
    365366                //broadcast sea level loads
    … …  
    443444        delete oceanareas;
    444445        xDelete<IssmDouble>(sealevelpercpu);
    445 
    446446}
    447447/*}}}*/
    … …  
    588588        IssmDouble* areae  = NULL;
    589589        int  nel;
     590        int* lids;
    590591        int  grdmodel=0;
    591592
    … …  
    603604        ElementCoordinatesx(&xxe,&yye,&zze,&areae,femmodel->elements);
    604605
     606
     607        /*Compute element ids, used to speed up computations in convolution phase:{{{*/
     608        lids=xNew<int>(femmodel->vertices->Size());
     609
     610        for(Object* & object : femmodel->elements->objects){
     611                Element*   element=xDynamicCast<Element*>(object);
     612                for(int i=0;i<3;i++){
     613                        lids[element->vertices[i]->lid]=element->lid;
     614                }
     615        }
     616
     617        /*}}}*/
     618
    605619        /*Run sealevel geometry routine in elements:*/
    606620        for(Object* & object : femmodel->elements->objects){
    607621                Element*   element=xDynamicCast<Element*>(object);
    608                 element->SealevelchangeGeometryInitial(xxe,yye,zze,areae);
     622                element->SealevelchangeGeometryInitial(xxe,yye,zze,areae,lids);
    609623        }
    610624
    … …  
    625639        xDelete<IssmDouble>(zze);
    626640        xDelete<IssmDouble>(areae);
     641        xDelete(lids);
    627642
    628643        return;
    … …  
    642657        int nel;
    643658        int  grdmodel=0;
     659        int isgrd=0;
    644660        SealevelGeometry* slgeom=NULL;
    645661
    646662        /*early return?:*/
    647663        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
    648         if(grdmodel==IvinsEnum) return NULL;
     664        femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
     665        if(grdmodel!=ElasticEnum || !isgrd) return NULL;
    649666
    650667        /*retrieve parameters:*/
    … …  
    701718
    702719}/*}}}*/
     720void slc_geometry_cleanup(SealevelGeometry* slgeom, FemModel* femmodel){  /*{{{*/
     721        int  grdmodel=0;
     722        int isgrd=0;
     723        int horiz=0;
     724
     725        /*early return?:*/
     726        femmodel->parameters->FindParam(&grdmodel,GrdModelEnum);
     727        femmodel->parameters->FindParam(&isgrd,SolidearthSettingsGRDEnum);
     728        femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     729        if(grdmodel!=ElasticEnum || !isgrd) return;
     730
     731        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     732                femmodel->inputs->DeleteInput(slgeom->AlphaIndexEnum(l));
     733                if (horiz) femmodel->inputs->DeleteInput(slgeom->AzimuthIndexEnum(l));
     734        }
     735
     736        delete slgeom;
     737} /*}}}*/
    703738
    704739/*subroutines:*/
    … …  
    756791        Vector<IssmDouble>* vsealevelloadsvolume=loads->vsealevelloads->Duplicate();
    757792        Vector<IssmDouble>* vsubsealevelloadsvolume=loads->vsubsealevelloads->Duplicate();
    758 
    759         vsealevelloadsvolume->PointwiseMult(loads->vsealevelloads,oceanareas);
    760         vsubsealevelloadsvolume->PointwiseMult(loads->vsubsealevelloads,suboceanareas);
    761793
    762794        vsealevelloadsvolume->Sum(&sealevelloadsaverage);
    … …  
    965997                /*free allocations:*/
    966998                xDelete<IssmDouble>(viscoustimes);
     999                if (rotation)   xDelete<IssmDouble>(viscouspolarmotion);
    9671000        }
    9681001

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

    r27122 r27131  
    289289syn keyword cConstant LoveInnerCoreBoundaryEnum
    290290syn keyword cConstant LoveComplexComputationEnum
    291 syn keyword cConstant LoveIntStepsPerLayerEnum
     291syn keyword cConstant LoveQuadPrecisionEnum
     292syn keyword cConstant LoveMinIntegrationStepsEnum
     293syn keyword cConstant LoveMaxIntegrationdrEnum
    292294syn keyword cConstant LoveKernelsEnum
    293295syn keyword cConstant LoveMu0Enum
    … …  
    301303syn keyword cConstant LoveUnderflowTolEnum
    302304syn keyword cConstant LovePostWidderThresholdEnum
     305syn keyword cConstant LoveDebugEnum
     306syn keyword cConstant LoveHypergeomNZEnum
     307syn keyword cConstant LoveHypergeomNAlphaEnum
    303308syn keyword cConstant MassFluxSegmentsEnum
    304309syn keyword cConstant MassFluxSegmentsPresentEnum
    … …  
    387392syn keyword cConstant SolidearthSettingsElasticEnum
    388393syn keyword cConstant SolidearthSettingsViscousEnum
     394syn keyword cConstant SolidearthSettingsSatelliteGraviEnum
     395syn keyword cConstant SolidearthSettingsDegreeAccuracyEnum
    389396syn keyword cConstant SealevelchangeGeometryDoneEnum
    390397syn keyword cConstant SealevelchangeViscousNumStepsEnum
    … …  
    409416syn keyword cConstant LoveTimeFreqEnum
    410417syn keyword cConstant LoveIsTimeEnum
     418syn keyword cConstant LoveHypergeomZEnum
     419syn keyword cConstant LoveHypergeomTable1Enum
     420syn keyword cConstant LoveHypergeomTable2Enum
    411421syn keyword cConstant SealevelchangeGSelfAttractionEnum
    412422syn keyword cConstant SealevelchangeGViscoElasticEnum
    … …  
    851861syn keyword cConstant SealevelEnum
    852862syn keyword cConstant SealevelGRDEnum
     863syn keyword cConstant SatGraviGRDEnum
    853864syn keyword cConstant SealevelBarystaticMaskEnum
    854865syn keyword cConstant SealevelBarystaticIceMaskEnum
    … …  
    890901syn keyword cConstant SealevelUNorthEsaEnum
    891902syn keyword cConstant SealevelchangeIndicesEnum
    892 syn keyword cConstant SealevelchangeGEnum
    893 syn keyword cConstant SealevelchangeGUEnum
    894 syn keyword cConstant SealevelchangeGEEnum
    895 syn keyword cConstant SealevelchangeGNEnum
     903syn keyword cConstant SealevelchangeAlphaIndexEnum
     904syn keyword cConstant SealevelchangeAzimuthIndexEnum
    896905syn keyword cConstant SealevelchangeGrotEnum
     906syn keyword cConstant SealevelchangeGSatGravirotEnum
    897907syn keyword cConstant SealevelchangeGUrotEnum
    898908syn keyword cConstant SealevelchangeGNrotEnum
    899909syn keyword cConstant SealevelchangeGErotEnum
    900 syn keyword cConstant SealevelchangeGsubelOceanEnum
    901 syn keyword cConstant SealevelchangeGUsubelOceanEnum
    902 syn keyword cConstant SealevelchangeGEsubelOceanEnum
    903 syn keyword cConstant SealevelchangeGNsubelOceanEnum
    904 syn keyword cConstant SealevelchangeGsubelIceEnum
    905 syn keyword cConstant SealevelchangeGUsubelIceEnum
    906 syn keyword cConstant SealevelchangeGEsubelIceEnum
    907 syn keyword cConstant SealevelchangeGNsubelIceEnum
    908 syn keyword cConstant SealevelchangeGsubelHydroEnum
    909 syn keyword cConstant SealevelchangeGUsubelHydroEnum
    910 syn keyword cConstant SealevelchangeGEsubelHydroEnum
    911 syn keyword cConstant SealevelchangeGNsubelHydroEnum
     910syn keyword cConstant SealevelchangeAlphaIndexOceanEnum
     911syn keyword cConstant SealevelchangeAlphaIndexIceEnum
     912syn keyword cConstant SealevelchangeAlphaIndexHydroEnum
     913syn keyword cConstant SealevelchangeAzimuthIndexOceanEnum
     914syn keyword cConstant SealevelchangeAzimuthIndexIceEnum
     915syn keyword cConstant SealevelchangeAzimuthIndexHydroEnum
    912916syn keyword cConstant SealevelchangeViscousRSLEnum
     917syn keyword cConstant SealevelchangeViscousSGEnum
    913918syn keyword cConstant SealevelchangeViscousUEnum
    914919syn keyword cConstant SealevelchangeViscousNEnum
    … …  
    12951300syn keyword cConstant DoubleArrayInputEnum
    12961301syn keyword cConstant ArrayInputEnum
     1302syn keyword cConstant IntArrayInputEnum
    12971303syn keyword cConstant DoubleExternalResultEnum
    12981304syn keyword cConstant DoubleMatArrayParamEnum
    … …  
    14101416syn keyword cConstant LovePMTF1tEnum
    14111417syn keyword cConstant LovePMTF2tEnum
     1418syn keyword cConstant LoveYiEnum
     1419syn keyword cConstant LoveRhsEnum
    14121420syn keyword cConstant LoveSolutionEnum
    14131421syn keyword cConstant MINIEnum
    … …  
    16241632syn keyword cType Cfsurfacesquare
    16251633syn keyword cType Channel
     1634syn keyword cType classes
    16261635syn keyword cType Constraint
    16271636syn keyword cType Constraints
    … …  
    16301639syn keyword cType ControlInput
    16311640syn keyword cType Covertree
     1641syn keyword cType DatasetInput
    16321642syn keyword cType DataSetParam
    1633 syn keyword cType DatasetInput
    16341643syn keyword cType Definition
    16351644syn keyword cType DependentObject
    … …  
    16441653syn keyword cType ElementInput
    16451654syn keyword cType ElementMatrix
     1655syn keyword cType Elements
    16461656syn keyword cType ElementVector
    1647 syn keyword cType Elements
    16481657syn keyword cType ExponentialVariogram
    16491658syn keyword cType ExternalResult
    … …  
    16521661syn keyword cType Friction
    16531662syn keyword cType Gauss
     1663syn keyword cType GaussianVariogram
     1664syn keyword cType gaussobjects
    16541665syn keyword cType GaussPenta
    16551666syn keyword cType GaussSeg
    16561667syn keyword cType GaussTetra
    16571668syn keyword cType GaussTria
    1658 syn keyword cType GaussianVariogram
    16591669syn keyword cType GenericExternalResult
    16601670syn keyword cType GenericOption
    … …  
    16651675syn keyword cType Input
    16661676syn keyword cType Inputs
     1677syn keyword cType IntArrayInput
    16671678syn keyword cType IntInput
    16681679syn keyword cType IntMatParam
    … …  
    16721683syn keyword cType IssmDirectApplicInterface
    16731684syn keyword cType IssmParallelDirectApplicInterface
     1685syn keyword cType krigingobjects
    16741686syn keyword cType Load
    16751687syn keyword cType Loads
    … …  
    16821694syn keyword cType Matice
    16831695syn keyword cType Matlitho
     1696syn keyword cType matrixobjects
    16841697syn keyword cType MatrixParam
    16851698syn keyword cType Misfit
    … …  
    16941707syn keyword cType Observations
    16951708syn keyword cType Option
     1709syn keyword cType Options
    16961710syn keyword cType OptionUtilities
    1697 syn keyword cType Options
    16981711syn keyword cType Param
    16991712syn keyword cType Parameters
    … …  
    17091722syn keyword cType Regionaloutput
    17101723syn keyword cType Results
     1724syn keyword cType Riftfront
    17111725syn keyword cType RiftStruct
    1712 syn keyword cType Riftfront
    17131726syn keyword cType SealevelGeometry
    17141727syn keyword cType Seg
    17151728syn keyword cType SegInput
     1729syn keyword cType Segment
    17161730syn keyword cType SegRef
    1717 syn keyword cType Segment
    17181731syn keyword cType SpcDynamic
    17191732syn keyword cType SpcStatic
    … …  
    17341747syn keyword cType Vertex
    17351748syn keyword cType Vertices
    1736 syn keyword cType classes
    1737 syn keyword cType gaussobjects
    1738 syn keyword cType krigingobjects
    1739 syn keyword cType matrixobjects
    17401749syn keyword cType AdjointBalancethickness2Analysis
    17411750syn keyword cType AdjointBalancethicknessAnalysis

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

    r27129 r27131  
    283283        LoveInnerCoreBoundaryEnum,
    284284        LoveComplexComputationEnum,
    285         LoveIntStepsPerLayerEnum,
     285        LoveQuadPrecisionEnum,
     286        LoveMinIntegrationStepsEnum,
     287        LoveMaxIntegrationdrEnum,
    286288        LoveKernelsEnum,
    287289        LoveMu0Enum,
    … …  
    295297        LoveUnderflowTolEnum,
    296298        LovePostWidderThresholdEnum,
     299        LoveDebugEnum,
     300        LoveHypergeomNZEnum,
     301        LoveHypergeomNAlphaEnum,
    297302        MassFluxSegmentsEnum,
    298303        MassFluxSegmentsPresentEnum,
    … …  
    381386        SolidearthSettingsElasticEnum,
    382387        SolidearthSettingsViscousEnum,
     388        SolidearthSettingsSatelliteGraviEnum,
     389        SolidearthSettingsDegreeAccuracyEnum,
    383390        SealevelchangeGeometryDoneEnum,
    384391        SealevelchangeViscousNumStepsEnum,
    … …  
    403410        LoveTimeFreqEnum,
    404411        LoveIsTimeEnum,
     412        LoveHypergeomZEnum,
     413        LoveHypergeomTable1Enum,
     414        LoveHypergeomTable2Enum,
    405415        SealevelchangeGSelfAttractionEnum,
    406416        SealevelchangeGViscoElasticEnum,
    … …  
    847857        SealevelEnum,
    848858        SealevelGRDEnum,
     859        SatGraviGRDEnum,
    849860        SealevelBarystaticMaskEnum,
    850861        SealevelBarystaticIceMaskEnum,
    … …  
    886897        SealevelUNorthEsaEnum,
    887898        SealevelchangeIndicesEnum,
    888         SealevelchangeGEnum,
    889         SealevelchangeGUEnum,
    890         SealevelchangeGEEnum,
    891         SealevelchangeGNEnum,
     899        SealevelchangeAlphaIndexEnum,
     900        SealevelchangeAzimuthIndexEnum,
    892901        SealevelchangeGrotEnum,
     902        SealevelchangeGSatGravirotEnum,
    893903        SealevelchangeGUrotEnum,
    894904        SealevelchangeGNrotEnum,
    895905        SealevelchangeGErotEnum,
    896         SealevelchangeGsubelOceanEnum,
    897         SealevelchangeGUsubelOceanEnum,
    898         SealevelchangeGEsubelOceanEnum,
    899         SealevelchangeGNsubelOceanEnum,
    900         SealevelchangeGsubelIceEnum,
    901         SealevelchangeGUsubelIceEnum,
    902         SealevelchangeGEsubelIceEnum,
    903         SealevelchangeGNsubelIceEnum,
    904         SealevelchangeGsubelHydroEnum,
    905         SealevelchangeGUsubelHydroEnum,
    906         SealevelchangeGEsubelHydroEnum,
    907         SealevelchangeGNsubelHydroEnum,
     906        SealevelchangeAlphaIndexOceanEnum,
     907        SealevelchangeAlphaIndexIceEnum,
     908        SealevelchangeAlphaIndexHydroEnum,
     909        SealevelchangeAzimuthIndexOceanEnum,
     910        SealevelchangeAzimuthIndexIceEnum,
     911        SealevelchangeAzimuthIndexHydroEnum,
    908912        SealevelchangeViscousRSLEnum,
     913        SealevelchangeViscousSGEnum,
    909914        SealevelchangeViscousUEnum,
    910915        SealevelchangeViscousNEnum,
    … …  
    14101415        LovePMTF1tEnum,
    14111416        LovePMTF2tEnum,
     1417        LoveYiEnum,
     1418        LoveRhsEnum,
    14121419        LoveSolutionEnum,
    14131420        MINIEnum,

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

    r27129 r27131  
    291291                case LoveInnerCoreBoundaryEnum : return "LoveInnerCoreBoundary";
    292292                case LoveComplexComputationEnum : return "LoveComplexComputation";
    293                 case LoveIntStepsPerLayerEnum : return "LoveIntStepsPerLayer";
     293                case LoveQuadPrecisionEnum : return "LoveQuadPrecision";
     294                case LoveMinIntegrationStepsEnum : return "LoveMinIntegrationSteps";
     295                case LoveMaxIntegrationdrEnum : return "LoveMaxIntegrationdr";
    294296                case LoveKernelsEnum : return "LoveKernels";
    295297                case LoveMu0Enum : return "LoveMu0";
    … …  
    303305                case LoveUnderflowTolEnum : return "LoveUnderflowTol";
    304306                case LovePostWidderThresholdEnum : return "LovePostWidderThreshold";
     307                case LoveDebugEnum : return "LoveDebug";
     308                case LoveHypergeomNZEnum : return "LoveHypergeomNZ";
     309                case LoveHypergeomNAlphaEnum : return "LoveHypergeomNAlpha";
    305310                case MassFluxSegmentsEnum : return "MassFluxSegments";
    306311                case MassFluxSegmentsPresentEnum : return "MassFluxSegmentsPresent";
    … …  
    389394                case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
    390395                case SolidearthSettingsViscousEnum : return "SolidearthSettingsViscous";
     396                case SolidearthSettingsSatelliteGraviEnum : return "SolidearthSettingsSatelliteGravi";
     397                case SolidearthSettingsDegreeAccuracyEnum : return "SolidearthSettingsDegreeAccuracy";
    391398                case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone";
    392399                case SealevelchangeViscousNumStepsEnum : return "SealevelchangeViscousNumSteps";
    … …  
    411418                case LoveTimeFreqEnum : return "LoveTimeFreq";
    412419                case LoveIsTimeEnum : return "LoveIsTime";
     420                case LoveHypergeomZEnum : return "LoveHypergeomZ";
     421                case LoveHypergeomTable1Enum : return "LoveHypergeomTable1";
     422                case LoveHypergeomTable2Enum : return "LoveHypergeomTable2";
    413423                case SealevelchangeGSelfAttractionEnum : return "SealevelchangeGSelfAttraction";
    414424                case SealevelchangeGViscoElasticEnum : return "SealevelchangeGViscoElastic";
    … …  
    853863                case SealevelEnum : return "Sealevel";
    854864                case SealevelGRDEnum : return "SealevelGRD";
     865                case SatGraviGRDEnum : return "SatGraviGRD";
    855866                case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask";
    856867                case SealevelBarystaticIceMaskEnum : return "SealevelBarystaticIceMask";
    … …  
    892903                case SealevelUNorthEsaEnum : return "SealevelUNorthEsa";
    893904                case SealevelchangeIndicesEnum : return "SealevelchangeIndices";
    894                 case SealevelchangeGEnum : return "SealevelchangeG";
    895                 case SealevelchangeGUEnum : return "SealevelchangeGU";
    896                 case SealevelchangeGEEnum : return "SealevelchangeGE";
    897                 case SealevelchangeGNEnum : return "SealevelchangeGN";
     905                case SealevelchangeAlphaIndexEnum : return "SealevelchangeAlphaIndex";
     906                case SealevelchangeAzimuthIndexEnum : return "SealevelchangeAzimuthIndex";
    898907                case SealevelchangeGrotEnum : return "SealevelchangeGrot";
     908                case SealevelchangeGSatGravirotEnum : return "SealevelchangeGSatGravirot";
    899909                case SealevelchangeGUrotEnum : return "SealevelchangeGUrot";
    900910                case SealevelchangeGNrotEnum : return "SealevelchangeGNrot";
    901911                case SealevelchangeGErotEnum : return "SealevelchangeGErot";
    902                 case SealevelchangeGsubelOceanEnum : return "SealevelchangeGsubelOcean";
    903                 case SealevelchangeGUsubelOceanEnum : return "SealevelchangeGUsubelOcean";
    904                 case SealevelchangeGEsubelOceanEnum : return "SealevelchangeGEsubelOcean";
    905                 case SealevelchangeGNsubelOceanEnum : return "SealevelchangeGNsubelOcean";
    906                 case SealevelchangeGsubelIceEnum : return "SealevelchangeGsubelIce";
    907                 case SealevelchangeGUsubelIceEnum : return "SealevelchangeGUsubelIce";
    908                 case SealevelchangeGEsubelIceEnum : return "SealevelchangeGEsubelIce";
    909                 case SealevelchangeGNsubelIceEnum : return "SealevelchangeGNsubelIce";
    910                 case SealevelchangeGsubelHydroEnum : return "SealevelchangeGsubelHydro";
    911                 case SealevelchangeGUsubelHydroEnum : return "SealevelchangeGUsubelHydro";
    912                 case SealevelchangeGEsubelHydroEnum : return "SealevelchangeGEsubelHydro";
    913                 case SealevelchangeGNsubelHydroEnum : return "SealevelchangeGNsubelHydro";
     912                case SealevelchangeAlphaIndexOceanEnum : return "SealevelchangeAlphaIndexOcean";
     913                case SealevelchangeAlphaIndexIceEnum : return "SealevelchangeAlphaIndexIce";
     914                case SealevelchangeAlphaIndexHydroEnum : return "SealevelchangeAlphaIndexHydro";
     915                case SealevelchangeAzimuthIndexOceanEnum : return "SealevelchangeAzimuthIndexOcean";
     916                case SealevelchangeAzimuthIndexIceEnum : return "SealevelchangeAzimuthIndexIce";
     917                case SealevelchangeAzimuthIndexHydroEnum : return "SealevelchangeAzimuthIndexHydro";
    914918                case SealevelchangeViscousRSLEnum : return "SealevelchangeViscousRSL";
     919                case SealevelchangeViscousSGEnum : return "SealevelchangeViscousSG";
    915920                case SealevelchangeViscousUEnum : return "SealevelchangeViscousU";
    916921                case SealevelchangeViscousNEnum : return "SealevelchangeViscousN";
    … …  
    14131418                case LovePMTF1tEnum : return "LovePMTF1t";
    14141419                case LovePMTF2tEnum : return "LovePMTF2t";
     1420                case LoveYiEnum : return "LoveYi";
     1421                case LoveRhsEnum : return "LoveRhs";
    14151422                case LoveSolutionEnum : return "LoveSolution";
    14161423                case MINIEnum : return "MINI";

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

    r27122 r27131  
    282282syn keyword juliaConstC LoveInnerCoreBoundaryEnum
    283283syn keyword juliaConstC LoveComplexComputationEnum
    284 syn keyword juliaConstC LoveIntStepsPerLayerEnum
     284syn keyword juliaConstC LoveQuadPrecisionEnum
     285syn keyword juliaConstC LoveMinIntegrationStepsEnum
     286syn keyword juliaConstC LoveMaxIntegrationdrEnum
    285287syn keyword juliaConstC LoveKernelsEnum
    286288syn keyword juliaConstC LoveMu0Enum
    … …  
    294296syn keyword juliaConstC LoveUnderflowTolEnum
    295297syn keyword juliaConstC LovePostWidderThresholdEnum
     298syn keyword juliaConstC LoveDebugEnum
     299syn keyword juliaConstC LoveHypergeomNZEnum
     300syn keyword juliaConstC LoveHypergeomNAlphaEnum
    296301syn keyword juliaConstC MassFluxSegmentsEnum
    297302syn keyword juliaConstC MassFluxSegmentsPresentEnum
    … …  
    380385syn keyword juliaConstC SolidearthSettingsElasticEnum
    381386syn keyword juliaConstC SolidearthSettingsViscousEnum
     387syn keyword juliaConstC SolidearthSettingsSatelliteGraviEnum
     388syn keyword juliaConstC SolidearthSettingsDegreeAccuracyEnum
    382389syn keyword juliaConstC SealevelchangeGeometryDoneEnum
    383390syn keyword juliaConstC SealevelchangeViscousNumStepsEnum
    … …  
    402409syn keyword juliaConstC LoveTimeFreqEnum
    403410syn keyword juliaConstC LoveIsTimeEnum
     411syn keyword juliaConstC LoveHypergeomZEnum
     412syn keyword juliaConstC LoveHypergeomTable1Enum
     413syn keyword juliaConstC LoveHypergeomTable2Enum
    404414syn keyword juliaConstC SealevelchangeGSelfAttractionEnum
    405415syn keyword juliaConstC SealevelchangeGViscoElasticEnum
    … …  
    844854syn keyword juliaConstC SealevelEnum
    845855syn keyword juliaConstC SealevelGRDEnum
     856syn keyword juliaConstC SatGraviGRDEnum
    846857syn keyword juliaConstC SealevelBarystaticMaskEnum
    847858syn keyword juliaConstC SealevelBarystaticIceMaskEnum
    … …  
    883894syn keyword juliaConstC SealevelUNorthEsaEnum
    884895syn keyword juliaConstC SealevelchangeIndicesEnum
    885 syn keyword juliaConstC SealevelchangeGEnum
    886 syn keyword juliaConstC SealevelchangeGUEnum
    887 syn keyword juliaConstC SealevelchangeGEEnum
    888 syn keyword juliaConstC SealevelchangeGNEnum
     896syn keyword juliaConstC SealevelchangeAlphaIndexEnum
     897syn keyword juliaConstC SealevelchangeAzimuthIndexEnum
    889898syn keyword juliaConstC SealevelchangeGrotEnum
     899syn keyword juliaConstC SealevelchangeGSatGravirotEnum
    890900syn keyword juliaConstC SealevelchangeGUrotEnum
    891901syn keyword juliaConstC SealevelchangeGNrotEnum
    892902syn keyword juliaConstC SealevelchangeGErotEnum
    893 syn keyword juliaConstC SealevelchangeGsubelOceanEnum
    894 syn keyword juliaConstC SealevelchangeGUsubelOceanEnum
    895 syn keyword juliaConstC SealevelchangeGEsubelOceanEnum
    896 syn keyword juliaConstC SealevelchangeGNsubelOceanEnum
    897 syn keyword juliaConstC SealevelchangeGsubelIceEnum
    898 syn keyword juliaConstC SealevelchangeGUsubelIceEnum
    899 syn keyword juliaConstC SealevelchangeGEsubelIceEnum
    900 syn keyword juliaConstC SealevelchangeGNsubelIceEnum
    901 syn keyword juliaConstC SealevelchangeGsubelHydroEnum
    902 syn keyword juliaConstC SealevelchangeGUsubelHydroEnum
    903 syn keyword juliaConstC SealevelchangeGEsubelHydroEnum
    904 syn keyword juliaConstC SealevelchangeGNsubelHydroEnum
     903syn keyword juliaConstC SealevelchangeAlphaIndexOceanEnum
     904syn keyword juliaConstC SealevelchangeAlphaIndexIceEnum
     905syn keyword juliaConstC SealevelchangeAlphaIndexHydroEnum
     906syn keyword juliaConstC SealevelchangeAzimuthIndexOceanEnum
     907syn keyword juliaConstC SealevelchangeAzimuthIndexIceEnum
     908syn keyword juliaConstC SealevelchangeAzimuthIndexHydroEnum
    905909syn keyword juliaConstC SealevelchangeViscousRSLEnum
     910syn keyword juliaConstC SealevelchangeViscousSGEnum
    906911syn keyword juliaConstC SealevelchangeViscousUEnum
    907912syn keyword juliaConstC SealevelchangeViscousNEnum
    … …  
    12881293syn keyword juliaConstC DoubleArrayInputEnum
    12891294syn keyword juliaConstC ArrayInputEnum
     1295syn keyword juliaConstC IntArrayInputEnum
    12901296syn keyword juliaConstC DoubleExternalResultEnum
    12911297syn keyword juliaConstC DoubleMatArrayParamEnum
    … …  
    14031409syn keyword juliaConstC LovePMTF1tEnum
    14041410syn keyword juliaConstC LovePMTF2tEnum
     1411syn keyword juliaConstC LoveYiEnum
     1412syn keyword juliaConstC LoveRhsEnum
    14051413syn keyword juliaConstC LoveSolutionEnum
    14061414syn keyword juliaConstC MINIEnum

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

    r27129 r27131  
    297297              else if (strcmp(name,"LoveInnerCoreBoundary")==0) return LoveInnerCoreBoundaryEnum;
    298298              else if (strcmp(name,"LoveComplexComputation")==0) return LoveComplexComputationEnum;
    299               else if (strcmp(name,"LoveIntStepsPerLayer")==0) return LoveIntStepsPerLayerEnum;
     299              else if (strcmp(name,"LoveQuadPrecision")==0) return LoveQuadPrecisionEnum;
     300              else if (strcmp(name,"LoveMinIntegrationSteps")==0) return LoveMinIntegrationStepsEnum;
     301              else if (strcmp(name,"LoveMaxIntegrationdr")==0) return LoveMaxIntegrationdrEnum;
    300302              else if (strcmp(name,"LoveKernels")==0) return LoveKernelsEnum;
    301303              else if (strcmp(name,"LoveMu0")==0) return LoveMu0Enum;
    … …  
    309311              else if (strcmp(name,"LoveUnderflowTol")==0) return LoveUnderflowTolEnum;
    310312              else if (strcmp(name,"LovePostWidderThreshold")==0) return LovePostWidderThresholdEnum;
     313              else if (strcmp(name,"LoveDebug")==0) return LoveDebugEnum;
     314              else if (strcmp(name,"LoveHypergeomNZ")==0) return LoveHypergeomNZEnum;
     315              else if (strcmp(name,"LoveHypergeomNAlpha")==0) return LoveHypergeomNAlphaEnum;
    311316              else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum;
    312317              else if (strcmp(name,"MassFluxSegmentsPresent")==0) return MassFluxSegmentsPresentEnum;
    … …  
    378383              else if (strcmp(name,"Modelname")==0) return ModelnameEnum;
    379384              else if (strcmp(name,"SamplingAlpha")==0) return SamplingAlphaEnum;
    380               else if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SamplingNumRequestedOutputs")==0) return SamplingNumRequestedOutputsEnum;
    381389              else if (strcmp(name,"SamplingRequestedOutputs")==0) return SamplingRequestedOutputsEnum;
    382390              else if (strcmp(name,"SamplingRobin")==0) return SamplingRobinEnum;
    383391              else if (strcmp(name,"SamplingSeed")==0) return SamplingSeedEnum;
    384392              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
     393              else if (strcmp(name,"SolidearthPartitionIce")==0) return SolidearthPartitionIceEnum;
    389394              else if (strcmp(name,"SolidearthPartitionHydro")==0) return SolidearthPartitionHydroEnum;
    390395              else if (strcmp(name,"SolidearthPartitionOcean")==0) return SolidearthPartitionOceanEnum;
    … …  
    398403              else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
    399404              else if (strcmp(name,"SolidearthSettingsViscous")==0) return SolidearthSettingsViscousEnum;
     405              else if (strcmp(name,"SolidearthSettingsSatelliteGravi")==0) return SolidearthSettingsSatelliteGraviEnum;
     406              else if (strcmp(name,"SolidearthSettingsDegreeAccuracy")==0) return SolidearthSettingsDegreeAccuracyEnum;
    400407              else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
    401408              else if (strcmp(name,"SealevelchangeViscousNumSteps")==0) return SealevelchangeViscousNumStepsEnum;
    … …  
    420427              else if (strcmp(name,"LoveTimeFreq")==0) return LoveTimeFreqEnum;
    421428              else if (strcmp(name,"LoveIsTime")==0) return LoveIsTimeEnum;
     429              else if (strcmp(name,"LoveHypergeomZ")==0) return LoveHypergeomZEnum;
     430              else if (strcmp(name,"LoveHypergeomTable1")==0) return LoveHypergeomTable1Enum;
     431              else if (strcmp(name,"LoveHypergeomTable2")==0) return LoveHypergeomTable2Enum;
    422432              else if (strcmp(name,"SealevelchangeGSelfAttraction")==0) return SealevelchangeGSelfAttractionEnum;
    423433              else if (strcmp(name,"SealevelchangeGViscoElastic")==0) return SealevelchangeGViscoElasticEnum;
    … …  
    496506              else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum;
    497507              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
    498               else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
    499512              else if (strcmp(name,"SmbIsmungsm")==0) return SmbIsmungsmEnum;
    500513              else if (strcmp(name,"SmbIsprecipscaled")==0) return SmbIsprecipscaledEnum;
    … …  
    506519              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    507520              else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
     521              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    512522              else if (strcmp(name,"SmbNumElevationBins")==0) return SmbNumElevationBinsEnum;
    513523              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    … …  
    619629              else if (strcmp(name,"Approximation")==0) return ApproximationEnum;
    620630              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
    621               else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum;
    622635              else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
    623636              else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum;
    … …  
    629642              else if (strcmp(name,"BasalforcingsFloatingiceMeltingRate")==0) return BasalforcingsFloatingiceMeltingRateEnum;
    630643              else if (strcmp(name,"BasalforcingsGeothermalflux")==0) return BasalforcingsGeothermalfluxEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
     644              else if (strcmp(name,"BasalforcingsGroundediceMeltingRate")==0) return BasalforcingsGroundediceMeltingRateEnum;
    635645              else if (strcmp(name,"BasalforcingsLinearBasinId")==0) return BasalforcingsLinearBasinIdEnum;
    636646              else if (strcmp(name,"BasalforcingsPerturbationMeltingRate")==0) return BasalforcingsPerturbationMeltingRateEnum;
    … …  
    742752              else if (strcmp(name,"EsaYmotion")==0) return EsaYmotionEnum;
    743753              else if (strcmp(name,"EtaDiff")==0) return EtaDiffEnum;
    744               else if (strcmp(name,"FlowequationBorderFS")==0) return FlowequationBorderFSEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"FlowequationBorderFS")==0) return FlowequationBorderFSEnum;
    745758              else if (strcmp(name,"FrictionAs")==0) return FrictionAsEnum;
    746759              else if (strcmp(name,"FrictionC")==0) return FrictionCEnum;
    … …  
    752765              else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
    753766              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
     767              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    758768              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    759769              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    … …  
    865875              else if (strcmp(name,"SampleOld")==0) return SampleOldEnum;
    866876              else if (strcmp(name,"SampleNoise")==0) return SampleNoiseEnum;
    867               else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
    868881              else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
    869882              else if (strcmp(name,"SamplingPhi")==0) return SamplingPhiEnum;
    … …  
    871884              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    872885              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
     886              else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum;
    873887              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    874888              else if (strcmp(name,"SealevelBarystaticIceMask")==0) return SealevelBarystaticIceMaskEnum;
    875889              else if (strcmp(name,"SealevelBarystaticIceWeights")==0) return SealevelBarystaticIceWeightsEnum;
    876890              else if (strcmp(name,"SealevelBarystaticIceArea")==0) return SealevelBarystaticIceAreaEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
     891              else if (strcmp(name,"SealevelBarystaticIceLatbar")==0) return SealevelBarystaticIceLatbarEnum;
    881892              else if (strcmp(name,"SealevelBarystaticIceLongbar")==0) return SealevelBarystaticIceLongbarEnum;
    882893              else if (strcmp(name,"SealevelBarystaticIceLoad")==0) return SealevelBarystaticIceLoadEnum;
    … …  
    913924              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    914925              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    915               else if (strcmp(name,"SealevelchangeG")==0) return SealevelchangeGEnum;
    916               else if (strcmp(name,"SealevelchangeGU")==0) return SealevelchangeGUEnum;
    917               else if (strcmp(name,"SealevelchangeGE")==0) return SealevelchangeGEEnum;
    918               else if (strcmp(name,"SealevelchangeGN")==0) return SealevelchangeGNEnum;
     926              else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
     927              else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
    919928              else if (strcmp(name,"SealevelchangeGrot")==0) return SealevelchangeGrotEnum;
     929              else if (strcmp(name,"SealevelchangeGSatGravirot")==0) return SealevelchangeGSatGravirotEnum;
    920930              else if (strcmp(name,"SealevelchangeGUrot")==0) return SealevelchangeGUrotEnum;
    921931              else if (strcmp(name,"SealevelchangeGNrot")==0) return SealevelchangeGNrotEnum;
    922932              else if (strcmp(name,"SealevelchangeGErot")==0) return SealevelchangeGErotEnum;
    923               else if (strcmp(name,"SealevelchangeGsubelOcean")==0) return SealevelchangeGsubelOceanEnum;
    924               else if (strcmp(name,"SealevelchangeGUsubelOcean")==0) return SealevelchangeGUsubelOceanEnum;
    925               else if (strcmp(name,"SealevelchangeGEsubelOcean")==0) return SealevelchangeGEsubelOceanEnum;
    926               else if (strcmp(name,"SealevelchangeGNsubelOcean")==0) return SealevelchangeGNsubelOceanEnum;
    927               else if (strcmp(name,"SealevelchangeGsubelIce")==0) return SealevelchangeGsubelIceEnum;
    928               else if (strcmp(name,"SealevelchangeGUsubelIce")==0) return SealevelchangeGUsubelIceEnum;
    929               else if (strcmp(name,"SealevelchangeGEsubelIce")==0) return SealevelchangeGEsubelIceEnum;
    930               else if (strcmp(name,"SealevelchangeGNsubelIce")==0) return SealevelchangeGNsubelIceEnum;
    931               else if (strcmp(name,"SealevelchangeGsubelHydro")==0) return SealevelchangeGsubelHydroEnum;
    932               else if (strcmp(name,"SealevelchangeGUsubelHydro")==0) return SealevelchangeGUsubelHydroEnum;
    933               else if (strcmp(name,"SealevelchangeGEsubelHydro")==0) return SealevelchangeGEsubelHydroEnum;
    934               else if (strcmp(name,"SealevelchangeGNsubelHydro")==0) return SealevelchangeGNsubelHydroEnum;
     933              else if (strcmp(name,"SealevelchangeAlphaIndexOcean")==0) return SealevelchangeAlphaIndexOceanEnum;
     934              else if (strcmp(name,"SealevelchangeAlphaIndexIce")==0) return SealevelchangeAlphaIndexIceEnum;
     935              else if (strcmp(name,"SealevelchangeAlphaIndexHydro")==0) return SealevelchangeAlphaIndexHydroEnum;
     936              else if (strcmp(name,"SealevelchangeAzimuthIndexOcean")==0) return SealevelchangeAzimuthIndexOceanEnum;
     937              else if (strcmp(name,"SealevelchangeAzimuthIndexIce")==0) return SealevelchangeAzimuthIndexIceEnum;
     938              else if (strcmp(name,"SealevelchangeAzimuthIndexHydro")==0) return SealevelchangeAzimuthIndexHydroEnum;
    935939              else if (strcmp(name,"SealevelchangeViscousRSL")==0) return SealevelchangeViscousRSLEnum;
     940              else if (strcmp(name,"SealevelchangeViscousSG")==0) return SealevelchangeViscousSGEnum;
    936941              else if (strcmp(name,"SealevelchangeViscousU")==0) return SealevelchangeViscousUEnum;
    937942              else if (strcmp(name,"SealevelchangeViscousN")==0) return SealevelchangeViscousNEnum;
    … …  
    993998              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    994999              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    995               else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    9961004              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    9971005              else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum;
    9981006              else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum;
    9991007              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
     1008              else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
    10041009              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
    10051010              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    … …  
    11161121              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    11171122              else if (strcmp(name,"Vx")==0) return VxEnum;
    1118               else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
    11191127              else if (strcmp(name,"VxObs")==0) return VxObsEnum;
    11201128              else if (strcmp(name,"VxShear")==0) return VxShearEnum;
    11211129              else if (strcmp(name,"VxSurface")==0) return VxSurfaceEnum;
    11221130              else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"VyBase")==0) return VyBaseEnum;
     1131              else if (strcmp(name,"VyBase")==0) return VyBaseEnum;
    11271132              else if (strcmp(name,"Vy")==0) return VyEnum;
    11281133              else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
    … …  
    12391244              else if (strcmp(name,"Outputdefinition91")==0) return Outputdefinition91Enum;
    12401245              else if (strcmp(name,"Outputdefinition92")==0) return Outputdefinition92Enum;
    1241               else if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Outputdefinition93")==0) return Outputdefinition93Enum;
    12421250              else if (strcmp(name,"Outputdefinition94")==0) return Outputdefinition94Enum;
    12431251              else if (strcmp(name,"Outputdefinition95")==0) return Outputdefinition95Enum;
    12441252              else if (strcmp(name,"Outputdefinition96")==0) return Outputdefinition96Enum;
    12451253              else if (strcmp(name,"Outputdefinition97")==0) return Outputdefinition97Enum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
     1254              else if (strcmp(name,"Outputdefinition98")==0) return Outputdefinition98Enum;
    12501255              else if (strcmp(name,"Outputdefinition99")==0) return Outputdefinition99Enum;
    12511256              else if (strcmp(name,"Outputdefinition9")==0) return Outputdefinition9Enum;
    … …  
    13621367              else if (strcmp(name,"FrontalForcingsRignotAutoregression")==0) return FrontalForcingsRignotAutoregressionEnum;
    13631368              else if (strcmp(name,"Fset")==0) return FsetEnum;
    1364               else if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"FullMeltOnPartiallyFloating")==0) return FullMeltOnPartiallyFloatingEnum;
    13651373              else if (strcmp(name,"GLheightadvectionAnalysis")==0) return GLheightadvectionAnalysisEnum;
    13661374              else if (strcmp(name,"GaussPenta")==0) return GaussPentaEnum;
    13671375              else if (strcmp(name,"GaussSeg")==0) return GaussSegEnum;
    13681376              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
     1377              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    13731378              else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    13741379              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    … …  
    14461451              else if (strcmp(name,"LovePMTF1t")==0) return LovePMTF1tEnum;
    14471452              else if (strcmp(name,"LovePMTF2t")==0) return LovePMTF2tEnum;
     1453              else if (strcmp(name,"LoveYi")==0) return LoveYiEnum;
     1454              else if (strcmp(name,"LoveRhs")==0) return LoveRhsEnum;
    14481455              else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    14491456              else if (strcmp(name,"MINI")==0) return MINIEnum;
    … …  
    14831490              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    14841491              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    1485               else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"MpiDense")==0) return MpiDenseEnum;
    14861496              else if (strcmp(name,"Mpi")==0) return MpiEnum;
    14871497              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    … …  
    14901500              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
    14911501              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
     1502              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    14961503              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    14971504              else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
    … …  
    16061613              else if (strcmp(name,"ThicknessErrorEstimator")==0) return ThicknessErrorEstimatorEnum;
    16071614              else if (strcmp(name,"TotalCalvingFluxLevelset")==0) return TotalCalvingFluxLevelsetEnum;
    1608               else if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"TotalCalvingMeltingFluxLevelset")==0) return TotalCalvingMeltingFluxLevelsetEnum;
    16091619              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    16101620              else if (strcmp(name,"TotalFloatingBmbScaled")==0) return TotalFloatingBmbScaledEnum;
    … …  
    16131623              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    16141624              else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
    1615          else stage=14;
    1616    }
    1617    if(stage==14){
    1618               if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
     1625              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    16191626              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
    16201627              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;

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

    r26358 r27131  
    2020                compute_bp_grd         = 0; %will GRD patterns for bottom pressures be computed?
    2121                degacc                 = 0; %degree increment for resolution of Green tables.
    22                 timeacc                = 0; %time step accuracy required to compute Green tables
     22                timeacc                = 1; %time step accuracy required to compute Green tables
    2323                horiz                  = 0; %compute horizontal deformation
    24                 grdmodel               = 0; %grd model (0 by default, 1 for (visco-)elastic, 2 for Ivins)
     24                grdmodel               = 1; %grd model (0 by default, 1 for (visco-)elastic, 2 for Ivins)
    2525                cross_section_shape    = 0; %cross section only used when grd model is Ivins
    2626        end
    … …  
    8989
    9090                        %no grd model by default:
    91                         self.grdmodel=0;
     91                        self.grdmodel=1;
    9292
    9393                end % }}}
    9494                function disp(self) % {{{
    9595                        disp(sprintf('   solidearth settings:'));
    96 
    97                         fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
    98                         fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
    99                         fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
    100                         fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
    101                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
    102                         fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
     96                        disp(sprintf('      core:'));
    10397                        fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
    104                         fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
    105                         fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
     98                        fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
     99                        fielddisplay(self,'runfrequency','How many time steps we let masstransport core accumulate changes before each run of the sealevelchange core (default: 1, i.e run slc every time step)');
     100                        disp(sprintf('      computational flags:'));
    106101                        fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
    107102                        fielddisplay(self,'elastic','enables elastic deformation from surface loading');
    108103                        fielddisplay(self,'viscous','enables viscous deformation from surface loading');
    109104                        fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
    110                         fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
     105                        fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
     106                        fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore. Used only for grdmodel=2 only');
     107                        disp(sprintf('      resolution:'));
     108                        fielddisplay(self,'degacc','spatial accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
    111109                        fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
    112                         fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
    113                         fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
     110                        disp(sprintf('      sea-level equation:'));
     111                        fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
     112                        fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
     113                        fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
     114                        fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
     115                        fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
     116                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
     117
    114118                end % }}}
    115119                function md = checkconsistency(self,md,solution,analyses) % {{{
    … …  
    133137                        if self.viscous==1 & self.elastic==0,
    134138                                error('solidearthsettings checkconsistency error message: need elastic on if viscous flag is set');
     139                        end
     140                        if self.rotation==1 & self.elastic==0,
     141                                error('solidearthsettings checkconsistency error message: need elastic on if rotation flag is set');
    135142                        end
    136143
    … …  
    164171                        WriteData(fid,prefix,'object',self,'fieldname','viscous','name','md.solidearth.settings.viscous','format','Boolean');
    165172                        WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.rotation','format','Boolean');
     173                        WriteData(fid,prefix,'object',self,'fieldname','rotation','name','md.solidearth.settings.satellitegravity','format','Boolean');
    166174                        WriteData(fid,prefix,'object',self,'fieldname','grdocean','name','md.solidearth.settings.grdocean','format','Boolean');
    167175                        WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','name','md.solidearth.settings.ocean_area_scaling','format','Boolean');
    … …  
    187195                        writejsdouble(fid,[modelname '.solidearth.settings.viscous'],self.viscous);
    188196                        writejsdouble(fid,[modelname '.solidearth.settings.rotation'],self.rotation);
    189                         writejsdouble(fid,[modelname '.solidearth.settings.grdocean'],self.rotation);
     197                        writejsdouble(fid,[modelname '.solidearth.settings.grdocean'],self.grdocean);
    190198                        writejsdouble(fid,[modelname '.solidearth.settings.ocean_area_scaling'],self.ocean_area_scaling);
    191199                        writejsdouble(fid,[modelname '.solidearth.settings.run_frequency'],self.run_frequency);
Note: See TracChangeset for help on using the changeset viewer.