Changeset 13119

Timestamp:

08/22/12 07:36:45 (13 years ago)

Author:

cborstad

Message:

NEW: damage inversion capability merged from branches/trunk-jpl-damage into trunk-jpl. This affects EVERYONE. You will need to initialize a new rheology parameter md.materials.rheology_Z in all your model runs. Check the updated .par files in test/Par for syntax.

Location:

issm/trunk-jpl

Files:

• src (modified) (1 prop)
• src/c/EnumDefinitions/EnumDefinitions.h (modified) (2 diffs)
• src/c/EnumDefinitions/EnumToModelField.cpp (copied) (copied from issm/branches/trunk-jpl-damage/src/c/EnumDefinitions/EnumToModelField.cpp )
• src/c/classes/objects/Elements/Penta.cpp (modified) (15 diffs)
• src/c/classes/objects/Elements/Tria.cpp (modified) (15 diffs)
• src/c/classes/objects/Elements/Tria.h (modified) (1 diff)
• src/c/classes/objects/Materials/Matice.cpp (modified) (11 diffs)
• src/c/classes/objects/Materials/Matice.h (modified) (1 diff)
• src/c/modules/EnumToStringx/EnumToStringx.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (modified) (2 diffs)
• src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (modified) (3 diffs)
• src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (modified) (1 diff)
• src/c/modules/Solverx/Solverx.cpp (modified) (1 diff)
• src/c/modules/StringToEnumx/StringToEnumx.cpp (modified) (5 diffs)
• src/m/classes/initialization.m (modified) (3 diffs)
• src/m/classes/inversion.m (modified) (1 diff)
• src/m/classes/materials.m (modified) (4 diffs)
• src/m/classes/model/model.m (modified) (3 diffs)
• src/m/contrib/hack/tres.m (modified) (1 diff)
• src/m/enum/EnumDefinitions.py (modified) (3 diffs)
• src/m/enum/MaximumNumberOfEnums.m (modified) (1 diff)
• src/m/enum/TemperatureBasalEnum.m (copied) (copied from issm/branches/trunk-jpl-damage/src/m/enum/TemperatureBasalEnum.m )
• src/m/enum/TemperatureSurfaceEnum.m (copied) (copied from issm/branches/trunk-jpl-damage/src/m/enum/TemperatureSurfaceEnum.m )
• test (modified) (1 prop)
• test/Archives/Archive270.mat (copied) (copied from issm/branches/trunk-jpl-damage/test/Archives/Archive270.mat )
• test/Archives/Archive272.mat (copied) (copied from issm/branches/trunk-jpl-damage/test/Archives/Archive272.mat )
• test/Archives/Archive274.mat (copied) (copied from issm/branches/trunk-jpl-damage/test/Archives/Archive274.mat )
• test/NightlyRun/IdToName.m (modified) (1 diff)
• test/NightlyRun/README (modified) (1 diff)
• test/NightlyRun/test270.m (copied) (copied from issm/branches/trunk-jpl-damage/test/NightlyRun/test270.m )
• test/NightlyRun/test272.m (copied) (copied from issm/branches/trunk-jpl-damage/test/NightlyRun/test272.m )
• test/NightlyRun/test274.m (copied) (copied from issm/branches/trunk-jpl-damage/test/NightlyRun/test274.m )
• test/Par/79North.par (modified) (1 diff)
• test/Par/ISMIPA.par (modified) (1 diff)
• test/Par/ISMIPB.par (modified) (1 diff)
• test/Par/ISMIPC.par (modified) (1 diff)
• test/Par/ISMIPD.par (modified) (1 diff)
• test/Par/ISMIPE.par (modified) (1 diff)
• test/Par/ISMIPF.par (modified) (1 diff)
• test/Par/Pig.par (modified) (1 diff)
• test/Par/RoundSheetEISMINT.par (modified) (1 diff)
• test/Par/RoundSheetShelf.par (modified) (1 diff)
• test/Par/RoundSheetStaticEISMINT.par (modified) (1 diff)
• test/Par/SquareEISMINT.par (modified) (1 diff)
• test/Par/SquareSheetConstrained.par (modified) (1 diff)
• test/Par/SquareSheetShelf.par (modified) (1 diff)
• test/Par/SquareShelf.par (modified) (1 diff)
• test/Par/SquareShelfConstrained.par (modified) (1 diff)
• test/Par/SquareThermal.par (modified) (1 diff)

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

@@ -100 +100 @@
         MaterialsRheologyLawEnum,
         MaterialsRheologyNEnum,
+        MaterialsRheologyZEnum,
+        MaterialsRheologyZbarEnum,
         MaterialsRhoIceEnum,
         MaterialsRhoWaterEnum,
@@ -361 +363 @@
         TemperatureOldEnum,
         TemperaturePicardEnum,
+        TemperatureSurfaceEnum,
+        TemperatureBasalEnum,
         ThicknessAbsMisfitEnum,
         TypeEnum,

issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp

@@ -1555 +1555 @@
         /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
         if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyBEnum);
