Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16033) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16034) @@ -4768,7 +4768,5 @@ for(i=0;ivalues[i]+=scalar*basis[i]; } - else{ - /*do nothing (water layer acts as insulation)*/ - } + else{ /*do nothing (water layer acts as insulation)*/ } } } @@ -5308,14 +5306,16 @@ IssmDouble xyz_list[NUMVERTICES][3]; IssmDouble xyz_list_tria[NUMVERTICES2D][3]; - IssmDouble heatflux, geothermalflux_value,kappa; + IssmDouble heatflux,kappa; IssmDouble vec_heatflux[3]; IssmDouble normal_base[3], d1enthalpy[3]; - IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES], meltingrate_enthalpy; + IssmDouble basalmeltingrate[NUMVERTICES], watercolumn[NUMVERTICES]; IssmDouble enthalpy[NUMVERTICES], enthalpyup; - IssmDouble pressure, pressureup; + IssmDouble pressure[NUMVERTICES], pressureup; IssmDouble temperature, waterfraction; IssmDouble latentheat, rho_ice; IssmDouble basalfriction, alpha2; - IssmDouble vx,vy,vz,dt,connectivity; + IssmDouble vx[NUMVERTICES],vy[NUMVERTICES],vz[NUMVERTICES]; + IssmDouble geothermalflux[NUMVERTICES]; + IssmDouble dt,meltingrate_enthalpy; Friction *friction = NULL; @@ -5330,12 +5330,13 @@ latentheat=matpar->GetLatentHeat(); rho_ice=this->matpar->GetRhoIce(); - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - Input* basalmeltingrate_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basalmeltingrate_input); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); + GetInputListOnVertices(&vx[0],VxEnum); + GetInputListOnVertices(&vy[0],VyEnum); + GetInputListOnVertices(&vz[0],VzEnum); + GetInputListOnVertices(&enthalpy[0],EnthalpyEnum); + GetInputListOnVertices(&pressure[0],PressureEnum); + GetInputListOnVertices(&watercolumn[0],WatercolumnEnum); + GetInputListOnVertices(&basalmeltingrate[0],BasalforcingsMeltingRateEnum); + GetInputListOnVertices(&geothermalflux[0],BasalforcingsGeothermalfluxEnum); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); for(i=0;iGaussVertex(iv); gaussup->GaussVertex(iv+3); - checkpositivethickness=true; - watercolumn_input->GetInputValue(&watercolumn[iv], gauss); - basalmeltingrate_input->GetInputValue(&basalmeltingrate[iv],gauss); - enthalpy_input->GetInputValue(&enthalpy[iv], gauss); - pressure_input->GetInputValue(&pressure, gauss); /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/ meltingrate_enthalpy=0.; - if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure))){ + if((watercolumn[iv]>0.) && (enthalpy[iv]PureIceEnthalpy(pressure[iv]))){ /*ensure that no ice is at TPureIceEnthalpy(pressure); + enthalpy[iv]=matpar->PureIceEnthalpy(pressure[iv]); //meltingrate_enthalpy=0.; } - else if(enthalpy[iv]PureIceEnthalpy(pressure)){ + else if(enthalpy[iv]PureIceEnthalpy(pressure[iv])){ /*cold base: set q*n=q_geo*n+frictionheating as Neumann BC in Penta::CreatePVectorEnthalpySheet*/ meltingrate_enthalpy=0.; @@ -5375,11 +5372,9 @@ /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */ istemperatelayer=false; - enthalpy_input->GetInputValue(&enthalpyup, gaussup); - pressure_input->GetInputValue(&pressureup, gaussup); - if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)) istemperatelayer=true; + if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true; if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; else{ enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]); for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; } @@ -5389,20 +5384,14 @@ heatflux=0.; for(i=0;i<3;i++) heatflux+=(vec_heatflux[i])*normal_base[i]; - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); /*basal friction*/ friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0)); - - matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure); - meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux_value))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent + basalfriction=alpha2*(pow(vx[iv],2.0)+pow(vy[iv],2.0)+pow(vz[iv],2.0)); + + matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]); + meltingrate_enthalpy=(basalfriction-(heatflux-geothermalflux[iv]))/((1-waterfraction)*latentheat*rho_ice); // m/yr water equivalent } /*Update water column, basal meltingrate*/ - connectivity=IssmDouble(vertices[iv]->Connectivity()); - meltingrate_enthalpy/=connectivity; // divide meltingrate_enthalpy by connectivity to address multiple iterations over vertex basalmeltingrate[iv]+=meltingrate_enthalpy; this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); @@ -5431,10 +5420,9 @@ IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES]; IssmDouble latentheat, dt; - IssmDouble connectivity; GaussPenta* gauss; - Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); /*Check wether enthalpy is activated*/ @@ -5453,7 +5441,5 @@ /*drain water fraction & update enthalpy*/ - // TODO: replace connectivity-wise draining by draining once per node per timestep - connectivity=IssmDouble(vertices[iv]->Connectivity()); - waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt)/connectivity; + waterfraction[iv]-=DrainageFunctionWaterfraction(waterfraction[iv], dt); matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]); }