Ignore:
Timestamp:
08/25/22 16:50:29 (3 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 27230

Location:
issm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

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

    r27035 r27232  
    376376        IssmDouble  time;
    377377        IssmDouble  coeff, indrate;
     378        IssmDouble  bed, bedrate = 1.0;
    378379
    379380        /*Retrieve all inputs and parameters we will need*/
    … …  
    382383        parameters->FindParam(&indrate,CalvingTestIndependentRateEnum,time);
    383384
     385        Input *bs_input = this->GetInput(BedEnum);_assert_(bs_input);
    384386        Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
    385387        Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
    … …  
    399401      lsf_slopey_input->GetInputValue(&dphidy,&gauss);
    400402
     403                bs_input->GetInputValue(&bed,&gauss);
     404                bedrate = (bed>0)?0.0:1.0;
     405
    401406      vel=sqrt(vx*vx + vy*vy) + 1e-14;
    402407      dphi=sqrt(dphidx*dphidx+dphidy*dphidy)+ 1e-14;
    403408
    404                 calvingratex[iv]= coeff*vx + indrate*dphidx/dphi;
    405                 calvingratey[iv]= coeff*vy + indrate*dphidy/dphi;
     409                calvingratex[iv]= coeff*vx + bedrate*indrate*dphidx/dphi;
     410                calvingratey[iv]= coeff*vy + bedrate*indrate*dphidy/dphi;
    406411                calvingrate[iv] = sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);
    407412        }
    … …  
    835840        IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
    836841        IssmDouble  sigma_vm[NUMVERTICES];
    837         IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
     842        IssmDouble  B, n;
    838843        IssmDouble  epse_2,groundedice,bed,sealevel;
    839         IssmDouble  arate, rho_ice, rho_water, thickness, paramX, Hab;
    840         int                     use_parameter=0;
    841         int                     nonlinear_law=0;
    842         IssmDouble  theta, alpha, midp, gamma;
     844        IssmDouble  arate, rho_ice, rho_water, thickness;
     845        int                     use_parameter=-1;
     846        IssmDouble  gamma, theta, alpha, xoffset, yoffset;
     847        IssmDouble  vel_lower, vel_upper, vrate, truncateVrate;
    843848
    844849        /* Get node coordinates and dof list: */
    … …  
    852857        Input *bs_input      = this->GetInput(BedEnum);                               _assert_(bs_input);
    853858        Input *H_input       = this->GetInput(ThicknessEnum);                         _assert_(H_input);
    854         Input *smax_fl_input = this->GetInput(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
    855         Input *smax_gr_input = this->GetInput(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
    856859        Input *n_input       = this->GetInput(MaterialsRheologyNEnum);                _assert_(n_input);
    857860        Input *sl_input      = this->GetInput(SealevelEnum);                          _assert_(sl_input);
    … …  
    864867        /* Use which parameter  */
    865868        this->FindParam(&use_parameter, CalvingUseParamEnum);
    866         this->FindParam(&theta, CalvingScaleThetaEnum);
    867         this->FindParam(&alpha, CalvingAmpAlphaEnum);
    868         this->FindParam(&midp, CalvingMidpointEnum);
    869         this->FindParam(&nonlinear_law, CalvingNonlinearLawEnum);
     869        this->FindParam(&theta, CalvingThetaEnum);
     870        this->FindParam(&alpha, CalvingAlphaEnum);
     871        this->FindParam(&xoffset, CalvingXoffsetEnum);
     872        this->FindParam(&yoffset, CalvingYoffsetEnum);
     873        this->FindParam(&vel_lower, CalvingVelLowerboundEnum);
     874        this->FindParam(&vel_upper, CalvingVelUpperboundEnum);
    870875
    871876        /* Start looping on the number of vertices: */
    … …  
    882887                bs_input->GetInputValue(&bed,gauss);
    883888                H_input->GetInputValue(&thickness,gauss);
    884                 smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
    885                 smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
    886889                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    887890                sl_input->GetInputValue(&sealevel,gauss);
    888891                arate_input->GetInputValue(&arate,gauss);
     892                vrate = 1.0;
     893                if (vel < vel_upper) vrate = vel / vel_upper;
    889894
    890895                /*Compute strain rate and viscosity: */
    … …  
    904909                sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));
    905910
    906                 /*Tensile stress threshold*/
    907                 if(groundedice<0)
    908                  sigma_max = sigma_max_floating;
    909                 else
    910                  sigma_max = sigma_max_grounded;
    911 
    912911                switch (use_parameter) {
    913912                        case 0:
    914                                 /* bed elevation */
    915                                 paramX = bed;
     913                                /* 0 Linear: f(x) = y_{o} + \alpha (x+x_{o}) */
     914                                gamma = yoffset = alpha * (bed+xoffset);
    916915                                break;
    917916                        case 1:
    918                                 /* Height above floatation */
    919                                 if (bed>sealevel)       paramX = 0.0;
    920                                 else paramX = thickness - (rho_water/rho_ice) * (sealevel-bed);
     917                                /* 1 tanh: f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
     918                                gamma = yoffset -  0.5*theta*tanh(alpha*(bed+xoffset));
    921919                                break;
    922920                        case 2:
    923                                 /* Thicknese */
    924                                 paramX = thickness;
     921                                /* 2 tanh(thicknes): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
     922                                gamma = yoffset -  0.5*theta*tanh(alpha*(-thickness+xoffset));
     923                                break;
     924                        case 3:
     925                                /* 3 tanh(normal vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
     926                                _error_("The normalized velocity is not supported yet!");
     927                                gamma = yoffset -  0.5*theta*tanh(alpha*(vel+xoffset));
    925928                                break;
    926929                        case 4:
    927                                 /* bed elevation+linear curve fitting */
    928                                 paramX = bed;
     930                                /* 4 tanh(truncated vel): f(x)=y_{o}-\frac{\theta}{2}\tanh(\alpha(x+x_{o})) */
     931                                truncateVrate = (min(vel_upper, max(vel_lower, vel))-vel_lower) / vel_upper;
     932                                gamma = 0.5*theta*(tanh(alpha*xoffset) - tanh(alpha*(truncateVrate+xoffset)));
    929933                                break;
    930934                        case -1:
    931                                 /* use nothing, just the arate*/
     935                                /* nothing, just the arate*/
     936                                gamma = 1;
    932937                                break;
    933938                        default:
    … …  
    935940                }
    936941
    937                 /* Compute the hyperbolic tangent function */
    938                 if ((use_parameter>-0.5) & (use_parameter<3)) {
    939                         gamma = 0.5*theta*(1.0-tanh((paramX+midp)/alpha))+(1.0-theta);
    940                         if (gamma<0.0) gamma =0.0;
    941                 }
    942                 else if (use_parameter>=3) {
    943                         gamma = alpha*paramX + theta;
    944                         if (gamma > 1.0) gamma = 1.0;
    945                         if (gamma < 0.0) gamma = 0.0;
    946                 }
    947                 else gamma = 1;
     942                /* set upper and lower bounds */
     943                if (gamma > 1.0) gamma = 1.0;
     944                if (gamma < 0.0) gamma = 0.0;
     945                if (bed >= sealevel) gamma = 0.0;
    948946
    949947                /*-------------------------------------------*/
    950                 if (nonlinear_law) {
    951                         /*This von Mises type has too strong positive feedback with vel included
    952                          * calvingrate[iv] = (arate+sigma_vm[iv]*vel/sigma_max)*gamma;
    953                          */
    954                         Hab = thickness - (rho_water/rho_ice) * (sealevel-bed);
    955                         if (Hab < 0.) Hab = 0.;
    956                         if (bed > sealevel) Hab = 0.;
    957 
    958                         calvingrate[iv] = (arate+Hab/sigma_max)*gamma;
    959                 }
     948                if (use_parameter < 3) {
     949                        calvingrate[iv] = arate*gamma*vrate;
     950                }
    960951                else {
    961952                        calvingrate[iv] = arate*gamma;
    962953                }
    963954        }
    964 
    965955        /*Add input*/
    966956        this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    … …  
    20182008        else for(int i=0;i<NUMVERTICES;i++)weights[i]=0;
    20192009
    2020         /*free ressources:*/
     2010        /*Free resources:*/
    20212011        delete gauss;
    20222012
    … …  
    39053895        this->AddInput(enum_type,values,this->element_type);
    39063896
    3907         /*Free ressources:*/
     3897        /*Free resources:*/
    39083898        xDelete<IssmDouble>(values);
    39093899        xDelete<int>(doflist);
    … …  
    60446034        pY->SetValues(gsize,indices,Y_values,ADD_VAL);
    60456035
    6046         /*free ressources:*/
     6036        /*Free resources:*/
    60476037        xDelete<int>(indices);
    60486038        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    … …  
    62066196        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
    62076197
    6208         /*free ressources:*/
     6198        /*Free resources:*/
    62096199        xDelete<int>(indices);
    62106200        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    … …  
    63136303        }
    63146304
    6315         /*Free ressources: */
     6305        /*Free resources: */
    63166306        xDelete<IssmDouble>(hes);
    63176307        xDelete<IssmDouble>(times);
    … …  
    63206310}
    63216311/*}}}*/
    6322 void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){ /*{{{*/
     6312void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/
    63236313
    63246314        /*Declarations:{{{*/
    … …  
    63426332
    63436333        #ifdef _HAVE_RESTRICT_
    6344         IssmDouble* __restrict__ G=NULL;
    6345         IssmDouble* __restrict__ GU=NULL;
    6346         IssmDouble* __restrict__ GN=NULL;
    6347         IssmDouble* __restrict__ GE=NULL;
    6348         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    63496334        IssmDouble* __restrict__ G_gravi_precomputed=NULL;
    6350         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
    6351         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    63526335        #else
    6353         IssmDouble* G=NULL;
    6354         IssmDouble* GU=NULL;
    6355         IssmDouble* GN=NULL;
    6356         IssmDouble* GE=NULL;
    6357         IssmDouble* G_viscoelastic_precomputed=NULL;
    63586336        IssmDouble* G_gravi_precomputed=NULL;
    6359         IssmDouble* U_viscoelastic_precomputed=NULL;
    6360         IssmDouble* H_viscoelastic_precomputed=NULL;
    63616337        #endif
    63626338
    … …  
    63646340        int index;
    63656341        int M;
     6342        IssmDouble degacc;
    63666343        IssmDouble doubleindex,lincoef;
    63676344
    … …  
    63806357        /*Rotational:*/
    63816358        #ifdef _HAVE_RESTRICT_
     6359        IssmDouble* __restrict__ Grot=NULL;
     6360        IssmDouble* __restrict__ GUrot=NULL;
     6361        IssmDouble* __restrict__ GNrot=NULL;
     6362        IssmDouble* __restrict__ GErot=NULL;
    63826363        IssmDouble* __restrict__ tide_love_h  = NULL;
    63836364        IssmDouble* __restrict__ tide_love_k  = NULL;
    … …  
    63866367        IssmDouble* __restrict__ LoveRotU     = NULL;
    63876368        IssmDouble* __restrict__ LoveRothoriz = NULL;
    6388         IssmDouble* __restrict__  Grot        = NULL;
    6389         IssmDouble* __restrict__ GUrot        = NULL;
    6390         IssmDouble* __restrict__ GNrot        = NULL;
    6391         IssmDouble* __restrict__ GErot        = NULL;
     6369        int* __restrict__ AplhaIndex   = NULL;
     6370        int* __restrict__ AzimuthIndex = NULL;
    63926371        #else
     6372        IssmDouble* Grot=NULL;
     6373        IssmDouble* GUrot=NULL;
     6374        IssmDouble* GNrot=NULL;
     6375        IssmDouble* GErot=NULL;
    63936376        IssmDouble* tide_love_h  = NULL;
    63946377        IssmDouble* tide_love_k  = NULL;
    … …  
    63976380        IssmDouble* LoveRotU     = NULL;
    63986381        IssmDouble* LoveRothoriz = NULL;
    6399         IssmDouble*  Grot        = NULL;
    6400         IssmDouble* GUrot        = NULL;
    6401         IssmDouble* GNrot        = NULL;
    6402         IssmDouble* GErot        = NULL;
     6382        int* AlphaIndex   = NULL;
     6383        int* AzimuthIndex = NULL;
    64036384        #endif
    64046385
    … …  
    64396420
    64406421        /*Recover precomputed green function kernels:{{{*/
    6441         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
    6442         parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
    6443 
     6422        parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
     6423        M=reCast<int,IssmDouble>(180.0/degacc+1.);
     6424
     6425        DoubleVecParam* parameter;
    64446426        if(computeelastic){
    6445                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    6446                 parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    6447 
    6448                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
    6449                 parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
    6450 
    6451                 if(horiz){
    6452                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    6453                         parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
    6454                 }
    6455 
    64566427                if(computerotation){
    64576428                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeTidalH2Enum)); _assert_(parameter);
    … …  
    64896460        }
    64906461
    6491         constant=3/rho_earth/planetarea;
    6492 
    6493         G=xNewZeroInit<IssmDouble>(3*nel*nt);
    6494         if(computeelastic){
    6495                 GU=xNewZeroInit<IssmDouble>(3*nel*nt);
    6496                 if(horiz){
    6497                         GN=xNewZeroInit<IssmDouble>(3*nel*nt);
    6498                         GE=xNewZeroInit<IssmDouble>(3*nel*nt);
    6499                 }
    6500         }
    6501 
    6502         for(int e=0;e<nel;e++){
    6503                 ratioe=constant*areae[e];
    6504                 for (int i=0;i<3;i++){
    6505                         IssmDouble alpha;
    6506                         IssmDouble delPhi,delLambda;
    6507                         /*recovers info for this element and vertex:*/
    6508                         IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
    6509                         IssmDouble longe= atan2(yye[e],xxe[e]);
    6510 
    6511                         lati=latitude[i];
    6512                         longi=longitude[i];
    6513 
    6514                         /*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
    6515                         delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    6516                         alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
    6517                         doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    6518                         index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    6519                         _assert_(index>=0 && index<M);
    6520 
    6521                         lincoef=doubleindex-index; //where between index and index+1 is alpha
    6522                         if (index==M-1){ //avoids out of bound case
    6523                                 index-=1;
    6524                                 lincoef=1;
    6525                         }
    6526 
    6527                         if(horiz){
    6528                                 /*Compute azimuths, both north and east components: */
    6529                                 x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    6530                                 if(lati==M_PI/2){
    6531                                         x=1e-12; y=1e-12;
     6462        AlphaIndex=xNewZeroInit<int>(3*nel);
     6463        if (horiz) AzimuthIndex=xNewZeroInit<int>(3*nel);
     6464        int intmax=pow(2,16)-1;
     6465
     6466
     6467        for (int i=0;i<3;i++){
     6468                if(lids[this->vertices[i]->lid]==this->lid){
     6469                        for(int e=0;e<nel;e++){
     6470                                IssmDouble alpha;
     6471                                IssmDouble delPhi,delLambda;
     6472                                /*recovers info for this element and vertex:*/
     6473                                IssmDouble late= asin(zze[e]/sqrt( pow(xxe[e],2.0)+ pow(yye[e],2.0)+ pow(zze[e],2.0)));
     6474                                IssmDouble longe= atan2(yye[e],xxe[e]);
     6475
     6476                                lati=latitude[i];
     6477                                longi=longitude[i];
     6478
     6479                                /*Computes alpha angle between centroid and current vertex, and indexes alpha in precomputed tables: */
     6480                                delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6481                                alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
     6482                                doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6483                                index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6484                                if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6485                                _assert_(index>=0 && index<M);
     6486
     6487                                if(horiz){
     6488                                        /*Compute azimuths*/
     6489                                        dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
     6490                                        dy=sin(longe-longi)*cos(late);
     6491                                        //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
     6492                                        AzimuthIndex[i*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    65326493                                }
    6533                                 if(lati==-M_PI/2){
    6534                                         x=1e-12; y=1e-12;
    6535                                 }
    6536                                 dx = xxe[e]-x; dy = yye[e]-y; dz = zze[e]-z;
    6537                                 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);
    6538                                 E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    6539                         }
    6540 
    6541                         for(int t=0;t<nt;t++){
    6542                                 int timeindex=i*nel*nt+e*nt+t;
    6543 
    6544                                 /*Rigid earth gravitational perturbation: */
    6545                                 G[timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef)
    6546                                              +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
    6547 
    6548                                 /*Find a way to interpolate precomputed Gkernels to our solution time stepping:*/
    6549                                 if(computeelastic){
    6550                                         G[timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6551                                                      +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
    6552                                 }
    6553                                 G[timeindex]=G[timeindex]*ratioe;
    6554 
    6555                                 /*Elastic components:*/
    6556                                 if(computeelastic){
    6557                                         GU[timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6558                                                                  +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6559                                         if(horiz){
    6560                                                 GN[timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6561                                                                                +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6562                                                 GE[timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    6563                                                                                +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    6564                                         }
    6565                                 }
    6566                         } //for(int t=0;t<nt;t++)
     6494                                AlphaIndex[i*nel+e]=index;
     6495                        }                       
    65676496                } //for (int i=0;i<3;i++)
    65686497        } //for(int e=0;e<nel;e++)
    65696498
    65706499        /*Add in inputs:*/
    6571         this->inputs->SetArrayInput(SealevelchangeGEnum,this->lid,G,nel*3*nt);
    6572         if(computeelastic){
    6573                 this->inputs->SetArrayInput(SealevelchangeGUEnum,this->lid,GU,nel*3*nt);
    6574                 if(horiz){
    6575                         this->inputs->SetArrayInput(SealevelchangeGNEnum,this->lid,GN,nel*3*nt);
    6576                         this->inputs->SetArrayInput(SealevelchangeGEEnum,this->lid,GE,nel*3*nt);
    6577                 }
    6578         }
     6500        this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3);
     6501        if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3);
     6502
    65796503        /*}}}*/
    65806504        /*Compute rotation kernel:{{{*/
    … …  
    65836507                LoveRotRSL  = xNewZeroInit<IssmDouble>(nt);
    65846508                LoveRotU    = xNewZeroInit<IssmDouble>(nt);
    6585                 LoveRothoriz= xNewZeroInit<IssmDouble>(nt);
     6509                if(horiz)LoveRothoriz= xNewZeroInit<IssmDouble>(nt);
    65866510                Grot        = xNewZeroInit<IssmDouble>(3*3*nt); //3 polar motion components * 3 vertices * number of time steps
    65876511                GUrot       = xNewZeroInit<IssmDouble>(3*3*nt);
    … …  
    66696593                        }
    66706594                }
     6595                /*Free resources:*/
     6596                xDelete<IssmDouble>(LoveRotRSL);
     6597                xDelete<IssmDouble>(LoveRotU);
     6598                if(horiz)xDelete<IssmDouble>(LoveRothoriz);
    66716599        }
    66726600        /*}}}*/
    … …  
    66906618        /*Free allocations:{{{*/
    66916619        #ifdef _HAVE_RESTRICT_
    6692         delete G;
    6693         if(computeelastic){
    6694                 delete GU;
    6695                 if(horiz){
    6696                         delete GN;
    6697                         delete GE;
    6698                 }
    6699                 if(computerotation){
    6700                         delete Grot;
    6701                         delete GUrot;
    6702                         if (horiz){
    6703                                 delete GNrot;
    6704                                 delete GErot;
    6705                         }
    6706                 }
    6707         }
     6620        delete AlphaIndex;
     6621        if(horiz) AzimuthIndex;
     6622
     6623        if(computerotation){
     6624                delete Grot;
     6625                delete GUrot;
     6626                if (horiz){
     6627                        delete GNrot;
     6628                        delete GErot;
     6629                }
     6630        }
     6631
    67086632        #else
    6709         xDelete(G);
    6710         if(computeelastic){
    6711                 xDelete(GU);
    6712                 if(horiz){
    6713                         xDelete(GN);
    6714                         xDelete(GE);
    6715                 }
    6716                 if(computerotation){
    6717                         xDelete(Grot);
    6718                         xDelete(GUrot);
    6719                         if (horiz){
    6720                                 xDelete(GNrot);
    6721                                 xDelete(GErot);
    6722                         }
     6633        xDelete<int>(AlphaIndex);
     6634        if(horiz){
     6635                xDelete<int>(AzimuthIndex);
     6636        }
     6637        if(computerotation){
     6638                xDelete<IssmDouble>(Grot);
     6639                xDelete<IssmDouble>(GUrot);
     6640                if (horiz){
     6641                        xDelete<IssmDouble>(GNrot);
     6642                        xDelete<IssmDouble>(GErot);
    67236643                }
    67246644        }
    67256645        #endif
     6646        /*}}}*/
     6647        return;
     6648
     6649}
     6650/*}}}*/
     6651
     6652void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
     6653
     6654        /*Declarations:{{{*/
     6655        int nel;
     6656        IssmDouble planetarea,planetradius;
     6657        IssmDouble constant,ratioe;
     6658        IssmDouble rho_earth;
     6659        IssmDouble lati,longi;
     6660        IssmDouble latitude[NUMVERTICES];
     6661        IssmDouble longitude[NUMVERTICES];
     6662        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
     6663        IssmDouble xyz_list[NUMVERTICES][3];
     6664
     6665        #ifdef _HAVE_RESTRICT_
     6666        int** __restrict__ AlphaIndex=NULL;
     6667        int** __restrict__ AzimIndex=NULL;
     6668
     6669        #else
     6670        int** AlphaIndex=NULL;
     6671        int** AzimIndex=NULL;
     6672        #endif
     6673
     6674        /*viscoelastic green function:*/
     6675        int index;
     6676        int M;
     6677        IssmDouble doubleindex,lincoef, degacc;
     6678
     6679        /*Computational flags:*/
     6680        bool computeselfattraction = false;
     6681        bool computeelastic = false;
     6682        bool computeviscous = false;
     6683        int  horiz;
     6684        int grd, grdmodel;
     6685
     6686        bool istime=true;
     6687        IssmDouble timeacc=0;
     6688        IssmDouble start_time,final_time;
     6689        int  nt,precomputednt;
     6690        int intmax=pow(2,16)-1;
     6691
     6692        /*}}}*/
     6693        /*Recover parameters:{{{ */
     6694        rho_earth=FindParam(MaterialsEarthDensityEnum);
     6695        this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
     6696        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     6697        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
     6698        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6699        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
     6700        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
     6701        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6702        this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
     6703        this->parameters->FindParam(&grdmodel,GrdModelEnum);
     6704        /*}}}*/
     6705
     6706        /*early return:*/
     6707        if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
     6708        if(!computeselfattraction)return;
     6709
     6710        /*Recover precomputed green function kernels:{{{*/
     6711        parameters->FindParam(&degacc,SolidearthSettingsDegreeAccuracyEnum);
     6712        M=reCast<int,IssmDouble>(180.0/degacc+1.);
     6713
     6714        /*}}}*/
     6715        /*Compute lat long of all vertices in the element:{{{*/
     6716        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     6717        for(int i=0;i<NUMVERTICES;i++){
     6718                latitude[i]= asin(xyz_list[i][2]/planetradius);
     6719                longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
     6720        }
     6721        /*}}}*/
     6722        /*Compute green functions:{{{ */
     6723
     6724        if(computeviscous){
     6725                this->parameters->FindParam(&istime,LoveIsTimeEnum);
     6726                if(!istime)_error_("Frequency love numbers not supported yet!");
     6727                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     6728                this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
     6729                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     6730                nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
     6731        }
     6732        else{
     6733                nt=1; //in elastic, or if we run only selfattraction, we need only one step
     6734        }
     6735        AlphaIndex=xNew<int*>(SLGEOM_NUMLOADS);
     6736        if(horiz) AzimIndex=xNew<int*>(SLGEOM_NUMLOADS);
     6737
     6738
     6739        //Allocate:
     6740        for(int l=0;l<SLGEOM_NUMLOADS;l++){
     6741                int nbar=slgeom->nbar[l];
     6742                AlphaIndex[l]=xNewZeroInit<int>(3*nbar);
     6743                if(horiz) AzimIndex[l]=xNewZeroInit<int>(3*nbar);
     6744
     6745
     6746                for (int i=0;i<3;i++){
     6747                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     6748                                for(int e=0;e<nbar;e++){
     6749                                        IssmDouble alpha;
     6750                                        IssmDouble delPhi,delLambda;
     6751                                        /*recover info for this element and vertex:*/
     6752                                        IssmDouble late= slgeom->latbarycentre[l][e];
     6753                                        IssmDouble longe= slgeom->longbarycentre[l][e];
     6754                                        late=late/180*M_PI;
     6755                                        longe=longe/180*M_PI;
     6756
     6757                                        lati=latitude[i];
     6758                                        longi=longitude[i];
     6759
     6760                                        if(horiz){
     6761                                                /*Compute azimuths*/
     6762                                                dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
     6763                                                dy=sin(longe-longi)*cos(late);
     6764                                                //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
     6765                                                AzimIndex[l][i*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6766                                        }
     6767
     6768                                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6769                                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6770                                        alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
     6771                                        doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6772                                        index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6773
     6774                                        if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6775                                        if (index==M-1){ //avoids out of bound case
     6776                                                index-=1;
     6777                                                lincoef=1;
     6778                                        }
     6779                                        AlphaIndex[l][i*nbar+e]=index;
     6780                                }
     6781                        }
     6782                }
     6783        }
     6784
     6785        /*Save all these arrayins for each element:*/
     6786        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6787                this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3);
     6788                if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3);
     6789        }
     6790        /*}}}*/
     6791        /*Free memory:{{{*/
     6792        for (int l=0;l<SLGEOM_NUMLOADS;l++){
     6793                xDelete<int>(AlphaIndex[l]);
     6794                if(horiz) xDelete<int>(AzimIndex[l]);
     6795        }
     6796        xDelete<int*>(AlphaIndex);
     6797        if(horiz) xDelete<int*>(AzimIndex);
     6798       
    67266799        /*}}}*/
    67276800        return;
    … …  
    69367009}
    69377010/*}}}*/
     7011void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
     7012
     7013        int nel;
     7014
     7015        /*Inputs:*/
     7016        IssmDouble I[NUMVERTICES];
     7017        IssmDouble W[NUMVERTICES];
     7018        IssmDouble BP[NUMVERTICES];
     7019        IssmDouble* areae=NULL;
     7020
     7021        /*output: */
     7022        IssmDouble bslcice=0;
     7023        IssmDouble bslchydro=0;
     7024        IssmDouble bslcbp=0;
     7025        IssmDouble BPavg=0;
     7026        IssmDouble Iavg=0;
     7027        IssmDouble Wavg=0;
     7028
     7029        /*ice properties: */
     7030        IssmDouble rho_ice,rho_water,rho_freshwater;
     7031
     7032        /*recover some parameters:*/
     7033        this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     7034        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     7035        this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
     7036        this->parameters->FindParam(&areae,&nel,AreaeEnum);
     7037
     7038        /*Retrieve inputs:*/
     7039        Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
     7040        Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
     7041        Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
     7042
     7043        for(int i=0;i<NUMVERTICES;i++){
     7044                Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid]*slgeom->LoadArea[SLGEOM_ICE][this->lid];
     7045                Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid]*slgeom->LoadArea[SLGEOM_WATER][this->lid];
     7046                BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]*slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     7047        }
     7048
     7049        /*convert from m^3 to kg:*/
     7050        Iavg*=rho_ice;
     7051        Wavg*=rho_freshwater;
     7052        BPavg*=rho_water;
     7053
     7054        #ifdef _ISSM_DEBUG_
     7055        this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
     7056        this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
     7057        this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
     7058        #endif
     7059
     7060        /*Compute barystatic component in kg:*/
     7061        // Note: Iavg, etc, already include partial area factor phi for subelement loading
     7062        bslcice =   -Iavg;
     7063        bslchydro = -Wavg;
     7064        bslcbp =    -BPavg;
     7065
     7066        _assert_(!xIsNan<IssmDouble>(bslcice));
     7067        _assert_(!xIsNan<IssmDouble>(bslchydro));
     7068        _assert_(!xIsNan<IssmDouble>(bslcbp));
     7069
     7070        /*Plug values into subelement load vector:*/
     7071        if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
     7072                int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
     7073                loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
     7074                Iavg=0; //avoid double counting centroid loads and subelement loads!
     7075        }
     7076        if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
     7077                int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
     7078                loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
     7079                Wavg=0;
     7080        }
     7081        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7082                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7083                loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL);
     7084                BPavg=0;
     7085        }
     7086        /*Plug remaining values into centroid load vector:*/
     7087        loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
     7088
     7089        /*Keep track of barystatic contributions:*/
     7090        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
     7091
     7092        /*Free ressources*/
     7093        xDelete<IssmDouble>(areae);
     7094
     7095}/*}}}*/
    69387096void       Tria::SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){ /*{{{*/
    69397097
    … …  
    70367194}
    70377195/*}}}*/
    7038 void       Tria::SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){ /*{{{*/
    7039 
    7040         /*Declarations:{{{*/
    7041         int nel;
    7042         IssmDouble planetarea,planetradius;
    7043         IssmDouble constant,ratioe;
    7044         IssmDouble rho_earth;
    7045         IssmDouble lati,longi;
    7046         IssmDouble latitude[NUMVERTICES];
    7047         IssmDouble longitude[NUMVERTICES];
    7048         IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    7049         IssmDouble xyz_list[NUMVERTICES][3];
    7050 
    7051         #ifdef _HAVE_RESTRICT_
    7052         IssmDouble* __restrict__ G_viscoelastic_precomputed=NULL;
    7053         IssmDouble* __restrict__ G_gravi_precomputed=NULL;
    7054         IssmDouble* __restrict__ U_viscoelastic_precomputed=NULL;
    7055         IssmDouble* __restrict__ H_viscoelastic_precomputed=NULL;
    7056         IssmDouble** __restrict__ Gsubel=NULL;
    7057         IssmDouble** __restrict__ GUsubel=NULL;
    7058         IssmDouble** __restrict__ GNsubel=NULL;
    7059         IssmDouble** __restrict__ GEsubel=NULL;
    7060 
    7061         #else
    7062         IssmDouble* G_viscoelastic_precomputed=NULL;
    7063         IssmDouble* G_gravi_precomputed=NULL;
    7064         IssmDouble* U_viscoelastic_precomputed=NULL;
    7065         IssmDouble* H_viscoelastic_precomputed=NULL;
    7066         IssmDouble** Gsubel=NULL;
    7067         IssmDouble** GUsubel=NULL;
    7068         IssmDouble** GNsubel=NULL;
    7069         IssmDouble** GEsubel=NULL;
    7070         #endif
    7071 
    7072         /*viscoelastic green function:*/
    7073         int index;
    7074         int M;
    7075         IssmDouble doubleindex,lincoef;
    7076 
    7077         /*Computational flags:*/
    7078         bool computeselfattraction = false;
    7079         bool computeelastic = false;
    7080         bool computeviscous = false;
    7081         int  horiz;
    7082         int grd, grdmodel;
    7083 
    7084         bool istime=true;
    7085         IssmDouble timeacc=0;
    7086         IssmDouble start_time,final_time;
    7087         int  nt,precomputednt;
    7088 
    7089         /*}}}*/
    7090         /*Recover parameters:{{{ */
    7091         rho_earth=FindParam(MaterialsEarthDensityEnum);
    7092         this->parameters->FindParam(&computeselfattraction,SolidearthSettingsSelfAttractionEnum);
    7093         this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
    7094         this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    7095         this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    7096         this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    7097         this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    7098         this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
    7099         this->parameters->FindParam(&grd,SolidearthSettingsGRDEnum);
    7100         this->parameters->FindParam(&grdmodel,GrdModelEnum);
    7101         /*}}}*/
    7102 
    7103         /*early return:*/
    7104         if (!grd || grdmodel!=ElasticEnum) return; //Love numbers won't be found in this case, return before loading them
    7105         if(!computeselfattraction)return;
    7106 
    7107         /*Recover precomputed green function kernels:{{{*/
    7108         DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGSelfAttractionEnum)); _assert_(parameter);
    7109         parameter->GetParameterValueByPointer((IssmDouble**)&G_gravi_precomputed,&M);
    7110 
    7111         if(computeelastic){
    7112                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    7113                 parameter->GetParameterValueByPointer((IssmDouble**)&G_viscoelastic_precomputed,NULL);
    7114 
    7115                 parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
    7116                 parameter->GetParameterValueByPointer((IssmDouble**)&U_viscoelastic_precomputed,NULL);
    7117 
    7118                 if(horiz){
    7119                         parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    7120                         parameter->GetParameterValueByPointer((IssmDouble**)&H_viscoelastic_precomputed,NULL);
    7121 
    7122                 }
    7123         }
    7124         /*}}}*/
    7125         /*Compute lat long of all vertices in the element:{{{*/
    7126         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7127         for(int i=0;i<NUMVERTICES;i++){
    7128                 latitude[i]= asin(xyz_list[i][2]/planetradius);
    7129                 longitude[i]= atan2(xyz_list[i][1],xyz_list[i][0]);
    7130         }
    7131         /*}}}*/
    7132         /*Compute green functions:{{{ */
    7133         constant=3/rho_earth/planetarea;
    7134 
    7135         if(computeviscous){
    7136                 this->parameters->FindParam(&istime,LoveIsTimeEnum);
    7137                 if(!istime)_error_("Frequency love numbers not supported yet!");
    7138                 this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
    7139                 this->parameters->FindParam(&start_time,TimesteppingStartTimeEnum);
    7140                 this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    7141                 nt=reCast<int,IssmDouble>((final_time-start_time)/timeacc)+1;
    7142         }
    7143         else{
    7144                 nt=1; //in elastic, or if we run only selfattraction, we need only one step
    7145         }
    7146         Gsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7147         if(computeelastic){
    7148                 GUsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7149                 if(horiz){
    7150                         GNsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7151                         GEsubel=xNew<IssmDouble*>(SLGEOM_NUMLOADS);
    7152                 }
    7153         }
    7154 
    7155         //Allocate:
    7156         for(int l=0;l<SLGEOM_NUMLOADS;l++){
    7157                 int nbar=slgeom->nbar[l];
    7158                 Gsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7159                 if(computeelastic){
    7160                         GUsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7161                         if(horiz){
    7162                                 GNsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7163                                 GEsubel[l]=xNewZeroInit<IssmDouble>(3*nbar*nt);
    7164                         }
    7165                 }
    7166 
    7167                 for(int e=0;e<nbar;e++){
    7168                         ratioe=constant*slgeom->area_subel[l][e];
    7169                         for (int i=0;i<3;i++){
    7170                                 IssmDouble alpha;
    7171                                 IssmDouble delPhi,delLambda;
    7172                                 /*recover info for this element and vertex:*/
    7173                                 IssmDouble late= slgeom->latbarycentre[l][e];
    7174                                 IssmDouble longe= slgeom->longbarycentre[l][e];
    7175                                 late=late/180*M_PI;
    7176                                 longe=longe/180*M_PI;
    7177 
    7178                                 lati=latitude[i];
    7179                                 longi=longitude[i];
    7180 
    7181                                 if(horiz){
    7182                                         /*Compute azimuths, both north and east components: */
    7183                                         x = xyz_list[i][0]; y = xyz_list[i][1]; z = xyz_list[i][2];
    7184                                         if(lati==90){
    7185                                                 x=1e-12; y=1e-12;
    7186                                         }
    7187                                         if(lati==-90){
    7188                                                 x=1e-12; y=1e-12;
    7189                                         }
    7190                                         IssmDouble xbar=planetradius*cos(late)*cos(longe);
    7191                                         IssmDouble ybar=planetradius*cos(late)*sin(longe);
    7192                                         IssmDouble zbar=planetradius*sin(late);
    7193 
    7194                                         dx = xbar-x; dy = ybar-y; dz = zbar-z;
    7195                                         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);
    7196                                         E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    7197                                 }
    7198 
    7199                                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    7200                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    7201                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
    7202                                 doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    7203                                 index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    7204 
    7205                                 lincoef=doubleindex-index; //where between index and index+1 is alpha
    7206                                 if (index==M-1){ //avoids out of bound case
    7207                                         index-=1;
    7208                                         lincoef=1;
    7209                                 }
    7210 
    7211                                 for(int t=0;t<nt;t++){
    7212                                         int timeindex=i*nbar*nt+e*nt+t;
    7213 
    7214                                         /*Rigid earth gravitational perturbation: */
    7215                                         Gsubel[l][timeindex]+=G_gravi_precomputed[index+0]*(1.-lincoef)
    7216                                                              +G_gravi_precomputed[index+1]*lincoef; //linear interpolation
    7217 
    7218                                         if(computeelastic){
    7219                                                 Gsubel[l][timeindex]+=G_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7220                                                                      +G_viscoelastic_precomputed[(index+1)*nt+t]*lincoef; //linear interpolation
    7221                                         }
    7222                                         Gsubel[l][timeindex]*=ratioe;
    7223 
    7224                                         /*Elastic components:*/
    7225                                         if(computeelastic){
    7226                                                 GUsubel[l][timeindex] =  ratioe *(U_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7227                                                                                  +U_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7228                                                 if(horiz){
    7229                                                         GNsubel[l][timeindex] = N_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7230                                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7231                                                         GEsubel[l][timeindex] = E_azim*ratioe *(H_viscoelastic_precomputed[(index+0)*nt+t]*(1.-lincoef)
    7232                                                                                               +H_viscoelastic_precomputed[(index+1)*nt+t]*lincoef);
    7233                                                 }
    7234                                         }
    7235                                 }
    7236                         }
    7237                 }
    7238         }
    7239 
    7240         /*Save all these arrayins for each element:*/
    7241         for (int l=0;l<SLGEOM_NUMLOADS;l++){
    7242                 this->inputs->SetArrayInput(slgeom->GEnum(l),this->lid,Gsubel[l],slgeom->nbar[l]*3*nt);
    7243                 if(computeelastic){
    7244                         this->inputs->SetArrayInput(slgeom->GUEnum(l),this->lid,GUsubel[l],slgeom->nbar[l]*3*nt);
    7245                         if(horiz){
    7246                                 this->inputs->SetArrayInput(slgeom->GNEnum(l),this->lid,GNsubel[l],slgeom->nbar[l]*3*nt);
    7247                                 this->inputs->SetArrayInput(slgeom->GEEnum(l),this->lid,GEsubel[l],slgeom->nbar[l]*3*nt);
    7248                         }
    7249                 }
    7250         }
    7251         /*}}}*/
    7252         /*Free memory:{{{*/
    7253         for (int l=0;l<SLGEOM_NUMLOADS;l++){
    7254                 xDelete<IssmDouble>(Gsubel[l]);
    7255                 if(computeelastic){
    7256                         xDelete<IssmDouble>(GUsubel[l]);
    7257                         if(horiz){
    7258                                 xDelete<IssmDouble>(GNsubel[l]);
    7259                                 xDelete<IssmDouble>(GEsubel[l]);
    7260                         }
    7261                 }
    7262         }
    7263         xDelete<IssmDouble*>(Gsubel);
    7264         if(computeelastic){
    7265                 xDelete<IssmDouble*>(GUsubel);
    7266                 if(horiz){
    7267                         xDelete<IssmDouble*>(GNsubel);
    7268                         xDelete<IssmDouble*>(GEsubel);
    7269                 }
    7270         }
    7271         /*}}}*/
    7272         return;
    7273 
    7274 }
    7275 /*}}}*/
    72767196void       Tria::SealevelchangeUpdateViscousFields(IssmDouble lincoeff, int newindex, int offset){ /*{{{*/
    72777197
    … …  
    72827202        IssmDouble* viscousE=NULL;
    72837203        int         viscousnumsteps;
    7284         int         dummy;
     7204        int         size;
    72857205        bool        viscous=false;
    72867206        int         horiz=0;
    … …  
    72927212                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    72937213
    7294                 this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&dummy);
    7295                 this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&dummy);
     7214                this->inputs->GetArrayPtr(SealevelchangeViscousRSLEnum,this->lid,&viscousRSL,&size);
     7215                this->inputs->GetArrayPtr(SealevelchangeViscousUEnum,this->lid,&viscousU,&size);
    72967216                if(horiz){
    7297                         this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&dummy);
    7298                         this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&dummy);
     7217                        this->inputs->GetArrayPtr(SealevelchangeViscousNEnum,this->lid,&viscousN,&size);
     7218                        this->inputs->GetArrayPtr(SealevelchangeViscousEEnum,this->lid,&viscousE,&size);
    72997219                }
    73007220
    … …  
    73077227                        }
    73087228                }
    7309 
    7310         }
    7311 
    7312 }
    7313 /*}}}*/
    7314 void       Tria::SealevelchangeBarystaticLoads(GrdLoads* loads,  BarystaticContributions* barycontrib, SealevelGeometry* slgeom){ /*{{{*/
    7315 
    7316         int nel;
    7317 
    7318         /*Inputs:*/
    7319         IssmDouble I[NUMVERTICES];
    7320         IssmDouble W[NUMVERTICES];
    7321         IssmDouble BP[NUMVERTICES];
    7322         IssmDouble* areae=NULL;
    7323 
    7324         /*output: */
    7325         IssmDouble bslcice=0;
    7326         IssmDouble bslchydro=0;
    7327         IssmDouble bslcbp=0;
    7328         IssmDouble BPavg=0;
    7329         IssmDouble Iavg=0;
    7330         IssmDouble Wavg=0;
    7331 
    7332         /*ice properties: */
    7333         IssmDouble rho_ice,rho_water,rho_freshwater;
    7334 
    7335         /*recover some parameters:*/
    7336         this->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     7229        }
     7230}
     7231/*}}}*/
     7232void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
     7233
     7234        IssmDouble oceanaverage=0;
     7235        IssmDouble oceanarea=0;
     7236        IssmDouble rho_water;
     7237
    73377238        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    7338         this->parameters->FindParam(&rho_freshwater,MaterialsRhoFreshwaterEnum);
    7339         this->parameters->FindParam(&areae,&nel,AreaeEnum);
    7340 
    7341         /*Retrieve inputs:*/
    7342         Element::GetInputListOnVertices(&I[0],DeltaIceThicknessEnum);
    7343         Element::GetInputListOnVertices(&W[0],DeltaTwsEnum);
    7344         Element::GetInputListOnVertices(&BP[0],DeltaBottomPressureEnum);
    7345 
     7239
     7240        /*retrieve ocean average and area:*/
    73467241        for(int i=0;i<NUMVERTICES;i++){
    7347                 Iavg+=I[i]*slgeom->LoadWeigths[SLGEOM_ICE][i][this->lid];
    7348                 Wavg+=W[i]*slgeom->LoadWeigths[SLGEOM_WATER][i][this->lid];
    7349                 BPavg+=BP[i]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    7350         }
    7351 
    7352         /*convert from m to kg/m^2:*/
    7353         Iavg*=rho_ice;
    7354         Wavg*=rho_freshwater;
    7355         BPavg*=rho_water;
    7356 
     7242                oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
     7243        }
    73577244        #ifdef _ISSM_DEBUG_
    7358         this->AddInput(SealevelBarystaticIceLoadEnum,&Iavg,P0Enum);
    7359         this->AddInput(SealevelBarystaticHydroLoadEnum,&Wavg,P0Enum);
    7360         this->AddInput(SealevelBarystaticBpLoadEnum,&BPavg,P0Enum);
     7245        this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    73617246        #endif
    7362 
    7363         /*Compute barystatic component in kg:*/
    7364         // Note: Iavg, etc, already include partial area factor phi for subelement loading
    7365         bslcice =   -slgeom->LoadArea[SLGEOM_ICE][this->lid]*Iavg;
    7366         bslchydro = -slgeom->LoadArea[SLGEOM_WATER][this->lid]*Wavg;
    7367         bslcbp =    -slgeom->LoadArea[SLGEOM_OCEAN][this->lid]*BPavg;
    7368 
    7369         _assert_(!xIsNan<IssmDouble>(bslcice));
    7370         _assert_(!xIsNan<IssmDouble>(bslchydro));
    7371         _assert_(!xIsNan<IssmDouble>(bslcbp));
    7372 
    7373         /*Plug values into subelement load vector:*/
    7374         if(slgeom->issubelement[SLGEOM_ICE][this->lid]){
    7375                 int intj=slgeom->subelementmapping[SLGEOM_ICE][this->lid];
    7376                 loads->vsubloads[SLGEOM_ICE]->SetValue(intj,Iavg,INS_VAL);
    7377                 Iavg=0; //avoid double counting centroid loads and subelement loads!
    7378         }
    7379         if(slgeom->issubelement[SLGEOM_WATER][this->lid]){
    7380                 int intj=slgeom->subelementmapping[SLGEOM_WATER][this->lid];
    7381                 loads->vsubloads[SLGEOM_WATER]->SetValue(intj,Wavg,INS_VAL);
    7382                 Wavg=0;
    7383         }
     7247        oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
     7248
     7249        /*add ocean average in the global sealevelloads vector:*/
    73847250        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    73857251                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7386                 loads->vsubloads[SLGEOM_OCEAN]->SetValue(intj,BPavg,INS_VAL);
    7387                 BPavg=0;
    7388         }
    7389         /*Plug remaining values into centroid load vector:*/
    7390         loads->vloads->SetValue(this->sid,Iavg+Wavg+BPavg,INS_VAL);
    7391 
    7392         /*Keep track of barystatic contributions:*/
    7393         barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
    7394 
     7252                loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL);
     7253                loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
     7254        }
     7255        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL);
     7256
     7257        /*add ocean area into a global oceanareas vector:*/
     7258        if(!loads->sealevelloads){
     7259                oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
     7260                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
     7261                        int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
     7262                        subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
     7263                }
     7264        }
    73957265}
    73967266/*}}}*/
    … …  
    73987268
    73997269        /*sal green function:*/
     7270        int* AlphaIndex=NULL;
     7271        int* AlphaIndexsub[SLGEOM_NUMLOADS];
    74007272        IssmDouble* G=NULL;
    74017273        IssmDouble* Grot=NULL;
    7402         IssmDouble* Gsub[SLGEOM_NUMLOADS];
     7274        DoubleVecParam* parameter;
    74037275        bool computefuture=false;
     7276        int spatial_component=0;
    74047277
    74057278        bool sal = false;
    … …  
    74177290
    74187291        if(sal){
    7419                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    7420                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
    7421                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    7422                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     7292                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     7293                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
     7294
     7295                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7296                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    74237297                if (rotation)   this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
    74247298
    7425                 this->SealevelchangeGxL(sealevelpercpu, G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    7426         }
    7427 
    7428         return;
    7429 } /*}}}*/
    7430 void       Tria::SealevelchangeOceanAverage(GrdLoads* loads, Vector<IssmDouble>* oceanareas, Vector<IssmDouble>* subelementoceanareas, IssmDouble* sealevelpercpu, SealevelGeometry* slgeom){ /*{{{*/
    7431 
    7432         IssmDouble oceanaverage=0;
    7433         IssmDouble oceanarea=0;
    7434         IssmDouble rho_water;
    7435 
    7436         this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
    7437 
    7438         /*retrieve ocean average and area:*/
    7439         for(int i=0;i<NUMVERTICES;i++){
    7440                 oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    7441         }
    7442         #ifdef _ISSM_DEBUG_
    7443         this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    7444         #endif
    7445         oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
    7446 
    7447         /*add ocean average in the global sealevelloads vector:*/
    7448         if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    7449                 int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7450                 loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water,INS_VAL);
    7451                 loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
    7452         }
    7453         else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    7454 
    7455         /*add ocean area into a global oceanareas vector:*/
    7456         if(!loads->sealevelloads){
    7457                 oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    7458                 if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    7459                         int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7460                         subelementoceanareas->SetValue(intj,oceanarea,INS_VAL);
    7461                 }
     7299                this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    74627300        }
    74637301
    … …  
    74737311        int nel,nbar;
    74747312        bool sal = false;
     7313        int* AlphaIndex=NULL;
     7314        int* AzimIndex=NULL;
     7315        int* AlphaIndexsub[SLGEOM_NUMLOADS];
     7316        int* AzimIndexsub[SLGEOM_NUMLOADS];
     7317        int spatial_component=0;
    74757318        IssmDouble* G=NULL;
    74767319        IssmDouble* GU=NULL;
    7477         IssmDouble* GE=NULL;
    7478         IssmDouble* GN=NULL;
     7320        IssmDouble* GH=NULL;
    74797321        IssmDouble* Grot=NULL;
    74807322        IssmDouble* GUrot=NULL;
    74817323        IssmDouble* GNrot=NULL;
    74827324        IssmDouble* GErot=NULL;
    7483         IssmDouble* Gsub[SLGEOM_NUMLOADS];
    7484         IssmDouble* GUsub[SLGEOM_NUMLOADS];
    7485         IssmDouble* GNsub[SLGEOM_NUMLOADS];
    7486         IssmDouble* GEsub[SLGEOM_NUMLOADS];
     7325
     7326        DoubleVecParam* parameter;
    74877327        bool computefuture=false;
    74887328
    … …  
    75037343
    75047344        if(sal){
    7505 
    7506                 this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    7507                 this->inputs->GetArrayPtr(SealevelchangeGsubelIceEnum,this->lid,&Gsub[SLGEOM_ICE],&size);
    7508                 this->inputs->GetArrayPtr(SealevelchangeGsubelHydroEnum,this->lid,&Gsub[SLGEOM_WATER],&size);
    7509                 this->inputs->GetArrayPtr(SealevelchangeGsubelOceanEnum,this->lid,&Gsub[SLGEOM_OCEAN],&size);
     7345                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7346                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
     7347
     7348                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
     7349                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    75107350
    75117351                if(elastic){
    7512                         this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
    7513                         this->inputs->GetArrayPtr(SealevelchangeGUsubelIceEnum,this->lid,&GUsub[SLGEOM_ICE],&size);
    7514                         this->inputs->GetArrayPtr(SealevelchangeGUsubelHydroEnum,this->lid,&GUsub[SLGEOM_WATER],&size);
    7515                         this->inputs->GetArrayPtr(SealevelchangeGUsubelOceanEnum,this->lid,&GUsub[SLGEOM_OCEAN],&size);
     7352                        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeUViscoElasticEnum)); _assert_(parameter);
     7353                        parameter->GetParameterValueByPointer((IssmDouble**)&GU,NULL);
     7354
    75167355                        if(horiz){
    7517                                 this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
    7518                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelIceEnum,this->lid,&GNsub[SLGEOM_ICE],&size);
    7519                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelHydroEnum,this->lid,&GNsub[SLGEOM_WATER],&size);
    7520                                 this->inputs->GetArrayPtr(SealevelchangeGNsubelOceanEnum,this->lid,&GNsub[SLGEOM_OCEAN],&size);
    7521 
    7522                                 this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
    7523                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelIceEnum,this->lid,&GEsub[SLGEOM_ICE],&size);
    7524                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelHydroEnum,this->lid,&GEsub[SLGEOM_WATER],&size);
    7525                                 this->inputs->GetArrayPtr(SealevelchangeGEsubelOceanEnum,this->lid,&GEsub[SLGEOM_OCEAN],&size);
     7356                                this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
     7357                                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
     7358
     7359                                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
     7360                                parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL);
    75267361                        }
    75277362                        if (rotation) {
    … …  
    75347369                        }
    75357370                }
    7536                 this->SealevelchangeGxL(&RSLGrd[0],G, Gsub, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
     7371
     7372                this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
    75377373
    75387374                if(elastic){
    7539                         this->SealevelchangeGxL(&UGrd[0],GU, GUsub, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
     7375                        this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
    75407376                        if(horiz){
    7541                                 this->SealevelchangeGxL(&NGrd[0],GN, GNsub, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
    7542                                 this->SealevelchangeGxL(&EGrd[0],GE, GEsub, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
     7377                                this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
     7378                                this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    75437379                        }
    75447380                }
    … …  
    75657401
    75667402} /*}}}*/
    7567 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) { /*{{{*/
     7403void       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) { /*{{{*/
    75687404
    75697405        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    75707406
    75717407        IssmDouble* grdfield=NULL;
    7572         int i,e,l,t,nbar;
     7408        int i,e,l,t,a, index, nbar;
     7409        bool rotation=false;
     7410        IssmDouble* Centroid_loads=NULL;
     7411        IssmDouble* Centroid_loads_copy=NULL;
     7412        IssmDouble* Subelement_loads[SLGEOM_NUMLOADS];
     7413        IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS];
     7414        IssmDouble* horiz_projection=NULL;
     7415        IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
     7416        int nt=1; //important, ensures there is a defined value if computeviscous is false
     7417
     7418        //viscous
    75737419        bool computeviscous=false;
    7574         bool rotation=false;
    7575         IssmDouble* viscousfield=NULL;
    7576         int nt=1; //important, ensures there is a defined value if computeviscous is false
    75777420        int viscousindex=0; //important
    75787421        int viscousnumsteps=1; //important
     7422        IssmDouble* viscousfield=NULL;
     7423        IssmDouble* grdfieldinterp=NULL;
     7424        IssmDouble* viscoustimes=NULL;
     7425        IssmDouble  final_time;
     7426        IssmDouble  lincoeff;
     7427        IssmDouble  timeacc;
    75797428
    75807429        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    … …  
    75837432                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    75847433                this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
     7434                this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
     7435                this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
     7436                this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
     7437                this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
    75857438                if(computefuture) {
    75867439                        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
    75877440                        if (nt>viscousnumsteps) nt=viscousnumsteps;
     7441                        grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps);
    75887442                }
    75897443                else nt=1;
    … …  
    76057459        }
    76067460
    7607 
    7608         if(loads->sealevelloads){ // general case: loads + sealevel loads
    7609                 for(i=0;i<NUMVERTICES;i++) {
    7610                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7611                         for (e=0;e<nel;e++){
    7612                                 for(t=0;t<nt;t++){
    7613                                         int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
    7614                                         grdfield[i*nt+t]+=G[index]*(loads->sealevelloads[e]+loads->loads[e]);
     7461        //Determine loads /*{{{*/
     7462        Centroid_loads=xNewZeroInit<IssmDouble>(nel);
     7463        for (e=0;e<nel;e++){
     7464                Centroid_loads[e]=loads->loads[e];
     7465        }
     7466        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7467                nbar=slgeom->nbar[l];
     7468                Subelement_loads[l]=xNewZeroInit<IssmDouble>(nbar);
     7469                for (e=0;e<nbar;e++){
     7470                        Subelement_loads[l][e]=(loads->subloads[l][e]);
     7471                }
     7472        }
     7473        if(loads->sealevelloads){
     7474                for (e=0;e<nel;e++){
     7475                        Centroid_loads[e]+=(loads->sealevelloads[e]);
     7476                }
     7477                nbar=slgeom->nbar[SLGEOM_OCEAN];
     7478                for (e=0;e<nbar;e++){
     7479                        Subelement_loads[SLGEOM_OCEAN][e]+=(loads->subsealevelloads[e]);
     7480                }
     7481        }
     7482
     7483        //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex
     7484        if (spatial_component!=0){
     7485                horiz_projection=xNewZeroInit<IssmDouble>(3*nel);
     7486                Centroid_loads_copy=xNewZeroInit<IssmDouble>(nel);
     7487                for (e=0;e<nel;e++){
     7488                        Centroid_loads_copy[e]=Centroid_loads[e];
     7489                }
     7490
     7491                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7492                        nbar=slgeom->nbar[l];
     7493                        Subelement_loads_copy[l]=xNewZeroInit<IssmDouble>(nbar);
     7494                        horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(3*nbar);
     7495                        for (e=0;e<nbar;e++){
     7496                                Subelement_loads_copy[l][e]=Subelement_loads[l][e];
     7497                        }
     7498                }
     7499        }
     7500        /*}}}*/
     7501
     7502        //Convolution
     7503        for(i=0;i<NUMVERTICES;i++) { /*{{{*/
     7504                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7505
     7506                if (spatial_component!=0){ //horizontals /*{{{*/
     7507                        //GxL needs to be projected on the right axis before summation into the grd field
     7508                        //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
     7509                        if (spatial_component==1){ //north
     7510                                for (e=0;e<nel;e++){
     7511                                        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
    76157512                                }
    7616                         }
     7513                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7514                                        nbar=slgeom->nbar[l];
     7515                                        for (e=0;e<nbar;e++){
     7516                                                horiz_projectionsub[l][i*nbar+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7517                                        }
     7518                                }
     7519                        }
     7520                        else if (spatial_component==2){ //east
     7521                                for (e=0;e<nel;e++){
     7522                                        horiz_projection[i*nel+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0);
     7523                                }
     7524                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7525                                        nbar=slgeom->nbar[l];
     7526                                        for (e=0;e<nbar;e++){
     7527                                                horiz_projectionsub[l][i*nbar+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7528                                        }
     7529                                }
     7530                        }
     7531                        for (e=0;e<nel;e++) Centroid_loads[e]=Centroid_loads_copy[e]*horiz_projection[i*nel+e];
    76177532                        for(l=0;l<SLGEOM_NUMLOADS;l++){
    76187533                                nbar=slgeom->nbar[l];
    76197534                                for (e=0;e<nbar;e++){
    7620                                         for(t=0;t<nt;t++){
    7621                                                 int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7622                                                 grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
    7623                                         }
     7535                                        Subelement_loads[l][e]=Subelement_loads_copy[l][e]*horiz_projectionsub[l][i*nbar+e];
    76247536                                }
    7625                                 if(l==SLGEOM_OCEAN){
    7626                                         for (e=0;e<nbar;e++){
    7627                                                 for(t=0;t<nt;t++){
    7628                                                         int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7629                                                         grdfield[i*nt+t]+=Gsub[l][index]*(loads->subsealevelloads[e]);
    7630                                                 }
    7631                                         }
     7537                        }
     7538                } /*}}}*/
     7539
     7540                for (e=0;e<nel;e++){
     7541                        for(t=0;t<nt;t++){
     7542                                a=AlphaIndex[i*nel+e];
     7543                                grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Centroid_loads[e];
     7544                        }
     7545                }
     7546                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7547                        nbar=slgeom->nbar[l];
     7548                        for (e=0;e<nbar;e++){
     7549                                for(t=0;t<nt;t++){
     7550                                        a=AlphaIndexsub[l][i*nbar+e];
     7551                                        grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Subelement_loads[l][e];
    76327552                                }
    76337553                        }
    76347554                }
    7635         }
    7636         else{  //this is the initial convolution where only loads are provided
    7637                 for(i=0;i<NUMVERTICES;i++) {
    7638                         if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7639                         for (e=0;e<nel;e++){
    7640                                 for(t=0;t<nt;t++){
    7641                                         int index=i*nel*viscousnumsteps+e*viscousnumsteps+t;
    7642                                         grdfield[i*nt+t]+=G[index]*(loads->loads[e]);
    7643                                 }
    7644                         }
    7645                         for(l=0;l<SLGEOM_NUMLOADS;l++){
    7646                                 nbar=slgeom->nbar[l];
    7647                                 for (e=0;e<nbar;e++){
    7648                                         for(t=0;t<nt;t++){
    7649                                                 int index=i*nbar*viscousnumsteps+e*viscousnumsteps+t;
    7650                                                 grdfield[i*nt+t]+=Gsub[l][index]*(loads->subloads[l][e]);
    7651                                         }
    7652                                 }
    7653                         }
    7654                 }
    7655         }
     7555        } /*}}}*/
     7556
     7557
    76567558
    76577559        if(computeviscous){ /*{{{*/
    … …  
    76617563                // 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
    76627564
    7663                 IssmDouble* grdfieldinterp=NULL;
    7664                 IssmDouble* viscoustimes=NULL;
    7665                 IssmDouble  final_time;
    7666                 IssmDouble  lincoeff;
    7667                 IssmDouble  timeacc;
    7668 
    7669                 this->parameters->FindParam(&viscoustimes,NULL,SealevelchangeViscousTimesEnum);
    7670                 this->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum);
    7671                 this->parameters->FindParam(&timeacc,SolidearthSettingsTimeAccEnum);
    7672                 this->inputs->GetArrayPtr(viscousenum,this->lid,&viscousfield,NULL);
    76737565                /* Map new grdfield generated by present-day loads onto viscous time vector*/
    76747566                if(computefuture){
    7675                         grdfieldinterp=xNewZeroInit<IssmDouble>(3*viscousnumsteps);
    76767567                        //viscousindex time and first time step of grdfield coincide, so just copy that value
    76777568                        for(int i=0;i<NUMVERTICES;i++){
    … …  
    76857576                                        for(int i=0;i<NUMVERTICES;i++){
    76867577                                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7687                                                 grdfieldinterp[i*viscousnumsteps+t]=  (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]+lincoeff*grdfield[i*nt+(t-viscousindex)];
     7578                                                grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]
     7579                                                                                         +lincoeff*grdfield[i*nt+(t-viscousindex)];
    76887580                                        }
    76897581                                }
    … …  
    76997591                /*update viscous stack with future deformation from present load: */
    77007592                if(computefuture){
    7701                         for(int t=viscousnumsteps-1;t>=viscousindex;t--){
     7593                        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
    77027594                                for(int i=0;i<NUMVERTICES;i++){
    77037595                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    77047596                                        //offset viscousfield to remove all deformation that has already been added
    7705                                         viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]-grdfieldinterp[i*viscousnumsteps+viscousindex]-viscousfield[i*viscousnumsteps+viscousindex];
     7597                                        viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]
     7598                                                                          -grdfieldinterp[i*viscousnumsteps+viscousindex]
     7599                                                                          -viscousfield[i*viscousnumsteps+viscousindex];
    77067600                                }
    77077601                        }
    77087602                        /*Save viscous stack now that we updated the values:*/
    77097603                        this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps);
    7710                 }
    7711 
    7712                 /*Free allocatoins:*/
    7713                 xDelete<IssmDouble>(viscoustimes);
     7604
     7605                        /*Free resources:*/
     7606                        xDelete<IssmDouble>(grdfieldinterp);
     7607                }
    77147608        }
    77157609        /*}}}*/
    … …  
    77267620                for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
    77277621        }
     7622        //free resources
     7623        xDelete<IssmDouble>(grdfield);
     7624        xDelete<IssmDouble>(Centroid_loads);
     7625        for(l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Subelement_loads[l]);
     7626        if (spatial_component!=0){
     7627                xDelete<IssmDouble>(horiz_projection);
     7628                xDelete<IssmDouble>(Centroid_loads_copy);
     7629                for(l=0;l<SLGEOM_NUMLOADS;l++) {
     7630                        xDelete<IssmDouble>(Subelement_loads_copy[l]);
     7631                        xDelete<IssmDouble>(horiz_projectionsub[l]);
     7632                }
     7633        }
     7634        if (computeviscous){
     7635                xDelete<IssmDouble>(viscoustimes);
     7636                if (computefuture){
     7637                        xDelete<IssmDouble>(grdfieldinterp);
     7638                }
     7639        }
    77287640
    77297641} /*}}}*/
    … …  
    77317643void       Tria::SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){ /*{{{*/
    77327644
     7645        offset*=slgeom->LoadArea[SLGEOM_OCEAN][this->lid]; //kg.m^-2 to kg
    77337646        if(slgeom->isoceanin[this->lid]){
    77347647                if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
Note: See TracChangeset for help on using the changeset viewer.