Changeset 27278

Timestamp:

09/10/22 08:38:49 (3 years ago)

Author:

vverjans

Message:

NEW: added possibility to use thermal forcing in basalforcingsbeckmanngoosse and ensures compatibility with thermal forcing of frontalforcingsrignot

Location:

issm/trunk-jpl/src

Files:

• c/analyses/FreeSurfaceBaseAnalysis.cpp (modified) (1 diff)
• c/analyses/LevelsetAnalysis.cpp (modified) (1 diff)
• c/analyses/MasstransportAnalysis.cpp (modified) (1 diff)
• c/analyses/StressbalanceAnalysis.cpp (modified) (1 diff)
• c/analyses/StressbalanceVerticalAnalysis.cpp (modified) (1 diff)
• c/classes/Elements/Element.cpp (modified) (5 diffs)
• c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• c/shared/Enum/Enum.vim (modified) (2 diffs)
• c/shared/Enum/EnumDefinitions.h (modified) (2 diffs)
• c/shared/Enum/EnumToStringx.cpp (modified) (2 diffs)
• c/shared/Enum/Enumjl.vim (modified) (2 diffs)
• c/shared/Enum/StringToEnumx.cpp (modified) (15 diffs)
• m/classes/basalforcingsbeckmanngoosse.m (modified) (5 diffs)

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

@@ -125 +125 @@
                         break;
                 case BeckmannGoosseFloatingMeltRateEnum:
+                        bool isthermalforcing;
+         iomodel->FindConstant(&isthermalforcing,"md.basalforcings.isthermalforcing");
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+         if(isthermalforcing==0){
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+         }
+         else{
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_thermalforcing",FrontalForcingsAndBasalforcingsThermalForcingEnum);
+         }
                         break;
                 default:

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

@@ -150 +150 @@
                 case FrontalForcingsRignotEnum:
          /*Retrieve thermal forcing only in the case of non-arma FrontalForcingsRignot*/
-         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsThermalForcingEnum);
+         iomodel->FetchDataToInput(inputs,elements,"md.frontalforcings.thermalforcing",FrontalForcingsAndBasalforcingsThermalForcingEnum);
          /* Do not break here, still retrieve basin_ID,subglacial_discharge, etc.*/
       case FrontalForcingsRignotarmaEnum:

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

@@ -215 +215 @@
                         break;
                 case BeckmannGoosseFloatingMeltRateEnum:
+                        bool isthermalforcing;
+                        iomodel->FindConstant(&isthermalforcing,"md.basalforcings.isthermalforcing");
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+                        if(isthermalforcing==0){
+                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+                        }
+                        else{
+                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_thermalforcing",FrontalForcingsAndBasalforcingsThermalForcingEnum);
+                        }
                         break;
                 case LinearFloatingMeltRatearmaEnum:

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

@@ -853 +853 @@
                                 break;
                         case BeckmannGoosseFloatingMeltRateEnum:
-                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
-                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
-                                iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+                                bool isthermalforcing;
+                                iomodel->FindConstant(&isthermalforcing,"md.basalforcings.isthermalforcing");
+           iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum);
+           if(isthermalforcing==0){
+              iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+              iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+           }
+           else{
+              iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_thermalforcing",FrontalForcingsAndBasalforcingsThermalForcingEnum);
+           }
                                 break;
                         default:

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

@@ -155 +155 @@
                         break;
                 case BeckmannGoosseFloatingMeltRateEnum:
+                        bool isthermalforcing;
+         iomodel->FindConstant(&isthermalforcing,"md.basalforcings.isthermalforcing");
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.meltrate_factor",BasalforcingsMeltrateFactorEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
-                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+         if(isthermalforcing==0){
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_salinity",BasalforcingsOceanSalinityEnum);
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_temp",BasalforcingsOceanTempEnum);
+         }
+         else{
+            iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.ocean_thermalforcing",FrontalForcingsAndBasalforcingsThermalForcingEnum);
+         }
                         break;
                 default:

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

