Changeset 24145

Timestamp:

09/16/19 00:00:54 (6 years ago)

Author:

rueckamp

Message:

NEW: Added BeckmannGoose parameterization for calculating basal melt rates of floating ice

Location:

issm/trunk-jpl

Files:

• src/c/analyses/FreeSurfaceBaseAnalysis.cpp (modified) (1 diff)
• src/c/analyses/MasstransportAnalysis.cpp (modified) (1 diff)
• src/c/analyses/StressbalanceAnalysis.cpp (modified) (1 diff)
• src/c/analyses/StressbalanceVerticalAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.h (modified) (1 diff)
• src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp (modified) (2 diffs)
• src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h (modified) (1 diff)
• src/c/modules/GeothermalFluxx/GeothermalFluxx.cpp (modified) (1 diff)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (2 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (9 diffs)
• src/c/shared/io/Marshalling/IoCodeConversions.cpp (modified) (1 diff)
• src/m/classes/basalforcingsbeckmanngoosse.m (added)
• test/Archives/Archive476.arch (added)
• test/NightlyRun/test476.m (added)

issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp ¶

@@ -110 +110 @@
                 case BasalforcingsIsmip6Enum:
                         iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
+                        break;
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
                         break;
                 default:

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp ¶

@@ -198 +198 @@
                         xDelete<IssmDouble*>(array3d);
                         }
+                        break;
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
                         break;
                 default:

issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp ¶

@@ -786 +786 @@
                                 iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
                                 break;
+                        case BeckmannGoosseFloatingMeltRateEnum:
+                                iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+                                iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+                                break;
                         default:
                                 _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");

issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp ¶

@@ -138 +138 @@
                         iomodel->FetchDataToInput(elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
                         break;
-
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+                        iomodel->FetchDataToInput(elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+                        break;
                 default:
                         _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");

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

@@ -2505 +2505 @@
         xDelete<IssmDouble>(values);
 
+}/*}}}*/
+void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/
+
+        int numvertices      = this->GetNumberOfVertices();
+        IssmDouble meltratefactor,T_f,ocean_heat_flux;
+        IssmDouble rho_water    = this->FindParam(MaterialsRhoSeawaterEnum);
+        IssmDouble rho_ice      = this->FindParam(MaterialsRhoIceEnum);
+        IssmDouble latentheat   = this->FindParam(MaterialsLatentheatEnum); 
+        IssmDouble mixed_layer_capacity = this->FindParam(MaterialsMixedLayerCapacityEnum);
+        IssmDouble thermal_exchange_vel = this->FindParam(MaterialsThermalExchangeVelocityEnum);
+
+        IssmDouble* base    = xNew<IssmDouble>(numvertices);
+        IssmDouble* values  = xNew<IssmDouble>(numvertices);
+        IssmDouble* oceansalinity   = xNew<IssmDouble>(numvertices);
+        IssmDouble* oceantemp       = xNew<IssmDouble>(numvertices);
+
+        this->GetInputListOnVertices(base,BaseEnum);
+        this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum);
+        this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum);
+        parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum);
+
+        Gauss* gauss=this->NewGauss();
+        for(int i=0;i<numvertices;i++){
+                T_f=(0.0939 - 0.057 * oceansalinity[i] + 7.64e-4 * base[i]); //degC
+
+                // compute ocean_heat_flux according to beckmann_goosse2003
+                // positive, if T_oc > T_ice ==> heat flux FROM ocean TO ice
+                ocean_heat_flux = meltratefactor * rho_water * mixed_layer_capacity * thermal_exchange_vel * (oceantemp[i] - T_f); // in W/m^2
+
+                // shelfbmassflux is positive if ice is freezing on; here it is always negative:
+                // same sign as ocean_heat_flux (positive if massflux FROM ice TO ocean)
+                values[i] = ocean_heat_flux / (latentheat * rho_ice); // m s-1
+        }
+
+        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
+        xDelete<IssmDouble>(base);
+        xDelete<IssmDouble>(values);
+        xDelete<IssmDouble>(oceantemp);
+        xDelete<IssmDouble>(oceansalinity);
+        delete gauss;
 }/*}}}*/
 void       Element::MungsmtpParameterization(void){/*{{{*/

issm/trunk-jpl/src/c/classes/Elements/Element.h ¶

@@ -144 +144 @@
                 void               MigrateGroundingLine(IssmDouble* sheet_ungrounding);
                 void               MismipFloatingiceMeltingRate();
+                void               BeckmannGoosseFloatingiceMeltingRate();
                 void               MungsmtpParameterization(void);
                 ElementMatrix*     NewElementMatrix(int approximation_enum=NoneApproximationEnum);
                 ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
                 ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
-      void               PicoUpdateBoxid(int* pmax_boxid_basin); 
+                void               PicoUpdateBoxid(int* pmax_boxid_basin); 
                 void               PicoUpdateBox(int loopboxid);
                 void               PicoComputeBasalMelt(); 

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp ¶

@@ -40 +40 @@
                         if(VerboseSolution())_printf0_("   call ISMIP 6 Floating melting rate module\n");
                         FloatingiceMeltingRateIsmip6x(femmodel);
+                        break;
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        if(VerboseSolution())_printf0_("        call BeckmannGoosse Floating melting rate module\n");
+                        BeckmannGoosseFloatingiceMeltingRatex(femmodel);
                         break;
                 default:
@@ -193 +197 @@
 }
 /*}}}*/
