Changeset 21550

Timestamp:

02/10/17 11:50:03 (8 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added enthalpy update in PDD scheme

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (3 diffs)

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

@@ -1993 +1993 @@
                                         desfac, s0t[iv], s0p[iv],rlaps,rlapslgm,TdiffTime,sealevTime,
                                         rho_water,rho_ice);
-        /*Get yearlytemperatures */
+                /*Get yearlytemperatures */
                 for(int month=0;month<12;month++) {
-                  yearlytemperatures[iv]=yearlytemperatures[iv]+(monthlytemperatures[iv*12+month]+273.15)*mavg; // Has to be in Kelvin
+                        yearlytemperatures[iv]=yearlytemperatures[iv]+(monthlytemperatures[iv*12+month]+273.15)*mavg; // Has to be in Kelvin
                 }
         }
@@ -2016 +2016 @@
         switch(this->ObjectEnum()){
                 case TriaEnum:  
-                        // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
-                        this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
+                        // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
+                        this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
                         this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&agd[0],P1Enum));
                         this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&accu[0],P1Enum));
@@ -2023 +2023 @@
                         break;
                 case PentaEnum:
-                        if(IsOnSurface()){
-                              GetInputListOnVertices(&s[0],TemperatureEnum);
-                              yearlytemperatures[0] = s[0];
-                              yearlytemperatures[1] = s[1];
-                              yearlytemperatures[2] = s[2];
-                              this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
-                        }
-                        this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum));
-                        this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
-                        this->InputExtrude(TemperaturePDDEnum,-1);
-                        this->InputExtrude(SmbMassBalanceEnum,-1);
-                        break;
+                        if(IsOnSurface()){
+                                /*Here, we want to change the BC of the thermal model, keep
+                                 * the temperatures as they are for the base of the penta and
+                                 * yse yearlytemperatures for the top*/
+                                GetInputListOnVertices(&s[0],TemperatureEnum);
+                                yearlytemperatures[0] = s[0];
+                                yearlytemperatures[1] = s[1];
+                                yearlytemperatures[2] = s[2];
+                                this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
+
+                                /*Convert that to enthalpy for the enthalpy model*/
+                                IssmDouble enthalpy[6];
+                                GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
+                                ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
+                                ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[3],0.,0.);
+                                ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[3],0.,0.);
+                                this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
+                        }
+                        this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&agd[0],P1Enum));
+                        this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
+                        this->InputExtrude(TemperaturePDDEnum,-1);
+                        this->InputExtrude(SmbMassBalanceEnum,-1);
+                        break;
                 case TetraEnum: 
-                        if(IsOnSurface()){
-                              GetInputListOnVertices(&s[0],TemperatureEnum);
-                              yearlytemperatures[0] = s[0];
-                              yearlytemperatures[1] = s[1];
-                              yearlytemperatures[2] = s[2];
-                              this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
-                        }
-                        this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&agd[0],P1Enum));
-                        this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
-                        this->InputExtrude(TemperaturePDDEnum,-1);
-                        this->InputExtrude(SmbMassBalanceEnum,-1);
-                        break;
+                        if(IsOnSurface()){
+                                GetInputListOnVertices(&s[0],TemperatureEnum);
+                                yearlytemperatures[0] = s[0];
+                                yearlytemperatures[1] = s[1];
+                                yearlytemperatures[2] = s[2];
+                                this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
+                        }
+                        this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&agd[0],P1Enum));
+                        this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
+                        this->InputExtrude(TemperaturePDDEnum,-1);
+                        this->InputExtrude(SmbMassBalanceEnum,-1);
+                        break;
                 default: _error_("Not implemented yet");
         }