Index: /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp (revision 21216) @@ -58,4 +58,16 @@ iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum); iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum); + iomodel->FetchDataToInput(elements,"md.smb.isInitialized",SmbIsInitializedEnum); + iomodel->FetchDataToInput(elements,"md.smb.isrestart",SmbIsrestartEnum); + iomodel->FetchDataToInput(elements,"md.smb.Dzrst",SmbDzrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Drst",SmbDrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Rerst",SmbRerstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Gdnrst",SmbGdnrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Gsprst",SmbGsprstEnum); + iomodel->FetchDataToInput(elements,"md.smb.ECrst",SmbECrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Wrst",SmbWrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Arst",SmbArstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Trst",SmbTrstEnum); + iomodel->FetchDataToInput(elements,"md.smb.Sizerst",SmbSizerstEnum); break; case SMBpddEnum: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 21216) @@ -1296,87 +1296,102 @@ /*}}}*/ void Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/ - - /*Intermediaries*/ - int i,t; - IssmDouble time; - - /*Branch on type of vector: nodal or elementary: */ - if(vector_type==1){ //nodal vector - - int numvertices = this->GetNumberOfVertices(); - int *vertexids = xNew(numvertices); - IssmDouble *values = xNew(numvertices); - - /*Recover vertices ids needed to initialize inputs*/ - _assert_(iomodel->elements); - for(i=0;i(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab - } - - /*Are we in transient or static? */ - if(M==iomodel->numberofvertices){ - for(i=0;iAddInput(vector_enum,values,P1Enum); - } - else if(M==iomodel->numberofvertices+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; - default: _error_("Not implemented yet"); - } - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - - xDelete(values); - xDelete(vertexids); - } - else if(vector_type==2){ //element vector - - IssmDouble value; - - /*Are we in transient or static? */ - if(M==iomodel->numberofelements){ - if (code==5){ //boolean - this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()]))); - } - else if (code==6){ //integer - this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()]))); - } - else if (code==7){ //IssmDouble - this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()])); - } - else _error_("could not recognize nature of vector from code " << code); - } - else if(M==iomodel->numberofelements+1){ - /*create transient input: */ - IssmDouble* times = xNew(N); - for(t=0;tSid()+t]; - switch(this->ObjectEnum()){ - case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; - case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; - case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; - default: _error_("Not implemented yet"); - } - } - this->inputs->AddInput(transientinput); - xDelete(times); - } - else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); - } - else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); + + /*Intermediaries*/ + int i,t; + IssmDouble time; + + /*Branch on type of vector: nodal or elementary: */ + if(vector_type==1){ //nodal vector + + int numvertices = this->GetNumberOfVertices(); + int *vertexids = xNew(numvertices); + IssmDouble *values = xNew(numvertices); + + /*Recover vertices ids needed to initialize inputs*/ + _assert_(iomodel->elements); + for(i=0;i(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab + } + + /*Are we in transient or static? */ + if(M==iomodel->numberofvertices){ + for(i=0;iAddInput(vector_enum,values,P1Enum); + } + else if(M==iomodel->numberofvertices+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + + xDelete(values); + xDelete(vertexids); + } + else if(vector_type==2){ //element vector + + IssmDouble value; + + /*Are we in transient or static? */ + if(M==iomodel->numberofelements){ + if (code==5){ //boolean + this->inputs->AddInput(new BoolInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==6){ //integer + this->inputs->AddInput(new IntInput(vector_enum,reCast(vector[this->Sid()]))); + } + else if (code==7){ //IssmDouble + this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()])); + } + else _error_("could not recognize nature of vector from code " << code); + } + else if(M==iomodel->numberofelements+1){ + /*create transient input: */ + IssmDouble* times = xNew(N); + for(t=0;tSid()+t]; + switch(this->ObjectEnum()){ + case TriaEnum: transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break; + case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break; + case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break; + default: _error_("Not implemented yet"); + } + } + this->inputs->AddInput(transientinput); + xDelete(times); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + else if(vector_type==3){ //element vector + + IssmDouble value; + + /*For right now we are static */ + if(M==iomodel->numberofelements){ + /*create transient input: */ + IssmDouble* layers = xNewZeroInit(N);; + for(t=0;tSid()+t]; + DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N); + this->inputs->AddInput(arrayinput); + xDelete(layers); + } + else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long"); + } + else _error_("Cannot add input for vector type " << vector_type << " (not supported)"); } /*}}}*/ @@ -2231,6 +2246,7 @@ /*Intermediary variables: {{{*/ - bool isinitialized=false; - IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; + IssmDouble isinitialized; + IssmDouble isrestart; + IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; IssmDouble Tmean; IssmDouble C; @@ -2251,7 +2267,8 @@ IssmDouble initMass; IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd; + IssmDouble sumdz_add; IssmDouble sumMass,dMass; bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; - IssmDouble init_scaling=1.0; + IssmDouble init_scaling; /*}}}*/ /*Output variables:{{{ */ @@ -2270,4 +2287,15 @@ IssmDouble R; IssmDouble mAdd; + IssmDouble dz_add; + + IssmDouble* dzrst=NULL; + IssmDouble* drst = NULL; + IssmDouble* rerst = NULL; + IssmDouble* gdnrst = NULL; + IssmDouble* gsprst = NULL; + IssmDouble* Wrst = NULL; + IssmDouble* arst = NULL; + IssmDouble* Trst = NULL; + int m; int count=0; @@ -2285,4 +2313,5 @@ parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ + parameters->FindParam(&yts,ConstantsYtsEnum); parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ parameters->FindParam(&aIdx,SmbAIdxEnum); @@ -2321,6 +2350,7 @@ Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input); Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input); - Input* isinitialized_input=this->inputs->GetInput(SmbIsInitializedEnum); - + Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input); + Input* isrestart_input=this->GetInput(SmbIsrestartEnum); _assert_(isrestart_input); + /*Retrieve input values:*/ Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); @@ -2336,31 +2366,77 @@ Tz_input->GetInputValue(&Tz,gauss); Vz_input->GetInputValue(&Vz,gauss); + isinitialized_input->GetInputValue(&isinitialized,gauss); + isrestart_input->GetInputValue(&isrestart,gauss); /*}}}*/ - /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ - if(!isinitialized_input){ - - if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n"); - GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); - //if(this->Sid()==1) for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;iAddInput(new BoolInput(SmbIsInitializedEnum,true)); - } - else{ + /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ + if(isinitialized){ + if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n"); + if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n"); + GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); + //if(this->Sid()==1) for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;i(m); for(int i=0;iAddInput(new DoubleInput(SmbIsInitializedEnum,false)); + } + else if(isrestart){ //Retrieve the snow properties from previous run + if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n"); + + /*Recover inputs: */ + DoubleArrayInput* dzrst_input= dynamic_cast(this->GetInput(SmbDzrstEnum)); _assert_(dzrst_input); + DoubleArrayInput* drst_input= dynamic_cast(this->GetInput(SmbDrstEnum));_assert_(drst_input); + DoubleArrayInput* rerst_input= dynamic_cast(this->GetInput(SmbRerstEnum));_assert_(rerst_input); + DoubleArrayInput* gdnrst_input= dynamic_cast(this->GetInput(SmbGdnrstEnum));_assert_(gdnrst_input); + DoubleArrayInput* gsprst_input= dynamic_cast(this->GetInput(SmbGsprstEnum));_assert_(gsprst_input); + DoubleInput* ECrst_input= dynamic_cast(this->GetInput(SmbECrstEnum));_assert_(ECrst_input); + DoubleArrayInput* Wrst_input= dynamic_cast(this->GetInput(SmbWrstEnum));_assert_(Wrst_input); + DoubleArrayInput* arst_input= dynamic_cast(this->GetInput(SmbArstEnum));_assert_(arst_input); + DoubleArrayInput* Trst_input= dynamic_cast(this->GetInput(SmbTrstEnum));_assert_(Trst_input); + + /*Recover arrays: */ + dzrst_input->GetValues(&dzrst,&m); + drst_input->GetValues(&drst,&m); + rerst_input->GetValues(&rerst,&m); + gdnrst_input->GetValues(&gdnrst,&m); + gsprst_input->GetValues(&gsprst,&m); + ECrst_input->GetInputValue(&EC); + Wrst_input->GetValues(&Wrst,&m); + arst_input->GetValues(&arst,&m); + Trst_input->GetValues(&Trst,&m); + + /*Retrive the correct value of m (without the zeroes at the end)*/ + Input* Sizerst_input=this->GetInput(SmbSizerstEnum); _assert_(Sizerst_input); + Sizerst_input->GetInputValue(&m); + + dz = xNewZeroInit(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;i(m);for(int i=0;iAddInput(new DoubleInput(SmbIsrestartEnum,false));//*CL* add + } + else{ + /*Recover inputs: */ DoubleArrayInput* dz_input= dynamic_cast(this->GetInput(SmbDzEnum)); _assert_(dz_input); @@ -2395,4 +2471,5 @@ // initialize cumulative variables sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; + sumdz_add=0; //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. @@ -2402,5 +2479,5 @@ /*Start loop: */ count=1; - for (t=time;t<=time+dt;t=t+smb_dt){ + for (t=time;tSid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n"); @@ -2442,5 +2519,5 @@ /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ - if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin,this->Sid()); + if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,this->Sid()); /*Allow non-melt densification and determine compaction [m]*/ @@ -2479,4 +2556,5 @@ sumP = P + sumP; sumEC = sumEC + EC; // evap (-)/cond(+) + sumdz_add=dz_add+sumdz_add; //*CL* /*Calculate total system mass:*/ @@ -2499,8 +2577,6 @@ /*increase counter:*/ count++; - - } //for (t=time;t<=time+dt;t=t+smb_dt) - - + } //for (t=time;tAddInput(new DoubleArrayInput(SmbDzEnum,dz,m)); @@ -2510,13 +2586,12 @@ this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m)); this->AddInput(new DoubleArrayInput(SmbTEnum,T,m)); - this->AddInput(new DoubleInput(SmbECEnum,EC)); + this->AddInput(new DoubleInput(SmbECEnum,sumEC/yts)); this->AddInput(new DoubleArrayInput(SmbWEnum,W,m)); this->AddInput(new DoubleArrayInput(SmbAEnum,a,m)); - this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/rho_water/dt)); - this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/rho_water/dt)); - this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/rho_water/dt)); - this->AddInput(new DoubleInput(SmbCondensationEnum,sumEC/rho_water/dt)); - - + this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/yts)); + this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/yts)); + this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/yts)); + this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts)); + this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts)); /*Free allocations:{{{*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 21216) @@ -1183,5 +1183,5 @@ *pm=m; } /*}}}*/ -void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, int sid){ /*{{{*/ +void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid){ /*{{{*/ //// MELT ROUTINE @@ -1218,4 +1218,24 @@ IssmDouble X; IssmDouble Wi; + + IssmDouble Ztot; + IssmDouble T_bot; + IssmDouble m_bot; + IssmDouble dz_bot; + IssmDouble d_bot; + IssmDouble W_bot; + IssmDouble a_bot; + IssmDouble re_bot; + IssmDouble gdn_bot; + IssmDouble gsp_bot; + bool top=false; + + IssmDouble* Zcum=NULL; + IssmDouble* dzMin2=NULL; + IssmDouble zY2=1.025; + bool lastCellFlag; + int X1=0; + int X2=0; + int D_size; int i; @@ -1223,4 +1243,5 @@ /*outputs:*/ IssmDouble mAdd; + IssmDouble dz_add; IssmDouble Rsum; IssmDouble* T=*pT; @@ -1264,4 +1285,5 @@ mAdd = 0; // mass added/removed to/from base of model [kg] addE = 0; // energy added/removed to/from base of model [J] + dz_add=0; // thickness of the layer added/removed to/from base of model [m] // calculate temperature excess above 0 deg C @@ -1301,36 +1323,37 @@ } - //// MELT, PERCOLATION AND REFREEZE - - // run melt algorithm if there is melt water or excess pore water - if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){ - - // _printf_(""MELT OCCURS"); - // check to see if thermal energy exceeds energy to melt entire cell - // if so redistribute temperature to lower cells (temperature surplus) - // (maximum T of snow before entire grid cell melts is a constant - // LF/CI = 159.1342) - surpT=xNew(n); for(int i=0;i 0 ){ - // _printf_("T Surplus"); - // calculate surplus energy - surpE=xNew(n); for(int i=0;i 0){ - // use surplus energy to increase the temperature of lower cell - T[i+1] = surpE[i] * m[i+1]/CI + T[i+1]; - surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342)); - surpE[i+1] = surpT[i+1] * CI / m[i+1]; - - // adjust current cell properties (again 159.1342 is the max T) - T[i] = CtoK + 159.1342; - surpE[i] = 0; - i = i + 1; - } - // recalculate temperature excess above 0 deg C - for(int i=0;i 0) || (cellsum(exsW,n) > 0)){ + // _printf_(""MELT OCCURS"); + // check to see if thermal energy exceeds energy to melt entire cell + // if so redistribute temperature to lower cells (temperature surplus) + // (maximum T of snow before entire grid cell melts is a constant + // LF/CI = 159.1342) + surpT=xNew(n); for(int i=0;i 0 ){ + // _printf_("T Surplus"); + // calculate surplus energy + surpE=xNew(n); for(int i=0;i 0){ + // use surplus energy to increase the temperature of lower cell + T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];//*CL* + + exsT[i+1] = fmax(0, T[i+1] - CtoK) + exsT[i+1]; + T[i+1] = fmin(CtoK, T[i+1]); + + surpT[i+1] = fmax(0, exsT[i+1] - 159.1342); + surpE[i+1] = surpT[i+1] * CI * m[i+1];//*CL* + + // adjust current cell properties (again 159.1342 is the max T) + exsT[i] = 159.1342; + surpE[i] = 0; + i = i + 1; + } + } // convert temperature excess to melt [kg] @@ -1469,35 +1492,92 @@ xDelete(R); } + + //Merging of cells as they are burried under snow. + Zcum=xNew(n); + dzMin2=xNew(n); + + Zcum[0]=dz[0]; // Compute a cumulative depth vector + + for (int i=1;i=0;i--){ - if(dz[i] dzMin * 2') - // adjust variables as a linearly weighted function of mass - IssmDouble m_new = m[i] + m[i+1]; - T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new; - a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new; - re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new; - gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new; - gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new; - - // merge with underlying grid cell and delete old cell - dz [i+1] = dz[i] + dz[i+1]; // combine cell depths - d[i+1] = m_new / dz[i+1]; // combine top densities - W[i+1] = W[i+1] + W[i]; // combine liquid water - m[i+1] = m_new; // combine top masses - - // set cell to 99999 for deletion - m[i] = 99999; - } - } + + //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell) + if(X==n-1){ + lastCellFlag = true; + X = X-1; + } + else{ + lastCellFlag = false; + } + + for (int i = 0; i<=X;i++){ + if (dz [i] < dzMin2[i]){ // merge top cells + // _printf_("dz > dzMin * 2') + // adjust variables as a linearly weighted function of mass + IssmDouble m_new = m[i] + m[i+1]; + T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new; + a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new; + re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new; + gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new; + gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new; + + // merge with underlying grid cell and delete old cell + dz [i+1] = dz[i] + dz[i+1]; // combine cell depths + d[i+1] = m_new / dz[i+1]; // combine top densities + W[i+1] = W[i+1] + W[i]; // combine liquid water + m[i+1] = m_new; // combine top masses + + // set cell to 99999 for deletion + m[i] = 99999; + } + } + + //If last cell has to be merged + if(lastCellFlag){ + //find closest cell to merge with + for(int i=n-2;i>=0;i--){ + if(m[i]!=99999){ + X2=i; + X1=n-1; + break; + } + } + + // adjust variables as a linearly weighted function of mass + IssmDouble m_new = m[X2] + m[X1]; + T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new; + a[X1] = (a[X2]*m[i] + a[X1]*m[X1]) / m_new; + re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new; + gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new; + gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new; + + // merge with underlying grid cell and delete old cell + dz [X1] = dz[X2] + dz[X1]; // combine cell depths + d[X1] = m_new / dz[X1]; // combine top densities + W[X1] = W[X1] + W[X2]; // combine liquid water + m[X1] = m_new; // combine top masses + + // set cell to 99999 for deletion + m[X2] = 99999; + } // delete combined cells @@ -1551,39 +1631,65 @@ } - //// CORRECT FOR TOTAL MODEL DEPTH - // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT - // INTERPRETATION - - // // calculate total model depth - // z = sum(dz); - // - // if (z < zMin){ // check if model is too shallow - // _printf_("z < zMin') - // // mass and energy to be added - // mAdd = m(end) + W(end); - // addE = T(end) * m(end) * CI; - // - // // add a grid cell of the same size and temperature to the bottom - // dz = [dz; dz(end)]; - // T = [T; T(end)]; - // W = [W; W(end)]; - // m = [m; m(end)]; - // d = [d; d(end)]; - // a = [a; a(end)]; - // re = [re; re(end)]; - // gdn = [gdn; gdn(end)]; - // gsp = [gsp; gsp(end)]; - // } - // else (if z > zMax){ // check if model is too deep - // _printf_("z > zMax') - // // mass and energy loss - // mAdd = -(m(end) + W(end)); - // addE = -(T(end) * m(end) * CI); - // - // // add a grid cell of the same size and temperature to the bottom - // dz(end) = []; T(end) = []; W(end) = []; m(end) = []; - // d(end) = []; a(end) = []; re(end) = []; gdn(end) = []; - // gsp(end) = []; - // } + //// CORRECT FOR TOTAL MODEL DEPTH + // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT + // INTERPRETATION + // // calculate total model depth + Ztot = cellsum(dz,n); + + if (Ztot < zMin){ +// printf("Total depth < zMin %f \n", Ztot); + // mass and energy to be added + mAdd = m[n-1]+W[n-1]; + addE = T[n-1]*m[n-1]*CI; + + // add a grid cell of the same size and temperature to the bottom + dz_bot=dz[n-1]; + T_bot=T[n-1]; + W_bot=W[n-1]; + m_bot=m[n-1]; + d_bot=d[n-1]; + a_bot=a[n-1]; + re_bot=re[n-1]; + gdn_bot=gdn[n-1]; + gsp_bot=gsp[n-1]; + + dz_add=dz_bot; + + newcell(&dz,dz_bot,top,n); + newcell(&T,T_bot,top,n); + newcell(&W,W_bot,top,n); + newcell(&m,m_bot,top,n); + newcell(&d,d_bot,top,n); + newcell(&a,a_bot,top,n); + newcell(&re,re_bot,top,n); + newcell(&gdn,gdn_bot,top,n); + newcell(&gsp,gsp_bot,top,n); + n=n+1; + } + else if (Ztot > zMax){ +// printf("Total depth > zMax %f \n", Ztot); + // mass and energy loss + mAdd = -(m[n-1]+W[n-1]); + addE = -(T[n-1]*m[n-1]*CI); + dz_add=-(dz[n-1]); + + // add a grid cell of the same size and temperature to the bottom + D_size=n-1; + D=xNew(D_size); + + for(int i=0;i(D); + + } //// CHECK FOR MASS AND ENERGY CONSERVATION @@ -1626,5 +1732,6 @@ *pR=Rsum; *pmAdd=mAdd; - + *pdz_add=dz_add; + *pT=T; *pd=d; Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 21215) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 21216) @@ -28,5 +28,5 @@ void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, int sid); void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid); -void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, int sid); +void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, int sid); void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble rho_ice,int m, int sid); void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, int sid); Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21215) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 21216) @@ -361,4 +361,15 @@ SmbRequestedOutputsEnum, SmbIsInitializedEnum, + SmbIsrestartEnum, + SmbDzrstEnum, + SmbDrstEnum, + SmbRerstEnum, + SmbGdnrstEnum, + SmbGsprstEnum, + SmbECrstEnum, + SmbWrstEnum, + SmbArstEnum, + SmbTrstEnum, + SmbSizerstEnum, /*SMBforcing*/ SMBforcingEnum, @@ -401,5 +412,4 @@ SmbGspEnum, SmbECEnum, - SmbCondensationEnum, SmbWEnum, SmbAEnum, @@ -413,4 +423,6 @@ SmbIsdensificationEnum, SmbIsturbulentfluxEnum, + SmbDz_addEnum, + SmbM_addEnum, /*SMBpdd*/ SMBpddEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 21216) @@ -364,4 +364,15 @@ case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; case SmbIsInitializedEnum : return "SmbIsInitialized"; + case SmbIsrestartEnum : return "SmbIsrestart"; + case SmbDzrstEnum : return "SmbDzrst"; + case SmbDrstEnum : return "SmbDrst"; + case SmbRerstEnum : return "SmbRerst"; + case SmbGdnrstEnum : return "SmbGdnrst"; + case SmbGsprstEnum : return "SmbGsprst"; + case SmbECrstEnum : return "SmbECrst"; + case SmbWrstEnum : return "SmbWrst"; + case SmbArstEnum : return "SmbArst"; + case SmbTrstEnum : return "SmbTrst"; + case SmbSizerstEnum : return "SmbSizerst"; case SMBforcingEnum : return "SMBforcing"; case SmbMassBalanceEnum : return "SmbMassBalance"; @@ -402,5 +413,4 @@ case SmbGspEnum : return "SmbGsp"; case SmbECEnum : return "SmbEC"; - case SmbCondensationEnum : return "SmbCondensation"; case SmbWEnum : return "SmbW"; case SmbAEnum : return "SmbA"; @@ -414,4 +424,6 @@ case SmbIsdensificationEnum : return "SmbIsdensification"; case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux"; + case SmbDz_addEnum : return "SmbDz_add"; + case SmbM_addEnum : return "SmbM_add"; case SMBpddEnum : return "SMBpdd"; case SmbDelta18oEnum : return "SmbDelta18o"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 21216) @@ -370,7 +370,21 @@ else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; + else if (strcmp(name,"SmbIsrestart")==0) return SmbIsrestartEnum; + else if (strcmp(name,"SmbDzrst")==0) return SmbDzrstEnum; + else if (strcmp(name,"SmbDrst")==0) return SmbDrstEnum; + else if (strcmp(name,"SmbRerst")==0) return SmbRerstEnum; + else if (strcmp(name,"SmbGdnrst")==0) return SmbGdnrstEnum; + else if (strcmp(name,"SmbGsprst")==0) return SmbGsprstEnum; + else if (strcmp(name,"SmbECrst")==0) return SmbECrstEnum; + else if (strcmp(name,"SmbWrst")==0) return SmbWrstEnum; + else if (strcmp(name,"SmbArst")==0) return SmbArstEnum; + else if (strcmp(name,"SmbTrst")==0) return SmbTrstEnum; + else if (strcmp(name,"SmbSizerst")==0) return SmbSizerstEnum; else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; - else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; @@ -383,8 +397,5 @@ else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum; else if (strcmp(name,"SmbC")==0) return SmbCEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbTz")==0) return SmbTzEnum; + else if (strcmp(name,"SmbTz")==0) return SmbTzEnum; else if (strcmp(name,"SmbVz")==0) return SmbVzEnum; else if (strcmp(name,"SmbDt")==0) return SmbDtEnum; @@ -411,5 +422,4 @@ else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum; else if (strcmp(name,"SmbEC")==0) return SmbECEnum; - else if (strcmp(name,"SmbCondensation")==0) return SmbCondensationEnum; else if (strcmp(name,"SmbW")==0) return SmbWEnum; else if (strcmp(name,"SmbA")==0) return SmbAEnum; @@ -423,4 +433,6 @@ else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum; + else if (strcmp(name,"SmbDz_add")==0) return SmbDz_addEnum; + else if (strcmp(name,"SmbM_add")==0) return SmbM_addEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum; @@ -494,5 +506,8 @@ else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; - else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"Vy")==0) return VyEnum; else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; @@ -506,8 +521,5 @@ else if (strcmp(name,"VzMesh")==0) return VzMeshEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; + else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; @@ -617,5 +629,8 @@ else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; - else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; @@ -629,8 +644,5 @@ else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum; else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; + else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum; else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum; else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; @@ -740,5 +752,8 @@ else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; - else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; @@ -752,8 +767,5 @@ else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum; else if (strcmp(name,"SealevelriseMaxiter")==0) return SealevelriseMaxiterEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SealevelriseReltol")==0) return SealevelriseReltolEnum; + else if (strcmp(name,"SealevelriseReltol")==0) return SealevelriseReltolEnum; else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum; else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum; @@ -863,5 +875,8 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; @@ -875,8 +890,5 @@ else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum; else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; + else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum; else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum; else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum; Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 21215) +++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 21216) @@ -577,26 +577,29 @@ void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){ /*{{{*/ - - IssmDouble* cell=NULL; - IssmDouble* dummy=NULL; - - /*recover pointer: */ - cell=*pcell; - - /*reallocate:*/ - dummy=xNew(m+1); - + + IssmDouble* cell=NULL; + IssmDouble* dummy=NULL; + + /*recover pointer: */ + cell=*pcell; + + /*reallocate:*/ + dummy=xNew(m+1); + /*copy data:*/ - for(int i=0;i(cell); cell=dummy; - - /*assign output pointer:*/ - *pcell=cell; + if(top){ + dummy[0]=newvalue; + for(int i=0;i(cell); cell=dummy; + + /*assign output pointer:*/ + *pcell=cell; } /*}}}*/ IssmDouble cellsum(IssmDouble* cell, int m){ /*{{{*/ Index: /issm/trunk-jpl/src/m/classes/SMBgemb.m =================================================================== --- /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 21215) +++ /issm/trunk-jpl/src/m/classes/SMBgemb.m (revision 21216) @@ -22,4 +22,6 @@ isdensification; isturbulentflux; + isInitialized = NaN; + isrestart = NaN; %inputs: @@ -36,4 +38,16 @@ Tz = NaN; %height above ground at which temperature (T) was sampled [m] Vz = NaN; %height above ground at which wind (V) eas sampled [m] + + % Initialization of snow properties w restart + Dzrst = NaN; %cell depth (m) + Drst = NaN; %snow density (kg m-3) + Rerst = NaN; %effective grain size (mm) + Gdnrst = NaN; %grain dendricity (0-1) + Gsprst = NaN; %grain sphericity (0-1) + ECrst = NaN; %evaporation/condensation (kg m-2) + Wrst = NaN; %Water content (kg m-2) + Arst = NaN; %albedo (0-1) + Trst = NaN; %snow temperature (K) + Sizerst = NaN; %Number of layers %settings: @@ -117,4 +131,6 @@ self.isdensification=1; self.isturbulentflux=1; + self.isInitialized = 1*ones(mesh.numberofelements,1); + self.isrestart = 0*ones(mesh.numberofelements,1); self.aIdx = 1; @@ -126,7 +142,6 @@ self.InitDensityScaling = 1.0; - he=sum(geometry.thickness(mesh.elements),2)/size(mesh.elements,2); - self.zMax=min(500,he); - self.zMin=min(30,he); + self.zMax=250*ones(mesh.numberofelements,1); + self.zMin=130*ones(mesh.numberofelements,1); self.zY = 1.10*ones(mesh.numberofelements,1); self.outputFreq = 30; @@ -139,4 +154,15 @@ self.t0dry = 30; self.K = 7; + + self.Dzrst=0.05*ones(mesh.numberofelements,1); + self.Drst=910.0*ones(mesh.numberofelements,1); + self.Rerst=2.5*ones(mesh.numberofelements,1); + self.Gdnrst=0.0*ones(mesh.numberofelements,1); + self.Gsprst=0.0*ones(mesh.numberofelements,1); + self.ECrst=0.0*ones(mesh.numberofelements,1); + self.Wrst=0.0*ones(mesh.numberofelements,1); + self.Arst=self.aSnow*ones(mesh.numberofelements,1); + self.Trst=273.15*ones(mesh.numberofelements,1); %*CL* Default value must be Tmean + self.Sizerst=200*ones(mesh.numberofelements,1); end % }}} @@ -152,6 +178,8 @@ md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]); md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); - - md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value + md = checkfield(md,'fieldname','smb.isInitialized','values',[0 1]); + md = checkfield(md,'fieldname','smb.isrestart','values',[0 1]); + + md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400); @@ -160,5 +188,5 @@ md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1); - md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value + md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0); md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000); @@ -194,4 +222,9 @@ end md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); + + % Check that isinitialization and isrestart have different values + if any(self.isInitialized==self.isrestart), + error('isInitialized and isrestart have the same value'); + end end % }}} @@ -232,4 +265,18 @@ '3: density and cloud amount [Greuell & Konzelmann, 1994]',... '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}); + + %snow properties init + fielddisplay(self,'Dzrst','Initial cel depth when restart [m]'); + fielddisplay(self,'Drst','Initial snow density when restart [kg m-3]'); + fielddisplay(self,'Rerst','Initial grain size when restart [mm]'); + fielddisplay(self,'Gdnrst','Initial grain dendricity when restart [-]'); + fielddisplay(self,'Gsprst','Initial grain sphericity when restart [-]'); + fielddisplay(self,'ECrst','Initial evaporation/condensation when restart [kg m-2]'); + fielddisplay(self,'Wrst','Initial snow water content when restart [kg m-2]'); + fielddisplay(self,'Arst','Initial albedo when restart [-]'); + fielddisplay(self,'Trst','Initial snow temperature when restart [K]'); + fielddisplay(self,'Sizerst','Initial number of layers when restart [K]'); + + %additional albedo parameters: switch self.aIdx @@ -279,6 +326,6 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); WriteData(fid,prefix,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); - WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); - + WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2); @@ -296,4 +343,7 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double'); + + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isInitialized','format','DoubleMat','mattype',2); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','isrestart','format','DoubleMat','mattype',2); WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double'); @@ -304,4 +354,16 @@ WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double'); WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double'); + + %snow properties init + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzrst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Drst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Rerst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnrst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gsprst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECrst','format','DoubleMat','mattype',2); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wrst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Arst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Trst','format','DoubleMat','mattype',3); + WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizerst','format','IntMat','mattype',2); %figure out dt from forcings: @@ -309,5 +371,11 @@ dtime=diff(time,1); dt=min(dtime); + WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts); + + % Check if smb_dt goes evenly into transient core time step + if (mod(md.timestepping.time_step,dt) >= 1e-10) + error('smb_dt/dt = %f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt); + end %process requested outputs Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 21215) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 21216) @@ -190,5 +190,5 @@ elseif strcmp(fieldname,'SmbRunoff'), field = field*yts; - elseif strcmp(fieldname,'SmbCondensation'), + elseif strcmp(fieldname,'SmbEC'), field = field*yts; elseif strcmp(fieldname,'SmbAccumulation'), @@ -196,4 +196,8 @@ elseif strcmp(fieldname,'SmbMelt'), field = field*yts; + elseif strcmp(fieldname,'SmbDz_add'), + field = field*yts; + elseif strcmp(fieldname,'SmbM_add'), + field = field*yts; elseif strcmp(fieldname,'CalvingCalvingrate'), field = field*yts; Index: /issm/trunk-jpl/src/py3/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/py3/enum/EnumDefinitions.py (revision 21215) +++ /issm/trunk-jpl/src/py3/enum/EnumDefinitions.py (revision 21216) @@ -445,4 +445,5 @@ def SmbEvaporationEnum(): return StringToEnum("SmbEvaporation")[0] def SmbRunoffEnum(): return StringToEnum("SmbRunoff")[0] +def SmbDz_addEnum(): return StringToEnum("SmbDz_add")[0] def SMBmeltcomponentsEnum(): return StringToEnum("SMBmeltcomponents")[0] def SmbMeltEnum(): return StringToEnum("SmbMelt")[0] Index: /issm/trunk-jpl/test/NightlyRun/test243.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test243.m (revision 21215) +++ /issm/trunk-jpl/test/NightlyRun/test243.m (revision 21216) @@ -35,7 +35,7 @@ %time stepping: -md.timestepping.start_time=1979; -md.timestepping.final_time=1980; -md.timestepping.time_step=.5; +md.timestepping.start_time=1965; +md.timestepping.final_time=1968; +md.timestepping.time_step=1/365.0; md.timestepping.interp_forcings=0;