Index: /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12676) +++ /issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h (revision 12677) @@ -11,5 +11,5 @@ enum definitions{ - /*Model fields {{{*/ + /*Model fields {{{1*/ AutodiffForwardEnum, AutodiffIsautodiffEnum, @@ -162,5 +162,13 @@ SurfaceforcingsMassBalanceEnum, SurfaceforcingsIspddEnum, + SurfaceforcingsIssmbgradientsEnum, SurfaceforcingsMonthlytemperaturesEnum, + SurfaceforcingsHcEnum, + SurfaceforcingsSmbPosMaxEnum, + SurfaceforcingsSmbPosMinEnum, + SurfaceforcingsAPosEnum, + SurfaceforcingsBPosEnum, + SurfaceforcingsANegEnum, + SurfaceforcingsBNegEnum, ThermalMaxiterEnum, ThermalPenaltyFactorEnum, @@ -183,5 +191,5 @@ TransientRequestedOutputsEnum, /*}}}*/ - /*Solutions and Analyses{{{ */ + /*Solutions and Analyses{{{1 */ SolutionTypeEnum, AnalysisTypeEnum, @@ -220,5 +228,5 @@ TransientSolutionEnum, /*}}}*/ - /*Approximations {{{*/ + /*Approximations {{{1*/ ApproximationEnum, HutterApproximationEnum, @@ -231,5 +239,5 @@ StokesApproximationEnum, /*}}}*/ - /*Datasets {{{*/ + /*Datasets {{{1*/ ConstraintsEnum, LoadsEnum, @@ -240,5 +248,5 @@ ResultsEnum, /*}}}*/ - /*Objects {{{*/ + /*Objects {{{1*/ BoolInputEnum, BoolParamEnum, @@ -289,5 +297,5 @@ VertexEnum, /*}}}*/ - /*Fill {{{*/ + /*Fill {{{1*/ AirEnum, IceEnum, @@ -295,10 +303,10 @@ WaterEnum, /*}}}*/ - /*Rift state {{{*/ + /*Rift state {{{1*/ ClosedEnum, FreeEnum, OpenEnum, /*}}}*/ - /*Inputs {{{*/ + /*Inputs {{{1*/ AdjointpEnum, AdjointxEnum, @@ -388,10 +396,10 @@ IceVolumeEnum, /*}}}*/ - /*Element Interpolations{{{*/ + /*Element Interpolations{{{1*/ P0Enum, P1Enum, P1DGEnum, /*}}}*/ - /*Results{{{*/ + /*Results{{{1*/ SaveResultsEnum, BoolElementResultEnum, @@ -414,5 +422,5 @@ WaterColumnOldEnum, /*}}}*/ - /*Responses{{{*/ + /*Responses{{{1*/ MinVelEnum, MaxVelEnum, @@ -427,18 +435,18 @@ MaxAbsVzEnum, /*}}}*/ - /*Convergence{{{*/ + /*Convergence{{{1*/ RelativeEnum, AbsoluteEnum, IncrementalEnum, /*}}}*/ - /*Grounding Line{{{*/ + /*Grounding Line{{{1*/ AgressiveMigrationEnum, NoneEnum, SoftMigrationEnum, /*}}}*/ - /*Solver{{{*/ + /*Solver{{{1*/ StokesSolverEnum, /*}}}*/ - /*Parameters{{{*/ + /*Parameters{{{1*/ AdjointEnum, ColinearEnum, @@ -468,14 +476,14 @@ VerboseEnum, /*}}}*/ - /*Interpolation {{{*/ + /*Interpolation {{{1*/ TriangleInterpEnum, BilinearInterpEnum, NearestInterpEnum, /*}}}*/ - /*Coordinate Systems{{{*/ + /*Coordinate Systems{{{1*/ XYEnum, XYZPEnum, /*}}}*/ - /*Options{{{*/ + /*Options{{{1*/ OptionEnum, OptionCellEnum, @@ -485,5 +493,5 @@ OptionLogicalEnum, /*}}}*/ - /*Rheology law (move too Material) {{{*/ + /*Rheology law (move too Material) {{{1*/ PatersonEnum, ArrheniusEnum, Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12676) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12677) @@ -284,4 +284,6 @@ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.h\ ./modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp\ + ./modules/SmbGradientsx/SmbGradientsx.h\ + ./modules/SmbGradientsx/SmbGradientsx.cpp\ ./modules/UpdateConstraintsx/UpdateConstraintsx.h\ ./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\ Index: /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 12677) @@ -167,5 +167,13 @@ case SurfaceforcingsMassBalanceEnum : return "SurfaceforcingsMassBalance"; case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd"; + case SurfaceforcingsIssmbgradientsEnum : return "SurfaceforcingsIssmbgradients"; case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures"; + case SurfaceforcingsHcEnum : return "SurfaceforcingsHc"; + case SurfaceforcingsSmbPosMaxEnum : return "SurfaceforcingsSmbPosMax"; + case SurfaceforcingsSmbPosMinEnum : return "SurfaceforcingsSmbPosMin"; + case SurfaceforcingsAPosEnum : return "SurfaceforcingsAPos"; + case SurfaceforcingsBPosEnum : return "SurfaceforcingsBPos"; + case SurfaceforcingsANegEnum : return "SurfaceforcingsANeg"; + case SurfaceforcingsBNegEnum : return "SurfaceforcingsBNeg"; case ThermalMaxiterEnum : return "ThermalMaxiter"; case ThermalPenaltyFactorEnum : return "ThermalPenaltyFactor"; @@ -467,5 +475,5 @@ len=strlen(EnumToStringx(enum_in)); string=xNew(len+1); - xMemCpy(string,EnumToStringx(enum_in),(len+1)); + memcpy(string,EnumToStringx(enum_in),(len+1)*sizeof(char)); /*Assign output pointer*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 12677) @@ -90,4 +90,5 @@ parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum)); parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum)); + parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIssmbgradientsEnum)); /*some parameters that did not come with the iomodel: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp (revision 12677) @@ -21,4 +21,5 @@ bool dakota_analysis; bool ispdd; + bool issmbgradients; /*Fetch data needed: */ @@ -29,4 +30,5 @@ iomodel->FetchData(1,MeshElementsEnum); iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum); + iomodel->Constant(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); /*Update elements: */ @@ -73,4 +75,13 @@ iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum); } + if(issmbgradients){ + iomodel->FetchDataToInput(elements,SurfaceforcingsHcEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMaxEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsSmbPosMinEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsAPosEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBPosEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsANegEnum); + iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum); + } else{ iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum); Index: /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.cpp (revision 12677) +++ /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.cpp (revision 12677) @@ -0,0 +1,27 @@ +/*!\file SmbGradientsx + * \brief: calculates SMB as function of local elevation + */ + +#include "./SmbGradientsx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../io/io.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void SmbGradientsx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters){ + +// void SmbGradientsx(hd,agd,ni){ +// INPUT parameters: ni: working size of arrays +// INPUT: surface elevation (m): hd(NA) +// OUTPUT: mass-balance (m/yr ice): agd(NA) + + int i; + + Element* element = NULL; + + for(i=0;iSize();i++){ + element=(Element*)elements->GetObjectByOffset(i); + element->SmbGradients(); + } +} Index: /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.h (revision 12677) +++ /issm/trunk-jpl/src/c/modules/SmbGradientsx/SmbGradientsx.h (revision 12677) @@ -0,0 +1,14 @@ +/*!\file: SmbGradientsx.h + * \brief header file for degree of freedoms distribution routines. + */ + +#ifndef _SMBGRADIENTSX_H +#define _SMBGRADIENTSX_H + +#include "../../Container/Container.h" +#include "../../objects/objects.h" + +/* local prototypes: */ +void SmbGradientsx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters); + +#endif /* _SMBGRADIENTSX_H*/ Index: /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 12677) @@ -171,5 +171,13 @@ else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; + else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum; + else if (strcmp(name,"SurfaceforcingsHc")==0) return SurfaceforcingsHcEnum; + else if (strcmp(name,"SurfaceforcingsSmbPosMax")==0) return SurfaceforcingsSmbPosMaxEnum; + else if (strcmp(name,"SurfaceforcingsSmbPosMin")==0) return SurfaceforcingsSmbPosMinEnum; + else if (strcmp(name,"SurfaceforcingsAPos")==0) return SurfaceforcingsAPosEnum; + else if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum; + else if (strcmp(name,"SurfaceforcingsANeg")==0) return SurfaceforcingsANegEnum; + else if (strcmp(name,"SurfaceforcingsBNeg")==0) return SurfaceforcingsBNegEnum; else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum; else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum; @@ -253,5 +261,8 @@ else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum; else if (strcmp(name,"DoubleVecParam")==0) return DoubleVecParamEnum; - else if (strcmp(name,"Element")==0) return ElementEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"Element")==0) return ElementEnum; else if (strcmp(name,"ElementResult")==0) return ElementResultEnum; else if (strcmp(name,"ExternalResult")==0) return ExternalResultEnum; @@ -261,8 +272,5 @@ else if (strcmp(name,"Input")==0) return InputEnum; else if (strcmp(name,"IntInput")==0) return IntInputEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"IntParam")==0) return IntParamEnum; + else if (strcmp(name,"IntParam")==0) return IntParamEnum; else if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum; else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum; @@ -376,5 +384,8 @@ else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum; else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum; - else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum; else if (strcmp(name,"StressTensor")==0) return StressTensorEnum; else if (strcmp(name,"StressTensorxx")==0) return StressTensorxxEnum; @@ -384,8 +395,5 @@ else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; + else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum; else if (strcmp(name,"P0")==0) return P0Enum; else if (strcmp(name,"P1")==0) return P1Enum; Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 12676) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 12677) @@ -106,4 +106,5 @@ #include "./RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h" #include "./Scotchx/Scotchx.h" +#include "./SmbGradientsx/SmbGradientsx.h" #include "./Solverx/Solverx.h" #include "./SpcNodesx/SpcNodesx.h" Index: /issm/trunk-jpl/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12676) +++ /issm/trunk-jpl/src/c/objects/Elements/Element.h (revision 12677) @@ -72,4 +72,5 @@ virtual void PotentialSheetUngrounding(Vector* potential_sheet_ungrounding)=0; virtual void PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm)=0; + virtual void SmbGradients()=0; virtual int UpdatePotentialSheetUngrounding(IssmDouble* potential_sheet_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf)=0; virtual void ResetCoordinateSystem()=0; Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.cpp (revision 12677) @@ -2738,4 +2738,59 @@ } /*}}}*/ +/*FUNCTION Penta::SmbGradients{{{*/ +void Penta::SmbGradients(void){ + + int i; + + // input + IssmDouble h[NUMVERTICES]; // ice thickness (m) + IssmDouble s[NUMVERTICES]; // surface elevation (m) + IssmDouble a_pos[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble b_pos[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble a_neg[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble b_neg[NUMVERTICES]; // Hs-SMB relation paremeter + IssmDouble Hc[NUMVERTICES]; // elevation of transition between accumulation regime and ablation regime + IssmDouble smb_pos_max[NUMVERTICES]; // maximum SMB value in the accumulation regime + IssmDouble smb_pos_min[NUMVERTICES]; // minimum SMB value in the accumulation regime + IssmDouble rho_water; // density of fresh water + IssmDouble rho_ice; // density of ice + + // output + IssmDouble smb[NUMVERTICES]; // surface mass balance (m/yr ice) + + /*Recover SmbGradients*/ + GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum); + GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum); + GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum); + GetInputListOnVertices(&a_pos[0],SurfaceforcingsAPosEnum); + GetInputListOnVertices(&b_pos[0],SurfaceforcingsBPosEnum); + GetInputListOnVertices(&a_neg[0],SurfaceforcingsANegEnum); + GetInputListOnVertices(&b_neg[0],SurfaceforcingsBNegEnum); + + /*Recover surface elevatio at vertices: */ + GetInputListOnVertices(&h[0],ThicknessEnum); + GetInputListOnVertices(&s[0],SurfaceEnum); + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + + // loop over all vertices + for(i=0;iHc[i]){ + smb[i]=a_pos[i]+b_pos[i]*s[i]; + if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} + if(smb[i]inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0])); +} +/*}}}*/ /*FUNCTION Penta::SurfaceArea {{{*/ IssmDouble Penta::SurfaceArea(void){ Index: /issm/trunk-jpl/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12676) +++ /issm/trunk-jpl/src/c/objects/Elements/Penta.h (revision 12677) @@ -110,4 +110,5 @@ void ProcessResultsUnits(void); void ResetCoordinateSystem(void); + void SmbGradients(); IssmDouble SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.cpp (revision 12677) @@ -2369,4 +2369,58 @@ } /*}}}*/ +/*FUNCTION Tria::SmbGradients{{{*/ +void Tria::SmbGradients(void){ + + int i; + + // input + IssmDouble h[NUMVERTICES]; // ice thickness (m) + IssmDouble s[NUMVERTICES]; // surface elevation (m) + IssmDouble a_pos[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble b_pos[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble a_neg[NUMVERTICES]; // Hs-SMB relation parameter + IssmDouble b_neg[NUMVERTICES]; // Hs-SMB relation paremeter + IssmDouble Hc[NUMVERTICES]; // elevation of transition between accumulation regime and ablation regime + IssmDouble smb_pos_max[NUMVERTICES]; // maximum SMB value in the accumulation regime + IssmDouble smb_pos_min[NUMVERTICES]; // minimum SMB value in the accumulation regime + IssmDouble rho_water; // density of fresh water + IssmDouble rho_ice; // density of ice + + // output + IssmDouble smb[NUMVERTICES]; // surface mass balance (m/yr ice) + + /*Recover SmbGradients*/ + GetInputListOnVertices(&Hc[0],SurfaceforcingsHcEnum); + GetInputListOnVertices(&smb_pos_max[0],SurfaceforcingsSmbPosMaxEnum); + GetInputListOnVertices(&smb_pos_min[0],SurfaceforcingsSmbPosMinEnum); + GetInputListOnVertices(&a_pos[0],SurfaceforcingsAPosEnum); + GetInputListOnVertices(&b_pos[0],SurfaceforcingsBPosEnum); + GetInputListOnVertices(&a_neg[0],SurfaceforcingsANegEnum); + GetInputListOnVertices(&b_neg[0],SurfaceforcingsBNegEnum); + + /*Recover surface elevatio at vertices: */ + GetInputListOnVertices(&h[0],ThicknessEnum); + GetInputListOnVertices(&s[0],SurfaceEnum); + + /*Get material parameters :*/ + rho_ice=matpar->GetRhoIce(); + rho_water=matpar->GetRhoFreshwater(); + + // loop over all vertices + for(i=0;iHc[i]){ + smb[i]=a_pos[i]+b_pos[i]*s[i]; + if(smb[i]>smb_pos_max[i]){smb[i]=smb_pos_max[i];} + if(smb[i]inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&smb[0])); +} +/*}}}*/ /*FUNCTION Tria::SetCurrentConfiguration {{{*/ void Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){ Index: /issm/trunk-jpl/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12676) +++ /issm/trunk-jpl/src/c/objects/Elements/Tria.h (revision 12677) @@ -110,4 +110,5 @@ void ProcessResultsUnits(void); void ResetCoordinateSystem(void){_error2_("not implemented yet");}; + void SmbGradients(); IssmDouble SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); Index: /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/shared/Numerics/UnitConversion.cpp (revision 12677) @@ -66,4 +66,10 @@ case SurfaceforcingsPrecipitationEnum: scale=yts;break; //m/yr case SurfaceforcingsMassBalanceEnum: scale=yts;break; //m/yr + case SurfaceforcingsSmbPosMaxEnum: scale=yts;break; //m/yr + case SurfaceforcingsSmbPosMinEnum: scale=yts;break; //m/yr + case SurfaceforcingsAPosEnum: scale=yts;break; //m/yr + case SurfaceforcingsBPosEnum: scale=yts;break; //m/yr + case SurfaceforcingsANegEnum: scale=yts;break; //m/yr + case SurfaceforcingsBNegEnum: scale=yts;break; //m/yr case MisfitEnum: scale=pow(yts,2);break; //(m/yr)^2 case MassFluxEnum: scale=pow((IssmDouble)10,-12)*yts;break; // (GigaTon/year) Index: /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 12676) +++ /issm/trunk-jpl/src/c/solutions/prognostic_core.cpp (revision 12677) @@ -17,4 +17,5 @@ bool save_results; bool ispdd; + bool issmbgradients; /*activate formulation: */ @@ -24,9 +25,15 @@ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); + femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); if(ispdd){ if(VerboseSolution()) _pprintLine_(" call positive degree day module"); PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - } + } + + if(issmbgradients){ + _printf_(VerboseSolution()," call smb gradients module\n"); + SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } if(VerboseSolution()) _pprintLine_(" call computational core"); Index: /issm/trunk-jpl/src/m/classes/surfaceforcings.m =================================================================== --- /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 12676) +++ /issm/trunk-jpl/src/m/classes/surfaceforcings.m (revision 12677) @@ -9,4 +9,12 @@ mass_balance = NaN; ispdd = 0; + issmbgradients = 0; + hc = NaN; + smb_pos_max = NaN; + smb_pos_min = NaN; + a_pos = NaN; + b_pos = NaN; + a_neg = NaN; + b_neg = NaN; monthlytemperatures = NaN; end @@ -24,4 +32,5 @@ %pdd method not used in default mode obj.ispdd=0; + obj.issmbgradients=0; end % }}} @@ -30,6 +39,15 @@ if ismember(PrognosticAnalysisEnum,analyses), md = checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0 1]); + checkfield(md,'surfaceforcings.issmbgradients','numel',1,'values',[0 1]); if(obj.ispdd) md = checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1); + elseif(obj.issmbgradients) + checkfield(md,'surfaceforcings.hc','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.smb_pos_max','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.smb_pos_min','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.a_pos','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.b_pos','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.a_neg','forcing',1,'NaN',1); + checkfield(md,'surfaceforcings.b_neg','forcing',1,'NaN',1); else md = checkfield(md,'surfaceforcings.mass_balance','forcing',1,'NaN',1); @@ -47,4 +65,12 @@ fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)'); fielddisplay(obj,'monthlytemperatures','monthly surface temperatures required if pdd is activated'); + fielddisplay(obj,'issmbgradients','is smb gradients method activated (0 or 1, default is 0)'); + fielddisplay(obj,'hc',' elevation of intersection between accumulation and ablation regime required if smb gradients is activated'); + fielddisplay(obj,'smb_pos_max',' maximum value of positive smb required if smb gradients is activated'); + fielddisplay(obj,'smb_pos_min',' minimum value of positive smb required if smb gradients is activated'); + fielddisplay(obj,'a_pos',' intercept of hs - smb regression line for accumulation regime required if smb gradients is activated'); + fielddisplay(obj,'b_pos',' slope of hs - smb regression line for accumulation regime required if smb gradients is activated'); + fielddisplay(obj,'a_neg',' intercept of hs - smb regression line for ablation regime required if smb gradients is activated'); + fielddisplay(obj,'b_neg',' slope of hs - smb regression line for ablation regime required if smb gradients is activated'); end % }}} @@ -56,4 +82,14 @@ WriteData(fid,'object',obj,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1); end + WriteData(fid,'object',obj,'fieldname','issmbgradients','format','Boolean'); + if obj.issmbgradients, + WriteData(fid,'object',obj,'fieldname','hc','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','smb_pos_max','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','smb_pos_min','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','a_pos','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','b_pos','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','a_neg','format','DoubleMat','mattype',1); + WriteData(fid,'object',obj,'fieldname','b_neg','format','DoubleMat','mattype',1); + end end % }}} Index: /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 12676) +++ /issm/trunk-jpl/src/m/enum/MaximumNumberOfEnums.m (revision 12677) @@ -9,3 +9,3 @@ % macro=MaximumNumberOfEnums() -macro=440; +macro=448;