Changeset 24505

Timestamp:

01/09/20 11:17:17 (5 years ago)

Author:

Eric.Larour

Message:

CHG: new dslmme class for multi-model ensembles.
New bottom pressure capabilities for dsl and dslmme classes.

Location:

issm/trunk-jpl/src

Files:

• c/analyses/SealevelriseAnalysis.cpp (modified) (3 diffs)
• c/classes/Elements/Tria.cpp (modified) (6 diffs)
• c/classes/Elements/Tria.h (modified) (1 diff)
• c/cores/sealevelrise_core.cpp (modified) (3 diffs)
• c/shared/Enum/Enum.vim (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (10 diffs)
• m/classes/dsl.m (modified) (6 diffs)
• m/classes/dslmme.m (added)

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

@@ -58 +58 @@
         /*dynamic sea level: */
         iomodel->FetchData(&dslmodel,"md.dsl.model");
-        if (dslmodel==1){
+        if (dslmodel==1){ /*standard dsl model:{{{*/
 
                 /*deal with global mean steric rate: */
@@ -94 +94 @@
 
                 /*deal with dynamic sea level fields: */
-
                 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_surface_height_change_above_geoid", DslSeaSurfaceHeightChangeAboveGeoidEnum);
                 iomodel->FetchDataToInput(inputs2,elements,"md.dsl.sea_water_pressure_change_at_sea_floor", DslSeaWaterPressureChangeAtSeaFloor);
                 
-                //TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,counter, times,N);
-        }
+        } /*}}}*/
+        else if (dslmodel==2){ /*multi-model ensemble dsl model:{{{*/
+        
+                /*variables:*/
+                int nummodels;
+                IssmDouble** pstr=NULL; 
+                IssmDouble*  str=NULL;
+                IssmDouble*  times = NULL;
+                int* pM = NULL;
+                int* pN = NULL;
+                int M,N;
+
+                /*deal with dsl.sea_surface_height_change_above_geoid {{{*/
+                iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.global_average_thermosteric_sea_level_change");
+
+                /*go through the mat array and create a dataset of transient inputs:*/
+                for (int i=0;i<nummodels;i++){
+
+                        M=pM[i];
+                        N=pN[i];
+                        str=pstr[i];
+
+                        //recover time vector: 
+                        times=xNew<IssmDouble>(N);
+                        for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
+
+                        TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum,i, times,N);
+                        
+                        for(int j=0;j<elements->Size();j++){
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+
+                                for(int t=0;t<N;t++){
+                                        switch(element->ObjectEnum()){
+                                                case TriaEnum:  transientinput->AddTriaTimeInput( t,1,&element->lid,&str[t],P0Enum); break;
+                                                case PentaEnum: transientinput->AddPentaTimeInput(t,1,&element->lid,&str[t],P0Enum); break;
+                                                default: _error_("Not implemented yet");
+                                        }
+                                }
+                        }
+                        xDelete<IssmDouble>(times);
+                }
+                /*Delete data:*/
+                for(int i=0;i<nummodels;i++){
+                        IssmDouble* str=pstr[i];
+                        xDelete<IssmDouble>(str);
+                }
+                xDelete<IssmDouble*>(pstr);
+                xDelete<int>(pM);
+                xDelete<int>(pN);
+                /*}}}*/
+                /*now do the same with the dsl.sea_surface_height_change_above_geoid field:{{{ */
+                iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_surface_height_change_above_geoid");
+
+                for (int i=0;i<nummodels;i++){
+                        M=pM[i];
+                        N=pN[i];
+                        str=pstr[i];
+                
+
+                        //recover time vector: 
+                        times=xNew<IssmDouble>(N);
+                        for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
+                        
+                        TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum,i, times,N);
+        
+                        for(int j=0;j<elements->Size();j++){
+
+                                /*Get the right transient input*/
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+
+                                /*Get values and lid list*/
+                                const int   numvertices = element->GetNumberOfVertices();
+                                int        *vertexlids = xNew<int>(numvertices);
+                                int        *vertexsids = xNew<int>(numvertices);
+
+
+                                /*Recover vertices ids needed to initialize inputs*/
+                                _assert_(iomodel->elements);
+                                for(int k=0;k<numvertices;k++){
+                                        vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
+                                        vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
+                                }
+
+                                //element->GetVerticesLidList(vertexlids);
+                                //element->GetVerticesSidList(vertexsids);
+                                IssmDouble* values=xNew<IssmDouble>(numvertices);
+
+                                for(int t=0;t<N;t++){
+                                        for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
+
+                                        switch(element->ObjectEnum()){
+                                                case TriaEnum:  transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
+                                                case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
+                                                default: _error_("Not implemented yet");
+                                        }
+                                }
+                                xDelete<IssmDouble>(values);
+                                xDelete<int>(vertexlids);
+                                xDelete<int>(vertexsids);
+                        }
+                        
+                        xDelete<IssmDouble>(times);
+                }
+                
+                /*Delete data:*/
+                for(int i=0;i<nummodels;i++){
+                        IssmDouble* str=pstr[i];
+                        xDelete<IssmDouble>(str);
+                }
+                xDelete<IssmDouble*>(pstr);
+                xDelete<int>(pM);
+                xDelete<int>(pN);
+                /*}}}*/
+                /*now do the same with the dsl.sea_water_pressure_change_at_sea_floor field:{{{ */
+                iomodel->FetchData(&pstr,&pM,&pN,&nummodels,"md.dsl.sea_water_pressure_change_at_sea_floor");
+
+                for (int i=0;i<nummodels;i++){
+                        M=pM[i];
+                        N=pN[i];
+                        str=pstr[i];
+
+                        //recover time vector: 
+                        times=xNew<IssmDouble>(N);
+                        for(int t=0;t<N;t++) times[t] = str[(M-1)*N+t];
+
+                        TransientInput2* transientinput=inputs2->SetDatasetTransientInput(DslSeaWaterPressureChangeAtSeaFloor,i, times,N);
+        
+                        for(int j=0;j<elements->Size();j++){
+
+                                /*Get the right transient input*/
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(j));
+
+                                /*Get values and lid list*/
+                                const int   numvertices = element->GetNumberOfVertices();
+                                int        *vertexlids = xNew<int>(numvertices);
+                                int        *vertexsids = xNew<int>(numvertices);
+                                
+                                
+                                /*Recover vertices ids needed to initialize inputs*/
+                                _assert_(iomodel->elements);
+                                for(int k=0;k<numvertices;k++){
+                                        vertexsids[k] =reCast<int>(iomodel->elements[numvertices*element->Sid()+k]); //ids for vertices are in the elements array from Matlab
+                                        vertexlids[k]=iomodel->my_vertices_lids[vertexsids[k]-1];
+                                }
+                                //element->GetVerticesLidList(vertexlids);
+                                //element->GetVerticesSidList(vertexsids);
+                                
+                                IssmDouble* values=xNew<IssmDouble>(numvertices);
+
+                                for(int t=0;t<N;t++){
+                                        for (int k=0;k<numvertices;k++)values[k]=str[N*vertexsids[k]+t];
+
+                                        switch(element->ObjectEnum()){
+                                                case TriaEnum:  transientinput->AddTriaTimeInput( t,numvertices,vertexlids,values,P1Enum); break;
+                                                case PentaEnum: transientinput->AddPentaTimeInput(t,numvertices,vertexlids,values,P1Enum); break;
+                                                default: _error_("Not implemented yet");
+                                        }
+                                }
+                                xDelete<IssmDouble>(values);
+                                xDelete<int>(vertexlids);
+                                xDelete<int>(vertexsids);
+                        }
+                        xDelete<IssmDouble>(times);
+                }
+                
+                /*Delete data:*/
+                for(int i=0;i<nummodels;i++){
+                        IssmDouble* str=pstr[i];
+                        xDelete<IssmDouble>(str);
+                }
+                xDelete<IssmDouble*>(pstr);
+                xDelete<int>(pM);
+                xDelete<int>(pN);
+                /*}}}*/
+
+        } /*}}}*/
         else _error_("Dsl model " << dslmodel << " not implemented yet!");
                         