@@ -86 +86 @@
          basinenum_type = FrontalForcingsBasinIdEnum;
          noiseenum_type = ThermalforcingARMANoiseEnum;
-         outenum_type   = FrontalForcingsThermalForcingEnum;
+         outenum_type   = FrontalForcingsAndBasalforcingsThermalForcingEnum;
          break;
                 case(BasalforcingsDeepwaterMeltingRatearmaEnum):
@@ -195 +195 @@
          basinenum_type = FrontalForcingsBasinIdEnum;
          noiseenum_type = ThermalforcingARMANoiseEnum;
-         outenum_type   = FrontalForcingsThermalForcingEnum;
+         outenum_type   = FrontalForcingsAndBasalforcingsThermalForcingEnum;
          break;
                 case(BasalforcingsDeepwaterMeltingRatearmaEnum):
@@ -2777 +2777 @@
 void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/
 
-        int numvertices      = this->GetNumberOfVertices();
-        IssmDouble T_f,ocean_heat_flux;
-        IssmDouble rho_water    = this->FindParam(MaterialsRhoSeawaterEnum);
-        IssmDouble rho_ice      = this->FindParam(MaterialsRhoIceEnum);
-        IssmDouble latentheat   = this->FindParam(MaterialsLatentheatEnum);
+        bool isthermalforcing;       
+        int numvertices                 = this->GetNumberOfVertices();
+        IssmDouble T_fr,T_forcing,ocean_heat_flux;
+        /*Material properties*/
+        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);
@@ -2789 +2791 @@
         IssmDouble oceansalinity [MAXVERTICES];
         IssmDouble oceantemp[MAXVERTICES];
+        IssmDouble oceanthermalforcing[MAXVERTICES];
         IssmDouble meltratefactor[MAXVERTICES];
 
+        /*Determine if we use temperature-and-salinity or thermal forcing*/
+        this->parameters->FindParam(&isthermalforcing,BasalforcingsIsThermalForcingEnum);
+        /*Retrieve Inputs*/
         this->GetInputListOnVertices(base,BaseEnum);
-        this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum);
-        this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum);
         this->GetInputListOnVertices(meltratefactor,BasalforcingsMeltrateFactorEnum);
+        if(isthermalforcing==false){
+                this->GetInputListOnVertices(oceansalinity,BasalforcingsOceanSalinityEnum);
+                this->GetInputListOnVertices(oceantemp,BasalforcingsOceanTempEnum);
+        }
+        else{
+                this->GetInputListOnVertices(oceanthermalforcing,FrontalForcingsAndBasalforcingsThermalForcingEnum);
+        }
 
         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[i] * 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
-        }
-
+                /*Case of temperature-and-salinity inputs*/
+                if(isthermalforcing==false){
+                        T_fr       = (0.0939-0.057*oceansalinity[i]+7.64e-4*base[i]); //freezing point [degC]
+                        T_forcing  = oceantemp[i]-T_fr; //thermal forcing [K]
+                }
+                /*Case of thermal forcing input*/
+                else{
+                        T_forcing  = oceanthermalforcing[i]; 
+                }
+                /*Compute heat flux from (Beckmann and Goosse, 2003) (>0 means heat flux from ocean to ice)*/
+                ocean_heat_flux = meltratefactor[i]*rho_water*mixed_layer_capacity*thermal_exchange_vel*T_forcing; // ocean-to-ice heat flux [W/m^2]
+                /*Save melt values [m/s]*/
+                values[i]       = ocean_heat_flux/(latentheat*rho_ice);
+        }
         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
         delete gauss;
@@ -3799 +3813 @@
         Input* bed_input = this->GetInput(BedEnum);                     _assert_(bed_input);
    Input* qsg_input = this->GetInput(FrontalForcingsSubglacialDischargeEnum);     _assert_(qsg_input);
