Changeset 21216


Ignore:
Timestamp:
09/21/16 17:27:00 (8 years ago)
Author:
langchar
Message:

BUG: debugging GEMB (time steps, energy balance, restart option, mergins cells)

Location:
issm/trunk-jpl
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r20690 r21216  
    5858                        iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum);
    5959                        iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum);
     60            iomodel->FetchDataToInput(elements,"md.smb.isInitialized",SmbIsInitializedEnum);
     61            iomodel->FetchDataToInput(elements,"md.smb.isrestart",SmbIsrestartEnum);
     62            iomodel->FetchDataToInput(elements,"md.smb.Dzrst",SmbDzrstEnum);
     63            iomodel->FetchDataToInput(elements,"md.smb.Drst",SmbDrstEnum);
     64            iomodel->FetchDataToInput(elements,"md.smb.Rerst",SmbRerstEnum);
     65            iomodel->FetchDataToInput(elements,"md.smb.Gdnrst",SmbGdnrstEnum);
     66            iomodel->FetchDataToInput(elements,"md.smb.Gsprst",SmbGsprstEnum);
     67            iomodel->FetchDataToInput(elements,"md.smb.ECrst",SmbECrstEnum);
     68            iomodel->FetchDataToInput(elements,"md.smb.Wrst",SmbWrstEnum);
     69            iomodel->FetchDataToInput(elements,"md.smb.Arst",SmbArstEnum);
     70            iomodel->FetchDataToInput(elements,"md.smb.Trst",SmbTrstEnum);
     71            iomodel->FetchDataToInput(elements,"md.smb.Sizerst",SmbSizerstEnum);
    6072                        break;
    6173                case SMBpddEnum:

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r21208 r21216  
    12961296/*}}}*/
    12971297void       Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/
    1298 
    1299         /*Intermediaries*/
    1300         int        i,t;
    1301         IssmDouble time;
    1302 
    1303         /*Branch on type of vector: nodal or elementary: */
    1304         if(vector_type==1){ //nodal vector
    1305 
    1306                 int         numvertices = this->GetNumberOfVertices();
    1307                 int        *vertexids   = xNew<int>(numvertices);
    1308                 IssmDouble *values      = xNew<IssmDouble>(numvertices);
    1309 
    1310                 /*Recover vertices ids needed to initialize inputs*/
    1311                 _assert_(iomodel->elements);
    1312                 for(i=0;i<numvertices;i++){
    1313                         vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    1314                 }
    1315 
    1316                 /*Are we in transient or static? */
    1317                 if(M==iomodel->numberofvertices){
    1318                         for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
    1319                         this->AddInput(vector_enum,values,P1Enum);
    1320                 }
    1321                 else if(M==iomodel->numberofvertices+1){
    1322                         /*create transient input: */
    1323                         IssmDouble* times = xNew<IssmDouble>(N);
    1324                         for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1325                         TransientInput* transientinput=new TransientInput(vector_enum,times,N);
    1326                         for(t=0;t<N;t++){
    1327                                 for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    1328                                 switch(this->ObjectEnum()){
    1329                                         case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
    1330                                         case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
    1331                                         case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
    1332                                         default: _error_("Not implemented yet");
    1333                                 }
    1334                         }
    1335                         this->inputs->AddInput(transientinput);
    1336                         xDelete<IssmDouble>(times);
    1337                 }
    1338                 else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    1339 
    1340                 xDelete<IssmDouble>(values);
    1341                 xDelete<int>(vertexids);
    1342         }
    1343         else if(vector_type==2){ //element vector
    1344 
    1345                 IssmDouble value;
    1346 
    1347                 /*Are we in transient or static? */
    1348                 if(M==iomodel->numberofelements){
    1349                         if (code==5){ //boolean
    1350                                 this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
    1351                         }
    1352                         else if (code==6){ //integer
    1353                                 this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
    1354                         }
    1355                         else if (code==7){ //IssmDouble
    1356                                 this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
    1357                         }
    1358                         else _error_("could not recognize nature of vector from code " << code);
    1359                 }
    1360                 else if(M==iomodel->numberofelements+1){
    1361                         /*create transient input: */
    1362                         IssmDouble* times = xNew<IssmDouble>(N);
    1363                         for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1364                         TransientInput* transientinput=new TransientInput(vector_enum,times,N);
    1365                         TriaInput* bof=NULL;
    1366                         for(t=0;t<N;t++){
    1367                                 value=vector[N*this->Sid()+t];
    1368                                 switch(this->ObjectEnum()){
    1369                                         case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
    1370                                         case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
    1371                                         case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
    1372                                         default: _error_("Not implemented yet");
    1373                                 }
    1374                         }
    1375                         this->inputs->AddInput(transientinput);
    1376                         xDelete<IssmDouble>(times);
    1377                 }
    1378                 else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    1379         }
    1380         else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
     1298   
     1299    /*Intermediaries*/
     1300    int        i,t;
     1301    IssmDouble time;
     1302   
     1303    /*Branch on type of vector: nodal or elementary: */
     1304    if(vector_type==1){ //nodal vector
     1305       
     1306        int         numvertices = this->GetNumberOfVertices();
     1307        int        *vertexids   = xNew<int>(numvertices);
     1308        IssmDouble *values      = xNew<IssmDouble>(numvertices);
     1309       
     1310        /*Recover vertices ids needed to initialize inputs*/
     1311        _assert_(iomodel->elements);
     1312        for(i=0;i<numvertices;i++){
     1313            vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1314        }
     1315       
     1316        /*Are we in transient or static? */
     1317        if(M==iomodel->numberofvertices){
     1318            for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
     1319            this->AddInput(vector_enum,values,P1Enum);
     1320        }
     1321        else if(M==iomodel->numberofvertices+1){
     1322            /*create transient input: */
     1323            IssmDouble* times = xNew<IssmDouble>(N);
     1324            for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1325            TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1326            for(t=0;t<N;t++){
     1327                for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1328                switch(this->ObjectEnum()){
     1329                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
     1330                    case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
     1331                    case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
     1332                    default: _error_("Not implemented yet");
     1333                }
     1334            }
     1335            this->inputs->AddInput(transientinput);
     1336            xDelete<IssmDouble>(times);
     1337        }
     1338        else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1339       
     1340        xDelete<IssmDouble>(values);
     1341        xDelete<int>(vertexids);
     1342    }
     1343    else if(vector_type==2){ //element vector
     1344       
     1345        IssmDouble value;
     1346       
     1347        /*Are we in transient or static? */
     1348        if(M==iomodel->numberofelements){
     1349            if (code==5){ //boolean
     1350                this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
     1351            }
     1352            else if (code==6){ //integer
     1353                this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
     1354            }
     1355            else if (code==7){ //IssmDouble
     1356                this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
     1357            }
     1358            else _error_("could not recognize nature of vector from code " << code);
     1359        }
     1360        else if(M==iomodel->numberofelements+1){
     1361            /*create transient input: */
     1362            IssmDouble* times = xNew<IssmDouble>(N);
     1363            for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1364            TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1365            TriaInput* bof=NULL;
     1366            for(t=0;t<N;t++){
     1367                value=vector[N*this->Sid()+t];
     1368                switch(this->ObjectEnum()){
     1369                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
     1370                    case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
     1371                    case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
     1372                    default: _error_("Not implemented yet");
     1373                }
     1374            }
     1375            this->inputs->AddInput(transientinput);
     1376            xDelete<IssmDouble>(times);
     1377        }
     1378        else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1379    }
     1380    else if(vector_type==3){ //element vector
     1381       
     1382        IssmDouble value;
     1383       
     1384        /*For right now we are static */
     1385        if(M==iomodel->numberofelements){
     1386            /*create transient input: */
     1387            IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
     1388            for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
     1389            DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N);
     1390            this->inputs->AddInput(arrayinput);
     1391            xDelete<IssmDouble>(layers);
     1392        }
     1393        else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1394    }
     1395    else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    13811396}
    13821397/*}}}*/
    … …  
    22312246
    22322247        /*Intermediary variables: {{{*/
    2233         bool       isinitialized=false;
    2234         IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
     2248    IssmDouble isinitialized;
     2249    IssmDouble isrestart;
     2250    IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
    22352251        IssmDouble Tmean;
    22362252        IssmDouble C;
    … …  
    22512267        IssmDouble initMass;
    22522268    IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
     2269    IssmDouble sumdz_add;
    22532270    IssmDouble sumMass,dMass;
    22542271        bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
    2255         IssmDouble init_scaling=1.0;
     2272        IssmDouble init_scaling;
    22562273        /*}}}*/
    22572274        /*Output variables:{{{ */
    … …  
    22702287        IssmDouble  R;
    22712288        IssmDouble  mAdd;
     2289    IssmDouble  dz_add;
     2290   
     2291    IssmDouble* dzrst=NULL;
     2292    IssmDouble* drst = NULL;
     2293    IssmDouble* rerst = NULL;
     2294    IssmDouble* gdnrst = NULL;
     2295    IssmDouble* gsprst = NULL;
     2296    IssmDouble* Wrst = NULL;
     2297    IssmDouble* arst = NULL;
     2298    IssmDouble* Trst = NULL;
     2299
    22722300        int         m;
    22732301        int         count=0;
    … …  
    22852313        parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
    22862314        parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
     2315    parameters->FindParam(&yts,ConstantsYtsEnum);
    22872316        parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
    22882317        parameters->FindParam(&aIdx,SmbAIdxEnum);
    … …  
    23212350        Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input);
    23222351        Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input);
    2323         Input* isinitialized_input=this->inputs->GetInput(SmbIsInitializedEnum);
    2324        
     2352    Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum);  _assert_(isinitialized_input);
     2353    Input* isrestart_input=this->GetInput(SmbIsrestartEnum);  _assert_(isrestart_input);
     2354   
    23252355        /*Retrieve input values:*/
    23262356        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
    … …  
    23362366        Tz_input->GetInputValue(&Tz,gauss);
    23372367        Vz_input->GetInputValue(&Vz,gauss);
     2368    isinitialized_input->GetInputValue(&isinitialized,gauss);
     2369    isrestart_input->GetInputValue(&isrestart,gauss);
    23382370        /*}}}*/
    23392371
    2340         /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
    2341         if(!isinitialized_input){
    2342 
    2343                 if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
    2344                 GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
    2345                 //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
    2346                 //dz[i] << "\n");
    2347 
    2348                 /*initialize profile variables:*/
    2349                 d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor
    2350                 re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5;         //set grain size to old snow [mm]
    2351                 gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0;         //set grain dentricity to old snow
    2352                 gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0;         //set grain sphericity to old snow
    2353                 EC = 0;                                                                //surface evaporation (-) condensation (+) [kg m-2]
    2354                 W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0;             //set water content to zero [kg m-2]
    2355                 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //set albedo equal to fresh snow [fraction]
    2356                 T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
    2357 
    2358                 //fixed lower temperatuer bounday condition - T is fixed
    2359                 T_bottom=T[m-1];
    2360 
    2361                 /*Flag the initialization:*/
    2362                 this->AddInput(new BoolInput(SmbIsInitializedEnum,true));
    2363         }
    2364         else{
     2372    /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
     2373    if(isinitialized){
     2374        if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
     2375        if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");
     2376        GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
     2377        //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
     2378        //dz[i] << "\n");
     2379        /*initialize profile variables:*/
     2380        d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor
     2381        re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5;         //set grain size to old snow [mm]
     2382        gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0;         //set grain dentricity to old snow
     2383        gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0;         //set grain sphericity to old snow
     2384        EC = 0;                                                                //surface evaporation (-) condensation (+) [kg m-2]
     2385        W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0;             //set water content to zero [kg m-2]
     2386        a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //set albedo equal to fresh snow [fraction]
     2387        T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
     2388       
     2389        //fixed lower temperature bounday condition - T is fixed
     2390        T_bottom=T[m-1];
     2391       
     2392        /*Flag the initialization:*/
     2393        this->AddInput(new DoubleInput(SmbIsInitializedEnum,false));
     2394    }
     2395    else if(isrestart){ //Retrieve the snow properties from previous run
     2396        if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
     2397       
     2398        /*Recover inputs: */
     2399        DoubleArrayInput* dzrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzrstEnum)); _assert_(dzrst_input);
     2400        DoubleArrayInput* drst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDrstEnum));_assert_(drst_input);
     2401        DoubleArrayInput* rerst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbRerstEnum));_assert_(rerst_input);
     2402        DoubleArrayInput* gdnrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnrstEnum));_assert_(gdnrst_input);
     2403        DoubleArrayInput* gsprst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGsprstEnum));_assert_(gsprst_input);
     2404        DoubleInput* ECrst_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECrstEnum));_assert_(ECrst_input);
     2405        DoubleArrayInput* Wrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWrstEnum));_assert_(Wrst_input);
     2406        DoubleArrayInput* arst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbArstEnum));_assert_(arst_input);
     2407        DoubleArrayInput* Trst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTrstEnum));_assert_(Trst_input);
     2408       
     2409        /*Recover arrays: */
     2410        dzrst_input->GetValues(&dzrst,&m);
     2411        drst_input->GetValues(&drst,&m);
     2412        rerst_input->GetValues(&rerst,&m);
     2413        gdnrst_input->GetValues(&gdnrst,&m);
     2414        gsprst_input->GetValues(&gsprst,&m);
     2415        ECrst_input->GetInputValue(&EC);
     2416        Wrst_input->GetValues(&Wrst,&m);
     2417        arst_input->GetValues(&arst,&m);
     2418        Trst_input->GetValues(&Trst,&m);
     2419       
     2420        /*Retrive the correct value of m (without the zeroes at the end)*/
     2421        Input* Sizerst_input=this->GetInput(SmbSizerstEnum); _assert_(Sizerst_input);
     2422        Sizerst_input->GetInputValue(&m);
     2423       
     2424        dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzrst[i];
     2425        d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=drst[i];
     2426        re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=rerst[i];
     2427        gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnrst[i];
     2428        gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gsprst[i];
     2429        W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wrst[i];
     2430        a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=arst[i];
     2431        T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Trst[i];
     2432       
     2433        //fixed lower temperatuer bounday condition - T is fixed
     2434        T_bottom=T[m-1];
     2435       
     2436        /*Flag the initialization:*/
     2437        this->AddInput(new DoubleInput(SmbIsrestartEnum,false));//*CL* add
     2438    }
     2439    else{
     2440 
    23652441                /*Recover inputs: */
    23662442                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
    … …  
    23952471    // initialize cumulative variables
    23962472    sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
     2473    sumdz_add=0;
    23972474
    23982475        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
    … …  
    24022479        /*Start loop: */
    24032480        count=1;
    2404         for (t=time;t<=time+dt;t=t+smb_dt){
     2481        for (t=time;t<time+dt;t=t+smb_dt){
    24052482
    24062483                if(VerboseSmb() && this->Sid()==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");
    … …  
    24422519                /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    24432520                 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    2444                 if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin,this->Sid());
     2521        if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,this->Sid());
    24452522
    24462523                /*Allow non-melt densification and determine compaction [m]*/
    … …  
    24792556        sumP = P +  sumP;
    24802557        sumEC = sumEC + EC;  // evap (-)/cond(+)
     2558        sumdz_add=dz_add+sumdz_add; //*CL*
    24812559
    24822560                /*Calculate total system mass:*/
    … …  
    24992577                /*increase counter:*/
    25002578                count++;
    2501 
    2502         } //for (t=time;t<=time+dt;t=t+smb_dt)
    2503 
    2504 
     2579        } //for (t=time;t<time+dt;t=t+smb_dt)
     2580   
    25052581        /*Save generated inputs: */
    25062582        this->AddInput(new DoubleArrayInput(SmbDzEnum,dz,m));
    … …  
    25102586        this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m));
    25112587        this->AddInput(new DoubleArrayInput(SmbTEnum,T,m));
    2512         this->AddInput(new DoubleInput(SmbECEnum,EC));
     2588        this->AddInput(new DoubleInput(SmbECEnum,sumEC/yts));
    25132589        this->AddInput(new DoubleArrayInput(SmbWEnum,W,m));
    25142590        this->AddInput(new DoubleArrayInput(SmbAEnum,a,m));
    2515         this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/rho_water/dt));
    2516         this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/rho_water/dt));
    2517         this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/rho_water/dt));
    2518         this->AddInput(new DoubleInput(SmbCondensationEnum,sumEC/rho_water/dt));
    2519 
    2520 
     2591    this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/yts));
     2592    this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/yts));
     2593    this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/yts));
     2594    this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts));
     2595    this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts));
    25212596
    25222597        /*Free allocations:{{{*/

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

    r19613 r21216  
    11831183        *pm=m;
    11841184} /*}}}*/
    1185 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){ /*{{{*/
     1185void 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){ /*{{{*/
    11861186
    11871187        //// MELT ROUTINE
    … …  
    12181218        IssmDouble X;
    12191219        IssmDouble Wi;
     1220   
     1221    IssmDouble Ztot;
     1222    IssmDouble T_bot;
     1223    IssmDouble m_bot;
     1224    IssmDouble dz_bot;
     1225    IssmDouble d_bot;
     1226    IssmDouble W_bot;
     1227    IssmDouble a_bot;
     1228    IssmDouble re_bot;
     1229    IssmDouble gdn_bot;
     1230    IssmDouble gsp_bot;
     1231    bool        top=false;
     1232   
     1233    IssmDouble* Zcum=NULL;
     1234    IssmDouble* dzMin2=NULL;
     1235    IssmDouble zY2=1.025;
     1236    bool lastCellFlag;
     1237    int X1=0;
     1238    int X2=0;
     1239   
    12201240        int        D_size;
    12211241        int         i;
    … …  
    12231243        /*outputs:*/
    12241244        IssmDouble  mAdd;
     1245    IssmDouble dz_add;
    12251246        IssmDouble  Rsum;
    12261247        IssmDouble* T=*pT;
    … …  
    12641285        mAdd = 0;       // mass added/removed to/from base of model [kg]
    12651286        addE = 0;       // energy added/removed to/from base of model [J]
     1287    dz_add=0;      // thickness of the layer added/removed to/from base of model [m]
    12661288
    12671289        // calculate temperature excess above 0 deg C
    … …  
    13011323        }
    13021324
    1303         //// MELT, PERCOLATION AND REFREEZE
    1304 
    1305         // run melt algorithm if there is melt water or excess pore water
    1306         if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
    1307                
    1308                 // _printf_(""MELT OCCURS");
    1309                 // check to see if thermal energy exceeds energy to melt entire cell
    1310                 // if so redistribute temperature to lower cells (temperature surplus)
    1311                 // (maximum T of snow before entire grid cell melts is a constant
    1312                 // LF/CI = 159.1342)
    1313                 surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0, exsT [i]- 159.1342);
    1314 
    1315                 if (cellsum(surpT,n) > 0 ){
    1316                         // _printf_("T Surplus");
    1317                         // calculate surplus energy
    1318                         surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI / m[i];
    1319                        
    1320                         int i = 0;
    1321                         while (cellsum(surpE,n) > 0){
    1322                                 // use surplus energy to increase the temperature of lower cell
    1323                                 T[i+1] = surpE[i] * m[i+1]/CI + T[i+1];
    1324                                 surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342));
    1325                                 surpE[i+1] = surpT[i+1] * CI / m[i+1];
    1326 
    1327                                 // adjust current cell properties (again 159.1342 is the max T)
    1328                                 T[i] = CtoK + 159.1342;
    1329                                 surpE[i] = 0;   
    1330                                 i = i + 1;
    1331                         }
    1332                         // recalculate temperature excess above 0 deg C
    1333                         for(int i=0;i<n;i++) exsT[i] = fmax(0, T[i] - CtoK);
    1334                 }
     1325    //// MELT, PERCOLATION AND REFREEZE
     1326   
     1327    // run melt algorithm if there is melt water or excess pore water
     1328    if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){
     1329        // _printf_(""MELT OCCURS");
     1330        // check to see if thermal energy exceeds energy to melt entire cell
     1331        // if so redistribute temperature to lower cells (temperature surplus)
     1332        // (maximum T of snow before entire grid cell melts is a constant
     1333        // LF/CI = 159.1342)
     1334        surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0, exsT [i]- 159.1342);
     1335       
     1336        if (cellsum(surpT,n) > 0 ){
     1337            // _printf_("T Surplus");
     1338            // calculate surplus energy
     1339            surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];//*CL*
     1340           
     1341            int i = 0;
     1342            while (cellsum(surpE,n) > 0){
     1343                // use surplus energy to increase the temperature of lower cell
     1344                T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];//*CL*
     1345               
     1346                exsT[i+1] = fmax(0, T[i+1] - CtoK) + exsT[i+1];
     1347                T[i+1] = fmin(CtoK, T[i+1]);
     1348               
     1349                surpT[i+1] = fmax(0, exsT[i+1] - 159.1342);
     1350                surpE[i+1] = surpT[i+1] * CI * m[i+1];//*CL*
     1351               
     1352                // adjust current cell properties (again 159.1342 is the max T)
     1353                exsT[i] = 159.1342;
     1354                surpE[i] = 0;   
     1355                i = i + 1;
     1356            }
     1357        }
    13351358
    13361359                // convert temperature excess to melt [kg]
    … …  
    14691492                xDelete<IssmDouble>(R);
    14701493        }
     1494   
     1495    //Merging of cells as they are burried under snow.
     1496    Zcum=xNew<IssmDouble>(n);
     1497    dzMin2=xNew<IssmDouble>(n);
     1498   
     1499    Zcum[0]=dz[0]; // Compute a cumulative depth vector
     1500   
     1501    for (int i=1;i<n;i++){
     1502        Zcum[i]=Zcum[i-1]+dz[i];
     1503    }
     1504   
     1505    for (int i=0;i<n;i++){
     1506        if (Zcum[i]<=zTop){
     1507            dzMin2[i]=dzMin;
     1508        }
     1509        else{
     1510            dzMin2[i]=zY2*dzMin2[i-1];
     1511        }
     1512    }
    14711513
    14721514        // check if depth is too small
    14731515        X = 0;
    14741516        for(int i=n-1;i>=0;i--){
    1475                 if(dz[i]<dzMin){
     1517         if(dz[i]<dzMin2[i]){
    14761518                        X=i;
    14771519                        break;
    14781520                }
    14791521        }
    1480 
    1481         for (int i = 0; i<=X;i++){
    1482                 if (dz [i] < dzMin){                               // merge top cells                                                                 
    1483                         //                                                                          _printf_("dz > dzMin * 2')
    1484                         // adjust variables as a linearly weighted function of mass
    1485                         IssmDouble m_new = m[i] + m[i+1];
    1486                         T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
    1487                         a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
    1488                         re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
    1489                         gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
    1490                         gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
    1491 
    1492                         // merge with underlying grid cell and delete old cell
    1493                         dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
    1494                         d[i+1] = m_new / dz[i+1];                   // combine top densities
    1495                         W[i+1] = W[i+1] + W[i];                     // combine liquid water
    1496                         m[i+1] = m_new;                             // combine top masses
    1497 
    1498                         // set cell to 99999 for deletion
    1499                         m[i] = 99999;
    1500                 }
    1501         }
     1522   
     1523    //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
     1524    if(X==n-1){
     1525        lastCellFlag = true;
     1526        X = X-1;
     1527    }
     1528    else{
     1529        lastCellFlag = false;
     1530    }
     1531
     1532    for (int i = 0; i<=X;i++){
     1533        if (dz [i] < dzMin2[i]){                               // merge top cells
     1534            //                                                                          _printf_("dz > dzMin * 2')
     1535            // adjust variables as a linearly weighted function of mass
     1536            IssmDouble m_new = m[i] + m[i+1];
     1537            T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
     1538            a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
     1539            re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
     1540            gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
     1541            gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
     1542           
     1543            // merge with underlying grid cell and delete old cell
     1544            dz [i+1] = dz[i] + dz[i+1];                 // combine cell depths
     1545            d[i+1] = m_new / dz[i+1];                   // combine top densities
     1546            W[i+1] = W[i+1] + W[i];                     // combine liquid water
     1547            m[i+1] = m_new;                             // combine top masses
     1548           
     1549            // set cell to 99999 for deletion
     1550            m[i] = 99999;
     1551        }
     1552    }
     1553
     1554    //If last cell has to be merged
     1555    if(lastCellFlag){
     1556        //find closest cell to merge with
     1557        for(int i=n-2;i>=0;i--){
     1558            if(m[i]!=99999){
     1559                X2=i;
     1560                X1=n-1;
     1561                break;
     1562            }
     1563        }
     1564       
     1565        // adjust variables as a linearly weighted function of mass
     1566        IssmDouble m_new = m[X2] + m[X1];
     1567        T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
     1568        a[X1] = (a[X2]*m[i] + a[X1]*m[X1]) / m_new;
     1569        re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
     1570        gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
     1571        gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
     1572       
     1573        // merge with underlying grid cell and delete old cell
     1574        dz [X1] = dz[X2] + dz[X1];                 // combine cell depths
     1575        d[X1] = m_new / dz[X1];                   // combine top densities
     1576        W[X1] = W[X1] + W[X2];                     // combine liquid water
     1577        m[X1] = m_new;                             // combine top masses
     1578       
     1579        // set cell to 99999 for deletion
     1580        m[X2] = 99999;
     1581    }
    15021582
    15031583        // delete combined cells
    … …  
    15511631        }
    15521632
    1553         //// CORRECT FOR TOTAL MODEL DEPTH
    1554         // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
    1555         // INTERPRETATION
    1556 
    1557         // // calculate total model depth
    1558         // z = sum(dz);
    1559         //
    1560         // if (z < zMin){ // check if model is too shallow                                                                     
    1561         //                                                                          _printf_("z < zMin')
    1562         //     // mass and energy to be added
    1563         //     mAdd = m(end) + W(end);
    1564         //     addE = T(end) * m(end) * CI;
    1565         //     
    1566         //     // add a grid cell of the same size and temperature to the bottom
    1567         //     dz = [dz; dz(end)];
    1568         //     T = [T; T(end)];
    1569         //     W = [W; W(end)];
    1570         //     m = [m; m(end)];
    1571         //     d = [d; d(end)];
    1572         //     a = [a; a(end)];
    1573         //     re = [re; re(end)];
    1574         //     gdn = [gdn; gdn(end)];
    1575         //     gsp = [gsp; gsp(end)];
    1576         // }
    1577         // else (if z > zMax){ // check if model is too deep                                                                                                                                       
    1578         //                                                                          _printf_("z > zMax')
    1579         //     // mass and energy loss
    1580         //     mAdd = -(m(end) + W(end));
    1581         //     addE = -(T(end) * m(end) * CI);
    1582         //     
    1583         //     // add a grid cell of the same size and temperature to the bottom
    1584         //     dz(end) = []; T(end) = []; W(end) = []; m(end) = [];
    1585         //     d(end) = []; a(end) = [];  re(end) = [];  gdn(end) = []; 
    1586         //     gsp(end) = [];
    1587         // }
     1633    //// CORRECT FOR TOTAL MODEL DEPTH
     1634    // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT
     1635    // INTERPRETATION
     1636    // // calculate total model depth
     1637    Ztot = cellsum(dz,n);
     1638   
     1639    if (Ztot < zMin){
     1640//        printf("Total depth < zMin %f \n", Ztot);
     1641        // mass and energy to be added
     1642        mAdd = m[n-1]+W[n-1];
     1643        addE = T[n-1]*m[n-1]*CI;
     1644       
     1645        // add a grid cell of the same size and temperature to the bottom
     1646        dz_bot=dz[n-1];
     1647        T_bot=T[n-1];
     1648        W_bot=W[n-1];
     1649        m_bot=m[n-1];
     1650        d_bot=d[n-1];
     1651        a_bot=a[n-1];
     1652        re_bot=re[n-1];
     1653        gdn_bot=gdn[n-1];
     1654        gsp_bot=gsp[n-1];
     1655       
     1656        dz_add=dz_bot;
     1657       
     1658        newcell(&dz,dz_bot,top,n);
     1659        newcell(&T,T_bot,top,n);
     1660        newcell(&W,W_bot,top,n);
     1661        newcell(&m,m_bot,top,n);
     1662        newcell(&d,d_bot,top,n);
     1663        newcell(&a,a_bot,top,n);
     1664        newcell(&re,re_bot,top,n);
     1665        newcell(&gdn,gdn_bot,top,n);
     1666        newcell(&gsp,gsp_bot,top,n);
     1667        n=n+1;
     1668    }
     1669    else if (Ztot > zMax){
     1670//        printf("Total depth > zMax %f \n", Ztot);
     1671        // mass and energy loss
     1672        mAdd = -(m[n-1]+W[n-1]);
     1673        addE = -(T[n-1]*m[n-1]*CI);
     1674        dz_add=-(dz[n-1]);
     1675       
     1676        // add a grid cell of the same size and temperature to the bottom
     1677        D_size=n-1;
     1678        D=xNew<int>(D_size);
     1679       
     1680        for(int i=0;i<n-1;i++) D[i]=i;
     1681        celldelete(&dz,n,D,D_size);
     1682        celldelete(&T,n,D,D_size);
     1683        celldelete(&W,n,D,D_size);
     1684        celldelete(&m,n,D,D_size);
     1685        celldelete(&d,n,D,D_size);
     1686        celldelete(&a,n,D,D_size);
     1687        celldelete(&re,n,D,D_size);
     1688        celldelete(&gdn,n,D,D_size);
     1689        celldelete(&gsp,n,D,D_size);
     1690        n=D_size;
     1691        xDelete<int>(D);
     1692       
     1693    }
    15881694
    15891695        //// CHECK FOR MASS AND ENERGY CONSERVATION
    … …  
    16261732        *pR=Rsum;
    16271733        *pmAdd=mAdd;
    1628        
     1734    *pdz_add=dz_add;
     1735
    16291736        *pT=T;
    16301737        *pd=d;

  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r19582 r21216  
    2828void 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);
    2929void 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);
    30 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);
     30void 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);
    3131void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,IssmDouble rho_ice,int m, int sid);
    3232void 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);

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r21208 r21216  
    361361        SmbRequestedOutputsEnum,
    362362        SmbIsInitializedEnum,
     363    SmbIsrestartEnum,
     364    SmbDzrstEnum,
     365    SmbDrstEnum,
     366    SmbRerstEnum,
     367    SmbGdnrstEnum,
     368    SmbGsprstEnum,
     369    SmbECrstEnum,
     370    SmbWrstEnum,
     371    SmbArstEnum,
     372    SmbTrstEnum,
     373    SmbSizerstEnum,
    363374        /*SMBforcing*/
    364375        SMBforcingEnum,
    … …  
    401412        SmbGspEnum,
    402413        SmbECEnum,
    403         SmbCondensationEnum,
    404414        SmbWEnum,
    405415        SmbAEnum,
    … …  
    413423        SmbIsdensificationEnum,
    414424        SmbIsturbulentfluxEnum,
     425    SmbDz_addEnum,
     426    SmbM_addEnum,
    415427        /*SMBpdd*/
    416428        SMBpddEnum,     

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r21103 r21216  
    364364                case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
    365365                case SmbIsInitializedEnum : return "SmbIsInitialized";
     366                case SmbIsrestartEnum : return "SmbIsrestart";
     367                case SmbDzrstEnum : return "SmbDzrst";
     368                case SmbDrstEnum : return "SmbDrst";
     369                case SmbRerstEnum : return "SmbRerst";
     370                case SmbGdnrstEnum : return "SmbGdnrst";
     371                case SmbGsprstEnum : return "SmbGsprst";
     372                case SmbECrstEnum : return "SmbECrst";
     373                case SmbWrstEnum : return "SmbWrst";
     374                case SmbArstEnum : return "SmbArst";
     375                case SmbTrstEnum : return "SmbTrst";
     376                case SmbSizerstEnum : return "SmbSizerst";
    366377                case SMBforcingEnum : return "SMBforcing";
    367378                case SmbMassBalanceEnum : return "SmbMassBalance";
    … …  
    402413                case SmbGspEnum : return "SmbGsp";
    403414                case SmbECEnum : return "SmbEC";
    404                 case SmbCondensationEnum : return "SmbCondensation";
    405415                case SmbWEnum : return "SmbW";
    406416                case SmbAEnum : return "SmbA";
    … …  
    414424                case SmbIsdensificationEnum : return "SmbIsdensification";
    415425                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
     426                case SmbDz_addEnum : return "SmbDz_add";
     427                case SmbM_addEnum : return "SmbM_add";
    416428                case SMBpddEnum : return "SMBpdd";
    417429                case SmbDelta18oEnum : return "SmbDelta18o";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r21103 r21216  
    370370              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    371371              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
     372              else if (strcmp(name,"SmbIsrestart")==0) return SmbIsrestartEnum;
     373              else if (strcmp(name,"SmbDzrst")==0) return SmbDzrstEnum;
     374              else if (strcmp(name,"SmbDrst")==0) return SmbDrstEnum;
     375              else if (strcmp(name,"SmbRerst")==0) return SmbRerstEnum;
     376              else if (strcmp(name,"SmbGdnrst")==0) return SmbGdnrstEnum;
     377              else if (strcmp(name,"SmbGsprst")==0) return SmbGsprstEnum;
     378              else if (strcmp(name,"SmbECrst")==0) return SmbECrstEnum;
     379              else if (strcmp(name,"SmbWrst")==0) return SmbWrstEnum;
     380              else if (strcmp(name,"SmbArst")==0) return SmbArstEnum;
     381              else if (strcmp(name,"SmbTrst")==0) return SmbTrstEnum;
     382              else if (strcmp(name,"SmbSizerst")==0) return SmbSizerstEnum;
    372383              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    373384              else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    374               else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
    375389              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    376390              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    … …  
    383397              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
    384398              else if (strcmp(name,"SmbC")==0) return SmbCEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
     399              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    389400              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    390401              else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
    … …  
    411422              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    412423              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
    413               else if (strcmp(name,"SmbCondensation")==0) return SmbCondensationEnum;
    414424              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
    415425              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    … …  
    423433              else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
    424434              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
     435              else if (strcmp(name,"SmbDz_add")==0) return SmbDz_addEnum;
     436              else if (strcmp(name,"SmbM_add")==0) return SmbM_addEnum;
    425437              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    426438              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    … …  
    494506              else if (strcmp(name,"Vx")==0) return VxEnum;
    495507              else if (strcmp(name,"VxPicard")==0) return VxPicardEnum;
    496               else if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"VyAverage")==0) return VyAverageEnum;
    497512              else if (strcmp(name,"Vy")==0) return VyEnum;
    498513              else if (strcmp(name,"VyPicard")==0) return VyPicardEnum;
    … …  
    506521              else if (strcmp(name,"VzMesh")==0) return VzMeshEnum;
    507522              else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
     523              else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
    512524              else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
    513525              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
    … …  
    617629              else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum;
    618630              else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum;
    619               else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum;
    620635              else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
    621636              else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
    … …  
    629644              else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
    630645              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
     646              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    635647              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    636648              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    … …  
    740752              else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
    741753              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    742               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
    743758              else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
    744759              else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum;
    … …  
    752767              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    753768              else if (strcmp(name,"SealevelriseMaxiter")==0) return SealevelriseMaxiterEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"SealevelriseReltol")==0) return SealevelriseReltolEnum;
     769              else if (strcmp(name,"SealevelriseReltol")==0) return SealevelriseReltolEnum;
    758770              else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum;
    759771              else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum;
    … …  
    863875              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    864876              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    865               else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    866881              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    867882              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    … …  
    875890              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    876891              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
     892              else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
    881893              else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
    882894              else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r20732 r21216  
    577577
    578578void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){  /*{{{*/
    579 
    580         IssmDouble* cell=NULL;
    581         IssmDouble* dummy=NULL;
    582 
    583         /*recover pointer: */
    584         cell=*pcell;
    585 
    586         /*reallocate:*/
    587         dummy=xNew<IssmDouble>(m+1);
    588 
     579   
     580    IssmDouble* cell=NULL;
     581    IssmDouble* dummy=NULL;
     582   
     583    /*recover pointer: */
     584    cell=*pcell;
     585   
     586    /*reallocate:*/
     587    dummy=xNew<IssmDouble>(m+1);
     588   
    589589        /*copy data:*/
    590         for(int i=0;i<m;i++)dummy[i+1]=cell[i];
    591 
    592         /*top or bottom layer:*/
    593         if(top) dummy[0]=newvalue;
    594         else dummy[m]=newvalue;
    595        
    596         /*delete and reassign: */
    597         xDelete<IssmDouble>(cell); cell=dummy;
    598 
    599         /*assign output pointer:*/
    600         *pcell=cell;
     590    if(top){
     591        dummy[0]=newvalue;
     592        for(int i=0;i<m;i++)dummy[i+1]=cell[i];
     593    }
     594    else{
     595        dummy[m]=newvalue;
     596        for(int i=0;i<m;i++)dummy[i]=cell[i];
     597    }
     598   
     599    /*delete and reassign: */
     600    xDelete<IssmDouble>(cell); cell=dummy;
     601   
     602    /*assign output pointer:*/
     603    *pcell=cell;
    601604} /*}}}*/
    602605IssmDouble  cellsum(IssmDouble* cell, int m){ /*{{{*/

  • issm/trunk-jpl/src/m/classes/SMBgemb.m

    r20902 r21216  
    2222                isdensification;
    2323                isturbulentflux;
     24        isInitialized = NaN;
     25        isrestart = NaN;     
    2426
    2527                %inputs:
    … …  
    3638                Tz    = NaN; %height above ground at which temperature (T) was sampled [m]
    3739                Vz    = NaN; %height above ground at which wind (V) eas sampled [m]
     40       
     41        % Initialization of snow properties w restart
     42        Dzrst = NaN; %cell depth (m)
     43        Drst = NaN; %snow density (kg m-3)
     44        Rerst = NaN; %effective grain size (mm)
     45        Gdnrst = NaN; %grain dendricity (0-1)
     46        Gsprst = NaN; %grain sphericity (0-1)
     47        ECrst = NaN; %evaporation/condensation (kg m-2)
     48        Wrst = NaN; %Water content (kg m-2)
     49        Arst = NaN; %albedo (0-1)
     50        Trst = NaN; %snow temperature (K)
     51        Sizerst = NaN; %Number of layers
    3852
    3953                %settings:
    … …  
    117131                self.isdensification=1;
    118132                self.isturbulentflux=1;
     133        self.isInitialized = 1*ones(mesh.numberofelements,1);
     134        self.isrestart = 0*ones(mesh.numberofelements,1);
    119135       
    120136                self.aIdx = 1;
    … …  
    126142                self.InitDensityScaling = 1.0;
    127143               
    128                 he=sum(geometry.thickness(mesh.elements),2)/size(mesh.elements,2);
    129                 self.zMax=min(500,he);
    130                 self.zMin=min(30,he);
     144        self.zMax=250*ones(mesh.numberofelements,1);
     145                self.zMin=130*ones(mesh.numberofelements,1);
    131146                self.zY = 1.10*ones(mesh.numberofelements,1);
    132147                self.outputFreq = 30;
    … …  
    139154                self.t0dry = 30;
    140155                self.K = 7;
     156       
     157        self.Dzrst=0.05*ones(mesh.numberofelements,1);
     158        self.Drst=910.0*ones(mesh.numberofelements,1);
     159        self.Rerst=2.5*ones(mesh.numberofelements,1);
     160        self.Gdnrst=0.0*ones(mesh.numberofelements,1);
     161        self.Gsprst=0.0*ones(mesh.numberofelements,1);
     162        self.ECrst=0.0*ones(mesh.numberofelements,1);
     163        self.Wrst=0.0*ones(mesh.numberofelements,1);
     164        self.Arst=self.aSnow*ones(mesh.numberofelements,1);
     165        self.Trst=273.15*ones(mesh.numberofelements,1); %*CL* Default value must be Tmean
     166        self.Sizerst=200*ones(mesh.numberofelements,1);
    141167
    142168                end % }}}
    … …  
    152178                        md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
    153179                        md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
    154 
    155                         md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
     180            md = checkfield(md,'fieldname','smb.isInitialized','values',[0 1]);
     181            md = checkfield(md,'fieldname','smb.isrestart','values',[0 1]);
     182
     183                        md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value
    156184                        md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<',45); %max 500 km/h
    157185                        md = checkfield(md,'fieldname','smb.dswrf','timeseries',1,'NaN',1,'Inf',1,'>=',0,'<=',1400);
    … …  
    160188                        md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1,'Inf',1);
    161189
    162                         md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>',273-60,'<',273+60); %60 celsius max value
     190                        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
    163191                        md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0);
    164192                        md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'Inf',1,'>=',0,'<=',5000);
    … …  
    194222                        end
    195223                        md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1);
     224           
     225            % Check that isinitialization and isrestart have different values
     226            if any(self.isInitialized==self.isrestart),
     227                error('isInitialized and isrestart have the same value');
     228            end
    196229
    197230                end % }}}
    … …  
    232265                                                                        '3: density and cloud amount [Greuell & Konzelmann, 1994]',...
    233266                                                                        '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'});
     267                               
     268            %snow properties init
     269            fielddisplay(self,'Dzrst','Initial cel depth when restart [m]');
     270            fielddisplay(self,'Drst','Initial snow density when restart [kg m-3]');
     271            fielddisplay(self,'Rerst','Initial grain size when restart [mm]');
     272            fielddisplay(self,'Gdnrst','Initial grain dendricity when restart [-]');
     273            fielddisplay(self,'Gsprst','Initial grain sphericity when restart [-]');
     274            fielddisplay(self,'ECrst','Initial evaporation/condensation when restart [kg m-2]');
     275            fielddisplay(self,'Wrst','Initial snow water content when restart [kg m-2]');
     276            fielddisplay(self,'Arst','Initial albedo when restart [-]');
     277            fielddisplay(self,'Trst','Initial snow temperature when restart [K]');
     278            fielddisplay(self,'Sizerst','Initial number of layers when restart [K]');
     279           
     280           
    234281                        %additional albedo parameters:
    235282                        switch self.aIdx
    … …  
    279326                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    280327                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    281                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
    282                        
     328                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);         
     329
    283330                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2);
    284331                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2);
    … …  
    296343                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer');
    297344                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double');
     345           
     346            WriteData(fid,prefix,'object',self,'class','smb','fieldname','isInitialized','format','DoubleMat','mattype',2);
     347            WriteData(fid,prefix,'object',self,'class','smb','fieldname','isrestart','format','DoubleMat','mattype',2);
    298348
    299349                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double');
    … …  
    304354                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','t0dry','format','Double');
    305355                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double');
     356           
     357            %snow properties init
     358            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzrst','format','DoubleMat','mattype',3);
     359            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Drst','format','DoubleMat','mattype',3);
     360            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Rerst','format','DoubleMat','mattype',3);
     361            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnrst','format','DoubleMat','mattype',3);
     362            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gsprst','format','DoubleMat','mattype',3);
     363            WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECrst','format','DoubleMat','mattype',2);
     364            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wrst','format','DoubleMat','mattype',3);
     365            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Arst','format','DoubleMat','mattype',3);
     366            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Trst','format','DoubleMat','mattype',3);
     367            WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizerst','format','IntMat','mattype',2);
    306368
    307369                        %figure out dt from forcings:
    … …  
    309371                        dtime=diff(time,1);
    310372                        dt=min(dtime);
     373           
    311374                        WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts);
     375           
     376            % Check if smb_dt goes evenly into transient core time step
     377            if (mod(md.timestepping.time_step,dt) >= 1e-10)
     378                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);
     379            end
    312380                       
    313381                        %process requested outputs

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r20900 r21216  
    190190        elseif strcmp(fieldname,'SmbRunoff'),
    191191                field = field*yts;
    192         elseif strcmp(fieldname,'SmbCondensation'),
     192        elseif strcmp(fieldname,'SmbEC'),
    193193                field = field*yts;
    194194        elseif strcmp(fieldname,'SmbAccumulation'),
    … …  
    196196        elseif strcmp(fieldname,'SmbMelt'),
    197197                field = field*yts;
     198    elseif strcmp(fieldname,'SmbDz_add'),
     199        field = field*yts;
     200    elseif strcmp(fieldname,'SmbM_add'),
     201        field = field*yts;
    198202        elseif strcmp(fieldname,'CalvingCalvingrate'),
    199203                field = field*yts;

  • issm/trunk-jpl/src/py3/enum/EnumDefinitions.py

    r20465 r21216  
    445445def SmbEvaporationEnum(): return StringToEnum("SmbEvaporation")[0]
    446446def SmbRunoffEnum(): return StringToEnum("SmbRunoff")[0]
     447def SmbDz_addEnum(): return StringToEnum("SmbDz_add")[0]
    447448def SMBmeltcomponentsEnum(): return StringToEnum("SMBmeltcomponents")[0]
    448449def SmbMeltEnum(): return StringToEnum("SmbMelt")[0]

  • issm/trunk-jpl/test/NightlyRun/test243.m

    r21056 r21216  
    3535
    3636%time stepping:
    37 md.timestepping.start_time=1979;
    38 md.timestepping.final_time=1980;
    39 md.timestepping.time_step=.5;
     37md.timestepping.start_time=1965;
     38md.timestepping.final_time=1968;
     39md.timestepping.time_step=1/365.0;
    4040md.timestepping.interp_forcings=0;
    4141
Note: See TracChangeset for help on using the changeset viewer.