+        else if (enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(MaterialsRheologyZEnum);
         else input=this->inputs->GetInput(enum_type);
         //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in penta->inputs"); why error out? if the requested input does not exist, we should still 
@@ -1667 +1668 @@
                                         break;
                                 case MaterialsRheologyBbarEnum:
+                                case MaterialsRheologyZbarEnum:
                                         /*Matice will take care of it*/ break;
                                 default:
@@ -1962 +1964 @@
 
                                 /*update input*/
-                                this->inputs->AddInput(new PentaP1Input(name,values));
+                                if (name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
+                                        matice->inputs->AddInput(new PentaP1Input(name,values));
+                                }
+                                else{
+                                        this->inputs->AddInput(new PentaP1Input(name,values));
+                                }
                                 return;
                                 break;
@@ -3230 +3237 @@
                         *presponse=this->matice->GetBbar();
                         break;
+                case MaterialsRheologyZbarEnum:
+                        *presponse=this->matice->GetZbar();
+                        break;
                 case VelEnum:
                         {
@@ -4460 +4470 @@
                 input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
         }
+        else if(enum_type==MaterialsRheologyZbarEnum){
+                if(!IsOnBed()) return;
+                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
+        }
+                
         else{
                 input=inputs->GetInput(enum_type);
@@ -4477 +4492 @@
         if(enum_type==MaterialsRheologyBbarEnum){
                 input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
+        }
+        else if(enum_type==MaterialsRheologyZbarEnum){
+                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
         }
         else{
@@ -4496 +4514 @@
         if(enum_type==MaterialsRheologyBbarEnum){
                 input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum);
+        }
+        else if(enum_type==MaterialsRheologyZbarEnum){
+                input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum);
         }
         else{
@@ -4540 +4561 @@
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
 
+        /*Depth Averaging Z*/
+        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
+
         /*Call Tria function*/
         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -4547 +4571 @@
         /*Delete B averaged*/
         this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+
+        /*Delete Z averaged*/
+        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
 
         /*clean up and return*/
@@ -5111 +5138 @@
                         input=(Input*)matice->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
                 }
+                else if(control_type[i]==MaterialsRheologyZbarEnum){
+                        if (!IsOnBed()) goto cleanup_and_return;
+                        input=(Input*)matice->inputs->GetInput(MaterialsRheologyZEnum); _assert_(input);
+                }
                 else{
                         input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
@@ -5123 +5154 @@
                 if(control_type[i]==MaterialsRheologyBbarEnum){
                         this->InputExtrude(MaterialsRheologyBEnum,MaterialsEnum);
+                }
+                else if(control_type[i]==MaterialsRheologyZbarEnum){
+                        this->InputExtrude(MaterialsRheologyZEnum,MaterialsEnum);
                 }
         }
@@ -6304 +6338 @@
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
 
+        /*Depth Averaging Z*/
+        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
+
         /*Call Tria function*/
         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -6311 +6348 @@
         /*Delete B averaged*/
         this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+
+        /*Delete Z averaged*/
+        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
 
         /*clean up and return*/
@@ -7860 +7900 @@
         this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
 
+        /*Depth Averaging Z*/
+        this->InputDepthAverageAtBase(MaterialsRheologyZEnum,MaterialsRheologyZbarEnum,MaterialsEnum);
+
         /*Call Tria function*/
         Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -7867 +7910 @@
         /*Delete B averaged*/
         this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+
+        /*Delete Z averaged*/
+        this->matice->inputs->DeleteInput(MaterialsRheologyZbarEnum);
 
         /*clean up and return*/

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

@@ -491 +491 @@
         /*Initialize Element matrix*/
         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
+
         GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);
 
@@ -1395 +1396 @@
 
         /*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
-        if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
+        if (enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum) input=this->matice->inputs->GetInput(enum_type);
         else input=this->inputs->GetInput(enum_type);
         //if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found in tria->inputs");
@@ -1501 +1502 @@
                                         break;
                                 case MaterialsRheologyBbarEnum:
+                                case MaterialsRheologyZbarEnum:
                                         /*Matice will take care of it*/ break;
                                 default:
@@ -1689 +1691 @@
 
                         /*update input*/
-                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
+                        if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum || name==MaterialsRheologyZEnum || name==MaterialsRheologyZbarEnum){
                                 matice->inputs->AddInput(new TriaP1Input(name,values));
                         }
@@ -1840 +1842 @@
                                 name==FrictionCoefficientEnum ||
                                 name==MaterialsRheologyBbarEnum ||
+                                name==MaterialsRheologyZbarEnum ||
                                 name==GradientEnum ||
                                 name==OldGradientEnum ||
@@ -2791 +2794 @@
                         *presponse=this->matice->GetBbar();
                         break;
