Changeset 25758

Timestamp:

11/16/20 13:22:37 (4 years ago)

Author:

jdquinn

Message:

CHG: Missing changes from merge from Eric’s branch

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (7 diffs)
• c/classes/Elements/Tria.cpp (modified) (13 diffs)
• m/classes/clusters/generic.py (modified) (2 diffs)
• m/classes/mask.m (modified) (1 diff)
• m/classes/mask.py (modified) (2 diffs)
• m/classes/model.m (modified) (9 diffs)
• m/classes/model.py (modified) (1 diff)
• m/classes/organizer.m (modified) (2 diffs)
• m/classes/solidearth.m (modified) (4 diffs)
• m/classes/solidearth.py (modified) (5 diffs)
• m/classes/solidearthsettings.m (modified) (6 diffs)
• m/classes/solidearthsettings.py (modified) (5 diffs)
• m/contrib/larour/glacier_inventory.m (modified) (2 diffs)
• m/contrib/larour/ismip6.m (modified) (6 diffs)
• m/contrib/larour/mme_autotime_correlation_matrix.m (modified) (1 diff)
• m/interp/averaging.m (modified) (1 diff)
• m/mesh/FixMesh.m (modified) (3 diffs)
• m/plot/processdatalatlong.m (modified) (1 diff)
• m/solve/solve.m (modified) (2 diffs)
• m/solve/solve.py (modified) (9 diffs)
• m/units/sletogt.m (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -175 +175 @@
         IssmDouble  planetradius=0;
         IssmDouble  planetarea=0;
-        bool            rigid=false;
         bool            elastic=false;
-        bool            rotation=false;
 
         int     numoutputs;
@@ -243 +241 @@
         }
 
-
         /*Deal with dsl multi-model ensembles: {{{*/
         iomodel->FetchData(&dslmodel,"md.dsl.model");
@@ -263 +260 @@
         } /*}}}*/
         /*Deal with elasticity {{{*/
-        iomodel->FetchData(&rigid,"md.solidearth.settings.rigid");
         iomodel->FetchData(&elastic,"md.solidearth.settings.elastic");
-        iomodel->FetchData(&rotation,"md.solidearth.settings.rotation");
-
-        if(elastic | rigid){
-                /*compute green functions for a range of angles*/
-                iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
-                M=reCast<int,IssmDouble>(180./degacc+1.);
-        }
-
-        /*love numbers: */
         if(elastic){
+                /*love numbers: */
                 iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h");
                 iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k");
@@ -281 +269 @@
                 iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk");
                 iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl");
+                parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
 
                 parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1));
@@ -289 +278 @@
                 parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1));
 
+                /*compute elastic green function for a range of angles*/
+                iomodel->FetchData(&degacc,"md.solidearth.settings.degacc");
+                M=reCast<int,IssmDouble>(180./degacc+1.);
+
                 // AD performance is sensitive to calls to ensurecontiguous.
                 // // Providing "t" will cause ensurecontiguous to be called.
                 #ifdef _HAVE_AD_
+                G_rigid=xNew<IssmDouble>(M,"t");
                 G_elastic=xNew<IssmDouble>(M,"t");
                 U_elastic=xNew<IssmDouble>(M,"t");
                 H_elastic=xNew<IssmDouble>(M,"t");
                 #else
+                G_rigid=xNew<IssmDouble>(M);
                 G_elastic=xNew<IssmDouble>(M);
                 U_elastic=xNew<IssmDouble>(M);
                 H_elastic=xNew<IssmDouble>(M);
                 #endif
-        }
-        if(rigid){
-                #ifdef _HAVE_AD_
-                G_rigid=xNew<IssmDouble>(M,"t");
-                #else
-                G_rigid=xNew<IssmDouble>(M);
-                #endif
-        }
-        
-        if(rotation)parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum));
-
-        if(rigid | elastic){
 
                 /*compute combined legendre + love number (elastic green function:*/
                 m=DetermineLocalSize(M,IssmComm::GetComm());
                 GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
-        }
-        if(elastic){
+                // AD performance is sensitive to calls to ensurecontiguous.
+                // // Providing "t" will cause ensurecontiguous to be called.
                 #ifdef _HAVE_AD_
                 G_elastic_local=xNew<IssmDouble>(m,"t");
+                G_rigid_local=xNew<IssmDouble>(m,"t");
                 U_elastic_local=xNew<IssmDouble>(m,"t");
                 H_elastic_local=xNew<IssmDouble>(m,"t");
                 #else
                 G_elastic_local=xNew<IssmDouble>(m);
+                G_rigid_local=xNew<IssmDouble>(m);
                 U_elastic_local=xNew<IssmDouble>(m);
                 H_elastic_local=xNew<IssmDouble>(m);
                 #endif
-        }
-        if(rigid){
-                #ifdef _HAVE_AD_
-                G_rigid_local=xNew<IssmDouble>(m,"t");
-                #else
-                G_rigid_local=xNew<IssmDouble>(m);
-                #endif
-        }
-
-        if(rigid){
+
                 for(int i=lower_row;i<upper_row;i++){
                         IssmDouble alpha,x;
                         alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
+
                         G_rigid_local[i-lower_row]= .5/sin(alpha/2.0);
-                }
-        }
-        if(elastic){
-                for(int i=lower_row;i<upper_row;i++){
-                        IssmDouble alpha,x;
-                        alpha= reCast<IssmDouble>(i)*degacc * PI / 180.0;
-
                         G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row];
                         U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row];
@@ -386 +356 @@
                         }
                 }
-        }
-        if(rigid){
 
                 /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/
@@ -401 +369 @@
                 /*All gather:*/
                 ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-                if(elastic){
-                        ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-                        ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-                        ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
-                }
-                
+                ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+                ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+                ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm());
+
                 /*free resources: */
                 xDelete<int>(recvcounts);
                 xDelete<int>(displs);
+                /*}}}*/
 
                 /*Avoid singularity at 0: */
                 G_rigid[0]=G_rigid[1];
+                parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
                 G_elastic[0]=G_elastic[1];
+                parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
                 U_elastic[0]=U_elastic[1];
+                parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
                 H_elastic[0]=H_elastic[1];
-
-                /*Save our precomputed tables into parameters*/
-                parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M));
-                if(elastic){
-                        parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M));
-                        parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M));
-                        parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
-                }
+                parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M));
 
                 /*free resources: */
+                xDelete<IssmDouble>(love_h);
+                xDelete<IssmDouble>(love_k);
+                xDelete<IssmDouble>(love_l);
+                xDelete<IssmDouble>(love_th);
+                xDelete<IssmDouble>(love_tk);
+                xDelete<IssmDouble>(love_tl);
                 xDelete<IssmDouble>(G_rigid);
                 xDelete<IssmDouble>(G_rigid_local);
