Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16321) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 16322) @@ -119,5 +119,5 @@ else{ enthalpy_core(femmodel); - PostprocessingEnthalpyx(femmodel); + PostprocessingEnthalpyx(femmodel); } #else Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16321) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16322) @@ -4000,6 +4000,6 @@ enthalpy_input->GetInputValue(&enthalpy, gauss); pressure_input->GetInputValue(&pressure, gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); - D_scalar_conduct=gauss->weight*Jdet*kappa; + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); _assert_(kappa>0.); + D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice; if(reCast(dt)) D_scalar_conduct=D_scalar_conduct*dt; @@ -4451,5 +4451,6 @@ scalar_def=phi/rho_ice*Jdet*gauss->weight; if(reCast(dt)) scalar_def=scalar_def*dt; - + + /*TODO: add -beta*laplace T_m(p)*/ for(i=0;ivalues[i]+=scalar_def*L[i]; @@ -4470,5 +4471,6 @@ enthalpypicard_input->GetInputValue(&enthalpypicard, gauss); kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); - tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); + kappa/=rho_ice; + tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa); for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]); @@ -4554,5 +4556,5 @@ IssmDouble Jdet2d,dt; IssmDouble rho_ice,heatcapacity,geothermalflux_value; - IssmDouble basalfriction,alpha2,vx,vy; + IssmDouble basalfriction,alpha2,vx,vy,vz; IssmDouble scalar,enthalpy,enthalpyup; IssmDouble pressure,pressureup; @@ -4600,9 +4602,27 @@ enthalpy_input->GetInputValue(&enthalpy,gauss); pressure_input->GetInputValue(&pressure,gauss); - if(enthalpy>=matpar->PureIceEnthalpy(pressure)){ + watercolumn_input->GetInputValue(&watercolumn,gauss); + + if((watercolumn==0.) && (enthalpyPureIceEnthalpy(pressure))){ + /* the above check is equivalent to + NOT ((watercolumn>0.) AND (enthalpyGetInputValue(&geothermalflux_value,gauss); + + 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)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + if(reCast(dt)) scalar=dt*scalar; + for(i=0;ivalues[i]+=scalar*basis[i]; + } + else if(enthalpy>=matpar->PureIceEnthalpy(pressure)){ + /* check positive thickness of temperate basal ice layer */ enthalpy_input->GetInputValue(&enthalpyup,gaussup); pressure_input->GetInputValue(&pressureup,gaussup); if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){ - // temperate ice has positive thickness: grad enthalpy*n=0. + // TODO: temperate ice has positive thickness: grad enthalpy*n=0. } else{ @@ -4610,20 +4630,6 @@ } } - else{ - watercolumn_input->GetInputValue(&watercolumn,gauss); - if(watercolumn==0.){ - /*add geothermal heat flux and basal friction*/ - geothermalflux_input->GetInputValue(&geothermalflux_value,gauss); - friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); - - scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); - if(reCast(dt)) scalar=dt*scalar; - - for(i=0;ivalues[i]+=scalar*basis[i]; - } - else{ /*do nothing (water layer acts as insulation)*/ } + else { + // base cold, but watercolumn positive. Set base to h_pmp. } } @@ -5011,6 +5017,7 @@ case LliboutryDuvalEnum: B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, - (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0, - material->GetN()); + (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0, + material->GetN(), + matpar->GetBeta(), matpar->GetReferenceTemperature(), matpar->GetHeatCapacity(), matpar->GetLatentHeat()); for(i=0;imaterial->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum)); @@ -5032,8 +5039,8 @@ /*Intermediary*/ - bool isenthalpy,isdynamicbasalspc,istemperatelayer; + bool isenthalpy,isdynamicbasalspc,setspc; int numindices, numindicesup; IssmDouble h_pmp,pressure, pressureup; - IssmDouble enthalpy, enthalpyup; + IssmDouble enthalpy, enthalpyup; int *indices = NULL, *indicesup = NULL; @@ -5045,42 +5052,37 @@ if(!isenthalpy) return; parameters->FindParam(&isdynamicbasalspc,ThermalIsdynamicbasalspcEnum); - isdynamicbasalspc = true; // TODO: remove before release if(!isdynamicbasalspc) return; /*Fetch indices of basal & surface nodes for this finite element*/ BasalNodeIndices(&numindices,&indices,this->element_type); - SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); - _assert_(numindices==numindicesup); + SurfaceNodeIndices(&numindicesup,&indicesup,this->element_type); + _assert_(numindices==numindicesup); /*Get parameters and inputs: */ Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); - Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); - - /*if there is a temperate layer of positive thickness, set enthalpy=h_pmp at that node*/ + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + + /*if there is a temperate layer of zero thickness, set spc enthalpy=h_pmp at that node*/ GaussPenta* gauss=new GaussPenta(); - GaussPenta* gaussup=new GaussPenta(); + GaussPenta* gaussup=new GaussPenta(); for(int i=0;iGaussNode(this->element_type,indices[i]); - gaussup->GaussNode(this->element_type,indicesup[i]); + gaussup->GaussNode(this->element_type,indicesup[i]); /*Check wether there is a temperate layer at the base or not */ - /*check if node is temperate, if not, continue*/ - enthalpy_input->GetInputValue(&enthalpy, gauss); + /*check if node is temperate, if not, continue*/ + enthalpy_input->GetInputValue(&enthalpy, gauss); pressure_input->GetInputValue(&pressure, gauss); - if (enthalpyPureIceEnthalpy(pressure)){ - continue; - } - /*check if upper node is temperate. - if yes, then we have a temperate layer of positive thickness and - reset the spc. - if not, apply dirichlet BC.*/ - enthalpy_input->GetInputValue(&enthalpyup, gaussup); - pressure_input->GetInputValue(&pressureup, gaussup); - istemperatelayer = false; - if (enthalpyup>=matpar->PureIceEnthalpy(pressureup)) - istemperatelayer=true; - - /*Add Dirichlet constraint to this node if there is no layer of temperate ice with positive thickness*/ - if(!istemperatelayer){ + setspc = false; + if (enthalpy>=matpar->PureIceEnthalpy(pressure)){ + /*check if upper node is temperate, too. + if yes, then we have a temperate layer of positive thickness and reset the spc. + if not, apply dirichlet BC.*/ + enthalpy_input->GetInputValue(&enthalpyup, gaussup); + pressure_input->GetInputValue(&pressureup, gaussup); + setspc=(enthalpyupPureIceEnthalpy(pressureup))?true:false; + } + + if (setspc) { /*Calculate enthalpy at pressure melting point */ h_pmp=matpar->PureIceEnthalpy(pressure); @@ -5088,13 +5090,15 @@ nodes[indices[i]]->ApplyConstraint(1,h_pmp); } - else { - /*remove spc*/ - nodes[indices[i]]->DofInFSet(0); - } + else { + /*remove spc*/ + nodes[indices[i]]->DofInFSet(0); + } } /*Free ressources:*/ xDelete(indices); - xDelete(indicesup); + xDelete(indicesup); + delete gauss; + delete gaussup; } /*}}}*/ @@ -5102,4 +5106,5 @@ void Penta::ComputeBasalMeltingrate(void){ /*Calculate the basal melt rates of the enthalpy model after Aschwanden 2012*/ + /* melting rate is positive when melting, negative when refreezing*/ /* Intermediaries */ @@ -5160,12 +5165,10 @@ /*ensure that no ice is at TPureIceEnthalpy(pressure[iv]); - //meltingrate_enthalpy=0.; } 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.; - checkpositivethickness=false; - } - else {/*do nothing, go to next check*/} + checkpositivethickness=false; // cold base, skip next test + } + else {/*we have a temperate base, go to next test*/} if(checkpositivethickness){ @@ -5176,5 +5179,5 @@ else{ enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy[iv],pressure[iv]); _assert_(kappa>0.); for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; } Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16321) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16322) @@ -378,8 +378,8 @@ IssmDouble Matpar::GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){ if(enthalpyelements->Size();i++){ + /*drain excess water fraction: */ + for (i=0;ielements->Size();i++){ element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); element->DrainWaterfraction(); } -*/ /*Update basal dirichlet BCs for enthalpy: */ Index: /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 16321) +++ /issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp (revision 16322) @@ -8,6 +8,6 @@ /* get ice stiffness B from enthalpy, pressure and flow law exponent*/ -IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n){ - /*Use parameterization for the rheology: Aschwanden 2012 +IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat){ + /*Use parameterization for the rheology: Grewe/Blatter 2009, Aschwanden 2012 * * A(H,p) = A0 exp(-Q/RT(H,p)), if H < H_s(p) @@ -22,5 +22,5 @@ * H_s(p) = c_i (Tpmp - Tref) * - * Tpmp = T - beta p; + * Tpmp = T - betaCC p; * * A0 constant of proportionality @@ -35,21 +35,22 @@ * Convert A to B : B = A^(-1/n) */ - /*Some physical constants (Aschwanden 2012)*/ - /*TODO: get those constants from model*/ - IssmDouble beta=7.9*pow(10.,-8.); - IssmDouble R=8.314; - IssmDouble heatcapacity=2009; // J/kg/K - IssmDouble Tref=253.15; - IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012 + /*check feasibility*/ + _assert_(pressure>0); + _assert_(n>0); + _assert_(betaCC>0); + _assert_(referencetemperature>0); + _assert_(heatcapacity>0); + _assert_(latentheat>0); + + /*Some physical constants*/ + IssmDouble R=8.314; /*Intermediaries*/ IssmDouble A,B,Tstar,Tpmp,H_sp,waterfraction; - _assert_(pressure>0); - _assert_(enthalpy>0); - Tpmp=273.15-beta*pressure; - H_sp=heatcapacity*(Tpmp - Tref); + Tpmp=273.15-betaCC*pressure; + H_sp=heatcapacity*(Tpmp - referencetemperature); if (enthalpy < H_sp){ - Tstar = Tref + enthalpy/heatcapacity - beta*pressure; + Tstar = referencetemperature + enthalpy/heatcapacity - betaCC*pressure; waterfraction = 0; } @@ -75,34 +76,36 @@ } -/*Get stiffness from temperature, waterfraction and depth*/ -IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){ - /*Use parameterization for the rheology: Greve and Blatter 2009 - * get enthalpy from temperature and water fraction, - * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */ - - IssmDouble rho_ice=910; // kg/m^3 - IssmDouble g=9.81; //kg*m/s^2 - IssmDouble heatcapacity=2009; // J/kg/K - IssmDouble Tref=253.15; - IssmDouble beta=7.9*pow(10.,-8.); - IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012 - IssmDouble Tstar, enthalpy, pressure, B; - _assert_(temperature>0); - _assert_(waterfraction>0); - _assert_(depth>0); - - /*get pressure*/ - pressure= rho_ice*g*depth; - Tstar=temperature-beta*pressure; // TODO: check whether plus or minus - /*get enthalpy*/ - if (Tstar < 273.15){ - enthalpy=heatcapacity*(Tstar - Tref); - } - else{ - enthalpy=heatcapacity*(273.15 - Tref) + waterfraction*latentheat; - } - - B=LliboutryDuval(enthalpy, pressure, n); - - return B; -} +// /*Get stiffness from temperature, waterfraction and depth*/ +// IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){ +// /*Use parameterization for the rheology: Greve and Blatter 2009 +// * get enthalpy from temperature and water fraction, +// * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */ +// +// /*TODO: update params from model*/ +// IssmDouble rho_ice=910; // kg/m^3 +// IssmDouble g=9.81; //kg*m/s^2 +// IssmDouble heatcapacity=2009; // J/kg/K +// IssmDouble referencetemperature=253.15; +// IssmDouble betaCC=7.9*pow(10.,-8.); +// IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012 +// +// IssmDouble Tstar, enthalpy, pressure, B; +// _assert_(temperature>0); +// _assert_(waterfraction>0); +// _assert_(depth>0); +// +// /*get pressure*/ +// pressure= rho_ice*g*depth; +// Tstar=temperature-betaCC*pressure; // TODO: check whether plus or minus +// /*get enthalpy*/ +// if (Tstar < 273.15){ +// enthalpy=heatcapacity*(Tstar - referencetemperature); +// } +// else{ +// enthalpy=heatcapacity*(273.15 - referencetemperature) + waterfraction*latentheat; +// } +// +// B=LliboutryDuval(enthalpy, pressure, n, betaCC, referencetemperature, heatcapacity, latentheat); +// +// return B; +// } Index: /issm/trunk-jpl/src/c/shared/Elements/elements.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 16321) +++ /issm/trunk-jpl/src/c/shared/Elements/elements.h (revision 16322) @@ -10,6 +10,6 @@ IssmDouble Paterson(IssmDouble temperature); IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n); -IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n); -IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); +IssmDouble LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure, IssmDouble n, IssmDouble betaCC, IssmDouble referencetemperature, IssmDouble heatcapacity, IssmDouble latentheat); +// IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n); IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures, IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, Index: /issm/trunk-jpl/src/m/classes/thermal.m =================================================================== --- /issm/trunk-jpl/src/m/classes/thermal.m (revision 16321) +++ /issm/trunk-jpl/src/m/classes/thermal.m (revision 16322) @@ -56,7 +56,10 @@ md = checkfield(md,'thermal.spctemperature(find(md.thermal.spctemperature(1:md.mesh.numberofvertices,:)~=NaN))','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*replicate(pos),'message','spctemperature should be below the adjusted melting point'); md = checkfield(md,'thermal.isenthalpy','numel',[1],'values',[0 1]); - md = checkfield(md,'thermal.isdynamicbasalspc','numel',[1],'values',[0 1]); - end - end % }}} + md = checkfield(md,'thermal.isdynamicbasalspc','numel', [1],'values',[0 1]); + if(md.thermal.isenthalpy) + md = checkfield(md,'thermal.isdynamicbasalspc','numel', [1],'values',[1], 'message',['for enthalpy run thermal.isdynamicbasalspc should be 1']); + end + end + end % }}} function disp(obj) % {{{ disp(sprintf(' Thermal solution parameters:')); Index: /issm/trunk-jpl/test/NightlyRun/test121.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test121.m (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test121.m (revision 16322) @@ -12,4 +12,5 @@ md.transient.isgroundingline=0; md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; md=solve(md,TransientSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test121.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test121.py (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test121.py (revision 16322) @@ -21,4 +21,5 @@ md.transient.isgroundingline=False md.thermal.isenthalpy=1 +md.thermal.isdynamicbasalspc=1 md=solve(md,TransientSolutionEnum()) Index: /issm/trunk-jpl/test/NightlyRun/test122.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test122.m (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test122.m (revision 16322) @@ -7,4 +7,5 @@ md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1); md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; md.thermal.stabilization=2; md.cluster=generic('name',oshostname(),'np',3); Index: /issm/trunk-jpl/test/NightlyRun/test122.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test122.py (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test122.py (revision 16322) @@ -16,4 +16,5 @@ md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,1)) md.thermal.isenthalpy=1 +md.thermal.isdynamicbasalspc=1 md.thermal.stabilization=2 md.cluster=generic('name',oshostname(),'np',3) Index: /issm/trunk-jpl/test/NightlyRun/test326.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test326.m (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test326.m (revision 16322) @@ -12,4 +12,5 @@ md.transient.isgroundingline=0; md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; md=solve(md,TransientSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test326.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test326.py (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test326.py (revision 16322) @@ -22,4 +22,5 @@ md.transient.isgroundingline=False md.thermal.isenthalpy=1 +md.thermal.isdynamicbasalspc=1 md=solve(md,TransientSolutionEnum()) Index: /issm/trunk-jpl/test/NightlyRun/test327.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test327.m (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test327.m (revision 16322) @@ -10,4 +10,5 @@ md.thermal.spctemperature(find(md.mesh.vertexonsurface))=272.; md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; md.basalforcings.geothermalflux(:)=5.; md=solve(md,TransientSolutionEnum()); Index: /issm/trunk-jpl/test/NightlyRun/test327.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test327.py (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test327.py (revision 16322) @@ -20,4 +20,5 @@ md.thermal.spctemperature[numpy.nonzero(md.mesh.vertexonsurface)[0]]=272. md.thermal.isenthalpy=1 +md.thermal.isdynamicbasalspc=1 md.basalforcings.geothermalflux[:]=5. md=solve(md,TransientSolutionEnum()) Index: /issm/trunk-jpl/test/NightlyRun/test431.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test431.m (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test431.m (revision 16322) @@ -7,4 +7,5 @@ md.timestepping.time_step=0.; md.thermal.isenthalpy=1; +md.thermal.isdynamicbasalspc=1; md.initialization.waterfraction=zeros(md.mesh.numberofvertices,1); md.initialization.watercolumn=zeros(md.mesh.numberofvertices,1); Index: /issm/trunk-jpl/test/NightlyRun/test431.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test431.py (revision 16321) +++ /issm/trunk-jpl/test/NightlyRun/test431.py (revision 16322) @@ -17,4 +17,5 @@ md.timestepping.time_step=0. md.thermal.isenthalpy=1 +md.thermal.isdynamicbasalspc=1 md.initialization.waterfraction=numpy.zeros((md.mesh.numberofvertices,1)) md.initialization.watercolumn=numpy.zeros((md.mesh.numberofvertices,1))