Changeset 18667

Timestamp:

10/23/14 07:11:42 (10 years ago)

Author:

jbondzio

Message:

CHG: For computation of basal melting rate, use enthalpypicard for steady state, if not converged yet. Update Enthalpy and water column only in case of transient runs.

File:

• issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (modified) (5 diffs)

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

@@ -933 +933 @@
 
         /* Intermediaries */
+        bool                    converged;
         const int   dim=3;
         int         i,is,state;
         int                     vertexdown,vertexup,numvertices,numsegments;
-        const int       enthalpy_enum=EnthalpyEnum;
+        int                     enthalpy_enum;
         IssmDouble  vec_heatflux[dim],normal_base[dim],d1enthalpy[dim],d1pressure[dim];
         IssmDouble  basalfriction,alpha2,geothermalflux,heatflux;
@@ -949 +950 @@
         element->GetVerticesCoordinates(&xyz_list);
         element->GetVerticesCoordinatesBase(&xyz_list_base);
+        element->GetInputValue(&converged,ConvergedEnum);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        element->FindParam(&yts, ConstantsYtsEnum);
+
+        if(dt==0. && !converged) enthalpy_enum=EnthalpyPicardEnum;
+        else enthalpy_enum=EnthalpyEnum;
+
         IssmDouble latentheat = element->GetMaterialParameter(MaterialsLatentheatEnum);
         IssmDouble rho_ice    = element->GetMaterialParameter(MaterialsRhoIceEnum);
         IssmDouble rho_water  = element->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
         IssmDouble beta          = element->GetMaterialParameter(MaterialsBetaEnum);
-        IssmDouble kappa                 = EnthalpyDiffusionParameterVolume(element,EnthalpyEnum);     _assert_(kappa>=0.);
+        IssmDouble kappa                 = EnthalpyDiffusionParameterVolume(element,enthalpy_enum);     _assert_(kappa>=0.);
         IssmDouble kappa_mix;
 
@@ -1002 +1010 @@
                         case 3:
                                 // temperate, thin melting base: set spc
-                                // enthalpies[vertexdown]=element->PureIceEnthalpy(pressures[vertexdown]); 
+                                enthalpies[vertexdown]=element->PureIceEnthalpy(pressures[vertexdown]); 
                                 enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],xyz_list,gauss);
                                 for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i];
@@ -1033 +1041 @@
 
         /******** UPDATE MELTINGRATES AND WATERCOLUMN **************//*{{{*/
-        element->FindParam(&dt,TimesteppingTimeStepEnum);
-        element->FindParam(&yts, ConstantsYtsEnum);
         for(is=0;is<numsegments;is++){
                 vertexdown = pairindices[is*2+0];
@@ -1062 +1068 @@
 
         /*feed updated variables back into model*/
-        element->AddInput(EnthalpyEnum,enthalpies,P1Enum); //TODO: distinguis for steadystate and transient run
-        element->AddInput(WatercolumnEnum,watercolumns,P1Enum);
+        if(dt!=0.){
+                element->AddInput(enthalpy_enum,enthalpies,P1Enum); //TODO: distinguis for steadystate and transient run
+                element->AddInput(WatercolumnEnum,watercolumns,P1Enum);
+        }
         element->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,P1Enum);