Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 27278) @@ -125,7 +125,14 @@ 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: Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 27278) @@ -150,5 +150,5 @@ 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: Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 27278) @@ -215,7 +215,14 @@ 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: Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 27278) @@ -853,7 +853,14 @@ 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: Index: /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp (revision 27278) @@ -155,7 +155,14 @@ 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: Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 27278) @@ -86,5 +86,5 @@ basinenum_type = FrontalForcingsBasinIdEnum; noiseenum_type = ThermalforcingARMANoiseEnum; - outenum_type = FrontalForcingsThermalForcingEnum; + outenum_type = FrontalForcingsAndBasalforcingsThermalForcingEnum; break; case(BasalforcingsDeepwaterMeltingRatearmaEnum): @@ -195,5 +195,5 @@ basinenum_type = FrontalForcingsBasinIdEnum; noiseenum_type = ThermalforcingARMANoiseEnum; - outenum_type = FrontalForcingsThermalForcingEnum; + outenum_type = FrontalForcingsAndBasalforcingsThermalForcingEnum; break; case(BasalforcingsDeepwaterMeltingRatearmaEnum): @@ -2777,9 +2777,11 @@ 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,24 +2791,36 @@ 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 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,5 +3813,5 @@ 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); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 27278) @@ -253,4 +253,5 @@ break; case BeckmannGoosseFloatingMeltRateEnum: + parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.isthermalforcing",BasalforcingsIsThermalForcingEnum)); break; case LinearFloatingMeltRatearmaEnum: Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27277) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 27278) @@ -92,4 +92,5 @@ syn keyword cConstant BasalforcingsARMAarlagcoefsEnum syn keyword cConstant BasalforcingsARMAmalagcoefsEnum +syn keyword cConstant BasalforcingsIsThermalForcingEnum syn keyword cConstant BasalforcingsPicoAverageOverturningEnum syn keyword cConstant BasalforcingsPicoAverageSalinityEnum @@ -775,5 +776,5 @@ syn keyword cConstant FrontalForcingsBasinIdEnum syn keyword cConstant FrontalForcingsSubglacialDischargeEnum -syn keyword cConstant FrontalForcingsThermalForcingEnum +syn keyword cConstant FrontalForcingsAndBasalforcingsThermalForcingEnum syn keyword cConstant GeometryHydrostaticRatioEnum syn keyword cConstant NGiaEnum Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27277) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 27278) @@ -86,4 +86,5 @@ BasalforcingsARMAarlagcoefsEnum, BasalforcingsARMAmalagcoefsEnum, + BasalforcingsIsThermalForcingEnum, BasalforcingsPicoAverageOverturningEnum, BasalforcingsPicoAverageSalinityEnum, @@ -771,5 +772,5 @@ FrontalForcingsBasinIdEnum, FrontalForcingsSubglacialDischargeEnum, - FrontalForcingsThermalForcingEnum, + FrontalForcingsAndBasalforcingsThermalForcingEnum, GeometryHydrostaticRatioEnum, NGiaEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 27278) @@ -94,4 +94,5 @@ case BasalforcingsARMAarlagcoefsEnum : return "BasalforcingsARMAarlagcoefs"; case BasalforcingsARMAmalagcoefsEnum : return "BasalforcingsARMAmalagcoefs"; + case BasalforcingsIsThermalForcingEnum : return "BasalforcingsIsThermalForcing"; case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning"; case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity"; @@ -777,5 +778,5 @@ case FrontalForcingsBasinIdEnum : return "FrontalForcingsBasinId"; case FrontalForcingsSubglacialDischargeEnum : return "FrontalForcingsSubglacialDischarge"; - case FrontalForcingsThermalForcingEnum : return "FrontalForcingsThermalForcing"; + case FrontalForcingsAndBasalforcingsThermalForcingEnum : return "FrontalForcingsAndBasalforcingsThermalForcing"; case GeometryHydrostaticRatioEnum : return "GeometryHydrostaticRatio"; case NGiaEnum : return "NGia"; Index: /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27277) +++ /issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim (revision 27278) @@ -85,4 +85,5 @@ syn keyword juliaConstC BasalforcingsARMAarlagcoefsEnum syn keyword juliaConstC BasalforcingsARMAmalagcoefsEnum +syn keyword juliaConstC BasalforcingsIsThermalForcingEnum syn keyword juliaConstC BasalforcingsPicoAverageOverturningEnum syn keyword juliaConstC BasalforcingsPicoAverageSalinityEnum @@ -768,5 +769,5 @@ syn keyword juliaConstC FrontalForcingsBasinIdEnum syn keyword juliaConstC FrontalForcingsSubglacialDischargeEnum -syn keyword juliaConstC FrontalForcingsThermalForcingEnum +syn keyword juliaConstC FrontalForcingsAndBasalforcingsThermalForcingEnum syn keyword juliaConstC GeometryHydrostaticRatioEnum syn keyword juliaConstC NGiaEnum Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27277) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 27278) @@ -94,4 +94,5 @@ 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,9 +137,9 @@ 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,9 +260,9 @@ 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,9 +383,9 @@ 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,9 +506,9 @@ 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,9 +629,9 @@ 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,9 +752,9 @@ 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,5 +796,5 @@ 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,9 +875,9 @@ 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,9 +998,9 @@ 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,9 +1121,9 @@ 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,9 +1244,9 @@ 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,9 +1367,9 @@ 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,9 +1490,9 @@ 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,9 +1613,9 @@ 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; Index: /issm/trunk-jpl/src/m/classes/basalforcingsbeckmanngoosse.m =================================================================== --- /issm/trunk-jpl/src/m/classes/basalforcingsbeckmanngoosse.m (revision 27277) +++ /issm/trunk-jpl/src/m/classes/basalforcingsbeckmanngoosse.m (revision 27278) @@ -11,4 +11,6 @@ ocean_temp = 0.; ocean_salinity = NaN; + ocean_thermalforcing = NaN; + isthermalforcing = 0; end methods @@ -20,4 +22,5 @@ fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean_temp [degC] ',''); fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean_salinity [psu] ',''); + fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','','
',' ocean_thermalforcing [K] ',''); end % }}} function self = extrude(self,md) % {{{ @@ -60,15 +63,42 @@ 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,6 +114,8 @@ 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,11 +131,14 @@ 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