Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21030) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 21031) @@ -3699,16 +3699,16 @@ int* indices=xNew(gsize); IssmDouble* values=xNew(gsize); + IssmDouble alpha; + IssmDouble delPhi,delLambda; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; + S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; /*Compute area of element:*/ @@ -3840,32 +3840,29 @@ int* indices=xNew(gsize); IssmDouble* values=xNewZeroInit(gsize); + IssmDouble alpha; + IssmDouble delPhi,delLambda; for(int i=0;i(alpha/PI*(M-1)); - G_elastic[i] += G_elastic_precomputed[index]; - } - } - - /*Add all components to the pSgo solution vectors:*/ - if(computerigid)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_rigid[i]; - if(computeelastic)values[i]+=3*rho_water/rho_earth*area/eartharea*I*G_elastic[i]; + + /*Compute alpha angle between centroid and current vertex : */ + 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))); + + //Rigid earth gravitational perturbation: + if(computerigid){ + G_rigid[i]=1.0/2.0/sin(alpha/2.0); + values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i]; + } + + //Elastic component (from Eq 17 in Adhikari et al, GMD 2015) + if(computeelastic){ + int index=reCast(alpha/PI*(M-1)); + G_elastic[i] += G_elastic_precomputed[index]; + values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i]; + } } @@ -3891,5 +3888,5 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble area; - IssmDouble I; //change in relative sea level or ice thickness + IssmDouble I, S; //change in relative ice thickness and sea level IssmDouble late,longe,re; IssmDouble lati,longi,ri; @@ -3922,6 +3919,6 @@ this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum); - /*early return if rigid or elastic not requested:*/ - if(!computerigid && !computeelastic) return; + /*early return if elastic not requested:*/ + if(!computeelastic) return; /*recover material parameters: */ @@ -3980,17 +3977,22 @@ z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0; - if(computeelastic){ + /*recover elastic Green's functions for displacement:*/ + DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter); + DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter); + U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); + H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); + + /*From Sg, recover water sea level rise:*/ + S=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - /*recover elastic Green's functions for displacement:*/ - DoubleVecParam* U_parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter); - DoubleVecParam* H_parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter); - U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M); - H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M); - - /*initialize: */ - U_elastic=xNewZeroInit(gsize); - N_elastic=xNewZeroInit(gsize); - E_elastic=xNewZeroInit(gsize); - } + /*Compute ice thickness change: */ + Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); + if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); + deltathickness_input->GetInputAverage(&I); + + /*initialize: */ + U_elastic=xNewZeroInit(gsize); + N_elastic=xNewZeroInit(gsize); + E_elastic=xNewZeroInit(gsize); int* indices=xNew(gsize); @@ -3998,67 +4000,48 @@ IssmDouble* N_values=xNewZeroInit(gsize); IssmDouble* E_values=xNewZeroInit(gsize); + IssmDouble alpha; + IssmDouble delPhi,delLambda; + IssmDouble dx, dy, dz, x, y, z; + IssmDouble N_azim, E_azim; for(int i=0;i(alpha/PI*(M-1)); - U_elastic[i] += U_elastic_precomputed[index]; - N_elastic[i] += H_elastic_precomputed[index]*N_azim; - E_elastic[i] += H_elastic_precomputed[index]*E_azim; - } + + /*Compute alpha angle between centroid and current vertex: */ + 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))); + + /*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); + + /*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]; + 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(computerigid){ - U_values[i]+=0; N_values[i]+=0; E_values[i]+=0; - } - if(computeelastic){ - - if(this->inputs->Max(MaskIceLevelsetEnum)<0){ - - /*Compute ice thickness change: */ - Input* deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum); - if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); - deltathickness_input->GetInputAverage(&I); - - U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; - N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; - E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; - } - else if(IsWaterInElement()) { - - /*From Sg, recover water sea level rise:*/ - I=0; for(int i=0;ivertices[i]->Sid()]/NUMVERTICES; - - U_values[i]+=3*rho_water/rho_earth*area/eartharea*I*U_elastic[i]; - N_values[i]+=3*rho_water/rho_earth*area/eartharea*I*N_elastic[i]; - E_values[i]+=3*rho_water/rho_earth*area/eartharea*I*E_elastic[i]; - } - } + if(this->inputs->Max(MaskIceLevelsetEnum)<0){ + U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i]; + N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i]; + E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i]; + } + else if(IsWaterInElement()) { + U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i]; + N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i]; + E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i]; + } } pUp->SetValues(gsize,indices,U_values,ADD_VAL); @@ -4069,9 +4052,5 @@ xDelete(indices); xDelete(U_values); xDelete(N_values); xDelete(E_values); - - /*Free ressources:*/ - if(computeelastic) { - xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); - } + xDelete(U_elastic); xDelete(N_elastic); xDelete(E_elastic); return;