+                case MaterialsRheologyZbarEnum:
+                        *presponse=this->matice->GetZbar();
+                        break;
                 case VelEnum:{
 
@@ -3402 +3408 @@
         for(int i=0;i<num_controls;i++){
 
-                if(control_type[i]==MaterialsRheologyBbarEnum){
+                if(control_type[i]==MaterialsRheologyBbarEnum || control_type[i]==MaterialsRheologyZbarEnum){
                         input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
                 }
@@ -3429 +3435 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3447 +3453 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3466 +3472 @@
         Input* input=NULL;
 
-        if(enum_type==MaterialsRheologyBbarEnum){
+        if(enum_type==MaterialsRheologyBbarEnum || enum_type==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(enum_type);
         }
@@ -3496 +3502 @@
                 case MaterialsRheologyBbarEnum:
                         GradjBMacAyeal(gradient,control_index);
+                        break;
+                case MaterialsRheologyZbarEnum:
+                        GradjZMacAyeal(gradient,control_index);
                         break;
                 case BalancethicknessThickeningRateEnum:
@@ -3583 +3592 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GradjZGradient{{{*/
+void  Tria::GradjZGradient(Vector* gradient,int weight_index,int control_index){
+
+        int        i,ig;
+        int        doflist1[NUMVERTICES];
+        IssmDouble     Jdet,weight;
+        IssmDouble     xyz_list[NUMVERTICES][3];
+        IssmDouble     dbasis[NDOF2][NUMVERTICES];
+        IssmDouble     dk[NDOF2]; 
+        IssmDouble     grade_g[NUMVERTICES]={0.0};
+        GaussTria  *gauss=NULL;
+
+        /*Retrieve all inputs we will be needing: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        GradientIndexing(&doflist1[0],control_index);
+        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
+        Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussTria(2);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
+                weights_input->GetInputValue(&weight,gauss,weight_index);
+
+                /*Build alpha_complement_list: */
+                rheologyz_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
+
+                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
+                for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
+        }
+        gradient->SetValues(NUMVERTICES,doflist1,grade_g,ADD_VAL);
+
+        /*Clean up and return*/
+        delete gauss;
+}
+/*}}}*/
 /*FUNCTION Tria::GradjBMacAyeal{{{*/
 void  Tria::GradjBMacAyeal(Vector* gradient,int control_index){
@@ -3624 +3673 @@
                 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
                 matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);
+
+                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
+                GetNodalFunctions(basis,gauss);
+
+                /*standard gradient dJ/dki*/
+                for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
+                                        (2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
+                                        )*Jdet*gauss->weight*basis[i];
+        }
+
+        gradient->SetValues(NUMVERTICES,doflist,grad,ADD_VAL);
+
+        /*clean-up*/
+        delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GradjZMacAyeal{{{*/
+void  Tria::GradjZMacAyeal(Vector* gradient,int control_index){
+
+        /*Intermediaries*/
+        int        i,ig;
+        int        doflist[NUMVERTICES];
+        IssmDouble     vx,vy,lambda,mu,thickness,Jdet;
+        IssmDouble     viscosity_complement;
+        IssmDouble     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dZ[NDOF2]; 
+        IssmDouble     xyz_list[NUMVERTICES][3];
+        IssmDouble     basis[3],epsilon[3];
+        IssmDouble     grad[NUMVERTICES]={0.0};
+        GaussTria *gauss = NULL;
+
+        /* Get node coordinates and dof list: */
+        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+        GradientIndexing(&doflist[0],control_index);
+
+        /*Retrieve all inputs*/
+        Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
+        Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
+        Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
+        Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
+        Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
+        Input* rheologyz_input=matice->inputs->GetInput(MaterialsRheologyZbarEnum); _assert_(rheologyz_input);
+
+        /* Start  looping on the number of gaussian points: */
+        gauss=new GaussTria(4);
+        for (ig=gauss->begin();ig<gauss->end();ig++){
+
+                gauss->GaussPoint(ig);
+
+                thickness_input->GetInputValue(&thickness,gauss);
+                rheologyz_input->GetInputDerivativeValue(&dZ[0],&xyz_list[0][0],gauss);
+                vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
+                vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
+                adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
+                adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);
+
+                this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+                matice->GetViscosityZComplement(&viscosity_complement,&epsilon[0]);
 
                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
@@ -5109 +5215 @@
 
         /*Get input (either in element or material)*/
-        if(control_enum==MaterialsRheologyBbarEnum){
+        if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
         }
@@ -5144 +5250 @@
         new_input = new TriaP1Input(control_enum,values);
 
-        if(control_enum==MaterialsRheologyBbarEnum){
+        if(control_enum==MaterialsRheologyBbarEnum || control_enum==MaterialsRheologyZbarEnum){
                 input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
         }

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

@@ -143 +143 @@
                 void   Gradj(Vector* gradient,int control_type,int control_index);
                 void   GradjBGradient(Vector* gradient,int weight_index,int control_index);
+                void   GradjZGradient(Vector* gradient,int weight_index,int control_index);
                 void   GradjBMacAyeal(Vector* gradient,int control_index);
+                void   GradjZMacAyeal(Vector* gradient,int control_index);
                 void   GradjDragMacAyeal(Vector* gradient,int control_index);
                 void   GradjDragStokes(Vector* gradient,int control_index);

issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp

@@ -172 +172 @@
 }
 /*}}}*/
+/*FUNCTION Matice::GetZ {{{*/
+IssmDouble Matice::GetZ(){
+
+        /*Output*/
+        IssmDouble Z;
+
+        inputs->GetInputAverage(&Z,MaterialsRheologyZEnum);
+        return Z;
+}
+/*}}}*/
+/*FUNCTION Matice::GetZbar {{{*/
+IssmDouble Matice::GetZbar(){
+
+        /*Output*/
+        IssmDouble Zbar;
+
+        inputs->GetInputAverage(&Zbar,MaterialsRheologyZbarEnum);
+        return Zbar;
+}
+/*}}}*/
 /*FUNCTION Matice::GetVectorFromInputs{{{*/
 void  Matice::GetVectorFromInputs(Vector* vector,int input_enum){
@@ -207 +227 @@
 void  Matice::GetViscosity2d(IssmDouble* pviscosity, IssmDouble* epsilon){
         /*From a string tensor and a material object, return viscosity, using Glen's flow law.
-                                                                                                    B
+                                                                                                   Z * B
           viscosity= -------------------------------------------------------------------
                                                   2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
@@ -226 +246 @@
         /*Intermediary: */
         IssmDouble A,e;
-        IssmDouble B,n;
+        IssmDouble Btmp,B,n,Z;
 
         /*Get B and n*/
-        B=GetBbar();
+        Btmp=GetBbar();
+        Z=GetZbar();
         n=GetN();
+        B=Z*Btmp;
 
         if (n==1){
@@ -292 +314 @@
         /*Intermediaries: */
         IssmDouble A,e;
-        IssmDouble B,n;
-
-        /*Get B and n*/
-        B=GetB();
+        IssmDouble B,n,Z;
+
+        /*Get B, Z and n*/
         n=GetN();
+        Z=GetZ();
+        B=Z*GetB();
 
         if (n==1){
@@ -362 +385 @@
         /*Intermediaries: */
         IssmDouble A,e;
-        IssmDouble B,n;
+        IssmDouble B,n,Z;
         IssmDouble eps0;
 
         /*Get B and n*/
         eps0=pow((IssmDouble)10,(IssmDouble)-27);
-        B=GetB();
         n=GetN();
+        Z=GetZ();
+        B=Z*GetB();
         
         if (n==1){
@@ -452 +476 @@
                 
                         viscosity_complement=1/(2*pow(A,e));
+                }
+        }
+        else{
+                viscosity_complement=4.5*pow((IssmDouble)10,(IssmDouble)17);
+        }
+
+        /*Checks in debugging mode*/
+        _assert_(B>0);
+        _assert_(n>0);
+        _assert_(viscosity_complement>0);
+                
+        /*Return: */
+        *pviscosity_complement=viscosity_complement;
+}
+/*}}}*/
+/*FUNCTION Matice::GetViscosityZComplement {{{*/
+void  Matice::GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* epsilon){
+        /*Return viscosity derivative for control method d(mu)/dZ: 
+         *
+         *                                                                                             B 
+         * dviscosity= -------------------------------------------------------------------
+         *                                2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n]
+         *
+         * If epsilon is NULL, it means this is the first time Gradjb is being run, and we 
+         * return mu20, initial viscosity.
+         */
+        
+        /*output: */
+        IssmDouble viscosity_complement;
+
+        /*input strain rate: */
+        IssmDouble exx,eyy,exy;
+
+        /*Intermediary value A and exponent e: */
+        IssmDouble A,e;
+        IssmDouble B,n;
+
+        /*Get B and n*/
+        B=GetBbar();
+        n=GetN();
+
+        if(epsilon){
+                exx=*(epsilon+0);
+                eyy=*(epsilon+1);
+                exy=*(epsilon+2);
+
+                /*Build viscosity: mu2=B/(2*A^e) */
+                A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy;
+                if(A==0){
+                        /*Maximum viscosity_complement for 0 shear areas: */
+                        viscosity_complement=2.25*pow((IssmDouble)10,(IssmDouble)17);
+                }
+                else{
+                        e=(n-1)/(2*n);
+                
+                        viscosity_complement=B/(2*pow(A,e));
                 }
         }
@@ -701 +781 @@
                 }
 
+                /*Get Z*/
+                if (iomodel->Data(MaterialsRheologyZEnum)) {
+                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[reCast<int,IssmDouble>(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
+                        this->inputs->AddInput(new TriaP1Input(MaterialsRheologyZbarEnum,nodeinputs));
+                }
+
                 /*Control Inputs*/
                 #ifdef _HAVE_CONTROL_
@@ -715 +801 @@
                                                 }
                                                 break;
+                                        case MaterialsRheologyZbarEnum:
+                                                if (iomodel->Data(MaterialsRheologyZEnum)){
+                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
+                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
+                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZbarEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                                }
+                                                break;
+
                                 }
                         }
@@ -741 +837 @@
                         for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyNEnum)[index];
                         this->inputs->AddInput(new PentaP1Input(MaterialsRheologyNEnum,nodeinputs));
+                }
+
+                /*Get Z*/
+                if (iomodel->Data(MaterialsRheologyZEnum)) {
+                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+i]-1)];
+                        this->inputs->AddInput(new PentaP1Input(MaterialsRheologyZEnum,nodeinputs));
                 }
 
@@ -757 +859 @@
                                                 }
                                                 break;
+                                        case MaterialsRheologyZbarEnum:
+                                                if (iomodel->Data(MaterialsRheologyZEnum)){
+                                                        _assert_(iomodel->Data(MaterialsRheologyZEnum));_assert_(iomodel->Data(InversionMinParametersEnum)); _assert_(iomodel->Data(InversionMaxParametersEnum)); 
+                                                        for(j=0;j<num_vertices;j++)nodeinputs[j]=iomodel->Data(MaterialsRheologyZEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)];
+                                                        for(j=0;j<num_vertices;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        for(j=0;j<num_vertices;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[int(iomodel->Data(MeshElementsEnum)[num_vertices*index+j]-1)*num_control_type+i];
+                                                        this->inputs->AddInput(new ControlInput(MaterialsRheologyZEnum,PentaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
+                                                }
+                                                break;
                                 }
                         }
@@ -775 +886 @@
                                 name==MaterialsRheologyBEnum ||
                                 name==MaterialsRheologyBbarEnum ||
-                                name==MaterialsRheologyNEnum
+                                name==MaterialsRheologyNEnum ||
+                                name==MaterialsRheologyZEnum ||
+                                name==MaterialsRheologyZbarEnum
                 ){
                 return true;

issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h

@@ -58 +58 @@
                 /*}}}*/
                 /*Matice Numerics: {{{*/
-                void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
-                void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
-                void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
-                void       GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
-                void       GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
-                void       GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
-                void       GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
+                void   SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin);
+                void   GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
+                void   GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
+                void   GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon);
+                void   GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
+                void   GetViscosityZComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon);
+                void   GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
+                void   GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon);
                 IssmDouble GetA();
                 IssmDouble GetB();
                 IssmDouble GetBbar();
                 IssmDouble GetN();
