Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24924) @@ -378,4 +378,5 @@ virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0; virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0; + virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0; virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24924) @@ -215,4 +215,5 @@ IssmDouble OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");}; void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24924) @@ -173,4 +173,5 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24924) @@ -179,4 +179,5 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24924) @@ -5564,25 +5564,5 @@ return; }/*}}}*/ -void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ - - /*Computational flags:*/ - int bp_compute_fingerprints= 0; - - /*some paramters first: */ - this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); - - if(!IsOceanInElement()){ - /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested - *bottom pressure fingerprints:*/ - if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); - } - //if(!IsIceInElement()){ - /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); - //} - -} -/*}}}*/ -void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ +void Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/ /*diverse:*/ @@ -5595,70 +5575,26 @@ IssmDouble late,longe,re; IssmDouble lati,longi,ri; - bool notfullygrounded=false; /*elastic green function:*/ - int index; - IssmDouble* G_elastic_precomputed=NULL; - IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* indices=NULL; int M; - - /*ice properties: */ - IssmDouble rho_ice,rho_water,rho_earth; - - /*constants:*/ - IssmDouble constant=0; - - /*Initialize eustatic component: do not skip this step :):*/ - IssmDouble eustatic = 0.; + IssmDouble* dummypointer=NULL; /*Computational flags:*/ bool computerigid = true; - bool computeelastic= true; - bool scaleoceanarea= false; - - /*early return if we are not on an ice cap:*/ - if(!IsIceOnlyInElement()){ - constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); - *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! - return; - } - - /*early return if we are fully floating:*/ - Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input); - if (gr_input->GetInputMax()<=0){ - constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); - *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! - return; - } - - /*If we are here, we are on ice that is fully grounded or half-way to floating:*/ - if ((gr_input->GetInputMin())<0){ - notfullygrounded=true; //used later on. - } - - /*Inform mask: */ - constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); - - /*recover material parameters: */ - rho_ice=FindParam(MaterialsRhoIceEnum); - rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_earth=FindParam(MaterialsEarthDensityEnum); - - /*recover love numbers and computational flags: */ + + /*recover computational flags: */ this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); - this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); - this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); - - /*recover elastic green function:*/ - DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); - - if(computeelastic){ - parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); - } + if(!computerigid) return; //we are running eustatic solution only. /*how many dofs are we working with here? */ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*recover size of indexed tables:*/ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&dummypointer,&M); + + /*allocate indices:*/ + indices=xNew(gsize); /* Where is the centroid of this element?:{{{*/ @@ -5700,4 +5636,123 @@ longe=longe/180*PI; /*}}}*/ + + for(int i=0;i(alpha/PI*reCast(M-1)); + } + + /*Add in inputs:*/ + this->inputs2->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); + + return; +} +/*}}}*/ +void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + + /*Computational flags:*/ + int bp_compute_fingerprints= 0; + + /*some paramters first: */ + this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); + + if(!IsOceanInElement()){ + /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested + *bottom pressure fingerprints:*/ + if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + } + //if(!IsIceInElement()){ + /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ + this->SealevelriseEustaticIce(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea); + //} + +} +/*}}}*/ +void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/ + + /*diverse:*/ + int gsize,dummy; + bool spherical=true; + IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble area; + IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) + IssmDouble rho; + IssmDouble late,longe,re; + IssmDouble lati,longi,ri; + bool notfullygrounded=false; + + /*elastic green function:*/ + int index; + IssmDouble* indices=NULL; + IssmDouble* G_elastic_precomputed=NULL; + IssmDouble* G_rigid_precomputed=NULL; + int M; + + /*ice properties: */ + IssmDouble rho_ice,rho_water,rho_earth; + + /*constants:*/ + IssmDouble constant=0; + + /*Initialize eustatic component: do not skip this step :):*/ + IssmDouble eustatic = 0.; + + /*Computational flags:*/ + bool computerigid = true; + bool computeelastic= true; + bool scaleoceanarea= false; + + /*early return if we are not on an ice cap:*/ + if(!IsIceOnlyInElement()){ + constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! + return; + } + + /*early return if we are fully floating:*/ + Input2* gr_input=this->GetInput2(MaskOceanLevelsetEnum); _assert_(gr_input); + if (gr_input->GetInputMax()<=0){ + constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! + return; + } + + /*If we are here, we are on ice that is fully grounded or half-way to floating:*/ + if ((gr_input->GetInputMin())<0){ + notfullygrounded=true; //used later on. + } + + /*Inform mask: */ + constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + + /*recover material parameters: */ + rho_ice=FindParam(MaterialsRhoIceEnum); + rho_water=FindParam(MaterialsRhoSeawaterEnum); + rho_earth=FindParam(MaterialsEarthDensityEnum); + + /*recover love numbers and computational flags: */ + this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); + + /*recover elastic green function:*/ + DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M); + + if(computeelastic){ + parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); + } + + /*how many dofs are we working with here? */ + this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + + /*retrieve indices:*/ + if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} /*Compute area of element. Scale it by grounded fraction if not fully grounded: */ @@ -5755,8 +5810,9 @@ /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */ - lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; + /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; delPhi=fabs(lati-late); delLambda=fabs(longi-longe); alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - index=reCast(alpha/PI*reCast(M-1)); + index=reCast(alpha/PI*reCast(M-1));*/ + index=reCast(indices[i]); //Elastic component (from Eq 17 in Adhikari et al, GMD 2015) @@ -5767,4 +5823,7 @@ } } + + /*Free ressources:*/ + if(computerigid)xDelete(indices); /*Assign output pointer:*/ @@ -5924,5 +5983,5 @@ /*diverse:*/ - int gsize; + int gsize,dummy; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; @@ -5940,4 +5999,5 @@ int M; int index; + IssmDouble* indices=NULL; IssmDouble alpha; IssmDouble delPhi,delLambda; @@ -5974,4 +6034,7 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + /*retrieve indices:*/ + if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} + /*From Sg_old, recover water sea level rise:*/ S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; @@ -5979,40 +6042,4 @@ /*Compute area of element:*/ area=GetAreaSpherical(); - - /* Where is the centroid of this element?:{{{*/ - ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); - - minlong=400; 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*PI; - longe=longe/180*PI; - /*}}}*/ /*recover rigid and elastic green functions:*/ @@ -6030,8 +6057,10 @@ /*Compute alpha angle between centroid and current vertex : */ - lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; + /*lati=latitude[i]/180*PI; longi=longitude[i]/180*PI; delPhi=fabs(lati-late); delLambda=fabs(longi-longe); alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); - index=reCast(alpha/PI*(M-1)); + index=reCast(alpha/PI*(M-1));*/ + + index=reCast(indices[i]); /*Rigid earth gravitational perturbation: */ @@ -6046,4 +6075,7 @@ } + /*Free ressources:*/ + if(computerigid)xDelete(indices); + return; @@ -6053,5 +6085,5 @@ /*diverse:*/ - int gsize; + int gsize,dummy; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; @@ -6069,4 +6101,6 @@ IssmDouble* H_elastic_precomputed = NULL; int M; + IssmDouble* indices=NULL; + int index; /*computation of Green functions:*/ @@ -6109,43 +6143,9 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); + /*retrieve indices:*/ + if(computerigid){this->inputs2->GetArray(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} + /*compute area of element:*/ area=GetAreaSpherical(); - - /*element centroid (spherical): */ - /* Where is the centroid of this element?:{{{*/ - ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical); - - minlong=400; 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*PI; - longe=longe/180*PI; - /*}}}*/ /*figure out gravity center of our element (Cartesian): */ @@ -6192,11 +6192,10 @@ for(int i=0;i(indices[i]); /*Compute azimuths, both north and east components: */ @@ -6215,5 +6214,4 @@ /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ - int index=reCast(alpha/PI*(M-1)); U_elastic[i] += U_elastic_precomputed[index]; if(horiz){ @@ -6240,4 +6238,5 @@ /*free ressources:*/ + if(computerigid)xDelete(indices); xDelete(U_values); xDelete(U_elastic); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24924) @@ -166,4 +166,5 @@ IssmDouble OceanAverage(IssmDouble* Sg); void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24924) @@ -4684,4 +4684,14 @@ #endif #ifdef _HAVE_SEALEVELRISE_ +void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ + + /*Run sealevelrie geometry routine in elements:*/ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseGeometry(latitude,longitude,radius); + } + +} +/*}}}*/ void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24923) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24924) @@ -163,4 +163,5 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 24923) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 24924) @@ -55,4 +55,5 @@ void geodetic_core(FemModel* femmodel); void steric_core(FemModel* femmodel); +void sealevelrise_core_geometry(FemModel* femmodel); Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea); Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24924) @@ -41,4 +41,7 @@ /*set SLR configuration: */ femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum); + + /*run geometry core:*/ + sealevelrise_core_geometry(femmodel); /*Run geodetic:*/ @@ -304,4 +307,30 @@ } /*}}}*/ +void sealevelrise_core_geometry(FemModel* femmodel) { /*{{{*/ + + /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */ + + /*parameters: */ + bool spherical=true; + IssmDouble *latitude = NULL; + IssmDouble *longitude = NULL; + IssmDouble *radius = NULL; + + /*Verbose: */ + if(VerboseSolution()) _printf0_(" computing geometrical offsets into precomputed Green tables \n"); + + /*first, recover lat,long and radius vectors from vertices: */ + VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical); + + /*call the FemModel geometry module: */ + femmodel->SealevelriseGeometry(latitude, longitude, radius); + + /*Free ressources:*/ + xDelete(latitude); + xDelete(longitude); + xDelete(radius); + + +}/*}}}*/ Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 24924) @@ -99,5 +99,5 @@ } - //if(isslr) sealevelrise_geometrycore(femmodel); + if(isslr) sealevelrise_core_geometry(femmodel); while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24923) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 24924) @@ -683,4 +683,5 @@ syn keyword cConstant SealevelriseSpcthicknessEnum syn keyword cConstant SealevelriseHydroRateEnum +syn keyword cConstant SealevelriseIndicesEnum syn keyword cConstant SedimentHeadEnum syn keyword cConstant SedimentHeadOldEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24923) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 24924) @@ -680,4 +680,5 @@ SealevelriseSpcthicknessEnum, SealevelriseHydroRateEnum, + SealevelriseIndicesEnum, SedimentHeadEnum, SedimentHeadOldEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 24924) @@ -685,4 +685,5 @@ case SealevelriseSpcthicknessEnum : return "SealevelriseSpcthickness"; case SealevelriseHydroRateEnum : return "SealevelriseHydroRate"; + case SealevelriseIndicesEnum : return "SealevelriseIndices"; case SedimentHeadEnum : return "SedimentHead"; case SedimentHeadOldEnum : return "SedimentHeadOld"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24923) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 24924) @@ -700,4 +700,5 @@ else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum; else if (strcmp(name,"SealevelriseHydroRate")==0) return SealevelriseHydroRateEnum; + else if (strcmp(name,"SealevelriseIndices")==0) return SealevelriseIndicesEnum; else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum; else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum; else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; - else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; + if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum; + else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum; else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum; else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum; else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum; - else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; else stage=8; } if(stage==8){ - if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; + if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum; + else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum; else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum; else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum; @@ -997,9 +998,9 @@ else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; - else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; else stage=9; } if(stage==9){ - if (strcmp(name,"BoolInput")==0) return BoolInputEnum; + if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; + else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum; else if (strcmp(name,"IntInput2")==0) return IntInput2Enum; @@ -1120,9 +1121,9 @@ else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum; else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum; - else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; + if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; + else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum; else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum; else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum; @@ -1243,9 +1244,9 @@ else if (strcmp(name,"PentaInput")==0) return PentaInputEnum; else if (strcmp(name,"Profiler")==0) return ProfilerEnum; - else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; + if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum; + else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum; else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum; else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;