Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24946) @@ -380,7 +380,7 @@ virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0; - virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; + virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0; virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0; - virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz)=0; + virtual void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24946) @@ -215,8 +215,8 @@ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_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 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24946) @@ -174,8 +174,8 @@ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24946) @@ -180,8 +180,8 @@ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");}; - void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");}; + void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; #endif Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24946) @@ -5578,10 +5578,10 @@ } /*}}}*/ -void Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/ - +void Tria::SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){ /*{{{*/ /*diverse:*/ int gsize; bool spherical=true; IssmDouble llr_list[NUMVERTICES][3]; + IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble area,eartharea; IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) @@ -5590,11 +5590,18 @@ IssmDouble lati,longi,ri; IssmDouble constant; + IssmDouble x_element,y_element,z_element,x,y,z,dx,dy,dz,N_azim,E_azim; IssmDouble* G=NULL; + IssmDouble* GU=NULL; + IssmDouble* GN=NULL; + IssmDouble* GE=NULL; IssmDouble* G_elastic_precomputed=NULL; IssmDouble* G_rigid_precomputed=NULL; + IssmDouble* U_elastic_precomputed=NULL; + IssmDouble* H_elastic_precomputed=NULL; /*elastic green function:*/ IssmDouble* indices=NULL; + int index; int M; @@ -5602,4 +5609,5 @@ bool computerigid = true; bool computeelastic = true; + int horiz; /*recover parameters: */ @@ -5609,4 +5617,5 @@ this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); + this->parameters->FindParam(&horiz,SealevelriseHorizEnum); /*recover precomputed green function kernels:*/ @@ -5616,11 +5625,28 @@ parameter = static_cast(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M); - + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + /*allocate indices:*/ indices=xNew(gsize); G=xNewZeroInit(gsize); + GU=xNewZeroInit(gsize); + if(horiz){ + GN=xNewZeroInit(gsize); + GE=xNewZeroInit(gsize); + } /*compute area:*/ area=GetAreaSpherical(); + + /*figure out gravity center of our element (Cartesian): */ + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; + y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; + z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; /* Where is the centroid of this element?:{{{*/ @@ -5674,13 +5700,33 @@ alpha=2.*asin(sqrt(pow(sin(delPhi/2),2)+cos(lati)*cos(late)*pow(sin(delLambda/2),2))); indices[i]=alpha/PI*reCast(M-1); - + index=reCast(indices[i]); + /*Rigid earth gravitational perturbation: */ if(computerigid){ - G[i]+=G_rigid_precomputed[reCast(indices[i])]; + G[i]+=G_rigid_precomputed[index]; } if(computeelastic){ - G[i]+=G_elastic_precomputed[reCast(indices[i])]; + G[i]+=G_elastic_precomputed[index]; } G[i]=G[i]*constant; + + /*Elastic components:*/ + GU[i] = constant * U_elastic_precomputed[index]; + if(horiz){ + /*Compute azimuths, both north and east components: */ + x = xx[i]; y = yy[i]; z = zz[i]; + if(latitude[i]==90){ + x=1e-12; y=1e-12; + } + if(latitude[i]==-90){ + x=1e-12; y=1e-12; + } + dx = x_element-x; dy = y_element-y; dz = z_element-z; + N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); + + GN[i] = constant*H_elastic_precomputed[index]*N_azim; + GE[i] = constant*H_elastic_precomputed[index]*E_azim; + } } @@ -5688,4 +5734,9 @@ this->inputs2->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize); this->inputs2->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize); + this->inputs2->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize); + if(horiz){ + this->inputs2->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize); + this->inputs2->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); + } this->inputs2->SetDoubleInput(AreaEnum,this->lid,area); @@ -5693,4 +5744,9 @@ xDelete(indices); xDelete(G); + xDelete(GU); + if(horiz){ + xDelete(GN); + xDelete(GE); + } return; @@ -6034,41 +6090,18 @@ } /*}}}*/ -void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){ /*{{{*/ +void Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/ /*diverse:*/ - int gsize,dummy; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble area,eartharea; + int gsize; IssmDouble I, S; //change in relative ice thickness and sea level - IssmDouble late,longe,re; - IssmDouble lati,longi,ri; - IssmDouble rho_ice,rho_water,rho_earth; - IssmDouble minlong=400; - IssmDouble maxlong=-20; + IssmDouble rho_ice,rho_water; + int horiz; /*precomputed elastic green functions:*/ - IssmDouble* U_elastic_precomputed = NULL; - IssmDouble* H_elastic_precomputed = NULL; - int M; - IssmDouble* indices=NULL; - int index; - - /*computation of Green functions:*/ - IssmDouble* U_elastic= NULL; - IssmDouble* N_elastic= NULL; - IssmDouble* E_elastic= NULL; - DoubleVecParam* U_parameter = NULL; - DoubleVecParam* H_parameter = NULL; - IssmDouble* U_values=NULL; - IssmDouble* N_values=NULL; - IssmDouble* E_values=NULL; - - /*optimization:*/ - bool store_green_functions=false; + IssmDouble* GU=NULL; + IssmDouble* GN=NULL; + IssmDouble* GE=NULL; /*computational flags:*/ - bool computerigid = true; bool computeelastic= true; @@ -6079,7 +6112,7 @@ if(masks->isfullyfloating[this->lid])return; - /*recover computational flags: */ - this->parameters->FindParam(&computerigid,SealevelriseRigidEnum); + /*recover parameters:*/ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); + this->parameters->FindParam(&horiz,SealevelriseHorizEnum); /*early return if elastic not requested:*/ @@ -6089,31 +6122,10 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_earth=FindParam(MaterialsEarthDensityEnum); - - /*how many dofs are we working with here? */ - this->parameters->FindParam(&gsize,MeshNumberofverticesEnum); - - /*recover earth area: */ - this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum); - - /*retrieve indices:*/ - if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);} - - /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ - this->GetInput2Value(&area,AreaEnum); - - /*figure out gravity center of our element (Cartesian): */ - IssmDouble x_element, y_element, z_element; - ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0; - y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0; - z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; /*recover elastic Green's functions for displacement:*/ - U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter); - U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize); if(horiz){ - H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter); - H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize); + this->inputs2->GetArrayPtr(SealevelriseGUEnum,this->lid,&GU,&gsize); } @@ -6126,79 +6138,21 @@ deltathickness_input->GetInputAverage(&I); - /*initialize: */ - U_elastic=xNewZeroInit(gsize); - if(horiz){ - N_elastic=xNewZeroInit(gsize); - E_elastic=xNewZeroInit(gsize); - } - - //int* indices=xNew(gsize); - //U_values=xNewZeroInit(gsize); - if(horiz){ - //N_values=xNewZeroInit(gsize); - //E_values=xNewZeroInit(gsize); - } - IssmDouble alpha; - IssmDouble delPhi,delLambda; - IssmDouble dx, dy, dz, x, y, z; - IssmDouble N_azim, E_azim; - - /*we are going to use the result of these masks a lot: */ - bool isiceonlyinelement=masks->isiceonly[this->lid]; - bool isoceaninelement=masks->isoceanin[this->lid]; - - for(int i=0;i(indices[i]); - - /*Compute azimuths, both north and east components: */ - x = xx[i]; y = yy[i]; z = zz[i]; - if(latitude[i]==90){ - x=1e-12; y=1e-12; - } - if(latitude[i]==-90){ - x=1e-12; y=1e-12; - } - dx = x_element-x; dy = y_element-y; dz = z_element-z; - if(horiz){ - N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5); - } - - /*Elastic component (from Eq 17 in Adhikari et al, GMD 2015): */ - U_elastic[i] += U_elastic_precomputed[index]; - if(horiz){ - N_elastic[i] += H_elastic_precomputed[index]*N_azim; - E_elastic[i] += H_elastic_precomputed[index]*E_azim; - } - - /*Add all components to the pUp solution vectors:*/ - if(isiceonlyinelement){ - Up[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; + if (masks->isiceonly[this->lid]){ + for(int i=0;iisoceanin[this->lid]){ + for(int i=0;i(U_values); - xDelete(U_elastic); - if(horiz){ - xDelete(N_values); xDelete(E_values); - xDelete(N_elastic); xDelete(E_elastic); } Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24945) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24946) @@ -166,10 +166,10 @@ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); void SetSealevelMasks(SealevelMasks* masks); - void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); + void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea); void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius); - void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz); + void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz); #endif /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24945) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24946) @@ -4673,10 +4673,10 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius) { /*{{{*/ +void FemModel::SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz) { /*{{{*/ /*Run sealevelrie geometry routine in elements:*/ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeometry(latitude,longitude,radius); + element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz); } @@ -4854,5 +4854,5 @@ } /*}}}*/ -void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop){/*{{{*/ /*serialized vectors:*/ @@ -4863,6 +4863,9 @@ IssmDouble* East = NULL; int* indices = NULL; - int gsize; - + int gsize; + int horiz; + + /*retrieve parameters:*/ + this->parameters->FindParam(&horiz,SealevelriseHorizEnum); /*Serialize vectors from previous iteration:*/ @@ -4880,5 +4883,5 @@ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz); + element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz); } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24945) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24946) @@ -164,8 +164,8 @@ #ifdef _HAVE_SEALEVELRISE_ void SealevelriseEustatic(Vector* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); - void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop); IssmDouble SealevelriseOceanAverage(Vector* Sg,SealevelMasks* masks, IssmDouble oceanarea); #endif