Changeset 27354 for issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp

Timestamp:

11/04/22 04:51:01 (2 years ago)

Author:

vverjans

Message:

NEW added autoregression and stochastic capabilities for subglacial discharge

File:

• issm/trunk-jpl/src/c/modules/FrontalForcingsx/FrontalForcingsx.cpp (modified) (2 diffs)

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

@@ -20 +20 @@
                 case FrontalForcingsRignotarmaEnum:
                         Thermalforcingarmax(femmodel);
+                        bool isdischargearma;
+                        femmodel->parameters->FindParam(&isdischargearma,FrontalForcingsIsDischargeARMAEnum);
+                        if(isdischargearma==true) Subglacialdischargearmax(femmodel);
                         /*Do not break here, call IcefrontAreax(),RignotMeltParameterizationx()*/
                 case FrontalForcingsRignotEnum:
@@ -105 +108 @@
    xDelete<IssmDouble>(monthtrends);
 }/*}}}*/
+void Subglacialdischargearmax(FemModel* femmodel){/*{{{*/
+
+        /*Get time parameters*/
+   IssmDouble time,dt,starttime,tstep_arma;
+   femmodel->parameters->FindParam(&time,TimeEnum);
+   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
+   femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
+   femmodel->parameters->FindParam(&tstep_arma,FrontalForcingsSdARMATimestepEnum);
+
+   /*Determine if this is a time step for the ARMA model*/
+   bool isstepforarma = false;
+
+   #ifndef _HAVE_AD_
+   if((fmod(time,tstep_arma)<fmod((time-dt),tstep_arma)) || (time<=starttime+dt) || tstep_arma==dt) isstepforarma = true;
+   #else
+   _error_("not implemented yet");
+   #endif
+
+   /*Load parameters*/
+        bool isstochastic;
+   bool isdischargestochastic = false;
+        int M,N,arorder,maorder,numbasins,numparams,numbreaks,my_rank;
+   femmodel->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
+   femmodel->parameters->FindParam(&numparams,FrontalForcingsSdNumberofParamsEnum);
+   femmodel->parameters->FindParam(&numbreaks,FrontalForcingsSdNumberofBreaksEnum);
+   femmodel->parameters->FindParam(&arorder,FrontalForcingsSdarOrderEnum);
+   femmodel->parameters->FindParam(&maorder,FrontalForcingsSdmaOrderEnum);
+   IssmDouble* datebreaks        = NULL;
+   IssmDouble* arlagcoefs        = NULL;
+   IssmDouble* malagcoefs        = NULL;
+   IssmDouble* monthlyfrac       = NULL;
+        IssmDouble* polyparams        = NULL;
+
+   femmodel->parameters->FindParam(&datebreaks,&M,&N,FrontalForcingsSddatebreaksEnum);            _assert_(M==numbasins); _assert_(N==max(numbreaks,1));        
+   femmodel->parameters->FindParam(&polyparams,&M,&N,FrontalForcingsSdpolyparamsEnum);            _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams);        
+   femmodel->parameters->FindParam(&arlagcoefs,&M,&N,FrontalForcingsSdarlagcoefsEnum);            _assert_(M==numbasins); _assert_(N==arorder);
+   femmodel->parameters->FindParam(&malagcoefs,&M,&N,FrontalForcingsSdmalagcoefsEnum);            _assert_(M==numbasins); _assert_(N==maorder);
+   femmodel->parameters->FindParam(&monthlyfrac,&M,&N,FrontalForcingsSdMonthlyFracEnum);          _assert_(M==numbasins); _assert_(N==12); 
+
+        femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum);
+        if(isstochastic){
+                int  numstochasticfields;
+                int* stochasticfields;
+                femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum);
+                femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields);
+                for(int i=0;i<numstochasticfields;i++){
+                        if(stochasticfields[i]==FrontalForcingsSubglacialDischargearmaEnum) isdischargestochastic = true;
+                }
+                xDelete<int>(stochasticfields);
+        }
+
+   /*Loop over each element to compute Subglacial Discharge at vertices*/
+   for(Object* &object:femmodel->elements->objects){
+      Element* element = xDynamicCast<Element*>(object);
+                /*Compute ARMA*/
+      element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,isdischargestochastic,FrontalForcingsSubglacialDischargearmaEnum);
+                /*Scale with monthly fractions*/
+                element->MonthlyFractionBasin(monthlyfrac,FrontalForcingsSubglacialDischargeEnum);
+        }
+
+   /*Cleanup*/
+   xDelete<IssmDouble>(arlagcoefs);
+   xDelete<IssmDouble>(malagcoefs);
+   xDelete<IssmDouble>(monthlyfrac);
+   xDelete<IssmDouble>(polyparams);
+   xDelete<IssmDouble>(datebreaks);
+}/*}}}*/
+