Index: /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26121) @@ -184,4 +184,5 @@ } /*New optimized code:*/ + ToolkitsOptionsFromAnalysis(parameters,SealevelchangeAnalysisEnum); //this is requested by the BarystaticContributions class inner vectors. BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel); parameters->AddObject(new GenericParam(barystaticcontributions,BarystaticContributionsEnum)); Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26121) @@ -26,4 +26,8 @@ cumice=new Vector(nice); cumice->Set(0); cumice->Assemble(); } + else{ + ice=new Vector(1); + cumice=new Vector(1); + } iomodel->FetchData(&nhydro,"md.solidearth.nparthydro"); @@ -33,4 +37,8 @@ cumhydro=new Vector(nhydro); cumhydro->Set(0); cumhydro->Assemble(); } + else{ + hydro=new Vector(1); + cumhydro=new Vector(1); + } iomodel->FetchData(&nocean,"md.solidearth.npartocean"); if(nocean){ @@ -38,4 +46,8 @@ ocean=new Vector(nocean); cumocean=new Vector(nocean); cumocean->Set(0); cumocean->Assemble(); + } + else{ + ocean=new Vector(1); + cumocean=new Vector(1); } @@ -90,4 +102,33 @@ } /*}}}*/ +void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/ + + int id; + if(nice){ + id=reCast(pice[eid]); + ice->SetValue(id,icevalue,ADD_VAL); + } + else{ + ice->SetValue(0,icevalue,ADD_VAL); + } + + if(nhydro){ + id=reCast(phydro[eid]); + hydro->SetValue(id,hydrovalue,ADD_VAL); + } + else{ + hydro->SetValue(0,hydrovalue,ADD_VAL); + } + + if(nocean){ + id=reCast(pocean[eid]); + ocean->SetValue(id,oceanvalue,ADD_VAL); + } + else{ + ocean->SetValue(0,oceanvalue,ADD_VAL); + } + + +} /*}}}*/ void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/ @@ -118,32 +159,23 @@ results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time)); - cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;icumhydro->ToMPISerial0(); for (int i=0;icumocean->ToMPISerial0(); for (int i=0;iAddResult(new GenericExternalResult(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time)); - results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time)); - results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time)); + if(nice){ + cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;iAddResult(new GenericExternalResult(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time)); + } + if(nhydro){ + cumhydro_serial=this->cumhydro->ToMPISerial0(); for (int i=0;iAddResult(new GenericExternalResult(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time)); + } + if(nocean){ + cumocean_serial=this->cumocean->ToMPISerial0(); for (int i=0;iAddResult(new GenericExternalResult(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time)); + } if(IssmComm::GetRank()==0){ - xDelete(cumice_serial); - xDelete(cumhydro_serial); - xDelete(cumocean_serial); + if(nice)xDelete(cumice_serial); + if(nhydro)xDelete(cumhydro_serial); + if(nocean)xDelete(cumocean_serial); } return; } /*}}}*/ -void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/ - - int id; - - id=reCast(pice[eid]); - ice->SetValue(id,icevalue,ADD_VAL); - - id=reCast(phydro[eid]); - hydro->SetValue(id,hydrovalue,ADD_VAL); - - id=reCast(pocean[eid]); - ocean->SetValue(id,oceanvalue,ADD_VAL); - -} /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/BarystaticContributions.h =================================================================== --- /issm/trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26121) @@ -17,19 +17,19 @@ public: - Vector* ice; //contributions to every ice partition + Vector* ice; //contributions to every ice partition (size nice x 1) Vector* cumice; //cumulated contributions to every ice partition int nice; //number of ice partitions - IssmDouble* pice; //ice partition + IssmDouble* pice; //ice partition (nel) - Vector* hydro; //contributions to every hydro partition + Vector* hydro; //contributions to every hydro partition (size nhydro x 1) Vector* cumhydro; //cumulated contributions to every hydro partition int nhydro; //number of hydro partitions - IssmDouble* phydro; //hydro partition + IssmDouble* phydro; //hydro partition (nel) - Vector* ocean; //contributions to every ocean partition + Vector* ocean; //contributions to every ocean partition (size nocean x 1) Vector* cumocean; //cumulated contributions to every ocean partition int nocean; //number of ocean partitions - IssmDouble* pocean; //ocean partition - + IssmDouble* pocean; //ocean partition (nel) + /*BarystaticContributions constructors, destructors :*/ BarystaticContributions(IoModel* iomodel ); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 26121) @@ -382,4 +382,5 @@ virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0; virtual void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; + virtual void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads)=0; virtual IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution,IssmDouble* partition,IssmDouble oceanarea)=0; virtual IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks,Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea)=0; @@ -393,6 +394,7 @@ virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0; virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; - virtual void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0; - virtual void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0; + virtual void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; + virtual void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; + virtual void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks)=0; virtual void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 26121) @@ -219,4 +219,5 @@ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; @@ -230,6 +231,7 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector,SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 26121) @@ -172,4 +172,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, @@ -185,6 +186,7 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26121) @@ -180,4 +180,5 @@ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; + void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){_error_("not implemented yet!");}; void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet!");}; @@ -191,6 +192,7 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){_error_("not implemented yet");}; void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26121) @@ -6252,5 +6252,6 @@ IssmDouble constant; IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim; - int sidlist[NUMVERTICES]; + int sidlist[NUMVERTICES]; + int sid; #ifdef _HAVE_RESTRICT_ @@ -6277,5 +6278,5 @@ IssmDouble* longes=NULL; #endif - + /*elastic green function:*/ int index; @@ -6285,5 +6286,18 @@ bool computerigid = false; bool computeelastic = false; + bool computerotation = false; int horiz; + + /*Rotational:*/ + IssmDouble* tide_love_h = NULL; + IssmDouble* tide_love_k = NULL; + IssmDouble* load_love_k = NULL; + IssmDouble tide_love_k2secular; + IssmDouble moi_e, moi_p, omega, g; + IssmDouble Grotm1[3]; + IssmDouble Grotm2[3]; + IssmDouble Grotm3[3]; + IssmDouble pre; + /*recover parameters: */ @@ -6291,8 +6305,19 @@ this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); + this->parameters->FindParam(&computerotation,SolidearthSettingsRotationEnum); this->parameters->FindParam(&nel,MeshNumberofelementsEnum); this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); this->parameters->FindParam(&planetradius,SolidearthPlanetRadiusEnum); this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + + if(computerotation){ + parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum); + parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum); + parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum); + parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum); + parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum); + parameters->FindParam(&moi_p,RotationalPolarMoiEnum); + parameters->FindParam(&omega,RotationalAngularVelocityEnum); + } /*compute area and add to inputs:*/ @@ -6349,5 +6374,5 @@ IssmDouble late=lates[e]; IssmDouble longe=longes[e]; - int sid=sidlist[i]; + sid=sidlist[i]; /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ @@ -6387,4 +6412,21 @@ } } + if(computerotation){ + + for(int i=0;i<3;i++){ + sid=sidlist[i]; + lati=latitude[sid]/180.*M_PI; + longi=longitude[sid]/180.*M_PI; + + pre=((1.0+tide_love_k[2]-tide_love_h[2])/9.81)*pow(omega*planetradius,2.0); + Grotm1[i]= - pre* 0.5*sin(2.*lati)*cos(longi); + Grotm2[i]= - pre* 0.5*sin(2.*lati)*sin(longi); + Grotm3[i]= - pre* (1/6.0 - 0.5*cos(2.0*lati)); + } + this->AddInput(SealevelGrotm1Enum,&Grotm1[0],P1Enum); + this->AddInput(SealevelGrotm2Enum,&Grotm2[0],P1Enum); + this->AddInput(SealevelGrotm3Enum,&Grotm3[0],P1Enum); + } + /*Add in inputs:*/ @@ -6491,7 +6533,7 @@ } - /*if we are an ice shelf, are we fully grounded or not? (used later):*/ if (masks->notfullygrounded[this->lid]) notfullygrounded=true; //used later on. + /*recover some parameters:*/ @@ -6606,5 +6648,5 @@ } /*}}}*/ -void Tria::SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ /*sal green function:*/ @@ -6612,11 +6654,27 @@ IssmDouble SealevelGRD[NUMVERTICES]={0,0,0}; IssmDouble oceanaverage,oceanarea=0; - + IssmDouble rho_water; + bool sal = false; + bool rotation= false; int size; int nel; + IssmDouble Grotm1[3]; + IssmDouble Grotm2[3]; + IssmDouble Grotm3[3]; this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); this->parameters->FindParam(&nel,MeshNumberofelementsEnum); + this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + + if(rotation){ + Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); + Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); + Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); + + for(int i=0;iinputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); @@ -6633,33 +6691,121 @@ /*add ocean average in the global sealevelloads vector:*/ - sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL); + sealevelloads->SetValue(this->sid,oceanaverage*rho_water,INS_VAL); oceanareas->SetValue(this->sid,oceanarea,INS_VAL); return; } /*}}}*/ -void Tria::SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/ - - bool converged=false; +void Tria::SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ + IssmDouble SealevelGRD[3]={0,0,0}; IssmDouble oceanaverage,oceanarea=0; + IssmDouble rho_water; int nel; bool sal = false; IssmDouble* G=NULL; int size; - + IssmDouble Grotm1[3]; + IssmDouble Grotm2[3]; + IssmDouble Grotm3[3]; + bool rotation= false; + this->parameters->FindParam(&nel,MeshNumberofelementsEnum); this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + this->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + + if(rotation){ + Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); + Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); + Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); + + for(int i=0;iinputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); + + for(int i=0;iSetValue(this->sid,oceanaverage*rho_water,INS_VAL); + +} /*}}}*/ +void Tria::SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks){ /*{{{*/ + + IssmDouble Sealevel[3]={0,0,0}; + IssmDouble SealevelRSL[3]={0,0,0}; + IssmDouble SealevelU[3]={0,0,0}; + IssmDouble SealevelN[3]={0,0,0}; + IssmDouble SealevelE[3]={0,0,0}; + int nel; + bool sal = false; + IssmDouble* G=NULL; + IssmDouble* GU=NULL; + IssmDouble* GE=NULL; + IssmDouble* GN=NULL; + int horiz; + int size; + IssmDouble Grotm1[3]; + IssmDouble Grotm2[3]; + IssmDouble Grotm3[3]; + bool rotation= false; + bool elastic=false; + + this->parameters->FindParam(&nel,MeshNumberofelementsEnum); + this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum); + this->parameters->FindParam(&elastic,SolidearthSettingsElasticEnum); + this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + + if(rotation){ + Element::GetInputListOnVertices(&Grotm1[0],SealevelGrotm1Enum); + Element::GetInputListOnVertices(&Grotm2[0],SealevelGrotm2Enum); + Element::GetInputListOnVertices(&Grotm3[0],SealevelGrotm3Enum); + + for(int i=0;iinputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); + if(elastic) this->inputs->GetArrayPtr(SealevelchangeGUEnum,this->lid,&GU,&size); + if(horiz && elastic){ + this->inputs->GetArrayPtr(SealevelchangeGEEnum,this->lid,&GE,&size); + this->inputs->GetArrayPtr(SealevelchangeGNEnum,this->lid,&GN,&size); + } for(int i=0;iSetValue(this->sid,oceanaverage,INS_VAL); + SealevelRSL[i]+=G[i*nel+e]*(sealevelloads[e]+loads[e]); + } + if(elastic){ + for (int e=0;eAddInput(SealevelGRDEnum,Sealevel,P1Enum); + this->AddInput(BedGRDEnum,SealevelU,P1Enum); + if(horiz){ + this->AddInput(BedEastGRDEnum,SealevelE,P1Enum); + this->AddInput(BedNorthGRDEnum,SealevelN,P1Enum); + } + } /*}}}*/ @@ -6733,4 +6879,71 @@ } /*}}}*/ +void Tria::SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads){/*{{{*/ + + IssmDouble S; + + /*Compute area of element:*/ + IssmDouble area,planetarea; + this->Element::GetInputValue(&area,AreaEnum); + + /*recover earth area: */ + this->parameters->FindParam(&planetarea,SolidearthPlanetAreaEnum); + + /*Compute lat,long,radius of elemental centroid: */ + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble late,longe,re; + /* Where is the centroid of this element?:{{{*/ + ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); + + IssmDouble minlong=400; + IssmDouble maxlong=-20; + for (int i=0;imaxlong)maxlong=llr_list[i][1]; + if(llr_list[i][1]180){ + if (llr_list[0][1]==0)llr_list[0][1]=360; + if (llr_list[1][1]==0)llr_list[1][1]=360; + if (llr_list[2][1]==0)llr_list[2][1]=360; + } + + // correction at the north pole + if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + //correction at the south pole + if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0; + if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0; + if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0; + + late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0; + longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0; + + late=90.-late; + if(longe>180.)longe=(longe-180.)-180.; + + late=late/180.*M_PI; + longe=longe/180.*M_PI; + /*}}}*/ + re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0; + + + /*recover total load: */ + S=sealevelloads[this->Sid()]+loads[this->Sid()]; + + /* Perturbation terms for moment of inertia (moi_list): + * computed analytically (see Wu & Peltier, eqs 10 & 32) + * also consistent with my GMD formulation! + * ALL in geographic coordinates + * */ + dI_list[0] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/planetarea; + dI_list[1] = -4*M_PI*(S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/planetarea; + dI_list[2] = +4*M_PI*(S*area)*pow(re,4)*(1-pow(sin(late),2))/planetarea; + return; +}/*}}}*/ #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 26121) @@ -163,4 +163,5 @@ #ifdef _HAVE_SEALEVELCHANGE_ void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); + void SealevelchangeMomentOfInertiaElement(IssmDouble* dI_list, IssmDouble* loads, IssmDouble* sealevelloads); void SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); @@ -175,6 +176,7 @@ void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks); - void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks); - void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks); + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); + void SealevelchangeDeformationConvolution(IssmDouble* sealevelloads, IssmDouble* loads, IssmDouble* rotationvector, SealevelMasks* masks); void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks); #endif Index: /issm/trunk-jpl/src/c/classes/Inputs/Input.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Input.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Inputs/Input.h (revision 26121) @@ -44,4 +44,5 @@ virtual void AXPY(Input* xinput,IssmDouble scalar){_error_("Not implemented yet");}; + virtual void Shift(IssmDouble scalar){_error_("Not implemented yet");}; virtual void PointWiseMult(Input* xinput){_error_("Not implemented yet");}; virtual void Pow(IssmDouble scale_factor){_error_("Not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Inputs/Inputs.cpp (revision 26121) @@ -247,5 +247,5 @@ } /*}}}*/ -void Inputs::AXPY(IssmDouble alpha, int xenum, int yenum, int zenum){/*{{{*/ +void Inputs::ZAXPY(IssmDouble alpha, int xenum, int yenum, int zenum){/*{{{*/ _assert_(this); @@ -268,4 +268,34 @@ /*AXPY: */ this->inputs[index_z]->AXPY(this->inputs[index_x],alpha); +} +/*}}}*/ +void Inputs::AXPY(IssmDouble alpha, int xenum, int yenum ){/*{{{*/ + + _assert_(this); + + /*Get indices from enums*/ + int index_x = EnumToIndex(xenum); + int index_y = EnumToIndex(yenum); + + /*Make sure that old one exists*/ + if(!this->inputs[index_x]) _error_("Input "<inputs[index_y]) _error_("Input "<inputs[index_y]->AXPY(this->inputs[index_x],alpha); +} +/*}}}*/ +void Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/ + + _assert_(this); + + /*Get indices from enums*/ + int index_x = EnumToIndex(xenum); + + /*Make sure that x exists*/ + if(!this->inputs[index_x]) _error_("Input "<inputs[index_x]->Shift(alpha); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Inputs/Inputs.h (revision 26121) @@ -46,5 +46,7 @@ int DeleteInput(int enum_type); void DuplicateInput(int original_enum,int new_enum); - void AXPY(IssmDouble alpha, int xenum, int yenum, int zenum); + void ZAXPY(IssmDouble alpha, int xenum, int yenum, int zenum); + void AXPY(IssmDouble alpha, int xenum, int yenum); + void Shift(int inputenum, IssmDouble alpha); void DeepEcho(void); void DeepEcho(int enum_in); Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp (revision 26121) @@ -384,4 +384,11 @@ } /*}}}*/ +void TriaInput::Shift(IssmDouble alpha){/*{{{*/ + + /*Carry out the shift operation:*/ + for(int i=0;iM*this->N;i++) this->values[i] +=alpha; + for(int i=0;iinterpolation);i++) this->element_values[i] += alpha; +} +/*}}}*/ void TriaInput::PointWiseMult(Input* xinput){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 26120) +++ /issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h (revision 26121) @@ -40,4 +40,5 @@ void Pow(IssmDouble scalar); void AXPY(Input* xinput,IssmDouble scalar); + void Shift(IssmDouble scalar); void PointWiseMult(Input* xinput); void Serve(int numindices,int* indices); Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 26121) @@ -246,11 +246,10 @@ /*From old and new thickness, create delta thickness and accumulate:*/ - femmodel->inputs->AXPY(-1, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum); - femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum,DummyEnum); - femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaTwsEnum); + femmodel->inputs->ZAXPY(-1, WaterColumnOldEnum,WatercolumnEnum,DeltaTwsEnum); + femmodel->inputs->AXPY(+1, DeltaTwsEnum,AccumulatedDeltaTwsEnum); /*compute total water column change between two sea-level solver time steps, ie. every frequency*dt:*/ if(count==frequency){ - femmodel->inputs->AXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum); + femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaTwsEnum,AccumulatedDeltaTwsEnum,DeltaTwsEnum); femmodel->inputs->DuplicateInput(AccumulatedDeltaTwsEnum,OldAccumulatedDeltaTwsEnum); } Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 26121) @@ -117,7 +117,6 @@ /*From old and new thickness, create delta thickness and accumulate:*/ - femmodel->inputs->AXPY(-1, ThicknessOldEnum,ThicknessEnum,DeltaIceThicknessEnum); - femmodel->inputs->AXPY(+1, DeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DummyEnum); - femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaIceThicknessEnum); + femmodel->inputs->ZAXPY(-1, ThicknessOldEnum,ThicknessEnum,DeltaIceThicknessEnum); + femmodel->inputs->AXPY(+1, DeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum); /*for Ivins deformation model, keep history of ice thickness changes inside TransientAccumulatedDeltaIceThicknessEnum:*/ @@ -154,5 +153,5 @@ /*compute total ice thickness change between two sea-level solver time steps, ie. every frequency*dt:*/ if(count==frequency){ - femmodel->inputs->AXPY(-1, OldAccumulatedDeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DeltaIceThicknessEnum); + femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaIceThicknessEnum,AccumulatedDeltaIceThicknessEnum,DeltaIceThicknessEnum); femmodel->inputs->DuplicateInput(AccumulatedDeltaIceThicknessEnum,OldAccumulatedDeltaIceThicknessEnum); } Index: /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/cores/oceantransport_core.cpp (revision 26121) @@ -75,14 +75,13 @@ /* From old and new bottom pressures, create delta bottom pressure, delta dsl and delta str. * Accumulate delta bottom pressure: */ - femmodel->inputs->AXPY(-1, BottomPressureOldEnum,BottomPressureEnum,DeltaBottomPressureEnum); - femmodel->inputs->AXPY(+1, DeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DummyEnum); - femmodel->inputs->DuplicateInput(DummyEnum,AccumulatedDeltaBottomPressureEnum); + femmodel->inputs->ZAXPY(-1, BottomPressureOldEnum,BottomPressureEnum,DeltaBottomPressureEnum); + femmodel->inputs->AXPY(+1, DeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum); - femmodel->inputs->AXPY(-1, DslOldEnum,DslEnum,DeltaDslEnum); - femmodel->inputs->AXPY(-1, StrOldEnum,StrEnum,DeltaStrEnum); + femmodel->inputs->ZAXPY(-1, DslOldEnum,DslEnum,DeltaDslEnum); + femmodel->inputs->ZAXPY(-1, StrOldEnum,StrEnum,DeltaStrEnum); /* Compute total bottom pressure change between two sea-level solver time steps, ie. every frequency*dt. */ if(count==frequency){ - femmodel->inputs->AXPY(-1, OldAccumulatedDeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DeltaBottomPressureEnum); + femmodel->inputs->ZAXPY(-1, OldAccumulatedDeltaBottomPressureEnum,AccumulatedDeltaBottomPressureEnum,DeltaBottomPressureEnum); femmodel->inputs->DuplicateInput(AccumulatedDeltaBottomPressureEnum,OldAccumulatedDeltaBottomPressureEnum); } Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26121) @@ -18,4 +18,5 @@ IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea); SealevelMasks* sealevel_masks(FemModel* femmodel); +void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads); /*}}}*/ @@ -398,4 +399,5 @@ SealevelMasks* masks=NULL; GenericParam* barycontribparam=NULL; + IssmDouble rotationaxismotionvector[3]; Vector* loads=NULL; @@ -413,4 +415,5 @@ int step; IssmDouble time; + bool converged=false; int modelid,earthid; @@ -418,4 +421,5 @@ int count,frequency,iscoupling; int grd=0; + int computesealevel=0; /*}}}*/ @@ -427,4 +431,5 @@ femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum); femmodel->parameters->FindParam(&count,SealevelchangeRunCountEnum); + femmodel->parameters->FindParam(&computesealevel,SolidearthSettingsComputesealevelchangeEnum); /*}}}*/ @@ -463,15 +468,25 @@ element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); } + loads->Assemble(); //broadcast loads - loads->Assemble(); allloads=loads->ToMPISerial(); + allloads=loads->ToMPISerial(); + + //compute rotation axis motion: + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,NULL); + + /*skip computation of sea level if requested, which means sea level loads should be zeroed */ + if(!computesealevel){ + allsealevelloads=xNewZeroInit(nel); + goto deformation; + } /*convolve loads:*/ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks); + element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,rotationaxismotionvector,masks); } sealevelloads->Assemble(); - + //Get ocean area: oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); @@ -479,10 +494,12 @@ //substract ocean average and barystatic contributionfrom sea level loads: - sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea)); + sealevelloads->Shift(barycontrib->Total()/oceanarea - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea)); //broadcast sea level loads allsealevelloads=sealevelloads->ToMPISerial(); - bool converged=false; + //compute rotation axis motion: + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); + for(;;){ @@ -492,15 +509,17 @@ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks); + element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,rotationaxismotionvector,masks); } sealevelloads->Assemble(); - //substract ocean average and barystatic contribution - sealevelloads->Shift(barycontrib->Total()/oceanarea/rho_water - SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea)); + sealevelloads->Shift(barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea)); //broadcast sea level loads allsealevelloads=sealevelloads->ToMPISerial(); + //compute rotation axis motion: + RotationAxisMotion(&rotationaxismotionvector[0],femmodel,allloads,allsealevelloads); + //convergence? slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs); @@ -508,22 +527,27 @@ } + deformation: + /*convolve loads and sea level loads to get the deformation:*/ for(Object* & object : femmodel->elements->objects){ Element* element = xDynamicCast(object); - //element->SealevelchangeDeformationConvolution(allsealevelloads, allloads,masks); + element->SealevelchangeDeformationConvolution(allsealevelloads, allloads, rotationaxismotionvector,masks); } /*Update bedrock motion and geoid:*/ - /*inputs->AXPY(SealevelEnum,SealevelGRD,1); - inputs->AXPY(BedEnum,BedGRD,1); + if(computesealevel){ + femmodel->inputs->Shift(SealevelGRDEnum,barycontrib->Total()/oceanarea- SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea)); + + /*cumulate barystatic contributions and save to results: */ + barycontrib->Cumulate(femmodel->parameters); + barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); + } + + femmodel->inputs->AXPY(SealevelEnum,SealevelGRDEnum,1); + femmodel->inputs->AXPY(BedEnum,BedGRDEnum,1); if(horiz){ - inputs->AXPY(BedrockeastEnum,BedEastGRD,1); - inputs->AXPY(BedrocknorthEnum,BedNorthGRD,1); - }*/ - - /*cumulate barystatic contributions and save to results: */ - barycontrib->Cumulate(femmodel->parameters); - barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); - + femmodel->inputs->AXPY(BedEastEnum,BedEastGRDEnum,1); + femmodel->inputs->AXPY(BedNorthEnum,BedNorthGRDEnum,1); + } } @@ -1216,2 +1240,48 @@ return sealevelloadsaverage/oceanarea; } /*}}}*/ +void RotationAxisMotion(IssmDouble* m, FemModel* femmodel,IssmDouble* loads, IssmDouble* sealeveloads){ /*{{{*/ + + IssmDouble moi_list[3]={0,0,0}; + IssmDouble moi_list_cpu[3]={0,0,0}; + IssmDouble* tide_love_h = NULL; + IssmDouble* tide_love_k = NULL; + IssmDouble* load_love_k = NULL; + IssmDouble tide_love_k2secular; + IssmDouble moi_e, moi_p; + IssmDouble m1, m2, m3; + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&load_love_k,NULL,NULL,LoadLoveKEnum); + femmodel->parameters->FindParam(&tide_love_h,NULL,NULL,TidalLoveHEnum); + femmodel->parameters->FindParam(&tide_love_k,NULL,NULL,TidalLoveKEnum); + femmodel->parameters->FindParam(&tide_love_k2secular,TidalLoveK2SecularEnum); + femmodel->parameters->FindParam(&moi_e,RotationalEquatorialMoiEnum); + femmodel->parameters->FindParam(&moi_p,RotationalPolarMoiEnum); + + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeMomentOfInertiaElement(&moi_list[0],loads,sealeveloads); + moi_list_cpu[0] += moi_list[0]; + moi_list_cpu[1] += moi_list[1]; + moi_list_cpu[2] += moi_list[2]; + } + + ISSM_MPI_Reduce (&moi_list_cpu[0],&moi_list[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + ISSM_MPI_Reduce (&moi_list_cpu[1],&moi_list[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + ISSM_MPI_Reduce (&moi_list_cpu[2],&moi_list[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&moi_list[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + + /*compute perturbation terms for angular velocity vector: */ + m1 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[0]; + m2 = 1/(1-tide_love_k[2]/tide_love_k2secular) * (1+load_love_k[2])/(moi_p-moi_e) * moi_list[1]; + m3 = -(1+load_love_k[2])/moi_p * moi_list[2]; // term associated with fluid number (3-order-of-magnitude smaller) is negelected + + /*Assign output pointers:*/ + m[0]=m1; + m[1]=m2; + m[2]=m3; +} /*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26120) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26121) @@ -353,4 +353,7 @@ syn keyword cConstant SolidearthSettingsElasticEnum syn keyword cConstant SealevelchangeGeometryDoneEnum +syn keyword cConstant SealevelGrotm1Enum +syn keyword cConstant SealevelGrotm2Enum +syn keyword cConstant SealevelGrotm3Enum syn keyword cConstant RotationalEquatorialMoiEnum syn keyword cConstant TidalLoveHEnum @@ -550,6 +553,9 @@ syn keyword cConstant BaseSlopeYEnum syn keyword cConstant BedEnum +syn keyword cConstant BedGRDEnum syn keyword cConstant BedEastEnum +syn keyword cConstant BedEastGRDEnum syn keyword cConstant BedNorthEnum +syn keyword cConstant BedNorthGRDEnum syn keyword cConstant BedSlopeXEnum syn keyword cConstant BedSlopeYEnum @@ -740,4 +746,5 @@ syn keyword cConstant SamplingKappaEnum syn keyword cConstant SealevelEnum +syn keyword cConstant SealevelGRDEnum syn keyword cConstant SealevelBarystaticMaskEnum syn keyword cConstant SealevelBarystaticOceanMaskEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26120) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26121) @@ -347,4 +347,7 @@ SolidearthSettingsElasticEnum, SealevelchangeGeometryDoneEnum, + SealevelGrotm1Enum, + SealevelGrotm2Enum, + SealevelGrotm3Enum, RotationalEquatorialMoiEnum, TidalLoveHEnum, @@ -546,6 +549,9 @@ BaseSlopeYEnum, BedEnum, + BedGRDEnum, BedEastEnum, + BedEastGRDEnum, BedNorthEnum, + BedNorthGRDEnum, BedSlopeXEnum, BedSlopeYEnum, @@ -736,4 +742,5 @@ SamplingKappaEnum, SealevelEnum, + SealevelGRDEnum, SealevelBarystaticMaskEnum, SealevelBarystaticOceanMaskEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26121) @@ -355,4 +355,7 @@ case SolidearthSettingsElasticEnum : return "SolidearthSettingsElastic"; case SealevelchangeGeometryDoneEnum : return "SealevelchangeGeometryDone"; + case SealevelGrotm1Enum : return "SealevelGrotm1"; + case SealevelGrotm2Enum : return "SealevelGrotm2"; + case SealevelGrotm3Enum : return "SealevelGrotm3"; case RotationalEquatorialMoiEnum : return "RotationalEquatorialMoi"; case TidalLoveHEnum : return "TidalLoveH"; @@ -552,6 +555,9 @@ case BaseSlopeYEnum : return "BaseSlopeY"; case BedEnum : return "Bed"; + case BedGRDEnum : return "BedGRD"; case BedEastEnum : return "BedEast"; + case BedEastGRDEnum : return "BedEastGRD"; case BedNorthEnum : return "BedNorth"; + case BedNorthGRDEnum : return "BedNorthGRD"; case BedSlopeXEnum : return "BedSlopeX"; case BedSlopeYEnum : return "BedSlopeY"; @@ -742,4 +748,5 @@ case SamplingKappaEnum : return "SamplingKappa"; case SealevelEnum : return "Sealevel"; + case SealevelGRDEnum : return "SealevelGRD"; case SealevelBarystaticMaskEnum : return "SealevelBarystaticMask"; case SealevelBarystaticOceanMaskEnum : return "SealevelBarystaticOceanMask"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26120) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26121) @@ -361,4 +361,7 @@ else if (strcmp(name,"SolidearthSettingsElastic")==0) return SolidearthSettingsElasticEnum; else if (strcmp(name,"SealevelchangeGeometryDone")==0) return SealevelchangeGeometryDoneEnum; + else if (strcmp(name,"SealevelGrotm1")==0) return SealevelGrotm1Enum; + else if (strcmp(name,"SealevelGrotm2")==0) return SealevelGrotm2Enum; + else if (strcmp(name,"SealevelGrotm3")==0) return SealevelGrotm3Enum; else if (strcmp(name,"RotationalEquatorialMoi")==0) return RotationalEquatorialMoiEnum; else if (strcmp(name,"TidalLoveH")==0) return TidalLoveHEnum; @@ -380,11 +383,11 @@ else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum; else if (strcmp(name,"RotationalPolarMoi")==0) return RotationalPolarMoiEnum; - else if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; - else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; - else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; + if (strcmp(name,"SolidearthSettingsReltol")==0) return SolidearthSettingsReltolEnum; + else if (strcmp(name,"SealevelchangeRequestedOutputs")==0) return SealevelchangeRequestedOutputsEnum; + else if (strcmp(name,"SolidearthSettingsRigid")==0) return SolidearthSettingsRigidEnum; + else if (strcmp(name,"SolidearthSettingsRotation")==0) return SolidearthSettingsRotationEnum; else if (strcmp(name,"SealevelchangeRunCount")==0) return SealevelchangeRunCountEnum; else if (strcmp(name,"SealevelchangeTransitions")==0) return SealevelchangeTransitionsEnum; @@ -503,11 +506,11 @@ else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum; else if (strcmp(name,"TransientIsgia")==0) return TransientIsgiaEnum; - else if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; - else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; - else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; + if (strcmp(name,"TransientIsgroundingline")==0) return TransientIsgroundinglineEnum; + else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum; + else if (strcmp(name,"TransientIsmasstransport")==0) return TransientIsmasstransportEnum; + else if (strcmp(name,"TransientIsoceantransport")==0) return TransientIsoceantransportEnum; else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum; else if (strcmp(name,"TransientIsoceancoupling")==0) return TransientIsoceancouplingEnum; @@ -564,6 +567,9 @@ else if (strcmp(name,"BaseSlopeY")==0) return BaseSlopeYEnum; else if (strcmp(name,"Bed")==0) return BedEnum; + else if (strcmp(name,"BedGRD")==0) return BedGRDEnum; else if (strcmp(name,"BedEast")==0) return BedEastEnum; + else if (strcmp(name,"BedEastGRD")==0) return BedEastGRDEnum; else if (strcmp(name,"BedNorth")==0) return BedNorthEnum; + else if (strcmp(name,"BedNorthGRD")==0) return BedNorthGRDEnum; else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum; else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum; @@ -623,5 +629,8 @@ else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; - else if (strcmp(name,"EplHead")==0) return EplHeadEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"EplHead")==0) return EplHeadEnum; else if (strcmp(name,"EplHeadOld")==0) return EplHeadOldEnum; else if (strcmp(name,"EplHeadSlopeX")==0) return EplHeadSlopeXEnum; @@ -629,8 +638,5 @@ else if (strcmp(name,"EplHeadSubstep")==0) return EplHeadSubstepEnum; else if (strcmp(name,"EplHeadTransient")==0) return EplHeadTransientEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; + else if (strcmp(name,"EsaEmotion")==0) return EsaEmotionEnum; else if (strcmp(name,"EsaNmotion")==0) return EsaNmotionEnum; else if (strcmp(name,"EsaRotationrate")==0) return EsaRotationrateEnum; @@ -746,5 +752,8 @@ else if (strcmp(name,"RadarAttenuationMacGregor")==0) return RadarAttenuationMacGregorEnum; else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum; - else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum; else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum; else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum; @@ -752,12 +761,10 @@ else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum; else if (strcmp(name,"RheologyBInitialguessMisfit")==0) return RheologyBInitialguessMisfitEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; + else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; else if (strcmp(name,"Sample")==0) return SampleEnum; else if (strcmp(name,"SamplingBeta")==0) return SamplingBetaEnum; else if (strcmp(name,"SamplingKappa")==0) return SamplingKappaEnum; else if (strcmp(name,"Sealevel")==0) return SealevelEnum; + else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum; else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum; else if (strcmp(name,"SealevelBarystaticOceanMask")==0) return SealevelBarystaticOceanMaskEnum; @@ -868,5 +875,8 @@ else if (strcmp(name,"SmbRunoff")==0) return SmbRunoffEnum; else if (strcmp(name,"SmbRunoffSubstep")==0) return SmbRunoffSubstepEnum; - else if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SmbRunoffTransient")==0) return SmbRunoffTransientEnum; else if (strcmp(name,"SmbS0gcm")==0) return SmbS0gcmEnum; else if (strcmp(name,"SmbS0p")==0) return SmbS0pEnum; @@ -875,8 +885,5 @@ else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; + else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; else if (strcmp(name,"SmbT")==0) return SmbTEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; @@ -991,5 +998,8 @@ else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum; else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum; - else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum; else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum; else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum; @@ -998,8 +1008,5 @@ else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; + else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; @@ -1114,5 +1121,8 @@ else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum; else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum; - else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum; else if (strcmp(name,"Channel")==0) return ChannelEnum; else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum; @@ -1121,8 +1131,5 @@ else if (strcmp(name,"Closed")==0) return ClosedEnum; else if (strcmp(name,"Colinear")==0) return ColinearEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Constraints")==0) return ConstraintsEnum; + else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; else if (strcmp(name,"Contact")==0) return ContactEnum; else if (strcmp(name,"Contour")==0) return ContourEnum; @@ -1237,5 +1244,8 @@ else if (strcmp(name,"Internal")==0) return InternalEnum; else if (strcmp(name,"Intersect")==0) return IntersectEnum; - else if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"InversionVzObs")==0) return InversionVzObsEnum; else if (strcmp(name,"J")==0) return JEnum; else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum; @@ -1244,8 +1254,5 @@ else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; + else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; @@ -1360,5 +1367,8 @@ else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; - else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum; else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; @@ -1367,8 +1377,5 @@ else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; + else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 26120) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 26121) @@ -143,5 +143,5 @@ if testagainst2002, % TODO: Check if the following works as expected: 'pos' is empty, so nothing is assigned to 'md.masstransport.spcthickness(pos)' - md.masstransport.spcthickness=zeros(md.mesh.numberofelements,1); + md.masstransport.spcthickness=zeros(md.mesh.numberofvertices,1); %antarctica late=sum(md.mesh.lat(md.mesh.elements),2)/3; @@ -150,5 +150,5 @@ ratio=0.225314032985172/0.193045366574523; %ratio= 1.276564103522540/.869956; - md.masstransport.spcthickness(pos)=-100*ratio; + md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100*ratio; else in_fileID=fopen('../Data/AIS_delH_trend.txt', 'r'); @@ -168,5 +168,5 @@ northpole=find_point(md.mesh.long,md.mesh.lat,0,90); delHAIS(northpole)=0; - md.masstransport.spcthickness=mean(delHAIS(md.mesh.elements),2)/100; + md.masstransport.spcthickness=delHAIS/100; end @@ -249,5 +249,5 @@ %}}} %slc loading/calibration: {{{ - md.masstransport.spcthickness=zeros(md.mesh.numberofelements,1); + md.masstransport.spcthickness=zeros(md.mesh.numberofvertices,1); if testagainst2002, @@ -258,4 +258,5 @@ pos=find(late > 70 & late < 80 & longe>-60 & longe<-30); ratio=.3823/.262344; + md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100*ratio; %md.masstransport.spcthickness(pos)=-100*ratio; @@ -292,14 +293,11 @@ delHGLAe=delHGLA(md.mesh.elements)*[1;1;1]/3; - pos=find(delHGISe); - md.masstransport.spcthickness(pos)=delHGISe(pos)/100; - pos=find(delHGLAe); - md.masstransport.spcthickness(pos)=delHGLAe(pos)/100; + pos=find(delHGIS); + md.masstransport.spcthickness(pos)= md.masstransport.spcthickness(pos)-delHGIS(pos)/100; + pos=find(delHGLA); + md.masstransport.spcthickness(pos)= md.masstransport.spcthickness(pos)-delHGLA(pos)/100; %adjust mask accordingly: pos=find(md.masstransport.spcthickness); - flags=zeros(md.mesh.numberofvertices,1); - flags(md.mesh.elements(pos,:))=1; - pos=find(flags); md.mask.ice_levelset(pos)=-1; md.mask.ocean_levelset(pos)=1; @@ -429,6 +427,6 @@ md.solidearth.settings.rotation=0; md.solidearth.requested_outputs= {'default',... - 'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate','SealevelchangeCumDeltathickness',... - 'SealevelNEsaRate', 'SealevelUEsaRate','SealevelStaticBarystaticMask','SealevelBarystaticOceanMask'}; + 'DeltaIceThickness','Sealevel','SealevelUGrd',... + 'SealevelchangeBarystaticMask','SealevelchangeBarystaticOceanMask'}; md=solve(md,'Transient'); Seustatic=md.results.TransientSolution.Sealevel;