-                if(elastic){
-                        xDelete<IssmDouble>(love_h);
-                        xDelete<IssmDouble>(love_k);
-                        xDelete<IssmDouble>(love_l);
-                        xDelete<IssmDouble>(love_th);
-                        xDelete<IssmDouble>(love_tk);
-                        xDelete<IssmDouble>(love_tl);
-                        xDelete<IssmDouble>(G_elastic);
-                        xDelete<IssmDouble>(G_elastic_local);
-                        xDelete<IssmDouble>(U_elastic);
-                        xDelete<IssmDouble>(U_elastic_local);
-                        xDelete<IssmDouble>(H_elastic);
-                        xDelete<IssmDouble>(H_elastic_local);
-                }
+                xDelete<IssmDouble>(G_elastic);
+                xDelete<IssmDouble>(G_elastic_local);
+                xDelete<IssmDouble>(U_elastic);
+                xDelete<IssmDouble>(U_elastic_local);
+                xDelete<IssmDouble>(H_elastic);
+                xDelete<IssmDouble>(H_elastic_local);
         }
-        /*}}}*/
 
         /*Indicate we have not yet run the Geometry Core module: */
         parameters->AddObject(new BoolParam(SealevelriseGeometryDoneEnum,false));
+        /*}}}*/
 
         /*Transitions:{{{ */

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -5446 +5446 @@
 
         /*Computational flags:*/
-        bool computerigid = false;
-        bool computeelastic = false;
+        bool computerigid = true;
+        bool computeelastic = true;
         int  horiz;
 
@@ -5463 +5463 @@
         this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
         this->AddInput(SealevelAreaEnum,&area,P0Enum);
-        if(!computerigid)return;
+        if(!computerigid & !computeelastic)return;
 
         /*recover precomputed green function kernels:*/
@@ -5469 +5469 @@
         parameter->GetParameterValueByPointer((IssmDouble**)&G_rigid_precomputed,&M);
 
+        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
+
+        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
+
+        parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
+        parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
+
         /*allocate indices:*/
         indices=xNew<IssmDouble>(gsize);
         G=xNewZeroInit<IssmDouble>(gsize);
-
-        if(computeelastic){
-                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
-                parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M);
-
-                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter);
-                parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M);
-
-                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter);
-                parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M);
-
-                /*allocate indices:*/
-                GU=xNewZeroInit<IssmDouble>(gsize);
-                if(horiz){
-                        GN=xNewZeroInit<IssmDouble>(gsize);
-                        GE=xNewZeroInit<IssmDouble>(gsize);
-                }
+        GU=xNewZeroInit<IssmDouble>(gsize);
+        if(horiz){
+                GN=xNewZeroInit<IssmDouble>(gsize);
+                GE=xNewZeroInit<IssmDouble>(gsize);
         }
 
@@ -5503 +5499 @@
         /*compute gravity center in lat,long: */
         late= asin(z_element/planetradius);
-        longe = atan2(y_element,x_element); 
-        /*}}}*/
+        longe = atan2(y_element,x_element);
 
         constant=3/rho_earth*area/planetarea;
+
         for(int i=0;i<gsize;i++){
                 IssmDouble alpha;
@@ -5519 +5515 @@
 
                 /*Rigid earth gravitational perturbation: */
-                G[i]+=G_rigid_precomputed[index];
-                if(computeelastic) G[i]+=G_elastic_precomputed[index];
+                if(computerigid){
+                        G[i]+=G_rigid_precomputed[index];
+                }
+                if(computeelastic){
+                        G[i]+=G_elastic_precomputed[index];
+                }
                 G[i]=G[i]*constant;
 
                 /*Elastic components:*/
-                if(computeelastic){
-                        GU[i] =  constant * U_elastic_precomputed[index];
-                        if(horiz){
-                                /*Compute azimuths, both north and east components: */
-                                x = xx[i]; y = yy[i]; z = zz[i];
-                                if(latitude[i]==90){
-                                        x=1e-12; y=1e-12;
-                                }
-                                if(latitude[i]==-90){
-                                        x=1e-12; y=1e-12;
-                                }
-                                dx = x_element-x; dy = y_element-y; dz = z_element-z;
-                                N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-                                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
-
-                                GN[i] = constant*H_elastic_precomputed[index]*N_azim;
-                                GE[i] = constant*H_elastic_precomputed[index]*E_azim;
-                        }
+                GU[i] =  constant * U_elastic_precomputed[index];
+                if(horiz){
+                        /*Compute azimuths, both north and east components: */
+                        x = xx[i]; y = yy[i]; z = zz[i];
+                        if(latitude[i]==90){
+                                x=1e-12; y=1e-12;
+                        }
+                        if(latitude[i]==-90){
+                                x=1e-12; y=1e-12;
+                        }
+                        dx = x_element-x; dy = y_element-y; dz = z_element-z;
+                        N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
+
+                        GN[i] = constant*H_elastic_precomputed[index]*N_azim;
+                        GE[i] = constant*H_elastic_precomputed[index]*E_azim;
                 }
         }
 
         /*Add in inputs:*/
-        this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
-        this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
-        if(computeelastic){
-                this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
-                if(horiz){
-                        this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
-                        this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
-                }
-        }
-        this->inputs->SetDoubleInput(AreaEnum,this->lid,area);
-        this->AddInput(SealevelAreaEnum,&area,P0Enum);
+    this->inputs->SetArrayInput(SealevelriseIndicesEnum,this->lid,indices,gsize);
+    this->inputs->SetArrayInput(SealevelriseGEnum,this->lid,G,gsize);
+    this->inputs->SetArrayInput(SealevelriseGUEnum,this->lid,GU,gsize);
+    if(horiz){
+                this->inputs->SetArrayInput(SealevelriseGNEnum,this->lid,GN,gsize);
+                this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize);
+        }
 
         /*Free allocations:*/
         #ifdef _HAVE_RESTRICT_
+        delete indices;
         delete G;