@@ -157 +330 @@
         parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
 
-                        
+        /*Deal with dsl multi-model ensembles: {{{*/
+        iomodel->FetchData(&dslmodel,"md.dsl.model");
+        parameters->AddObject(iomodel->CopyConstantObject("md.dsl.compute_fingerprints",DslComputeFingerprintsEnum));
+        if(dslmodel==2){
+                int modelid; 
+                int nummodels;
+
+                parameters->AddObject(iomodel->CopyConstantObject("md.dsl.modelid",DslModelidEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.dsl.nummodels",DslNummodelsEnum));
+                iomodel->FetchData(&modelid,"md.dsl.modelid");
+                iomodel->FetchData(&nummodels,"md.dsl.nummodels");
+
+                /*quick checks: */
+                if(nummodels<=0)_error_("dslmme object in  md.dsl field should contain at least 1 ensemble model!");
+                if(modelid<=0 || modelid>nummodels)_error_("modelid field in dslmme object of md.dsl field should be between 1 and the number of ensemble runs!");
+        } /*}}}*/
         /*Deal with elasticity {{{*/
         iomodel->FetchData(&elastic,"md.slr.elastic");

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

@@ -5555 +5555 @@
 void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
 
+        /*Computational flags:*/
+        int bp_compute_fingerprints= 0;
+        
+        /*some paramters first: */
+        this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
+        
+        if(!IsOceanInElement()){
+                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested 
+                 *bottom pressure fingerprints:*/
+                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
+        }
+        if(!IsIceInElement()){
+                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
+                this->SealevelriseEustaticIce(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
+        }
+
+}
+/*}}}*/
+void    Tria::SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
+
         /*diverse:*/
         int gsize;
@@ -5583 +5603 @@
         bool computeelastic= true;
         bool scaleoceanarea= false;
-
+        
         /*early return if we are not on an ice cap:*/
         if(!IsIceOnlyInElement()){
@@ -5591 +5611 @@
         }
 
-        /*early return if we are fully floating: */
+        /*early return if we are fully floating:*/
         Input2* gr_input=this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gr_input);
         if (gr_input->GetInputMax()<=0){
@@ -5599 +5619 @@
         }
 
-        /*If we are here, we are on ice that is fully grounded or half-way to floating: */
+        /*If we are here, we are on ice that is fully grounded or half-way to floating:*/
         if ((gr_input->GetInputMin())<0){
                 notfullygrounded=true; //used later on.
@@ -5677 +5697 @@
         }
 
-        /*Compute ice thickness change: */
+        /*Compute ice thickness: */
         Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
         if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
@@ -5740 +5760 @@
                         /*Add all components to the pSgi or pSgo solution vectors:*/
                         values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
+                }
+                pSgi->SetValues(gsize,indices,values,ADD_VAL);
+
+                /*free ressources:*/
+                xDelete<IssmDouble>(values);
+                xDelete<int>(indices);
+        }
+
+        /*Assign output pointer:*/
+        _assert_(!xIsNan<IssmDouble>(eustatic));
+        _assert_(!xIsInf<IssmDouble>(eustatic));
+        *peustatic=eustatic;
+        return;
+}
+/*}}}*/
+void    Tria::SealevelriseEustaticBottomPressure(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
+
+        /*diverse:*/
+        int gsize;
+        bool spherical=true;
+        IssmDouble llr_list[NUMVERTICES][3];
+        IssmDouble area;
+        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
+        IssmDouble rho;
+        IssmDouble late,longe,re;
+        IssmDouble lati,longi,ri;
+
+        /*elastic green function:*/
+        IssmDouble* G_elastic_precomputed=NULL;
+        int         M;
+
+        /*ice properties: */
+        IssmDouble rho_water,rho_earth;
+
+        /*constants:*/
+        IssmDouble constant=0;
+
+        /*Initialize eustatic component: do not skip this step :):*/
+        IssmDouble eustatic = 0.;
+
+        /*Computational flags:*/
+        bool computerigid = true;
+        bool computeelastic= true;
+        bool scaleoceanarea= false;
+        bool bp_compute_fingerprints= false;
+        
+        /*we are here to compute fingerprints originating fromn bottom pressure loads:*/
+
+        /*Inform mask: */
+        constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
+
+        /*recover material parameters: */
+        rho_water=FindParam(MaterialsRhoFreshwaterEnum);
+        rho_earth=FindParam(MaterialsEarthDensityEnum);
+
+        /*recover love numbers and computational flags: */
+        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
+        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
+        this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
+
+        /*recover elastic green function:*/
+        if(computeelastic){
+                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum));
+                _assert_(parameter);
+                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
+        }
+
+        /*how many dofs are we working with here? */
+        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
+
+        /* Where is the centroid of this element?:{{{*/
+
+        /*retrieve coordinates: */
+        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
+
+        IssmDouble minlong=400;
+        IssmDouble maxlong=-20;
+        for (int i=0;i<NUMVERTICES;i++){
+                llr_list[i][0]=(90-llr_list[i][0]);
+                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
+                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
+                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
+        }
+        if(minlong==0 && maxlong>180){
+                if (llr_list[0][1]==0)llr_list[0][1]=360;
+                if (llr_list[1][1]==0)llr_list[1][1]=360;
+                if (llr_list[2][1]==0)llr_list[2][1]=360;
+        }
+
+        // correction at the north pole
+        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
+        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
+        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
+
+        //correction at the south pole
+        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
+        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
+        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
+
+        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
+        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
+
+        late=90-late;
+        if(longe>180)longe=(longe-180)-180;
+
+        late=late/180*PI;
+        longe=longe/180*PI;
+        /*}}}*/
+
+        /*Compute area of element. For now, we dont do partially grounded elements:*/
+        area=GetAreaSpherical();
+
+        /*Compute bottom pressure change: */
+        Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
+        if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");
+
+        /*If we are fully grounded, take the average over the element: */
+        bottompressure_change_input->GetInputAverage(&I);
+
+        /*Compute eustatic compoent:*/
+        _assert_(oceanarea>0.);
+        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
+
+        /*We do not need to add the bottom pressure component to the eustatic value: */
+        eustatic += 0;
+
+        if(computeelastic | computerigid){
+                int* indices=xNew<int>(gsize);
+                IssmDouble* values=xNew<IssmDouble>(gsize);
+                IssmDouble alpha;
+                IssmDouble delPhi,delLambda;
+                for(int i=0;i<gsize;i++){
+                        indices[i]=i;
+
+                        IssmDouble G_rigid=0;  //do not remove =0!
+                        IssmDouble G_elastic=0;  //do not remove =0!
+
+                        /*Compute alpha angle between centroid and current vertex : */
+                        lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;
+
+                   delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
+                        alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));
+
+                        //Rigid earth gravitational perturbation:
+                        if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);
+
+                        //Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
+                        if(computeelastic){
+                                int index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
+                                G_elastic += G_elastic_precomputed[index];
+                        }
+
+                        /*Add all components to the pSgi or pSgo solution vectors:*/
+                        values[i]=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
                 }
                 pSgi->SetValues(gsize,indices,values,ADD_VAL);

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

