Changeset 26121


Ignore:
Timestamp:
03/19/21 18:55:34 (4 years ago)
Author:
Eric.Larour
Message:

CHG: fully fleshed out implementation of the new sea level convolution
based on element load distribution across the cluster.

Location:
issm/trunk-jpl
Files:
23 edited

Legend:

Unmodified
Added
Removed

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

    r26110 r26121  
    184184        }
    185185        /*New optimized code:*/
     186        ToolkitsOptionsFromAnalysis(parameters,SealevelchangeAnalysisEnum); //this is requested by the BarystaticContributions class inner vectors.
    186187        BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel);
    187188        parameters->AddObject(new GenericParam<BarystaticContributions*>(barystaticcontributions,BarystaticContributionsEnum));

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

    r26110 r26121  
    2626                cumice=new Vector<IssmDouble>(nice); cumice->Set(0); cumice->Assemble();
    2727        }
     28        else{
     29                ice=new Vector<IssmDouble>(1);
     30                cumice=new Vector<IssmDouble>(1);
     31        }
    2832
    2933        iomodel->FetchData(&nhydro,"md.solidearth.nparthydro");
    … …  
    3337                cumhydro=new Vector<IssmDouble>(nhydro); cumhydro->Set(0); cumhydro->Assemble();
    3438        }
     39        else{
     40                hydro=new Vector<IssmDouble>(1);
     41                cumhydro=new Vector<IssmDouble>(1);
     42        }
    3543        iomodel->FetchData(&nocean,"md.solidearth.npartocean");
    3644        if(nocean){
    … …  
    3846                ocean=new Vector<IssmDouble>(nocean);
    3947                cumocean=new Vector<IssmDouble>(nocean); cumocean->Set(0); cumocean->Assemble();
     48        }
     49        else{
     50                ocean=new Vector<IssmDouble>(1);
     51                cumocean=new Vector<IssmDouble>(1);
    4052        }
    4153
    … …  
    90102
    91103} /*}}}*/
     104void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
     105       
     106        int id;
     107        if(nice){
     108                id=reCast<int>(pice[eid]);
     109                ice->SetValue(id,icevalue,ADD_VAL);
     110        }
     111        else{
     112                ice->SetValue(0,icevalue,ADD_VAL);
     113        }
     114
     115        if(nhydro){
     116                id=reCast<int>(phydro[eid]);
     117                hydro->SetValue(id,hydrovalue,ADD_VAL);
     118        }
     119        else{
     120                hydro->SetValue(0,hydrovalue,ADD_VAL);
     121        }
     122
     123        if(nocean){
     124                id=reCast<int>(pocean[eid]);
     125                ocean->SetValue(id,oceanvalue,ADD_VAL);
     126        }
     127        else{
     128                ocean->SetValue(0,oceanvalue,ADD_VAL);
     129        }
     130               
     131
     132} /*}}}*/
    92133void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/
    93134
    … …  
    118159        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time));
    119160
    120         cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;i<nice;i++)cumice_serial[i]=cumice_serial[i]/oceanarea/rho_water;
    121         cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;i<nhydro;i++)cumhydro_serial[i]=cumhydro_serial[i]/oceanarea/rho_water;
    122         cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;i<nocean;i++)cumocean_serial[i]=cumocean_serial[i]/oceanarea/rho_water;
    123        
    124         results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time));
    125         results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time));
    126         results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time));
     161        if(nice){
     162                cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;i<nice;i++)cumice_serial[i]=cumice_serial[i]/oceanarea/rho_water;
     163                results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time));
     164        }
     165        if(nhydro){
     166                cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;i<nhydro;i++)cumhydro_serial[i]=cumhydro_serial[i]/oceanarea/rho_water;
     167                results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time));
     168        }
     169        if(nocean){
     170                cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;i<nocean;i++)cumocean_serial[i]=cumocean_serial[i]/oceanarea/rho_water;
     171                results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time));
     172        }
    127173
    128174        if(IssmComm::GetRank()==0){
    129                 xDelete<IssmDouble>(cumice_serial);
    130                 xDelete<IssmDouble>(cumhydro_serial);
    131                 xDelete<IssmDouble>(cumocean_serial);
     175                if(nice)xDelete<IssmDouble>(cumice_serial);
     176                if(nhydro)xDelete<IssmDouble>(cumhydro_serial);
     177                if(nocean)xDelete<IssmDouble>(cumocean_serial);
    132178        }
    133179        return;
    134180
    135181} /*}}}*/
    136 void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/
    137        
    138         int id;
    139                
    140         id=reCast<int>(pice[eid]);
    141         ice->SetValue(id,icevalue,ADD_VAL);
    142 
    143         id=reCast<int>(phydro[eid]);
    144         hydro->SetValue(id,hydrovalue,ADD_VAL);
    145 
    146         id=reCast<int>(pocean[eid]);
    147         ocean->SetValue(id,oceanvalue,ADD_VAL);
    148 
    149 } /*}}}*/

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

    r26110 r26121  
    1717        public:
    1818
    19                 Vector<IssmDouble>* ice;  //contributions to every ice partition
     19                Vector<IssmDouble>* ice;  //contributions to every ice partition (size nice x 1)
    2020                Vector<IssmDouble>* cumice;  //cumulated contributions to every ice partition
    2121                int                 nice; //number of ice partitions
    22                 IssmDouble*         pice; //ice partition
     22                IssmDouble*         pice; //ice partition (nel)
    2323
    24                 Vector<IssmDouble>* hydro;  //contributions to every hydro partition
     24                Vector<IssmDouble>* hydro;  //contributions to every hydro partition (size nhydro x 1)
    2525                Vector<IssmDouble>* cumhydro;  //cumulated contributions to every hydro partition
    2626                int                 nhydro; //number of hydro partitions
    27                 IssmDouble*         phydro; //hydro partition
     27                IssmDouble*         phydro; //hydro partition (nel)
    2828
    29                 Vector<IssmDouble>* ocean;  //contributions to every ocean partition
     29                Vector<IssmDouble>* ocean;  //contributions to every ocean partition (size nocean x 1)
    3030                Vector<IssmDouble>* cumocean;  //cumulated contributions to every ocean partition
    3131                int                 nocean; //number of ocean partitions
    32                 IssmDouble*         pocean; //ocean partition
    33 
     32                IssmDouble*         pocean; //ocean partition (nel)
     33 
    3434                /*BarystaticContributions constructors, destructors :*/
    3535                BarystaticContributions(IoModel* iomodel );

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

    r26110 r26121  
    382382                virtual IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0;
    383383                virtual void          SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0;
     384                virtual void          SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0;
    384385                virtual IssmDouble    SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector<IssmDouble>* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0;
    385386                virtual IssmDouble    SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0;
    … …  
    393394                virtual void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0;
    394395                virtual void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0;
    395                 virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0;
    396                 virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0;
     396                virtual void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0;
     397                virtual void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0;
     398                virtual void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,  SealevelMasks* masks)=0;
    397399                virtual void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0;
    398400                #endif

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

    r26110 r26121  
    219219                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    220220                void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
     221                void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    221222                void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    222223                                        IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");};
    … …  
    230231                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
    231232                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    232                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
    233                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     233                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
     234                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,SealevelMasks* masks){_error_("not implemented yet");};
     235                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    234236                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
    235237                #endif

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

    r26110 r26121  
    172172#ifdef _HAVE_SEALEVELCHANGE_
    173173                void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
     174                void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    174175                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    175176                void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    … …  
    185186                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
    186187                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    187                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
    188                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     188                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
     189                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
     190                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    189191                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
    190192

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

    r26110 r26121  
    180180                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    181181                void    SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
     182                void    SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");};
    182183                void    SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,
    183184                                        IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");};
    … …  
    191192                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");};
    192193                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");};
    193                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
    194                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");};
     194                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
     195                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
     196                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");};
    195197                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");};
    196198

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

    r26110 r26121  
    62526252        IssmDouble constant;
    62536253        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    6254         int sidlist[NUMVERTICES];
     6254        int sidlist[NUMVERTICES];
     6255        int sid;
    62556256
    62566257        #ifdef _HAVE_RESTRICT_
    … …  
    62776278        IssmDouble* longes=NULL;
    62786279        #endif
    6279 
     6280       
    62806281        /*elastic green function:*/
    62816282        int index;
    … …  
    62856286        bool computerigid = false;
    62866287        bool computeelastic = false;
     6288        bool computerotation = false;
    62876289        int  horiz;
     6290
     6291        /*Rotational:*/
     6292        IssmDouble* tide_love_h  = NULL;
     6293        IssmDouble* tide_love_k  = NULL;
     6294        IssmDouble* load_love_k  = NULL;
     6295        IssmDouble  tide_love_k2secular;
     6296        IssmDouble  moi_e, moi_p, omega, g;
     6297        IssmDouble Grotm1[3];
     6298        IssmDouble Grotm2[3];
     6299        IssmDouble Grotm3[3];
     6300        IssmDouble pre;
     6301
    62886302
    62896303        /*recover parameters: */
    … …  
    62916305        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
    62926306        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
     6307        this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum);
    62936308        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    62946309        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
    62956310        this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum);
    62966311        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6312
     6313        if(computerotation){
     6314                parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
     6315                parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
     6316                parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
     6317                parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
     6318                parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
     6319                parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
     6320                parameters->FindParam(&omega,RotationalAngularVelocityEnum);
     6321        }
    62976322
    62986323        /*compute area and add to inputs:*/
    … …  
    63496374                        IssmDouble late=lates[e];
    63506375                        IssmDouble longe=longes[e];
    6351                         int sid=sidlist[i];
     6376                        sid=sidlist[i];
    63526377
    63536378                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    … …  
    63876412                }
    63886413        }
     6414        if(computerotation){
     6415
     6416                for(int i=0;i<3;i++){
     6417                        sid=sidlist[i];
     6418                        lati=latitude[sid]/180.*M_PI;
     6419                        longi=longitude[sid]/180.*M_PI;
     6420
     6421                        pre=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*planetradius,2.0);
     6422                        Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi);
     6423                        Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi);
     6424                        Grotm3[i]= - pre* (1/6.0 - 0.5*cos(2.0*lati));
     6425                }
     6426                this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum);
     6427                this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum);
     6428                this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum);
     6429        }
     6430
    63896431
    63906432        /*Add in inputs:*/
    … …  
    64916533
    64926534        }
    6493        
    64946535        /*if we are an ice shelf, are we fully grounded or not? (used later):*/
    64956536        if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on.
     6537       
    64966538
    64976539        /*recover some parameters:*/
    … …  
    66066648}
    66076649/*}}}*/
    6608 void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
     6650void       Tria::SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
    66096651
    66106652        /*sal green function:*/
    … …  
    66126654        IssmDouble SealevelGRD[NUMVERTICES]={0,0,0};
    66136655        IssmDouble oceanaverage,oceanarea=0;
    6614 
     6656        IssmDouble rho_water;
     6657       
    66156658        bool sal = false;
     6659        bool rotation= false;
    66166660        int  size;
    66176661        int  nel;
     6662        IssmDouble Grotm1[3];
     6663        IssmDouble Grotm2[3];
     6664        IssmDouble Grotm3[3];
    66186665       
    66196666        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
     6667        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    66206668        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6669        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     6670
     6671        if(rotation){
     6672                Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
     6673                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
     6674                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
     6675               
     6676                for(int i=0;i<NUMVERTICES;i++) SealevelGRD[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
     6677        }
     6678
    66216679        if(sal){
    66226680                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
    … …  
    66336691       
    66346692        /*add ocean average in the global sealevelloads vector:*/
    6635         sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
     6693        sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
    66366694        oceanareas->SetValue(this->sid,oceanarea,INS_VAL);
    66376695       
    66386696        return;
    66396697} /*}}}*/
    6640 void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/
    6641 
    6642         bool converged=false;
     6698void       Tria::SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
     6699
    66436700        IssmDouble SealevelGRD[3]={0,0,0};
    66446701        IssmDouble oceanaverage,oceanarea=0;
     6702        IssmDouble rho_water;
    66456703        int nel;
    66466704        bool sal = false;
    66476705        IssmDouble* G=NULL;
    66486706        int size;
    6649        
     6707        IssmDouble Grotm1[3];
     6708        IssmDouble Grotm2[3];
     6709        IssmDouble Grotm3[3];
     6710        bool rotation= false;
     6711
    66506712        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
    66516713        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
     6714        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     6715        this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum);
     6716       
     6717        if(rotation){
     6718                Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
     6719                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
     6720                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
     6721               
     6722                for(int i=0;i<NUMVERTICES;i++) SealevelGRD[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
     6723        }       
     6724
     6725        if(sal){
     6726                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
     6727
     6728                for(int i=0;i<NUMVERTICES;i++) {
     6729                        for (int e=0;e<nel;e++){
     6730                                SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
     6731                        }
     6732                }
     6733        }
     6734        OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
     6735        newsealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL);
     6736
     6737} /*}}}*/
     6738void       Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/
     6739
     6740        IssmDouble Sealevel[3]={0,0,0};
     6741        IssmDouble SealevelRSL[3]={0,0,0};
     6742        IssmDouble SealevelU[3]={0,0,0};
     6743        IssmDouble SealevelN[3]={0,0,0};
     6744        IssmDouble SealevelE[3]={0,0,0};
     6745        int nel;
     6746        bool sal = false;
     6747        IssmDouble* G=NULL;
     6748        IssmDouble* GU=NULL;
     6749        IssmDouble* GE=NULL;
     6750        IssmDouble* GN=NULL;
     6751        int horiz;
     6752        int size;
     6753        IssmDouble Grotm1[3];
     6754        IssmDouble Grotm2[3];
     6755        IssmDouble Grotm3[3];
     6756        bool rotation= false;
     6757        bool elastic=false;
     6758
     6759        this->parameters->FindParam(&nel,MeshNumberofelementsEnum);
     6760        this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum);
     6761        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
     6762        this->parameters->FindParam(&elastic,SolidearthSettingsElasticEnum);
     6763        this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
     6764
     6765        if(rotation){
     6766                Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum);
     6767                Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum);
     6768                Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum);
     6769               
     6770                for(int i=0;i<NUMVERTICES;i++) SealevelRSL[i]+=Grotm1[i]*rotationvector[0]+Grotm2[i]*rotationvector[1]+Grotm3[i]*rotationvector[2];
     6771        }
    66526772               
    66536773        if(sal){
    66546774                this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size);
     6775                if(elastic) this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size);
     6776                if(horiz && elastic){
     6777                        this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size);
     6778                        this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size);
     6779                }
    66556780
    66566781                for(int i=0;i<NUMVERTICES;i++) {
    66576782                        for (int e=0;e<nel;e++){
    6658                                 SealevelGRD[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
    6659                         }
    6660                 }
    6661         }
    6662         OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks);
    6663         newsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
     6783                                SealevelRSL[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]);
     6784                        }
     6785                        if(elastic){
     6786                                for (int e=0;e<nel;e++){
     6787                                        SealevelU[i]+=GU[i*nel+e]*(sealevelloads[e]+loads[e]);
     6788                                }
     6789                        }
     6790                        if(horiz && elastic){
     6791                                for (int e=0;e<nel;e++){
     6792                                        SealevelN[i]+=GN[i*nel+e]*(sealevelloads[e]+loads[e]);
     6793                                        SealevelE[i]+=GE[i*nel+e]*(sealevelloads[e]+loads[e]);
     6794                                }
     6795                        }
     6796                }
     6797        }
     6798
     6799        /*Create geoid: */
     6800        for(int i=0;i<NUMVERTICES;i++)Sealevel[i]=SealevelU[i]+SealevelRSL[i];
     6801       
     6802        /*Create inputs*/
     6803        this->AddInput(SealevelGRDEnum,Sealevel,P1Enum);
     6804        this->AddInput(BedGRDEnum,SealevelU,P1Enum);
     6805        if(horiz){
     6806                this->AddInput(BedEastGRDEnum,SealevelE,P1Enum);
     6807                this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum);
     6808        }
     6809
    66646810
    66656811} /*}}}*/
    … …  
    67336879}
    67346880/*}}}*/
     6881void       Tria::SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/
     6882               
     6883        IssmDouble S;
     6884
     6885        /*Compute area of element:*/
     6886        IssmDouble area,planetarea;
     6887        this->Element::GetInputValue(&area,AreaEnum);
     6888
     6889        /*recover earth area: */
     6890        this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum);
     6891
     6892        /*Compute lat,long,radius of elemental centroid: */
     6893        bool spherical=true;
     6894        IssmDouble llr_list[NUMVERTICES][3];
     6895        IssmDouble late,longe,re;
     6896        /* Where is the centroid of this element?:{{{*/
     6897        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     6898
     6899        IssmDouble minlong=400;
     6900        IssmDouble maxlong=-20;
     6901        for (int i=0;i<NUMVERTICES;i++){
     6902                llr_list[i][0]=(90-llr_list[i][0]);
     6903                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     6904                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     6905                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     6906        }
     6907        if(minlong==0 && maxlong>180){
     6908                if (llr_list[0][1]==0)llr_list[0][1]=360;
     6909                if (llr_list[1][1]==0)llr_list[1][1]=360;
     6910                if (llr_list[2][1]==0)llr_list[2][1]=360;
     6911        }
     6912
     6913        // correction at the north pole
     6914        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     6915        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     6916        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     6917
     6918        //correction at the south pole
     6919        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     6920        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     6921        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     6922
     6923        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     6924        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     6925
     6926        late=90.-late;
     6927        if(longe>180.)longe=(longe-180.)-180.;
     6928
     6929        late=late/180.*M_PI;
     6930        longe=longe/180.*M_PI;
     6931        /*}}}*/
     6932        re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
     6933
     6934
     6935        /*recover total load: */
     6936        S=sealevelloads[this->Sid()]+loads[this->Sid()];
     6937
     6938        /* Perturbation terms for moment of inertia (moi_list):
     6939         * computed analytically (see Wu & Peltier, eqs 10 & 32)
     6940         * also consistent with my GMD formulation!
     6941         * ALL in geographic coordinates
     6942         * */
     6943        dI_list[0] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea;
     6944        dI_list[1] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea;
     6945        dI_list[2] = +4*M_PI*(S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea;
     6946        return;
     6947}/*}}}*/
    67356948#endif
    67366949

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

    r26110 r26121  
    163163                #ifdef _HAVE_SEALEVELCHANGE_
    164164                void       SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
     165                void       SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads);
    165166                void       SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    166167                IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector<IssmDouble>* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea);
    … …  
    175176                void       SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze);
    176177                void       SealevelchangeBarystaticLoads(Vector<IssmDouble>* loads, BarystaticContributions* barycontrib, SealevelMasks* masks);
    177                 void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, SealevelMasks* masks);
    178                 void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks);
     178                void       SealevelchangeInitialConvolution(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
     179                void       SealevelchangeOceanConvolution(Vector<IssmDouble>* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
     180                void       SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks);
    179181                void       OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks);
    180182                #endif

  • issm/trunk-jpl/src/c/classes/Inputs/Input.h

    r25507 r26121  
    4444
    4545                virtual void   AXPY(Input* xinput,IssmDouble scalar){_error_("Not implemented yet");};
     46                virtual void   Shift(IssmDouble scalar){_error_("Not implemented yet");};
    4647                virtual void   PointWiseMult(Input* xinput){_error_("Not implemented yet");};
    4748                virtual void   Pow(IssmDouble scale_factor){_error_("Not implemented yet");};

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

    r26099 r26121  
    247247}
    248248/*}}}*/
    249 void Inputs::AXPY(IssmDouble alpha, int xenum, int yenum, int zenum){/*{{{*/
     249void Inputs::ZAXPY(IssmDouble alpha, int xenum, int yenum, int zenum){/*{{{*/
    250250
    251251        _assert_(this);
    … …  
    268268        /*AXPY: */
    269269        this->inputs[index_z]->AXPY(this->inputs[index_x],alpha);
     270}
     271/*}}}*/
     272void Inputs::AXPY(IssmDouble alpha, int xenum, int yenum ){/*{{{*/
     273
     274        _assert_(this);
     275
     276        /*Get indices from enums*/
     277        int index_x = EnumToIndex(xenum);
     278        int index_y = EnumToIndex(yenum);
     279
     280        /*Make sure that old one exists*/
     281        if(!this->inputs[index_x]) _error_("Input "<<EnumToStringx(xenum)<<" not found");
     282        if(!this->inputs[index_y]) _error_("Input "<<EnumToStringx(yenum)<<" not found");
     283
     284        /*AXPY: */
     285        this->inputs[index_y]->AXPY(this->inputs[index_x],alpha);
     286}
     287/*}}}*/
     288void     Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/
     289
     290        _assert_(this);
     291
     292        /*Get indices from enums*/
     293        int index_x = EnumToIndex(xenum);
     294
     295        /*Make sure that x exists*/
     296        if(!this->inputs[index_x]) _error_("Input "<<EnumToStringx(xenum)<<" not found");
     297
     298        /*Shift: */
     299        this->inputs[index_x]->Shift(alpha);
    270300}
    271301/*}}}*/

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

    r26099 r26121  
    4646                int      DeleteInput(int enum_type);
    4747                void     DuplicateInput(int original_enum,int new_enum);
    48                 void     AXPY(IssmDouble alpha, int xenum, int yenum, int zenum);
     48                void     ZAXPY(IssmDouble alpha, int xenum, int yenum, int zenum);
     49                void     AXPY(IssmDouble alpha, int xenum, int yenum);
     50                void     Shift(int inputenum, IssmDouble alpha);
    4951                void     DeepEcho(void);
    5052                void     DeepEcho(int enum_in);

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

    r25508 r26121  
    384384}
    385385/*}}}*/
     386void TriaInput::Shift(IssmDouble alpha){/*{{{*/
     387
     388        /*Carry out the shift operation:*/
     389        for(int i=0;i<this->M*this->N;i++) this->values[i] +=alpha;
     390        for(int i=0;i<TriaRef::NumberofNodes(this->interpolation);i++) this->element_values[i] += alpha;
     391}
     392/*}}}*/
    386393void TriaInput::PointWiseMult(Input* xinput){/*{{{*/
    387394

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r25508 r26121  
    4040                void Pow(IssmDouble scalar);
    4141                void AXPY(Input* xinput,IssmDouble scalar);
     42                void Shift(IssmDouble scalar);
    4243                void PointWiseMult(Input* xinput);
    4344                void Serve(int numindices,int* indices);

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

    r26075 r26121  
    246246
    247247        /*From old and new thickness, create delta thickness  and accumulate:*/
    248         femmodel->inputs->AXPY(-1, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum);
    249         femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum,DummyEnum);
    250         femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaTwsEnum);
     248        femmodel->inputs->ZAXPY(-1, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum);
     249        femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum);
    251250
    252251        /*compute total water column change between two sea-level solver time steps, ie. every frequency*dt:*/
    253252        if(count==frequency){
    254                 femmodel->inputs->AXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum);
     253                femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum);
    255254                femmodel->inputs->DuplicateInput(AccumulatedDeltaTwsEnum,OldAccumulatedDeltaTwsEnum);
    256255        }

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

    r26075 r26121  
    117117
    118118        /*From old and new thickness, create delta thickness  and accumulate:*/
    119         femmodel->inputs->AXPY(-1, ThicknessOldEnum,ThicknessEnum,DeltaIceThicknessEnum);
    120         femmodel->inputs->AXPY(+1, DeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DummyEnum);
    121         femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaIceThicknessEnum);
     119        femmodel->inputs->ZAXPY(-1, ThicknessOldEnum,ThicknessEnum,DeltaIceThicknessEnum);
     120        femmodel->inputs->AXPY(+1, DeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum);
    122121
    123122        /*for Ivins deformation model, keep history of ice thickness changes inside TransientAccumulatedDeltaIceThicknessEnum:*/
    … …  
    154153        /*compute total ice thickness change between two sea-level solver time steps, ie. every frequency*dt:*/
    155154        if(count==frequency){
    156                 femmodel->inputs->AXPY(-1, OldAccumulatedDeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DeltaIceThicknessEnum);
     155                femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DeltaIceThicknessEnum);
    157156                femmodel->inputs->DuplicateInput(AccumulatedDeltaIceThicknessEnum,OldAccumulatedDeltaIceThicknessEnum);
    158157        }

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

    r26099 r26121  
    7575        /* From old and new bottom pressures, create delta bottom pressure, delta dsl and delta str.
    7676         * Accumulate delta bottom pressure: */
    77         femmodel->inputs->AXPY(-1, BottomPressureOldEnum,BottomPressureEnum,DeltaBottomPressureEnum);
    78         femmodel->inputs->AXPY(+1, DeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DummyEnum);
    79         femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaBottomPressureEnum);
     77        femmodel->inputs->ZAXPY(-1, BottomPressureOldEnum,BottomPressureEnum,DeltaBottomPressureEnum);
     78        femmodel->inputs->AXPY(+1, DeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum);
    8079
    81         femmodel->inputs->AXPY(-1, DslOldEnum,DslEnum,DeltaDslEnum);
    82         femmodel->inputs->AXPY(-1, StrOldEnum,StrEnum,DeltaStrEnum);
     80        femmodel->inputs->ZAXPY(-1, DslOldEnum,DslEnum,DeltaDslEnum);
     81        femmodel->inputs->ZAXPY(-1, StrOldEnum,StrEnum,DeltaStrEnum);
    8382
    8483        /* Compute total bottom pressure change between two sea-level solver time steps, ie. every frequency*dt. */
    8584        if(count==frequency){
    86                 femmodel->inputs->AXPY(-1, OldAccumulatedDeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DeltaBottomPressureEnum);
     85                femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DeltaBottomPressureEnum);
    8786                femmodel->inputs->DuplicateInput(AccumulatedDeltaBottomPressureEnum,OldAccumulatedDeltaBottomPressureEnum);
    8887        }

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

    r26119 r26121  
    1818IssmDouble SealevelloadsOceanAverage(Vector<IssmDouble>* sealevelloads,Vector<IssmDouble>* oceanareas, IssmDouble oceanarea);
    1919SealevelMasks* sealevel_masks(FemModel* femmodel);
     20void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads);
    2021/*}}}*/
    2122
    … …  
    398399        SealevelMasks* masks=NULL;
    399400        GenericParam<BarystaticContributions*>* barycontribparam=NULL;
     401        IssmDouble rotationaxismotionvector[3];
    400402       
    401403        Vector<IssmDouble>*    loads=NULL;
    … …  
    413415        int                  step;
    414416        IssmDouble           time;
     417        bool converged=false;
    415418
    416419        int  modelid,earthid;
    … …  
    418421        int  count,frequency,iscoupling;
    419422        int  grd=0;
     423        int computesealevel=0;
    420424        /*}}}*/
    421425
    … …  
    427431        femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
    428432        femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum);
     433        femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum);
    429434        /*}}}*/
    430435
    … …  
    463468                element->SealevelchangeBarystaticLoads(loads, barycontrib,masks);
    464469        }
     470        loads->Assemble();
    465471
    466472        //broadcast loads
    467         loads->Assemble(); allloads=loads->ToMPISerial();
     473        allloads=loads->ToMPISerial();
     474
     475        //compute rotation axis motion:
     476        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL);
     477
     478        /*skip computation of sea level if requested, which means sea level loads should be zeroed */
     479        if(!computesealevel){
     480                allsealevelloads=xNewZeroInit<IssmDouble>(nel);
     481                goto deformation;
     482        }
    468483
    469484        /*convolve loads:*/
    470485        for(Object* & object : femmodel->elements->objects){
    471486                Element* element = xDynamicCast<Element*>(object);
    472                 element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks);
     487                element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,rotationaxismotionvector,masks);
    473488        }
    474489        sealevelloads->Assemble();
    475        
     490
    476491        //Get ocean area:
    477492        oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.);
    … …  
    479494
    480495        //substract ocean average and barystatic contributionfrom sea level loads:
    481         sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
     496        sealevelloads->Shift(barycontrib->Total()/oceanarea - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
    482497       
    483498        //broadcast sea level loads
    484499        allsealevelloads=sealevelloads->ToMPISerial();
    485500
    486         bool converged=false;
     501        //compute rotation axis motion:
     502        RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
     503
    487504        for(;;){
    488505                       
    … …  
    492509                for(Object* & object : femmodel->elements->objects){
    493510                        Element* element = xDynamicCast<Element*>(object);
    494                         element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks);
     511                        element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,rotationaxismotionvector,masks);
    495512                }
    496513                sealevelloads->Assemble();
    497514       
    498                
    499515                //substract ocean average and barystatic contribution
    500                 sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
     516                sealevelloads->Shift(barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
    501517               
    502518                //broadcast sea level loads
    503519                allsealevelloads=sealevelloads->ToMPISerial();
    504520
     521                //compute rotation axis motion:
     522                RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads);
     523
    505524                //convergence?
    506525                slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs);
    … …  
    508527        }
    509528
     529        deformation:
     530
    510531        /*convolve loads and sea level loads to get the deformation:*/
    511532        for(Object* & object : femmodel->elements->objects){
    512533                Element* element = xDynamicCast<Element*>(object);
    513                 //element->SealevelchangeDeformationConvolution(allsealevelloads, allloads,masks);
     534                element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector,masks);
    514535        }
    515536
    516537        /*Update bedrock motion and geoid:*/
    517         /*inputs->AXPY(SealevelEnum,SealevelGRD,1);
    518         inputs->AXPY(BedEnum,BedGRD,1);
     538        if(computesealevel){
     539                femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea));
     540
     541                /*cumulate barystatic contributions and save to results: */
     542                barycontrib->Cumulate(femmodel->parameters);
     543                barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
     544        }
     545               
     546        femmodel->inputs->AXPY(SealevelEnum,SealevelGRDEnum,1);
     547        femmodel->inputs->AXPY(BedEnum,BedGRDEnum,1);
    519548        if(horiz){
    520                 inputs->AXPY(BedrockeastEnum,BedEastGRD,1);
    521                 inputs->AXPY(BedrocknorthEnum,BedNorthGRD,1);
    522         }*/
    523 
    524         /*cumulate barystatic contributions and save to results: */
    525         barycontrib->Cumulate(femmodel->parameters);
    526         barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea);
    527 
     549                femmodel->inputs->AXPY(BedEastEnum,BedEastGRDEnum,1);
     550                femmodel->inputs->AXPY(BedNorthEnum,BedNorthGRDEnum,1);
     551        }
    528552
    529553}
    … …  
    12161240        return sealevelloadsaverage/oceanarea;
    12171241} /*}}}*/
     1242void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads){ /*{{{*/
     1243
     1244        IssmDouble  moi_list[3]={0,0,0};
     1245        IssmDouble  moi_list_cpu[3]={0,0,0};
     1246        IssmDouble*     tide_love_h  = NULL;
     1247        IssmDouble*     tide_love_k  = NULL;
     1248        IssmDouble*     load_love_k  = NULL;
     1249        IssmDouble  tide_love_k2secular;
     1250        IssmDouble  moi_e, moi_p;
     1251        IssmDouble      m1, m2, m3;
     1252
     1253        /*retrieve parameters: */
     1254        femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum);
     1255        femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum);
     1256        femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum);
     1257        femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum);
     1258        femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum);
     1259        femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum);
     1260
     1261        for(Object* & object : femmodel->elements->objects){
     1262                Element* element = xDynamicCast<Element*>(object);
     1263                element->SealevelchangeMomentOfInertiaElement(&moi_list[0],loads,sealeveloads);
     1264                moi_list_cpu[0] += moi_list[0];
     1265                moi_list_cpu[1] += moi_list[1];
     1266                moi_list_cpu[2] += moi_list[2];
     1267        }
     1268
     1269        ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1270        ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1271       
     1272        ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1273        ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1274       
     1275        ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     1276        ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     1277
     1278        /*compute perturbation terms for angular velocity vector: */
     1279        m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0];
     1280        m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1];
     1281        m3 = -(1+load_love_k[2])/moi_p * moi_list[2];   // term associated with fluid number (3-order-of-magnitude smaller) is negelected
     1282
     1283        /*Assign output pointers:*/
     1284        m[0]=m1;
     1285        m[1]=m2;
     1286        m[2]=m3;
     1287} /*}}}*/

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

    r26118 r26121  
    353353syn keyword cConstant SolidearthSettingsElasticEnum
    354354syn keyword cConstant SealevelchangeGeometryDoneEnum
     355syn keyword cConstant SealevelGrotm1Enum
     356syn keyword cConstant SealevelGrotm2Enum
     357syn keyword cConstant SealevelGrotm3Enum
    355358syn keyword cConstant RotationalEquatorialMoiEnum
    356359syn keyword cConstant TidalLoveHEnum
    … …  
    550553syn keyword cConstant BaseSlopeYEnum
    551554syn keyword cConstant BedEnum
     555syn keyword cConstant BedGRDEnum
    552556syn keyword cConstant BedEastEnum
     557syn keyword cConstant BedEastGRDEnum
    553558syn keyword cConstant BedNorthEnum
     559syn keyword cConstant BedNorthGRDEnum
    554560syn keyword cConstant BedSlopeXEnum
    555561syn keyword cConstant BedSlopeYEnum
    … …  
    740746syn keyword cConstant SamplingKappaEnum
    741747syn keyword cConstant SealevelEnum
     748syn keyword cConstant SealevelGRDEnum
    742749syn keyword cConstant SealevelBarystaticMaskEnum
    743750syn keyword cConstant SealevelBarystaticOceanMaskEnum

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

    r26118 r26121  
    347347        SolidearthSettingsElasticEnum,
    348348        SealevelchangeGeometryDoneEnum,
     349        SealevelGrotm1Enum,
     350        SealevelGrotm2Enum,
     351        SealevelGrotm3Enum,
    349352        RotationalEquatorialMoiEnum,
    350353        TidalLoveHEnum,
    … …  
    546549        BaseSlopeYEnum,
    547550        BedEnum,
     551        BedGRDEnum,
    548552        BedEastEnum,
     553        BedEastGRDEnum,
    549554        BedNorthEnum,
     555        BedNorthGRDEnum,
    550556        BedSlopeXEnum,
    551557        BedSlopeYEnum,
    … …  
    736742        SamplingKappaEnum,
    737743        SealevelEnum,
     744        SealevelGRDEnum,
    738745        SealevelBarystaticMaskEnum,
    739746        SealevelBarystaticOceanMaskEnum,

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

    r26118 r26121  
    355355                case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic";
    356356                case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone";
     357                case SealevelGrotm1Enum : return "SealevelGrotm1";
     358                case SealevelGrotm2Enum : return "SealevelGrotm2";
     359                case SealevelGrotm3Enum : return "SealevelGrotm3";
    357360                case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi";
    358361                case TidalLoveHEnum : return "TidalLoveH";
    … …  
    552555                case BaseSlopeYEnum : return "BaseSlopeY";
    553556                case BedEnum : return "Bed";
     557                case BedGRDEnum : return "BedGRD";
    554558                case BedEastEnum : return "BedEast";
     559                case BedEastGRDEnum : return "BedEastGRD";
    555560                case BedNorthEnum : return "BedNorth";
     561                case BedNorthGRDEnum : return "BedNorthGRD";
    556562                case BedSlopeXEnum : return "BedSlopeX";
    557563                case BedSlopeYEnum : return "BedSlopeY";
    … …  
    742748                case SamplingKappaEnum : return "SamplingKappa";
    743749                case SealevelEnum : return "Sealevel";
     750                case SealevelGRDEnum : return "SealevelGRD";
    744751                case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask";
    745752                case SealevelBarystaticOceanMaskEnum : return "SealevelBarystaticOceanMask";

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

    r26118 r26121  
    361361              else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum;
    362362              else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum;
     363              else if (strcmp(name,"SealevelGrotm1")==0) return SealevelGrotm1Enum;
     364              else if (strcmp(name,"SealevelGrotm2")==0) return SealevelGrotm2Enum;
     365              else if (strcmp(name,"SealevelGrotm3")==0) return SealevelGrotm3Enum;
    363366              else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum;
    364367              else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum;
    … …  
    380383              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
    381384              else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum;
    382               else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
    383               else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
    384               else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
     388              if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum;
     389              else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum;
     390              else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum;
     391              else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum;
    389392              else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum;
    390393              else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum;
    … …  
    503506              else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
    504507              else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum;
    505               else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
    506               else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
    507               else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
     511              if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum;
     512              else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
     513              else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum;
     514              else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum;
    512515              else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
    513516              else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum;
    … …  
    564567              else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum;
    565568              else if (strcmp(name,"Bed")==0) return BedEnum;
     569              else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
    566570              else if (strcmp(name,"BedEast")==0) return BedEastEnum;
     571              else if (strcmp(name,"BedEastGRD")==0) return BedEastGRDEnum;
    567572              else if (strcmp(name,"BedNorth")==0) return BedNorthEnum;
     573              else if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum;
    568574              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    569575              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    … …  
    623629              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    624630              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    625               else if (strcmp(name,"EplHead")==0) return EplHeadEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"EplHead")==0) return EplHeadEnum;
    626635              else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum;
    627636              else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum;
    … …  
    629638              else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum;
    630639              else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
     640              else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum;
    635641              else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum;
    636642              else if (strcmp(name,"EsaRotationrate")==0) return EsaRotationrateEnum;
    … …  
    746752              else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum;
    747753              else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
    748               else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
    749758              else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
    750759              else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
    … …  
    752761              else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
    753762              else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
     763              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
    758764              else if (strcmp(name,"Sample")==0) return SampleEnum;
    759765              else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum;
    760766              else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum;
    761767              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
     768              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    762769              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    763770              else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum;
    … …  
    868875              else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum;
    869876              else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum;
    870               else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum;
    871881              else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum;
    872882              else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum;
    … …  
    875885              else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum;
    876886              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
     887              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    881888              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    882889              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    … …  
    991998              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    992999              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    993               else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    9941004              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    9951005              else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
    … …  
    9981008              else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum;
    9991009              else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
     1010              else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum;
    10041011              else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum;
    10051012              else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum;
    … …  
    11141121              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
    11151122              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
    1116               else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
    11171127              else if (strcmp(name,"Channel")==0) return ChannelEnum;
    11181128              else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
    … …  
    11211131              else if (strcmp(name,"Closed")==0) return ClosedEnum;
    11221132              else if (strcmp(name,"Colinear")==0) return ColinearEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
     1133              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    11271134              else if (strcmp(name,"Contact")==0) return ContactEnum;
    11281135              else if (strcmp(name,"Contour")==0) return ContourEnum;
    … …  
    12371244              else if (strcmp(name,"Internal")==0) return InternalEnum;
    12381245              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    1239               else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum;
    12401250              else if (strcmp(name,"J")==0) return JEnum;
    12411251              else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
    … …  
    12441254              else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum;
    12451255              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
     1256              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    12501257              else if (strcmp(name,"LambdaS")==0) return LambdaSEnum;
    12511258              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    … …  
    13601367              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    13611368              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
    1362               else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
    13631373              else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
    13641374              else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;
    … …  
    13671377              else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum;
    13681378              else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
     1379              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    13731380              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    13741381              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;

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

    r26110 r26121  
    143143                if testagainst2002,
    144144                        % TODO: Check if the following works as expected: 'pos' is empty, so nothing is assigned to 'md.masstransport.spcthickness(pos)'
    145                         md.masstransport.spcthickness=zeros(md.mesh.numberofelements,1);
     145                        md.masstransport.spcthickness=zeros(md.mesh.numberofvertices,1);
    146146                        %antarctica
    147147                        late=sum(md.mesh.lat(md.mesh.elements),2)/3;
    … …  
    150150                        ratio=0.225314032985172/0.193045366574523;
    151151                        %ratio=   1.276564103522540/.869956;
    152                         md.masstransport.spcthickness(pos)=-100*ratio;
     152                        md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100*ratio;
    153153                else
    154154                        in_fileID=fopen('../Data/AIS_delH_trend.txt', 'r');
    … …  
    168168                        northpole=find_point(md.mesh.long,md.mesh.lat,0,90);
    169169                        delHAIS(northpole)=0;
    170                         md.masstransport.spcthickness=mean(delHAIS(md.mesh.elements),2)/100;
     170                        md.masstransport.spcthickness=delHAIS/100;
    171171                end
    172172
    … …  
    249249        %}}}
    250250        %slc loading/calibration:  {{{
    251         md.masstransport.spcthickness=zeros(md.mesh.numberofelements,1);
     251        md.masstransport.spcthickness=zeros(md.mesh.numberofvertices,1);
    252252
    253253        if testagainst2002,
    … …  
    258258                pos=find(late > 70 &  late < 80 & longe>-60 & longe<-30);
    259259                ratio=.3823/.262344;
     260                md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100*ratio;
    260261                %md.masstransport.spcthickness(pos)=-100*ratio;
    261262
    … …  
    292293                delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3;
    293294
    294                 pos=find(delHGISe);
    295                 md.masstransport.spcthickness(pos)=delHGISe(pos)/100;
    296                 pos=find(delHGLAe);
    297                 md.masstransport.spcthickness(pos)=delHGLAe(pos)/100;
     295                pos=find(delHGIS);
     296                md.masstransport.spcthickness(pos)= md.masstransport.spcthickness(pos)-delHGIS(pos)/100;
     297                pos=find(delHGLA);
     298                md.masstransport.spcthickness(pos)= md.masstransport.spcthickness(pos)-delHGLA(pos)/100;
    298299
    299300                %adjust mask accordingly:
    300301                pos=find(md.masstransport.spcthickness);
    301                 flags=zeros(md.mesh.numberofvertices,1);
    302                 flags(md.mesh.elements(pos,:))=1;
    303                 pos=find(flags);
    304302                md.mask.ice_levelset(pos)=-1;
    305303                md.mask.ocean_levelset(pos)=1;
    … …  
    429427md.solidearth.settings.rotation=0;
    430428md.solidearth.requested_outputs= {'default',...
    431         'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate','SealevelchangeCumDeltathickness',...
    432         'SealevelNEsaRate', 'SealevelUEsaRate','SealevelStaticBarystaticMask','SealevelBarystaticOceanMask'};
     429        'DeltaIceThickness','Sealevel','SealevelUGrd',...
     430        'SealevelchangeBarystaticMask','SealevelchangeBarystaticOceanMask'};
    433431md=solve(md,'Transient');
    434432Seustatic=md.results.TransientSolution.Sealevel;
Note: See TracChangeset for help on using the changeset viewer.