-        if(computeelastic){
-                delete GU;
-                if(horiz){
-                        delete GN;
-                        delete GE;
-                }
+        delete GU;
+        if(horiz){
+                delete GN;
+                delete GE;
         }
         #else
+        xDelete(indices);
         xDelete(G);
-        if(computeelastic){
-                xDelete(GU);
-                if(horiz){
-                        xDelete(GN);
-                        xDelete(GE);
-                }
+        xDelete(GU);
+        if(horiz){
+                xDelete(GN);
+                xDelete(GE);
         }
         #endif
@@ -5591 +5583 @@
         bool notfullygrounded=false;
         bool scaleoceanarea= false;
-        bool computerigid= false;
+        bool computeelastic= true;
         int  glfraction=1;
 
@@ -5637 +5629 @@
         rho_ice=FindParam(MaterialsRhoIceEnum);
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
-        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
+        this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum);
 
         /*retrieve precomputed G:*/
-        if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
+        if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
         /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5695 +5688 @@
         _assert_(!xIsNan<IssmDouble>(bslrice));
 
-        if(computerigid){
+        if(computeelastic){
                 /*convert from m to kg/m^2:*/
                 I=I*rho_ice*phi;
@@ -5721 +5714 @@
         bool notfullygrounded=false;
         bool scaleoceanarea= false;
-        bool computerigid= false;
+        bool computeelastic= true;
 
         /*elastic green function:*/
@@ -5751 +5744 @@
         rho_water=FindParam(MaterialsRhoSeawaterEnum);
         rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
-        this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum);
+        this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum);
         this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
 
         /*retrieve precomputed G:*/
-        if(computerigid)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
+        if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
 
         /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
@@ -5771 +5764 @@
         _assert_(!xIsNan<IssmDouble>(bslrhydro));
 
-        if(computerigid){
+        if(computeelastic){
                 /*convert from m to kg/m^2:*/
                 W=W*rho_freshwater*phi;
@@ -5795 +5788 @@
         IssmDouble BP;  //change in bottom pressure (Farrel and Clarke, Equ. 4)
         IssmDouble constant;
-        bool computeelastic= false;
+        bool computeelastic= true;
 
         /*elastic green function:*/
@@ -5897 +5890 @@
 
         /*computational flags:*/
-        bool computeelastic= false;
+        bool computeelastic= true;
 
         /*early return if we are not on the ocean or on an ice cap:*/

issm/trunk-jpl/src/m/classes/clusters/generic.py

@@ -1 @@
-import numpy as np
 from socket import gethostname
 from subprocess import call
-from IssmConfig import IssmConfig
-from issmdir import issmdir
-from pairoptions import pairoptions
-from issmssh import issmssh
-from issmscpin import issmscpin
-from issmscpout import issmscpout
-from MatlabFuncs import ispc
+
+import numpy as np
+
 try:
     from generic_settings import generic_settings
 except ImportError:
     print('Warning generic_settings.py not found, default will be used')
+from MatlabFuncs import ispc
+from IssmConfig import IssmConfig
+from issmdir import issmdir
+from issmscpin import issmscpin
+from issmscpout import issmscpout
+from issmssh import issmssh
+from pairoptions import pairoptions
 
 
 class generic(object):
-    """
-    GENERIC cluster class definition
-
-       Usage:
-          cluster = generic('name', 'astrid', 'np', 3)
-          cluster = generic('name', gethostname(), 'np', 3, 'login', 'username')
+    """GENERIC cluster class definition
+
+    Usage:
+        cluster = generic('name', 'astrid', 'np', 3)
+        cluster = generic('name', gethostname(), 'np', 3, 'login', 'username')
     """
 
@@ -50 +51 @@
         except NameError:
             print("generic_settings.py not found, using default settings")
-        # else:
-        #     raise
 
         #OK get other fields

issm/trunk-jpl/src/m/classes/mask.m

@@ -114 +114 @@
                                 WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','name','md.mask.ice_levelset','format','MatArray','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
                         else
-                                WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1);
+                                WriteData(fid,prefix,'object',self,'fieldname','ice_levelset','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
                         end
                 end % }}}

issm/trunk-jpl/src/m/classes/mask.py

@@ -8 +8 @@
 
 class mask(object):
-    '''
-    MASK class definition
+    """MASK class definition
 
-        Usage:
-            mask = mask()
-    '''
+    Usage:
+        mask = mask()
+    """
+    
     def __init__(self):  # {{{
         self.ice_levelset = float('NaN')
@@ -91 +91 @@
 
     def marshall(self, prefix, md, fid):  # {{{
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_levelset', 'format', 'DoubleMat', 'mattype', 1)
-        WriteData(fid, prefix, 'object', self, 'fieldname', 'ice_levelset', 'format', 'DoubleMat', 'mattype', 1)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_levelset', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'ice_levelset', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
     # }}}

issm/trunk-jpl/src/m/classes/model.m

@@ -192 +192 @@
         methods
                 function md = model(varargin) % {{{
-                
 
                         switch nargin
                                 case 0
-                                        md=setdefaultparameters(md,'earth');
+                                        md=setdefaultparameters(md);
+                                case 1
+                                        error('model constructor not supported yet');
+
                                 otherwise
-                                        options=pairoptions(varargin{:});
-                                        planet=getfieldvalue(options,'planet','earth');
-                                        md=setdefaultparameters(md,planet);
-                                end
+                                        error('model constructor error message: 0 of 1 argument only in input.');
+                                end
+
                 end
                 %}}}
@@ -424 +425 @@
                                 if isobject(md.outputdefinition.definitions{i})
                                         %get subfields
-                                        solutionsubfields=fieldnames(md.outputdefinition.definitions{i});
+                                        solutionsubfields=fields(md.outputdefinition.definitions{i});
                                         for j=1:length(solutionsubfields),
                                                 field=md.outputdefinition.definitions{i}.(solutionsubfields{j});
@@ -548 +549 @@
 
                         %loop over model fields
-                        model_fields=fieldnames(md1);
+                        model_fields=fields(md1);
                         for i=1:length(model_fields),
                                 %get field
@@ -554 +555 @@
                                 fieldsize=size(field);
                                 if isobject(field), %recursive call
-                                        object_fields=fieldnames(md1.(model_fields{i}));
+                                        object_fields=fields(md1.(model_fields{i}));
                                         for j=1:length(object_fields),
                                                 %get field
@@ -721 +722 @@
                         if isstruct(md1.results),
                                 md2.results=struct();
-                                solutionfields=fieldnames(md1.results);
+                                solutionfields=fields(md1.results);
                                 for i=1:length(solutionfields),
                                         if isstruct(md1.results.(solutionfields{i}))
@@ -728 +729 @@
                                                 for p=1:length(md1.results.(solutionfields{i}))
                                                     current = md1.results.(solutionfields{i})(p);
-                                                    solutionsubfields=fieldnames(current);
+                                                    solutionsubfields=fields(current);
                                                     for j=1:length(solutionsubfields),
                                                         field=md1.results.(solutionfields{i})(p).(solutionsubfields{j});
@@ -757 +758 @@
                                 if isobject(md1.outputdefinition.definitions{i})
                                         %get subfields
-                                        solutionsubfields=fieldnames(md1.outputdefinition.definitions{i});
+                                        solutionsubfields=fields(md1.outputdefinition.definitions{i});
                                         for j=1:length(solutionsubfields),
                                                 field=md1.outputdefinition.definitions{i}.(solutionsubfields{j});
@@ -1287 +1288 @@
                         end
                 end% }}}
-                function md = setdefaultparameters(md,planet) % {{{
+                function md = setdefaultparameters(md) % {{{
 
                         %initialize subclasses
@@ -1299 +1300 @@
                         md.friction         = friction();
                         md.rifts            = rifts();
-                        md.solidearth       = solidearth(planet);
+                        md.solidearth       = solidearth();
                         md.dsl              = dsl();
                         md.timestepping     = timestepping();

issm/trunk-jpl/src/m/classes/model.py

@@ -173 +173 @@
 
     def __repr__(obj):  #{{{
+        # TODO:
+        # - Convert all formatting to calls to <string>.format (see any 
+        #   already converted <class>.__repr__ method for examples)
+        #
+
         #print "Here %s the number: %d" % ("is", 37)
-        string = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties")
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ..."))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("constants", "[%s %s]" % ("1x1", obj.constants.__class__.__name__), "physical constants"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("smb", "[%s %s]" % ("1x1", obj.smb.__class__.__name__), "surface mass balance"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("basalforcings", "[%s %s]" % ("1x1", obj.basalforcings.__class__.__name__), "bed forcings"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("materials", "[%s %s]" % ("1x1", obj.materials.__class__.__name__), "material properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("damage", "[%s %s]" % ("1x1", obj.damage.__class__.__name__), "damage propagation laws"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("friction", "[%s %s]" % ("1x1", obj.friction.__class__.__name__), "basal friction / drag properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("flowequation", "[%s %s]" % ("1x1", obj.flowequation.__class__.__name__), "flow equations"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("timestepping", "[%s %s]" % ("1x1", obj.timestepping.__class__.__name__), "time stepping for transient models"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("initialization", "[%s %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("rifts", "[%s %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("solidearth", "[%s %s]" % ("1x1", obj.solidearth.__class__.__name__), "solidearth inputs and settings"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("dsl", "[%s %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("debug", "[%s %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("verbose", "[%s %s]" % ("1x1", obj.verbose.__class__.__name__), "verbosity level in solve"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("settings", "[%s %s]" % ("1x1", obj.settings.__class__.__name__), "settings properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("toolkits", "[%s %s]" % ("1x1", obj.toolkits.__class__.__name__), "PETSc options for each solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("cluster", "[%s %s]" % ("1x1", obj.cluster.__class__.__name__), "cluster parameters (number of cpus...)"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("balancethickness", "[%s %s]" % ("1x1", obj.balancethickness.__class__.__name__), "parameters for balancethickness solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("stressbalance", "[%s %s]" % ("1x1", obj.stressbalance.__class__.__name__), "parameters for stressbalance solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("groundingline", "[%s %s]" % ("1x1", obj.groundingline.__class__.__name__), "parameters for groundingline solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("hydrology", "[%s %s]" % ("1x1", obj.hydrology.__class__.__name__), "parameters for hydrology solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("masstransport", "[%s %s]" % ("1x1", obj.masstransport.__class__.__name__), "parameters for masstransport solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("thermal", "[%s %s]" % ("1x1", obj.thermal.__class__.__name__), "parameters for thermal solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("steadystate", "[%s %s]" % ("1x1", obj.steadystate.__class__.__name__), "parameters for steadystate solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("transient", "[%s %s]" % ("1x1", obj.transient.__class__.__name__), "parameters for transient solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution"))
-        string = "%s\n%s" % (string, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution"))
-        string = "%s\n%s" % (string, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("qmu", "[%s %s]" % ("1x1", obj.qmu.__class__.__name__), "dakota properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("amr", "[%s %s]" % ("1x1", obj.amr.__class__.__name__), "adaptive mesh refinement properties"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("outputdefinition", "[%s %s]" % ("1x1", obj.outputdefinition.__class__.__name__), "output definition"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("results", "[%s %s]" % ("1x1", obj.results.__class__.__name__), "model results"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay"))
-        string = "%s\n%s" % (string, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields"))
-        return string
+        s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties")
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ..."))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("constants", "[%s %s]" % ("1x1", obj.constants.__class__.__name__), "physical constants"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("smb", "[%s %s]" % ("1x1", obj.smb.__class__.__name__), "surface mass balance"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("basalforcings", "[%s %s]" % ("1x1", obj.basalforcings.__class__.__name__), "bed forcings"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("materials", "[%s %s]" % ("1x1", obj.materials.__class__.__name__), "material properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("damage", "[%s %s]" % ("1x1", obj.damage.__class__.__name__), "damage propagation laws"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("friction", "[%s %s]" % ("1x1", obj.friction.__class__.__name__), "basal friction / drag properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("flowequation", "[%s %s]" % ("1x1", obj.flowequation.__class__.__name__), "flow equations"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("timestepping", "[%s %s]" % ("1x1", obj.timestepping.__class__.__name__), "time stepping for transient models"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("initialization", "[%s %s]" % ("1x1", obj.initialization.__class__.__name__), "initial guess / state"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("rifts", "[%s %s]" % ("1x1", obj.rifts.__class__.__name__), "rifts properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("solidearth", "[%s %s]" % ("1x1", obj.solidearth.__class__.__name__), "solidearth inputs and settings"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("dsl", "[%s %s]" % ("1x1", obj.dsl.__class__.__name__), "dynamic sea level"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("debug", "[%s %s]" % ("1x1", obj.debug.__class__.__name__), "debugging tools (valgrind, gprof)"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("verbose", "[%s %s]" % ("1x1", obj.verbose.__class__.__name__), "verbosity level in solve"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("settings", "[%s %s]" % ("1x1", obj.settings.__class__.__name__), "settings properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("toolkits", "[%s %s]" % ("1x1", obj.toolkits.__class__.__name__), "PETSc options for each solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("cluster", "[%s %s]" % ("1x1", obj.cluster.__class__.__name__), "cluster parameters (number of cpus...)"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("balancethickness", "[%s %s]" % ("1x1", obj.balancethickness.__class__.__name__), "parameters for balancethickness solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("stressbalance", "[%s %s]" % ("1x1", obj.stressbalance.__class__.__name__), "parameters for stressbalance solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("groundingline", "[%s %s]" % ("1x1", obj.groundingline.__class__.__name__), "parameters for groundingline solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("hydrology", "[%s %s]" % ("1x1", obj.hydrology.__class__.__name__), "parameters for hydrology solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("masstransport", "[%s %s]" % ("1x1", obj.masstransport.__class__.__name__), "parameters for masstransport solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("thermal", "[%s %s]" % ("1x1", obj.thermal.__class__.__name__), "parameters for thermal solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("steadystate", "[%s %s]" % ("1x1", obj.steadystate.__class__.__name__), "parameters for steadystate solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("transient", "[%s %s]" % ("1x1", obj.transient.__class__.__name__), "parameters for transient solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("levelset", "[%s %s]" % ("1x1", obj.levelset.__class__.__name__), "parameters for moving boundaries (level - set method)"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("calving", "[%s %s]" % ("1x1", obj.calving.__class__.__name__), "parameters for calving"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("frontalforcings", "[%s %s]" % ("1x1", obj.frontalforcings.__class__.__name__), "parameters for frontalforcings"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("gia", "[%s %s]" % ("1x1", obj.gia.__class__.__name__), "parameters for gia solution"))
+        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("esa", "[%s %s]" % ("1x1", obj.esa.__class__.__name__), "parameters for elastic adjustment solution"))
+        s = "%s\n%s" % (s, '%19s: %-22s -- %s' % ("love", "[%s %s]" % ("1x1", obj.love.__class__.__name__), "parameters for love solution"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("autodiff", "[%s %s]" % ("1x1", obj.autodiff.__class__.__name__), "automatic differentiation parameters"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("inversion", "[%s %s]" % ("1x1", obj.inversion.__class__.__name__), "parameters for inverse methods"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("qmu", "[%s %s]" % ("1x1", obj.qmu.__class__.__name__), "dakota properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("amr", "[%s %s]" % ("1x1", obj.amr.__class__.__name__), "adaptive mesh refinement properties"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("outputdefinition", "[%s %s]" % ("1x1", obj.outputdefinition.__class__.__name__), "output definition"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("results", "[%s %s]" % ("1x1", obj.results.__class__.__name__), "model results"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay"))
+        s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields"))
+        return s
     # }}}
 

issm/trunk-jpl/src/m/classes/organizer.m

@@ -156 +156 @@
                         error(['Could not find ' path ]);
                 end%}}}
+                function loaddatanoprefix(org,string),% {{{
+
+                        %Get model path
+                        if ~ischar(string), error('argument provided is not a string'); end
+                        path=[org.repository '/' string];
+
+                        %figure out if the data is there, otherwise, we have to use the default path supplied by user.
+                        if exist(path,'file'),
+                                path=path;
+                        elseif exist([path '.mat'],'file'),
+                                path=[path '.mat'];
+                        else
+                                error(['Could not find ' path ]);
+                        end
+                        if exist(path,'file')
+                                evalin('caller',['load -mat ' path]);
+                                return;
+                        end
+
+                        %If we are here, the data has not been found. 
+                        error(['Could not find ' path ]);
+                end%}}}
                 function bool=perform(org,varargin) % {{{
 
@@ -254 +276 @@
                         eval(['save(''' name '''' variables ',''-v7.3'');']);
                 end%}}}
+                function savedatanoprefix(org,varargin) % {{{
+
+                        %check
+                        if (org.currentstep==0), error('Cannot save data because organizer (org) is empty! Make sure you did not skip any perform call'); end
+                        if (org.currentstep>length(org.steps)), error('Cannot save data because organizer (org) is not up to date!'); end
+
+                        name=[org.repository '/' org.steps(org.currentstep).string ];
+                        disp(['saving data in: ' name]);
+
+                        %Skip if requested
+                        if org.skipio,
+                                disp(['WARNING: Skipping saving ' name]);
+                                return;
+                        end
+
+                        %check that md is a model
+                        if (org.currentstep>length(org.steps)), error(['organizer error message: element with id ' num2str(org.currentstep) ' not found']); end
+
+                        %list of variable names: 
+                        variables='';
+                        for i=2:nargin, 
+                                variables=[variables ',' '''' inputname(i) ''''];
+                                eval([inputname(i) '= varargin{' num2str(i-1) '};']);
+                        end
+                        eval(['save(''' name '''' variables ',''-v7.3'');']);
+                end%}}}
         end
 end

issm/trunk-jpl/src/m/classes/solidearth.m

@@ -14 +14 @@
                 requested_outputs      = {};
                 transitions            = {};
-                partitionice              = [];
-                partitionhydro             = [];
+                partitionice           = [];
+                partitionhydro         = [];
         end
         methods
@@ -21 +21 @@
                         switch nargin
                                 case 0
-                                        self=setdefaultparameters(self,earth);
-                                case 1
-                                        self=setdefaultparameters(self,varargin{:});
+                                        self=setdefaultparameters(self);
                                 otherwise
-                                        error('solidearth constructor error message: zero or one argument  only!');
+                                        error('constructor not supported');
                         end
                 end % }}}
-                function self = setdefaultparameters(self,planet) % {{{
+                function self = setdefaultparameters(self) % {{{
                 
                 %output default:
@@ -35 +33 @@
                 %transitions should be a cell array of vectors: 
                 self.transitions={};
-
+                
                 %no partitions requested for barystatic contribution: 
                 self.partitionice=[];
@@ -41 +39 @@
 
                 %earth radius
-                self.planetradius= planetradius(planet);
+                self.planetradius= planetradius('earth');
                 
                 end % }}}

issm/trunk-jpl/src/m/classes/solidearth.py

@@ -1 +1 @@
 import numpy as np
+
 from checkfield import checkfield
 from fielddisplay import fielddisplay
@@ -19 +20 @@
 
     def __init__(self, *args):  #{{{
-        self.sealevel = np.nan
-        self.settings = solidearthsettings()
-        self.surfaceload = surfaceload()
-        self.lovenumbers = lovenumbers()
-        self.rotational = rotational()
-        self.planetradius = planetradius('earth')
-        self.requested_outputs = ['default']
-        self.transitions = []
+        self.sealevel           = np.nan
+        self.settings           = solidearthsettings()
+        self.surfaceload        = surfaceload()
+        self.lovenumbers        = lovenumbers()
+        self.rotational         = rotational()
+        self.planetradius       = planetradius('earth')
+        self.requested_outputs  = []
+        self.transitions        = []
+        self.partitionice       = []
+        self.partitionhydro     = []
 
         nargin = len(args)
@@ -50 +53 @@
 
     def setdefaultparameters(self):  # {{{
-        return self
+        # Default output
+        self.requested_outputs = ['default']
+
+        # Transitions should be a list
+        self.transitions = []
+
+        # No partitions requested for barystatic contribution
+        self.partitionice = []
+        self.partitionhydro = []
+
+        # Earth radius
+        self.planetradius = planetradius('earth')
     #}}}
 
@@ -58 +72 @@
         md = checkfield(md, 'fieldname', 'solidearth.sealevel', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
         md = checkfield(md, 'fieldname', 'solidearth.requested_outputs', 'stringrow', 1)
+
         self.settings.checkconsistency(md, solution, analyses)
         self.surfaceload.checkconsistency(md, solution, analyses)
@@ -73 +88 @@
         WriteData(fid, prefix, 'object', self, 'fieldname', 'planetradius', 'format', 'Double')
         WriteData(fid, prefix, 'object', self, 'fieldname', 'transitions', 'format', 'MatArray')
+
+        if len(self.partitionice):
+            npartice = np.max(self.partitionice) + 2
+        else:
+            npartice = 0
+
+        if len(self.partitionhydro):
+            nparthydro = np.max(self.partitionhydro) + 2
+        else:
+            nparthydro = 0
+
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionice', 'mattype', 1, 'format', 'DoubleMat');
+        WriteData(fid, prefix, 'data', npartice, 'format', 'Integer', 'name', 'md.solidearth.npartice');
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'partitionhydro', 'mattype', 1, 'format', 'DoubleMat');
+        WriteData(fid, prefix, 'data', nparthydro,'format', 'Integer', 'name','md.solidearth.nparthydro');
 
         self.settings.marshall(prefix, md, fid)

issm/trunk-jpl/src/m/classes/solidearthsettings.m

@@ -17 +17 @@
                 degacc                 = 0; %degree increment for resolution of Green tables
                 horiz                  = 0; %compute horizontal deformation
+                glfraction             = 1; %barystatic contribution full or fractional (default fractional)
         end
         methods
@@ -48 +49 @@
                 %how many time steps we skip before we run solidearthsettings solver during transient
                 self.runfrequency=1;
+                
+                %fractional contribution: 
+                self.glfraction=1;
         
                 %horizontal displacement?  (not by default)
@@ -63 +67 @@
                         md = checkfield(md,'fieldname','solidearth.settings.degacc','size',[1 1],'>=',1e-10);
                         md = checkfield(md,'fieldname','solidearth.settings.horiz','NaN',1,'Inf',1,'values',[0 1]);
-
-                        %checks on computational flags: 
-                        if self.elastic==1 & self.rigid==0,
-                                error('solidearthsettings checkconsistency error message: need rigid on if elastic flag is set');
-                        end
+                        md = checkfield(md,'fieldname','solidearth.settings.glfraction','values',[0 1]);
 
                         %a coupler to a planet model is provided. 
@@ -96 +96 @@
                         fielddisplay(self,'rotation','earth rotational potential perturbation');
                         fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
+                        fielddisplay(self,'glfraction','contribute fractionally (default, 1) to barystatic sea level');
                 end % }}}
                 function marshall(self,prefix,md,fid) % {{{
@@ -109 +110 @@
                         WriteData(fid,prefix,'object',self,'fieldname','horiz','name','md.solidearth.settings.horiz','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','computesealevelchange','name','md.solidearth.settings.computesealevelchange','format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','glfraction','name','md.solidearth.settings.glfraction','format','Integer');
                 end % }}}
                 function savemodeljs(self,fid,modelname) % {{{
@@ -121 +123 @@
                         writejsdouble(fid,[modelname '.slr.settings.run_frequency'],self.run_frequency);
                         writejsdouble(fid,[modelname '.slr.settings.degacc'],self.degacc);
+                        writejsdouble(fid,[modelname '.slr.settings.glfraction'],self.glfraction);
                 end % }}}
                 function self = extrude(self,md) % {{{

issm/trunk-jpl/src/m/classes/solidearthsettings.py

@@ -14 +14 @@
 
     def __init__(self, *args): #{{{
+        self.reltol                 = 0
+        self.abstol                 = 0
+        self.maxiter                = 0
+        self.rigid                  = 0
+        self.elastic                = 0
+        self.rotation               = 0
+        self.ocean_area_scaling     = 0
+        self.runfrequency           = 1 # How many time steps we skip before we run grd_core
+        self.computesealevelchange  = 0 # Will grd_core compute sea level?
+        self.degacc                 = 0 # Degree increment for resolution of Green tables
+        self.horiz                  = 0 # Compute horizontal displacement?
+        self.glfraction             = 1 # Barystatic contribution: full or fractional (default: fractional)
+
         nargin = len(args)
 
@@ -33 +46 @@
         s += '{}\n'.format(fielddisplay(self, 'elastic', 'elastic earth graviational potential perturbation'))
         s += '{}\n'.format(fielddisplay(self, 'degacc', 'accuracy (default .01 deg) for numerical discretization of the Green\'s functions'))
+        s += '{}\n'.format(fielddisplay(self, 'glfraction', 'contribute fractionally (default, 1) to barystatic sea level'))
         return s
     #}}}
@@ -59 +73 @@
         # Horizontal displacement? (not on by default)
         self.horiz = 0
-        return self
+
+        # Fractional contribution
+        self.glfraction = 1
     #}}}
 
@@ -71 +87 @@
         md = checkfield(md, 'fieldname', 'solidearth.settings.degacc', 'size', [1], '>=', 1e-10)
         md = checkfield(md, 'fieldname', 'solidearth.settings.horiz', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
+        md = checkfield(md, 'fieldname', 'solidearth.settings.glfraction', 'values', [0, 1])
 
         # A coupler to planet model is provided
@@ -88 +105 @@
     def marshall(self, prefix, md, fid): #{{{
         WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double')
-        WriteData(fid,prefix,'object', self,'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
-        WriteData(fid,prefix,'object', self,'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
-        WriteData(fid,prefix,'object', self,'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean')
-        WriteData(fid,prefix,'object', self,'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
-        WriteData(fid,prefix,'object', self,'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
-        WriteData(fid,prefix,'object', self,'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean')
-        WriteData(fid,prefix,'object', self,'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
-        WriteData(fid,prefix,'object', self,'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
-        WriteData(fid,prefix,'object', self,'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
-        WriteData(fid,prefix,'object', self,'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'rigid', 'name', 'md.solidearth.settings.rigid', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname', 'computesealevelchange', 'name', 'md.solidearth.settings.computesealevelchange', 'format', 'Integer')
+        WriteData(fid, prefix, 'object', self, 'fieldname','glfraction', 'name', 'md.solidearth.settings.glfraction', 'format', 'Integer')
     #}}}
 

issm/trunk-jpl/src/m/contrib/larour/glacier_inventory.m

@@ -183 +183 @@
                         %Go through O2 regions: 
                         for i=subsetregions,
+                        %for i=33,
                                 string=self.boxes(i).RGI_CODE; 
                                 disp(['progressing with region ' num2str(i) ' ' string]);
@@ -208 +209 @@
                                                 case 19, radius=60;
                                                 case 32, radius=60;
-                                                case 33, radius=5;
+                                                case 33, radius=10;
                                                 case 41, radius=75;
                                                 case 42, radius=45;

issm/trunk-jpl/src/m/contrib/larour/ismip6.m

@@ -13 +13 @@
                 directories    = {};   %directories where the files are
                 experiments    = {};   %names of experiments
+                base           = {};   %placeholder for base
+                surface        = {};   %placeholder for surface
                 thickness      = {};   %placeholder for thicknesses
                 deltathickness = {};   %placeholder for delta thicknesses
+                deltathicknessvaf = {};   %placeholder for delta thicknesses above floatation
+                deltathicknesshal = {};   %placeholder for delta thicknesses halosteric origins
+                deltathicknessbar = {};   %placeholder for delta thicknesses halosteric origins
+                thicknesscorrection={};  
                 icemask        = {};   %placeholder for ice masks
                 oceanmask      = {};   %placeholder for ocean masks
@@ -20 +26 @@
                 timestart      = {};   %placeholder for times
                 calendar      = {};   %placeholder for times
+                di            = {}; %ice densities
         end
         methods
@@ -146 +153 @@
                         end
 
+                        if ~exist('timestart','var'), timestart=2015; end
+                        if ~exist('calendar','var'), calendar=0; end
+
                 end % }}}
                 function info=readinfo(self,experiment,field) % {{{
@@ -203 +213 @@
 
                                 %map onto mesh: correct only for thicknesses
-                                hg=ismip2mesh_correction.*(ismip2mesh*ht) ;
+                                if strcmpi(field,'lithk') | strcmpi(field,'orog') | strcmpi(field,'base'),
+                                        hg=ismip2mesh_correction.*(ismip2mesh*ht) ;
+                                        %hg=ismip2mesh*ht ;
+                                else
+                                        hg=ismip2mesh*ht ;
+                                end
 
                                 %keep field:
@@ -209 +224 @@
                                         pos=find(isnan(hg)); hg(pos)=0;
                                         self.thickness{i}=hg; 
+                                end
+                                if strcmpi(field,'orog'),
+                                        pos=find(isnan(hg)); hg(pos)=0;
+                                        self.surface{i}=hg; 
+                                end
+                                if strcmpi(field,'base'),
+                                        pos=find(isnan(hg)); hg(pos)=0;
+                                        self.base{i}=hg; 
                                 end
                                 if strcmpi(field,'sftgif'),
@@ -220 +243 @@
                                 end
                                 if strcmpi(field,'sftgrf'),
-                                        hge=-ones(md.mesh.numberofvertices,size(hg,2));
+                                        hgv=-ones(md.mesh.numberofvertices,size(hg,2));
                                         for j=1:size(hg,2),
                                                 hgj=hg(:,j);
-                                                pos=find(hgj>0);
-                                                hge(md.mesh.elements(pos,:),j)=1;
-                                        end
-                                        self.oceanmask{i}=hge; 
+                                                pos=find(hgj>.99); %we want fully grounded
+                                                %pos=find(hgj>0); %we want slightly grounded
+                                                hgv(md.mesh.elements(pos,:),j)=1;
+                                        end
+                                        self.oceanmask{i}=hgv; 
                                 end
 

issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m

@@ -1 @@
-function matrix=mme_time_correlation_matrix(mme,varargin); 
+function matrix=mme_autotime_correlation_matrix(mme,type)
 
-        if nargin==2,
-                type=varargin{1};
-        elseif nargin==3,
-                mme2=varargin{1};
-                type=varargin{2};
-        else 
-                error('mme_time_correlation_matrix usage error: 2 or 3 arguments only allowed!');
+
+        %Out of a multi model ensemble (nsamples x nsteps) of runs, build 
+        %a temporal correlation matrix (of size nsteps x nsteps)
+
+        nsamples=size(mme,1);
+        nsteps=size(mme,2);
+
+        %initialize with 1 in the diagonal: 
+        matrix=eye(nsteps,nsteps);
+
+        %go through time steps, and fill up the top part. 
+        for i=1:nsteps,
+                for j=i+1:nsteps,
+                        matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type);
+                        matrix(j,i)=matrix(i,j);
+                end
         end
 
-        if nargin==2,
-
-                %Out of a multi model ensemble (nsamples x nsteps) of runs, build 
-                %a temporal correlation matrix (of size nsteps x nsteps)
-
-                nsamples=size(mme,1);
-                nsteps=size(mme,2);
-
-                %initialize with 1 in the diagonal: 
-                matrix=eye(nsteps,nsteps);
-
-                %go through time steps, and fill up the top part. 
-                for i=1:nsteps,
-                        for j=i+1:nsteps,
-                                matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type);
-                                matrix(j,i)=matrix(i,j);
-                        end
-                end
-        else
-
-                %Same kind of computations, except it's not autocorrelation:
-                nsamples=size(mme,1); nsamples2=size(mme2,1);
-                nsteps=size(mme,2); nsteps2=size(mme2,2);
-
-                if nsteps2~=nsteps,
-                        error('number of time steps from both sample matrices should be identical!');
-                end
-                if nsamples2~=nsamples,
-                        error('number of samples from both sample matrices should be identical!');
-                end
-
-                %initialize with 1 in the diagonal: 
-                matrix=zeros(nsteps,nsteps);
-
-                %go through time steps, and fill up the top part. 
-                for i=1:nsteps,
-                        for j=i:nsteps,
-                                matrix(i,j)=corr(mme(:,i),mme2(:,j),'Type',type);
-                                matrix(j,i)=matrix(i,j);
-                        end
-                end
-        end

issm/trunk-jpl/src/m/interp/averaging.m

@@ -64 +64 @@
 else
         rep=3;
-        areas=GetAreas(index,md.mesh.x2d,md.mesh.y2d);
+        if isa(md.mesh,'mesh3dsurface'),
+                areas=GetAreas3DTria(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.z);
+        else
+                areas=GetAreas(index,md.mesh.x,md.mesh.y);
+        end
 end
 summation=1/rep*ones(rep,1);

issm/trunk-jpl/src/m/mesh/FixMesh.m

@@ -1 @@
-function  [index2 x2 y2 value2]=FixMesh(index,x,y,value)
+function  [index2 x2 y2 value2 newpos]=FixMesh(index,x,y,value)
 % FIXMESH - FixMesh fix mesh with broken triangles, orphan vertices, etc ...
 %
@@ -15 +15 @@
 y2=y;
 value2=value;
+newpos=1:length(x);
 
 %First, look for orphan vertices, and take them out.
@@ -29 +30 @@
         y2(orphan)=[];
         value2(orphan)=[];
+        newpos(orphan)=[];
 
         %now, the index:

issm/trunk-jpl/src/m/plot/processdatalatlong.m

@@ -22 +22 @@
                 %interpolate data: 
                 extradata=griddata(x0,y0,data,xextra,yextra,'nearest');
-
                 data=[data; extradata];
         elseif length(data)==length(md.mesh.elements),
-                error('processdatalatlong error message: coord ''latlong'' case not covered for element data ');
                 datatype=1;
         end

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

@@ -4 +4 @@
 %   Usage:
 %      md=solve(md,solutionstring,varargin)
-%      where varargin is a lit of paired arguments of string OR enums
 %
-%   solution types available comprise:
-%      - 'Stressbalance'      or 'sb'
-%      - 'Masstransport'      or 'mt'
-%      - 'Thermal'            or 'th'
-%      - 'Steadystate'        or 'ss'
-%      - 'Transient'          or 'tr'
-%      - 'Balancethickness'   or 'mc'
-%      - 'Balancevelocity'    or 'bv'
-%      - 'BedSlope'           or 'bsl'
-%      - 'SurfaceSlope'       or 'ssl'
-%      - 'Hydrology'          or 'hy'
-%      - 'DamageEvolution'    or 'da'
-%      - 'Gia'                or 'gia'
-%      - 'Love'               or 'lv'
-%      - 'Esa'                or 'esa'
-%      - 'Sealevelrise'       or 'slr'
+%   where varargin is a lit of paired arguments of string OR enums
 %
-%  extra options:
-%      - loadonly    : does not solve. only load results
-%      - runtimename : true or false (default is true), makes name unique
-%      - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on consistency of model
-%      - restart: 'directory name (relative to the execution directory) where the restart file is located.
-%      - outbinread  : if 0, download the outbin but do not process is (md.results is not updated)
+%   Solution types available comprise:
+%   - 'Stressbalance'      or 'sb'
+%   - 'Masstransport'      or 'mt'
+%   - 'Thermal'            or 'th'
+%   - 'Steadystate'        or 'ss'
+%   - 'Transient'          or 'tr'
+%   - 'Balancethickness'   or 'mc'
+%   - 'Balancevelocity'    or 'bv'
+%   - 'BedSlope'           or 'bsl'
+%   - 'SurfaceSlope'       or 'ssl'
+%   - 'Hydrology'          or 'hy'
+%   - 'DamageEvolution'    or 'da'
+%   - 'Gia'                or 'gia'
+%   - 'Love'               or 'lv'
+%   - 'Esa'                or 'esa'
+%   - 'Sealevelrise'       or 'slr'
+%
+%   Extra options:
+%   - loadonly         : does not solve. only load results
+%   - runtimename      : true or false (default is true), makes name unique
+%   - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on 
+%                        consistency of model
+%   - restart          : 'directory name (relative to the execution directory) 
+%                        where the restart file is located
 %
 %   Examples:
@@ -111 +113 @@
 end
 
-%if running qmu analysis, some preprocessing of dakota files using models fields needs to be carried out. 
+%if running QMU analysis, some preprocessing of Dakota files using models fields needs to be carried out. 
 if md.qmu.isdakota,
         md=preqmu(md,options);

issm/trunk-jpl/src/m/solve/solve.py

@@ -1 @@
-import datetime
+from datetime import datetime
 import os
 import shutil
@@ -16 +16 @@
     Usage:
         md = solve(md, solutionstring, varargin)
-        where varargin is a list of paired arguments of string OR enums
+    
+    where varargin is a list of paired arguments of string OR enums
 
-        solution types available comprise:
-            - 'Stressbalance'      or 'sb'
-            - 'Masstransport'      or 'mt'
-            - 'Thermal'            or 'th'
-            - 'Steadystate'        or 'ss'
-            - 'Transient'          or 'tr'
-            - 'Balancethickness'   or 'mc'
-            - 'Balancevelocity'    or 'bv'
-            - 'BedSlope'           or 'bsl'
-            - 'SurfaceSlope'       or 'ssl'
-            - 'Hydrology'          or 'hy'
-            - 'DamageEvolution'    or 'da'
-            - 'Gia'                or 'gia'
-            - 'Esa'                or 'esa'
-            - 'Sealevelrise'       or 'slr'
-            - 'Love'               or 'lv'
+    Solution types available comprise:
+    - 'Stressbalance'      or 'sb'
+    - 'Masstransport'      or 'mt'
+    - 'Thermal'            or 'th'
+    - 'Steadystate'        or 'ss'
+    - 'Transient'          or 'tr'
+    - 'Balancethickness'   or 'mc'
+    - 'Balancevelocity'    or 'bv'
+    - 'BedSlope'           or 'bsl'
+    - 'SurfaceSlope'       or 'ssl'
+    - 'Hydrology'          or 'hy'
+    - 'DamageEvolution'    or 'da'
+    - 'Gia'                or 'gia'
+    - 'Esa'                or 'esa'
+    - 'Sealevelrise'       or 'slr'
+    - 'Love'               or 'lv'
 
-        extra options:
-            - loadonly          : does not solve. only load results
-            - checkconsistency  : 'yes' or 'no' (default is 'yes'), ensures 
-                                  checks on consistency of model
-            - restart           : 'directory name (relative to the execution 
-                                  directory) where the restart file is located
+    Extra options:
+    - loadonly         : does not solve. only load results
+    - runtimename      : true or false (default is true), makes name unique
+    - checkconsistency : 'yes' or 'no' (default is 'yes'), ensures checks on 
+                         consistency of model
+    - restart          : directory name (relative to the execution directory) 
+                         where the restart file is located
 
     Examples:
@@ -47 +49 @@
     """
 
-    #recover and process solve options
+    # Recover and process solve options
     if solutionstring.lower() == 'sb' or solutionstring.lower() == 'stressbalance':
         solutionstring = 'StressbalanceSolution'
@@ -82 +84 @@
     options = pairoptions('solutionstring', solutionstring, *args)
 
-    #recover some fields
+    # Recover some fields
     md.private.solution = solutionstring
     cluster = md.cluster
@@ -90 +92 @@
         batch = 0
 
-    #check model consistency
+    # Check model consistency
     if options.getfieldvalue('checkconsistency', 'yes') == 'yes':
         if md.verbose.solution:
@@ -96 +98 @@
         ismodelselfconsistent(md)
 
-    #First, build a runtime name that is unique
+    # First, build a runtime name that is unique
     restart = options.getfieldvalue('restart', '')
     if restart == 1:
-        pass  #do nothing
+        pass  # do nothing
     else:
         if restart:
@@ -105 +107 @@
         else:
             if options.getfieldvalue('runtimename', True):
-                c = datetime.datetime.now()
+                c = datetime.now()
                 md.private.runtimename = "%s-%02i-%02i-%04i-%02i-%02i-%02i-%i" % (md.miscellaneous.name, c.month, c.day, c.year, c.hour, c.minute, c.second, os.getpid())
             else:
                 md.private.runtimename = md.miscellaneous.name
 
-    #if running qmu analysis, some preprocessing of dakota files using models
-    #fields needs to be carried out.
+    # If running QMU analysis, some preprocessing of Dakota files using models 
+    # fields needs to be carried out
     if md.qmu.isdakota:
         md = preqmu(md, options)
 
-    #Do we load results only?
+    # Do we load results only?
     if options.getfieldvalue('loadonly', False):
         md = loadresultsfromcluster(md)
         return md
 
-    #Write all input files
-    marshall(md)  # bin file
-    md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits')  # toolkits file
-    cluster.BuildQueueScript(md.private.runtimename, md.miscellaneous.name, md.private.solution, md.settings.io_gather, md.debug.valgrind, md.debug.gprof, md.qmu.isdakota, md.transient.isoceancoupling)  # queue file
+    # Write all input files
+    marshall(md) # bin file
+    md.toolkits.ToolkitsFile(md.miscellaneous.name + '.toolkits') # toolkits file
+    cluster.BuildQueueScript(md.private.runtimename, md.miscellaneous.name, md.private.solution, md.settings.io_gather, md.debug.valgrind, md.debug.gprof, md.qmu.isdakota, md.transient.isoceancoupling) # queue file
 
-    #Stop here if batch mode
+    # Stop here if batch mode
     if options.getfieldvalue('batch', 'no') == 'yes':
         print('batch mode requested: not launching job interactively')
@@ -131 +133 @@
         return md
 
-    #Upload all required files:
+    # Upload all required files:
     modelname = md.miscellaneous.name
     filelist = [modelname + '.bin ', modelname + '.toolkits ', modelname + '.queue ']
@@ -140 +142 @@
         cluster.UploadQueueJob(md.miscellaneous.name, md.private.runtimename, filelist)
 
-    #Launch job
+    # Launch job
     cluster.LaunchQueueJob(md.miscellaneous.name, md.private.runtimename, filelist, restart, batch)
 
-    #wait on lock
+    # Wait on lock
     if md.settings.waitonlock > 0:
-        #we wait for the done file
         islock = waitonlock(md)
-        if islock == 0:  #no results to be loaded
+        if islock == 0:  # no results to be loaded
             print('The results must be loaded manually with md = loadresultsfromcluster(md).')
-        else:  #load results
+        else:  # load results
             if md.verbose.solution:
                 print('loading results from cluster')
             md = loadresultsfromcluster(md)
+
     return md

issm/trunk-jpl/src/m/units/sletogt.m

@@ -1 +1 @@
 function conversionfactor=sletogt()
-        
-        conversionfactor=361.9;  %361.9 Gigatons to 1 mm sea-level equivalent.
+
+        rho_water=1023;
+        conversionfactor=rho_water/1000*361.9; %361.9 Gigatons to 1 mm sea-level equivalent.