-                bool       IsInput(int name);
+                IssmDouble GetZ();
+                IssmDouble GetZbar();
+                bool   IsInput(int name);
                 /*}}}*/
 };

issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

@@ -105 +105 @@
                 case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
                 case MaterialsRheologyNEnum : return "MaterialsRheologyN";
+                case MaterialsRheologyZEnum : return "MaterialsRheologyZ";
+                case MaterialsRheologyZbarEnum : return "MaterialsRheologyZbar";
                 case MaterialsRhoIceEnum : return "MaterialsRhoIce";
                 case MaterialsRhoWaterEnum : return "MaterialsRhoWater";
@@ -352 +354 @@
                 case TemperatureOldEnum : return "TemperatureOld";
                 case TemperaturePicardEnum : return "TemperaturePicard";
+                case TemperatureSurfaceEnum : return "TemperatureSurface";
+                case TemperatureBasalEnum : return "TemperatureBasal";
                 case ThicknessAbsMisfitEnum : return "ThicknessAbsMisfit";
                 case TypeEnum : return "Type";

issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

@@ -48 +48 @@
                         case FrictionCoefficientEnum: iomodel->FetchData(1,FrictionCoefficientEnum); break;
                         case MaterialsRheologyBbarEnum:    iomodel->FetchData(1,MaterialsRheologyBEnum); break;