@@ -167 +167 @@
                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea); 
                 void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
+                void    SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
+                void    SealevelriseEustaticBottomPressure(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
                 void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
                 void    SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);

issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp

@@ -8 +8 @@
 #include "../classes/Inputs2/TriaInput2.h"
 #include "../classes/Inputs2/TransientInput2.h"
+#include "../classes/Inputs2/DatasetInput2.h"
 #include "../shared/shared.h"
 #include "../modules/modules.h"
@@ -583 +584 @@
                 femmodel->inputs2->AddInput(tria_input_copy);
         }
+        else if(dslmodel==2){
+        
+                int modelid;
+                
+                /*Recover modelid:*/
+                femmodel->parameters->FindParam(&modelid,DslModelidEnum);
+                modelid--; //from matlab. 
+                
+                /*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/
+                DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslSeaSurfaceHeightChangeAboveGeoidEnum);
+                
+                /*Go find the modelid'th transient input:*/
+                TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
+                
+                /*Plug into DslDynamicRate input: */
+                Input2* tria_input_copy=tria_input->copy();
+                tria_input_copy->ChangeEnum(DslDynamicRateEnum);
+                femmodel->inputs2->AddInput(tria_input_copy);
+        }
         else _error_("not implemented yet");
 
@@ -603 +623 @@
                 TransientInput2* transient_input  = femmodel->inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
                 TriaInput2* tria_input=transient_input->GetTriaInput(time);
