Changeset 27285


Ignore:
Timestamp:
09/17/22 06:51:49 (3 years ago)
Author:
vverjans
Message:

NEW: added monthly effect capabilities to frontalforcingsrignotarma

Location:
issm/trunk-jpl/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

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

    r27284 r27285  
    256256         iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.malag_coefs");
    257257         parameters->AddObject(new DoubleMatParam(FrontalForcingsARMAmalagcoefsEnum,transparam,M,N));
     258         xDelete<IssmDouble>(transparam);
     259                        iomodel->FetchData(&transparam,&M,&N,"md.frontalforcings.monthly_effects");
     260         parameters->AddObject(new DoubleMatParam(ThermalForcingMonthlyEffectsEnum,transparam,M,N));
    258261         xDelete<IssmDouble>(transparam);
    259262                        /*Do not break here, generic FrontalForcingsRignot parameters still to be retrieved*/

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

    r27284 r27285  
    27752775        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
    27762776
     2777}/*}}}*/
     2778void       Element::MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type){/*{{{*/
     2779       
     2780        /*Variable declaration*/
     2781   const int numvertices = this->GetNumberOfVertices();
     2782        int basinid,mindex,basinenum_type,varenum_type;
     2783   IssmDouble time,fracyear,yts;
     2784   IssmDouble monthsteps[12] = {0.,1./12,2./12,3./12,4./12,5./12,6./12,7./12,8./12,9./12,10./12,11./12};
     2785        IssmDouble* monthlyeff_b = xNew<IssmDouble>(12);
     2786   IssmDouble* varlist      = xNew<IssmDouble>(numvertices);
     2787
     2788        /*Get field-specific enums*/
     2789   switch(enum_type){
     2790      case(FrontalForcingsRignotarmaEnum):
     2791         basinenum_type = FrontalForcingsBasinIdEnum;
     2792         varenum_type   = ThermalForcingEnum;
     2793         break;
     2794        }
     2795       
     2796        /*Evaluate the month index*/
     2797        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2798        this->parameters->FindParam(&time,TimeEnum);
     2799        fracyear = time/yts-floor(time/yts);
     2800        for(int i=1;i<12;i++){
     2801                if(fracyear>=monthsteps[i-1]) mindex = i-1;
     2802        }
     2803        if(fracyear>=monthsteps[11])     mindex = 11;
     2804
     2805        /*Get basin-specific parameters of the element*/
     2806   this->GetInputValue(&basinid,basinenum_type);
     2807        for(int ii=0;ii<12;ii++){
     2808                monthlyeff_b[ii] = monthlyeff[basinid*12+ii];
     2809        }
     2810
     2811        /*Retrieve values non-adjusted for monthly effects*/
     2812   this->GetInputListOnVertices(varlist,varenum_type);
     2813
     2814        /*Adjust values using monthly effects*/
     2815        for(int v=0;v<numvertices;v++){
     2816                varlist[v] = varlist[v]+monthlyeff_b[mindex];
     2817        }
     2818
     2819        /*Add input to element*/
     2820   this->AddInput(varenum_type,varlist,P1Enum);
     2821
     2822        /*Clean-up*/
     2823   xDelete<IssmDouble>(monthlyeff_b);
     2824        xDelete<IssmDouble>(varlist);
    27772825}/*}}}*/
    27782826void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/

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

    r27260 r27285  
    164164                void               MigrateGroundingLine(IssmDouble* sheet_ungrounding);
    165165                void               MismipFloatingiceMeltingRate();
     166                void               MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type);
    166167                void               BeckmannGoosseFloatingiceMeltingRate();
    167168                void               MungsmtpParameterization(void);

  • issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

    r27261 r27285  
    5959   IssmDouble* arlagcoefs    = NULL;
    6060   IssmDouble* malagcoefs    = NULL;
     61   IssmDouble* monthlyeff    = NULL;
    6162
    6263        femmodel->parameters->FindParam(&tinit_arma,FrontalForcingsARMAInitialTimeEnum);
    … …  
    6566   femmodel->parameters->FindParam(&arlagcoefs,&M,&Narlagcoefs,FrontalForcingsARMAarlagcoefsEnum);  _assert_(M==numbasins); _assert_(Narlagcoefs==arorder);
    6667   femmodel->parameters->FindParam(&malagcoefs,&M,&Nmalagcoefs,FrontalForcingsARMAmalagcoefsEnum);  _assert_(M==numbasins); _assert_(Nmalagcoefs==maorder);
     68   femmodel->parameters->FindParam(&monthlyeff,&M,&N,ThermalForcingMonthlyEffectsEnum);             _assert_(M==numbasins); _assert_(N==12);
    6769
    6870        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
    … …  
    8385   for(Object* &object:femmodel->elements->objects){
    8486      Element* element = xDynamicCast<Element*>(object);
     87                /*Compute ARMA*/
    8588      element->ArmaProcess(isstepforarma,arorder,maorder,telapsed_arma,tstep_arma,termconstant,trend,arlagcoefs,malagcoefs,istfstochastic,FrontalForcingsRignotarmaEnum);
    86    }
     89                /*Compute monthly effects*/
     90                element->MonthlyEffectBasin(monthlyeff,FrontalForcingsRignotarmaEnum);
     91        }
    8792
    8893   /*Cleanup*/
    … …  
    9196   xDelete<IssmDouble>(arlagcoefs);
    9297   xDelete<IssmDouble>(malagcoefs);
     98   xDelete<IssmDouble>(monthlyeff);
    9399}/*}}}*/

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

    r27284 r27285  
    564564syn keyword cConstant StressbalanceRiftPenaltyThresholdEnum
    565565syn keyword cConstant StressbalanceShelfDampeningEnum
     566syn keyword cConstant ThermalForcingMonthlyEffectsEnum
    566567syn keyword cConstant ThermalIsdrainicecolumnEnum
    567568syn keyword cConstant ThermalIsdynamicbasalspcEnum

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

    r27284 r27285  
    558558        StressbalanceRiftPenaltyThresholdEnum,
    559559        StressbalanceShelfDampeningEnum,
     560   ThermalForcingMonthlyEffectsEnum,
    560561        ThermalIsdrainicecolumnEnum,
    561562        ThermalIsdynamicbasalspcEnum,

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

    r27284 r27285  
    566566                case StressbalanceRiftPenaltyThresholdEnum : return "StressbalanceRiftPenaltyThreshold";
    567567                case StressbalanceShelfDampeningEnum : return "StressbalanceShelfDampening";
     568                case ThermalForcingMonthlyEffectsEnum : return "ThermalForcingMonthlyEffects";
    568569                case ThermalIsdrainicecolumnEnum : return "ThermalIsdrainicecolumn";
    569570                case ThermalIsdynamicbasalspcEnum : return "ThermalIsdynamicbasalspc";

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

    r27284 r27285  
    557557syn keyword juliaConstC StressbalanceRiftPenaltyThresholdEnum
    558558syn keyword juliaConstC StressbalanceShelfDampeningEnum
     559syn keyword juliaConstC ThermalForcingMonthlyEffectsEnum
    559560syn keyword juliaConstC ThermalIsdrainicecolumnEnum
    560561syn keyword juliaConstC ThermalIsdynamicbasalspcEnum

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

    r27284 r27285  
    578578              else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
    579579              else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
     580              else if (strcmp(name,"ThermalForcingMonthlyEffects")==0) return ThermalForcingMonthlyEffectsEnum;
    580581              else if (strcmp(name,"ThermalIsdrainicecolumn")==0) return ThermalIsdrainicecolumnEnum;
    581582              else if (strcmp(name,"ThermalIsdynamicbasalspc")==0) return ThermalIsdynamicbasalspcEnum;
    … …  
    628629              else if (strcmp(name,"Xxe")==0) return XxeEnum;
    629630              else if (strcmp(name,"Yye")==0) return YyeEnum;
    630               else if (strcmp(name,"Zze")==0) return ZzeEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Areae")==0) return AreaeEnum;
     634              if (strcmp(name,"Zze")==0) return ZzeEnum;
     635              else if (strcmp(name,"Areae")==0) return AreaeEnum;
    635636              else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
    636637              else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum;
    … …  
    751752              else if (strcmp(name,"DrivingStressX")==0) return DrivingStressXEnum;
    752753              else if (strcmp(name,"DrivingStressY")==0) return DrivingStressYEnum;
    753               else if (strcmp(name,"Dummy")==0) return DummyEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
     757              if (strcmp(name,"Dummy")==0) return DummyEnum;
     758              else if (strcmp(name,"EffectivePressure")==0) return EffectivePressureEnum;
    758759              else if (strcmp(name,"EffectivePressureSubstep")==0) return EffectivePressureSubstepEnum;
    759760              else if (strcmp(name,"EffectivePressureTransient")==0) return EffectivePressureTransientEnum;
    … …  
    874875              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    875876              else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
    876               else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
     880              if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
     881              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    881882              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    882883              else if (strcmp(name,"Node")==0) return NodeEnum;
    … …  
    997998              else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
    998999              else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
    999               else if (strcmp(name,"SmbD")==0) return SmbDEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
     1003              if (strcmp(name,"SmbD")==0) return SmbDEnum;
     1004              else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
    10041005              else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
    10051006              else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
    … …  
    11201121              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    11211122              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    1122               else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
     1126              if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
     1127              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    11271128              else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
    11281129              else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
    … …  
    12431244              else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum;
    12441245              else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
    1245               else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
     1249              if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
     1250              else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
    12501251              else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
    12511252              else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
    … …  
    13661367              else if (strcmp(name,"Element")==0) return ElementEnum;
    13671368              else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
    1368               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
     1372              if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
     1373              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
    13731374              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
    13741375              else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
    … …  
    14891490              else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
    14901491              else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
    1491               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
     1495              if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
     1496              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
    14961497              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
    14971498              else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
    … …  
    16121613              else if (strcmp(name,"SmoothAnalysis")==0) return SmoothAnalysisEnum;
    16131614              else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
    1614               else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
     1618              if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
     1619              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
    16191620              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
    16201621              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;

  • issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.m

    r27260 r27285  
    1515                arlag_coefs          = NaN;
    1616                malag_coefs          = NaN;
     17                monthly_effects      = NaN;
    1718                basin_id             = NaN;
    1819                subglacial_discharge = NaN;
    … …  
    6465                        md = checkfield(md,'fieldname','frontalforcings.arlag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ar_order]);
    6566                        md = checkfield(md,'fieldname','frontalforcings.malag_coefs','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,md.frontalforcings.ma_order]);
     67                        if(any(~isnan(md.frontalforcings.monthly_effects)))
     68                                md = checkfield(md,'fieldname','frontalforcings.monthly_effects','NaN',1,'Inf',1,'size',[md.frontalforcings.num_basins,12]);
     69                                if(md.timestepping.time_step>=1)
     70                                        error('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1');
     71                                end
     72                        end
    6673
    6774                end % }}}
    … …  
    7986         fielddisplay(self,'arlag_coefs','basin-specific vectors of AR lag coefficients [unitless]');
    8087         fielddisplay(self,'malag_coefs','basin-specific vectors of MA lag coefficients [unitless]');
     88         fielddisplay(self,'monthly_effects','basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]');
    8189                end % }}}
    8290                function marshall(self,prefix,md,fid) % {{{
    8391                        yts=md.constants.yts;
     92                        if(any(isnan(md.frontalforcings.monthly_effects))) %monthly effects not provided, set to 0
     93                                meffects = zeros(md.frontalforcings.num_basins,12);
     94                        else
     95                                meffects = md.frontalforcings.monthly_effects;
     96                        end
     97
    8498                        WriteData(fid,prefix,'name','md.frontalforcings.parameterization','data',3,'format','Integer');
    8599                        WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','num_basins','format','Integer');
    … …  
    94108         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','arlag_coefs','format','DoubleMat','name','md.frontalforcings.arlag_coefs','yts',yts);
    95109         WriteData(fid,prefix,'object',self,'class','frontalforcings','fieldname','malag_coefs','format','DoubleMat','name','md.frontalforcings.malag_coefs','yts',yts);
     110         WriteData(fid,prefix,'data',meffects,'name','md.frontalforcings.monthly_effects','format','DoubleMat');
    96111                end % }}}
    97112        end

  • issm/trunk-jpl/src/m/classes/frontalforcingsrignotarma.py

    r27260 r27285  
    2222        self.arma_timestep = 0
    2323        self.arlag_coefs = np.nan
     24        self.malag_coefs = np.nan
     25        self.monthly_effects = np.nan
    2426        self.basin_id = np.nan
    2527        self.subglacial_discharge = np.nan
    … …  
    4345        s += '{}\n'.format(fielddisplay(self, 'arlag_coefs', 'basin-specific vectors of AR lag coefficients [unitless]'))
    4446        s += '{}\n'.format(fielddisplay(self, 'malag_coefs', 'basin-specific vectors of MA lag coefficients [unitless]'))
     47        s += '{}\n'.format(fielddisplay(self, 'monthly_effects', 'basin-specific monthly values of TF added at corresponding month (default: all 0) [°C]'))
    4548        return s
    4649    #}}}
    … …  
    7477        md = checkfield(md, 'fieldname', 'frontalforcings.arlag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ar_order])
    7578        md = checkfield(md, 'fieldname', 'frontalforcings.malag_coefs', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, md.frontalforcings.ma_order])
     79        if(np.any(np.isnan(md.frontalforcings.monthly_effects)==False)):
     80            md = checkfield(md, 'fieldname', 'frontalforcings.monthly_effects', 'NaN', 1, 'Inf', 1, 'size', [md.frontalforcings.num_basins, 12])
     81            if(md.timestepping.time_step>=1):
     82                raise RuntimeError('md.frontalforcings.monthly_effects are provided but md.timestepping.time_step>=1')
    7683        return md
    7784    # }}}
    … …  
    8491    def marshall(self, prefix, md, fid):  # {{{
    8592        yts = md.constants.yts
     93        if(np.any(np.isnan(md.frontalforcings.monthly_effects))): #monthly effects not provided, set to 0
     94            meffects = np.zeros((md.frontalforcings.num_basins,12))
     95        else:
     96            meffects = 1*md.frontalforcings.monthly_effects
    8697        WriteData(fid, prefix, 'name', 'md.frontalforcings.parameterization', 'data', 3, 'format', 'Integer')
    8798        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'num_basins', 'format', 'Integer')
    … …  
    96107        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'arlag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.arlag_coefs', 'yts', yts)
    97108        WriteData(fid, prefix, 'object', self, 'class', 'frontalforcings', 'fieldname', 'malag_coefs', 'format', 'DoubleMat', 'name', 'md.frontalforcings.malag_coefs', 'yts', yts)
     109        WriteData(fid, prefix, 'data', meffects, 'name', 'md.frontalforcings.monthly_effects', 'format', 'DoubleMat')
    98110    # }}}
Note: See TracChangeset for help on using the changeset viewer.