Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 20155) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 20156) @@ -673,10 +673,11 @@ IssmDouble* newthickness = xNew(numnodes); IssmDouble* deltathickness = xNew(numnodes); - IssmDouble* newbed = xNew(numnodes); + IssmDouble* newbase = xNew(numnodes); IssmDouble* newsurface = xNew(numnodes); IssmDouble* oldthickness = xNew(numnodes); - IssmDouble* oldbed = xNew(numnodes); + IssmDouble* oldbase = xNew(numnodes); IssmDouble* oldsurface = xNew(numnodes); IssmDouble* phi = xNew(numnodes); + IssmDouble* sealevel = xNew(numnodes); /*Use the dof list to index into the solution vector: */ @@ -689,12 +690,12 @@ } - /*Get previous bed, thickness and surface*/ - basalelement->GetInputListOnNodes(&oldbed[0],BaseEnum); + /*Get previous base, thickness, surfac and current sealevel:*/ + basalelement->GetInputListOnNodes(&oldbase[0],BaseEnum); basalelement->GetInputListOnNodes(&oldsurface[0],SurfaceEnum); basalelement->GetInputListOnNodes(&oldthickness[0],ThicknessEnum); basalelement->GetInputListOnNodes(&phi[0],MaskGroundediceLevelsetEnum); - - - /*What is the delta thickness, useful for Sea-level rise:*/ + basalelement->GetInputListOnNodes(&sealevel[0],SealevelEnum); + + /*What is the delta thickness forcing the sea-level rise core: */ for(i=0;i0.){ - newsurface[i] = oldbed[i]+newthickness[i]; //surface = oldbed + newthickness - newbed[i] = oldbed[i]; //same bed: do nothing - } - else{ //this is an ice shelf + if (phi[i]>0.){ //this is an ice sheet: just add thickness to base. + newsurface[i] = oldbase[i]+newthickness[i]; //surface = oldbase + newthickness + newbase[i] = oldbase[i]; //same base: do nothing + } + else{ //this is an ice shelf: hydrostatic equilibrium*/ if(hydroadjustment==AbsoluteEnum){ - newsurface[i] = newthickness[i]*(1-rho_ice/rho_water); - newbed[i] = newthickness[i]*(-rho_ice/rho_water); + newsurface[i] = newthickness[i]*(1-rho_ice/rho_water)+sealevel[i]; + newbase[i] = newthickness[i]*(-rho_ice/rho_water)+sealevel[i]; } else if(hydroadjustment==IncrementalEnum){ - newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH - newbed[i] = oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH + newsurface[i] = oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i])+sealevel[i]; //surface = oldsurface + (1-di) * dH + newbase[i] = oldbase[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i])+sealevel[i]; //base = oldbed + di * dH } else _error_("Hydrostatic adjustment " << hydroadjustment << " (" << EnumToStringx(hydroadjustment) << ") not supported yet"); @@ -727,13 +727,13 @@ element->AddBasalInput(SealevelriseDeltathicknessEnum,deltathickness,P1Enum); element->AddBasalInput(SurfaceEnum,newsurface,P1Enum); - element->AddBasalInput(BaseEnum,newbed,P1Enum); + element->AddBasalInput(BaseEnum,newbase,P1Enum); /*Free ressources:*/ xDelete(newthickness); - xDelete(newbed); + xDelete(newbase); xDelete(newsurface); xDelete(oldthickness); xDelete(deltathickness); - xDelete(oldbed); + xDelete(oldbase); xDelete(oldsurface); xDelete(phi); Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 20155) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 20156) @@ -1610,5 +1610,5 @@ /*Intermediaries*/ int dim,domaintype; - IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure; + IssmDouble Jdet,thickness,base,sealevel,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble *xyz_list = NULL; @@ -1635,4 +1635,5 @@ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); @@ -1648,10 +1649,11 @@ gauss->GaussPoint(ig); thickness_input->GetInputValue(&thickness,gauss); - base_input->GetInputValue(&bed,gauss); + sealevel_input->GetInputValue(&sealevel,gauss); + base_input->GetInputValue(&base,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); - surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level - base_under_water = min(0.,bed); // 0 if the bottom of the glacier is above water level + surface_under_water = min(0.,thickness+base-sealevel); // 0 if the top of the glacier is above water level + base_under_water = min(0.,base-sealevel); // 0 if the bottom of the glacier is above water level water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water); ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; @@ -2070,5 +2072,5 @@ /*Intermediaries*/ - IssmDouble Jdet,thickness,bed,water_pressure,ice_pressure; + IssmDouble Jdet,thickness,bed,sealevel,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble *xyz_list = NULL; @@ -2086,4 +2088,5 @@ Input* thickness_input = element->GetInput(ThicknessEnum); _assert_(thickness_input); Input* base_input = element->GetInput(BaseEnum); _assert_(base_input); + Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); @@ -2100,9 +2103,10 @@ thickness_input->GetInputValue(&thickness,gauss); base_input->GetInputValue(&bed,gauss); + sealevel_input->GetInputValue(&sealevel,gauss); element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); element->NodalFunctions(basis,gauss); - surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level - base_under_water = min(0.,bed); // 0 if the bottom of the glacier is above water level + surface_under_water = min(0.,thickness+bed-sealevel); // 0 if the top of the glacier is above water level + base_under_water = min(0.,bed-sealevel); // 0 if the bottom of the glacier is above water level water_pressure = 1.0/2.0*gravity*rho_water*(surface_under_water*surface_under_water - base_under_water*base_under_water); ice_pressure = 1.0/2.0*gravity*rho_ice*thickness*thickness; @@ -2584,5 +2588,5 @@ /*Intermediaries*/ int dim; - IssmDouble Jdet,surface,z,water_pressure,ice_pressure; + IssmDouble Jdet,surface,sealevel,z,water_pressure,ice_pressure; IssmDouble surface_under_water,base_under_water,pressure; IssmDouble* xyz_list = NULL; @@ -2604,4 +2608,5 @@ /*Retrieve all inputs and parameters*/ Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); @@ -2622,4 +2627,5 @@ gauss->GaussPoint(ig); surface_input->GetInputValue(&surface,gauss); + sealevel_input->GetInputValue(&sealevel,gauss); if(dim==3) z=element->GetZcoord(xyz_list,gauss); else z=element->GetYcoord(xyz_list,gauss); @@ -2627,5 +2633,5 @@ element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss); - water_pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level + water_pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level ice_pressure = rho_ice*gravity*(surface-z); pressure = ice_pressure + water_pressure; @@ -3811,5 +3817,5 @@ /*Intermediaries*/ int i,dim; - IssmDouble Jdet,pressure,surface,z; + IssmDouble Jdet,pressure,surface,sealevel,z; IssmDouble normal[3]; IssmDouble *xyz_list = NULL; @@ -3843,4 +3849,5 @@ element->NormalSection(&normal[0],xyz_list_front); Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input); + Input* sealevel_input = element->GetInput(SealevelEnum); _assert_(sealevel_input); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); @@ -3859,7 +3866,8 @@ element->NodalFunctionsVelocity(vbasis,gauss); surface_input->GetInputValue(&surface,gauss); + sealevel_input->GetInputValue(&sealevel,gauss); if(dim==3) z=element->GetZcoord(xyz_list,gauss); else z=element->GetYcoord(xyz_list,gauss); - pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level + pressure = rho_water*gravity*min(0.,z-sealevel);//0 if the gaussian point is above water level for (int i=0;i(numvertices); IssmDouble* r = xNew(numvertices); + IssmDouble* sl = xNew(numvertices); /*Recover info at the vertices: */ @@ -1661,4 +1662,5 @@ GetInputListOnVertices(&b[0],BaseEnum); GetInputListOnVertices(&r[0],BedEnum); + GetInputListOnVertices(&sl[0],SealevelEnum); GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum); rho_water = matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -1684,14 +1686,14 @@ /*Change only if AggressiveMigration or if the ice sheet is in contact with the ocean*/ else{ // phi>0 - bed_hydro=-density*h[i]; + bed_hydro=-density*h[i]+sl[i]; if (bed_hydro>r[i]){ /*Unground only if the element is connected to the ice shelf*/ if(migration_style==AggressiveMigrationEnum || migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){ - s[i] = (1-density)*h[i]; - b[i] = -density*h[i]; + s[i] = (1-density)*h[i]+sl[i]; + b[i] = -density*h[i]+sl[i]; } else if(migration_style==SoftMigrationEnum && phi_ungrounding[vertices[i]->Pid()]<0.){ - s[i] = (1-density)*h[i]; - b[i] = -density*h[i]; + s[i] = (1-density)*h[i]+sl[i]; + b[i] = -density*h[i]+sl[i]; } else{ @@ -1705,10 +1707,10 @@ for(i=0;iPid()]<0.){ - phi[i]=h[i]+r[i]/density; + phi[i]=h[i]+(r[i]-sl[i])/density; } } - else if(migration_style!=ContactEnum) phi[i]=h[i]+r[i]/density; + else if(migration_style!=ContactEnum) phi[i]=h[i]+(r[i]-sl[i])/density; else{ /*do nothing*/ @@ -1727,4 +1729,5 @@ xDelete(b); xDelete(s); + xDelete(sl); xDelete(h); Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 20155) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 20156) @@ -67,5 +67,5 @@ IssmDouble Neff; IssmDouble drag_coefficient; - IssmDouble bed,thickness; + IssmDouble bed,thickness,sealevel; IssmDouble alpha_complement; @@ -75,4 +75,5 @@ element->GetInputValue(&thickness, gauss,ThicknessEnum); element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -85,5 +86,5 @@ //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); + Neff=gravity*(rho_ice*thickness+rho_water*(bed-sealevel)); if(Neff<0)Neff=0; @@ -231,5 +232,5 @@ IssmDouble drag_p, drag_q; IssmDouble Neff; - IssmDouble thickness,base,bed,floatation_thickness; + IssmDouble thickness,base,bed,floatation_thickness,sealevel; IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient,drag_coefficient_coulomb; @@ -241,4 +242,5 @@ element->GetInputValue(&thickness, gauss,ThicknessEnum); element->GetInputValue(&base, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&bed, gauss,BedEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); @@ -253,5 +255,5 @@ //From base and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*base); + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); if(Neff<0)Neff=0; @@ -407,5 +409,5 @@ IssmDouble drag_p, drag_q; IssmDouble Neff; - IssmDouble thickness,bed; + IssmDouble thickness,base,sealevel; IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient; @@ -416,5 +418,6 @@ element->GetInputValue(&drag_q,FrictionQEnum); element->GetInputValue(&thickness, gauss,ThicknessEnum); - element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&base, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); @@ -426,6 +429,6 @@ s=1./drag_p; - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); + //From base and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel)); if(Neff<0)Neff=0; @@ -467,5 +470,5 @@ IssmDouble drag_p, drag_q; IssmDouble Neff,F; - IssmDouble thickness,bed; + IssmDouble thickness,bed,sealevel; IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient,water_layer; @@ -478,4 +481,5 @@ element->GetInputValue(&thickness, gauss,ThicknessEnum); element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); element->GetInputValue(&water_layer, gauss,FrictionWaterLayerEnum); @@ -490,5 +494,5 @@ //From bed and thickness, compute effective pressure when drag is viscous: if(bed>0) bed=0; - if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*bed; + if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(bed-sealevel); else if(water_layer>0) Neff=gravity*rho_ice*thickness*F; else _error_("negative water layer thickness"); @@ -599,5 +603,5 @@ IssmDouble Neff; IssmDouble drag_coefficient; - IssmDouble bed,thickness,head; + IssmDouble bed,thickness,head,sealevel; IssmDouble alpha2; @@ -606,4 +610,5 @@ element->GetInputValue(&bed, gauss,BaseEnum); element->GetInputValue(&head, gauss,HydrologyHeadEnum); + element->GetInputValue(&sealevel, gauss,SealevelEnum); element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum); @@ -613,5 +618,5 @@ //From bed and thickness, compute effective pressure when drag is viscous: pressure_ice = rho_ice*gravity*thickness; - pressure_water = rho_water*gravity*(head-bed); + pressure_water = rho_water*gravity*(head-bed+sealevel); Neff=pressure_ice-pressure_water; if(Neff<0.) Neff=0.;