Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24967) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24968) @@ -5751,19 +5751,12 @@ void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ - /*Computational flags:*/ int bp_compute_fingerprints= 0; - - /*some paramters first: */ + + /*Compute bottom pressure contribution from ocean if requested:*/ this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); - - if(!masks->isoceanin[this->lid]){ - /*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,oceanarea); - } - //if(!IsIceInElement()){ - /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/ - this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea); - //} + if(bp_compute_fingerprints)this->SealevelriseBottomPressure(Sgi,masks); + + /*Compute eustatic ice contribution to sea level rise: */ + this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea); } @@ -5777,4 +5770,5 @@ IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) bool notfullygrounded=false; + bool scaleoceanarea= false; /*elastic green function:*/ @@ -5782,5 +5776,5 @@ /*ice properties: */ - IssmDouble rho_ice,rho_water,rho_earth; + IssmDouble rho_ice,rho_water; /*constants:*/ @@ -5789,9 +5783,4 @@ /*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:*/ @@ -5827,5 +5816,5 @@ rho_water=FindParam(MaterialsRhoSeawaterEnum); - /*recover love numbers and computational flags: */ + /*recover ocean area scaling: */ this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum); @@ -5836,18 +5825,18 @@ this->GetInput2Value(&area,AreaEnum); - /*factor to reduce area if we are not fully grounded:*/ - if(notfullygrounded){ /*{{{*/ + /*Compute fraction of the element that is grounded: */ + if(notfullygrounded){ IssmDouble xyz_list[NUMVERTICES][3]; ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias - } - /*}}}*/ - - /*Compute ice thickness: */ + } + else phi=1.0; + + /*Retrieve ice thickness at vertices: */ Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum); if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!"); - /*If we are fully grounded, take the average over the element: {{{*/ + /*/Average ice thickness over grounded area of the element only: {{{*/ if(!notfullygrounded)deltathickness_input->GetInputAverage(&I); else{ @@ -5876,5 +5865,5 @@ /*}}}*/ - /*Compute eustatic compoent:*/ + /*Compute eustatic component:*/ _assert_(oceanarea>0.); if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 @@ -5893,152 +5882,45 @@ } /*}}}*/ -void Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble oceanarea){ /*{{{*/ +void Tria::SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ /*diverse:*/ int gsize; - bool spherical=true; - IssmDouble llr_list[NUMVERTICES][3]; - IssmDouble area,planetarea; - IssmDouble I; //change in ice thickness or water level(Farrel and Clarke, Equ. 4) - IssmDouble rho; - IssmDouble late,longe,re; - IssmDouble lati,longi,ri; + IssmDouble area; + IssmDouble BP; //change in bottom pressure (Farrel and Clarke, Equ. 4) /*elastic green function:*/ - IssmDouble* G_elastic_precomputed=NULL; - int M; - - /*ice properties: */ - IssmDouble 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; - bool bp_compute_fingerprints= false; - - /*we are here to compute fingerprints originating fromn bottom pressure loads:*/ + IssmDouble* G=NULL; + + /*water properties: */ + IssmDouble rho_water; + + /*we are here to compute fingerprints originating fromn bottom pressure changes:*/ + if(!masks->isoceanin[this->lid])return; /*Inform mask: */ + #ifdef _ISSM_DEBUG_ constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum); + #endif /*recover material parameters: */ rho_water=FindParam(MaterialsRhoSeawaterEnum); - rho_earth=FindParam(MaterialsEarthDensityEnum); - - /*recover earth area: */ - this->parameters->FindParam(&planetarea,SealevelPlanetAreaEnum); - - /*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:*/ - if(computeelastic){ - DoubleVecParam* 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); - - /* Where is the centroid of this element?:{{{*/ - - /*retrieve coordinates: */ - ::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*PI; - longe=longe/180*PI; - /*}}}*/ + + /*retrieve precomputed G:*/ + this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ this->GetInput2Value(&area,AreaEnum); - /*Compute bottom pressure change: */ + /*Retrieve bottom pressure change and average over the element: */ Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor); if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!"); - - /*If we are fully grounded, take the average over the element: */ - bottompressure_change_input->GetInputAverage(&I); - - /*Compute eustatic compoent:*/ - _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - - /*We do not need to add the bottom pressure component to the eustatic value: */ - eustatic += 0; - IssmDouble* latitude=NULL; //NOT GOING TO WORK! - IssmDouble* longitude=NULL; //NOT GOING TO WORK! - - if(computeelastic | computerigid){ - IssmDouble alpha; - IssmDouble delPhi,delLambda; - for(int i=0;i(alpha/PI*reCast(M-1)); - G_elastic += G_elastic_precomputed[index]; - } - - /*Add all components to the pSgi or pSgo solution vectors:*/ - Sgi[i]+=3*rho_water/rho_earth*area/planetarea*I*(G_rigid+G_elastic); - } - } - - /*Assign output pointer:*/ - _assert_(!xIsNan(eustatic)); - _assert_(!xIsInf(eustatic)); - *peustatic=eustatic; + bottompressure_change_input->GetInputAverage(&BP); + + /*convert from m to kg/m^2:*/ + BP=BP*rho_water; + + /*convolve:*/ + for(int i=0;i