-   Input* TF_input  = this->GetInput(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
+   Input* TF_input  = this->GetInput(FrontalForcingsAndBasalforcingsThermalForcingEnum);          _assert_(TF_input);
 
         this->FindParam(&yts, ConstantsYtsEnum);

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

@@ -253 +253 @@
                         break;
                 case BeckmannGoosseFloatingMeltRateEnum:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isthermalforcing",BasalforcingsIsThermalForcingEnum));
                         break;
                 case LinearFloatingMeltRatearmaEnum:

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

@@ -92 +92 @@
 syn keyword cConstant BasalforcingsARMAarlagcoefsEnum
 syn keyword cConstant BasalforcingsARMAmalagcoefsEnum
+syn keyword cConstant BasalforcingsIsThermalForcingEnum
 syn keyword cConstant BasalforcingsPicoAverageOverturningEnum
 syn keyword cConstant BasalforcingsPicoAverageSalinityEnum
@@ -775 +776 @@
 syn keyword cConstant FrontalForcingsBasinIdEnum
 syn keyword cConstant FrontalForcingsSubglacialDischargeEnum
-syn keyword cConstant FrontalForcingsThermalForcingEnum
+syn keyword cConstant FrontalForcingsAndBasalforcingsThermalForcingEnum
 syn keyword cConstant GeometryHydrostaticRatioEnum
 syn keyword cConstant NGiaEnum

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

@@ -86 +86 @@
         BasalforcingsARMAarlagcoefsEnum,
         BasalforcingsARMAmalagcoefsEnum,
+        BasalforcingsIsThermalForcingEnum,
         BasalforcingsPicoAverageOverturningEnum,
         BasalforcingsPicoAverageSalinityEnum,
@@ -771 +772 @@
         FrontalForcingsBasinIdEnum,
         FrontalForcingsSubglacialDischargeEnum,
-        FrontalForcingsThermalForcingEnum,
+        FrontalForcingsAndBasalforcingsThermalForcingEnum,
         GeometryHydrostaticRatioEnum,
         NGiaEnum,

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

@@ -94 +94 @@
                 case BasalforcingsARMAarlagcoefsEnum : return "BasalforcingsARMAarlagcoefs";
                 case BasalforcingsARMAmalagcoefsEnum : return "BasalforcingsARMAmalagcoefs";
+                case BasalforcingsIsThermalForcingEnum : return "BasalforcingsIsThermalForcing";
                 case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning";
                 case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity";
@@ -777 +778 @@
                 case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId";
                 case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge";
-                case FrontalForcingsThermalForcingEnum : return "FrontalForcingsThermalForcing";
+                case FrontalForcingsAndBasalforcingsThermalForcingEnum : return "FrontalForcingsAndBasalforcingsThermalForcing";
                 case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio";
                 case NGiaEnum : return "NGia";

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

@@ -85 +85 @@
 syn keyword juliaConstC BasalforcingsARMAarlagcoefsEnum
 syn keyword juliaConstC BasalforcingsARMAmalagcoefsEnum
+syn keyword juliaConstC BasalforcingsIsThermalForcingEnum
 syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum
 syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum
@@ -768 +769 @@
 syn keyword juliaConstC FrontalForcingsBasinIdEnum
 syn keyword juliaConstC FrontalForcingsSubglacialDischargeEnum
-syn keyword juliaConstC FrontalForcingsThermalForcingEnum
+syn keyword juliaConstC FrontalForcingsAndBasalforcingsThermalForcingEnum
 syn keyword juliaConstC GeometryHydrostaticRatioEnum
 syn keyword juliaConstC NGiaEnum

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

@@ -94 +94 @@
               else if (strcmp(name,"BasalforcingsARMAarlagcoefs")==0) return BasalforcingsARMAarlagcoefsEnum;
               else if (strcmp(name,"BasalforcingsARMAmalagcoefs")==0) return BasalforcingsARMAmalagcoefsEnum;
+              else if (strcmp(name,"BasalforcingsIsThermalForcing")==0) return BasalforcingsIsThermalForcingEnum;
               else if (strcmp(name,"BasalforcingsPicoAverageOverturning")==0) return BasalforcingsPicoAverageOverturningEnum;
               else if (strcmp(name,"BasalforcingsPicoAverageSalinity")==0) return BasalforcingsPicoAverageSalinityEnum;
@@ -136 +137 @@
               else if (strcmp(name,"ConstantsYts")==0) return ConstantsYtsEnum;
               else if (strcmp(name,"ControlInputSizeM")==0) return ControlInputSizeMEnum;
-              else if (strcmp(name,"ControlInputSizeN")==0) return ControlInputSizeNEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"ControlInputInterpolation")==0) return ControlInputInterpolationEnum;
+              if (strcmp(name,"ControlInputSizeN")==0) return ControlInputSizeNEnum;
+              else if (strcmp(name,"ControlInputInterpolation")==0) return ControlInputInterpolationEnum;
               else if (strcmp(name,"CumBslc")==0) return CumBslcEnum;
               else if (strcmp(name,"CumBslcIce")==0) return CumBslcIceEnum;