+                        case MaterialsRheologyZbarEnum:    iomodel->FetchData(1,MaterialsRheologyZEnum); break;
                         default: _error_("Control " << EnumToStringx(reCast<int,IssmDouble>(iomodel->Data(InversionControlParametersEnum)[i])) << " not implemented yet");
                 }
@@ -66 +67 @@
         
         /*Free data: */
-        iomodel->DeleteData(1+4+5,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum);
+        iomodel->DeleteData(1+4+6,MeshElementsEnum,InversionControlParametersEnum,InversionCostFunctionsCoefficientsEnum,InversionMinParametersEnum,InversionMaxParametersEnum,BalancethicknessThickeningRateEnum,VxEnum,VyEnum,FrictionCoefficientEnum,MaterialsRheologyBEnum,MaterialsRheologyZEnum);
 }

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

@@ -45 +45 @@
         
         /*Fetch data needed: */
-        iomodel->FetchData(4,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum);
+        iomodel->FetchData(5,MeshElementsEnum,MeshElementconnectivityEnum,MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum);
         #ifdef _HAVE_3D_
         if(dim==3)iomodel->FetchData(2,MeshUpperelementsEnum,MeshLowerelementsEnum);
@@ -67 +67 @@
         
         /*Free data: */
-        iomodel->DeleteData(9,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum,
-                                MaterialsRheologyBEnum,MaterialsRheologyNEnum,InversionControlParametersEnum,InversionMinParametersEnum,InversionMaxParametersEnum);
+        iomodel->DeleteData(10,MeshElementsEnum,MeshElementconnectivityEnum,MeshUpperelementsEnum,MeshLowerelementsEnum,
+                                MaterialsRheologyBEnum,MaterialsRheologyNEnum,MaterialsRheologyZEnum,InversionControlParametersEnum,InversionMinParametersEnum,
+                                InversionMaxParametersEnum);
 
-        /*Add new constrant material property tgo materials, at the end: */
+        /*Add new constant material property to materials, at the end: */
         materials->AddObject(new Matpar(numberofelements+1,iomodel));//put it at the end of the materials
         