+                Input2* tria_input_copy=tria_input->copy();
+                tria_input_copy->ChangeEnum(DslStericRateEnum);
+                femmodel->inputs2->AddInput(tria_input_copy);
+        }
+        else if (dslmodel==2){
+                int modelid;
+                
+                /*Recover modelid:*/
+                femmodel->parameters->FindParam(&modelid,DslModelidEnum);
+                
+                modelid--; //from matlab. 
+                
+                /*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/
+                DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslGlobalAverageThermostericSeaLevelChangeEnum);
+                
+                /*Go find the modelid'th transient input:*/
+                TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
+                
+                /*Plug into DslStericRate input: */
                 Input2* tria_input_copy=tria_input->copy();
                 tria_input_copy->ChangeEnum(DslStericRateEnum);

issm/trunk-jpl/src/c/shared/Enum/Enum.vim

@@ -130 +130 @@
 syn keyword cConstant DomainTypeEnum
 syn keyword cConstant DslModelEnum
+syn keyword cConstant DslModelidEnum
+syn keyword cConstant DslNummodelsEnum
+syn keyword cConstant DslComputeFingerprintsEnum
 syn keyword cConstant EarthIdEnum
 syn keyword cConstant EplZigZagCounterEnum

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

@@ -124 +124 @@
         DomainTypeEnum,
         DslModelEnum,
+        DslModelidEnum,
+        DslNummodelsEnum,
+        DslComputeFingerprintsEnum,
         EarthIdEnum,
         EplZigZagCounterEnum,

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

@@ -132 +132 @@
                 case DomainTypeEnum : return "DomainType";
                 case DslModelEnum : return "DslModel";
+                case DslModelidEnum : return "DslModelid";
+                case DslNummodelsEnum : return "DslNummodels";
+                case DslComputeFingerprintsEnum : return "DslComputeFingerprints";
                 case EarthIdEnum : return "EarthId";
                 case EplZigZagCounterEnum : return "EplZigZagCounter";

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

@@ -132 +132 @@
               else if (strcmp(name,"DomainType")==0) return DomainTypeEnum;
               else if (strcmp(name,"DslModel")==0) return DslModelEnum;
+              else if (strcmp(name,"DslModelid")==0) return DslModelidEnum;
+              else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum;
+              else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum;
               else if (strcmp(name,"EarthId")==0) return EarthIdEnum;
               else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum;
-              else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum;
-              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
-              else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
+              if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum;
+              else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum;
+              else if (strcmp(name,"EsaRequestedOutputs")==0) return EsaRequestedOutputsEnum;
+              else if (strcmp(name,"EsaUElastic")==0) return EsaUElasticEnum;
               else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum;
               else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum;
@@ -257 +260 @@
               else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
               else if (strcmp(name,"MaterialsLithosphereDensity")==0) return MaterialsLithosphereDensityEnum;
-              else if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum;
-              else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum;
-              else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
+              if (strcmp(name,"MaterialsLithosphereShearModulus")==0) return MaterialsLithosphereShearModulusEnum;
+              else if (strcmp(name,"MaterialsMantleDensity")==0) return MaterialsMantleDensityEnum;
+              else if (strcmp(name,"MaterialsMantleShearModulus")==0) return MaterialsMantleShearModulusEnum;
+              else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
               else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
               else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
@@ -380 +383 @@
               else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum;
               else if (strcmp(name,"SmbRunoffgrad")==0) return SmbRunoffgradEnum;
-              else if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum;
-              else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
-              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
+              if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum;
+              else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
+              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
+              else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
               else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
               else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
@@ -503 +506 @@
               else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
               else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
-              else if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
-              else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum;
-              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum;
+              if (strcmp(name,"CalvingFluxLevelset")==0) return CalvingFluxLevelsetEnum;
+              else if (strcmp(name,"CalvingMeltingFluxLevelset")==0) return CalvingMeltingFluxLevelsetEnum;
+              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
+              else if (strcmp(name,"CrevasseDepth")==0) return CrevasseDepthEnum;
               else if (strcmp(name,"DamageD")==0) return DamageDEnum;
               else if (strcmp(name,"DamageDOld")==0) return DamageDOldEnum;
@@ -626 +629 @@
               else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum;
               else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum;
-              else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
-              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
-              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
+              if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum;
+              else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum;
+              else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum;
+              else if (strcmp(name,"MaskGroundediceLevelset")==0) return MaskGroundediceLevelsetEnum;
               else if (strcmp(name,"MaskIceLevelset")==0) return MaskIceLevelsetEnum;
               else if (strcmp(name,"MaskLandLevelset")==0) return MaskLandLevelsetEnum;
@@ -749 +752 @@
               else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
               else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
-              else if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
-              else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum;
-              else if (strcmp(name,"SmbP")==0) return SmbPEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum;
+              if (strcmp(name,"SmbNetSW")==0) return SmbNetSWEnum;
+              else if (strcmp(name,"SmbPAir")==0) return SmbPAirEnum;
+              else if (strcmp(name,"SmbP")==0) return SmbPEnum;
+              else if (strcmp(name,"SmbPddfacIce")==0) return SmbPddfacIceEnum;
               else if (strcmp(name,"SmbPddfacSnow")==0) return SmbPddfacSnowEnum;
               else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
@@ -872 +875 @@
               else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
               else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
-              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
-              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
-              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
+              if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
+              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
+              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
+              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
               else if (strcmp(name,"Outputdefinition29")==0) return Outputdefinition29Enum;
               else if (strcmp(name,"Outputdefinition2")==0) return Outputdefinition2Enum;
@@ -995 +998 @@
               else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
               else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
-              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
-              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
-              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
+              if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
+              else if (strcmp(name,"Cfdragcoeffabsgrad")==0) return CfdragcoeffabsgradEnum;
+              else if (strcmp(name,"Cfsurfacelogvel")==0) return CfsurfacelogvelEnum;
+              else if (strcmp(name,"Cfsurfacesquare")==0) return CfsurfacesquareEnum;
               else if (strcmp(name,"Channel")==0) return ChannelEnum;
               else if (strcmp(name,"ChannelArea")==0) return ChannelAreaEnum;
@@ -1118 +1121 @@
               else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
               else if (strcmp(name,"SegInput2")==0) return SegInput2Enum;
-              else if (strcmp(name,"TriaInput2")==0) return TriaInput2Enum;
-              else if (strcmp(name,"PentaInput2")==0) return PentaInput2Enum;
-              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
+              if (strcmp(name,"TriaInput2")==0) return TriaInput2Enum;
+              else if (strcmp(name,"PentaInput2")==0) return PentaInput2Enum;
+              else if (strcmp(name,"IntMatExternalResult")==0) return IntMatExternalResultEnum;
+              else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
               else if (strcmp(name,"IntParam")==0) return IntParamEnum;
               else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
@@ -1241 +1244 @@
               else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
               else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
-              else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
-              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
-              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
+              if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum;
+              else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum;
+              else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum;
+              else if (strcmp(name,"SMBgradientscomponents")==0) return SMBgradientscomponentsEnum;
               else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum;
               else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum;

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

@@ -10 +10 @@
                 sea_surface_height_change_above_geoid; %corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)
                 sea_water_pressure_change_at_sea_floor; %corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in mm/yr equivalent, not in Pa/yr!) for each ensemble