@@ -259 +260 @@
               else if (strcmp(name,"HydrologydcTransferFlag")==0) return HydrologydcTransferFlagEnum;
               else if (strcmp(name,"HydrologydcUnconfinedFlag")==0) return HydrologydcUnconfinedFlagEnum;
-              else if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"IcecapToEarthComm")==0) return IcecapToEarthCommEnum;
+              if (strcmp(name,"HydrologyshreveStabilization")==0) return HydrologyshreveStabilizationEnum;
+              else if (strcmp(name,"IcecapToEarthComm")==0) return IcecapToEarthCommEnum;
               else if (strcmp(name,"Index")==0) return IndexEnum;
               else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum;
@@ -382 +383 @@
               else if (strcmp(name,"QmuVariablePartitions")==0) return QmuVariablePartitionsEnum;
               else if (strcmp(name,"QmuVariablePartitionsNpart")==0) return QmuVariablePartitionsNpartEnum;
-              else if (strcmp(name,"QmuVariablePartitionsNt")==0) return QmuVariablePartitionsNtEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum;
+              if (strcmp(name,"QmuVariablePartitionsNt")==0) return QmuVariablePartitionsNtEnum;
+              else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum;
               else if (strcmp(name,"QmuResponsePartitionsNpart")==0) return QmuResponsePartitionsNpartEnum;
               else if (strcmp(name,"QmuStatistics")==0) return QmuStatisticsEnum;
@@ -505 +506 @@
               else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
               else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
-              else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
+              if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
+              else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
               else if (strcmp(name,"Smb")==0) return SmbEnum;
               else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum;
@@ -628 +629 @@
               else if (strcmp(name,"Yye")==0) return YyeEnum;
               else if (strcmp(name,"Zze")==0) return ZzeEnum;
-              else if (strcmp(name,"Areae")==0) return AreaeEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
+              if (strcmp(name,"Areae")==0) return AreaeEnum;
+              else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
               else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
               else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum;
@@ -751 +752 @@
               else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum;
               else if (strcmp(name,"Dummy")==0) return DummyEnum;
-              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"EffectivePressureSubstep")==0) return EffectivePressureSubstepEnum;
+              if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
+              else if (strcmp(name,"EffectivePressureSubstep")==0) return EffectivePressureSubstepEnum;
               else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum;
               else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum;
@@ -795 +796 @@
               else if (strcmp(name,"FrontalForcingsBasinId")==0) return FrontalForcingsBasinIdEnum;
               else if (strcmp(name,"FrontalForcingsSubglacialDischarge")==0) return FrontalForcingsSubglacialDischargeEnum;