@@ -81 +82 @@
         for (i=0;i<numberofvertices;i++){
 
-                /*vertices and nodes (same number, as we are running continuous galerkin formulation: */
+                /*vertices and nodes (same number, as we are running continuous galerkin formulation): */
                 if(iomodel->my_vertices[i]){
                         

issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp

@@ -61 +61 @@
         iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
         iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
+        iomodel->FetchDataToInput(elements,MaterialsRheologyZEnum);
         iomodel->FetchDataToInput(elements,VxEnum);
         iomodel->FetchDataToInput(elements,VyEnum);

issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp

@@ -47 +47 @@
         }
 
-        /*Assign output pointers: */
+        /*Assign output pointers:*/
         *puf=uf;
 }

issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

@@ -106 +106 @@
               else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
               else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
+              else if (strcmp(name,"MaterialsRheologyZ")==0) return MaterialsRheologyZEnum;
+              else if (strcmp(name,"MaterialsRheologyZbar")==0) return MaterialsRheologyZbarEnum;
               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
               else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
@@ -136 +138 @@
               else if (strcmp(name,"PrognosticHydrostaticAdjustment")==0) return PrognosticHydrostaticAdjustmentEnum;
               else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
-              else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
-              else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
+              if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
+              else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
+              else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
               else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
               else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
@@ -259 +261 @@
               else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
               else if (strcmp(name,"Contour")==0) return ContourEnum;
-              else if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
-              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
+              if (strcmp(name,"ControlInput")==0) return ControlInputEnum;
+              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
+              else if (strcmp(name,"DofIndexing")==0) return DofIndexingEnum;
               else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
               else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
@@ -359 +361 @@
               else if (strcmp(name,"TemperatureOld")==0) return TemperatureOldEnum;
               else if (strcmp(name,"TemperaturePicard")==0) return TemperaturePicardEnum;
+              else if (strcmp(name,"TemperatureSurface")==0) return TemperatureSurfaceEnum;
+              else if (strcmp(name,"TemperatureBasal")==0) return TemperatureBasalEnum;
               else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
               else if (strcmp(name,"Type")==0) return TypeEnum;
@@ -380 +384 @@
               else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum;
               else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum;
-              else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
+         else stage=4;
+   }
+   if(stage==4){
+              if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
               else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
               else if (strcmp(name,"StepResponses")==0) return StepResponsesEnum;
               else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum;
-         else stage=4;
-   }
-   if(stage==4){
-              if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
+              else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum;
               else if (strcmp(name,"DragCoefficientAbsGradient")==0) return DragCoefficientAbsGradientEnum;
               else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;

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

@@ -12 +12 @@
                 pressure      = NaN;
                 temperature   = NaN;
+                surfacetemp   = NaN;
+                basaltemp     = NaN;
                 watercolumn   = NaN;
                 waterfraction = NaN;
@@ -68 +70 @@
                         fielddisplay(obj,'pressure','pressure field');
                         fielddisplay(obj,'temperature','temperature in Kelvins');
+                        fielddisplay(obj,'surfacetemp','surface temperature in Kelvins');
+                        fielddisplay(obj,'basaltemp','basal temperature in Kelvins');
                         fielddisplay(obj,'watercolumn','thickness of subglacial water');
                         fielddisplay(obj,'waterfraction','fraction of water in the ice');
@@ -73 +77 @@
                 end % }}}
                 function marshall(obj,fid) % {{{
-                        WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum());
-                        WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum());
-                        WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum());
-                        WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum());
-                        WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum());
-                        WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum());
-                        WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum());
+                        WriteData(fid,'data',obj.vx,'format','DoubleMat','mattype',1,'enum',VxEnum);
+                        WriteData(fid,'data',obj.vy,'format','DoubleMat','mattype',1,'enum',VyEnum);
+                        WriteData(fid,'data',obj.vz,'format','DoubleMat','mattype',1,'enum',VzEnum);
+                        WriteData(fid,'data',obj.pressure,'format','DoubleMat','mattype',1,'enum',PressureEnum);
+                        WriteData(fid,'data',obj.temperature,'format','DoubleMat','mattype',1,'enum',TemperatureEnum);
+                        WriteData(fid,'data',obj.surfacetemp,'format','DoubleMat','mattype',1,'enum',TemperatureSurfaceEnum);
+                        WriteData(fid,'data',obj.basaltemp,'format','DoubleMat','mattype',1,'enum',TemperatureBasalEnum);
+                        WriteData(fid,'data',obj.watercolumn,'format','DoubleMat','mattype',1,'enum',WatercolumnEnum);
+                        WriteData(fid,'data',obj.waterfraction,'format','DoubleMat','mattype',1,'enum',WaterfractionEnum);
                 end % }}}
         end

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

@@ -83 +83 @@
                         num_costfunc=size(md.inversion.cost_functions,2);
 
