Index: /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 24009) +++ /issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp (revision 24010) @@ -17,5 +17,7 @@ /*First fetch data: */ + if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); ::CreateNodes(nodes,iomodel,GLheightadvectionAnalysisEnum,P1Enum); + iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface"); }/*}}}*/ int GLheightadvectionAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ @@ -56,7 +58,10 @@ ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*/ + /* Spawn basal element */ + if(!element->IsOnBase()) return NULL; + Element* basalelement = element->SpawnBasalElement(); + /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - //if(!element->IsFloating()) return NULL; + if(!basalelement->IsIceInElement()) return NULL; /*Intermediaries */ @@ -65,8 +70,11 @@ IssmDouble Jdet,D_scalar,onboundary; IssmDouble vel,vx,vy; - IssmDouble* xyz_list = NULL; + IssmDouble* xyz_list = NULL; + Input* vx_input = NULL; + Input* vy_input = NULL; + Input* bc_input = NULL; /*Get problem dimension*/ - element->FindParam(&domaintype,DomainTypeEnum); + basalelement->FindParam(&domaintype,DomainTypeEnum); switch(domaintype){ case Domain2DhorizontalEnum: dim = 2; break; @@ -76,8 +84,8 @@ /*Fetch number of nodes and dof for this finite element*/ - int numnodes = element->GetNumberOfNodes(); + int numnodes = basalelement->GetNumberOfNodes(); /*Initialize Element vector and other vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(); + ElementMatrix* Ke = basalelement->NewElementMatrix(); IssmDouble* basis = xNew(numnodes); IssmDouble* dbasis = xNew(dim*numnodes); @@ -86,22 +94,29 @@ /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinates(&xyz_list); - //Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input); - //Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input); - Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); - Input* bc_input=element->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input); - - IssmDouble h = element->CharacteristicLength(); + basalelement->GetVerticesCoordinates(&xyz_list); + switch(domaintype){ + case Domain2DhorizontalEnum: + vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input); + vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input); + bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input); + break; + case Domain3DEnum: + vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input); + vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input); + bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input); + break; + default: _error_("mesh "<CharacteristicLength(); /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGauss(4); + Gauss* gauss=basalelement->NewGauss(4); for(int ig=gauss->begin();igend();ig++){ gauss->GaussPoint(ig); - element->JacobianDeterminant(&Jdet,xyz_list,gauss); - element->NodalFunctions(basis,gauss); - element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); - GetBprime(Bprime,element,dim,xyz_list,gauss); + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + GetBprime(Bprime,basalelement,dim,xyz_list,gauss); /*Get velocity*/ @@ -168,9 +183,9 @@ xDelete(Bprime); delete gauss; + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; return Ke; }/*}}}*/ ElementVector* GLheightadvectionAnalysis::CreatePVector(Element* element){/*{{{*/ - return NULL; }/*}}}*/ @@ -203,5 +218,4 @@ for(int i=0;ielements->Size();i++){ Element *element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - int numnodes = element->GetNumberOfNodes(); IssmDouble *mask = xNew(numnodes); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24010) @@ -2157,4 +2157,61 @@ } /*}}}*/ +void Element::Ismip6FloatingiceMeltingRate(){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int numvertices = this->GetNumberOfVertices(); + int basinid,num_basins,M,N; + IssmDouble tf,gamma0,base,delta_t_basin,mean_tf_basin,absval; + IssmDouble basalmeltrate[numvertices]; + bool islocal; + IssmDouble* delta_t = NULL; + IssmDouble* mean_tf = NULL; + IssmDouble* depths = NULL; + + /*Get variables*/ + IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); + IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum); + IssmDouble lf = this->FindParam(MaterialsLatentheatEnum); + IssmDouble cp = this->FindParam(MaterialsMixedLayerCapacityEnum); + + /*Hard code sea water density to be consistent with ISMIP6 documentation*/ + rhow = 1028.; + + /* Get parameters and inputs */ + this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); + this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum); + this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum); + this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum); _assert_(M==num_basins); + this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum); + if(!islocal) { + this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins); + } + Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); + delta_t_basin = delta_t[basinid]; + if(!islocal) mean_tf_basin = mean_tf[basinid]; + + /*Compute melt rate for Local and Nonlocal parameterizations*/ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + tf_input->GetInputValue(&tf,gauss); + if(!islocal) { + absval = mean_tf_basin+delta_t_basin; + if (absval<0) {absval = -1.*absval;} + basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;} + else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);} + } + + /*Return basal melt rate*/ + this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum); + + /*Cleanup and return*/ + delete gauss; + xDelete(delta_t); + xDelete(mean_tf); + xDelete(depths); + +}/*}}}*/ bool Element::IsWaterInElement(){/*{{{*/ return (this->inputs->Max(MaskOceanLevelsetEnum)>0.); @@ -2520,4 +2577,280 @@ } /*}}}*/ +void Element::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int basin_id; + IssmDouble dist_gl,dist_cf; + + inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum); + IssmDouble boxid_max=reCast(max_boxid_basin_list[basin_id])+1.; + + Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input); + Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); + + /*Get dist_gl and dist_cf at center of element*/ + Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); + dist_gl_input->GetInputValue(&dist_gl,gauss); + dist_cf_input->GetInputValue(&dist_cf,gauss); + delete gauss; + + /*Ensure values are positive for floating ice*/ + dist_gl = fabs(dist_gl); + dist_cf = fabs(dist_cf); + + /*Compute relative distance to grounding line*/ + IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf); + + /*Assign box numbers based on rel_dist_gl*/ + int boxid = -1; + for(IssmDouble i=0.;i=lowbound && rel_dist_gl<=highbound){ + boxid=reCast(i); + break; + } + } + if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!"); + + this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); + +}/*}}}*/ +void Element::PicoUpdateBox(int loopboxid){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int boxid; + this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); + if(loopboxid!=boxid) return; + + int NUMVERTICES = this->GetNumberOfVertices(); + int basinid, maxbox, num_basins, numnodes, M; + IssmDouble gamma_T, overturning_coeff, thickness; + IssmDouble pressure, T_star,p_coeff, q_coeff; + bool isplume; + + /*Get variables*/ + IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); + IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum); + IssmDouble earth_grav = this->FindParam(ConstantsGEnum); + IssmDouble rho_star = 1033.; // kg/m^3 + IssmDouble nu = rhoi/rhow; + IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum); + IssmDouble c_p_ocean = this->FindParam(MaterialsMixedLayerCapacityEnum); + IssmDouble lambda = latentheat/c_p_ocean; + IssmDouble a = -0.0572; // K/psu + IssmDouble b = 0.0788 + this->FindParam(MaterialsMeltingpointEnum); //K + IssmDouble c = 7.77e-4; + IssmDouble alpha = 7.5e-5; // 1/K + IssmDouble Beta = 7.7e-4; // K + + /* Get non-box-specific parameters and inputs */ + this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum); + this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); + this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); + this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); + this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum); + Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); + _assert_(basinid<=num_basins); + + IssmDouble* boxareas = NULL; + this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); + _assert_(M==num_basins*maxbox); + + IssmDouble area_boxi = boxareas[basinid*maxbox+boxid]; + IssmDouble g1 = area_boxi*gamma_T; + + IssmDouble basalmeltrates_shelf[NUMVERTICES]; + IssmDouble potential_pressure_melting_point[NUMVERTICES]; + IssmDouble Tocs[NUMVERTICES]; + IssmDouble Socs[NUMVERTICES]; + + /* First box calculations */ + if(boxid==0){ + /* Get box1 parameters and inputs */ + IssmDouble time, toc_farocean, soc_farocean; + this->parameters->FindParam(&time,TimeEnum); + this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); + this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum); + IssmDouble s1 = soc_farocean/(nu*lambda); + IssmDouble overturnings[NUMVERTICES]; + + /* Start looping on the number of verticies and calculate ocean vars */ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + thickness_input->GetInputValue(&thickness,gauss); + pressure = (rhoi*earth_grav*1e-4)*thickness; + T_star = a*soc_farocean+b-c*pressure-toc_farocean; + p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha)); + q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha))); + + /* To avoid negatives under the square root */ + if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff; + + Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff)); + Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]); + potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; + if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); + overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i])); + } + + if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum); + + /*Cleanup and return*/ + delete gauss; + } + + /* Subsequent box calculations */ + else { + /* Get subsequent box parameters and inputs */ + IssmDouble* toc_weighted_avg = NULL; + IssmDouble* soc_weighted_avg = NULL; + IssmDouble* overturning_weighted_avg = NULL; + this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum); + this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum); + this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum); + IssmDouble mean_toc = toc_weighted_avg[basinid]; + IssmDouble mean_soc = soc_weighted_avg[basinid]; + IssmDouble mean_overturning = overturning_weighted_avg[basinid]; + IssmDouble g2 = g1/(nu*lambda); + + /* Start looping on the number of verticies and calculate ocean vars */ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + thickness_input->GetInputValue(&thickness,gauss); + pressure = (rhoi*earth_grav*1.e-4)*thickness; + T_star = a*mean_soc+b-c*pressure-mean_toc; + Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc)); + Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); + potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; + if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); + } + + if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); + this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); + + /*Cleanup and return*/ + xDelete(toc_weighted_avg); + xDelete(soc_weighted_avg); + xDelete(overturning_weighted_avg); + delete gauss; + } + + /*Cleanup and return*/ + xDelete(boxareas); + +}/*}}}*/ +void Element::PicoComputeBasalMelt(){/*{{{*/ + + if(!this->IsIceInElement() || !this->IsFloating()) return; + + int NUMVERTICES = this->GetNumberOfVertices(); + IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0; + IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat; + IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey; + + /*Get variables*/ + E0 = 3.6e-2; //Entrainment coefficient + Cd = 2.5e-3; //Drag coefficient + CdT = 1.1e-3; //Turbulent heat exchange coefficient + YT = CdT/sqrt(Cd); //Heat exchange coefficient + lam1 = -5.73e-2; //Freezing point-salinity coefficient (degrees C) + lam2 = 8.32e-2; //Freezing point offset (degrees C) + lam3 = 7.61e-4; //Freezing point-depth coefficient (K m-1) + M0 = 10.; //Melt-rate parameter (m yr-1 C-2) + CdTS0 = 6e-4; //Heat exchange parameter + y1 = 0.545; //Heat exchange parameter + y2 = 3.5e-5; //Heat exchange parameter + x0 = 0.56; //Dimentionless scaling factor + + /*Define arrays*/ + IssmDouble basalmeltrates_shelf[NUMVERTICES]; //Basal melt-rate + + /*Polynomial coefficients*/ + IssmDouble p[12]; + p[0] = 0.1371330075095435; + p[1] = 5.527656234709359E1; + p[2] = -8.951812433987858E2; + p[3] = 8.927093637594877E3; + p[4] = -5.563863123811898E4; + p[5] = 2.218596970948727E5; + p[6] = -5.820015295669482E5; + p[7] = 1.015475347943186E6; + p[8] = -1.166290429178556E6; + p[9] = 8.466870335320488E5; + p[10] = -3.520598035764990E5; + p[11] = 6.387953795485420E4; + + /*Get inputs*/ + Input* zgl_input = this->GetInput(GroundinglineHeightEnum); _assert_(zgl_input); + Input* toc_input = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(toc_input); + Input* soc_input = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(soc_input); + Input* base_input = this->GetInput(BaseEnum); _assert_(base_input); + Input* baseslopex_input = this->GetInput(BaseSlopeXEnum); _assert_(baseslopex_input); + Input* baseslopey_input = this->GetInput(BaseSlopeYEnum); _assert_(baseslopey_input); + this->FindParam(&yts, ConstantsYtsEnum); + + /*Loop over the number of vertices in this element*/ + Gauss* gauss=this->NewGauss(); + for(int i=0;iGaussVertex(i); + + /*Get inputs*/ + zgl_input->GetInputValue(&zgl,gauss); + toc_input->GetInputValue(&Toc,gauss); //(K) + soc_input->GetInputValue(&Soc,gauss); + base_input->GetInputValue(&z_base,gauss); + baseslopex_input->GetInputValue(&slopex,gauss); + baseslopey_input->GetInputValue(&slopey,gauss); + + /*Compute ice shelf base slope (radians)*/ + alpha = atan(sqrt(slopex*slopex + slopey*slopey)); + if(alpha>=M_PI) alpha = M_PI - 0.001; //ensure sin(alpha) > 0 for meltrate calculations + + /*Make necessary conversions*/ + Toc = Toc-273.15; + if(zgl>z_base) zgl=z_base; + + /*Low bound for Toc to ensure X_hat is between 0 and 1*/ + if(TocAddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); + + /*Cleanup and return*/ + delete gauss; + +}/*}}}*/ void Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 24010) @@ -134,4 +134,5 @@ bool IsIceInElement(); bool IsLandInElement(); + void Ismip6FloatingiceMeltingRate(); bool IsWaterInElement(); void LinearFloatingiceMeltingRate(); @@ -145,4 +146,7 @@ ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); + void PicoUpdateBoxid(int* pmax_boxid_basin); + void PicoUpdateBox(int loopboxid); + void PicoComputeBasalMelt(); void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac); void PositiveDegreeDaySicopolis(bool isfirnwarming); @@ -257,5 +261,4 @@ virtual bool IsFaceOnBoundary(void)=0; virtual bool IsIcefront(void)=0; - virtual void Ismip6FloatingiceMeltingRate(void)=0; virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; virtual bool IsOnBase()=0; @@ -299,7 +302,4 @@ virtual int NumberofNodesPressure(void)=0; virtual int NumberofNodesVelocity(void)=0; - virtual void PicoUpdateBoxid(int* pmax_boxid_basin)=0; - virtual void PicoUpdateBox(int loopboxid)=0; - virtual void PicoComputeBasalMelt(void)=0; virtual void PotentialUngrounding(Vector* potential_sheet_ungrounding)=0; virtual int PressureInterpolation()=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24010) @@ -95,5 +95,4 @@ IssmDouble IceVolume(bool scaled); IssmDouble IceVolumeAboveFloatation(bool scaled); - void Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");}; void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputExtrude(int enum_type,int start); @@ -144,7 +143,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; - void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; - void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24010) @@ -88,5 +88,4 @@ bool IsFaceOnBoundary(void){_error_("not implemented yet");}; bool IsIcefront(void); - void Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");}; bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; bool IsOnBase(){_error_("not implemented yet");}; @@ -130,7 +129,4 @@ int NumberofNodesPressure(void){_error_("not implemented yet");}; int NumberofNodesVelocity(void){_error_("not implemented yet");}; - void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; - void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; - void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 24010) @@ -88,5 +88,4 @@ IssmDouble IceVolume(bool scaled){_error_("not implemented yet");}; IssmDouble IceVolumeAboveFloatation(bool scaled){_error_("not implemented yet");}; - void Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");}; bool IsFaceOnBoundary(void){_error_("not implemented yet");}; bool IsIcefront(void); @@ -138,7 +137,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");}; - void PicoUpdateBox(int loopboxid){_error_("not implemented yet");}; - void PicoComputeBasalMelt(void){_error_("not implemented yet");}; void PotentialUngrounding(Vector* potential_sheet_ungrounding){_error_("not implemented yet");}; int PressureInterpolation(void); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 24010) @@ -2726,60 +2726,4 @@ return isicefront; }/*}}}*/ -void Tria::Ismip6FloatingiceMeltingRate(void){/*{{{*/ - - if(!this->IsIceInElement() || !this->IsFloating()) return; - - int basinid,num_basins,M,N; - IssmDouble tf,gamma0,base,delta_t_basin,mean_tf_basin,absval; - IssmDouble basalmeltrate[NUMVERTICES]; - bool islocal; - IssmDouble* delta_t = NULL; - IssmDouble* mean_tf = NULL; - IssmDouble* depths = NULL; - - /*Get variables*/ - IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); - IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum); - IssmDouble lf = this->FindParam(MaterialsLatentheatEnum); - IssmDouble cp = this->FindParam(MaterialsMixedLayerCapacityEnum); - - /*Hard code sea water density to be consistent with ISMIP6 documentation*/ - rhow = 1028.; - - /* Get parameters and inputs */ - this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); - this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum); - this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum); - this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum); _assert_(M==num_basins); - this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum); - if(!islocal) { - this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins); - } - Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); - delta_t_basin = delta_t[basinid]; - if(!islocal) mean_tf_basin = mean_tf[basinid]; - - /* Compute melt rate */ - Gauss* gauss=this->NewGauss(); - for(int i=0;iGaussVertex(i); - tf_input->GetInputValue(&tf,gauss); - if(!islocal) { - absval = mean_tf_basin+delta_t_basin; - if (absval<0) {absval = -1.*absval;} - basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;} - else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);} - } - - /*Return basal melt rate*/ - this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum); - - /*Cleanup and return*/ - delete gauss; - xDelete(delta_t); - xDelete(mean_tf); - xDelete(depths); - -}/*}}}*/ bool Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/ @@ -3269,276 +3213,4 @@ } /*}}}*/ -void Tria::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/ - - if(!this->IsIceInElement() || !this->IsFloating()) return; - - int basin_id; - IssmDouble dist_gl,dist_cf; - - inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum); - IssmDouble boxid_max=reCast(max_boxid_basin_list[basin_id])+1.; - - Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input); - Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input); - - /*Get dist_gl and dist_cf at center of element*/ - Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); - dist_gl_input->GetInputValue(&dist_gl,gauss); - dist_cf_input->GetInputValue(&dist_cf,gauss); - delete gauss; - - /*Ensure values are positive for floating ice*/ - dist_gl = fabs(dist_gl); - dist_cf = fabs(dist_cf); - - /*Compute relative distance to grounding line*/ - IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf); - - /*Assign box numbers based on rel_dist_gl*/ - int boxid = -1; - for(IssmDouble i=0.;i=lowbound && rel_dist_gl<=highbound){ - boxid=reCast(i); - break; - } - } - if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!"); - - this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid)); -}/*}}}*/ -void Tria::PicoUpdateBox(int loopboxid){/*{{{*/ - - if(!this->IsIceInElement() || !this->IsFloating()) return; - - int boxid; - this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); - if(loopboxid!=boxid) return; - - int basinid, maxbox, num_basins, numnodes, M; - IssmDouble gamma_T, overturning_coeff, thickness; - IssmDouble pressure, T_star,p_coeff, q_coeff; - bool isplume; - - /*Get variables*/ - IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum); - IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum); - IssmDouble earth_grav = this->FindParam(ConstantsGEnum); - IssmDouble rho_star = 1033.; // kg/m^3 - IssmDouble nu = rhoi/rhow; - IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum); - IssmDouble c_p_ocean = this->FindParam(MaterialsMixedLayerCapacityEnum); - IssmDouble lambda = latentheat/c_p_ocean; - IssmDouble a = -0.0572; // K/psu - IssmDouble b = 0.0788 + this->FindParam(MaterialsMeltingpointEnum); //K - IssmDouble c = 7.77e-4; - IssmDouble alpha = 7.5e-5; // 1/K - IssmDouble Beta = 7.7e-4; // K - - /* Get non-box-specific parameters and inputs */ - this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum); - this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum); - this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum); - this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); - this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum); - Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input); - _assert_(basinid<=num_basins); - - IssmDouble* boxareas = NULL; - this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum); - _assert_(M==num_basins*maxbox); - - IssmDouble area_boxi = boxareas[basinid*maxbox+boxid]; - IssmDouble g1 = area_boxi*gamma_T; - - IssmDouble basalmeltrates_shelf[NUMVERTICES]; - IssmDouble potential_pressure_melting_point[NUMVERTICES]; - IssmDouble Tocs[NUMVERTICES]; - IssmDouble Socs[NUMVERTICES]; - - /* First box calculations */ - if(boxid==0){ - /* Get box1 parameters and inputs */ - IssmDouble time, toc_farocean, soc_farocean; - this->parameters->FindParam(&time,TimeEnum); - this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum); - this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum); - IssmDouble s1 = soc_farocean/(nu*lambda); - IssmDouble overturnings[NUMVERTICES]; - - /* Start looping on the number of verticies and calculate ocean vars */ - Gauss* gauss=this->NewGauss(); - for(int i=0;iGaussVertex(i); - thickness_input->GetInputValue(&thickness,gauss); - pressure = (rhoi*earth_grav*1e-4)*thickness; - T_star = a*soc_farocean+b-c*pressure-toc_farocean; - p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha)); - q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha))); - - /* To avoid negatives under the square root */ - if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff; - - Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff)); - Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]); - potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; - if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); - overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i])); - } - - if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum); - - /*Cleanup and return*/ - delete gauss; - } - - /* Subsequent box calculations */ - else { - /* Get subsequent box parameters and inputs */ - IssmDouble* toc_weighted_avg = NULL; - IssmDouble* soc_weighted_avg = NULL; - IssmDouble* overturning_weighted_avg = NULL; - this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum); - this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum); - this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum); - IssmDouble mean_toc = toc_weighted_avg[basinid]; - IssmDouble mean_soc = soc_weighted_avg[basinid]; - IssmDouble mean_overturning = overturning_weighted_avg[basinid]; - IssmDouble g2 = g1/(nu*lambda); - - /* Start looping on the number of verticies and calculate ocean vars */ - Gauss* gauss=this->NewGauss(); - for(int i=0;iGaussVertex(i); - thickness_input->GetInputValue(&thickness,gauss); - pressure = (rhoi*earth_grav*1.e-4)*thickness; - T_star = a*mean_soc+b-c*pressure-mean_toc; - Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc)); - Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda)); - potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure; - if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]); - } - - if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum); - this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum); - - /*Cleanup and return*/ - xDelete(toc_weighted_avg); - xDelete(soc_weighted_avg); - xDelete(overturning_weighted_avg); - delete gauss; - } - - /*Cleanup and return*/ - xDelete(boxareas); -} -/*}}}*/ -void Tria::PicoComputeBasalMelt(void){/*{{{*/ - - if(!this->IsIceInElement() || !this->IsFloating()) return; - - IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0; - IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat; - IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey; - - /*Get variables*/ - E0 = 3.6e-2; //Entrainment coefficient - Cd = 2.5e-3; //Drag coefficient - CdT = 1.1e-3; //Turbulent heat exchange coefficient - YT = CdT/sqrt(Cd); //Heat exchange coefficient - lam1 = -5.73e-2; //Freezing point-salinity coefficient (degrees C) - lam2 = 8.32e-2; //Freezing point offset (degrees C) - lam3 = 7.61e-4; //Freezing point-depth coefficient (K m-1) - M0 = 10.; //Melt-rate parameter (m yr-1 C-2) - CdTS0 = 6e-4; //Heat exchange parameter - y1 = 0.545; //Heat exchange parameter - y2 = 3.5e-5; //Heat exchange parameter - x0 = 0.56; //Dimentionless scaling factor - - /*Define arrays*/ - IssmDouble basalmeltrates_shelf[NUMVERTICES]; //Basal melt-rate - - /*Polynomial coefficients*/ - IssmDouble p[12]; - p[0] = 0.1371330075095435; - p[1] = 5.527656234709359E1; - p[2] = -8.951812433987858E2; - p[3] = 8.927093637594877E3; - p[4] = -5.563863123811898E4; - p[5] = 2.218596970948727E5; - p[6] = -5.820015295669482E5; - p[7] = 1.015475347943186E6; - p[8] = -1.166290429178556E6; - p[9] = 8.466870335320488E5; - p[10] = -3.520598035764990E5; - p[11] = 6.387953795485420E4; - - /*Get inputs*/ - Input* zgl_input = this->GetInput(GroundinglineHeightEnum); _assert_(zgl_input); - Input* toc_input = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(toc_input); - Input* soc_input = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(soc_input); - Input* base_input = this->GetInput(BaseEnum); _assert_(base_input); - Input* baseslopex_input = this->GetInput(BaseSlopeXEnum); _assert_(baseslopex_input); - Input* baseslopey_input = this->GetInput(BaseSlopeYEnum); _assert_(baseslopey_input); - this->FindParam(&yts, ConstantsYtsEnum); - - /*Loop over the number of vertices in this element*/ - Gauss* gauss=this->NewGauss(); - for(int i=0;iGaussVertex(i); - - /*Get inputs*/ - zgl_input->GetInputValue(&zgl,gauss); - toc_input->GetInputValue(&Toc,gauss); //(K) - soc_input->GetInputValue(&Soc,gauss); - base_input->GetInputValue(&z_base,gauss); - baseslopex_input->GetInputValue(&slopex,gauss); - baseslopey_input->GetInputValue(&slopey,gauss); - - /*Compute ice shelf base slope (radians)*/ - alpha = atan(sqrt(slopex*slopex + slopey*slopey)); - if(alpha>=M_PI) alpha = M_PI - 0.001; //ensure sin(alpha) > 0 for meltrate calculations - - /*Make necessary conversions*/ - Toc = Toc-273.15; - if(zgl>z_base) zgl=z_base; - - /*Low bound for Toc to ensure X_hat is between 0 and 1*/ - if(TocAddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum); - - /*Cleanup and return*/ - delete gauss; -}/*}}}*/ void Tria::PotentialUngrounding(Vector* potential_ungrounding){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24009) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 24010) @@ -98,5 +98,4 @@ IssmDouble IcefrontMassFluxLevelset(bool scaled); IssmDouble GroundinglineMassFlux(bool scaled); - void Ismip6FloatingiceMeltingRate(void); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/}; @@ -116,7 +115,4 @@ int NumberofNodesPressure(void); int NumberofNodesVelocity(void); - void PicoUpdateBoxid(int* pmax_boxid_basin); - void PicoUpdateBox(int loopboxid); - void PicoComputeBasalMelt(void); void PotentialUngrounding(Vector* potential_sheet_ungrounding); int PressureInterpolation(); Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 24009) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (revision 24010) @@ -12,20 +12,20 @@ bool isplume; - /*First, reset all melt to 0 */ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - int numvertices = element->GetNumberOfVertices(); - IssmDouble* values = xNewZeroInit(numvertices); - element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); - xDelete(values); - } + /*First, reset all melt to 0 */ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + int numvertices = element->GetNumberOfVertices(); + IssmDouble* values = xNewZeroInit(numvertices); + element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum); + xDelete(values); + } /*PICO melt rate parameterization (Reese et al., 2018)*/ - femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); - UpdateBoxIdsPico(femmodel); - ComputeBoxAreasPico(femmodel); - for(int i=0;iparameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum); + UpdateBoxIdsPico(femmodel); + ComputeBoxAreasPico(femmodel); + for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; - element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); - element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - boxareas[basinid*maxbox+boxid]+=element->GetHorizontalSurfaceArea(); + if(!element->IsOnBase()) continue; + Element* basalelement = element->SpawnBasalElement(); + if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue; + basalelement->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum); + basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea(); + basalelement->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; } @@ -140,5 +144,5 @@ void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/ - int num_basins, basinid, maxbox, M; + int num_basins, basinid, maxbox, M, domaintype; IssmDouble area, toc, soc, overturning; IssmDouble* boxareas=NULL; @@ -157,20 +161,24 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; + if(!element->IsOnBase()) continue; + Element* basalelement = element->SpawnBasalElement(); + if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue; int el_boxid; - element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum); + basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum); if(el_boxid!=boxid) continue; - Input* tocs_input=element->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); - Input* socs_input=element->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input); - - element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); + Input* tocs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); + Input* socs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input); + + basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0); tocs_input->GetInputValue(&toc,gauss); socs_input->GetInputValue(&soc,gauss); delete gauss; - area=element->GetHorizontalSurfaceArea(); + area=basalelement->GetHorizontalSurfaceArea(); toc_weighted_avg[basinid]+=toc*area; soc_weighted_avg[basinid]+=soc*area; + basalelement->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; } @@ -194,17 +202,21 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; + if(!element->IsOnBase()) continue; + Element* basalelement = element->SpawnBasalElement(); + if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue; int el_boxid; - element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum); + basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum); if(el_boxid!=boxid) continue; - Input* overturnings_input=element->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input); - - element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); - Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); + Input* overturnings_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input); + + basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum); + Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0); overturnings_input->GetInputValue(&overturning,gauss); delete gauss; - area=element->GetHorizontalSurfaceArea(); + area=basalelement->GetHorizontalSurfaceArea(); overturning_weighted_avg[basinid]+=overturning*area; + basalelement->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; } Index: /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 24009) +++ /issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (revision 24010) @@ -73,5 +73,5 @@ void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/ - int num_basins, basinid,num_depths; + int num_basins, basinid, num_depths, domaintype; IssmDouble area, tf, base, time; bool islocal; @@ -154,13 +154,19 @@ for(int i=0;ielements->Size();i++){ Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - if(!element->IsIceInElement() || !element->IsFloating()) continue; - Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); - element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); - Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0); + if(!element->IsOnBase()) continue; + /*Spawn basal element if on base to compute element area*/ + Element* basalelement = element->SpawnBasalElement(); + if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue; + Input* tf_input=basalelement->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input); + basalelement->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum); + Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0); tf_input->GetInputValue(&tf,gauss); delete gauss; - area=element->GetHorizontalSurfaceArea(); + area=basalelement->GetHorizontalSurfaceArea(); tf_weighted_avg[basinid]+=tf*area; areas_summed[basinid] +=area; + /*Delete spawned element if we are in 3D*/ + basalelement->FindParam(&domaintype,DomainTypeEnum); + if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; } @@ -174,4 +180,10 @@ } + /*Compute meltrates*/ + for(int i=0;ielements->Size();i++){ + Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + element->Ismip6FloatingiceMeltingRate(); + } + /*Cleanup and return */ xDelete(tf_weighted_avg); @@ -180,10 +192,4 @@ xDelete(areas_summed_cpu); xDelete(tf_depths); - - /*Compute meltrates*/ - for(int i=0;ielements->Size();i++){ - Element* element=xDynamicCast(femmodel->elements->GetObjectByOffset(i)); - element->Ismip6FloatingiceMeltingRate(); - } } /*}}}*/ Index: /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m (revision 24009) +++ /issm/trunk-jpl/src/m/classes/basalforcingsismip6.m (revision 24010) @@ -19,6 +19,7 @@ function self = extrude(self,md) % {{{ self.basin_id=project3d(md,'vector',self.basin_id,'type','element','layer',1); - self.tf=project3d(md,'vector',self.tf,'type','element','layer',1); - self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1); + %self.tf=project3d(md,'vector',self.tf,'type','element','layer',1); + %self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1); + self.tf=project3d(md,'vector',self.tf,'type','node'); self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','element','layer',1); %bedrock only gets geothermal flux self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1); Index: /issm/trunk-jpl/src/m/classes/basalforcingspico.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcingspico.m (revision 24009) +++ /issm/trunk-jpl/src/m/classes/basalforcingspico.m (revision 24010) @@ -13,5 +13,5 @@ farocean_temperature = NaN; farocean_salinity = NaN; - isplume = NaN; + isplume = 0; geothermalflux = NaN; groundedice_melting_rate = NaN; @@ -68,5 +68,5 @@ md = checkfield(md,'fieldname','basalforcings.overturning_coeff','numel',1,'NaN',1,'Inf',1,'>',0); md = checkfield(md,'fieldname','basalforcings.gamma_T','numel',1,'NaN',1,'Inf',1,'>',0); - md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]); + md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins+1 NaN]); md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]); md = checkfield(md,'fieldname','basalforcings.isplume','values',[0 1]);