-              else if (strcmp(name,"FrontalForcingsThermalForcing")==0) return FrontalForcingsThermalForcingEnum;
+              else if (strcmp(name,"FrontalForcingsAndBasalforcingsThermalForcing")==0) return FrontalForcingsAndBasalforcingsThermalForcingEnum;
               else if (strcmp(name,"GeometryHydrostaticRatio")==0) return GeometryHydrostaticRatioEnum;
               else if (strcmp(name,"NGia")==0) return NGiaEnum;
@@ -874 +875 @@
               else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
               else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
-              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
+              if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
+              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
               else if (strcmp(name,"Node")==0) return NodeEnum;
               else if (strcmp(name,"OmegaAbsGradient")==0) return OmegaAbsGradientEnum;
@@ -997 +998 @@
               else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
               else if (strcmp(name,"SmbD")==0) return SmbDEnum;
-              else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
+              if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
+              else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
               else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
               else if (strcmp(name,"SmbDailydsradiation")==0) return SmbDailydsradiationEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
-              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
+              if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
+              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
               else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
               else if (strcmp(name,"SurfaceSlopeX")==0) return SurfaceSlopeXEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
               else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
-              else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
+              if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
+              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
               else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
               else if (strcmp(name,"Outputdefinition74")==0) return Outputdefinition74Enum;
@@ -1366 +1367 @@
               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
               else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
-              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
+              if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
+              else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
               else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum;
               else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum;
@@ -1489 +1490 @@
               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 stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
+              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;
@@ -1612 +1613 @@
               else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
               else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
-              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
          else stage=14;
    }
    if(stage==14){
-              if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
+              if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
+              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
               else if (strcmp(name,"Sset")==0) return SsetEnum;
               else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;

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

@@ -11 +11 @@
                 ocean_temp                = 0.;
                 ocean_salinity            = NaN;
+                ocean_thermalforcing      = NaN;
+                isthermalforcing          = 0;
         end
         methods
@@ -20 +22 @@
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="ocean_temp" type="',class(self.ocean_temp),'" default="',num2str(self.ocean_temp),'">','     <section name="basalforcings" />','     <help> ocean_temp [degC] </help>','</parameter>');
                         fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="ocean_salinity" type="',class(self.ocean_salinity),'" default="',num2str(self.ocean_salinity),'">','     <section name="basalforcings" />','     <help> ocean_salinity [psu] </help>','</parameter>');
+                        fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','<parameter key ="ocean_thermalforcing" type="',class(self.ocean_thermalforcing),'" default="',num2str(self.ocean_thermalforcing),'">','     <section name="basalforcings" />','     <help> ocean_thermalforcing [K] </help>','</parameter>');
                 end % }}}
                 function self = extrude(self,md) % {{{
@@ -60 +63 @@
                 function md = checkconsistency(self,md,solution,analyses) % {{{
 
+                        if(self.isthermalforcing==0)
+                                if(any(~isnan(self.ocean_thermalforcing)))
+                                        error('basalforcingsbeckmanngoosse has isthermalforcing=0 but has specified entries for ocean_thermalforcing (should be unused)');
+                                end
+                        elseif(self.isthermalforcing==1)
+                                if(any(self.ocean_temp~=0) || any(~isnan(self.ocean_salinity)))
+                                        error('basalforcingsbeckmanngoosse has isthermalforcing=1 but has specified entries for ocean_temp or ocean_salinity (should be unused)');
+                                end
+                                if(strcmp(class(md.frontalforcings),'frontalforcingsrignot') && ~isequal(md.frontalforcings.thermalforcing,self.ocean_thermalforcing))
+                                        error('basalforcing.ocean_thermalforcing is not consistent with frontalforcings.thermalforcing')
+                                end
+                                if(strcmp(class(md.frontalforcings),'frontalforcingsrignotarma') && any(~isnan(self.ocean_thermalforcing)))
+                                        error('basalforcings has specified entries for ocean_thermalforcing, but frontalforcingsrignotarma already calculates thermalforcing internally')
+                                end
+                        end
+
+
                         if ismember('MasstransportAnalysis',analyses) & ~(solution=='TransientSolution' & md.transient.ismasstransport==0),
                                 md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
-                                md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1);
-                                md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1);  
                                 md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'universal',1,'NaN',1,'Inf',1);
+                                md = checkfield(md,'fieldname','basalforcings.isthermalforcing','values',[0 1]);
+                                if(self.isthermalforcing==0)
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1);
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1);  
+                                elseif(self.isthermalforcing==1 && ~strcmp(class(md.frontalforcings),'frontalforcingsrignotarma'))
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_thermalforcing','NaN',1,'Inf',1,'timeseries',1);    
+                                end
                         end
                         if ismember('BalancethicknessAnalysis',analyses),
                                 md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
-                                md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1);
-                                md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1);
                                 md = checkfield(md,'fieldname','basalforcings.meltrate_factor','>=',0,'universal',1,'NaN',1,'Inf',1);
