Changeset 27956


Ignore:
Timestamp:
10/10/23 15:27:28 (18 months ago)
Author:
dlcheng
Message:

CHG: Collecting dlcheng-ASE branch changes to misfit calculation. Changes work with r27296 and not from the current revision r27954.

Location:
issm/branches/trunk-dlcheng-ASE/src
Files:
2 added
25 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/Makefile.am

    r27944 r27956  
    8181        ./classes/Numberedcostfunction.cpp \
    8282        ./classes/Misfit.cpp \
     83        ./classes/MisfitAnnual.cpp \
    8384        ./classes/Cfsurfacesquare.cpp \
    8485        ./classes/Cfsurfacesquaretransient.cpp \

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.cpp

    r27940 r27956  
    3030extern "C" void run_semic_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
    3131                        IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out);
    32 
    33 extern "C" void run_semic_transient_(int *nx, int *ntime, int *nloop,
    34                         IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in,
    35                         IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
    36                         IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *qmr_in,
    37                         IssmDouble *tstic,
    38                         IssmDouble *hcrit, IssmDouble *rcrit,
    39                         IssmDouble *mask, IssmDouble *hice, IssmDouble *hsnow,
    40                         IssmDouble *albedo_in, IssmDouble *albedo_snow_in,
    41                         int *alb_scheme, IssmDouble *alb_smax, IssmDouble *alb_smin, IssmDouble *albi, IssmDouble *albl,
    42                         IssmDouble *Tamp,
    43                         IssmDouble *tmin, IssmDouble *tmax, IssmDouble *tmid, IssmDouble *mcrit, IssmDouble *wcrit, IssmDouble *tau_a, IssmDouble* tau_f, IssmDouble *afac, bool *verbose,
    44                         IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *smbi_out, IssmDouble *smbs_out, IssmDouble *saccu_out, IssmDouble *smelt_out, IssmDouble *refr_out, IssmDouble *albedo_out, IssmDouble *albedo_snow_out, IssmDouble *hsnow_out, IssmDouble *hice_out, IssmDouble *qmr_out);
    4532#endif
    4633// _HAVE_SEMIC_
    … …  
    5744        this->parameters = NULL;
    5845        this->element_type_list=NULL;
     46        this->accumulator_values=NULL;
    5947}/*}}}*/
    6048Element::~Element(){/*{{{*/
    … …  
    7462        return false;
    7563}/*}}}*/
    76 void       Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type){/*{{{*/
     64void       Element::ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/
    7765   const int numvertices = this->GetNumberOfVertices();
    78         int         numperiods = numbreaks+1;
    79    int         basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type,indperiod;
    80    IssmDouble  time,dt,starttime,noiseterm;
     66   int         basinid,M,N,arenum_type,maenum_type,basinenum_type,noiseenum_type,outenum_type;
     67   IssmDouble  time,dt,starttime,termconstant_basin,trend_basin,noiseterm;
    8168   IssmDouble* arlagcoefs_basin     = xNew<IssmDouble>(arorder);
    8269   IssmDouble* malagcoefs_basin     = xNew<IssmDouble>(maorder);
    83    IssmDouble* datebreaks_basin     = xNew<IssmDouble>(numbreaks);
    84    IssmDouble* polyparams_basin     = xNew<IssmDouble>(numperiods*numparams);
    8570   IssmDouble* varlist              = xNew<IssmDouble>(numvertices);
    86    IssmDouble* sumpoly              = xNewZeroInit<IssmDouble>(arorder+1);
    8771   IssmDouble* valuesautoregression = NULL;
    8872   IssmDouble* valuesmovingaverage  = NULL;
    … …  
    11296         outenum_type   = BasalforcingsSpatialDeepwaterMeltingRateEnum;
    11397         break;
    114                 case(FrontalForcingsSubglacialDischargearmaEnum):
    115          arenum_type    = SubglacialdischargeValuesAutoregressionEnum;
    116          maenum_type    = SubglacialdischargeValuesMovingaverageEnum;
    117          basinenum_type = FrontalForcingsBasinIdEnum;
    118          noiseenum_type = SubglacialdischargeARMANoiseEnum;
    119          outenum_type   = FrontalForcingsSubglacialDischargeEnum;
    120          break;
    121                 case(HydrologyarmapwEnum):
    122          arenum_type    = WaterPressureValuesAutoregressionEnum;
    123          maenum_type    = WaterPressureValuesMovingaverageEnum;
    124          basinenum_type = HydrologyBasinsIdEnum;
    125          noiseenum_type = FrictionWaterPressureNoiseEnum;
    126          outenum_type   = WaterPressureArmaPerturbationEnum;
    127          break;
    128         }
     98   }
    12999
    130100        /*Get time parameters*/
    … …  
    135105   /*Get basin coefficients*/
    136106   this->GetInputValue(&basinid,basinenum_type);
    137         for(int i=0;i<arorder;i++) arlagcoefs_basin[i]   = arlagcoefs[basinid*arorder+i];
    138         for(int i=0;i<maorder;i++) malagcoefs_basin[i]   = malagcoefs[basinid*maorder+i];
    139         for(int i=0;i<numparams;i++){
    140                 for(int j=0;j<numperiods;j++) polyparams_basin[i*numperiods+j] = polyparams[basinid*numparams*numperiods+i*numperiods+j];
    141         }
    142         if(numbreaks>0){
    143                 for(int i=0;i<numbreaks;i++) datebreaks_basin[i] = datebreaks[basinid*numbreaks+i];
    144         }
    145 
    146         /*Compute terms from polynomial parameters from arorder steps back to present*/
    147         IssmDouble telapsed_break;
    148         IssmDouble tatstep;
    149         for(int s=0;s<arorder+1;s++){
    150                 tatstep = time-s*tstep_arma;
    151                 if(numbreaks>0){
    152                         /*Find index of tatstep compared to the breakpoints*/
    153                         indperiod = 0;
    154                         for(int i=0;i<numbreaks;i++){
    155                                 if(tatstep>=datebreaks_basin[i]) indperiod = i+1;
    156                         }
    157                         /*Compute polynomial with parameters of indperiod*/
    158                         if(indperiod==0) telapsed_break = tatstep-starttime;
    159                         else             telapsed_break = tatstep-datebreaks_basin[indperiod-1];
    160                         for(int j=0;j<numparams;j++)   sumpoly[s] = sumpoly[s]+polyparams_basin[indperiod+j*numperiods]*pow(telapsed_break,j);
    161                 }
    162                 else for(int j=0;j<numparams;j++) sumpoly[s] = sumpoly[s]+polyparams_basin[j*numperiods]*pow(tatstep-starttime,j);
    163         }
     107        for(int ii=0;ii<arorder;ii++) arlagcoefs_basin[ii] = arlagcoefs[basinid*arorder+ii];
     108        for(int ii=0;ii<maorder;ii++) malagcoefs_basin[ii] = malagcoefs[basinid*maorder+ii];
     109   termconstant_basin   = termconstant[basinid];
     110   trend_basin          = trend[basinid];
    164111
    165112        /*Initialze autoregressive and moving-average values at first time step*/
    … …  
    167114                IssmDouble* initvaluesautoregression = xNewZeroInit<IssmDouble>(numvertices*arorder);
    168115                IssmDouble* initvaluesmovingaverage  = xNewZeroInit<IssmDouble>(numvertices*maorder);
    169                 for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=polyparams_basin[0];
     116                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=termconstant_basin;
    170117      this->inputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder);
    171118      this->inputs->SetArrayInput(maenum_type,this->lid,initvaluesmovingaverage,numvertices*maorder);
    … …  
    195142         /*Compute autoregressive term*/
    196143         IssmDouble autoregressionterm=0.;
    197          for(int ii=0;ii<arorder;ii++) autoregressionterm += arlagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-sumpoly[ii+1]);
    198                         /*Compute moving-average term*/
     144         for(int ii=0;ii<arorder;ii++) autoregressionterm += arlagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-(termconstant_basin+trend_basin*(telapsed_arma-(ii+1)*tstep_arma)));
     145         /*Compute moving-average term*/
    199146         IssmDouble movingaverageterm=0.;
    200147         for(int ii=0;ii<maorder;ii++) movingaverageterm  += malagcoefs_basin[ii]*valuesmovingaverage[v+ii*numvertices];
    201148
    202149                        /*Stochastic variable value*/
    203          varlist[v] = sumpoly[0]+autoregressionterm+movingaverageterm+noiseterm;
    204      
    205                         /*Impose zero-bound*/
    206                         if(outenum_type == ThermalForcingEnum || outenum_type == FrontalForcingsSubglacialDischargeEnum) varlist[v] = max(varlist[v],0.0);
    207 
    208                 }
    209 
     150         varlist[v] = termconstant_basin+trend_basin*telapsed_arma+autoregressionterm+movingaverageterm+noiseterm;
     151      }
    210152      /*Update autoregression and moving-average values*/
    211153      IssmDouble* temparrayar = xNew<IssmDouble>(numvertices*arorder);
    … …  
    228170   xDelete<IssmDouble>(arlagcoefs_basin);
    229171   xDelete<IssmDouble>(malagcoefs_basin);
    230    xDelete<IssmDouble>(datebreaks_basin);
    231    xDelete<IssmDouble>(polyparams_basin);
    232    xDelete<IssmDouble>(sumpoly);
    233172   xDelete<IssmDouble>(varlist);
    234173   xDelete<IssmDouble>(valuesautoregression);
    235174   xDelete<IssmDouble>(valuesmovingaverage);
    236 }
    237 
    238 /*}}}*/
     175}/*}}}*/
     176void       Element::Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type){/*{{{*/
     177
     178   const int numvertices = this->GetNumberOfVertices();
     179   int         basinid,M,N,arenum_type,basinenum_type,noiseenum_type,outenum_type;
     180   IssmDouble  time,dt,starttime,termconstant_basin,trend_basin,noiseterm;
     181   IssmDouble* lagcoefs_basin  = xNew<IssmDouble>(arorder);
     182   IssmDouble* varlist         = xNew<IssmDouble>(numvertices);
     183   IssmDouble* valuesautoregression = NULL;
     184   Input*      noiseterm_input      = NULL;
     185
     186   /*Get field-specific enums*/
     187   switch(enum_type){
     188      case(SMBarmaEnum):
     189         arenum_type    = SmbValuesAutoregressionEnum;
     190         basinenum_type = SmbBasinsIdEnum;
     191         noiseenum_type = SmbARMANoiseEnum;
     192         outenum_type   = SmbMassBalanceEnum;
     193         break;
     194      case(FrontalForcingsRignotarmaEnum):
     195         arenum_type    = ThermalforcingValuesAutoregressionEnum;
     196         basinenum_type = FrontalForcingsBasinIdEnum;
     197         noiseenum_type = ThermalforcingARMANoiseEnum;
     198         outenum_type   = ThermalForcingEnum;
     199         break;
     200                case(BasalforcingsDeepwaterMeltingRatearmaEnum):
     201         arenum_type    = BasalforcingsDeepwaterMeltingRateValuesAutoregressionEnum;
     202         basinenum_type = BasalforcingsLinearBasinIdEnum;
     203         noiseenum_type = BasalforcingsDeepwaterMeltingRateNoiseEnum;
     204         outenum_type   = BasalforcingsSpatialDeepwaterMeltingRateEnum;
     205         break;
     206   }
     207
     208        /*Get time parameters*/
     209   this->parameters->FindParam(&time,TimeEnum);
     210   this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     211   this->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     212
     213   /*Get basin coefficients*/
     214   this->GetInputValue(&basinid,basinenum_type);
     215        for(int ii=0;ii<arorder;ii++) lagcoefs_basin[ii] = lagcoefs[basinid*arorder+ii];
     216   termconstant_basin   = termconstant[basinid];
     217   trend_basin   = trend[basinid];
     218
     219        /*Initialze autoregressive values at first time step*/
     220        if(time<=starttime+dt){
     221                IssmDouble* initvaluesautoregression = xNewZeroInit<IssmDouble>(numvertices*arorder);
     222                for(int i=0;i<numvertices*arorder;i++) initvaluesautoregression[i]=termconstant_basin;
     223      this->inputs->SetArrayInput(arenum_type,this->lid,initvaluesautoregression,numvertices*arorder);
     224      xDelete<IssmDouble>(initvaluesautoregression);
     225        }
     226
     227   /*Get noise and autoregressive terms*/
     228        if(isfieldstochastic){
     229      noiseterm_input = this->GetInput(noiseenum_type);
     230      Gauss* gauss = this->NewGauss();
     231      noiseterm_input->GetInputValue(&noiseterm,gauss);
     232      delete gauss;
     233   }
     234   else noiseterm = 0.0;
     235   this->inputs->GetArray(arenum_type,this->lid,&valuesautoregression,&M);
     236
     237        /*If not AR model timestep: take the old values of variable*/
     238   if(isstepforar==false){
     239      for(int i=0;i<numvertices;i++) varlist[i]=valuesautoregression[i];
     240   }
     241   /*If AR model timestep: compute variable values on vertices from AR*/
     242   else{
     243      for(int v=0;v<numvertices;v++){
     244
     245         /*Compute autoregressive term*/
     246         IssmDouble autoregressionterm=0.;
     247         for(int ii=0;ii<arorder;ii++) autoregressionterm += lagcoefs_basin[ii]*(valuesautoregression[v+ii*numvertices]-(termconstant_basin+trend_basin*(telapsed_ar-(ii+1)*tstep_ar)));
     248
     249         /*Stochastic variable value*/
     250         varlist[v] = termconstant_basin+trend_basin*telapsed_ar+autoregressionterm+noiseterm;
     251      }
     252      /*Update autoregression values*/
     253      IssmDouble* temparray = xNew<IssmDouble>(numvertices*arorder);
     254      /*Assign newest values and shift older values*/
     255      for(int i=0;i<numvertices;i++) temparray[i] = varlist[i];
     256      for(int i=0;i<(arorder-1)*numvertices;i++) temparray[i+numvertices] = valuesautoregression[i];
     257      this->inputs->SetArrayInput(arenum_type,this->lid,temparray,numvertices*arorder);
     258      xDelete<IssmDouble>(temparray);
     259   }
     260
     261   /*Add input to element*/
     262   this->AddInput(outenum_type,varlist,P1Enum);
     263
     264   /*Cleanup*/
     265   xDelete<IssmDouble>(lagcoefs_basin);
     266   xDelete<IssmDouble>(varlist);
     267   xDelete<IssmDouble>(valuesautoregression);
     268}/*}}}*/
    239269void       Element::BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation){/*{{{*/
    240270
    … …  
    264294         values[i]=deepwatermelt[i]*alpha+(1.-alpha)*upperwatermelt[basinid];
    265295      }
    266                 if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in melt");
    267                 if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in melt");
    268                 if(fabs(values[i])>1.e+10) _error_("melt exceeds 1.e+10");
    269296   }
    270297
    … …  
    944971
    945972        /*Get material parameters :*/
    946         IssmDouble rho_water = this->FindParam(MaterialsRhoFreshwaterEnum);
     973        IssmDouble rho_water = this->FindParam(MaterialsRhoSeawaterEnum);
    947974        IssmDouble rho_ice   = this->FindParam(MaterialsRhoIceEnum);
    948975
    … …  
    19031930                /*Are we in transient or static? */
    19041931                if(M==1){
    1905                         if(N!=1) _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
    1906                         _assert_(N==1);
     1932                        values[0]=vector[0];
    19071933                        this->SetElementInput(inputs,vector_enum,vector[0]);
    19081934                }
    1909 
    19101935                else if(M==iomodel->numberofvertices){
    1911                         if(N!=1) _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
    19121936                        for(int i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    19131937                        this->SetElementInput(inputs,NUM_VERTICES,vertexlids,values,vector_enum);
    … …  
    19481972                        }
    19491973                        else{
    1950                                 _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
     1974                                _error_("Patch interpolation not supported yet");
    19511975                        }
    19521976                        xDelete<IssmDouble>(evalues);
    … …  
    19541978                }
    19551979                else{
    1956                         _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
     1980                        _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    19571981                }
    19581982        }
    19591983        else if(vector_type==2){ //element vector
    19601984
     1985                IssmDouble value;
     1986
    19611987                /*Are we in transient or static? */
    1962                 if(M==1){
    1963                         if(N!=1) _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
    1964                         this->SetElementInput(inputs,vector_enum,vector[0]);
    1965                 }
    1966                 else if(M==2){
    1967                         /*create transient input: */
    1968                         IssmDouble* times = xNew<IssmDouble>(N);
    1969                         for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1970 
    1971                         inputs->SetTransientInput(vector_enum,times,N);
    1972                         TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
    1973 
    1974                         for(int t=0;t<N;t++){
    1975                                 IssmDouble value=vector[t]; //values are on the first line, times are on the second line
    1976                                 switch(this->ObjectEnum()){
    1977                                         case TriaEnum:  transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break;
    1978                                         case PentaEnum: transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum); break;
    1979                                         default: _error_("Not implemented yet");
    1980                                 }
    1981                         }
    1982                         xDelete<IssmDouble>(times);
    1983                 }
    1984                 else if(M==iomodel->numberofelements){
    1985                         if(N!=1) _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
     1988                if(M==iomodel->numberofelements){
    19861989                        if (code==5){ //boolean
    19871990                                this->SetBoolInput(inputs,vector_enum,reCast<bool>(vector[this->Sid()]));
    … …  
    20022005                        TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
    20032006                        for(int t=0;t<N;t++){
    2004                                 IssmDouble value=vector[N*this->Sid()+t];
     2007                                value=vector[N*this->Sid()+t];
    20052008                                switch(this->ObjectEnum()){
    20062009                                        case TriaEnum:  transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break;
    … …  
    20112014                        xDelete<IssmDouble>(times);
    20122015                }
    2013 
    2014                 else{
    2015                         _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
    2016                 }
     2016                else if(M==1 || M==2){
     2017                        /*create transient input: */
     2018                        IssmDouble* times = xNew<IssmDouble>(N);
     2019                        if(M==1)times[0]=0;
     2020                        if(M==2)for(int t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     2021
     2022                        inputs->SetTransientInput(vector_enum,times,N);
     2023                        TransientInput* transientinput = inputs->GetTransientInput(vector_enum);
     2024
     2025                        for(int t=0;t<N;t++){
     2026                                value=vector[t]; //values are on the first line, times are on the second line
     2027                                switch(this->ObjectEnum()){
     2028                                        case TriaEnum:  transientinput->AddTriaTimeInput( t,1,&(this->lid),&value,P0Enum); break;
     2029                                        case PentaEnum: transientinput->AddPentaTimeInput(t,1,&(this->lid),&value,P0Enum); break;
     2030                                        default: _error_("Not implemented yet");
     2031                                }
     2032                        }
     2033                        xDelete<IssmDouble>(times);
     2034                }
     2035                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
    20172036        }
    20182037        else if(vector_type==3){ //Double array matrix
    … …  
    20252044                        xDelete<IssmDouble>(layers);
    20262045                }
    2027                 else{
    2028                         _error_("Size of Input "<<EnumToStringx(vector_enum)<<" "<<M<<"x"<<N<<" not supported");
    2029                 }
    2030         }
    2031         else{
    2032                 _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    2033         }
     2046                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     2047        }
     2048        else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    20342049}
    20352050/*}}}*/
    … …  
    21592174        else _error_("not currently supported type of M and N attempted");
    21602175}/*}}}*/
    2161 void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int input_id){/*{{{*/
     2176void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/
    21622177        /*enum_type: the name of the DatasetInput (eg Outputdefinition1)
    21632178         * vector: information being stored (eg observations)
    21642179         * vector_type: is if by element or by vertex
    21652180         * input_enum: is the name of the vector being stored
     2181         * code: what type of data is in the vector (booleans, ints, doubles)
    21662182         */
    21672183
    … …  
    22222238                /*Are we in transient or static? */
    22232239                if(M==iomodel->numberofelements){
    2224                         _error_("not implemented");
     2240                        if (code==5){ //boolean
     2241                                _error_("not implemented");
     2242                                //datasetinput->AddInput(new BoolInput(input_enum,reCast<bool>(vector[this->Sid()])),input_id);
     2243                        }
     2244                        else if (code==6){ //integer
     2245                                _error_("not implemented");
     2246                                //datasetinput->AddInput(new IntInput(input_enum,reCast<int>(vector[this->Sid()])),input_id);
     2247                        }
     2248                        else if (code==7){ //IssmDouble
     2249                                _error_("not implemented");
     2250                                //datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id);
     2251                        }
     2252                        else _error_("could not recognize nature of vector from code " << code);
    22252253                }
    22262254                else if(M==iomodel->numberofelements+1){
    22272255                        _error_("not supported");
     2256                        ///*create transient input: */
     2257                        //IssmDouble* times = xNew<IssmDouble>(N);
     2258                        //for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     2259                        //TransientInput* transientinput=new TransientInput(input_enum,times,N);
     2260                        //TriaInput* bof=NULL;
     2261                        //for(t=0;t<N;t++){
     2262                        //      value=vector[N*this->Sid()+t];
     2263                        //      switch(this->ObjectEnum()){
     2264                        //              case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break;
     2265                        //              case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break;
     2266                        //              case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break;
     2267                        //              default: _error_("Not implemented yet");
     2268                        //      }
     2269                        //}
     2270                        //xDelete<IssmDouble>(times);
    22282271                }
    22292272                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    … …  
    23572400        Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
    23582401        return (input->GetInputMax()<=0.);
    2359 }
    2360 /*}}}*/
    2361 bool       Element::IsAllMinThicknessInElement(){/*{{{*/
    2362 
    2363         IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum);
    2364 
    2365         Input* input=this->GetInput(ThicknessEnum); _assert_(input);
    2366         if(input->GetInputMax()<=(minthickness+0.00000001)){
    2367                 return true;
    2368         }
    2369         else{
    2370                 return false;
    2371         }
    23722402}
    23732403/*}}}*/
    … …  
    24412471   const int numvertices = this->GetNumberOfVertices();
    24422472   bool isadjustsmb = false;
    2443         int basinid,bb1,bb2,mindex;
    2444         IssmDouble ela,refelevation_b,time,dt,fracyear,yts;
    2445    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};
     2473        int basinid,bb1,bb2;
     2474   IssmDouble ela,refelevation_b;
    24462475   IssmDouble* surfacelist  = xNew<IssmDouble>(numvertices);
    24472476   IssmDouble* smblist      = xNew<IssmDouble>(numvertices);
    2448    /* numelevbins values of lapse rates at current month */
     2477   /* numelevbins values of lapse rates */
    24492478        IssmDouble* lapserates_b = xNew<IssmDouble>(numelevbins);
    2450    /* (numelevbins-1) limits between elevation bins at current month (be cautious with indexing) */
     2479   /* (numelevbins-1) limits between elevation bins (be cautious with indexing) */
    24512480        IssmDouble* elevbins_b   = xNew<IssmDouble>(numelevbins-1);
    2452 
    2453         /*Find month of current time step*/
    2454         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    2455    this->parameters->FindParam(&time,TimeEnum);
    2456    this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    2457    fracyear     = time/yts-floor(time/yts);
    2458    for(int i=1;i<12;i++){
    2459                 if(fracyear>=monthsteps[i-1]) mindex = i-1;
    2460         }
    2461    if(fracyear>=monthsteps[11]) mindex = 11;
    24622481
    24632482   /*Retrieve SMB values non-adjusted for SMB lapse rate*/
    … …  
    24682487   this->GetInputValue(&basinid,SmbBasinsIdEnum);
    24692488   refelevation_b = refelevation[basinid];
    2470         /*Retrieve bins and laps rates for this basin at this month*/
    2471         for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)*12+mindex*(numelevbins-1)+ii];
     2489        for(int ii=0;ii<(numelevbins-1);ii++) elevbins_b[ii] = elevbins[basinid*(numelevbins-1)+ii];
    24722490        for(int ii=0;ii<numelevbins;ii++){
    2473                 lapserates_b[ii] = lapserates[basinid*numelevbins*12+mindex*numelevbins+ii];
     2491                lapserates_b[ii] = lapserates[basinid*numelevbins+ii];
    24742492                if(lapserates_b[ii]!=0) isadjustsmb=true;
    24752493        }
    2476        
    24772494        /*Adjust SMB if any lapse rate value is non-zero*/
    24782495        if(isadjustsmb){
    2479 
    2480                 _assert_(dt<yts);
     2496       
    24812497           for(int v=0;v<numvertices;v++){
    24822498              /*Find elevation bin of Reference elevation and of Vertex*/
    2483                         bb1 = 0;
    2484                         bb2 = 0;
    24852499                        for(int ii=0;ii<(numelevbins-1);ii++){
    2486                                 if(surfacelist[v]>=elevbins_b[ii]) bb1 = ii+1;
    2487                                 if(refelevation_b>=elevbins_b[ii]) bb2 = ii+1;
     2500                                if(surfacelist[v]<=elevbins_b[ii]) bb1 = ii;   
     2501                                if(refelevation_b<=elevbins_b[ii]) bb2 = ii;
    24882502                        }
    2489 
     2503                        /*Check for elevations above highest bin limit */
     2504                        if(surfacelist[v]>elevbins_b[numelevbins-1-1]) bb1 = numelevbins-1;
     2505                        if(refelevation_b>elevbins_b[numelevbins-1-1]) bb2 = numelevbins-1;
    24902506                        /*Vertex and Reference elevation in same elevation bin*/
    24912507                        if(bb1==bb2){
    … …  
    25222538        IssmDouble deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
    25232539        IssmDouble base[MAXVERTICES];
    2524         IssmDouble perturbation[MAXVERTICES];
    25252540        IssmDouble values[MAXVERTICES];
    25262541        IssmDouble time;
    … …  
    25342549
    25352550        this->GetInputListOnVertices(&base[0],BaseEnum);
    2536         this->GetInputListOnVertices(&perturbation[0],BasalforcingsPerturbationMeltingRateEnum);
    25372551        for(int i=0;i<NUM_VERTICES;i++){
    25382552                if(base[i]>=upperwaterel){
    … …  
    25462560                        values[i]=deepwatermelt*alpha+(1.-alpha)*upperwatermelt;
    25472561                }
    2548 
    2549                 /*Add perturbation*/
    2550                 values[i] += perturbation[i];
    25512562        }
    25522563
    … …  
    27662777
    27672778}/*}}}*/
    2768 void       Element::MonthlyFactorBasin(IssmDouble* monthlyfac,int enum_type){/*{{{*/
    2769        
    2770         /*Variable declaration*/
    2771         bool ratevariable;
    2772    const int numvertices = this->GetNumberOfVertices();
    2773         int basinid,mindex,mindexnext,basinenum_type,varenum_type,indperiod;
    2774    IssmDouble time,dt,fracyear,fracyearnext,fracmonth,fracmonthnext,yts;
    2775    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};
    2776    IssmDouble* monthlyfac_b   = xNew<IssmDouble>(12);
    2777    IssmDouble* monthlyrate_b  = xNew<IssmDouble>(12);
    2778         IssmDouble* fracdtinmonth  = xNew<IssmDouble>(12);
    2779         IssmDouble* rateinmonth    = xNew<IssmDouble>(numvertices*12);
    2780         IssmDouble* varlistinput   = xNew<IssmDouble>(numvertices);
    2781         IssmDouble* varlist        = xNewZeroInit<IssmDouble>(numvertices);
    2782 
    2783         /*Get field-specific enums*/
    2784    switch(enum_type){
    2785       case(FrontalForcingsSubglacialDischargearmaEnum):
    2786          basinenum_type = FrontalForcingsBasinIdEnum;
    2787          varenum_type   = FrontalForcingsSubglacialDischargeEnum;
    2788          ratevariable   = true;
    2789                         break;
    2790                 case(HydrologyarmapwEnum):
    2791          basinenum_type = HydrologyBasinsIdEnum;
    2792          varenum_type   = FrictionWaterPressureEnum;
    2793          ratevariable   = false;
    2794                         break;
    2795         }
    2796        
    2797         /*Evaluate the month index now and at (now-timestepjump)*/
    2798         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    2799         this->parameters->FindParam(&time,TimeEnum);
    2800    this->parameters->FindParam(&dt,TimesteppingTimeStepEnum); _assert_(dt<yts);
    2801         fracyear     = time/yts-floor(time/yts);
    2802         fracyearnext = (time+dt)/yts-floor((time+dt)/yts);
    2803         for(int i=1;i<12;i++){
    2804                 if(fracyear>=monthsteps[i-1])     mindex     = i-1;
    2805                 if(fracyearnext>=monthsteps[i-1]) mindexnext = i-1;
    2806         }
    2807         if(fracyear>=monthsteps[11])         mindex     = 11;
    2808         if(fracyearnext>=monthsteps[11])     mindexnext = 11;
    2809 
    2810         /*Calculate fraction of the time step spent in each month*/
    2811         for(int i=0;i<12;i++){
    2812                 if(mindex<i && mindexnext>i)                            fracdtinmonth[i] = 1.0/dt*yts/12.0;
    2813                 else if(mindex<i && mindexnext<i && mindexnext<mindex)  fracdtinmonth[i] = 1.0/dt*yts/12.0;
    2814                 else if(mindex>i && mindexnext<mindex && mindexnext>i)  fracdtinmonth[i] = 1.0/dt*yts/12.0;
    2815                 else if(mindex>i && mindexnext<mindex && mindexnext==i) fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-monthsteps[i]);
    2816                 else if(mindex==i && mindexnext==i)                     fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-fracyear);
    2817                 else if(mindex==i && mindexnext!=mindex)                fracdtinmonth[i] = 1.0/dt*yts*(1.0/12-(fracyear-monthsteps[i]));
    2818                 else if(mindexnext==i && mindex!=mindexnext)            fracdtinmonth[i] = 1.0/dt*yts*(fracyearnext-monthsteps[i]);
    2819                 else                                                      fracdtinmonth[i] = 0.0;
    2820         }
    2821 
    2822         /*Get basin-specific parameters of the element*/
    2823    this->GetInputValue(&basinid,basinenum_type);
    2824         for(int i=0;i<12;i++) monthlyfac_b[i]   = monthlyfac[basinid*12+i];
    2825 
    2826         /*Retrieve input*/
    2827         this->GetInputListOnVertices(varlistinput,varenum_type);
    2828 
    2829         /*Calculate monthly rate for each month and weight-average it for application over dt*/
    2830         for(int v=0;v<numvertices;v++){
    2831                 for(int i=0;i<12;i++){
    2832                         if(ratevariable){
    2833                                 rateinmonth[v*12+i] = varlistinput[v]*monthlyfac_b[i]*12;
    2834                                 varlist[v]          = varlist[v]+fracdtinmonth[i]*rateinmonth[v*12+i];
    2835                         }
    2836                         else varlist[v]       = varlist[v]+fracdtinmonth[i]*monthlyfac_b[i]*varlistinput[v];
    2837                 }
    2838         }
    2839         /*Update input*/
    2840    this->AddInput(varenum_type,varlist,P1DGEnum);
    2841 
    2842         /*Clean-up*/
    2843         xDelete<IssmDouble>(fracdtinmonth);
    2844         xDelete<IssmDouble>(rateinmonth);
    2845         xDelete<IssmDouble>(monthlyfac_b);
    2846         xDelete<IssmDouble>(monthlyrate_b);
    2847         xDelete<IssmDouble>(varlist);
    2848         xDelete<IssmDouble>(varlistinput);
    2849 }/*}}}*/
    2850 void       Element::MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type){/*{{{*/
     2779void       Element::MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type){/*{{{*/
    28512780       
    28522781        /*Variable declaration*/
    28532782   const int numvertices = this->GetNumberOfVertices();
    2854         int basinid,mindex,basinenum_type,varenum_type,indperiod;
    2855         int numperiods = nummonthbreaks+1;
     2783        int basinid,mindex,basinenum_type,varenum_type;
    28562784   IssmDouble time,fracyear,yts;
    28572785   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};
    2858    IssmDouble* datebreaks_b  = xNew<IssmDouble>(nummonthbreaks);
    2859         IssmDouble* intercepts_b  = xNew<IssmDouble>(numperiods*12);
    2860         IssmDouble* trends_b      = xNew<IssmDouble>(numperiods*12);
    2861         IssmDouble* varlist       = xNew<IssmDouble>(numvertices);
     2786        IssmDouble* monthlyeff_b = xNew<IssmDouble>(12);
     2787   IssmDouble* varlist      = xNew<IssmDouble>(numvertices);
    28622788
    28632789        /*Get field-specific enums*/
    … …  
    28802806        /*Get basin-specific parameters of the element*/
    28812807   this->GetInputValue(&basinid,basinenum_type);
    2882         if(nummonthbreaks>0){
    2883       for(int i=0;i<nummonthbreaks;i++) datebreaks_b[i] = monthlydatebreaks[basinid*nummonthbreaks+i];
    2884    }
    2885         for(int i=0;i<numperiods;i++){
    2886                 intercepts_b[i]  = monthlyintercepts[basinid*12*numperiods+mindex+12*i];
    2887                 trends_b[i]      = monthlytrends[basinid*12*numperiods+mindex+12*i];
    2888         }
    2889 
    2890         /*Compute piecewise-linear function*/
    2891         IssmDouble telapsed_break,piecewiselin;
    2892         if(nummonthbreaks>0){
    2893                 /*Find index of time compared to the breakpoints*/
    2894                 indperiod = 0;
    2895                 for(int i=0;i<nummonthbreaks;i++){
    2896                         if(time>datebreaks_b[i]) indperiod = i+1;
    2897                 }
    2898                 /*Compute intercept+trend terms with parameters of indperiod*/
    2899       if(indperiod==0) telapsed_break = time;
    2900       else             telapsed_break = time-datebreaks_b[indperiod-1];
    2901       piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*telapsed_break;
    2902         }
    2903    else piecewiselin = intercepts_b[indperiod]+trends_b[indperiod]*time;
     2808        for(int ii=0;ii<12;ii++){
     2809                monthlyeff_b[ii] = monthlyeff[basinid*12+ii];
     2810        }
    29042811
    29052812        /*Retrieve values non-adjusted for monthly effects*/
    … …  
    29072814
    29082815        /*Adjust values using monthly effects*/
    2909         for(int v=0;v<numvertices;v++) varlist[v] = varlist[v]+piecewiselin;
     2816        for(int v=0;v<numvertices;v++){
     2817                varlist[v] = varlist[v]+monthlyeff_b[mindex];
     2818        }
    29102819
    29112820        /*Add input to element*/
    … …  
    29132822
    29142823        /*Clean-up*/
    2915    xDelete<IssmDouble>(datebreaks_b);
    2916    xDelete<IssmDouble>(intercepts_b);
    2917    xDelete<IssmDouble>(trends_b);
     2824   xDelete<IssmDouble>(monthlyeff_b);
    29182825        xDelete<IssmDouble>(varlist);
    2919 }
    2920 /*}}}*/
     2826}/*}}}*/
    29212827void       Element::BeckmannGoosseFloatingiceMeltingRate(){/*{{{*/
    29222828
    … …  
    37413647}
    37423648/*}}}*/
    3743 void       Element::SmbDebrisEvatt(){/*{{{*/
    3744 
    3745         const int NUM_VERTICES          = this->GetNumberOfVertices();
    3746         const int NUM_VERTICES_DAYS_PER_YEAR  = NUM_VERTICES * 365; // 365 FIXME
    3747 
    3748         int             i,vertexlids[MAXVERTICES];;
    3749         IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);
    3750         IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);
    3751         IssmDouble* summermelt=xNew<IssmDouble>(NUM_VERTICES);
    3752         IssmDouble* albedo=xNew<IssmDouble>(NUM_VERTICES);
    3753         IssmDouble* summeralbedo=xNew<IssmDouble>(NUM_VERTICES);
    3754         IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES);
    3755        
    3756         // climate inputs
    3757         IssmDouble* temperature=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3758         IssmDouble* precip=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3759         IssmDouble* lw=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3760         IssmDouble* sw=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3761         IssmDouble* wind=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3762         IssmDouble* humidity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    3763         IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
    3764         IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3765         IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3766         IssmDouble* lw_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3767         IssmDouble* sw_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3768         IssmDouble* wind_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3769         IssmDouble* humidity_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3770 
    3771         IssmDouble* surface=xNew<IssmDouble>(NUM_VERTICES);
    3772         IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);
    3773         IssmDouble* snowheight=xNew<IssmDouble>(NUM_VERTICES);
    3774         IssmDouble* debriscover=xNew<IssmDouble>(NUM_VERTICES);
    3775         IssmDouble rho_water,rho_ice,Tf,debris,debris_here;
    3776         IssmDouble qlaps,rlaps,dsgrad,dlgrad,windspeedgrad,humiditygrad,Tm;
    3777         IssmDouble inv_twelve=1./365.;
    3778         IssmDouble time,yts,time_yr,lambda;
    3779         IssmDouble DailyMelt,CleanIceDailyMelt, CumDailyMelt=0,CleanIceMelt,CumDailySummerMelt=0;
    3780         IssmDouble MeanAlbedo=0, MeanSummerAlbedo=0;
    3781         bool isdebris,isAnderson,iscryokarst;
    3782         this->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
    3783         this->parameters->FindParam(&isAnderson,SmbDebrisIsAndersonEnum);
    3784         this->parameters->FindParam(&iscryokarst,SmbDebrisIsCryokarstEnum);
    3785         IssmDouble PhiD=0.,p;
    3786         IssmDouble icealbedo=this->FindParam(SmbIcealbedoEnum);
    3787         IssmDouble snowalbedo=this->FindParam(SmbSnowalbedoEnum);
    3788         IssmDouble debrisalbedo=this->FindParam(SmbDebrisalbedoEnum);
    3789         IssmDouble Lm=this->FindParam(MaterialsLatentheatEnum);
    3790         IssmDouble D0=this->FindParam(SmbDebrisAndersonD0Enum);
    3791         int step;
    3792         this->FindParam(&step,StepEnum);
    3793 
    3794         // cryokarst
    3795         int dim=1,domaintype;
    3796         this->parameters->FindParam(&domaintype,DomainTypeEnum);
    3797         if(domaintype!=Domain2DverticalEnum){
    3798                         dim=2;
    3799         }
    3800         IssmDouble taud_plus=110e3, taud_minus=60e3;
    3801         IssmDouble taud, slope, gravity, taudx, taudy;
    3802         this->parameters->FindParam(&gravity,ConstantsGEnum);
    3803         IssmDouble* slopex         = xNew<IssmDouble>(NUM_VERTICES);
    3804         IssmDouble* slopey         = xNew<IssmDouble>(NUM_VERTICES);
    3805         IssmDouble* icethickness   = xNew<IssmDouble>(NUM_VERTICES);
    3806 
    3807         /*Get material parameters :*/
    3808         rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
    3809         rho_ice=this->FindParam(MaterialsRhoIceEnum);
    3810         IssmDouble sconv=(rho_water/rho_ice);
    3811         Tf=this->FindParam(MaterialsMeltingpointEnum);
    3812 
    3813         /*Get parameters for height corrections*/
    3814         qlaps=this->FindParam(SmbDesfacEnum); // comment MR; on alpine galciers we dont have the desertification effect
    3815         rlaps=this->FindParam(SmbRlapsEnum);
    3816         dsgrad=this->FindParam(SmbSWgradEnum);
    3817         dlgrad=this->FindParam(SmbLWgradEnum);
    3818         windspeedgrad=this->FindParam(SmbWindspeedgradEnum);
    3819         humiditygrad=this->FindParam(SmbHumiditygradEnum);
    3820 
    3821         /* Get time */
    3822         this->parameters->FindParam(&time,TimeEnum);
    3823         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    3824         time_yr=floor(time/yts)*yts;
    3825 
    3826         /*Get inputs*/
    3827         DatasetInput* tempday     =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(tempday);
    3828         DatasetInput* precipday   =this->GetDatasetInput(SmbPrecipitationEnum);       _assert_(precipday);
    3829         DatasetInput* lwday       =this->GetDatasetInput(SmbMonthlydlradiationEnum); _assert_(lwday);
    3830         DatasetInput* swday       =this->GetDatasetInput(SmbMonthlydsradiationEnum);       _assert_(swday);
    3831         DatasetInput* windday     =this->GetDatasetInput(SmbMonthlywindspeedEnum); _assert_(windday);
    3832         DatasetInput* humidityday =this->GetDatasetInput(SmbMonthlyairhumidityEnum); _assert_(humidityday);
    3833 
    3834         /*loop over vertices: */
    3835         Gauss* gauss=this->NewGauss();
    3836         for(int month=0;month<365;month++){
    3837                 for(int iv=0;iv<NUM_VERTICES;iv++){
    3838                         gauss->GaussVertex(iv);
    3839                         tempday->GetInputValue(&temperature[iv*365+month],gauss,month);
    3840                         temperature[iv*365+month]=temperature[iv*365+month]-Tf; // conversion from Kelvin to celcius for PDD module
    3841                         precipday->GetInputValue(&precip[iv*365+month],gauss,month);
    3842                         precip[iv*365+month]=precip[iv*365+month]*yts; // from m/s to m/a
    3843                         lwday->GetInputValue(&lw[iv*365+month],gauss,month);
    3844                         swday->GetInputValue(&sw[iv*365+month],gauss,month);
    3845                         windday->GetInputValue(&wind[iv*365+month],gauss,month);
    3846                         humidityday->GetInputValue(&humidity[iv*365+month],gauss,month);
    3847                 }
    3848         }
    3849 
    3850         /*Recover info at the vertices: */
    3851         GetInputListOnVertices(&surface[0],SurfaceEnum);
    3852         GetInputListOnVertices(&s0t[0],SmbS0tEnum);
    3853         GetInputListOnVertices(&snowheight[0],SmbSnowheightEnum);
    3854         GetInputListOnVertices(&debriscover[0],DebrisThicknessEnum);
    3855         GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum);
    3856         GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum);
    3857         GetInputListOnVertices(&lw_ampl[0],SmbDsradiationAnomalyEnum);
    3858         GetInputListOnVertices(&sw_ampl[0],SmbDlradiationAnomalyEnum);
    3859         GetInputListOnVertices(&wind_ampl[0],SmbWindspeedAnomalyEnum);
    3860         GetInputListOnVertices(&humidity_ampl[0],SmbAirhumidityAnomalyEnum);
    3861         if(iscryokarst){
    3862                 GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
    3863                 GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    3864                 if(dim==2){
    3865                         GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum);
    3866                 }
    3867                 taudx=rho_ice*gravity*icethickness[i]*slopex[i];
    3868                 if(dim==2) taudy=rho_ice*gravity*icethickness[i]*slopey[i];
    3869                 taud=sqrt(taudx*taudx+taudy*taudy);
    3870         }
    3871         IssmDouble Alphaeff,Alphaeff_cleanice;
    3872 
    3873         /*measure the surface mass balance*/
    3874         for (int iv = 0; iv<NUM_VERTICES; iv++){
    3875 
    3876                 IssmDouble st=(surface[iv]-s0t[iv])/1000.;
    3877 
    3878                 int ismb_end=1;
    3879                 if(isdebris & !isAnderson) ismb_end=2;
    3880                 for (int ismb=0;ismb<ismb_end;ismb++){
    3881                         if(ismb==0){
    3882                                 // calc a reference smb to identify accum and melt region; debris only develops in ablation area
    3883                                 debris=0.;
    3884                                 PhiD=0.;
    3885                                 if(isAnderson) debris_here=debriscover[iv]; // store debris for later
    3886                         }else{
    3887                                 // debris only develops in ablation area
    3888                                 /*if((accu[iv]/yts-CleanIceMelt)<(-1e-2)/yts){
    3889                                         debris=debriscover[iv];
    3890                                 }else{
    3891                                         debris=0.;
    3892                                 }*/
    3893                                 debris=0.;
    3894                                 if(debris<=0.) debris=0.;
    3895                                 if(isdebris) PhiD=FindParam(DebrisPackingFractionEnum);
    3896                                 CumDailyMelt=0;
    3897                                 CumDailySummerMelt=0;
    3898                                 debris_here=debriscover[iv];
    3899                         }
    3900 
    3901                         /* Now run the debris part */
    3902 
    3903                         // Climate inputs
    3904                         IssmDouble Tm;          // C air temperature
    3905                         IssmDouble In;          // Wm^-2 incoming long wave
    3906                         IssmDouble Q;           // Wm^-2 incoming short wave
    3907                         IssmDouble Um;          // ms^-1 measured wind speed
    3908                         IssmDouble Humidity;    // relative humidity
    3909                         IssmDouble P;           // precip
    3910 
    3911                         // other parameters
    3912                         IssmDouble Qh=0.006;   // kg m^-3      saturated humidity level // not used
    3913                         IssmDouble Qm=0.8*Qh;  // kg m^-3      measured humiditiy level // not used
    3914                         IssmDouble Rhoaa=1.22; // kgm^-3       air densitiy
    3915                         IssmDouble K=0.585;    // Wm^-1K^-1    thermal conductivity          0.585
    3916                         IssmDouble Xr=0.01;    // ms^-1        surface roughness             0.01
    3917                         IssmDouble Ustar=0.16; // ms^-1        friction velocity             0.16
    3918                         IssmDouble Ca=1000;    // jkg^-1K^-1   specific heat capacity of air
    3919                         IssmDouble Lv=2.50E+06;// jkg^-1K^-1   latent heat of evaporation
    3920                         IssmDouble Eps=0.95;   //              thermal emissivity
    3921                         IssmDouble Sigma=5.67E-08;// Wm^-2K^-4    Stefan Boltzmann constant
    3922                         IssmDouble Gamma=180.;    // m^-1         wind speed attenuation        234
    3923                
    3924                         // Calculate effective albedo
    3925                         IssmDouble Alphaeff,Alphaeff_cleanice;
    3926                         IssmDouble mean_ela,delta=2000;
    3927                        
    3928                         // compute cleanice albedo based on previous SMB distribution
    3929                         //if(step==1){
    3930                                 mean_ela=3000; //FIXME
    3931                         //}else{
    3932                         //        mean_ela=FindParam(SmbMeanElaEnum);
    3933                         //}
    3934                         Alphaeff_cleanice=icealbedo+(snowalbedo-icealbedo)*(1+tanh(PI*(surface[iv]-mean_ela)/delta))/2;
    3935                         Alphaeff=Alphaeff_cleanice; // will be updated below
    3936 
    3937                        
    3938                         accu[iv]=0.;
    3939                         for (int iday=0;iday<365;iday++) {
    3940 
    3941                                 Tm=temperature[iv*365+iday]-st*rlaps;//+t_ampl[iv];//+(rand()%10-5)/5;
    3942                                 In=lw[iv*365+iday]-st*dlgrad+lw_ampl[iv];
    3943                                 Q=sw[iv*365+iday]+st*dsgrad+sw_ampl[iv];
    3944                                 Humidity=humidity[iv*365+iday]-st*humiditygrad+humidity_ampl[iv];
    3945                                 Um=wind[iv*365+iday]-st*windspeedgrad+wind_ampl[iv];
    3946                                 P=(qlaps*st*precip[iv*365+iday]+precip[iv*365+iday]+p_ampl[iv])*sconv/365.; // convert precip from w.e. -> i.e
    3947 
    3948                                 /*Partition of precip in solid and liquid parts */
    3949                                 IssmDouble temp_plus=1;
    3950                                 IssmDouble temp_minus=-1.;
    3951                                 IssmDouble frac_solid;
    3952                                 if(Tm>=temp_plus){
    3953                                         frac_solid=0;
    3954                                 }else if(Tm<=temp_minus){
    3955                                         frac_solid=1;
    3956                                 }else{
    3957                                         frac_solid=1*(1-cos(PI*(temp_plus-Tm)/(temp_plus-temp_minus)))/2;
    3958                                 }
    3959 
    3960                                 /*Get yearly temperatures and accumulation */
    3961                                 yearlytemperatures[iv]=yearlytemperatures[iv]+((temperature[iv*365+iday]-rlaps*st+Tf+t_ampl[iv]))/365; // Has to be in Kelvin
    3962                                 accu[iv]=accu[iv]+P*frac_solid;
    3963                                 if(yearlytemperatures[iv]>Tf) yearlytemperatures[iv]=Tf;
    3964 
    3965                                 CleanIceDailyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    3966                                         Q*(1.-Alphaeff)+
    3967                                         (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+
    3968                                         ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    3969                                         K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    3970                                         rho_ice*Lm*Ustar)/(((Um
    3971                                         -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)));
    3972                                 if(CleanIceDailyMelt<0) CleanIceDailyMelt=0.;
    3973                                 DailyMelt=CleanIceDailyMelt;
    3974 
    3975                                 if(ismb==1){
    3976 
    3977                                         //snowheight[iv]=snowheight[iv]+(P-CleanIceDailyMelt*yts/365);
    3978                                         IssmDouble sn_prev;
    3979                                         sn_prev=snowheight[iv];
    3980                                         snowheight[iv]=sn_prev+(-CleanIceDailyMelt*yts/365);//P
    3981                                        
    3982                                         if(snowheight[iv]<=0) snowheight[iv]=0.;
    3983                                         if(snowheight[iv]<=0.0001){
    3984                                                 p=debris_here*PhiD/(2*0.2*0.01); //Eq. 51 from Evatt et al 2015 without source term g*t
    3985                                                 if(p>1.) p=1.;
    3986                                                 if(p>=0.999){
    3987                                                         Alphaeff=debrisalbedo;
    3988                                                 } else {
    3989                                                         Alphaeff=Alphaeff_cleanice+p*(debrisalbedo-Alphaeff_cleanice);
    3990                                                 }
    3991                                                 debris=debris_here;
    3992                                                 DailyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    3993                                                         Q*(1.-Alphaeff)+
    3994                                                         (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+
    3995                                                         ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    3996                                                         K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    3997                                                         rho_ice*Lm*Ustar)/(((Um-2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)));
    3998                                                 if(DailyMelt<0) DailyMelt=0.;
    3999                                                 MeanSummerAlbedo=MeanSummerAlbedo+Alphaeff;
    4000                                                 CumDailySummerMelt=CumDailySummerMelt+DailyMelt/365;
    4001                                         }
    4002                                 }
    4003                                 CumDailyMelt=CumDailyMelt+DailyMelt/365;
    4004                         }
    4005                         MeanAlbedo=MeanAlbedo+Alphaeff;
    4006                         if(ismb==0) CleanIceMelt=CumDailyMelt;
    4007                 }
    4008 
    4009                 if(iscryokarst){
    4010                         if(taud>=taud_plus){
    4011                                 lambda=0;
    4012                         }else if(taud>=taud_minus & taud<taud_plus){
    4013                                 lambda=0.1*(1-cos(PI*(taud_plus-taud)/(taud_plus-taud_minus)))/2;
    4014                         }else if(taud<taud_minus){
    4015                                 lambda=0.1;
    4016                         }
    4017                 }
    4018 
    4019                 // update values
    4020                 melt[iv]=CumDailyMelt; // is already in m/s
    4021                 accu[iv]=accu[iv]/yts;
    4022                 if(isAnderson){
    4023                         smb[iv]=(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4024                         if(iscryokarst){
    4025                                 smb[iv]=lambda*(accu[iv]-melt[iv])+(1-lambda)*(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4026                         }else{
    4027                                 smb[iv]=(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4028                         }
    4029                 }else{
    4030                         if(iscryokarst){
    4031                                 smb[iv]=lambda*(accu[iv]-CleanIceMelt)+(1-lambda)*(accu[iv]-melt[iv]);
    4032                         }else{
    4033                                 smb[iv]=(accu[iv]-melt[iv]);
    4034                         }
    4035                 }
    4036                 albedo[iv]=MeanAlbedo;
    4037                 summeralbedo[iv]=MeanSummerAlbedo;
    4038                 summermelt[iv]=CumDailySummerMelt;
    4039         }
    4040 
    4041         this->AddInput(SmbMassBalanceEnum,smb,P1Enum);
    4042         this->AddInput(SmbAccumulationEnum,accu,P1Enum);
    4043         this->AddInput(SmbMeltEnum,melt,P1Enum);
    4044         this->AddInput(SmbSummerMeltEnum,summermelt,P1Enum);
    4045         this->AddInput(SmbSnowheightEnum,snowheight,P1Enum);
    4046         this->AddInput(SmbAlbedoEnum,albedo,P1Enum);
    4047         this->AddInput(SmbSummerAlbedoEnum,summeralbedo,P1Enum);
    4048         this->AddInput(TemperaturePDDEnum,yearlytemperatures,P1Enum); // TemperaturePDD is wrong here, but don't want to create new Enum ...
    4049 
    4050         /*clean-up*/
    4051         xDelete<IssmDouble>(temperature);
    4052         xDelete<IssmDouble>(precip);
    4053         xDelete<IssmDouble>(lw);
    4054         xDelete<IssmDouble>(sw);
    4055         xDelete<IssmDouble>(wind);
    4056         xDelete<IssmDouble>(humidity);
    4057         xDelete<IssmDouble>(smb);
    4058         xDelete<IssmDouble>(surface);
    4059         xDelete<IssmDouble>(melt);
    4060         xDelete<IssmDouble>(summermelt);
    4061         xDelete<IssmDouble>(albedo);
    4062         xDelete<IssmDouble>(summeralbedo);
    4063         xDelete<IssmDouble>(accu);
    4064         xDelete<IssmDouble>(yearlytemperatures);
    4065         xDelete<IssmDouble>(s0t);
    4066         xDelete<IssmDouble>(snowheight);
    4067         xDelete<IssmDouble>(debriscover);
    4068         xDelete<IssmDouble>(t_ampl);
    4069         xDelete<IssmDouble>(p_ampl);
    4070         xDelete<IssmDouble>(lw_ampl);
    4071         xDelete<IssmDouble>(sw_ampl);
    4072         xDelete<IssmDouble>(humidity_ampl);
    4073         xDelete<IssmDouble>(wind_ampl);
    4074         xDelete<IssmDouble>(slopex);
    4075         xDelete<IssmDouble>(slopey);
    4076         xDelete<IssmDouble>(icethickness);
    4077 }
    4078 /*}}}*/
    4079 
    4080 
    4081 
    40823649void       Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/
    40833650
    … …  
    41363703                                                                                                                  this->CalvingRateVonmises();
    41373704                                                                                                                  break;
    4138                                                                                                           case CalvingVonmisesADEnum:
    4139                                                                                                                   this->CalvingRateVonmisesAD();
    4140                                                                                                                   break;
    41413705                                                                                                          case CalvingCrevasseDepthEnum:
    41423706                                                                                                                  this->CalvingCrevasseDepth();
    … …  
    41443708                                                                                                          case CalvingParameterizationEnum:
    41453709                                                                                                                  this->CalvingRateParameterization();
    4146                                                                                                                   break;
    4147                                                                                                           case CalvingCalvingMIPEnum:
    4148                                                                                                                   this->CalvingRateCalvingMIP();
    41493710                                                                                                                  break;
    41503711                                                                                                          case CalvingTestEnum:
    … …  
    43133874
    43143875   /* Start looping on the number of vertices: */
    4315         Gauss* gauss=this->NewGauss();
     3876   GaussTria gauss;
    43163877   for(int iv=0;iv<numvertices;iv++){
    4317                 gauss->GaussVertex(iv);
     3878      gauss.GaussVertex(iv);
    43183879
    43193880      /* Get variables */
    4320       bed_input->GetInputValue(&bed,gauss);
    4321       qsg_input->GetInputValue(&qsg,gauss);
    4322                 TF_input->GetInputValue(&TF,gauss);
     3881      bed_input->GetInputValue(&bed,&gauss);
     3882      qsg_input->GetInputValue(&qsg,&gauss);
     3883      TF_input->GetInputValue(&TF,&gauss);
    43233884
    43243885      if(basin_icefront_area[basinid]==0.) meltrates[iv]=0.;
    … …  
    43293890         /* calculate melt rates */
    43303891         meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
    4331                 }
     3892      }
    43323893
    43333894      if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
    … …  
    43393900
    43403901   /*Cleanup and return*/
    4341    delete gauss;
    4342         xDelete<IssmDouble>(basin_icefront_area);
     3902   xDelete<IssmDouble>(basin_icefront_area);
    43433903}
    43443904/*}}}*/
    … …  
    45594119        xDelete<IssmDouble>(st);
    45604120        xDelete<IssmDouble>(s0gcm);
    4561 }
    4562 /*}}}*/
    4563 void       Element::SmbSemicTransient(){/*{{{*/
    4564 
    4565         bool isverbose=VerboseSmb();
    4566         if(isverbose && this->Sid()==0){
    4567                 _printf0_("smb core: initialize.\n");
    4568         }
    4569         /*only compute SMB at the surface: */
    4570         if (!IsOnSurface()) return;
    4571 
    4572         const int NUM_VERTICES                                  = this->GetNumberOfVertices();
    4573 
    4574         // daily forcing inputs
    4575         IssmDouble* dailyrainfall   =xNew<IssmDouble>(NUM_VERTICES);
    4576         IssmDouble* dailysnowfall   =xNew<IssmDouble>(NUM_VERTICES);
    4577         IssmDouble* dailydlradiation=xNew<IssmDouble>(NUM_VERTICES);
    4578         IssmDouble* dailydsradiation=xNew<IssmDouble>(NUM_VERTICES);
    4579         IssmDouble* dailywindspeed  =xNew<IssmDouble>(NUM_VERTICES);
    4580         IssmDouble* dailypressure   =xNew<IssmDouble>(NUM_VERTICES);
    4581         IssmDouble* dailyairdensity =xNew<IssmDouble>(NUM_VERTICES);
    4582         IssmDouble* dailyairhumidity=xNew<IssmDouble>(NUM_VERTICES);
    4583         IssmDouble* dailytemperature=xNew<IssmDouble>(NUM_VERTICES);
    4584 
    4585         // inputs: geometry
    4586         IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
    4587         IssmDouble* s0gcm=xNew<IssmDouble>(NUM_VERTICES);
    4588         IssmDouble* st=xNew<IssmDouble>(NUM_VERTICES);
    4589 
    4590         // inputs
    4591         IssmDouble* tsurf_in        =xNew<IssmDouble>(NUM_VERTICES);
    4592         IssmDouble* mask_in         =xNew<IssmDouble>(NUM_VERTICES);
    4593         IssmDouble* Tamp_in         =xNew<IssmDouble>(NUM_VERTICES);
    4594         IssmDouble* albedo_in       =xNew<IssmDouble>(NUM_VERTICES);
    4595         IssmDouble* albedo_snow_in  =xNew<IssmDouble>(NUM_VERTICES);
    4596         IssmDouble* hice_in         =xNew<IssmDouble>(NUM_VERTICES);
    4597         IssmDouble* hsnow_in        =xNew<IssmDouble>(NUM_VERTICES);
    4598         IssmDouble* qmr_in          =xNew<IssmDouble>(NUM_VERTICES);
    4599 
    4600         // outputs
    4601         IssmDouble* tsurf_out  =xNew<IssmDouble>(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4602         IssmDouble* smb_out    =xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4603         IssmDouble* smbi_out   =xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4604         IssmDouble* smbs_out   =xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4605         IssmDouble* saccu_out  =xNew<IssmDouble>(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4606         IssmDouble* smelt_out  =xNew<IssmDouble>(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4607         IssmDouble* refr_out  =xNew<IssmDouble>(NUM_VERTICES); memset(refr_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4608         IssmDouble* albedo_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4609         IssmDouble* albedo_snow_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_snow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4610         IssmDouble* hsnow_out   =xNew<IssmDouble>(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4611         IssmDouble* hice_out    =xNew<IssmDouble>(NUM_VERTICES); memset(hice_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4612         IssmDouble* qmr_out     =xNew<IssmDouble>(NUM_VERTICES); memset(qmr_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    4613 
    4614         IssmDouble rho_water,rho_ice,desfac,desfacElev,rlaps,rdl;
    4615         IssmDouble alb_smax, alb_smin, albi, albl;
    4616         IssmDouble hcrit, rcrit; // parameters for ? and refreezing.
    4617         int alb_scheme;
    4618         // albedo parameters - slatter
    4619         IssmDouble tmin, tmax;
    4620         // albedo parameters - isba
    4621         IssmDouble mcrit, tau_a, tau_f, wcrit;
    4622         // albedo parameters - alex
    4623         IssmDouble tmid, afac;
    4624 
    4625         IssmDouble tstart, time,yts,time_yr,dt;
    4626 
    4627         /* Get time: */
    4628         this->parameters->FindParam(&time,TimeEnum);
    4629         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4630         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    4631         this->parameters->FindParam(&tstart,TimesteppingStartTimeEnum);
    4632         time_yr=floor(time/yts)*yts;
    4633         //dt = dt * yts;
    4634 
    4635         /*Get material parameters :*/
    4636         rho_water=this->FindParam(MaterialsRhoFreshwaterEnum);
    4637         rho_ice=this->FindParam(MaterialsRhoIceEnum);
    4638         desfac=this->FindParam(SmbDesfacEnum);
    4639         desfacElev=this->FindParam(SmbDesfacElevEnum);
    4640         rlaps=this->FindParam(SmbRlapsEnum);
    4641         rdl=this->FindParam(SmbRdlEnum);
    4642 
    4643         this->FindParam(&alb_scheme,SmbAlbedoSchemeEnum);
    4644         this->FindParam(&hcrit,SmbSemicHcritEnum);
    4645         this->FindParam(&rcrit,SmbSemicRcritEnum);
    4646         alb_smax=this->FindParam(SmbAlbedoSnowMaxEnum);
    4647         alb_smin=this->FindParam(SmbAlbedoSnowMinEnum);
    4648         albi=this->FindParam(SmbAlbedoIceEnum);
    4649         albl=this->FindParam(SmbAlbedoLandEnum);
    4650 
    4651         // albedo parameters
    4652         this->FindParam(&tmid,SmbSemicTmidEnum);
    4653         this->FindParam(&tmin,SmbSemicTminEnum);
    4654         this->FindParam(&tmax,SmbSemicTmaxEnum);
    4655         this->FindParam(&mcrit,SmbSemicMcritEnum);
    4656         this->FindParam(&wcrit,SmbSemicWcritEnum);
    4657         this->FindParam(&tau_a,SmbSemicTauAEnum);
    4658         this->FindParam(&tau_f,SmbSemicTauFEnum);
    4659         this->FindParam(&afac,SmbSemicAfacEnum);
    4660 
    4661         /* Recover info at the vertices: */
    4662         GetInputListOnVertices(&s[0],SurfaceEnum);
    4663         GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum);
    4664 
    4665         if(isverbose && this->Sid()==0){
    4666                 _printf0_("smb core: allocate inputs.\n");
    4667                 _printf0_("smb core: time_yr  : " << time_yr/yts <<"\n");
    4668                 _printf0_("smb core: time     : " << time <<"\n");
    4669                 _printf0_("smb core: dt       : " << dt <<"\n");
    4670         }
    4671         /* loop over vertices and days */
    4672         Gauss* gauss=this->NewGauss();
    4673         /* Retrieve inputs: */
    4674         Input* dailysnowfall_input    = this->GetInput(SmbDailysnowfallEnum,time); _assert_(dailysnowfall_input);
    4675         Input* dailyrainfall_input    = this->GetInput(SmbDailyrainfallEnum,time); _assert_(dailyrainfall_input);
    4676         Input* dailydlradiation_input = this->GetInput(SmbDailydlradiationEnum,time); _assert_(dailydlradiation_input);
    4677         Input* dailydsradiation_input = this->GetInput(SmbDailydsradiationEnum,time); _assert_(dailydsradiation_input);
    4678         Input* dailywindspeed_input   = this->GetInput(SmbDailywindspeedEnum,time); _assert_(dailywindspeed_input);
    4679         Input* dailypressure_input    = this->GetInput(SmbDailypressureEnum,time); _assert_(dailypressure_input);
    4680         Input* dailyairdensity_input  = this->GetInput(SmbDailyairdensityEnum,time); _assert_(dailyairdensity_input);
    4681         Input* dailyairhumidity_input = this->GetInput(SmbDailyairhumidityEnum,time); _assert_(dailyairhumidity_input);
    4682         Input* dailytemperature_input = this->GetInput(SmbDailytemperatureEnum,time); _assert_(dailytemperature_input);
    4683 
    4684         /*temporal Enum depending on time*/
    4685         int enum_temp       =TemperatureSEMICEnum;
    4686         int enum_hice       =SmbHIceEnum;
    4687         int enum_hsnow      =SmbHSnowEnum;
    4688         int enum_albedo     =SmbAlbedoEnum;
    4689         int enum_albedo_snow=SmbAlbedoSnowEnum;
    4690         int enum_qmr        =SmbSemicQmrEnum;
    4691         if (tstart+dt == time) {
    4692                 /* Load inital value at first time step*/
    4693                 enum_temp=TemperatureEnum;
    4694                 enum_hice=SmbHIceInitEnum;
    4695                 enum_hsnow=SmbHSnowInitEnum;
    4696                 enum_albedo=SmbAlbedoInitEnum;
    4697                 enum_albedo_snow=SmbAlbedoSnowInitEnum;
    4698                 enum_qmr        =SmbSemicQmrInitEnum;
    4699         }
    4700         //if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n");
    4701         Input* tsurf_input       = this->GetInput(enum_temp); _assert_(tsurf_in);
    4702         //if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n");
    4703         Input* mask_input        = this->GetInput(SmbMaskEnum); _assert_(mask_input);
    4704         //if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n");
    4705         Input* Tamp_input        = this->GetInput(SmbTampEnum); _assert_(Tamp_input);
    4706         //if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n");
    4707         Input* albedo_input      = this->GetInput(enum_albedo); _assert_(albedo_input);
    4708         Input* albedo_snow_input = this->GetInput(enum_albedo_snow); _assert_(albedo_snow_input);
    4709         Input* hice_input        = this->GetInput(enum_hice); _assert_(hice_input);
    4710         Input* hsnow_input       = this->GetInput(enum_hsnow); _assert_(hsnow_input);
    4711         Input* qmr_input         = this->GetInput(enum_qmr); _assert_(qmr_input);
    4712 
    4713         if(isverbose && this->Sid()==0)_printf0_("smb core: assign inputs done....\n");
    4714         for(int iv=0;iv<NUM_VERTICES;iv++){
    4715                 gauss->GaussVertex(iv);
    4716                 /* get forcing */
    4717                 dailyrainfall_input->GetInputValue(&dailyrainfall[iv],gauss);
    4718                 dailysnowfall_input->GetInputValue(&dailysnowfall[iv],gauss);
    4719                 dailydlradiation_input->GetInputValue(&dailydlradiation[iv],gauss);
    4720                 dailydsradiation_input->GetInputValue(&dailydsradiation[iv],gauss);
    4721                 dailywindspeed_input->GetInputValue(&dailywindspeed[iv],gauss);
    4722                 dailypressure_input->GetInputValue(&dailypressure[iv],gauss);
    4723                 dailyairdensity_input->GetInputValue(&dailyairdensity[iv],gauss);
    4724                 dailyairhumidity_input->GetInputValue(&dailyairhumidity[iv],gauss);
    4725                 dailytemperature_input->GetInputValue(&dailytemperature[iv],gauss);
    4726                 tsurf_input->GetInputValue(&tsurf_in[iv],gauss);
    4727 
    4728                 /* Get Albedo information */
    4729                 albedo_input->GetInputValue(&albedo_in[iv],gauss);
    4730                 albedo_snow_input->GetInputValue(&albedo_snow_in[iv],gauss);
    4731                 mask_input->GetInputValue(&mask_in[iv],gauss);
    4732                 Tamp_input->GetInputValue(&Tamp_in[iv],gauss);
    4733 
    4734                 hsnow_input->GetInputValue(&hsnow_in[iv],gauss);
    4735                 hice_input->GetInputValue(&hice_in[iv],gauss);
    4736                 qmr_input->GetInputValue(&qmr_in[iv],gauss);
    4737 
    4738                 /* Surface temperature correction */
    4739                 st[iv]=(s[iv]-s0gcm[iv])/1000.;
    4740                 dailytemperature[iv]=dailytemperature[iv]-rlaps *st[iv];
    4741 
    4742                 /* Precipitation correction (Vizcaino et al. 2010) */
    4743                 if (s0gcm[iv] < desfacElev) {
    4744                         dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev));
    4745                         dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev));
    4746                 }else{
    4747                         dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv]));
    4748                         dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv]));
    4749                 }
    4750 
    4751                 /* downward longwave radiation correction (Marty et al. 2002) */
    4752                 st[iv]=(s[iv]-s0gcm[iv])/1000.;
    4753                 dailydlradiation[iv]=dailydlradiation[iv]+rdl*st[iv];
    4754         }
    4755         if(isverbose && this->Sid()==0){
    4756                 _printf0_("smb core: assign tsurf_in        :" << tsurf_in[0] << "\n");
    4757                 _printf0_("smb core: assign dailytemperature:" << dailytemperature[0] << "\n");
    4758                 _printf0_("smb core: assign hsnow           :" << hsnow_in[0] << "\n");
    4759                 _printf0_("smb core: assign hice            :" << hice_in[0] << "\n");
    4760                 _printf0_("smb core: assign mask            :" << mask_in[0] << "\n");
    4761                 _printf0_("smb core: assign Tamp            :" << Tamp_in[0] << "\n");
    4762                 _printf0_("smb core: assign albedo          :" << albedo_in[0] << "\n");
    4763                 _printf0_("smb core: assign albedo_snow     :" << albedo_snow_in[0] << "\n");
    4764                 _printf0_("smb core: assign albedo_scheme   :" << alb_scheme  << "\n");
    4765                 _printf0_("smb core: assign qmr             :" << qmr_in[0]  << "\n");
    4766         }
    4767 
    4768         if(isverbose && this->Sid()==0)_printf0_("smb core: call run_semic_transient module.\n");
    4769         /* call semic */
    4770         int nx=NUM_VERTICES, ntime=1, nloop=1;
    4771         bool semic_verbose=false; //VerboseSmb();
    4772         run_semic_transient_(&nx, &ntime, &nloop,
    4773                         dailysnowfall,  dailyrainfall, dailydsradiation, dailydlradiation,
    4774                         dailywindspeed, dailypressure, dailyairdensity,  dailyairhumidity, dailytemperature, tsurf_in, qmr_in,
    4775                         &dt,
    4776                         &hcrit, &rcrit,
    4777                         mask_in, hice_in, hsnow_in,
    4778                         albedo_in, albedo_snow_in,
    4779                         &alb_scheme, &alb_smax, &alb_smin, &albi, &albl,
    4780                         Tamp_in,
    4781                         &tmin, &tmax, &tmid, &mcrit, &wcrit, &tau_a, &tau_f, &afac, &semic_verbose,
    4782                         tsurf_out, smb_out, smbi_out, smbs_out, saccu_out, smelt_out, refr_out, albedo_out, albedo_snow_out, hsnow_out, hice_out, qmr_out);
    4783 
    4784         for (int iv = 0; iv<NUM_VERTICES; iv++){
    4785                 /*
    4786                  unit conversion: water -> ice
    4787                  w.e. : water equivalenet.
    4788                  */
    4789                 smb_out[iv]  = smb_out[iv]*yts;  // w.e. m/sec -> m/yr
    4790                 smbi_out[iv] = smbi_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
    4791                 smbs_out[iv] = smbs_out[iv]*yts; // w.e. m/sec -> m/yr
    4792                 saccu_out[iv] = saccu_out[iv]*yts; // w.e. m/sec -> m/yr
    4793                 smelt_out[iv] = smelt_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
    4794                 refr_out[iv]  = refr_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
    4795         }
    4796 
    4797         if(isverbose && this->Sid()==0){
    4798                 _printf0_("smb core: tsurf_out " << tsurf_out[0] << " " << tsurf_out[1] << " " << tsurf_out[2] << "\n");
    4799                 _printf0_("smb core: hice_out  " << hice_out[0] << " " <<  hice_out[1] << " " << hice_out[2] << "\n");
    4800                 _printf0_("smb core: hsnow_out " << hsnow_out[0] << "\n");
    4801                 _printf0_("smb core: smb_ice   " << smbi_out[0]*yts << "\n");
    4802                 _printf0_("smb core: smb_ice   " << albedo_out[0] <<" "<<albedo_out[1] << " " << albedo_out[2] << "\n");
    4803         }
    4804 
    4805         switch(this->ObjectEnum()){
    4806                 case TriaEnum:
    4807                         this->AddInput(TemperatureSEMICEnum,  &tsurf_out[0],P1DGEnum);
    4808                         // SMBout = SMB_ice + SMB_snow values.
    4809                         //this->AddInput(SmbMassBalanceTotalEnum,&smb_out[0],P1DGEnum);
    4810                         // water equivalent SMB ice to ice equivalent.
    4811                         this->AddInput(SmbMassBalanceEnum,     &smbi_out[0],P1DGEnum);
    4812                         this->AddInput(SmbMassBalanceIceEnum,  &smbi_out[0],P1DGEnum);
    4813                         this->AddInput(SmbMassBalanceSnowEnum, &smbs_out[0],P1DGEnum);
    4814                         this->AddInput(SmbMassBalanceSemicEnum,&smb_out[0],P1DGEnum);
    4815                         //this->AddInput(SmbMassBalanceSnowEnum,&smbs_out[0],P1DGEnum);
    4816                         // saccu - accumulation of snow.
    4817                         this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1DGEnum);
    4818                         // smelt
    4819                         this->AddInput(SmbMeltEnum,        &smelt_out[0],P1DGEnum);
    4820                         this->AddInput(SmbRefreezeEnum,    &refr_out[0],P1DGEnum);
    4821                         this->AddInput(SmbAlbedoEnum,      &albedo_out[0],P1DGEnum);
    4822                         this->AddInput(SmbAlbedoSnowEnum,  &albedo_snow_out[0],P1DGEnum);
    4823                         this->AddInput(SmbHSnowEnum,       &hsnow_out[0],P1DGEnum);
    4824                         this->AddInput(SmbHIceEnum,        &hice_out[0],P1DGEnum);
    4825                         this->AddInput(SmbSemicQmrEnum,    &qmr_out[0],P1DGEnum);
    4826                         break;
    4827                 case PentaEnum:
    4828                         // TODO
    4829                         break;
    4830                 case TetraEnum:
    4831                         // TODO
    4832                         break;
    4833                 default: _error_("Not implemented yet");
    4834         }
    4835 
    4836         /*clean-up {{{*/
    4837         delete gauss;
    4838         xDelete<IssmDouble>(dailysnowfall);
    4839         xDelete<IssmDouble>(dailyrainfall);
    4840         xDelete<IssmDouble>(dailydlradiation);
    4841         xDelete<IssmDouble>(dailydsradiation);
    4842         xDelete<IssmDouble>(dailywindspeed);
    4843         xDelete<IssmDouble>(dailypressure);
    4844         xDelete<IssmDouble>(dailyairdensity);
    4845         xDelete<IssmDouble>(dailyairhumidity);
    4846         xDelete<IssmDouble>(dailypressure);
    4847         xDelete<IssmDouble>(dailytemperature);
    4848 
    4849         /*for outputs*/
    4850         xDelete<IssmDouble>(tsurf_out);
    4851         xDelete<IssmDouble>(smb_out);
    4852         xDelete<IssmDouble>(smbi_out);
    4853         xDelete<IssmDouble>(smbs_out);
    4854         xDelete<IssmDouble>(saccu_out);
    4855         xDelete<IssmDouble>(smelt_out);
    4856         xDelete<IssmDouble>(refr_out);
    4857         xDelete<IssmDouble>(albedo_out);
    4858         xDelete<IssmDouble>(albedo_snow_out);
    4859         xDelete<IssmDouble>(hsnow_out);
    4860         xDelete<IssmDouble>(hice_out);
    4861         xDelete<IssmDouble>(qmr_out);
    4862 
    4863         /*for inputs*/
    4864         xDelete<IssmDouble>(hsnow_in);
    4865         xDelete<IssmDouble>(hice_in);
    4866         xDelete<IssmDouble>(mask_in);
    4867         xDelete<IssmDouble>(Tamp_in);
    4868         xDelete<IssmDouble>(albedo_in);
    4869         xDelete<IssmDouble>(albedo_snow_in);
    4870         xDelete<IssmDouble>(tsurf_in);
    4871         xDelete<IssmDouble>(qmr_in);
    4872 
    4873         /* for inputs:geometry */
    4874         xDelete<IssmDouble>(s);
    4875         xDelete<IssmDouble>(st);
    4876         xDelete<IssmDouble>(s0gcm);
    4877         /*}}}*/
    48784121}
    48794122/*}}}*/
    … …  
    54204663}
    54214664/*}}}*/
    5422 void       Element::SubglacialWaterPressure(int output_enum){/*{{{*/
    5423 
    5424         bool ispwHydroArma;
    5425    int M;
    5426    int numvertices = this->GetNumberOfVertices();
    5427    IssmDouble p_water[numvertices];
    5428    IssmDouble* perturbationvalues = xNew<IssmDouble>(numvertices);
    5429    Gauss* gauss=this->NewGauss();
    5430    Friction* friction = new Friction(this);
    5431    /*Calculate subglacial water pressure*/
    5432    for(int i=0;i<numvertices;i++){
    5433          gauss->GaussVertex(i);
    5434          p_water[i] = friction->SubglacialWaterPressure(gauss);
    5435    }
    5436    /*Add perturbation in water pressure if HydrologyIsWaterPressureArmaEnum is true*/
    5437    this->parameters->FindParam(&ispwHydroArma,HydrologyIsWaterPressureArmaEnum);
    5438    if(ispwHydroArma){
    5439       this->GetInputListOnVertices(perturbationvalues,WaterPressureArmaPerturbationEnum);
    5440                 for(int i=0;i<numvertices;i++) p_water[i] = p_water[i]+perturbationvalues[i];
    5441    }
    5442    /*Save*/
    5443    this->AddInput(output_enum,p_water,P1DGEnum);
    5444    /*Clean-up*/
    5445    delete gauss;
    5446    delete friction;
    5447    xDelete<IssmDouble>(perturbationvalues);
    5448 
    5449 }/*}}}*/
    54504665void       Element::StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    54514666
    … …  
    57554970}
    57564971/*}}}*/
    5757 IssmDouble Element::TotalSmbMelt(IssmDouble* mask, bool scaled){/*{{{*/
    5758 
    5759         /*Retrieve values of the mask defining the element: */
    5760         for(int i=0;i<this->GetNumberOfVertices();i++){
    5761                 if(mask[this->vertices[i]->Sid()]<=0.){
    5762                         return 0.;
    5763                 }
    5764         }
    5765 
    5766         /*Return: */
    5767         return this->TotalSmbMelt(scaled);
    5768 }
    5769 /*}}}*/
    5770 IssmDouble Element::TotalSmbRefreeze(IssmDouble* mask, bool scaled){/*{{{*/
    5771 
    5772         /*Retrieve values of the mask defining the element: */
    5773         for(int i=0;i<this->GetNumberOfVertices();i++){
    5774                 if(mask[this->vertices[i]->Sid()]<=0.){
    5775                         return 0.;
    5776                 }
    5777         }
    5778 
    5779         /*Return: */
    5780         return this->TotalSmbRefreeze(scaled);
    5781 }
    5782 /*}}}*/
    57834972void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
    57844973

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Element.h

    r27856 r27956  
    5757
    5858                int* element_type_list;
     59                IssmDouble* accumulator_values;
    5960                int  element_type;
    6061
    … …  
    6869                /*bool               AnyActive(void);*/
    6970                bool               AnyFSet(void);
    70       void               ArmaProcess_pre18Oct2022(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type);
    71       void               ArmaProcess(bool isstepforarma,int arorder,int maorder,int numparams,int numbreaks,IssmDouble tstep_arma,IssmDouble* polyparams,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,IssmDouble* datebreaks,bool isfieldstochastic,int enum_type);
     71      void               ArmaProcess(bool isstepforarma,int arorder,int maorder,IssmDouble telapsed_arma,IssmDouble tstep_arma,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* arlagcoefs,IssmDouble* malagcoefs,bool isfieldstochastic,int enum_type);
     72      void               Autoregression(bool isstepforar,int arorder,IssmDouble telapsed_ar,IssmDouble tstep_ar,IssmDouble* termconstant,IssmDouble* trend,IssmDouble* lagcoefs,bool isfieldstochastic,int enum_type);
    7273                void               BasinLinearFloatingiceMeltingRate(IssmDouble* deepwaterel,IssmDouble* upperwatermelt,IssmDouble* upperwaterel,IssmDouble* perturbation);
    7374                void               CalvingSetZeroRate(void);
    … …  
    140141                void               InputCreateP1FromConstant(Inputs* inputs,IoModel* iomodel,IssmDouble value,int vector_enum);
    141142                void               ControlInputCreate(IssmDouble* doublearray,IssmDouble* independents_min,IssmDouble* independents_max,Inputs*inputs,IoModel* iomodel,int M,int N,IssmDouble scale,int input_enum,int id);
    142                 void                                     DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int input_enum);
     143                void                                     DatasetInputAdd(int enum_type,IssmDouble* vector,Inputs* inputs,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum);
    143144                void               InputUpdateFromConstant(IssmDouble constant, int name);
    144145                void               InputUpdateFromConstant(int constant, int name);
    … …  
    155156                bool               IsOceanInElement();
    156157                bool               IsOceanOnlyInElement();
    157                 bool               IsAllMinThicknessInElement();
    158158                bool               IsLandInElement();
    159159                void               Ismip6FloatingiceMeltingRate();
    … …  
    165165                void               MigrateGroundingLine(IssmDouble* sheet_ungrounding);
    166166                void               MismipFloatingiceMeltingRate();
    167                 void               MonthlyFactorBasin(IssmDouble* monthlyfac,int enum_type);
    168                 void               MonthlyPiecewiseLinearEffectBasin(int nummonthbreaks,IssmDouble* monthlyintercepts,IssmDouble* monthlytrends,IssmDouble* monthlydatebreaks,int enum_type);
     167                void               MonthlyEffectBasin(IssmDouble* monthlyeff, int enum_type);
    169168                void               BeckmannGoosseFloatingiceMeltingRate();
    170169                void               MungsmtpParameterization(void);
    … …  
    177176                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    178177                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
    179                 void               SmbDebrisEvatt();
    180178                void               RignotMeltParameterization();
    181179                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);
    … …  
    188186                void               SetIntInput(Inputs* inputs,int enum_in,int value);
    189187                void               SmbSemic();
    190                 void               SmbSemicTransient();
    191188                int                Sid();
    192189                void               SmbGemb(IssmDouble timeinputs, int count);
    … …  
    199196                void               StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input);
    200197                void               StressMaxPrincipalCreateInput(void);
    201                 void               SubglacialWaterPressure(int output_enum);
    202198                IssmDouble         TotalFloatingBmb(IssmDouble* mask, bool scaled);
    203199                IssmDouble         TotalGroundedBmb(IssmDouble* mask, bool scaled);
    204200                IssmDouble         TotalSmb(IssmDouble* mask, bool scaled);
    205                 IssmDouble         TotalSmbMelt(IssmDouble* mask, bool scaled);
    206                 IssmDouble         TotalSmbRefreeze(IssmDouble* mask, bool scaled);
    207201                void               TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum);
    208202                void               TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
    … …  
    241235                virtual void             BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){_error_("not implemented yet");};
    242236                virtual void       CalvingRateParameterization(void){_error_("not implemented yet");};
    243                 virtual void       CalvingRateCalvingMIP(void){_error_("not implemented yet");};
    244237                virtual void       CalvingRateVonmises(void){_error_("not implemented yet");};
    245                 virtual void       CalvingRateVonmisesAD(void){_error_("not implemented yet");};
    246238                virtual void       CalvingRateTest(void){_error_("not implemented yet");};
    247239                virtual void       CalvingCrevasseDepth(void){_error_("not implemented yet");};
    … …  
    332324                virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
    333325                virtual IssmDouble Misfit(int modelenum,int observationenum,int weightsenum)=0;
     326                virtual void       MisfitAccumulate(int modelenum,IssmDouble dt)=0;
     327                virtual IssmDouble MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt)=0;
    334328                virtual IssmDouble MisfitArea(int weightsenum)=0;
    335329                virtual void       MovingFrontalVelocity(void){_error_("not implemented yet");};
    … …  
    390384                virtual IssmDouble TotalGroundedBmb(bool scaled)=0;
    391385                virtual IssmDouble TotalSmb(bool scaled)=0;
    392                 virtual IssmDouble TotalSmbMelt(bool scaled)=0;
    393                 virtual IssmDouble TotalSmbRefreeze(bool scaled)=0;
    394386                virtual void       Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
    395387                virtual void       UpdateConstraintsExtrudeFromBase(void)=0;
    … …  
    417409                virtual void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom)=0;
    418410                virtual void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom)=0;
    419                 virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount)=0;
     411                virtual void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids)=0;
    420412                virtual void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae)=0;
    421413                virtual void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae)=0;

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.cpp

    r27909 r27956  
    177177                        }
    178178                }
    179                 else if(interpolation_enum==P0Enum){
    180                         Penta* penta=this;
    181                         for(;;){
    182                                 penta->AddInput(input_enum,&values[0],interpolation_enum);
    183                                 if (penta->IsOnSurface()) break;
    184                                 penta=penta->GetUpperPenta(); _assert_(penta->Id()!=this->id);
    185                         }
    186                 }
    187                 else _error_("Interpolation "<<EnumToStringx(interpolation_enum)<<" not implemented yet");
     179                else _error_("not implemented yet");
    188180        }
    189181
    … …  
    717709        IssmDouble  xyz_list[NUMVERTICES][3];
    718710        IssmDouble  viscosity;
    719         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     711        IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
    720712        IssmDouble  tau_xx[NUMVERTICES];
    721713        IssmDouble      tau_yy[NUMVERTICES];
    … …  
    730722
    731723        /*Retrieve all inputs we will be needing: */
    732         Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
    733         Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
    734         Input* vz_input=this->GetInput(VzEnum); _assert_(vz_input);
     724        Input* vx_input=this->GetInput(VxEnum);             _assert_(vx_input);
     725        Input* vy_input=this->GetInput(VyEnum);             _assert_(vy_input);
     726        Input* vz_input=this->GetInput(VzEnum);             _assert_(vz_input);
    735727
    736728        /* Start looping on the number of vertices: */
    … …  
    14511443IssmDouble Penta::GetIcefrontArea(){/*{{{*/
    14521444
    1453         /*We need to be on base and cross the levelset*/
     1445        IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
     1446        IssmDouble      Haverage,frontarea;
     1447        IssmDouble  x1,y1,x2,y2,distance;
     1448        IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
     1449        int* indices=NULL;
     1450        IssmDouble* H=NULL;;
     1451        int nrfrontbed,numiceverts;
     1452
     1453        if(!this->IsOnBase()) return 0;
    14541454        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    1455         if(!this->IsOnBase()) return 0;
    1456 
    1457         /*Spawn Tria element from the base of the Penta: */
    1458         Tria* tria=(Tria*)SpawnTria(0,1,2);
    1459         IssmDouble frontarea = tria->GetIcefrontArea();
    1460         delete tria->material; delete tria;
    1461 
     1455
     1456        /*Retrieve all inputs and parameters*/
     1457        Element::GetInputListOnVertices(&bed[0],BedEnum);
     1458        Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
     1459        Element::GetInputListOnVertices(&bases[0],BaseEnum);
     1460        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     1461
     1462        nrfrontbed=0;
     1463        for(int i=0;i<NUMVERTICES2D;i++){
     1464                /*Find if bed<0*/
     1465                if(bed[i]<0.) nrfrontbed++;
     1466        }
     1467
     1468        if(nrfrontbed==3){
     1469                /*2. Find coordinates of where levelset crosses 0*/
     1470                int         numiceverts;
     1471                IssmDouble  s[2],x[2],y[2];
     1472                this->GetLevelsetIntersectionBase(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
     1473                _assert_(numiceverts);
     1474
     1475                /*3 Write coordinates*/
     1476                IssmDouble  xyz_list[NUMVERTICES][3];
     1477                ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     1478                int counter = 0;
     1479                if((numiceverts>0) && (numiceverts<NUMVERTICES2D)){
     1480                        for(int i=0;i<numiceverts;i++){
     1481                                for(int n=numiceverts;n<NUMVERTICES2D;n++){ // iterate over no-ice vertices
     1482                                        x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
     1483                                        y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
     1484                                        counter++;
     1485                                }
     1486                        }
     1487                }
     1488                else if(numiceverts==NUMVERTICES2D){ //NUMVERTICES ice vertices: calving front lies on element edge
     1489
     1490                        for(int i=0;i<NUMVERTICES2D;i++){
     1491                                if(lsf[indices[i]]==0.){
     1492                                        x[counter]=xyz_list[indices[i]][0];
     1493                                        y[counter]=xyz_list[indices[i]][1];
     1494                                        counter++;
     1495                                }
     1496                                if(counter==2) break;
     1497                        }
     1498                        if(counter==1){
     1499                                /*We actually have only 1 vertex on levelset, write a single point as a segment*/
     1500                                x[counter]=x[0];
     1501                                y[counter]=y[0];
     1502                                counter++;
     1503                        }
     1504                }
     1505                else{
     1506                        _error_("not sure what's going on here...");
     1507                }
     1508                x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
     1509                distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
     1510
     1511                int numthk=numiceverts+2;
     1512                H=xNew<IssmDouble>(numthk);
     1513                for(int iv=0;iv<NUMVERTICES2D;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
     1514
     1515                switch(numiceverts){
     1516                        case 1: // average over triangle
     1517                                H[0]=Haux[0];
     1518                                H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
     1519                                H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
     1520                                Haverage=(H[1]+H[2])/2;
     1521                                break;
     1522                        case 2: // average over quadrangle
     1523                                H[0]=Haux[0];
     1524                                H[1]=Haux[1];
     1525                                H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
     1526                                H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
     1527                                Haverage=(H[2]+H[3])/2;
     1528                                break;
     1529                        default:
     1530                                _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
     1531                                break;
     1532                }
     1533                frontarea=distance*Haverage;
     1534        }
     1535        else return 0;
     1536
     1537        xDelete<int>(indices);
     1538        xDelete<IssmDouble>(H);
     1539
     1540        _assert_(frontarea>0);
    14621541        return frontarea;
    1463 }/*}}}*/
     1542}
     1543/*}}}*/
    14641544void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    14651545
    … …  
    34113491        }
    34123492
     3493        /*Get out if this is not an element input*/
     3494        if(!IsInputEnum(control_enum)) return;
     3495
    34133496        /*Prepare index list*/
    34143497        ElementInput* input=this->inputs->GetControlInputData(control_enum,"value");   _assert_(input);
    … …  
    42704353        /*Return: */
    42714354        return Total_Smb;
    4272 }
    4273 /*}}}*/
    4274 IssmDouble Penta::TotalSmbMelt(bool scaled){/*{{{*/
    4275 
    4276         /*The smbmelt[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    4277         IssmDouble base,smbmelt,rho_ice,scalefactor;
    4278         IssmDouble Total_SmbMelt=0;
    4279         IssmDouble lsf[NUMVERTICES];
    4280         IssmDouble xyz_list[NUMVERTICES][3];
    4281 
    4282         /*Get material parameters :*/
    4283         rho_ice=FindParam(MaterialsRhoIceEnum);
    4284 
    4285         if(!IsIceInElement() || !IsOnSurface()) return 0.;
    4286 
    4287         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4288 
    4289         /*First calculate the area of the base (cross section triangle)
    4290          * http://en.wikipedia.org/wiki/Triangle
    4291          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    4292         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    4293 
    4294         /*Now get the average SMB over the element*/
    4295         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    4296    if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
    4297                 /*Partially ice-covered element*/
    4298                 bool mainlyice;
    4299       int point;
    4300       IssmDouble* smbmelt_vertices = xNew<IssmDouble>(NUMVERTICES);
    4301                 IssmDouble  weights[NUMVERTICES2D];
    4302                 IssmDouble  lsf2d[NUMVERTICES2D];
    4303       IssmDouble f1,f2,phi;
    4304       Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum);
    4305                 for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
    4306                 GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
    4307                 smbmelt = 0.0;
    4308                 for(int i=0;i<NUMVERTICES2D;i++) smbmelt += weights[i]*smbmelt_vertices[i];
    4309 
    4310                 if(scaled==true){
    4311          IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
    4312          Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    4313          /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
    4314          scalefactor = 0.0;
    4315          for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    4316          xDelete<IssmDouble>(scalefactor_vertices);
    4317                 }
    4318                 else scalefactor = 1.0;
    4319 
    4320                 /*Cleanup*/
    4321       xDelete<IssmDouble>(smbmelt_vertices);
    4322         }
    4323 
    4324         else{
    4325                 /*Fully ice-covered element*/
    4326                 Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input);
    4327                 smbmelt_input->GetInputAverage(&smbmelt);
    4328 
    4329                 if(scaled==true){
    4330                         Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    4331                         scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
    4332                 }
    4333                 else scalefactor=1.0;
    4334         }
    4335 
    4336         Total_SmbMelt=rho_ice*base*smbmelt*scalefactor;// smbmelt on element in kg s-1
    4337 
    4338         /*Return: */
    4339         return Total_SmbMelt;
    4340 }
    4341 /*}}}*/
    4342 IssmDouble Penta::TotalSmbRefreeze(bool scaled){/*{{{*/
    4343 
    4344         /*The smbrefreeze[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
    4345         IssmDouble base,smbrefreeze,rho_ice,scalefactor;
    4346         IssmDouble Total_SmbRefreeze=0;
    4347         IssmDouble lsf[NUMVERTICES];
    4348         IssmDouble xyz_list[NUMVERTICES][3];
    4349 
    4350         /*Get material parameters :*/
    4351         rho_ice=FindParam(MaterialsRhoIceEnum);
    4352 
    4353         if(!IsIceInElement() || !IsOnSurface()) return 0.;
    4354 
    4355         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4356 
    4357         /*First calculate the area of the base (cross section triangle)
    4358          * http://en.wikipedia.org/wiki/Triangle
    4359          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    4360         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
    4361 
    4362         /*Now get the average SMB over the element*/
    4363         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    4364    if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
    4365                 /*Partially ice-covered element*/
    4366                 bool mainlyice;
    4367       int point;
    4368       IssmDouble* smbrefreeze_vertices = xNew<IssmDouble>(NUMVERTICES);
    4369                 IssmDouble  weights[NUMVERTICES2D];
    4370                 IssmDouble  lsf2d[NUMVERTICES2D];
    4371       IssmDouble f1,f2,phi;
    4372       Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum);
    4373                 for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
    4374                 GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
    4375                 smbrefreeze = 0.0;
    4376                 for(int i=0;i<NUMVERTICES2D;i++) smbrefreeze += weights[i]*smbrefreeze_vertices[i];
    4377 
    4378                 if(scaled==true){
    4379          IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
    4380          Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    4381          /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
    4382          scalefactor = 0.0;
    4383          for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    4384          xDelete<IssmDouble>(scalefactor_vertices);
    4385                 }
    4386                 else scalefactor = 1.0;
    4387 
    4388                 /*Cleanup*/
    4389       xDelete<IssmDouble>(smbrefreeze_vertices);
    4390         }
    4391 
    4392         else{
    4393                 /*Fully ice-covered element*/
    4394                 Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input);
    4395                 smbrefreeze_input->GetInputAverage(&smbrefreeze);
    4396 
    4397                 if(scaled==true){
    4398                         Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    4399                         scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
    4400                 }
    4401                 else scalefactor=1.0;
    4402         }
    4403 
    4404         Total_SmbRefreeze=rho_ice*base*smbrefreeze*scalefactor;// smbrefreeze on element in kg s-1
    4405 
    4406         /*Return: */
    4407         return Total_SmbRefreeze;
    44084355}
    44094356/*}}}*/

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Penta.h

    r27850 r27956  
    6464                void           ComputeSigmaNN(){_error_("not implemented yet");};
    6565                void           ComputeStressTensor();
    66                 //void           ComputeMeanEla(IssmDouble* paltitude, int* pcounter);
    6766                void           Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs* inputsin);
    6867                void           ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp);
    … …  
    138137                IssmDouble     MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
    139138                IssmDouble     MinEdgeLength(IssmDouble* xyz_list);
    140                 IssmDouble     Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
     139        IssmDouble     Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
     140                void               MisfitAccumulate(int modelenum,IssmDouble dt){_error_("not implemented yet");};
     141                IssmDouble         MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt){_error_("not implemented yet");};
    141142                IssmDouble     MisfitArea(int weightsenum){_error_("not implemented yet");};
    142143                void           MovingFrontalVelocity(void);
    … …  
    201202                IssmDouble     TotalGroundedBmb(bool scaled);
    202203                IssmDouble     TotalSmb(bool scaled);
    203                 IssmDouble     TotalSmbMelt(bool scaled);
    204                 IssmDouble     TotalSmbRefreeze(bool scaled);
    205204                void           Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    206205                void           UpdateConstraintsExtrudeFromBase(void);
    … …  
    230229                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    231230                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    232                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");};
     231                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
    233232                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    234233                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/PentaRef.cpp

    r27720 r27956  
    307307                        basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
    308308                        return;
    309                         #ifndef _HAVE_AD_
    310309                case P2xP1Enum:
    311310                        /*Corner nodes*/
    … …  
    460459                        basis[14]=gauss->coord3*(-8./3.)*(zeta-1.0)*zeta*(zeta+0.5)*(zeta+1.);
    461460                        return;
    462                         #endif
    463461                default:
    464462                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    … …  
    573571                        dbasis[NUMNODESP1b*2+6] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta;
    574572                        return;
    575                         #ifndef _HAVE_AD_
    576573                case P2xP1Enum:
    577574                        /*Nodal function 1*/
    … …  
    10601057
    10611058                        return;
    1062                         #endif
    10631059                default:
    10641060                        _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Seg.h

    r27852 r27956  
    108108                IssmDouble  MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    109109                IssmDouble  Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
     110                void            MisfitAccumulate(int modelenum,IssmDouble dt){_error_("not implemented yet");};
     111                IssmDouble      MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt){_error_("not implemented yet");};
    110112                IssmDouble  MisfitArea(int weightsenum){_error_("not implemented yet");};
    111113                Gauss*      NewGauss(void);
    … …  
    158160                IssmDouble  TotalGroundedBmb(bool scaled){_error_("not implemented yet");};
    159161                IssmDouble  TotalSmb(bool scaled){_error_("not implemented yet");};
    160                 IssmDouble  TotalSmbMelt(bool scaled){_error_("not implemented yet");};
    161                 IssmDouble  TotalSmbRefreeze(bool scaled){_error_("not implemented yet");};
    162162                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
    163163                void        UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
    … …  
    182182                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    183183                void       SealevelchangeShift(GrdLoads* loads, IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    184                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids,int* vcount){_error_("not implemented yet");};
     184                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
    185185                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    186186                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tetra.h

    r27852 r27956  
    115115                IssmDouble  MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
    116116                IssmDouble  Misfit(int modelenum,int observationenum,int weightsenum){_error_("not implemented yet");};
     117                void            MisfitAccumulate(int modelenum,IssmDouble dt){_error_("not implemented yet");};
     118                IssmDouble      MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt){_error_("not implemented yet");};
    117119                IssmDouble  MisfitArea(int weightsenum){_error_("not implemented yet");};
    118120                Gauss*      NewGauss(void);
    … …  
    167169                IssmDouble  TotalGroundedBmb(bool scaled){_error_("not implemented yet");};
    168170                IssmDouble  TotalSmb(bool scaled){_error_("not implemented yet");};
    169                 IssmDouble  TotalSmbMelt(bool scaled){_error_("not implemented yet");};
    170                 IssmDouble  TotalSmbRefreeze(bool scaled){_error_("not implemented yet");};
    171171                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    172172                void        UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};
    … …  
    189189                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom){_error_("not implemented yet");};
    190190                void       SealevelchangeShift(GrdLoads* loads,  IssmDouble offset, SealevelGeometry* slgeom){_error_("not implemented yet");};
    191                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount){_error_("not implemented yet");};
     191                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){_error_("not implemented yet");};
    192192                void       SealevelchangeGeometryCentroidLoads(SealevelGeometry* slgeom, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae){_error_("not implemented yet");};
    193193                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae){_error_("not implemented yet");};

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tria.cpp

    r27909 r27956  
    2929#define NUMVERTICES   3
    3030#define NUMVERTICES1D 2
    31 //#define MICI          0 //1 = DeConto & Pollard, 2 = Anna Crawford DOMINOS
     31//#define MICI          1 //1 = DeConto & Pollard, 2 = DOMINOS
    3232
    3333/*Constructors/destructor/copy*/
    … …  
    4949                this->vertices = NULL;
    5050                this->material = NULL;
     51                this->accumulator_values = xNew<IssmDouble>(NUMVERTICES);
     52                for(int i=0;i<NUMVERTICES;i++) this->accumulator_values[i] = 0;
    5153                if(nummodels>0){
    5254                        this->element_type_list=xNew<int>(nummodels);
    … …  
    105107        else tria->element_type_list = NULL;
    106108        tria->element_type=this->element_type;
     109        tria->accumulator_values=this->accumulator_values;
    107110        tria->numanalyses=nanalyses;
    108111
    … …  
    345348                smax_fl_input->GetInputValue(&sigma_max_floating,&gauss);
    346349                smax_gr_input->GetInputValue(&sigma_max_grounded,&gauss);
    347                 sl_input->GetInputValue(&sealevel,&gauss);
    348 
    349                 /*Tensile stress threshold*/
    350                 if(groundedice<0)
    351                  sigma_max = sigma_max_floating;
    352                 else
    353                  sigma_max = sigma_max_grounded;
    354 
    355                 /*Assign values*/
    356                 if(bed>sealevel){
    357                         calvingrate[iv] = 0.;
    358                 }
    359                 else{
    360                         calvingrate[iv] = sqrt(vx*vx+vy*vy)*sigma_vm/sigma_max;
    361                 }
    362         }
    363 
    364         /*Add input*/
    365         this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    366    this->CalvingRateToVector();
    367 }
    368 /*}}}*/
    369 void       Tria::CalvingRateVonmisesAD(){/*{{{*/
    370 
    371         /*First, compute Von Mises Stress*/
    372         this->ComputeSigmaVM();
    373 
    374         /*Now compute calving rate*/
    375         IssmDouble  calvingrate[NUMVERTICES];
    376         IssmDouble  sigma_vm,vx,vy;
    377         IssmDouble  sigma_max,sigma_max_floating,sigma_max_grounded,n;
    378         IssmDouble  groundedice,bed,sealevel;
    379         int M;
    380         int basinid;
    381         IssmDouble* sigma_max_floating_basin=NULL;
    382         IssmDouble* sigma_max_grounded_basin=NULL;
    383 
    384         /*Retrieve all inputs and parameters we will need*/
    385         Input* vx_input       = this->GetInput(VxEnum); _assert_(vx_input);
    386         Input* vy_input       = this->GetInput(VyEnum); _assert_(vy_input);
    387         Input* gr_input       = this->GetInput(MaskOceanLevelsetEnum); _assert_(gr_input);
    388         Input* bs_input       = this->GetInput(BaseEnum);                    _assert_(bs_input);
    389         Input* sl_input       = this->GetInput(SealevelEnum); _assert_(sl_input);
    390         Input* sigma_vm_input = this->GetInput(SigmaVMEnum); _assert_(sigma_vm_input);
    391 
    392         this->Element::GetInputValue(&basinid,CalvingBasinIdEnum);
    393 
    394         parameters->FindParam(&sigma_max_floating_basin,&M,CalvingADStressThresholdFloatingiceEnum);
    395         parameters->FindParam(&sigma_max_grounded_basin,&M,CalvingADStressThresholdGroundediceEnum);
    396 
    397         sigma_max_floating = sigma_max_floating_basin[basinid];
    398         sigma_max_grounded = sigma_max_grounded_basin[basinid];
    399 
    400         /* Start looping on the number of vertices: */
    401         GaussTria gauss;
    402         for(int iv=0;iv<NUMVERTICES;iv++){
    403                 gauss.GaussVertex(iv);
    404 
    405                 /*Get velocity components and thickness*/
    406                 sigma_vm_input->GetInputValue(&sigma_vm,&gauss);
    407                 vx_input->GetInputValue(&vx,&gauss);
    408                 vy_input->GetInputValue(&vy,&gauss);
    409                 gr_input->GetInputValue(&groundedice,&gauss);
    410                 bs_input->GetInputValue(&bed,&gauss);
    411350                sl_input->GetInputValue(&sealevel,&gauss);
    412351
    … …  
    633572        IssmDouble ds, db, da, dt, dw, r, R;
    634573
    635         if(!IsIceInElement()){
    636                 for (int iv=0;iv<NUMVERTICES;iv++) calvingrate[iv]=0.;
    637                 this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    638                 this->CalvingRateToVector();
    639                 return;
    640         }
    641 
    642574        /*Retrieve all inputs and parameters we will need*/
    643575        IssmDouble rc        = FindParam(CalvingRcEnum);
    … …  
    689621
    690622         /*3. "Additional" crevasse opening*/
    691                         vel = sqrt(vx*vx + vy*vy);
    692                         da = H* max(0., log(vel*365*24*3600/1600.))/log(1.2);
     623         vel = sqrt(vx*vx + vy*vy)/(365.25*24*3600);
     624         da = H* max(0., log(vel/1600.))/log(1.2);
    693625
    694626         /*4. deal with shallow ice*/
    … …  
    698630         dw = 0.;
    699631         R = smb*365.25*24*3600; //convert from m/s to m/yr
    700                         if(R>1.5 && R<=3.){
     632         if(R>1.5 && R<=3.){
    701633            dw = 4*1.5*(R - 1.5);
    702634         }
    … …  
    707639         /*Total calving rate*/
    708640         r = (ds+db+da+dt+dw)/H;
    709                         //if(this->Id()==1){
    710                         //      printf("rc = %g\n",rc);
    711                         //      printf("ds = %g\n",ds);
    712                         //}
    713641                        calvingrate[iv]= mig_max * max(0., min(1., (r - rc)/(1 - rc))); //P&DC: mig_max = 3000 m/yr
    714642                        _assert_(!xIsNan<IssmDouble>(calvingrate[iv]));
    … …  
    11201048}
    11211049/*}}}*/
    1122 void       Tria::CalvingRateCalvingMIP(){/*{{{*/
    1123 
    1124         IssmDouble  calvingrate[NUMVERTICES];
    1125         IssmDouble  calvingratex[NUMVERTICES];
    1126         IssmDouble  calvingratey[NUMVERTICES];
    1127         int                     experiment = 1;  /* exp:1 by default */
    1128         int         dim, domaintype;
    1129         IssmDouble      vx, vy, vel, c, wrate;
    1130         IssmDouble  time, groundedice, yts;
    1131 
    1132         /*Get problem dimension and whether there is moving front or not*/
    1133         this->FindParam(&domaintype,DomainTypeEnum);
    1134         this->FindParam(&time,TimeEnum);
    1135         this->FindParam(&yts,ConstantsYtsEnum);
    1136 
    1137         switch(domaintype){
    1138                 case Domain2DverticalEnum:   dim = 1; break;
    1139                 case Domain2DhorizontalEnum: dim = 2; break;
    1140                 case Domain3DEnum:           dim = 2; break;
    1141                 default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    1142         }
    1143         if(dim==1) _error_("not implemented in 1D...");
    1144 
    1145         /*Retrieve all inputs and parameters we will need*/
    1146         Input *vx_input      = this->GetInput(VxEnum);                                _assert_(vx_input);
    1147         Input *vy_input      = this->GetInput(VyEnum);                                _assert_(vy_input);
    1148         Input *wrate_input   = this->GetInput(CalvingAblationrateEnum);               _assert_(wrate_input);
    1149         Input* gr_input      = this->GetInput(MaskOceanLevelsetEnum);                                           _assert_(gr_input);
    1150 
    1151         /* Use which experiment: use existing Enum */
    1152         this->FindParam(&experiment, CalvingUseParamEnum);
    1153 
    1154         /* Start looping on the number of vertices: */
    1155         GaussTria gauss;
    1156         for(int iv=0;iv<NUMVERTICES;iv++){
    1157                 gauss.GaussVertex(iv);
    1158 
    1159                 /*Get velocity components */
    1160                 vx_input->GetInputValue(&vx,&gauss);
    1161                 vy_input->GetInputValue(&vy,&gauss);
    1162                 vel=sqrt(vx*vx+vy*vy)+1.e-14;
    1163 
    1164                 /* no calving for grounded ice in EXP4 */
    1165                 gr_input->GetInputValue(&groundedice,&gauss);
    1166 
    1167                 switch (experiment) {
    1168                         case 1:
    1169                         case 3:
    1170                                 /* Exp 1 and 3: set c=v-wrate, wrate=0, so that w=0 */
    1171                                 wrate = 0.0;
    1172                                 break;
    1173                         case 2:
    1174                                 /* Exp 2: set c=v-wrate(given)*/
    1175                                 wrate_input->GetInputValue(&wrate,&gauss);
    1176                                 break;
    1177                         case 4:
    1178                                 /* Exp 4: set c=v-wrate(given), for the first 500 years, then c=0 for the second 500 years*/
    1179                                 if((groundedice<0) && (time<=500.0*yts)) {
    1180                                 //      wrate_input->GetInputValue(&wrate,&gauss);
    1181                                         wrate = -750*sin(2.0*M_PI*time/yts/1000)/yts;  // m/a -> m/s
    1182                                 }
    1183                                 else {
    1184                                         /* no calving on the grounded ice*/
    1185                                         wrate = vel;
    1186                                 }
    1187                                 break;
    1188                         default:
    1189                                 _error_("The experiment is not supported yet!");
    1190                 }
    1191 
    1192                 calvingrate[iv] = vel - wrate;
    1193                 calvingratex[iv] = vx - wrate*vx/vel;
    1194                 calvingratey[iv] = vy - wrate*vy/vel;
    1195         }
    1196         /*Add input*/
    1197         this->AddInput(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
    1198         this->AddInput(CalvingratexEnum,&calvingratex[0],P1DGEnum);
    1199         this->AddInput(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
    1200 }
    1201 /*}}}*/
    12021050IssmDouble Tria::CharacteristicLength(void){/*{{{*/
    12031051
    … …  
    23532201        //compute sea level load weights
    23542202        this->GetFractionGeometry(loadweights,&phi,&point1,&fraction1,&fraction2,&istrapeze1,levelset);
    2355 
    2356         //failsafe for phi so small that GetFractionGeometry returns 0 
    2357         if (phi==0) phi=1e-16;
    2358 
    23592203        for (int i=0;i<NUMVERTICES;i++) loadweights[i]/=phi;
    23602204        this->GetBarycenterFromLevelset(platbar,plongbar, phi, fraction1, fraction2, late, longe, point1,istrapeze1,planetradius);
    … …  
    27722616IssmDouble Tria::GetIcefrontArea(){/*{{{*/
    27732617
    2774         IssmDouble  bed[NUMVERTICES];
     2618        IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
    27752619        IssmDouble      Haverage,frontarea;
    27762620        IssmDouble  x1,y1,x2,y2,distance;
    27772621        IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
    27782622        int* indices=NULL;
    2779 
    2780         /*Return if no ice front present*/
     2623        IssmDouble* H=NULL;;
     2624        int nrfrontbed,numiceverts;
     2625
    27812626        if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;
    2782         //if(!this->IsIcefront()) return 0.;
    27832627
    27842628        /*Retrieve all inputs and parameters*/
    … …  
    27882632        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    27892633
    2790         /*Only continue if all 3 vertices are below sea level*/
    2791         for(int i=0;i<NUMVERTICES;i++) if(bed[i]>=0.) return 0.;
    2792 
    2793         /*2. Find coordinates of where levelset crosses 0*/
    2794         int         numiceverts;
    2795         IssmDouble  s[2],x[2],y[2];
    2796         this->GetLevelsetIntersection(&indices, &numiceverts, &s[0],MaskIceLevelsetEnum,0.);
    2797         _assert_(numiceverts);
    2798         if(numiceverts>2){
    2799                 Input* ls_input = this->GetInput(MaskIceLevelsetEnum);
    2800                 ls_input->Echo();
    2801         }
    2802 
    2803         /*3 Write coordinates*/
    2804         IssmDouble  xyz_list[NUMVERTICES][3];
    2805         ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
    2806         int counter = 0;
    2807         if((numiceverts>0) && (numiceverts<NUMVERTICES)){
    2808                 for(int i=0;i<numiceverts;i++){
    2809                         for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
    2810                                 x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
    2811                                 y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
     2634        nrfrontbed=0;
     2635        for(int i=0;i<NUMVERTICES;i++){
     2636                /*Find if bed<0*/
     2637                if(bed[i]<0.) nrfrontbed++;
     2638        }
     2639
     2640        if(nrfrontbed==3){
     2641                /*2. Find coordinates of where levelset crosses 0*/
     2642                int         numiceverts;
     2643                IssmDouble  s[2],x[2],y[2];
     2644                this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
     2645                _assert_(numiceverts);
     2646
     2647                /*3 Write coordinates*/
     2648                IssmDouble  xyz_list[NUMVERTICES][3];
     2649                ::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
     2650                int counter = 0;
     2651                if((numiceverts>0) && (numiceverts<NUMVERTICES)){
     2652                        for(int i=0;i<numiceverts;i++){
     2653                                for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
     2654                                        x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
     2655                                        y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
     2656                                        counter++;
     2657                                }
     2658                        }
     2659                }
     2660                else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
     2661
     2662                        for(int i=0;i<NUMVERTICES;i++){
     2663                                if(lsf[indices[i]]==0.){
     2664                                        x[counter]=xyz_list[indices[i]][0];
     2665                                        y[counter]=xyz_list[indices[i]][1];
     2666                                        counter++;
     2667                                }
     2668                                if(counter==2) break;
     2669                        }
     2670                        if(counter==1){
     2671                                /*We actually have only 1 vertex on levelset, write a single point as a segment*/
     2672                                x[counter]=x[0];
     2673                                y[counter]=y[0];
    28122674                                counter++;
    28132675                        }
    28142676                }
    2815         }
    2816         else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge
    2817 
    2818                 for(int i=0;i<NUMVERTICES;i++){
    2819                         if(lsf[indices[i]]==0.){
    2820                                 x[counter]=xyz_list[indices[i]][0];
    2821                                 y[counter]=xyz_list[indices[i]][1];
    2822                                 counter++;
    2823                         }
    2824                         if(counter==2) break;
    2825                 }
    2826                 if(counter==1){
    2827                         /*We actually have only 1 vertex on levelset, write a single point as a segment*/
    2828                         x[counter]=x[0];
    2829                         y[counter]=y[0];
    2830                         counter++;
    2831                 }
    2832         }
    2833         else{
    2834                 _error_("not sure what's going on here...");
    2835         }
    2836         x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
    2837         distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
    2838         if(distance<1e-3) return 0.;
    2839 
    2840         IssmDouble H[4];
    2841         for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
     2677                else{
     2678                        _error_("not sure what's going on here...");
     2679                }
     2680                x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
     2681                distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));
     2682
     2683                int numthk=numiceverts+2;
     2684                H=xNew<IssmDouble>(numthk);
     2685                for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice
     2686
     2687                switch(numiceverts){
     2688                        case 1: // average over triangle
     2689                                H[0]=Haux[0];
     2690                                H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
     2691                                H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
     2692                                Haverage=(H[1]+H[2])/2;
     2693                                break;
     2694                        case 2: // average over quadrangle
     2695                                H[0]=Haux[0];
     2696                                H[1]=Haux[1];
     2697                                H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
     2698                                H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
     2699                                Haverage=(H[2]+H[3])/2;
     2700                                break;
     2701                        default:
     2702                                _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
     2703                                break;
     2704                }
     2705                frontarea=distance*Haverage;
     2706        }
     2707        else return 0;
     2708
    28422709        xDelete<int>(indices);
    2843 
    2844         switch(numiceverts){
    2845                 case 1: // average over triangle
    2846                         H[0]=Haux[0];
    2847                         H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
    2848                         H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
    2849                         Haverage=(H[1]+H[2])/2;
    2850                         break;
    2851                 case 2: // average over quadrangle
    2852                         H[0]=Haux[0];
    2853                         H[1]=Haux[1];
    2854                         H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
    2855                         H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
    2856                         Haverage=(H[2]+H[3])/2;
    2857                         break;
    2858                 case 3:
    2859                         if(counter==1) distance = 0; //front has 0 width on this element because levelset is 0 at a single vertex
    2860                         else if(counter==2){ //two vertices with levelset=0: averaging ice front depth over both
    2861                                 Haverage = 0;
    2862                                 for(int i=0;i<NUMVERTICES;i++){
    2863                                         if(lsf[indices[i]]==0.) Haverage -= Haux[indices[i]]/2;
    2864                                         if(Haverage<Haux[indices[i]]/2-1e-3) break; //done with the two vertices
    2865                                 }
    2866                         }
    2867                         break;
    2868                 default:
    2869                         _error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
    2870                         break;
    2871         }
    2872         frontarea=distance*Haverage;
     2710        xDelete<IssmDouble>(H);
    28732711
    28742712        _assert_(frontarea>0);
    … …  
    44564294        IssmDouble Jdet;
    44574295        IssmDouble Jelem = 0;
     4296        IssmDouble Jelem1 = 0;
     4297        IssmDouble Jelem2 = 0;
     4298        IssmDouble Jelem3 = 0;
     4299        IssmDouble scaling = 0;
    44584300        IssmDouble xyz_list[NUMVERTICES][3];
    44594301        GaussTria *gauss = NULL;
    … …  
    44614303        /*If on water, return 0: */
    44624304        if(!IsIceInElement())return 0;
     4305
    44634306
    44644307        /*Retrieve all inputs we will be needing: */
    … …  
    44814324
    44824325                /*compute misfit between model and observation */
    4483                 Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
    4484         }
     4326                scaling=Jdet*weight*weight*gauss->weight;
     4327                Jelem+=0.5*(model-observation)*(model-observation);
     4328                //Jelem+=1e-19*log(fabs((model+1e-16)/(observation+1e-16)))*log(fabs((model+1e-16)/(observation+1e-16)));
     4329        }
     4330        Jelem*=scaling;
     4331        //Jelem1*=scaling;
     4332        //Jelem2=Jelem+Jelem1;
     4333    //_printf0_("   Tria.cpp:4327 Jtotal: "<<Jelem2<<"Jabs:"<<Jelem<<" Jrel:"<<Jelem1<<" model-obs:"<<(model-observation)<<" model:"<<model<<" obs:"<<observation<<"\n");
     4334
     4335        /* clean up and Return: */
     4336        delete gauss;
     4337        return Jelem;
     4338}
     4339/*}}}*/
     4340void Tria::MisfitAccumulate(int modelenum, IssmDouble dt){/*{{{*/
     4341        /*Intermediaries*/
     4342        IssmDouble model;
     4343        IssmDouble xyz_list[NUMVERTICES][3];
     4344        GaussTria *gauss = NULL;
     4345        int i = 0;
     4346
     4347        /*If on water, return 0: */
     4348        if(!IsIceInElement())return;
     4349
     4350        /*Retrieve all inputs we will be needing: */
     4351        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     4352        Input* model_input=this->GetInput(modelenum);   _assert_(model_input);
     4353
     4354        /* Start  looping on the number of gaussian points: */
     4355        gauss=new GaussTria(2);
     4356        while(gauss->next()){
     4357                /*Get parameters at gauss point*/
     4358                model_input->GetInputValue(&model,gauss);
     4359                if (model != model) {
     4360                        _printf0_("   Tria.cpp:4358 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " dt: " << dt << "\n");
     4361                } else {
     4362                        //_printf0_("   Tria.cpp:4361 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " dt: " << dt << "\n");
     4363                        this->accumulator_values[i++] += model * dt;
     4364                }
     4365        }
     4366        delete gauss;
     4367}
     4368/*}}}*/
     4369IssmDouble Tria::MisfitAnnual(int modelenum,int observationenum,int weightsenum, IssmDouble annual_dt){/*{{{*/
     4370        /*Intermediaries*/
     4371        IssmDouble model,observation,weight;
     4372        IssmDouble Jdet;
     4373        IssmDouble Jelem = 0;
     4374        IssmDouble Jelem1 = 0;
     4375        IssmDouble Jelem2 = 0;
     4376        IssmDouble Jelem3 = 0;
     4377        IssmDouble scaling = 0;
     4378        IssmDouble xyz_list[NUMVERTICES][3];
     4379        GaussTria *gauss = NULL;
     4380        int i = 0;
     4381
     4382        /*If on water, return 0: */
     4383        if(!IsIceInElement())return 0;
     4384
     4385        /*Retrieve all inputs we will be needing: */
     4386        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     4387        Input* model_input=this->GetInput(modelenum);   _assert_(model_input);
     4388        Input* observation_input=this->GetInput(observationenum);_assert_(observation_input);
     4389        Input* weights_input     =this->GetInput(weightsenum);     _assert_(weights_input);
     4390
     4391        /* Start  looping on the number of gaussian points: */
     4392        gauss=new GaussTria(2);
     4393        while(gauss->next()){
     4394
     4395                /* Get Jacobian determinant: */
     4396                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     4397
     4398                /*Get parameters at gauss point*/
     4399                observation_input->GetInputValue(&observation,gauss);
     4400                weights_input->GetInputValue(&weight,gauss);
     4401
     4402                /*compute misfit between model and observation */
     4403                scaling=Jdet*weight*gauss->weight;
     4404                model = this->accumulator_values[i] / annual_dt;
     4405                if (model != model ) {
     4406                        _printf0_("   Tria.cpp:4405 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " annual_dt: " << annual_dt << "\n");
     4407                        model = 0;
     4408                } else {
     4409                        //_printf0_("   Tria.cpp:4407 accumalator_value: " << this->accumulator_values[i] << " model: " << model << " annual_dt: " << annual_dt << "\n");
     4410                }
     4411                this->accumulator_values[i++] = 0;
     4412                Jelem+=0.5*(model-observation)*(model-observation);
     4413                //Jelem+=1e-19*log(fabs((model+1e-16)/(observation+1e-16)))*log(fabs((model+1e-16)/(observation+1e-16)));
     4414        }
     4415        Jelem*=scaling;
     4416        //Jelem1*=scaling;
     4417        //Jelem2=Jelem+Jelem1;
     4418    //_printf0_("   Tria.cpp:4327 Jtotal: "<<Jelem2<<"Jabs:"<<Jelem<<" Jrel:"<<Jelem1<<" model-obs:"<<(model-observation)<<" model:"<<model<<" obs:"<<observation<<"\n");
     4419    //_printf0_("   Tria.cpp:4327 Jtotal: " << Jelem << " model-obs:"<<(model-observation)<<" model:"<<model<<" obs:"<<observation<<"\n");
    44854420
    44864421        /* clean up and Return: */
    … …  
    45644499                case DefaultCalvingEnum:
    45654500                case CalvingVonmisesEnum:
    4566                 case CalvingVonmisesADEnum:
    45674501                case CalvingLevermannEnum:
    45684502                case CalvingPollardEnum:
    45694503                case CalvingTestEnum:
    45704504                case CalvingParameterizationEnum:
    4571                 case CalvingCalvingMIPEnum:
    45724505                        calvingratex_input=this->GetInput(CalvingratexEnum); _assert_(calvingratex_input);
    45734506                        calvingratey_input=this->GetInput(CalvingrateyEnum); _assert_(calvingratey_input);
    … …  
    46104543                        case DefaultCalvingEnum:
    46114544                        case CalvingVonmisesEnum:
    4612                         case CalvingVonmisesADEnum:
    46134545                        case CalvingTestEnum:
    46144546                        case CalvingParameterizationEnum:
    46154547                        case CalvingLevermannEnum:
    46164548                        case CalvingPollardEnum:
    4617                         case CalvingCalvingMIPEnum:
    46184549                                calvingratex_input->GetInputValue(&c[0],&gauss);
    46194550                                calvingratey_input->GetInputValue(&c[1],&gauss);
    … …  
    47294660
    47304661                /*Do we assume that the calving front does not move if MICI is not engaged?*/
    4731                 bool regrowth = false;
    4732                 bool apply_as_retreat = true;
    4733                 if(!regrowth){
    4734                         movingfrontvx[iv] = 0.;
    4735                         movingfrontvy[iv] = 0.;
    4736                 }
    4737 
     4662                movingfrontvx[iv] = 0.;
     4663                movingfrontvy[iv] = 0.;
    47384664                //movingfrontvx[iv] = -2000./(365*24*3600.)*dlsf[0]/norm_dlsf;
    47394665                //movingfrontvy[iv] = -2000./(365*24*3600.)*dlsf[1]/norm_dlsf;
    … …  
    47474673                }
    47484674                else if (MICI==2 && Hc>135. && bed<0. && fabs(ls)<100.e3){ // Crawford et all
    4749 
    4750                         /*if 1: RETREAT rate
    4751                          *if 0: calving rate*/
    4752                         if(0) v[0]=0.; v[1]=0.;
    47534675
    47544676                        /*5C Bn (worst case scenario)*/
    … …  
    47604682                        movingfrontvx[iv] = v[0] -C*dlsf[0]/norm_dlsf;
    47614683                        movingfrontvy[iv] = v[1] -C*dlsf[1]/norm_dlsf;
    4762 
    4763                         /*disable regrowth if calving rate is too low*/
    4764                         if(!regrowth && C<vel){
    4765                                 movingfrontvx[iv] = 0.;
    4766                                 movingfrontvy[iv] = 0.;
    4767                         }
    47684684                }
    47694685        }
    … …  
    51355051                }
    51365052        }
     5053
     5054        /*Get out if this is not an element input*/
     5055        if(!IsInputEnum(control_enum)) return;
    51375056
    51385057        /*Get list of ids for this element and this control*/
    … …  
    58875806        /*Return: */
    58885807        return Total_Smb;
    5889 }
    5890 /*}}}*/
    5891 IssmDouble Tria::TotalSmbMelt(bool scaled){/*{{{*/
    5892 
    5893         /*The smbmelt[kg yr-1] of one element is area[m2] * smbmelt [kg m^-2 yr^-1]*/
    5894         IssmDouble base,smbmelt,rho_ice,scalefactor;
    5895         IssmDouble Total_Melt=0;
    5896         IssmDouble lsf[NUMVERTICES];
    5897         IssmDouble xyz_list[NUMVERTICES][3];
    5898 
    5899         /*Get material parameters :*/
    5900         rho_ice=FindParam(MaterialsRhoIceEnum);
    5901 
    5902    if(!IsIceInElement())return 0;
    5903 
    5904         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5905 
    5906         /*First calculate the area of the base (cross section triangle)
    5907          * http://en.wikipedia.org/wiki/Triangle
    5908          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    5909         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
    5910        
    5911         /*Now get the average SMB over the element*/
    5912         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    5913         if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
    5914                 /*Partially ice-covered element*/
    5915                 bool mainlyice;
    5916       int point;
    5917       IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
    5918       IssmDouble* smbmelt_vertices  = xNew<IssmDouble>(NUMVERTICES);
    5919       IssmDouble f1,f2,phi;
    5920                
    5921                 Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum);
    5922                 GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
    5923                 smbmelt = 0.0;
    5924                 for(int i=0;i<NUMVERTICES;i++) smbmelt += weights[i]*smbmelt_vertices[i];
    5925        
    5926                 if(scaled==true){
    5927          IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
    5928          Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    5929          scalefactor = 0.0;
    5930          for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    5931          xDelete<IssmDouble>(scalefactor_vertices);
    5932       }
    5933                 else scalefactor = 1.0;
    5934 
    5935                 /*Cleanup*/
    5936       xDelete<IssmDouble>(weights);
    5937       xDelete<IssmDouble>(smbmelt_vertices);
    5938         }
    5939         else{
    5940                 /*Fully ice-covered element*/
    5941                 Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input);
    5942                 smbmelt_input->GetInputAverage(&smbmelt);   // average smbmelt on element in m ice s-1
    5943 
    5944                 if(scaled==true){
    5945                         Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    5946                         scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
    5947                 }
    5948                 else scalefactor=1.0;
    5949         }
    5950        
    5951    Total_Melt=rho_ice*base*smbmelt*scalefactor; // smbmelt on element in kg s-1
    5952 
    5953         /*Return: */
    5954         return Total_Melt;
    5955 }
    5956 /*}}}*/
    5957 IssmDouble Tria::TotalSmbRefreeze(bool scaled){/*{{{*/
    5958 
    5959         /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
    5960         IssmDouble base,smbrefreeze,rho_ice,scalefactor;
    5961         IssmDouble Total_Refreeze=0;
    5962         IssmDouble lsf[NUMVERTICES];
    5963         IssmDouble xyz_list[NUMVERTICES][3];
    5964 
    5965         /*Get material parameters :*/
    5966         rho_ice=FindParam(MaterialsRhoIceEnum);
    5967 
    5968    if(!IsIceInElement())return 0;
    5969 
    5970         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5971 
    5972         /*First calculate the area of the base (cross section triangle)
    5973          * http://en.wikipedia.org/wiki/Triangle
    5974          * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
    5975         base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
    5976        
    5977         /*Now get the average SMB over the element*/
    5978         Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
    5979         if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
    5980                 /*Partially ice-covered element*/
    5981                 bool mainlyice;
    5982       int point;
    5983       IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
    5984       IssmDouble* smbrefreeze_vertices  = xNew<IssmDouble>(NUMVERTICES);
    5985       IssmDouble f1,f2,phi;
    5986                
    5987                 Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum);
    5988                 GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
    5989                 smbrefreeze = 0.0;
    5990                 for(int i=0;i<NUMVERTICES;i++) smbrefreeze += weights[i]*smbrefreeze_vertices[i];
    5991        
    5992                 if(scaled==true){
    5993          IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
    5994          Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
    5995          scalefactor = 0.0;
    5996          for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
    5997          xDelete<IssmDouble>(scalefactor_vertices);
    5998       }
    5999                 else scalefactor = 1.0;
    6000 
    6001                 /*Cleanup*/
    6002       xDelete<IssmDouble>(weights);
    6003       xDelete<IssmDouble>(smbrefreeze_vertices);
    6004         }
    6005         else{
    6006                 /*Fully ice-covered element*/
    6007                 Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input);
    6008                 smbrefreeze_input->GetInputAverage(&smbrefreeze);   // average smbrefreeze on element in m ice s-1
    6009 
    6010                 if(scaled==true){
    6011                         Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
    6012                         scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
    6013                 }
    6014                 else scalefactor=1.0;
    6015         }
    6016        
    6017    Total_Refreeze=rho_ice*base*smbrefreeze*scalefactor; // smbrefreeze on element in kg s-1
    6018 
    6019         /*Return: */
    6020         return Total_Refreeze;
    60215808}
    60225809/*}}}*/
    … …  
    67186505}
    67196506/*}}}*/
    6720 void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* n_activevertices){ /*{{{*/
     6507void       Tria::SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids){ /*{{{*/
    67216508
    67226509        /*Declarations:{{{*/
    … …  
    68686655        }
    68696656
    6870         AlphaIndex=xNewZeroInit<int>(n_activevertices[this->lid]*nel);
    6871         if (horiz) AzimuthIndex=xNewZeroInit<int>(n_activevertices[this->lid]*nel);
     6657        AlphaIndex=xNewZeroInit<int>(3*nel);
     6658        if (horiz) AzimuthIndex=xNewZeroInit<int>(3*nel);
    68726659        int intmax=pow(2,16)-1;
    68736660
    6874         int* activevertices=xNew<int>(n_activevertices[this->lid]);
    6875        
    6876         int av=0;
    68776661
    68786662        for (int i=0;i<3;i++){
    68796663                if(lids[this->vertices[i]->lid]==this->lid){
    6880                         activevertices[av]=i;
    68816664                        for(int e=0;e<nel;e++){
    68826665                                IssmDouble alpha;
    … …  
    69026685                                        dy=sin(longe-longi)*cos(late);
    69036686                                        //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
    6904                                         AzimuthIndex[av*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6687                                        AzimuthIndex[i*nel+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    69056688                                }
    6906                                 AlphaIndex[av*nel+e]=index;
    6907                         }
    6908                         av+=1;                 
     6689                                AlphaIndex[i*nel+e]=index;
     6690                        }                       
    69096691                } //for (int i=0;i<3;i++)
    69106692        } //for(int e=0;e<nel;e++)
    69116693
    69126694        /*Add in inputs:*/
    6913         this->inputs->SetIntArrayInput(SealevelchangeConvolutionVerticesEnum,this->lid,activevertices,n_activevertices[this->lid]);
    6914         this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*n_activevertices[this->lid]);
    6915         if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*n_activevertices[this->lid]);
     6695        this->inputs->SetIntArrayInput(SealevelchangeAlphaIndexEnum,this->lid,AlphaIndex,nel*3);
     6696        if(horiz) this->inputs->SetIntArrayInput(SealevelchangeAzimuthIndexEnum,this->lid,AzimuthIndex,nel*3);
    69166697
    69176698        /*}}}*/
    … …  
    70326813        /*Free allocations:{{{*/
    70336814        #ifdef _HAVE_RESTRICT_
    7034         delete activevertices;
    70356815        delete AlphaIndex;
    70366816        if(horiz) AzimuthIndex;
    … …  
    70466826
    70476827        #else
    7048         xDelete<int>(activevertices);
    70496828        xDelete<int>(AlphaIndex);
    70506829        if(horiz){
    … …  
    70786857        IssmDouble x,y,z,dx,dy,dz,N_azim,E_azim;
    70796858        IssmDouble xyz_list[NUMVERTICES][3];
    7080         int* activevertices = NULL;
    7081         int n_activevertices, av;
    70826859
    70836860        #ifdef _HAVE_RESTRICT_
    … …  
    71546931        if(horiz) AzimIndex=xNew<int*>(SLGEOM_NUMLOADS);
    71556932
    7156         this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
    7157         // 0<=n_activevertices<=3 is the number of vertices this element is in charge of computing fields in during the sea level convolutions
    7158         // activevertices contains the vertex indices (1,2 and/or 3) in case debugging is required, they are supposed to appear in the same order as slgeom->lids
    71596933
    71606934        //Allocate:
    71616935        for(int l=0;l<SLGEOM_NUMLOADS;l++){
    71626936                int nbar=slgeom->nbar[l];
    7163                 AlphaIndex[l]=xNewZeroInit<int>(n_activevertices*nbar);
    7164                 if(horiz) AzimIndex[l]=xNewZeroInit<int>(n_activevertices*nbar);
    7165 
    7166                 //av=0;
    7167                 //for (int i=0;i<3;i++){
    7168                 for (int av=0;av<n_activevertices;av++){
    7169                         //if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    7170                         int i=activevertices[av];
    7171                         for(int e=0;e<nbar;e++){
    7172                                 IssmDouble alpha;
    7173                                 IssmDouble delPhi,delLambda;
    7174                                 /*recover info for this element and vertex:*/
    7175                                 IssmDouble late= slgeom->latbarycentre[l][e];
    7176                                 IssmDouble longe= slgeom->longbarycentre[l][e];
    7177                                 late=late/180*M_PI;
    7178                                 longe=longe/180*M_PI;
    7179                                 lati=latitude[i];
    7180                                 longi=longitude[i];
     6937                AlphaIndex[l]=xNewZeroInit<int>(3*nbar);
     6938                if(horiz) AzimIndex[l]=xNewZeroInit<int>(3*nbar);
     6939
     6940
     6941                for (int i=0;i<3;i++){
     6942                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     6943                                for(int e=0;e<nbar;e++){
     6944                                        IssmDouble alpha;
     6945                                        IssmDouble delPhi,delLambda;
     6946                                        /*recover info for this element and vertex:*/
     6947                                        IssmDouble late= slgeom->latbarycentre[l][e];
     6948                                        IssmDouble longe= slgeom->longbarycentre[l][e];
     6949                                        late=late/180*M_PI;
     6950                                        longe=longe/180*M_PI;
     6951
     6952                                        lati=latitude[i];
     6953                                        longi=longitude[i];
     6954
    71816955                                        if(horiz){
    7182                                         /*Compute azimuths*/
     6956                                                /*Compute azimuths*/
    71836957                                                dx=cos(lati)*sin(late)-sin(lati)*cos(late)*cos(longe-longi);
    71846958                                                dy=sin(longe-longi)*cos(late);
    71856959                                                //angle between horiz motion and North, remapped from a double on [0,2*pi] to a int [0,intmax]
    7186                                                 AzimIndex[l][av*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
     6960                                                AzimIndex[l][i*nbar+e]=reCast<int,IssmDouble>(intmax*(atan2(dy,dx)/2/M_PI));
    71876961                                        }
    71886962
    7189                                 /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
    7190                                 delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
    7191                                 alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
    7192                                 doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
    7193                                 index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
    7194 
    7195                                 if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
    7196                                 if (index==M-1){ //avoids out of bound case
    7197                                         index-=1;
    7198                                         lincoef=1;
     6963                                        /*Compute alpha angle between centroid and current vertex and index into precomputed tables: */
     6964                                        delPhi=fabs(lati-late); delLambda=fabs(longi-longe); if (delLambda>M_PI)delLambda=2*M_PI-delLambda;
     6965                                        alpha=2.*asin(sqrt(pow(sin(delPhi/2.0),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2.0),2.0)));
     6966                                        doubleindex=alpha/M_PI*reCast<IssmDouble,int>(M-1); //maps 0<alpha<PI on [0:M-1]
     6967                                        index=reCast<int,IssmDouble>(doubleindex); //truncates doubleindex to integer part
     6968
     6969                                        if ((doubleindex-index)>=0.5) index+=1; //nearest neighbour
     6970                                        if (index==M-1){ //avoids out of bound case
     6971                                                index-=1;
     6972                                                lincoef=1;
     6973                                        }
     6974                                        AlphaIndex[l][i*nbar+e]=index;
    71996975                                }
    7200                                 AlphaIndex[l][av*nbar+e]=index;
    7201                         //}
    7202                         //av+=1;
    7203                         }
    7204 
     6976                        }
    72056977                }
    72066978        }
    … …  
    72086980        /*Save all these arrayins for each element:*/
    72096981        for (int l=0;l<SLGEOM_NUMLOADS;l++){
    7210                 this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*n_activevertices);
    7211                 if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*n_activevertices);
     6982                this->inputs->SetIntArrayInput(slgeom->AlphaIndexEnum(l),this->lid,AlphaIndex[l],slgeom->nbar[l]*3);
     6983                if(horiz) this->inputs->SetIntArrayInput(slgeom->AzimuthIndexEnum(l),this->lid,AzimIndex[l],slgeom->nbar[l]*3);
    72126984        }
    72136985        /*}}}*/
    … …  
    75137285        barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp);
    75147286
    7515         /*Free resources*/
     7287        /*Free ressources*/
    75167288        xDelete<IssmDouble>(areae);
    75177289
    … …  
    75557327
    75567328                for(int i=0;i<NUMVERTICES;i++) slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid]=loadweightsocean[i];
    7557 
    75587329                #ifdef _ISSM_DEBUG_ /*{{{*/
    75597330                /*Inform mask: */
    … …  
    76667437                oceanaverage+=sealevelpercpu[this->vertices[i]->lid]*slgeom->LoadWeigths[SLGEOM_OCEAN][i][this->lid];
    76677438        }
    7668 
     7439        #ifdef _ISSM_DEBUG_
     7440        this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
     7441        #endif
    76697442        oceanarea=slgeom->LoadArea[SLGEOM_OCEAN][this->lid];
    7670         oceanaverage*=rho_water*oceanarea;
    76717443
    76727444        /*add ocean average in the global sealevelloads vector:*/
    76737445        if(slgeom->issubelement[SLGEOM_OCEAN][this->lid]){
    76747446                int intj=slgeom->subelementmapping[SLGEOM_OCEAN][this->lid];
    7675                 loads->vsubsealevelloads->SetValue(intj,oceanaverage,INS_VAL);
     7447                loads->vsubsealevelloads->SetValue(intj,oceanaverage*rho_water*oceanarea,INS_VAL);
    76767448                loads->vsealevelloads->SetValue(this->sid,0.,INS_VAL);
    76777449        }
    7678         else loads->vsealevelloads->SetValue(this->sid,oceanaverage,INS_VAL);
    7679 
    7680         #ifdef _ISSM_DEBUG_
    7681         this->AddInput(SealevelBarystaticOceanLoadEnum,&oceanaverage,P0Enum);
    7682         #endif
     7450        else loads->vsealevelloads->SetValue(this->sid,oceanaverage*rho_water*oceanarea,INS_VAL);
    76837451
    76847452        /*add ocean area into a global oceanareas vector:*/
    … …  
    76997467        IssmDouble* G=NULL;
    77007468        IssmDouble* Grot=NULL;
    7701         IssmDouble* rslfield=NULL;
    77027469        DoubleVecParam* parameter;
    77037470        bool computefuture=false;
     7471        int spatial_component=0;
    77047472
    77057473        bool sal = false;
    … …  
    77207488                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    77217489
     7490                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7491                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    77227492                if (rotation)   this->inputs->GetArrayPtr(SealevelchangeGrotEnum,this->lid,&Grot,&size);
    77237493
    7724                 rslfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=false);
    7725                 this->SealevelchangeCollectGrdfield(sealevelpercpu,rslfield,slgeom,nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    7726 
     7494                this->SealevelchangeGxL(sealevelpercpu, spatial_component=0,AlphaIndex, AlphaIndexsub, NULL, NULL, G, Grot, loads, polarmotionvector, slgeom, nel,percpu=true,SealevelchangeViscousRSLEnum,computefuture=false);
    77277495        }
    77287496
    … …  
    77387506        int nel,nbar;
    77397507        bool sal = false;
     7508        int* AlphaIndex=NULL;
     7509        int* AzimIndex=NULL;
     7510        int* AlphaIndexsub[SLGEOM_NUMLOADS];
     7511        int* AzimIndexsub[SLGEOM_NUMLOADS];
    77407512        int spatial_component=0;
    77417513        IssmDouble* G=NULL;
    … …  
    77467518        IssmDouble* GNrot=NULL;
    77477519        IssmDouble* GErot=NULL;
    7748         IssmDouble* grdfield=NULL;
    77497520
    77507521        DoubleVecParam* parameter;
    … …  
    77677538
    77687539        if(sal){
     7540                this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
     7541                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
     7542
    77697543                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeGViscoElasticEnum)); _assert_(parameter);
    77707544                parameter->GetParameterValueByPointer((IssmDouble**)&G,NULL);
    … …  
    77757549
    77767550                        if(horiz){
     7551                                this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
     7552                                for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
     7553
    77777554                                parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelchangeHViscoElasticEnum)); _assert_(parameter);
    77787555                                parameter->GetParameterValueByPointer((IssmDouble**)&GH,NULL);
    … …  
    77877564                        }
    77887565                }
    7789                 //Relative sea level convolution
    7790                 grdfield=this->SealevelchangeGxL(G,Grot,loads,polarmotionvector,slgeom,nel,computefuture=true);
    7791                 this->SealevelchangeCollectGrdfield(&RSLGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
     7566
     7567                this->SealevelchangeGxL(&RSLGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL,G, Grot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousRSLEnum,computefuture=true);
    77927568
    77937569                if(elastic){
    7794                         //Bedrock Uplift
    7795                         grdfield=this->SealevelchangeGxL(GU,GUrot,loads,polarmotionvector,slgeom,nel,computefuture=true);
    7796                         this->SealevelchangeCollectGrdfield(&UGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
    7797 
     7570                        this->SealevelchangeGxL(&UGrd[0],spatial_component=0, AlphaIndex, AlphaIndexsub,NULL, NULL, GU, GUrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousUEnum,computefuture=true);
    77987571                        if(horiz){
    7799                                 //Bedrock North displacement
    7800                                 grdfield=this->SealevelchangeHorizGxL(spatial_component=1,GH,GNrot,loads,polarmotionvector,slgeom,nel,computefuture=true);
    7801                                 this->SealevelchangeCollectGrdfield(&NGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
    7802 
    7803                                 //Bedrock East displacement
    7804                                 grdfield=this->SealevelchangeHorizGxL(spatial_component=2,GH,GErot,loads,polarmotionvector,slgeom,nel,computefuture=true);
    7805                                 this->SealevelchangeCollectGrdfield(&EGrd[0],grdfield,slgeom,nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
     7572                                this->SealevelchangeGxL(&NGrd[0],spatial_component=1,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GNrot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousNEnum,computefuture=true);
     7573                                this->SealevelchangeGxL(&EGrd[0],spatial_component=2,AlphaIndex, AlphaIndexsub,AzimIndex,AzimIndexsub,GH, GErot, loads, polarmotionvector, slgeom, nel,percpu=false,SealevelchangeViscousEEnum,computefuture=true);
    78067574                        }
    78077575                }
    … …  
    78287596
    78297597} /*}}}*/
    7830 IssmDouble*       Tria::SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/
     7598void       Tria::SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
    78317599
    78327600        //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    7833         int* AlphaIndex=NULL;
    7834         int* AlphaIndexsub[SLGEOM_NUMLOADS];
    7835         int* activevertices=NULL;
     7601
    78367602        IssmDouble* grdfield=NULL;
    7837         int i,e,l,t,a, index, nbar, size, av,ae,b,c;
     7603        int i,e,l,t,a, index, nbar;
    78387604        bool rotation=false;
     7605        IssmDouble* Centroid_loads=NULL;
     7606        IssmDouble* Centroid_loads_copy=NULL;
     7607        IssmDouble* Subelement_loads[SLGEOM_NUMLOADS];
     7608        IssmDouble* Subelement_loads_copy[SLGEOM_NUMLOADS];
     7609        IssmDouble* horiz_projection=NULL;
     7610        IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
    78397611        int nt=1; //important, ensures there is a defined value if computeviscous is false
    7840         int n_activevertices=0;
    78417612
    78427613        //viscous
    … …  
    78447615        int viscousindex=0; //important
    78457616        int viscousnumsteps=1; //important
    7846 
    7847         this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    7848         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    7849 
    7850         //Get green functions indexing & geometry
    7851         this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices); //the order in which the vertices appear here should be the same as in slgeom->lids
    7852         this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
    7853         for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    7854 
    7855         if(computeviscous){
    7856                 this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    7857                 this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
    7858                 if(computefuture) {
    7859                         nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity
    7860                         if (nt>viscousnumsteps) nt=viscousnumsteps;
    7861                 }
    7862                 else nt=1;
    7863         }
    7864         //allocate
    7865         grdfield=xNewZeroInit<IssmDouble>(3*nt);
    7866 
    7867         //early return
    7868         if (n_activevertices==0) return grdfield;
    7869 
    7870         if(rotation){ //add rotational feedback
    7871                 for(av=0;av<n_activevertices;av++) { //vertices
    7872                         i=activevertices[av];
    7873                         //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7874                         b=i*nt;
    7875                         for (int m=0;m<3;m++){ //polar motion components
    7876                                 for(t=0;t<nt;t++){ //time
    7877                                         int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
    7878                                         grdfield[b+t]+=Grot[index]*polarmotionvector[m];
    7879                                 }
    7880                         }
    7881                 }
    7882         }
    7883 
    7884         //Convolution
    7885         for(av=0;av<n_activevertices;av++) { /*{{{*/
    7886                 i=activevertices[av];
    7887                 //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7888                 b=i*nt;
    7889                 c=av*nel;
    7890                 for(ae=0;ae<loads->nactiveloads;ae++){
    7891                         e=loads->combined_loads_index[ae];
    7892                         a=AlphaIndex[c+e]*viscousnumsteps;
    7893                         for(t=0;t<nt;t++){
    7894                                 grdfield[b+t]+=G[a+t]*loads->combined_loads[ae];
    7895                         }
    7896                 }
    7897                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    7898                         nbar=slgeom->nbar[l];
    7899                         c=av*nbar;
    7900                         for (ae=0;ae<loads->nactivesubloads[l];ae++){
    7901                                 e=loads->combined_subloads_index[l][ae];
    7902                                 a=AlphaIndexsub[l][c+e]*viscousnumsteps;
    7903                                 for(t=0;t<nt;t++){
    7904                                         grdfield[b+t]+=G[a+t]*loads->combined_subloads[l][ae];
    7905                                 }
    7906                         }
    7907                 }
    7908                 //av+=1;
    7909         } /*}}}*/
    7910 
    7911         return grdfield;
    7912 
    7913 } /*}}}*/
    7914 IssmDouble*       Tria::SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture) { /*{{{*/
    7915 
    7916         //This function performs the actual convolution between Green functions and surface Loads for a particular grd field
    7917         int* AlphaIndex=NULL;
    7918         int* AzimIndex=NULL;
    7919         int* AlphaIndexsub[SLGEOM_NUMLOADS];
    7920         int* AzimIndexsub[SLGEOM_NUMLOADS];
    7921         int* activevertices = NULL;
    7922         IssmDouble* grdfield=NULL;
    7923         int i,e,l,t,a,b,c, index, nbar, av, ae,n_activevertices, size;
    7924         bool rotation=false;
    7925         IssmDouble* projected_loads=NULL;
    7926         IssmDouble* projected_subloads[SLGEOM_NUMLOADS];
    7927         IssmDouble* horiz_projection=NULL;
    7928         IssmDouble* horiz_projectionsub[SLGEOM_NUMLOADS];
    7929         int nt=1; //important, ensures there is a defined value if computeviscous is false
    7930 
    7931         //viscous
    7932         bool computeviscous=false;
    7933         int viscousindex=0; //important
    7934         int viscousnumsteps=1; //important
    7935 
    7936         //Get green functions indexing & geometry
    7937         this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
    7938         this->inputs->GetIntArrayPtr(SealevelchangeAlphaIndexEnum,this->lid,&AlphaIndex,&size);
    7939         for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AlphaIndexEnum(l),this->lid,&AlphaIndexsub[l],&size);
    7940         this->inputs->GetIntArrayPtr(SealevelchangeAzimuthIndexEnum,this->lid,&AzimIndex,&size);
    7941         for (int l=0;l<SLGEOM_NUMLOADS;l++) this->inputs->GetIntArrayPtr(slgeom->AzimuthIndexEnum(l),this->lid,&AzimIndexsub[l],&size);
    7942 
    7943         //First, figure out how many time steps to compute grdfield for
    7944         this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    7945         this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    7946         if(computeviscous){
    7947                 this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    7948                 this->parameters->FindParam(&viscousindex,SealevelchangeViscousIndexEnum);
    7949                 if(computefuture) {
    7950                         nt=viscousnumsteps-viscousindex+2; //number of time steps remaining to reach final_time, +1 is sufficient with no adaptative time stepping, +2 necessary otherwise; we assume the safe choice here for the sake of simplicity
    7951                         if (nt>viscousnumsteps) nt=viscousnumsteps;
    7952                 }
    7953                 else nt=1;
    7954         }
    7955         //allocate
    7956         grdfield=xNewZeroInit<IssmDouble>(3*nt);
    7957         if (n_activevertices==0) return grdfield;
    7958 
    7959         if(rotation){ //add rotational feedback
    7960                 for(av=0;av<n_activevertices;av++) { //vertices
    7961                         i=activevertices[av];
    7962                         //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7963                         for (int m=0;m<3;m++){ //polar motion components
    7964                                 for(t=0;t<nt;t++){ //time
    7965                                         int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
    7966                                         grdfield[i*nt+t]+=Grot[index]*polarmotionvector[m];
    7967                                 }
    7968                         }
    7969                         //}
    7970                 }
    7971         }
    7972 
    7973         //Initialize projection vectors
    7974         horiz_projection=xNewZeroInit<IssmDouble>(loads->nactiveloads);
    7975         projected_loads=xNewZeroInit<IssmDouble>(loads->nactiveloads);
    7976         for(l=0;l<SLGEOM_NUMLOADS;l++){
    7977                 //nbar=slgeom->nbar[l];
    7978                 projected_subloads[l]=xNewZeroInit<IssmDouble>(loads->nactivesubloads[l]);
    7979                 horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(loads->nactivesubloads[l]);
    7980         }
    7981 
    7982 
    7983         //Convolution
    7984         //av=0;
    7985         for(av=0;av<n_activevertices;av++) { //vertices
    7986                 i=activevertices[av];
    7987                 //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    7988                 b=i*nt;
    7989 
    7990                 //GxL needs to be projected on the right axis before summation into the grd field
    7991                 //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
    7992 
    7993                 //get projection
    7994                 if (spatial_component==1){ //north
    7995                         for(ae=0;ae<loads->nactiveloads;ae++){
    7996                                 e=loads->combined_loads_index[ae];
    7997                                 horiz_projection[ae]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[av*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
    7998                         }
    7999                         for(l=0;l<SLGEOM_NUMLOADS;l++){
    8000                                 nbar=slgeom->nbar[l];
    8001                                 for(ae=0;ae<loads->nactivesubloads[l];ae++){
    8002                                         e=loads->combined_subloads_index[l][ae];
    8003                                         horiz_projectionsub[l][ae]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][av*nbar+e])/65535.0);
    8004                                 }
    8005                         }
    8006                 }
    8007                 else if (spatial_component==2){ //east
    8008                         for(ae=0;ae<loads->nactiveloads;ae++){
    8009                                 e=loads->combined_loads_index[ae];
    8010                                 horiz_projection[ae]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[av*nel+e])/65535.0);
    8011                         }
    8012                         for(l=0;l<SLGEOM_NUMLOADS;l++){
    8013                                 nbar=slgeom->nbar[l];
    8014                                 for(ae=0;ae<loads->nactivesubloads[l];ae++){
    8015                                         e=loads->combined_subloads_index[l][ae];
    8016                                         horiz_projectionsub[l][ae]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][av*nbar+e])/65535.0);
    8017                                 }
    8018                         }
    8019                 }
    8020 
    8021                 //project load in the right direction
    8022                 for (ae=0;ae<loads->nactiveloads;ae++){
    8023                         projected_loads[ae]=loads->combined_loads[ae]*horiz_projection[ae];
    8024                 }
    8025                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    8026                         nbar=slgeom->nbar[l];
    8027                         for(ae=0;ae<loads->nactivesubloads[l];ae++){
    8028                                 projected_subloads[l][ae]=loads->combined_subloads[l][ae]*horiz_projectionsub[l][ae];
    8029                         }
    8030                 }
    8031 
    8032                 //do the convolution
    8033                 c=av*nel;
    8034                 for(ae=0;ae<loads->nactiveloads;ae++){
    8035                         e=loads->combined_loads_index[ae];
    8036                         a=AlphaIndex[c+e]*viscousnumsteps;
    8037                         for(t=0;t<nt;t++){
    8038                                 grdfield[b+t]+=G[a+t]*projected_loads[ae];
    8039                         }
    8040                 }
    8041                 for(l=0;l<SLGEOM_NUMLOADS;l++){
    8042                         nbar=slgeom->nbar[l];
    8043                         c=av*nbar;
    8044                         for(ae=0;ae<loads->nactivesubloads[l];ae++){
    8045                                 e=loads->combined_subloads_index[l][ae];
    8046                                 a=AlphaIndexsub[l][c+e]*viscousnumsteps;
    8047                                 for(t=0;t<nt;t++){
    8048                                         grdfield[b+t]+=G[a+t]*projected_subloads[l][ae];
    8049                                 }
    8050                         }
    8051                 }
    8052                 //av+=1;
    8053         } /*}}}*/
    8054 
    8055 
    8056         //free resources
    8057         xDelete<IssmDouble>(horiz_projection);
    8058         xDelete<IssmDouble>(projected_loads);
    8059         for(l=0;l<SLGEOM_NUMLOADS;l++) {
    8060                 xDelete<IssmDouble>(projected_subloads[l]);
    8061                 xDelete<IssmDouble>(horiz_projectionsub[l]);
    8062         }
    8063         return grdfield;
    8064 
    8065 } /*}}}*/
    8066 
    8067 
    8068 void       Tria::SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture) { /*{{{*/
    8069 
    8070         //This function aligns grdfield with the requested output format: in a size 3 vector or in a size numberofvertices vector
    8071         // if compute viscous is on, we also interpolate the field timewise given the current timestepping as well as collect viscous deformation and update the viscous deformation time series for future time steps
    8072         int i,e,l,t,a, index, nbar, av, n_activevertices;
    8073         int nt=1;
    8074 
    8075 
    8076         //viscous
    8077         bool computeviscous=false;
    8078         int viscousindex=0; //important
    8079         int viscousnumsteps=1; //important
    8080         int* activevertices = NULL;
    80817617        IssmDouble* viscousfield=NULL;
    80827618        IssmDouble* grdfieldinterp=NULL;
    … …  
    80867622        IssmDouble  timeacc;
    80877623
    8088         //parameters & initialization
    80897624        this->parameters->FindParam(&computeviscous,SolidearthSettingsViscousEnum);
    8090         this->inputs->GetIntArrayPtr(SealevelchangeConvolutionVerticesEnum,this->lid,&activevertices,&n_activevertices);
    8091 
     7625        this->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
    80927626        if(computeviscous){
    80937627                this->parameters->FindParam(&viscousnumsteps,SealevelchangeViscousNumStepsEnum);
    … …  
    81047638                else nt=1;
    81057639        }
    8106        
    8107 
    8108         if(!computeviscous){ /*{{{*/
    8109                 /*elastic or self attraction only case
    8110                   store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
    8111                 if(percpu){
    8112                         for(av=0;av<n_activevertices;av++) { //vertices
    8113                                 i=activevertices[av];
    8114                                 //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    8115                                 grdfieldout[this->vertices[i]->lid]=grdfield[i];
    8116                                 //}
    8117                         }
    8118                 }
    8119                 else{
    8120                         for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i];
    8121                 }
    8122                 //free resources
    8123                 xDelete<IssmDouble>(grdfield);
    8124                 return;
    8125         }
    8126         else { //viscous case
     7640        //allocate
     7641        grdfield=xNewZeroInit<IssmDouble>(3*nt);
     7642
     7643        if(rotation){ //add rotational feedback
     7644                for(int i=0;i<NUMVERTICES;i++){ //vertices
     7645                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
     7646                                for (int m=0;m<3;m++){ //polar motion components
     7647                                        for(t=0;t<nt;t++){ //time
     7648                                                int index=m*3*viscousnumsteps+i*viscousnumsteps+t;
     7649                                                grdfield[i*nt+t]+=Grot[index]*polarmotionvector[m];
     7650                                        }
     7651                                }
     7652                        }
     7653                }
     7654        }
     7655
     7656        //Determine loads /*{{{*/
     7657        Centroid_loads=xNewZeroInit<IssmDouble>(nel);
     7658        for (e=0;e<nel;e++){
     7659                Centroid_loads[e]=loads->loads[e];
     7660        }
     7661        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7662                nbar=slgeom->nbar[l];
     7663                Subelement_loads[l]=xNewZeroInit<IssmDouble>(nbar);
     7664                for (e=0;e<nbar;e++){
     7665                        Subelement_loads[l][e]=(loads->subloads[l][e]);
     7666                }
     7667        }
     7668        if(loads->sealevelloads){
     7669                for (e=0;e<nel;e++){
     7670                        Centroid_loads[e]+=(loads->sealevelloads[e]);
     7671                }
     7672                nbar=slgeom->nbar[SLGEOM_OCEAN];
     7673                for (e=0;e<nbar;e++){
     7674                        Subelement_loads[SLGEOM_OCEAN][e]+=(loads->subsealevelloads[e]);
     7675                }
     7676        }
     7677
     7678        //Copy loads if dealing with a horizontal component: the result will need to be projected against the North or East axis for each vertex
     7679        if (spatial_component!=0){
     7680                horiz_projection=xNewZeroInit<IssmDouble>(3*nel);
     7681                Centroid_loads_copy=xNewZeroInit<IssmDouble>(nel);
     7682                for (e=0;e<nel;e++){
     7683                        Centroid_loads_copy[e]=Centroid_loads[e];
     7684                }
     7685
     7686                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7687                        nbar=slgeom->nbar[l];
     7688                        Subelement_loads_copy[l]=xNewZeroInit<IssmDouble>(nbar);
     7689                        horiz_projectionsub[l]=xNewZeroInit<IssmDouble>(3*nbar);
     7690                        for (e=0;e<nbar;e++){
     7691                                Subelement_loads_copy[l][e]=Subelement_loads[l][e];
     7692                        }
     7693                }
     7694        }
     7695        /*}}}*/
     7696
     7697        //Convolution
     7698        for(i=0;i<NUMVERTICES;i++) { /*{{{*/
     7699                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7700
     7701                if (spatial_component!=0){ //horizontals /*{{{*/
     7702                        //GxL needs to be projected on the right axis before summation into the grd field
     7703                        //here we apply the projection scalar to the load prior to the actual convolution loop for more efficiency
     7704                        if (spatial_component==1){ //north
     7705                                for (e=0;e<nel;e++){
     7706                                        horiz_projection[i*nel+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0); // 65535=2^16-1 = max value of 16 bits unsigned int
     7707                                }
     7708                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7709                                        nbar=slgeom->nbar[l];
     7710                                        for (e=0;e<nbar;e++){
     7711                                                horiz_projectionsub[l][i*nbar+e]=cos(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7712                                        }
     7713                                }
     7714                        }
     7715                        else if (spatial_component==2){ //east
     7716                                for (e=0;e<nel;e++){
     7717                                        horiz_projection[i*nel+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndex[i*nel+e])/65535.0);
     7718                                }
     7719                                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7720                                        nbar=slgeom->nbar[l];
     7721                                        for (e=0;e<nbar;e++){
     7722                                                horiz_projectionsub[l][i*nbar+e]=sin(2.0*M_PI*reCast<IssmDouble,int>(AzimIndexsub[l][i*nbar+e])/65535.0);;
     7723                                        }
     7724                                }
     7725                        }
     7726                        for (e=0;e<nel;e++) Centroid_loads[e]=Centroid_loads_copy[e]*horiz_projection[i*nel+e];
     7727                        for(l=0;l<SLGEOM_NUMLOADS;l++){
     7728                                nbar=slgeom->nbar[l];
     7729                                for (e=0;e<nbar;e++){
     7730                                        Subelement_loads[l][e]=Subelement_loads_copy[l][e]*horiz_projectionsub[l][i*nbar+e];
     7731                                }
     7732                        }
     7733                } /*}}}*/
     7734
     7735                for (e=0;e<nel;e++){
     7736                        for(t=0;t<nt;t++){
     7737                                a=AlphaIndex[i*nel+e];
     7738                                grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Centroid_loads[e];
     7739                        }
     7740                }
     7741                for(l=0;l<SLGEOM_NUMLOADS;l++){
     7742                        nbar=slgeom->nbar[l];
     7743                        for (e=0;e<nbar;e++){
     7744                                for(t=0;t<nt;t++){
     7745                                        a=AlphaIndexsub[l][i*nbar+e];
     7746                                        grdfield[i*nt+t]+=G[a*viscousnumsteps+t]*Subelement_loads[l][e];
     7747                                }
     7748                        }
     7749                }
     7750        } /*}}}*/
     7751
     7752
     7753
     7754        if(computeviscous){ /*{{{*/
    81277755                // we need to do up to 3 things (* = only if computefuture)
    8128                 // 1: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time]
    8129                 // 2*: add new grdfield contribution to the viscous stack for future time steps
     7756                // 1*: add new grdfield contribution to the viscous stack for future time steps
     7757                // 2: collect viscous grdfield from past loads due at present-day and add it to grdfield[current_time]
    81307758                // 3*: subtract from viscous stack the grdfield that has already been accounted for so we don't add it again at the next time step
    8131 
    8132                 /*update grdfield at present time using viscous stack at present time: */
    8133                 for(av=0;av<n_activevertices;av++) { //vertices
    8134                         i=activevertices[av];
    8135                         //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    8136                         grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex];
    8137                 }
    8138 
    81397759
    81407760                /* Map new grdfield generated by present-day loads onto viscous time vector*/
    81417761                if(computefuture){
    81427762                        //viscousindex time and first time step of grdfield coincide, so just copy that value
    8143                         for(av=0;av<n_activevertices;av++) { //vertices
    8144                                 i=activevertices[av];
    8145                                 //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    8146                                 grdfieldinterp[i*viscousnumsteps+viscousindex]=grdfield[i*nt+0];
     7763                        for(int i=0;i<NUMVERTICES;i++){
     7764                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7765                                grdfieldinterp[i*viscousnumsteps+viscousindex]=  grdfield[i*nt+0];
    81477766                        }
    81487767                        if(viscoustimes[viscousindex]<final_time){
    81497768                                //And interpolate the rest of the points in the future
    81507769                                lincoeff=(viscoustimes[viscousindex+1]-viscoustimes[viscousindex])/timeacc;
    8151                                 for(av=0;av<n_activevertices;av++) { //vertices
    8152                                         i=activevertices[av];
    8153                                         //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
    8154                                         int i_time1= i*nt-viscousindex;
    8155                                         int i_time2= i*viscousnumsteps;
    8156                                         for(int t=viscousindex+1;t<viscousnumsteps;t++){
    8157                                                 grdfieldinterp[i_time2+t] = (1-lincoeff)*grdfield[i_time1+t-1]
    8158                                                                           +    lincoeff *grdfield[i_time1+t]
    8159                                                                           +          viscousfield[i_time2+t];
    8160                                                 /*update viscous stack with future deformation from present load: */
    8161                                                 viscousfield[i_time2+t]=grdfieldinterp[i_time2+t]
    8162                                                                        -grdfieldinterp[i_time2+viscousindex];
     7770                                for(int t=viscousindex+1;t<viscousnumsteps;t++){
     7771                                        for(int i=0;i<NUMVERTICES;i++){
     7772                                                if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7773                                                grdfieldinterp[i*viscousnumsteps+t] = (1-lincoeff)*grdfield[i*nt+(t-viscousindex-1)]
     7774                                                                                         +lincoeff*grdfield[i*nt+(t-viscousindex)];
    81637775                                        }
    81647776                                }
    81657777                        }
     7778                }
     7779
     7780                /*update grdfield at present time using viscous stack at present time: */
     7781                for(int i=0;i<NUMVERTICES;i++){
     7782                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7783                        grdfield[i*nt+0]+=viscousfield[i*viscousnumsteps+viscousindex];
     7784                }
     7785
     7786                /*update viscous stack with future deformation from present load: */
     7787                if(computefuture){
     7788                        for(int t=viscousnumsteps-1;t>=viscousindex;t--){ //we need to go backwards so as not to zero out viscousfield[i*viscousnumsteps+viscousindex] until the end
     7789                                for(int i=0;i<NUMVERTICES;i++){
     7790                                        if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7791                                        //offset viscousfield to remove all deformation that has already been added
     7792                                        viscousfield[i*viscousnumsteps+t]+=grdfieldinterp[i*viscousnumsteps+t]
     7793                                                                          -grdfieldinterp[i*viscousnumsteps+viscousindex]
     7794                                                                          -viscousfield[i*viscousnumsteps+viscousindex];
     7795                                }
     7796                        }
    81667797                        /*Save viscous stack now that we updated the values:*/
    81677798                        this->inputs->SetArrayInput(viscousenum,this->lid,viscousfield,3*viscousnumsteps);
    8168                 }
    8169 
    8170                 /*store values computed on vertices*/
    8171                 if(percpu){
    8172                         for(av=0;av<n_activevertices;av++) { //vertices
    8173                                 i=activevertices[av];
    8174                                 //if(slgeom->lids[this->vertices[i]->lid]!=this->lid)continue;
     7799
     7800                        /*Free resources:*/
     7801                        xDelete<IssmDouble>(grdfieldinterp);
     7802                }
     7803        }
     7804        /*}}}*/
     7805
     7806        /*store values computed on vertices, but don't repeat the computation if another element already computed it!:*/
     7807        if(percpu){
     7808                for(i=0;i<NUMVERTICES;i++){
     7809                        if(slgeom->lids[this->vertices[i]->lid]==this->lid){
    81757810                                grdfieldout[this->vertices[i]->lid]=grdfield[i*nt+0];
    8176                                 //}
    8177                         }
    8178                 }
    8179                 else{
    8180                         for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
    8181                 }
    8182                 //free resources
    8183                 xDelete<IssmDouble>(grdfield);
     7811                        }
     7812                }
     7813        }
     7814        else{
     7815                for(i=0;i<NUMVERTICES;i++) grdfieldout[i]=grdfield[i*nt+0];
     7816        }
     7817        //free resources
     7818        xDelete<IssmDouble>(grdfield);
     7819        xDelete<IssmDouble>(Centroid_loads);
     7820        for(l=0;l<SLGEOM_NUMLOADS;l++) xDelete<IssmDouble>(Subelement_loads[l]);
     7821        if (spatial_component!=0){
     7822                xDelete<IssmDouble>(horiz_projection);
     7823                xDelete<IssmDouble>(Centroid_loads_copy);
     7824                for(l=0;l<SLGEOM_NUMLOADS;l++) {
     7825                        xDelete<IssmDouble>(Subelement_loads_copy[l]);
     7826                        xDelete<IssmDouble>(horiz_projectionsub[l]);
     7827                }
     7828        }
     7829        if (computeviscous){
    81847830                xDelete<IssmDouble>(viscoustimes);
    81857831                if (computefuture){
    81867832                        xDelete<IssmDouble>(grdfieldinterp);
    81877833                }
    8188                 /*}}}*/
    8189         }
     7834        }
     7835
    81907836} /*}}}*/
    81917837

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Elements/Tria.h

    r27852 r27956  
    5555                void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
    5656                void                    CalvingRateVonmises();
    57                 void                    CalvingRateVonmisesAD();
    5857                void                    CalvingRateTest();
    5958                void        CalvingCrevasseDepth();
    … …  
    6362                void                    CalvingMeltingFluxLevelset();
    6463                void                    CalvingRateParameterization();
    65                 void                    CalvingRateCalvingMIP();
    6664                IssmDouble  CharacteristicLength(void);
    6765                void        ComputeBasalStress(void);
    … …  
    128126                void        MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
    129127                IssmDouble  Misfit(int modelenum,int observationenum,int weightsenum);
     128                void            MisfitAccumulate(int modelenum,IssmDouble dt);
     129                IssmDouble  MisfitAnnual(int modelenum,int observationenum,int weightsenum,IssmDouble annual_dt);
    130130                IssmDouble  MisfitArea(int weightsenum);
    131131                int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
    … …  
    158158                IssmDouble  TotalGroundedBmb(bool scaled);
    159159                IssmDouble  TotalSmb(bool scaled);
    160                 IssmDouble  TotalSmbMelt(bool scaled);
    161                 IssmDouble  TotalSmbRefreeze(bool scaled);
    162160                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    163161                int         UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
    … …  
    176174                #ifdef _HAVE_SEALEVELCHANGE_
    177175                void       GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y);
    178                 void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids, int* vcount);
     176                void       SealevelchangeGeometryInitial(IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze, IssmDouble* areae, int* lids);
    179177                void       SealevelchangeGeometrySubElementKernel(SealevelGeometry* slgeom);
    180178                void       SealevelchangeGeometrySubElementLoads(SealevelGeometry* slgeom, IssmDouble* areae);
    … …  
    249247                void           UpdateConstraintsExtrudeFromBase(void);
    250248                void           UpdateConstraintsExtrudeFromTop(void);
    251                 IssmDouble*    SealevelchangeGxL(IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture);
    252                 IssmDouble*    SealevelchangeHorizGxL(int spatial_component, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool computefuture);
    253                 void           SealevelchangeCollectGrdfield(IssmDouble* grdfieldout, IssmDouble* grdfield, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
     249                void           SealevelchangeGxL(IssmDouble* grdfieldout, int spatial_component, int* AlphaIndex, int** AlphaIndexsub, int* AzimIndex, int**AzimIndexsub, IssmDouble* G, IssmDouble* Grot, GrdLoads* loads, IssmDouble* polarmotionvector, SealevelGeometry* slgeom, int nel, bool percpu, int viscousenum, bool computefuture);
    254250                /*}}}*/
    255251

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ArrayInput.cpp

    r27308 r27956  
    4444        ArrayInput* output = new ArrayInput(this->numberofelements_local);
    4545
     46        output->N = xNew<int>(this->numberofelements_local);
    4647        xMemCpy<int>(output->N,this->N,this->numberofelements_local);
    4748
     49        output->values = xNew<IssmDouble*>(this->numberofelements_local);
    4850        for(int i=0;i<this->numberofelements_local;i++){
    4951                if(this->values[i]){

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.cpp

    r27703 r27956  
    227227}
    228228/*}}}*/
    229 void ControlInput::AverageAndReplace(void){/*{{{*/
    230         this->values->AverageAndReplace();
    231 }
    232 /*}}}*/
    233229TriaInput* ControlInput::GetTriaInput(){/*{{{*/
    234230

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/ControlInput.h

    r27703 r27956  
    4848                TriaInput* GetTriaInput();
    4949                PentaInput* GetPentaInput();
    50                 void AverageAndReplace(void);
    5150};
    5251#endif  /* _CONTROLINPUT_H */

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Input.h

    r27691 r27956  
    4848                virtual void   Pow(IssmDouble scale_factor){_error_("Not implemented yet");};
    4949                virtual void   Scale(IssmDouble scale_factor){_error_("Not implemented yet");};
    50                 virtual void   AverageAndReplace(void){_error_("Not implemented yet");};
    5150
    5251                virtual int  GetResultArraySize(void){_error_("Not implemented yet");};

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.cpp

    r27691 r27956  
    302302}
    303303/*}}}*/
    304 void Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/
     304void     Inputs::Shift(int xenum, IssmDouble alpha){/*{{{*/
    305305
    306306        _assert_(this);
    … …  
    314314        /*Shift: */
    315315        this->inputs[index_x]->Shift(alpha);
    316 }
    317 /*}}}*/
    318 void Inputs::AverageAndReplace(int inputenum){/*{{{*/
    319 
    320         _assert_(this);
    321 
    322         /*Get indices from enums*/
    323         int index = EnumToIndex(inputenum);
    324         if(!this->inputs[index]) _error_("Input "<<EnumToStringx(inputenum)<<" not found");
    325 
    326         this->inputs[index]->AverageAndReplace();
    327316}
    328317/*}}}*/

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/Inputs.h

    r27691 r27956  
    5050                void     AXPY(IssmDouble alpha, int xenum, int yenum);
    5151                void     Shift(int inputenum, IssmDouble alpha);
    52                 void     AverageAndReplace(int inputenum);
    5352                void     DeepEcho(void);
    5453                void     DeepEcho(int enum_in);

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/IntArrayInput.cpp

    r27308 r27956  
    4444        IntArrayInput* output = new IntArrayInput(this->numberofelements_local);
    4545
     46        output->N = xNew<int>(this->numberofelements_local);
    4647        xMemCpy<int>(output->N,this->N,this->numberofelements_local);
    4748
     49        output->values = xNew<int*>(this->numberofelements_local);
    4850        for(int i=0;i<this->numberofelements_local;i++){
    4951                if(this->values[i]){

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TransientInput.cpp

    r27822 r27956  
    435435        /*First, recover current time from parameters: */
    436436        bool linear_interp,average,cycle;
    437         int  timestepping;
    438437        IssmDouble dt;
    439438        this->parameters->FindParam(&linear_interp,TimesteppingInterpForcingEnum);
    … …  
    441440        this->parameters->FindParam(&cycle,TimesteppingCycleForcingEnum);
    442441        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
    443         this->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
     442
     443        /*Change input time if we cycle through the forcing*/
     444        IssmDouble time0 = this->timesteps[0];
     445        IssmDouble time1 = this->timesteps[this->numtimesteps - 1];
     446
     447        /*We need the end time to be the last timestep that would be taken*/
     448        /* i.e., the case where GEMB has time stamps (finer timestep) after the last timestep */
     449        IssmDouble nsteps = reCast<int,IssmDouble>(time1/dt);
     450        if (reCast<IssmDouble>(nsteps)<time1/dt) nsteps=nsteps+1;
     451        time1 = nsteps*dt;
    444452
    445453        if(cycle){
    446 
    447                 /*Change input time if we cycle through the forcing*/
    448                 IssmDouble time0 = this->timesteps[0];
    449                 IssmDouble time1 = this->timesteps[this->numtimesteps - 1];
    450 
    451                 if(timestepping!=AdaptiveTimesteppingEnum){
    452                         /*We need the end time to be the last timestep that would be taken*/
    453                         /* i.e., the case where GEMB has time stamps (finer timestep) after the last timestep */
    454                         /* warning: this assumes dt = constant!!*/
    455                         IssmDouble nsteps = reCast<int,IssmDouble>(time1/dt);
    456                         if (reCast<IssmDouble>(nsteps)<time1/dt) nsteps=nsteps+1;
    457                         time1 = nsteps*dt;
    458                 }
    459454
    460455                /*See by how many intervals we have to offset time*/

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TriaInput.cpp

    r27703 r27956  
    406406}
    407407/*}}}*/
    408 void TriaInput::AverageAndReplace(void){/*{{{*/
    409 
    410         if(this->M!=this->numberofelements_local) _error_("not implemented for P1");
    411 
    412         /*Get local sum and local size*/
    413         IssmDouble sum  = 0.;
    414         int        weight;
    415         for(int i=0;i<this->M*this->N;i++) sum += this->values[i];
    416         weight = this->M*this->N;
    417 
    418         /*Get sum across all procs*/
    419         IssmDouble all_sum;
    420         int        all_weight;
    421         ISSM_MPI_Allreduce((void*)&sum,(void*)&all_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    422         ISSM_MPI_Allreduce((void*)&weight,(void*)&all_weight,1,ISSM_MPI_INT,ISSM_MPI_SUM,IssmComm::GetComm());
    423 
    424         /*Divide by number of procs*/
    425         IssmDouble newvalue = all_sum/reCast<IssmPDouble>(all_weight);
    426 
    427         /*Now replace existing input*/
    428         this->Reset(P0Enum);
    429         for(int i=0;i<this->M*this->N;i++) this->values[i] = newvalue;
    430 }
    431 /*}}}*/
    432408
    433409/*Object functions*/

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/Inputs/TriaInput.h

    r27691 r27956  
    4141                void AXPY(Input* xinput,IssmDouble scalar);
    4242                void Shift(IssmDouble scalar);
    43                 void AverageAndReplace(void);
    4443                void PointWiseMult(Input* xinput);
    4544                void Serve(int numindices,int* indices);

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/classes/classes.h

    r27920 r27956  
    1818#include "./Massfluxatgate.h"
    1919#include "./Misfit.h"
     20#include "./MisfitAnnual.h"
    2021#include "./SealevelGeometry.h"
    2122#include "./GrdLoads.h"
    … …  
    2425#include "./Numberedcostfunction.h"
    2526#include "./Cfsurfacesquare.h"
    26 #include "./Cfsurfacesquaretransient.h"
    2727#include "./Cfdragcoeffabsgrad.h"
    28 #include "./Cfdragcoeffabsgradtransient.h"
    29 #include "./Cfrheologybbarabsgrad.h"
    30 #include "./Cfrheologybbarabsgradtransient.h"
    3128#include "./Cfsurfacelogvel.h"
    3229#include "./Cflevelsetmisfit.h"
    … …  
    9390#include "./Params/GenericParam.h"
    9491#include "./Params/BoolParam.h"
    95 #include "./Params/ControlParam.h"
    9692#include "./Params/DoubleMatParam.h"
    9793#include "./Params/DoubleTransientMatParam.h"

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

    r27920 r27956  
    154154                                /*}}}*/
    155155                        }
     156                        else if (output_definition_enums[i]==MisfitannualEnum){
     157                                /*Deal with misfits: {{{*/
     158
     159                                /*misfit variables: */
     160                                int          nummisfits;
     161                                char**       misfit_name_s                                              = NULL;   
     162                                char**           misfit_definitionstring_s              = NULL;   
     163                                char**       misfit_model_string_s                      = NULL;
     164                                IssmDouble** misfit_observation_s                       = NULL;
     165                                char**           misfit_observation_string_s    = NULL;
     166                                int*         misfit_observation_M_s                     = NULL;
     167                                int*         misfit_observation_N_s                     = NULL;
     168                                int*         misfit_local_s                                     = NULL;
     169                                char**       misfit_timeinterpolation_s = NULL;
     170                                IssmDouble** misfit_weights_s                                   = NULL;
     171                                int*         misfit_weights_M_s                         = NULL;
     172                                int*         misfit_weights_N_s                         = NULL;
     173                                char**       misfit_weights_string_s            = NULL;
     174
     175                                /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/misfit.m): */
     176                                iomodel->FetchMultipleData(&misfit_name_s,&nummisfits,                                                        "md.misfitannual.name");
     177                                iomodel->FetchMultipleData(&misfit_definitionstring_s,&nummisfits,                                            "md.misfitannual.definitionstring");
     178                                iomodel->FetchMultipleData(&misfit_model_string_s,&nummisfits,                                                "md.misfitannual.model_string");
     179                                iomodel->FetchMultipleData(&misfit_observation_s,&misfit_observation_M_s,&misfit_observation_N_s,&nummisfits, "md.misfitannual.observation");
     180                                iomodel->FetchMultipleData(&misfit_observation_string_s,&nummisfits,                                          "md.misfitannual.observation_string");
     181                                iomodel->FetchMultipleData(&misfit_timeinterpolation_s,&nummisfits,                                           "md.misfitannual.timeinterpolation");
     182                                iomodel->FetchMultipleData(&misfit_local_s,&nummisfits,                                                       "md.misfitannual.local");
     183                                iomodel->FetchMultipleData(&misfit_weights_s,&misfit_weights_M_s,&misfit_weights_N_s,&nummisfits,             "md.misfitannual.weights");
     184                                iomodel->FetchMultipleData(&misfit_weights_string_s,&nummisfits,                                              "md.misfitannual.weights_string");
     185
     186                                for(j=0;j<nummisfits;j++){
     187
     188                                        int obs_vector_type=0;
     189                                        if ((misfit_observation_M_s[j]==iomodel->numberofvertices) || (misfit_observation_M_s[j]==iomodel->numberofvertices+1)){
     190                                                obs_vector_type=1;
     191                                        }
     192                                        else if ((misfit_observation_M_s[j]==iomodel->numberofelements) || (misfit_observation_M_s[j]==iomodel->numberofelements+1)){
     193                                                obs_vector_type=2;
     194                                        }
     195                                        else
     196                                         _error_("misfit observation size not supported yet");
     197
     198                                        int weight_vector_type=0;
     199                                        if ((misfit_weights_M_s[j]==iomodel->numberofvertices) || (misfit_weights_M_s[j]==iomodel->numberofvertices+1)){
     200                                                weight_vector_type=1;
     201                                        }
     202                                        else if ((misfit_weights_M_s[j]==iomodel->numberofelements) || (misfit_weights_M_s[j]==iomodel->numberofelements+1)){
     203                                                weight_vector_type=2;
     204                                        }
     205                                        else
     206                                         _error_("misfit weight size not supported yet");
     207
     208                                        /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
     209                                        output_definitions->AddObject(new MisfitAnnual(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
     210
     211                                        /*Now, for this particular misfit object, make sure we plug into the elements: the observation, and the weights.*/
     212                                        for(Object* & object : elements->objects){
     213                                                Element* element=xDynamicCast<Element*>(object);
     214                                                element->InputCreate(misfit_observation_s[j],inputs,iomodel,misfit_observation_M_s[j],misfit_observation_N_s[j],obs_vector_type,StringToEnumx(misfit_observation_string_s[j]),7);
     215                                                element->InputCreate(misfit_weights_s[j],inputs,iomodel,misfit_weights_M_s[j],misfit_weights_N_s[j],weight_vector_type,StringToEnumx(misfit_weights_string_s[j]),7);
     216                                        }
     217
     218                                }
     219
     220                                /*Free resources:*/
     221                                for(j=0;j<nummisfits;j++){
     222                                        char* string=NULL;
     223                                        IssmDouble* matrix = NULL;
     224                                        string = misfit_definitionstring_s[j];          xDelete<char>(string);
     225                                        string = misfit_observation_string_s[j];        xDelete<char>(string);
     226                                        string = misfit_model_string_s[j];                      xDelete<char>(string);
     227                                        string = misfit_weights_string_s[j];            xDelete<char>(string);
     228                                        string = misfit_name_s[j];    xDelete<char>(string);
     229                                        string = misfit_timeinterpolation_s[j];    xDelete<char>(string);
     230                                        matrix = misfit_observation_s[j]; xDelete<IssmDouble>(matrix);
     231                                        matrix = misfit_weights_s[j]; xDelete<IssmDouble>(matrix);
     232                                }
     233                                xDelete<char*>(misfit_name_s);
     234                                xDelete<char*>(misfit_model_string_s);
     235                                xDelete<char*>(misfit_definitionstring_s);
     236                                xDelete<IssmDouble*>(misfit_observation_s);
     237                                xDelete<char*>(misfit_observation_string_s);
     238                                xDelete<int>(misfit_observation_M_s);
     239                                xDelete<int>(misfit_observation_N_s);
     240                                xDelete<int>(misfit_local_s);
     241                                xDelete<char*>(misfit_timeinterpolation_s);
     242                                xDelete<IssmDouble*>(misfit_weights_s);
     243                                xDelete<int>(misfit_weights_M_s);
     244                                xDelete<int>(misfit_weights_N_s);
     245                                xDelete<char*>(misfit_weights_string_s);
     246                                /*}}}*/
     247                        }
    156248                        else if (output_definition_enums[i]==CfsurfacesquareEnum){
    157249                                /*Deal with cfsurfacesquare: {{{*/
    … …  
    205297
    206298                                        /*First create a cfsurfacesquare object for that specific string (cfsurfacesquare_model_string_s[j]):*/
    207                                         output_definitions->AddObject(new Cfsurfacesquare(cfsurfacesquare_name_s[j],StringToEnumx(cfsurfacesquare_definitionstring_s[j]),StringToEnumx(cfsurfacesquare_model_string_s[j]),cfsurfacesquare_datatime_s[j]));
     299                                        output_definitions->AddObject(new Cfsurfacesquare(cfsurfacesquare_name_s[j],StringToEnumx(cfsurfacesquare_definitionstring_s[j]),StringToEnumx(cfsurfacesquare_model_string_s[j]),StringToEnumx(cfsurfacesquare_observation_string_s[j]),StringToEnumx(cfsurfacesquare_weights_string_s[j]),cfsurfacesquare_datatime_s[j],false));
    208300
    209301                                        /*Now, for this particular cfsurfacesquare object, make sure we plug into the elements: the observation, and the weights.*/
    210302                                        for(Object* & object : elements->objects){
    211303                                                Element* element=xDynamicCast<Element*>(object);
    212                                                 element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j],inputs,iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),SurfaceObservationEnum);
    213                                                 element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j],inputs,iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),WeightsSurfaceObservationEnum);
     304                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j],inputs,iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum);
     305                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j],inputs,iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    214306
    215307                                        }
    … …  
    244336                                /*}}}*/
    245337                        }
    246                         else if (output_definition_enums[i]==CfsurfacesquaretransientEnum){
    247                                 /*Deal with cfsurfacesquaretransient: {{{*/
    248 
    249                                 /*cfsurfacesquaretransient variables: */
    250                                 int          num_cfsurfacesquaretransients,test;
    251                                 char       **cfssqt_name_s                = NULL;
    252                                 char       **cfssqt_definitionstring_s    = NULL;
    253                                 char       **cfssqt_model_string_s        = NULL;
    254                                 IssmDouble **cfssqt_observations_s        = NULL;
    255                                 int         *cfssqt_observations_M_s      = NULL;
    256                                 int         *cfssqt_observations_N_s      = NULL;
    257                                 IssmDouble **cfssqt_weights_s             = NULL;
    258                                 int         *cfssqt_weights_M_s           = NULL;
    259                                 int         *cfssqt_weights_N_s           = NULL;
    260 
    261                                 /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfsurfacesquaretransient.m): */
    262                                 iomodel->FetchMultipleData(&cfssqt_name_s,&num_cfsurfacesquaretransients,"md.cfsurfacesquaretransient.name");
    263                                 iomodel->FetchMultipleData(&cfssqt_definitionstring_s,&test,"md.cfsurfacesquaretransient.definitionstring"); _assert_(test==num_cfsurfacesquaretransients);
    264                                 iomodel->FetchMultipleData(&cfssqt_model_string_s,&test,"md.cfsurfacesquaretransient.model_string"); _assert_(test==num_cfsurfacesquaretransients);
    265                                 iomodel->FetchMultipleData(&cfssqt_observations_s,&cfssqt_observations_M_s,&cfssqt_observations_N_s,&test, "md.cfsurfacesquaretransient.observations"); _assert_(test==num_cfsurfacesquaretransients);
    266                                 iomodel->FetchMultipleData(&cfssqt_weights_s,&cfssqt_weights_M_s,&cfssqt_weights_N_s, &test,"md.cfsurfacesquaretransient.weights"); _assert_(test==num_cfsurfacesquaretransients);
    267 
    268                                 for(j=0;j<num_cfsurfacesquaretransients;j++){
    269 
    270                /*Check that we can use P1 inputs*/
    271                                         if (cfssqt_observations_M_s[j]!=(iomodel->numberofvertices+1)) _error_("observations should be a P1 time series");
    272                if (cfssqt_weights_M_s[j]!=iomodel->numberofvertices+1)        _error_("weights should be a P1 time series");
    273                                         _assert_(cfssqt_observations_N_s[j]>0);
    274 
    275                                         /*extract data times from last row of observations*/
    276                                         IssmDouble *datatimes = xNew<IssmDouble>(cfssqt_observations_N_s[j]);
    277                                         for(int k=0;k<cfssqt_observations_N_s[j];k++) datatimes[k] = (cfssqt_observations_s[j])[cfssqt_observations_N_s[j]*(cfssqt_weights_M_s[j]-1)+k];
    278 
    279                                         /*First create a cfsurfacesquaretransient object for that specific string (cfssqt_model_string_s[j]):*/
    280                                         output_definitions->AddObject(new Cfsurfacesquaretransient(cfssqt_name_s[j], StringToEnumx(cfssqt_definitionstring_s[j]), StringToEnumx(cfssqt_model_string_s[j]), cfssqt_observations_N_s[j],datatimes ));
    281                                         xDelete<IssmDouble>(datatimes);
    282 
    283                                         /*Now, for this particular cfsurfacesquaretransient object, make sure we plug into the elements: the observation, and the weights.*/
    284                                         for(Object* & object : elements->objects){
    285                                                 Element* element=xDynamicCast<Element*>(object);
    286                                                 element->DatasetInputAdd(StringToEnumx(cfssqt_definitionstring_s[j]),cfssqt_observations_s[j],inputs,iomodel,cfssqt_observations_M_s[j],cfssqt_observations_N_s[j],1,SurfaceObservationEnum,SurfaceObservationEnum);
    287                                                 element->DatasetInputAdd(StringToEnumx(cfssqt_definitionstring_s[j]),cfssqt_weights_s[j],inputs,iomodel,cfssqt_weights_M_s[j],cfssqt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum);
    288 
    289                                         }
    290                                 }
    291 
    292                                 /*Free resources:*/
    293                                 for(j=0;j<num_cfsurfacesquaretransients;j++){
    294                                         char* string=NULL;
    295                                         IssmDouble* matrix = NULL;
    296                                         string = cfssqt_definitionstring_s[j];          xDelete<char>(string);
    297                                         string = cfssqt_model_string_s[j];                      xDelete<char>(string);
    298                                         string = cfssqt_name_s[j];    xDelete<char>(string);
    299                                         matrix = cfssqt_observations_s[j]; xDelete<IssmDouble>(matrix);
    300                                         matrix = cfssqt_weights_s[j]; xDelete<IssmDouble>(matrix);
    301                                 }
    302                                 xDelete<char*>(cfssqt_name_s);
    303                                 xDelete<char*>(cfssqt_model_string_s);
    304                                 xDelete<char*>(cfssqt_definitionstring_s);
    305                                 xDelete<IssmDouble*>(cfssqt_observations_s);
    306                                 xDelete<int>(cfssqt_observations_M_s);
    307                                 xDelete<int>(cfssqt_observations_N_s);
    308                                 xDelete<IssmDouble*>(cfssqt_weights_s);
    309                                 xDelete<int>(cfssqt_weights_M_s);
    310                                 xDelete<int>(cfssqt_weights_N_s);
    311                                 /*}}}*/
    312                         }
    313338                        else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){
    314339                                /*Deal with cfdragcoeffabsgrad: {{{*/
    … …  
    342367
    343368                                        /*First create a cfdragcoeffabsgrad object for that specific string (cfdragcoeffabsgrad_model_string_s[j]):*/
    344                                         output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j])));
     369                                        output_definitions->AddObject(new Cfdragcoeffabsgrad(cfdragcoeffabsgrad_name_s[j],StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),false));
    345370
    346371                                        /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/
    … …  
    349374                                                Element* element=xDynamicCast<Element*>(object);
    350375
    351                                                 element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j],inputs,iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),WeightsSurfaceObservationEnum);
     376                                                element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j],inputs,iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    352377
    353378                                        }
    … …  
    373398                                /*}}}*/
    374399                        }
    375                         else if (output_definition_enums[i]==CfdragcoeffabsgradtransientEnum){
    376                                 /*Deal with cfdragcoeffabsgradtransient: {{{*/
    377 
    378                                 /*cfdragcoeffabsgrad variables: */
    379                                 int          num_cfdragcoeffabsgradtransients, test;
    380                                 char**       cfdraggradt_name_s                                         = NULL;   
    381                                 char**           cfdraggradt_definitionstring_s         = NULL;   
    382                                 IssmDouble** cfdraggradt_weights_s                                      = NULL;
    383                                 int*         cfdraggradt_weights_M_s                            = NULL;
    384                                 int*         cfdraggradt_weights_N_s                            = NULL;
    385 
    386                                 /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfdragcoeffabsgradtransient.m): */
    387                                 iomodel->FetchMultipleData(&cfdraggradt_name_s,&num_cfdragcoeffabsgradtransients,                                                        "md.cfdragcoeffabsgradtransient.name");
    388                                 iomodel->FetchMultipleData(&cfdraggradt_definitionstring_s,&num_cfdragcoeffabsgradtransients,                                            "md.cfdragcoeffabsgradtransient.definitionstring");
    389                                 iomodel->FetchMultipleData(&cfdraggradt_weights_s,&cfdraggradt_weights_M_s,&cfdraggradt_weights_N_s,&test,             "md.cfdragcoeffabsgradtransient.weights");
    390                                        
    391                                 for(j=0;j<num_cfdragcoeffabsgradtransients;j++){
    392                
    393                                         /*Check that we can use P1 inputs*/
    394                                         if (cfdraggradt_weights_M_s[j]!=iomodel->numberofvertices+1)  _error_("weights should be a P1 time series");
    395                                        
    396                                         /*extract data times from last row of observations*/
    397                                         IssmDouble *datatimes = xNew<IssmDouble>(cfdraggradt_weights_N_s[j]);
    398                                         for(int k=0;k<cfdraggradt_weights_N_s[j];k++) datatimes[k] = (cfdraggradt_weights_s[j])[cfdraggradt_weights_N_s[j]*(cfdraggradt_weights_M_s[j]-1)+k];
    399 
    400                                          /*First create a cfdragcoeffabsgradtransient object for that specific string:*/
    401                                         output_definitions->AddObject(new Cfdragcoeffabsgradtransient(cfdraggradt_name_s[j],StringToEnumx(cfdraggradt_definitionstring_s[j]), cfdraggradt_weights_N_s[j], datatimes));
    402 
    403                                         /*Now, for this particular cfdragcoeffabsgrad object, make sure we plug into the elements: the observation, and the weights.*/
    404                                         for(Object* & object : elements->objects){
    405 
    406                                                 Element* element=xDynamicCast<Element*>(object);
    407 
    408                                                 element->DatasetInputAdd(StringToEnumx(cfdraggradt_definitionstring_s[j]),cfdraggradt_weights_s[j],inputs,iomodel,cfdraggradt_weights_M_s[j],cfdraggradt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum);
    409 
    410                                         }
    411                                 }
    412 
    413                                 /*Free resources:*/
    414                                 for(j=0;j<num_cfdragcoeffabsgradtransients;j++){
    415                                         char* string=NULL;
    416                                         IssmDouble* matrix = NULL;
    417 
    418                                         string = cfdraggradt_definitionstring_s[j];             xDelete<char>(string);
    419                                         string = cfdraggradt_name_s[j];    xDelete<char>(string);
    420                                         matrix = cfdraggradt_weights_s[j]; xDelete<IssmDouble>(matrix);
    421                                 }
    422                                 xDelete<char*>(cfdraggradt_name_s);
    423                                 xDelete<char*>(cfdraggradt_definitionstring_s);
    424                                 xDelete<IssmDouble*>(cfdraggradt_weights_s);
    425                                 xDelete<int>(cfdraggradt_weights_M_s);
    426                                 xDelete<int>(cfdraggradt_weights_N_s);
    427                                 /*}}}*/
    428                         }
    429                         else if (output_definition_enums[i]==CfrheologybbarabsgradEnum){
    430                                 /*Deal with cfrheologybbarabsgrad: {{{*/
    431 
    432                                 /*cfrheologybbarabsgrad variables: */
    433                                 int          num_cfrheologybbarabsgrads;
    434                                 char**       cfrheologybbarabsgrad_name_s                = NULL;
    435                                 char**       cfrheologybbarabsgrad_definitionstring_s    = NULL;
    436                                 IssmDouble** cfrheologybbarabsgrad_weights_s             = NULL;
    437                                 int*         cfrheologybbarabsgrad_weights_M_s           = NULL;
    438                                 int*         cfrheologybbarabsgrad_weights_N_s           = NULL;
    439                                 char**       cfrheologybbarabsgrad_weights_string_s      = NULL;
    440 
    441                                 /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfrheologybbarabsgrad.m): */
    442                                 iomodel->FetchMultipleData(&cfrheologybbarabsgrad_name_s,&num_cfrheologybbarabsgrads,                                                        "md.cfrheologybbarabsgrad.name");
    443                                 iomodel->FetchMultipleData(&cfrheologybbarabsgrad_definitionstring_s,&num_cfrheologybbarabsgrads,                                            "md.cfrheologybbarabsgrad.definitionstring");
    444                                 iomodel->FetchMultipleData(&cfrheologybbarabsgrad_weights_s,&cfrheologybbarabsgrad_weights_M_s,&cfrheologybbarabsgrad_weights_N_s,&num_cfrheologybbarabsgrads,             "md.cfrheologybbarabsgrad.weights");
    445                                 iomodel->FetchMultipleData(&cfrheologybbarabsgrad_weights_string_s,&num_cfrheologybbarabsgrads,                                              "md.cfrheologybbarabsgrad.weights_string");
    446 
    447                                 for(j=0;j<num_cfrheologybbarabsgrads;j++){
    448 
    449                                         int weight_vector_type=0;
    450                                         if ((cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofvertices) || (cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofvertices+1)){
    451                                                 weight_vector_type=1;
    452                                         }
    453                                         else if ((cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofelements) || (cfrheologybbarabsgrad_weights_M_s[j]==iomodel->numberofelements+1)){
    454                                                 weight_vector_type=2;
    455                                         }
    456                                         else
    457                                          _error_("cfrheologybbarabsgrad weight size not supported yet");
    458 
    459                                         /*First create a cfrheologybbarabsgrad object for that specific string (cfrheologybbarabsgrad_model_string_s[j]):*/
    460                                         output_definitions->AddObject(new Cfrheologybbarabsgrad(cfrheologybbarabsgrad_name_s[j],StringToEnumx(cfrheologybbarabsgrad_definitionstring_s[j]),StringToEnumx(cfrheologybbarabsgrad_weights_string_s[j])));
    461 
    462                                         /*Now, for this particular cfrheologybbarabsgrad object, make sure we plug into the elements: the observation, and the weights.*/
    463                                         for(Object* & object : elements->objects){
    464 
    465                                                 Element* element=xDynamicCast<Element*>(object);
    466 
    467                                                 element->DatasetInputAdd(StringToEnumx(cfrheologybbarabsgrad_definitionstring_s[j]),cfrheologybbarabsgrad_weights_s[j],inputs,iomodel,cfrheologybbarabsgrad_weights_M_s[j],cfrheologybbarabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfrheologybbarabsgrad_weights_string_s[j]),WeightsSurfaceObservationEnum);
    468 
    469                                         }
    470 
    471                                 }
    472                                     /*Free resources:*/
    473             for(j=0;j<num_cfrheologybbarabsgrads;j++){
    474                char* string=NULL;
    475                IssmDouble* matrix = NULL;
    476 
    477                string = cfrheologybbarabsgrad_definitionstring_s[j];    xDelete<char>(string);
    478                string = cfrheologybbarabsgrad_weights_string_s[j];      xDelete<char>(string);
    479                string = cfrheologybbarabsgrad_name_s[j];    xDelete<char>(string);
    480                matrix = cfrheologybbarabsgrad_weights_s[j]; xDelete<IssmDouble>(matrix);
    481             }
    482             xDelete<char*>(cfrheologybbarabsgrad_name_s);
    483             xDelete<char*>(cfrheologybbarabsgrad_definitionstring_s);
    484             xDelete<IssmDouble*>(cfrheologybbarabsgrad_weights_s);
    485             xDelete<int>(cfrheologybbarabsgrad_weights_M_s);
    486             xDelete<int>(cfrheologybbarabsgrad_weights_N_s);
    487             xDelete<char*>(cfrheologybbarabsgrad_weights_string_s);
    488             /*}}}*/
    489          }
    490                         else if (output_definition_enums[i]==CfrheologybbarabsgradtransientEnum){
    491                                 /*Deal with cfrheologybbarabsgradtransient: {{{*/
    492 
    493                                 /*cfrheologybbarabsgrad variables: */
    494                                 int          num_cfrheologybbarabsgradtransients, test;
    495                                 char**       cfrheogradt_name_s                = NULL;
    496                                 char**       cfrheogradt_definitionstring_s    = NULL;
    497                                 IssmDouble** cfrheogradt_weights_s             = NULL;
    498                                 int*         cfrheogradt_weights_M_s           = NULL;
    499                                 int*         cfrheogradt_weights_N_s           = NULL;
    500                                 char**       cfrheogradt_weights_string_s      = NULL;
    501 
    502                                 /*Fetch name, model_string, observation, observation_string, etc ... (see src/m/classes/cfrheologybbarabsgradtransient.m): */
    503                                 iomodel->FetchMultipleData(&cfrheogradt_name_s,&num_cfrheologybbarabsgradtransients,                                                        "md.cfrheologybbarabsgradtransient.name");
    504                                 iomodel->FetchMultipleData(&cfrheogradt_definitionstring_s,&num_cfrheologybbarabsgradtransients,                                            "md.cfrheologybbarabsgradtransient.definitionstring");
    505                                 iomodel->FetchMultipleData(&cfrheogradt_weights_s,&cfrheogradt_weights_M_s,&cfrheogradt_weights_N_s,&test,             "md.cfrheologybbarabsgradtransient.weights");
    506 
    507                                 for(j=0;j<num_cfrheologybbarabsgradtransients;j++){
    508 
    509                                         if (cfrheogradt_weights_M_s[j]!=iomodel->numberofvertices+1) _error_("weights should be a P1 time series");
    510                                        
    511                                         /*extract data times from last row of observations*/
    512                                         IssmDouble *datatimes = xNew<IssmDouble>(cfrheogradt_weights_N_s[j]);
    513                                         for(int k=0;k<cfrheogradt_weights_N_s[j];k++) datatimes[k] = (cfrheogradt_weights_s[j])[cfrheogradt_weights_N_s[j]*(cfrheogradt_weights_M_s[j]-1)+k];
    514 
    515                                         /*First create a cfrheologybbarabsgradtransient object for that specific string:*/
    516                                         output_definitions->AddObject(new Cfrheologybbarabsgradtransient(cfrheogradt_name_s[j],StringToEnumx(cfrheogradt_definitionstring_s[j]), cfrheogradt_weights_N_s[j], datatimes));
    517 
    518                                         /*Now, for this particular cfrheologybbarabsgrad object, make sure we plug into the elements: the observation, and the weights.*/
    519                                         for(Object* & object : elements->objects){
    520 
    521                                                 Element* element=xDynamicCast<Element*>(object);
    522 
    523                                                 element->DatasetInputAdd(StringToEnumx(cfrheogradt_definitionstring_s[j]),cfrheogradt_weights_s[j],inputs,iomodel,cfrheogradt_weights_M_s[j],cfrheogradt_weights_N_s[j],1,WeightsSurfaceObservationEnum,WeightsSurfaceObservationEnum);
    524 
    525                                         }
    526                                 }
    527                                
    528                                 /*Free resources:*/
    529             for(j=0;j<num_cfrheologybbarabsgradtransients;j++){
    530                char* string=NULL;
    531                IssmDouble* matrix = NULL;
    532 
    533                string = cfrheogradt_definitionstring_s[j];    xDelete<char>(string);
    534                string = cfrheogradt_name_s[j];    xDelete<char>(string);
    535                matrix = cfrheogradt_weights_s[j]; xDelete<IssmDouble>(matrix);
    536             }
    537             xDelete<char*>(cfrheogradt_name_s);
    538             xDelete<char*>(cfrheogradt_definitionstring_s);
    539             xDelete<IssmDouble*>(cfrheogradt_weights_s);
    540             xDelete<int>(cfrheogradt_weights_M_s);
    541             xDelete<int>(cfrheogradt_weights_N_s);
    542             /*}}}*/
    543          }
    544400                        else if (output_definition_enums[i]==CfsurfacelogvelEnum){
    545401                                /*Deal with cfsurfacelogvel: {{{*/
    … …  
    595451
    596452                                        /*First create a cfsurfacelogvel object for that specific string (cfsurfacelogvel_modeltring[j]):*/
    597                                         output_definitions->AddObject(new Cfsurfacelogvel(cfsurfacelogvel_name[j],StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_datatime[j]));
     453                                        output_definitions->AddObject(new Cfsurfacelogvel(cfsurfacelogvel_name[j],StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_datatime[j],false));
    598454
    599455                                        /*Now, for this particular cfsurfacelogvel object, make sure we plug into the elements: the observation, and the weights.*/
    … …  
    602458                                                Element* element=xDynamicCast<Element*>(object);
    603459
    604                                                 element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),VxObsEnum);
    605                                                         element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),VyObsEnum);
    606                                                 element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j],inputs,iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),WeightsSurfaceObservationEnum);
     460                                                element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum);
     461                                                        element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j],inputs,iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum);
     462                                                element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j],inputs,iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum);
    607463
    608464                                        }
    … …  
    689545
    690546                                        /*First create a cflevelsetmisfit object for that specific string (cflevelsetmisfit_model_string_s[j]):*/
    691                                         output_definitions->AddObject(new Cflevelsetmisfit(cflevelsetmisfit_name_s[j],StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),StringToEnumx(cflevelsetmisfit_model_string_s[j]),cflevelsetmisfit_datatime_s[j]));
     547                                        output_definitions->AddObject(new Cflevelsetmisfit(cflevelsetmisfit_name_s[j],StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),StringToEnumx(cflevelsetmisfit_model_string_s[j]),StringToEnumx(cflevelsetmisfit_observation_string_s[j]),StringToEnumx(cflevelsetmisfit_weights_string_s[j]),cflevelsetmisfit_datatime_s[j],false));
    692548
    693549                                        /*Now, for this particular cflevelsetmisfit object, make sure we plug into the elements: the observation, and the weights.*/
    694550                                        for(Object* & object : elements->objects){
    695551                                                Element* element=xDynamicCast<Element*>(object);
    696                                                 element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_observation_s[j],inputs,iomodel,cflevelsetmisfit_observation_M_s[j],cflevelsetmisfit_observation_N_s[j],obs_vector_type,StringToEnumx(cflevelsetmisfit_observation_string_s[j]),LevelsetObservationEnum);
    697                                                 element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_weights_s[j],inputs,iomodel,cflevelsetmisfit_weights_M_s[j],cflevelsetmisfit_weights_N_s[j],weight_vector_type,StringToEnumx(cflevelsetmisfit_weights_string_s[j]),WeightsLevelsetObservationEnum);
     552                                                element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_observation_s[j],inputs,iomodel,cflevelsetmisfit_observation_M_s[j],cflevelsetmisfit_observation_N_s[j],obs_vector_type,StringToEnumx(cflevelsetmisfit_observation_string_s[j]),7,LevelsetObservationEnum);
     553                                                element->DatasetInputAdd(StringToEnumx(cflevelsetmisfit_definitionstring_s[j]),cflevelsetmisfit_weights_s[j],inputs,iomodel,cflevelsetmisfit_weights_M_s[j],cflevelsetmisfit_weights_N_s[j],weight_vector_type,StringToEnumx(cflevelsetmisfit_weights_string_s[j]),7,WeightsLevelsetObservationEnum);
    698554                                        }
    699555                                }

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/m/boundaryconditions/SetMarineIceSheetBC.m

    r26352 r27956  
    3333else
    3434        %Guess where the ice front is
    35         pos=find(sum(md.mask.ocean_levelset(md.mesh.elements)<0.,2) >0.);
    36         vertexonfloatingice=zeros(md.mesh.numberofvertices,1);
    37         vertexonfloatingice(md.mesh.elements(pos,:))=1.;
    38         vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice);
     35        %pos=find(sum(md.mask.ocean_levelset(md.mesh.elements)<0.,2) >0.);
     36        %vertexonfloatingice=zeros(md.mesh.numberofvertices,1);
     37        %vertexonfloatingice(md.mesh.elements(pos,:))=1.;
     38        %vertexonicefront=double(md.mesh.vertexonboundary & vertexonfloatingice);
     39
     40    only_ocean_levelset=-double(md.mask.ocean_levelset<=0. & md.mask.ice_levelset>=0.);
     41    only_ocean_levelset_sum=sum(only_ocean_levelset(md.mesh.elements)<0.,2);
     42    pos=find(only_ocean_levelset_sum >0. & only_ocean_levelset_sum < 3.);
     43    vertexonicefront=zeros(md.mesh.numberofvertices,1);
     44    vertexonicefront(md.mesh.elements(pos,:))=1.;
    3945end
    4046pos=find(md.mesh.vertexonboundary & ~vertexonicefront);
    … …  
    5965        error('mesh type not supported yet');
    6066end
    61 segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0;
    62 segments=find(sum(segmentsfront,2)~=numbernodesfront);
     67%segmentsfront=md.mask.ice_levelset(md.mesh.segments(:,1:numbernodesfront))==0;
     68%segments=find(sum(segmentsfront,2)~=numbernodesfront);
    6369%Find all nodes for these segments and spc them
    64 pos=md.mesh.segments(segments,1:numbernodesfront);
     70%pos=md.mesh.segments(segments,1:numbernodesfront);
     71
     72numbernodesfront=3;
     73elementsfront=md.mask.ice_levelset(md.mesh.elements(:,1:numbernodesfront))==0;
     74elements=find(sum(elementsfront,2)>0);
     75pos=md.mesh.elements(elements,1:numbernodesfront);
     76
    6577md.stressbalance.spcvx(pos(:))=0;
    6678md.stressbalance.spcvy(pos(:))=0;

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/m/exp/exptool.m

    r27950 r27956  
    1616%      - markersize (default=7)
    1717%      - markeredgecolor (default='r')
    18 %      - nofigurecopy (default=0) do not copy current figure, this is needed on some platform to avoid an offset in the figure
     18%      - nofigurecopy (default=0) do not copy current figure, this is needed on some plateform to avoid an offset in the figure
    1919%
    2020%   Usage:
    … …  
    139139disableDefaultInteractivity(gca); %disables the built-in interactions for the specified axes
    140140
    141 %Build backup structure for do and redo
     141%Build backup structre for do and redo
    142142backup=cell(1,3);
    143143backup{1,1}=A;

  • TabularUnified issm/branches/trunk-dlcheng-ASE/src/m/solve/waitonlock.m

    r27683 r27956  
    5858                command = [command ' "[ -f ' lockfilename ' ] && [ -f ' logfilename ' ]" 2>/dev/null'];
    5959        end
     60    disp(command)
     61    disp(num2str(cluster.port))
    6062end
    6163
Note: See TracChangeset for help on using the changeset viewer.