-                        md = checkfield(md,'inversion.iscontrol','values',[0 1]);
-                        md = checkfield(md,'inversion.tao','values',[0 1]);
-                        md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
-                        md = checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'});
-                        md = checkfield(md,'inversion.nsteps','numel',[1],'>=',1);
-                        md = checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
-                        md = checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]);
-                        md = checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:505]);
-                        md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
-                        md = checkfield(md,'inversion.gradient_only','values',[0 1]);
-                        md = checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]);
-                        md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
-                        md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
+                        checkfield(md,'inversion.iscontrol','values',[0 1]);
+                        checkfield(md,'inversion.tao','values',[0 1]);
+                        checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
+                        checkfield(md,'inversion.control_parameters','cell',1,'values',{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy'});
+                        checkfield(md,'inversion.nsteps','numel',1,'>=',1);
+                        checkfield(md,'inversion.maxiter_per_step','size',[md.inversion.nsteps 1],'>=',0);
+                        checkfield(md,'inversion.step_threshold','size',[md.inversion.nsteps 1]);
+                        checkfield(md,'inversion.cost_functions','size',[md.inversion.nsteps num_costfunc],'values',[101:105 201 501:503]);
+                        checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
+                        checkfield(md,'inversion.gradient_only','values',[0 1]);
+                        checkfield(md,'inversion.gradient_scaling','size',[md.inversion.nsteps num_controls]);
+                        checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
+                        checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
 
                         if solution==BalancethicknessSolutionEnum()

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

@@ -19 +19 @@
                 rheology_B   = NaN;
                 rheology_n   = NaN;
+                rheology_Z   = NaN;
                 rheology_law = '';
         end
@@ -76 +77 @@
                         md = checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
                         md = checkfield(md,'materials.rheology_n','>',0,'size',[md.mesh.numberofelements 1]);
+                        md = checkfield(md,'materials.rheology_Z','>',0,'size',[md.mesh.numberofvertices 1]);
                         md = checkfield(md,'materials.rheology_law','values',{'None' 'Paterson' 'Arrhenius'});
                 end % }}}
@@ -94 +96 @@
                         fielddisplay(obj,'rheology_B','flow law parameter [Pa/s^(1/n)]');
                         fielddisplay(obj,'rheology_n','Glen''s flow law exponent');
+                        fielddisplay(obj,'rheology_Z','rheology multiplier');
                         fielddisplay(obj,'rheology_law','law for the temperature dependance of the rheology: ''None'', ''Paterson'' or ''Arrhenius''');
                 end % }}}
@@ -110 +113 @@
                         WriteData(fid,'object',obj,'fieldname','rheology_B','format','DoubleMat','mattype',1);
                         WriteData(fid,'object',obj,'fieldname','rheology_n','format','DoubleMat','mattype',2);
+                        WriteData(fid,'object',obj,'fieldname','rheology_Z','format','DoubleMat','mattype',1);
                         WriteData(fid,'data',StringToEnum(obj.rheology_law),'enum',MaterialsRheologyLawEnum(),'format','Integer');
                 end % }}}

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

@@ -168 +168 @@
                          md.materials.rheology_B=DepthAverage(md,md.materials.rheology_B);
                          md.materials.rheology_n=project2d(md,md.materials.rheology_n,1);
+                         md.materials.rheology_Z=DepthAverage(md,md.materials.rheology_Z);
 
                          %special for thermal modeling:
@@ -716 +717 @@
                          md.materials.rheology_B=project3d(md,'vector',md.materials.rheology_B,'type','node');
                          md.materials.rheology_n=project3d(md,'vector',md.materials.rheology_n,'type','element');
+                         md.materials.rheology_Z=project3d(md,'vector',md.materials.rheology_Z,'type','node');
 
                          %parameters
@@ -829 +831 @@
                          if isfield(structmd,'B'), md.materials.rheology_B=structmd.B; end
                          if isfield(structmd,'n'), md.materials.rheology_n=structmd.n; end
+                         if isfield(structmd,'Z'), md.materials.rheology_Z=structmd.Z; end
                          if isfield(structmd,'rheology_B'), md.materials.rheology_B=structmd.rheology_B; end
                          if isfield(structmd,'rheology_n'), md.materials.rheology_n=structmd.rheology_n; end
+                         if isfield(structmd,'rheology_Z'), md.materials.rheology_Z=structmd.rheology_Z; end
                          if isfield(structmd,'elementoniceshelf'), md.mask.elementonfloatingice=structmd.elementoniceshelf; end
                          if isfield(structmd,'elementonicesheet'), md.mask.elementongroundedice=structmd.elementonicesheet; end

issm/trunk-jpl/src/m/contrib/hack/tres.m

@@ -105 +105 @@
                 case MaterialsRheologyBEnum(), string='rheology_B'; return
                 case MaterialsRheologyBbarEnum(), string='rheology_B'; return
+                case MaterialsRheologyZEnum(), string='rheology_Z'; return
+                case MaterialsRheologyZbarEnum(), string='rheology_Z'; return
                 case BalancethicknessThickeningRateEnum(), string='dhdt'; return
                 case VxEnum(), string='vx'; return

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -889 +889 @@
         return StringToEnum('MaterialsRheologyN')[0]
 