+void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
+
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                element->BeckmannGoosseFloatingiceMeltingRate();
+        }
+}
+/*}}}*/

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.h ¶

@@ -15 +15 @@
 void MismipFloatingiceMeltingRatex(FemModel* femmodel);
 void FloatingiceMeltingRateIsmip6x(FemModel* femmodel);
+void BeckmannGoosseFloatingiceMeltingRatex(FemModel* femmodel);
 
 #endif  /* _FloatingiceMeltingRatex_H*/

issm/trunk-jpl/src/c/modules/GeothermalFluxx/GeothermalFluxx.cpp ¶

@@ -28 +28 @@
                         MantlePlumeGeothermalFluxx(femmodel);
                         break;
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        /*Nothing to be done*/
+                        break;
                 default:
                         _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp ¶

@@ -238 +238 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.num_basins",BasalforcingsIsmip6NumBasinsEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_0",BasalforcingsIsmip6Gamma0Enum));
-                   parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum));       
+                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.islocal",BasalforcingsIsmip6IsLocalEnum));  
                         iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.delta_t");
                         parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6DeltaTEnum,transparam,N));
@@ -245 +245 @@
                         parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6TfDepthsEnum,transparam,N));
                         xDelete<IssmDouble>(transparam);
+                        break;
+                case BeckmannGoosseFloatingMeltRateEnum:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum));
                         break;
                 default:

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

@@ -457 +457 @@
         BasalforcingsIsmip6TfShelfEnum,
         BasalforcingsIsmip6MeltAnomalyEnum,
+        BasalforcingsOceanSalinityEnum,
+        BasalforcingsOceanTempEnum,
         BasalforcingsPicoBasinIdEnum,
         BasalforcingsPicoBoxIdEnum,
@@ -929 +931 @@
         BasalforcingsIsmip6Enum,
         BasalforcingsPicoEnum,
+        BeckmannGoosseFloatingMeltRateEnum,     
         BedSlopeSolutionEnum,
         BoolExternalResultEnum,

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

@@ -463 +463 @@
                 case BasalforcingsIsmip6TfShelfEnum : return "BasalforcingsIsmip6TfShelf";
                 case BasalforcingsIsmip6MeltAnomalyEnum : return "BasalforcingsIsmip6MeltAnomaly";
+                case BasalforcingsOceanSalinityEnum : return "BasalforcingsOceanSalinity";
+                case BasalforcingsOceanTempEnum : return "BasalforcingsOceanTemp";
                 case BasalforcingsPicoBasinIdEnum : return "BasalforcingsPicoBasinId";
                 case BasalforcingsPicoBoxIdEnum : return "BasalforcingsPicoBoxId";
@@ -933 +935 @@
                 case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6";
                 case BasalforcingsPicoEnum : return "BasalforcingsPico";
+                case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate";
                 case BedSlopeSolutionEnum : return "BedSlopeSolution";
                 case BoolExternalResultEnum : return "BoolExternalResult";

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

@@ -472 +472 @@
               else if (strcmp(name,"BasalforcingsIsmip6TfShelf")==0) return BasalforcingsIsmip6TfShelfEnum;
               else if (strcmp(name,"BasalforcingsIsmip6MeltAnomaly")==0) return BasalforcingsIsmip6MeltAnomalyEnum;