+                                md = checkfield(md,'fieldname','basalforcings.isthermalforcing','values',[0 1]);
+                                if(self.isthermalforcing==0)
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_temp','NaN',1,'Inf',1,'timeseries',1);
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_salinity','NaN',1,'Inf',1,'timeseries',1);  
+                                elseif(self.isthermalforcing==1 && ~strcmp(class(md.frontalforcings),'frontalforcingsrignotarma'))
+                                        md = checkfield(md,'fieldname','basalforcings.ocean_thermalforcing','NaN',1,'Inf',1,'timeseries',1);    
+                                end
                         end
                         if ismember('ThermalAnalysis',analyses) & ~(solution=='TransientSolution' & md.transient.isthermal==0),
@@ -84 +114 @@
                         fielddisplay(self,'geothermalflux','geothermal heat flux (W/m^2)');
                         fielddisplay(self,'meltrate_factor','Melt-rate rate factor');
+                        fielddisplay(self,'isthermalforcing','whether to use ocean_temp and ocean_salinity (0) or ocean_thermalforcing (1) (default:0)');
                         fielddisplay(self,'ocean_temp','ocean temperature (degC)');
-                        fielddisplay(self,'ocean_salinity','ocean ocean_salinity (psu)');
+                        fielddisplay(self,'ocean_salinity','ocean salinity (psu)');
+                        fielddisplay(self,'ocean_thermalforcing','ocean thermalforcing (K)');
 
                 end % }}}
@@ -99 +131 @@
 
                         WriteData(fid,prefix,'name','md.basalforcings.model','data',8,'format','Integer');
-                        WriteData(fid,prefix,'data',floatingice_melting_rate,'format','DoubleMat','name','md.basalforcings.floatingice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1)
                         WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','name','md.basalforcings.groundedice_melting_rate','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1)
                         WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','name','md.basalforcings.geothermalflux','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
                         WriteData(fid,prefix,'object',self,'fieldname','meltrate_factor','format','DoubleMat','mattype',1,'name','md.basalforcings.meltrate_factor');
-                        WriteData(fid,prefix,'object',self,'fieldname','ocean_temp','format','DoubleMat','name','md.basalforcings.ocean_temp','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
-                        WriteData(fid,prefix,'object',self,'fieldname','ocean_salinity','format','DoubleMat','name','md.basalforcings.ocean_salinity','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
-
+                        WriteData(fid,prefix,'object',self,'fieldname','isthermalforcing','format','Boolean');
+                        if(self.isthermalforcing==0)
+                                WriteData(fid,prefix,'object',self,'fieldname','ocean_temp','format','DoubleMat','name','md.basalforcings.ocean_temp','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
+                                WriteData(fid,prefix,'object',self,'fieldname','ocean_salinity','format','DoubleMat','name','md.basalforcings.ocean_salinity','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
+                        elseif(self.isthermalforcing==1)
+                                WriteData(fid,prefix,'object',self,'fieldname','ocean_thermalforcing','format','DoubleMat','name','md.basalforcings.ocean_thermalforcing','mattype',1,'timeserieslength',md.mesh.numberofvertices+1);
+                        end
                 end % }}}
         end