+def MaterialsRheologyZEnum():
+        """
+        MATERIALSRHEOLOGYZENUM - Enum of MaterialsRheologyZ
+
+           Usage:
+              macro=MaterialsRheologyZEnum()
+        """
+
+        return StringToEnum('MaterialsRheologyZ')[0]
+
+def MaterialsRheologyZbarEnum():
+        """
+        MATERIALSRHEOLOGYZBARENUM - Enum of MaterialsRheologyZbar
+
+           Usage:
+              macro=MaterialsRheologyZbarEnum()
+        """
+
+        return StringToEnum('MaterialsRheologyZbar')[0]
+
 def MaterialsRhoIceEnum():
         """
@@ -3359 +3379 @@
         return StringToEnum('TemperaturePicard')[0]
 
+def TemperatureSurfaceEnum():
+        """
+        TEMPERATURESURFACEENUM - Enum of TemperatureSurface
+
+           Usage:
+              macro=TemperatureSurfaceEnum()
+        """
+
+        return StringToEnum('TemperatureSurface')[0]
+
+def TemperatureBasalEnum():
+        """
+        TEMPERATUREBASALENUM - Enum of TemperatureBasal
+
+           Usage:
+              macro=TemperatureBasalEnum()
+        """
+
+        return StringToEnum('TemperatureBasal')[0]
+
 def ThicknessAbsMisfitEnum():
         """
@@ -4607 +4647 @@
         """
 
-        return 459
-
+        return 463
+

issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m

@@ -9 +9 @@
 %      macro=MaximumNumberOfEnums()
 
-macro=459;
+macro=463;

issm/trunk-jpl/test/NightlyRun/IdToName.m

@@ -65 +65 @@
         case 236, name='SquareShelfTranIspddIsdeltaM2d';
         case 237, name='SquareShelfTranIspddIsdeltaM3d';
+        case 270, name='SquareShelfDiagM2dDamage';
+        case 272, name='SquareShelfCMZM2dDamage';
+        case 274, name='SquareShelfDiagM2dDamageRift';
         case 301, name='SquareSheetConstrainedDiagM2d';
         case 302, name='SquareSheetConstrainedDiagH2d';

issm/trunk-jpl/test/NightlyRun/README

@@ -31 +31 @@
 6: 79North
 Add the id and testname in IdToName.m (incresing order)
-We try not to create to many .par and .exp files, so try to use the existing ones as much as possible.
+We try not to create too many .par and .exp files, so try to use the existing ones as much as possible.
 To modify some characteristics, do it in the testxxx.m file.
 The testxxx_nightly.m is used to define the parameters you want to check in the nightlyruns.

issm/trunk-jpl/test/Par/79North.par

@@ -13 +13 @@
 md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_B=paterson(md.initialization.temperature);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
 

issm/trunk-jpl/test/Par/ISMIPA.par

@@ -17 +17 @@
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      boundary conditions for diagnostic model');

issm/trunk-jpl/test/Par/ISMIPB.par

@@ -17 +17 @@
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      boundary conditions for diagnostic model');

issm/trunk-jpl/test/Par/ISMIPC.par

@@ -18 +18 @@
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      boundary conditions for diagnostic model: ');

issm/trunk-jpl/test/Par/ISMIPD.par

@@ -17 +17 @@
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      boundary conditions for diagnostic model: ');

issm/trunk-jpl/test/Par/ISMIPE.par

@@ -26 +26 @@
 md.materials.rheology_B=6.8067*10^7*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      boundary conditions for diagnostic model: ');

issm/trunk-jpl/test/Par/ISMIPF.par

@@ -16 +16 @@
 md.materials.rheology_B=1.4734*10^14*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=1*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 md.materials.rheology_law='None';
 

issm/trunk-jpl/test/Par/Pig.par

@@ -17 +17 @@
 md.initialization.temperature=(273-20)*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_B=paterson(md.initialization.temperature);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
 md.initialization.temperature=md.initialization.temperature;

issm/trunk-jpl/test/Par/RoundSheetEISMINT.par

@@ -20 +20 @@
 md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution 
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      creating surface mass balance');

issm/trunk-jpl/test/Par/RoundSheetShelf.par

@@ -60 +60 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 %Surface mass balance and basal melting

issm/trunk-jpl/test/Par/RoundSheetStaticEISMINT.par

@@ -25 +25 @@
 md.materials.rheology_B=6.81*10^(7)*ones(md.mesh.numberofvertices,1); %to have the same B as the analytical solution 
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      creating surface mass balance');

issm/trunk-jpl/test/Par/SquareEISMINT.par

@@ -27 +27 @@
 md.materials.rheology_B=1.7687*10^8*ones(md.mesh.numberofvertices,1);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      creating surface mass balance');

issm/trunk-jpl/test/Par/SquareSheetConstrained.par

@@ -22 +22 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 %Friction

issm/trunk-jpl/test/Par/SquareSheetShelf.par

@@ -25 +25 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 %Accumulation and melting

issm/trunk-jpl/test/Par/SquareShelf.par

@@ -22 +22 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 %Friction

issm/trunk-jpl/test/Par/SquareShelfConstrained.par

@@ -22 +22 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 %Surface mass balance and basal melting

issm/trunk-jpl/test/Par/SquareThermal.par

@@ -28 +28 @@
 md.materials.rheology_B=paterson(md.initialization.temperature);
 md.materials.rheology_n=3*ones(md.mesh.numberofelements,1);
+md.materials.rheology_Z=ones(md.mesh.numberofvertices,1);
 
 disp('      creating surface mass balance');