+                compute_fingerprints; %do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid
 
         end
@@ -32 +33 @@
                         self.sea_surface_height_change_above_geoid=NaN;
                         self.sea_water_pressure_change_at_sea_floor=NaN;
+                        self.compute_fingerprints=0;
 
                 end % }}}
@@ -38 +40 @@
                         %Early return
                         if ~ismember('SealevelriseAnalysis',analyses), return; end
+                        if (strcmp(solution,'TransientSolution') & md.transient.isslr == 0), return; end
                         md = checkfield(md,'fieldname','dsl.global_average_thermosteric_sea_level_change','NaN',1,'Inf',1);
                         md = checkfield(md,'fieldname','dsl.sea_surface_height_change_above_geoid','NaN',1,'Inf',1,'timeseries',1);
                         md = checkfield(md,'fieldname','dsl.sea_water_pressure_change_at_sea_floor','NaN',1,'Inf',1,'timeseries',1);
+                        md = checkfield(md,'fieldname','dsl.compute_fingerprints','NaN',1,'Inf',1,'values',[0,1]);
+                        if self.compute_fingerprints,
+                                %check geodetic flag of slr is on: 
+                                if md.slr.geodetic==0,
+                                        error('DSL checkconsistency error message: if bottom pressure fingerprints computations are requested, slr class should have geodetic flag on');
+                                end
+                        end
 
                 end % }}}
@@ -49 +59 @@
                         fielddisplay(self,'sea_surface_height_change_above_geoid','corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable rate (mm/yr)');
                         fielddisplay(self,'sea_water_pressure_change_at_sea_floor','corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable rate (in Pa/yr)');
+                        fielddisplay(self,'compute_fingerprints','%do we use the sea water pressure change to compute fingerprints and correct sea_surface_height_change_above_geoid');
 
                 end % }}}
@@ -54 +65 @@
 
                         WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
+                        WriteData(fid,prefix,'object',self,'fieldname','compute_fingerprints','format','Integer');
                         WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level_change','format','DoubleMat','mattype',1,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
                         WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_change_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts,'scale',1e-3/md.constants.yts);
@@ -62 +74 @@
                 
                         writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level_change'],self.global_average_thermosteric_sea_level_change);
+                        writejs1Darray(fid,[modelname '.dsl.compute_fingerprints'],self.compute_fingerprints);
                         writejs1Darray(fid,[modelname '.dsl.sea_surface_height_change_above_geoid'],self.sea_surface_height_change_above_geoid);
                         writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_change_at_sea_floor'],self.sea_water_pressure_change_at_sea_floor);