+              else if (strcmp(name,"BasalforcingsOceanSalinity")==0) return BasalforcingsOceanSalinityEnum;
+              else if (strcmp(name,"BasalforcingsOceanTemp")==0) return BasalforcingsOceanTempEnum;
               else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
               else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
@@ -504 +506 @@
               else if (strcmp(name,"DamageF")==0) return DamageFEnum;
               else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
-              else if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
-              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
+              if (strcmp(name,"DeviatoricStresseffective")==0) return DeviatoricStresseffectiveEnum;
+              else if (strcmp(name,"DeviatoricStressxx")==0) return DeviatoricStressxxEnum;
+              else if (strcmp(name,"DeviatoricStressxy")==0) return DeviatoricStressxyEnum;
               else if (strcmp(name,"DeviatoricStressxz")==0) return DeviatoricStressxzEnum;
               else if (strcmp(name,"DeviatoricStressyy")==0) return DeviatoricStressyyEnum;
@@ -627 +629 @@
               else if (strcmp(name,"MaterialsRheologyEs")==0) return MaterialsRheologyEsEnum;
               else if (strcmp(name,"MaterialsRheologyEsbar")==0) return MaterialsRheologyEsbarEnum;
-              else if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
-              else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
+              if (strcmp(name,"MaterialsRheologyN")==0) return MaterialsRheologyNEnum;
+              else if (strcmp(name,"MeshScaleFactor")==0) return MeshScaleFactorEnum;
+              else if (strcmp(name,"MeshVertexonbase")==0) return MeshVertexonbaseEnum;
               else if (strcmp(name,"MeshVertexonboundary")==0) return MeshVertexonboundaryEnum;
               else if (strcmp(name,"MeshVertexonsurface")==0) return MeshVertexonsurfaceEnum;
@@ -750 +752 @@
               else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
               else if (strcmp(name,"SmbTemperaturesReconstructed")==0) return SmbTemperaturesReconstructedEnum;
-              else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
-              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
+              if (strcmp(name,"SmbTini")==0) return SmbTiniEnum;
+              else if (strcmp(name,"SmbTmean")==0) return SmbTmeanEnum;
+              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
               else if (strcmp(name,"SmbV")==0) return SmbVEnum;
               else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
@@ -873 +875 @@
               else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum;
               else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum;
-              else if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
-              else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
+              if (strcmp(name,"Outputdefinition52")==0) return Outputdefinition52Enum;
+              else if (strcmp(name,"Outputdefinition53")==0) return Outputdefinition53Enum;
+              else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
               else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
               else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
@@ -954 +956 @@
               else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
               else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
+              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
               else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
@@ -995 +998 @@
               else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
               else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
-              else if (strcmp(name,"Divergence")==0) return DivergenceEnum;
-              else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
-              else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
+              if (strcmp(name,"Divergence")==0) return DivergenceEnum;
+              else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum;
+              else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
+              else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
               else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
               else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
@@ -1118 +1121 @@
               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
-              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
-              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
-              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Materials")==0) return MaterialsEnum;
+              if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
+              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
+              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
               else if (strcmp(name,"Matestar")==0) return MatestarEnum;
               else if (strcmp(name,"Matice")==0) return MaticeEnum;
@@ -1241 +1244 @@
               else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
               else if (strcmp(name,"StressbalanceConvergenceNumSteps")==0) return StressbalanceConvergenceNumStepsEnum;
-              else if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
-              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
-              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
+              if (strcmp(name,"StressbalanceSIAAnalysis")==0) return StressbalanceSIAAnalysisEnum;
+              else if (strcmp(name,"StressbalanceSolution")==0) return StressbalanceSolutionEnum;
+              else if (strcmp(name,"StressbalanceVerticalAnalysis")==0) return StressbalanceVerticalAnalysisEnum;
+              else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
               else if (strcmp(name,"StringParam")==0) return StringParamEnum;

issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp ¶

@@ -198 +198 @@
                 case 6: return SpatialLinearFloatingMeltRateEnum;
                 case 7: return BasalforcingsIsmip6Enum;
+                case 8: return BeckmannGoosseFloatingMeltRateEnum;
                 default: _error_("Marshalled Basal Forcings code \""<<enum_in<<"\" not supported yet");
         }