Changeset 27491


Ignore:
Timestamp:
01/02/23 00:15:34 (2 years ago)
Author:
rueckamp
Message:

CHG: MinThickness Mask and minor changes for Debris Analysis

Location:
issm/trunk-jpl/src/c
Files:
13 edited

Legend:

Unmodified
Added
Removed

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

    r27304 r27491  
    1010
    1111#define FINITEELEMENT P1Enum
     12#define EPS 1e-14
    1213
    1314/*Model processing*/
    … …  
    5960
    6061        if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    61         ::CreateNodes(nodes,iomodel,DebrisAnalysisEnum,FINITEELEMENT);
     62        ::CreateNodes(nodes,iomodel,DebrisAnalysisEnum,FINITEELEMENT,isamr);
    6263        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    6364
    … …  
    127128void           DebrisAnalysis::Core(FemModel* femmodel){/*{{{*/
    128129
    129 
    130         Element* element= NULL;
    131         for(Object* & object : femmodel->elements->objects){
    132                 element=xDynamicCast<Element*>(object);
    133 
    134 
    135                 int numvertices = element->GetNumberOfNodes();
    136 
    137                 IssmDouble* vx = xNew<IssmDouble>(numvertices);
    138                 IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
    139                 IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
    140                 IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
    141                 IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
    142 
    143                 element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
    144                 element->GetInputListOnVertices(&vx[0],VxEnum);
    145                 element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
    146                 element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
    147                 element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    148 
    149                 IssmDouble slope,rad2deg=180./M_PI; //=57.2958
    150                 IssmDouble vslipx,rhod=1900.;
    151                 IssmDouble gravity=element->FindParam(ConstantsGEnum);
    152                 IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
    153                 IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
    154 
    155                 int step;
    156                 IssmDouble dt, maxv;
    157                 IssmDouble yts=31536000.;
    158                 femmodel->parameters->FindParam(&step,StepEnum);
    159                 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    160 
    161                 bool isminthicknessinelement=false;
    162                 for(int i=0;i<numvertices;i++){
    163                         if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
    164                 }
    165                 if(isminthicknessinelement){
    166                         //do nothing
    167                         for(int i=0;i<numvertices;i++){
    168                                 if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;                         
    169                         }
    170                 }else{
    171                         for(int i=0;i<numvertices;i++){
    172                                 //if(onsurface[i]>.5){
    173                                 slope=fabs(slopex[i]);
    174                                 if((atan(slope)*rad2deg)>25.){
    175                                         //if(debristhickness[i]>0.01){
    176                                         vslipx=slope_threshold/yts;
    177                                         //maxv=10.0/2./dt;
    178                                         //vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts;
    179                                         vx[i]=vx[i]+vslipx;
    180                                         //debristhickness[i]=debristhickness[i];
    181                                         //if(vx[i]>maxv) vx[i]=maxv;
    182                                         //}
    183                                 }
    184                                 //}
    185                         }
    186                 }
    187                 //if(step%100==0)   
    188                 element->AddInput(VxDebrisEnum,vx,P1Enum);
    189                 //element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    190 
    191                 /* Free resources */
    192                 xDelete<IssmDouble>(debristhickness);
    193                 xDelete<IssmDouble>(icethickness);
    194                 xDelete<IssmDouble>(vx);
    195                 xDelete<IssmDouble>(slopex);
    196                 xDelete<IssmDouble>(onsurface);
    197         }
    198 
    199         //if(step%7==0)
    200130        //PreProcessing(femmodel);
    201131        //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
    202132        //extrudefromtop_core(femmodel);
    203 
    204133        femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);       
    205134        solutionsequence_linear(femmodel);
    206 
    207         int step;
    208         femmodel->parameters->FindParam(&step,StepEnum);
    209         //if(step%7==0) PreProcessing(femmodel);
    210135        PostProcessing(femmodel);
    211136
    … …  
    230155        IssmDouble Jdet,D_scalar,dt,h;
    231156        IssmDouble vel,vx,vy,dvxdx,dvydy;
    232         IssmDouble xi,tau;
     157        IssmDouble yts=31536000.;
     158        IssmDouble tau;
    233159        IssmDouble dvx[2],dvy[2];
    234160        Element*    topelement = NULL;
    … …  
    271197        Input* vy_input=NULL;
    272198        if(dim>1){vy_input = topelement->GetInput(VyDebrisEnum); _assert_(vy_input);}
    273         h = topelement->CharacteristicLength();
     199        h=topelement->CharacteristicLength();
    274200
    275201        /* Start  looping on the number of gaussian points: */
    … …  
    302228                                }
    303229                        }
    304                 }
    305                 else{
     230                }else{
    306231                        dvxdx=dvx[0];
    307232                        for(int i=0;i<numnodes;i++){
    … …  
    313238                }
    314239
    315                 /*******************************************************************/
    316                 /* Diffusion */
    317                 bool isdisplacement=false;
    318                 int step;
    319                 topelement->FindParam(&step,StepEnum);
    320                 IssmDouble slope_threshold;
    321                 topelement->FindParam(&slope_threshold,DebrisRemovalSlopeThresholdEnum);
    322                 IssmDouble kappa,f,smb,debristhickness,slopex;
    323                 IssmDouble Diff,fraction,M=1,C;
    324                 IssmDouble rad2deg=180./M_PI;
    325                 Diff=3.2/3e7;
    326                 Input* slopex_input=topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    327                 Input* smb_input=topelement->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    328                 Input* debristhickness_input=topelement->GetInput(DebrisThicknessEnum); _assert_(debristhickness_input);
    329 
    330                 if(isdisplacement){
    331 
    332                         slopex_input->GetInputValue(&slopex, gauss);
    333                         smb_input->GetInputValue(&smb, gauss);
    334                         debristhickness_input->GetInputValue(&debristhickness, gauss);
    335                         if((atan(fabs(slopex))*rad2deg)>30.){
    336                                 f=1.;
    337                         }else{
    338                                 f=0.;
    339                         }
    340                         //f=1;
    341                         //kappa=-5.6e16*smb*debristhickness*f;
    342                         //kappa=debristhickness/h*4e9*f;
    343                         //kappa=14.2809e8*f; // 25°
    344                         kappa=slope_threshold*1e8*f;
    345                         if(dim==2){
    346                                 for(int i=0;i<numnodes;i++){
    347                                         for(int j=0;j<numnodes;j++){
    348                                                 Ke->values[i*numnodes+j] +=  D_scalar*kappa*(
    349                                                                 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
    350                                                                 );
    351                                         }
    352                                 }
    353                         }else{
    354                                 for(int i=0;i<numnodes;i++){
    355                                         for(int j=0;j<numnodes;j++){
    356                                                 Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i]);
    357                                         }
    358                                 }
    359                         }               
    360                 }
    361 
    362                 /*******************************************************************/               
    363 
    364240                IssmDouble rho;
    365241                if(FINITEELEMENT==P1Enum){
    366                         rho=2;
     242                        rho=2.;
    367243                }else if(FINITEELEMENT==P2Enum){
    368244                        rho=4.;
    369245                }
    370                 if(stabilization==2){
     246
     247                for(int i=0;i<(dim*dim);i++) D[i]=0.;
     248                if(stabilization==1){
     249                        /*SSA*/
     250                        if(dim==1){
     251                                vx_input->GetInputValue(&vx,gauss);
     252                                D[0]=h/rho*fabs(vx);
     253                        }else{
     254                                vx_input->GetInputValue(&vx,gauss);
     255                                vy_input->GetInputValue(&vy,gauss);
     256                                vel=sqrt(vx*vx+vy*vy);
     257                                D[0*dim+0]=h/rho*fabs(vx);
     258                                D[1*dim+1]=h/rho*fabs(vy);
     259                        }
     260                }else if(stabilization==2){ 
    371261                        /*Streamline upwinding*/
    372262                        if(dim==1){
    373263                                vx_input->GetInputValue(&vx,gauss);
    374                                 vel=fabs(vx)+1.e-10;
     264                                vel=fabs(vx)+EPS;
    375265                                D[0] = h/(rho*vel)*vx*vx;
    376                         }
    377                         else{
    378                                 vx_input->GetInputAverage(&vx);
    379                                 vy_input->GetInputAverage(&vy);
    380                                 vel=sqrt(vx*vx+vy*vy)+1.e-10;
     266                        }else{
     267                                vx_input->GetInputValue(&vx,gauss);
     268                                vy_input->GetInputValue(&vy,gauss);
     269                                vel=sqrt(vx*vx+vy*vy)+EPS;
    381270                                D[0*dim+0]=h/(rho*vel)*vx*vx;
    382271                                D[1*dim+0]=h/(rho*vel)*vy*vx;
    383272                                D[0*dim+1]=h/(rho*vel)*vx*vy;
    384273                                D[1*dim+1]=h/(rho*vel)*vy*vy;
    385                         }
    386                 }
    387                 else if(stabilization==1){ 
    388                         /*SSA*/
    389                         if(dim==1){
    390                                 vx_input->GetInputAverage(&vx);
    391                                 D[0]=h/rho*fabs(vx);
    392                         }
    393                         else{
    394                                 vx_input->GetInputAverage(&vx);
    395                                 vy_input->GetInputAverage(&vy);
    396                                 D[0*dim+0]=h/rho*fabs(vx);
    397                                 D[1*dim+1]=h/rho*fabs(vy);
    398                         }
    399                 }
    400                 else if(stabilization==3){ 
     274                        }               
     275                }else if(stabilization==3){ 
    401276                        /*SUPG*/
    402277                        if(dim==1){
    403                                 vx_input->GetInputAverage(&vx);
    404                                 tau=h/(rho*fabs(vx)+1e-10);
    405                         }
    406                         else{
    407                                 vx_input->GetInputAverage(&vx);
    408                                 vy_input->GetInputAverage(&vy);
    409                                 tau=1*h/(rho*sqrt(vx*vx+vy*vy)+1e-10);
    410                         }
    411                 }
     278                                vx_input->GetInputValue(&vx,gauss);
     279                                tau=h/(rho*max(fabs(vx),EPS));
     280                        }else{
     281                                vx_input->GetInputValue(&vx,gauss);
     282                                vy_input->GetInputValue(&vy,gauss);
     283                                tau=h/(rho*sqrt(vx*vx+vy*vy)+EPS);
     284                        }
     285                }
     286
    412287                if(stabilization==1 || stabilization==2){
    413288                        for(int i=0;i<dim*dim;i++) D[i]=D_scalar*D[i];
    … …  
    417292                                                Ke->values[i*numnodes+j] += (
    418293                                                                dbasis[0*numnodes+i] *(D[0*dim+0]*dbasis[0*numnodes+j] + D[0*dim+1]*dbasis[1*numnodes+j]) +
    419                                                                 dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j])
    420                                                                 );
     294                                                                dbasis[1*numnodes+i] *(D[1*dim+0]*dbasis[0*numnodes+j] + D[1*dim+1]*dbasis[1*numnodes+j]));
    421295                                        }   
    422296                                }
    423                         }
    424                         else{
     297                        }else{
    425298                                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += dbasis[0*numnodes+i]*D[0]*dbasis[0*numnodes+j];
    426299                        }
    427300                }else if(stabilization==3){
    428                         /*Mass matrix - part 2*/
    429301                        if(dim==1){
     302                                /*Mass matrix - part 2*/
    430303                                for(int i=0;i<numnodes;i++){
    431304                                        for(int j=0;j<numnodes;j++){
    … …  
    454327                                }
    455328
     329
    456330                                /*Advection matrix - part 2, B*/
    457331                                for(int i=0;i<numnodes;i++){
    … …  
    464338                                for(int i=0;i<numnodes;i++){
    465339                                        for(int j=0;j<numnodes;j++){
    466                                                 Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx);
     340                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);;
     341                                        }
     342                                }
     343                        }else if(dim==2){
     344                                /*Mass matrix - part 2*/
     345                                for(int i=0;i<numnodes;i++){
     346                                        for(int j=0;j<numnodes;j++){
     347                                                Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     348                                        }
     349                                }
     350                                /*Mass matrix - part 3*/
     351                                for(int i=0;i<numnodes;i++){
     352                                        for(int j=0;j<numnodes;j++){
     353                                                Ke->values[i*numnodes+j]+=gauss->weight*Jdet*tau*basis[j]*(basis[i]*dvxdx+basis[i]*dvydy);
     354                                        }
     355                                }
     356
     357                                /*Advection matrix - part 2, A*/
     358                                for(int i=0;i<numnodes;i++){
     359                                        for(int j=0;j<numnodes;j++){
     360                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     361                                        }
     362                                }
     363                                /*Advection matrix - part 3, A*/
     364                                for(int i=0;i<numnodes;i++){
     365                                        for(int j=0;j<numnodes;j++){
     366                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(vx*dbasis[0*numnodes+j]+vy*dbasis[1*numnodes+j])*(basis[i]*dvxdx+basis[i]*dvydy);
     367                                        }
     368                                }
     369
     370                                /*Advection matrix - part 2, B*/
     371                                for(int i=0;i<numnodes;i++){
     372                                        for(int j=0;j<numnodes;j++){
     373                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     374                                        }
     375                                }
     376                                /*Advection matrix - part 3, B*/
     377                                for(int i=0;i<numnodes;i++){
     378                                        for(int j=0;j<numnodes;j++){
     379                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx+basis[j]*dvydy)*(basis[i]*dvxdx+basis[i]*dvydy);
    467380                                        }
    468381                                }
    … …  
    487400        int     stabilization,dim,domaintype;
    488401        IssmDouble  Jdet,dt;
    489         IssmDouble  smb,thickness;
    490         IssmDouble  vx,vy,vel,dvxdx,dvydy,xi,h,tau,pf;
     402        IssmDouble  smb,thickness,psi;
     403        IssmDouble  vx,vy,vel,dvxdx,dvydy,h,tau,pf;
     404        IssmDouble yts=31536000.;
    491405        IssmDouble  dvx[2],dvy[2];
    492406        IssmDouble* xyz_list = NULL;
    … …  
    533447                vy_input=topelement->GetInput(VyDebrisEnum); _assert_(vy_input);
    534448        }
     449        h=topelement->CharacteristicLength();
    535450
    536451        IssmDouble rho;
    … …  
    550465                smb_input->GetInputValue(&smb,gauss);
    551466                thickness_input->GetInputValue(&thickness,gauss);
    552 
    553467                if(smb>0.){
    554                         for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness-0.*dt*smb*pf)*basis[i];
    555                 } else {
    556                         for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness-dt*smb*pf)*basis[i]; // take the negative of melt, because it is a debris production term here
    557                 }
     468                        psi=thickness-0.*dt*smb*pf;
     469                }else{
     470                        psi=thickness-dt*smb*pf;
     471                }
     472
     473                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*psi*basis[i];
    558474
    559475                if(stabilization==3){
    560476                        /*SUPG*/
    561477                        topelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     478                        vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     479                        dvxdx=dvx[0];
    562480                        if(dim==1){
    563 
    564481                                vx_input->GetInputValue(&vx,gauss);
    565                                 vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    566                                 dvxdx=dvx[0];
    567                                 tau=h/(rho*fabs(vx)+1e-10);
    568                                 IssmDouble psi;
    569                                 if(smb>0.){
    570                                         psi=thickness;
    571                                 } else {
    572                                         psi=thickness-dt*smb*pf;
    573                                 }
     482                                tau=h/(rho*max(fabs(vx),EPS));
     483
    574484                                /*Force vector - part 2*/
    575485                                for(int i=0;i<numnodes;i++){
    576                                         pe->values[i]+=Jdet*gauss->weight*psi*tau*(vx*dbasis[0*numnodes+i]);
     486                                        pe->values[i]+=Jdet*gauss->weight*psi*tau*vx*dbasis[0*numnodes+i];
    577487                                }
    578488                                /*Force vector - part 3*/
    579489                                for(int i=0;i<numnodes;i++){
    580                                         pe->values[i]+=Jdet*gauss->weight*psi*tau*(basis[i]*dvxdx);
    581                                 }
    582 
     490                                        pe->values[i]+=Jdet*gauss->weight*psi*tau*basis[i]*dvxdx;
     491                                }
     492
     493                        }else if(dim==2){
     494                                vx_input->GetInputValue(&vx,gauss);
     495                                vy_input->GetInputValue(&vy,gauss);
     496                                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     497                                vel=sqrt(vx*vx+vy*vy);
     498                                dvydy=dvy[1];
     499                                tau=h/(rho*vel+EPS);
     500
     501                                /*Force vector - part 2*/
     502                                for(int i=0;i<numnodes;i++){
     503                                        pe->values[i]+=Jdet*gauss->weight*psi*tau*(vx*dbasis[0*numnodes+i]+vy*dbasis[1*numnodes+i]);
     504                                }
     505                                /*Force vector - part 3*/
     506                                for(int i=0;i<numnodes;i++){
     507                                        pe->values[i]+=Jdet*gauss->weight*psi*tau*(basis[i]*dvxdx+basis[i]*dvydy);
     508                                }
    583509                        }
    584510                }
    … …  
    600526}/*}}}*/
    601527void           DebrisAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    602         //element->InputUpdateFromSolutionOneDof(solution,DebrisThicknessEnum);
    603         int*         ddoflist=NULL;
     528
     529        int *ddoflist=NULL;
    604530
    605531        int numnodes = element->GetNumberOfNodes();
    606         IssmDouble* thickness     = xNew<IssmDouble>(numnodes);
    607         IssmDouble* thicknessold  = xNew<IssmDouble>(numnodes);
    608532        IssmDouble* newthickness  = xNew<IssmDouble>(numnodes);
    609         IssmDouble* icethickness  = xNew<IssmDouble>(numnodes);
    610         IssmDouble* bedslopex     = xNew<IssmDouble>(numnodes);
    611         IssmDouble* surfaceslopex = xNew<IssmDouble>(numnodes);
    612533
    613534        /*Use the dof list to index into the solution vector: */
    614535        IssmDouble minthickness = element->FindParam(DebrisMinThicknessEnum);
    615         IssmDouble iceminthickness = element->FindParam(MasstransportMinThicknessEnum);
    616         element->GetInputListOnVertices(&thickness[0],DebrisThicknessEnum);   
    617         element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    618         element->GetInputListOnVertices(&bedslopex[0],BedSlopeXEnum);
    619         element->GetInputListOnVertices(&surfaceslopex[0],SurfaceSlopeXEnum);
    620536        element->GetDofListLocal(&ddoflist,NoneApproximationEnum,GsetEnum);
    621537
    … …  
    625541                if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
    626542
    627                 /* check for thickness<minthickness */
    628                 if(thickness[i]<minthickness) newthickness[i]=minthickness;
    629 
    630                 /* Carlos model sets all values below Hmin to zero */
    631                 if(icethickness[i]<=(iceminthickness+0.0001) & fabs((fabs(surfaceslopex[i])-fabs(bedslopex[i])))<1e-3 ) newthickness[i]=0;
    632                 //if(icethickness[i]<=(iceminthickness+0.01)) newthickness[i]=0;
    633         }
    634 
    635         /* update inputs */
     543                // check for thickness<minthickness
     544                if(newthickness[i]<minthickness) newthickness[i]=minthickness;
     545        }
     546
     547        // update inputs
    636548        element->AddInput(DebrisThicknessEnum,newthickness,P1Enum);
    637549
    638         /* Free resources */
     550        // Free resources
    639551        xDelete<IssmDouble>(newthickness);
    640         xDelete<IssmDouble>(thickness);
    641         xDelete<IssmDouble>(icethickness);
    642         xDelete<IssmDouble>(bedslopex);
    643         xDelete<IssmDouble>(surfaceslopex);
    644552        xDelete<int>(ddoflist);
     553        //*/
    645554}/*}}}*/
    646555void           DebrisAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    647         //        return;
    648         SetActiveNodesLSMx(femmodel);
     556        //SetActiveNodesLSMx(femmodel);
     557
     558        // Update active elements based on ice levelset and ocean levelset*/
     559        GetMaskOfIceVerticesLSMx(femmodel,false,true); //FIXME?
     560        SetActiveNodesLSMx(femmodel,false,true); //FIXME?
     561
     562        /*Constrain all nodes that are grounded and unconstrain the ones that float*/
     563        for(Object* & object : femmodel->elements->objects){
     564                Element    *element  = xDynamicCast<Element*>(object);
     565                int         numnodes  = element->GetNumberOfNodes();
     566                IssmDouble *mask      = xNew<IssmDouble>(numnodes);
     567                IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
     568
     569                element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
     570                element->GetInputListOnNodes(&ls_active[0],DebrisMaskNodeActivationEnum);
     571
     572                for(int in=0;in<numnodes;in++){
     573                        Node* node=element->GetNode(in);
     574                        if(mask[in]>0. && ls_active[in]==1.){
     575                                // Do nothing
     576                                node->Activate(); //Not sure if we need this!
     577                        }
     578                        else{
     579                                IssmDouble phi=0;
     580                                node->Deactivate();// Not sure if we need this
     581                                node->ApplyConstraint(0,phi);
     582                        }
     583                }
     584                xDelete<IssmDouble>(mask);
     585                xDelete<IssmDouble>(ls_active);
     586        }
     587        //*/
     588        return;
    649589}/*}}}*/
    650590void           DebrisAnalysis::PostProcessing(FemModel* femmodel){/*{{{*/
    … …  
    676616                        IssmDouble* slopey         = xNew<IssmDouble>(numnodes);
    677617                        IssmDouble* onsurface      = xNew<IssmDouble>(numnodes);
     618                        IssmDouble* ls_active      = xNew<IssmDouble>(numnodes);
    678619                        element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
    679620                        element->GetInputListOnNodes(icethickness,ThicknessEnum);
    680621                        element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
     622                        element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
     623
    681624                        dim=1;
    682625                        element->GetInputListOnNodes(slopex,SurfaceSlopeXEnum);
    … …  
    688631                        bool isminthicknessinelement=false;
    689632                        bool remove_debris=false;
     633                        bool isactive=false;
    690634
    691635                        IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);                       
    … …  
    693637                        switch(removalmodel){
    694638                                case 1:{
    695                                         IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
    696 
    697                                         for(k=0; k<numnodes;k++){
    698                                                 if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
    699                                         }
    700                                         isminthicknessinelement=true;
    701                                         if(isminthicknessinelement){
    702                                                 for(k=0; k<numnodes;k++){
    703                                                         if(onsurface[k]>0.5){
    704                                                                 slope=fabs(slopex[k]);
    705                                                                 if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
    706                                                                 if((atan(slope)*rad2deg)>slope_threshold) debristhickness[k]=remove_debris=true;
    707                                                         }
    708                                                 }
    709                                                 if(remove_debris){
    710                                                         for(k=0; k<numnodes;k++){
    711                                                                 if(icethickness[k]<=(iceminthickness+0.01)) debristhickness[k]=0.;
    712                                                         }
    713                                                 }
    714                                         }
    715                                         //int finite_element = element->GetElementType();
    716                                         //element->AddInput(DebrisThicknessEnum,debristhickness,FINITEELEMENT);
    717                                         element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    718 
    719                                         xDelete<IssmDouble>(debristhickness);
    720                                         xDelete<IssmDouble>(icethickness);
    721                                         xDelete<IssmDouble>(slopex);
    722                                         xDelete<IssmDouble>(slopey);
    723                                         break;
    724                                                 }
     639                                               IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
     640                                               int kk=0;
     641                                               for(k=0; k<numnodes;k++){
     642                                                       if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
     643                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
     644                                               }
     645                                               isminthicknessinelement=true;
     646                                               if(kk<numnodes && isminthicknessinelement){
     647                                                       for(k=0; k<numnodes;k++){
     648                                                               slope=fabs(slopex[k]);
     649                                                               if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
     650                                                               //slope_mean=slope_mean+slope;
     651                                                               if((atan(slope)*rad2deg)>slope_threshold) remove_debris=true;
     652                                                               //if((atan(slope)*rad2deg)>slope_threshold) debristhickness[k]=0.;
     653                                                       }
     654                                                       //if((atan(slope_mean)*rad2deg)>slope_threshold) remove_debris=true;
     655                                                       if(remove_debris){
     656                                                               for(k=0; k<numnodes;k++){
     657                                                                       debristhickness[k]=0.;
     658                                                               }
     659                                                       }
     660                                               }
     661                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
     662
     663                                               xDelete<IssmDouble>(debristhickness);
     664                                               xDelete<IssmDouble>(icethickness);
     665                                               xDelete<IssmDouble>(slopex);
     666                                               xDelete<IssmDouble>(slopey);
     667                                               break;
     668                                      }
    725669                                case 2:{
    726                                         IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
    727                                         IssmDouble gravity = element->FindParam(ConstantsGEnum);
    728                                         IssmDouble stress,rhod=1900.;
    729 
    730                                         for(k=0; k<numnodes;k++){
    731                                                 if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
    732                                         }
    733                                         isminthicknessinelement=true;
    734                                         if(isminthicknessinelement){
    735                                                 //stress=0;
    736                                                 int kk=0;
    737                                                 for(k=0; k<numnodes;k++){
    738                                                         if(onsurface[k]>0.5){
    739                                                                 slope=fabs(slopex[k]);
    740                                                                 if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
    741                                                                 stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
    742                                                                 if(stress>stress_threshold) debristhickness[k]=0.;
    743                                                                 //kk++;
    744                                                         }
    745                                                 }
    746                                                 /*if((stress/double(kk))>stress_threshold) remove_debris=true;
    747                                                   if(remove_debris){
    748                                                   for(k=0; k<numnodes;k++){
    749                                                   debristhickness[k]=0.;
    750                                                   }
    751                                                   }*/
    752                                         }
    753                                         //int finite_element = element->GetElementType();
    754                                         //element->AddInput(DebrisThicknessEnum,debristhickness,FINITEELEMENT);
    755                                         element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    756 
    757                                         xDelete<IssmDouble>(debristhickness);
    758                                         xDelete<IssmDouble>(icethickness);
    759                                         xDelete<IssmDouble>(slopex);
    760                                         xDelete<IssmDouble>(slopey);
    761                                         break;
    762                                                  }
     670                                               IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
     671                                               IssmDouble gravity = element->FindParam(ConstantsGEnum);
     672                                               IssmDouble stress,rhod=1900.;
     673                                               int kk=0;
     674                                               for(k=0; k<numnodes;k++){
     675                                                       if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
     676                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
     677                                               }
     678                                               isminthicknessinelement=true;
     679                                               if(kk<numnodes && isminthicknessinelement){
     680                                                       //stress=0;
     681                                                       IssmDouble stress_sum=0.;
     682                                                       for(k=0; k<numnodes;k++){
     683                                                               slope=fabs(slopex[k]);
     684                                                               if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
     685                                                               stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
     686                                                               //stress_sum=stress_sum+stress;
     687                                                               if(stress>stress_threshold) remove_debris=true;
     688                                                               //if(stress>stress_threshold) debristhickness[k]=0.;
     689                                                       }
     690                                                       //if((stress_sum/double(kk))>stress_threshold) remove_debris=true;
     691                                                       if(remove_debris){
     692                                                               for(k=0; k<numnodes;k++){
     693                                                                       debristhickness[k]=0.;
     694                                                               }
     695                                                       }
     696                                               }
     697                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
     698
     699                                               xDelete<IssmDouble>(debristhickness);
     700                                               xDelete<IssmDouble>(icethickness);
     701                                               xDelete<IssmDouble>(slopex);
     702                                               xDelete<IssmDouble>(slopey);
     703                                               xDelete<IssmDouble>(ls_active);
     704                                               break;
     705                                       }
    763706                                default: _error_("removalmodel "<<EnumToStringx(removalmodel)<<" not implemented yet");
    764707                        }
    … …  
    767710
    768711}/*}}}*/
    769 void           DebrisAnalysis::PreProcessing(FemModel* femmodel){/*{{{*/
     712void DebrisAnalysis::PreProcessing(FemModel* femmodel){/*{{{*/
    770713
    771714        if(VerboseSolution()) _printf0_("   Debris preprocessing\n");
    772715
    773         /*Intermediaries*/
    774         bool isdebrisdisplacement=true;
    775         int displacementmodel=1;
    776         int k,numnodes;
    777         int domaintype,dim;
    778         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    779         //femmodel->parameters->FindParam(&displacementmodel,DebrisDisplacementmodelEnum);
    780716        Element* element= NULL;
    781         IssmDouble *xyz_list = NULL;
    782         //IssmDouble top_normal[2];
    783 
    784         if(isdebrisdisplacement){
    785 
    786                 //if(displacementmodel==0){
    787                 // no displacement, do nothing
    788                 //}else{
    789                 // deterministic or random displacement
    790 
    791                 for(Object* & object : femmodel->elements->objects){
    792                         element=xDynamicCast<Element*>(object);
    793 
    794                         numnodes=element->GetNumberOfNodes();
    795                         IssmDouble* debristhickness= xNew<IssmDouble>(numnodes);
    796                         IssmDouble* icethickness= xNew<IssmDouble>(numnodes);
    797                         IssmDouble* slopex         = xNew<IssmDouble>(numnodes);
    798                         IssmDouble* slopey         = xNew<IssmDouble>(numnodes);
    799                         IssmDouble* vx             = xNew<IssmDouble>(numnodes);
    800                         IssmDouble* vy             = xNew<IssmDouble>(numnodes);
    801                         IssmDouble* surface        = xNew<IssmDouble>(numnodes);
    802                         IssmDouble* onsurface      = xNew<IssmDouble>(numnodes);
    803                         element->GetInputListOnNodes(vx,VxDebrisEnum);
    804                         element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
    805                         element->GetInputListOnNodes(icethickness,ThicknessEnum);
    806                         element->GetInputListOnNodes(surface,SurfaceEnum);
    807                         element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
    808                         element->GetVerticesCoordinates(&xyz_list);
    809 
    810                         dim=1;
    811                         element->GetInputListOnNodes(slopex,SurfaceSlopeXEnum);
    812                         if(domaintype!=Domain2DverticalEnum){
    813                                 element->GetInputListOnNodes(slopey,SurfaceSlopeYEnum);
    814                                 element->GetInputListOnNodes(vy,VyDebrisEnum);
    815                                 dim=2;
    816                         }
    817                         IssmDouble slope,rad2deg=180./M_PI; //=57.2958
    818                         IssmDouble h=10.,f;
    819                         IssmDouble debrissum;
    820                         IssmDouble dt = element->FindParam(TimesteppingTimeStepEnum);
    821                         bool displacedebris;
    822 
    823                         switch(displacementmodel){
    824                                 case 1:{
    825                                         IssmDouble slope_threshold = element->FindParam(DebrisRemovalSlopeThresholdEnum);
    826                                         IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
    827 
    828                                         bool isminthicknessinelement=false;
    829                                         for(k=0; k<numnodes;k++){
    830                                                 if(icethickness[k]<=(iceminthickness+0.01)) isminthicknessinelement=true;
    831                                         }
    832                                         if(isminthicknessinelement){
    833                                                 //do nothing
    834                                         }else{
    835                                                 debrissum=0.;
    836                                                 int test;
    837                                                 test=0;
    838                                                 for(k=0; k<numnodes;k++){
    839                                                         if(onsurface[k]>.5){
    840                                                                 slope=pow(slopex[k]*slopex[k],0.5);
    841                                                                 if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
    842                                                                 if((atan(slope)*rad2deg)>30.){
    843                                                                         f=0.5;
    844                                                                         test=test+0;
    845                                                                         debrissum=debrissum+debristhickness[k]*f;
    846                                                                         //displacedebris=true;
    847                                                                         //if(debristhickness[k]>1.e-6) vx[k]=vx[k]+10./31536000.;//*vx[k]/pow(pow(vx[k],2),0.5);
    848                                                                         //debristhickness[k]=debristhickness[k]*(1.-f);
    849                                                                 }
    850                                                         }
    851                                                 }
    852                                                 if(test>1){
    853                                                         test=test;
    854                                                 }else{
    855                                                         test=1;
    856                                                 }
    857                                                 //if(displacedebris){
    858                                                 int index=-1;
    859                                                 IssmDouble min=1e14;
    860                                                 for(k=0; k<numnodes;k++){
    861                                                         if(onsurface[k]>.5){
    862                                                                 if(surface[k]<min){
    863                                                                         index=k;
    864                                                                         min=surface[k];
    865                                                                 }
    866                                                         }
    867                                                 }
    868                                                 for(k=0; k<numnodes;k++){
    869                                                         if(onsurface[k]>.5){
    870                                                                 if(k==index){
    871                                                                         debristhickness[k]=debristhickness[k]+debrissum;
    872                                                                 }else{
    873                                                                         debristhickness[k]=debristhickness[k]-debrissum;
    874                                                                         if(debristhickness[k]<=0) debristhickness[k]=0;
    875                                                                 }
    876                                                                 //if(debristhickness[k]>10.) debristhickness[k]=10.;
    877                                                         }
    878                                                 }
    879                                                 //}
    880                                         }
    881 
    882                                         int finite_element = element->GetElementType();
    883                                         //element->AddInput(DebrisThicknessEnum,debristhickness,finite_element);
    884                                         //element->AddInput(VxDebrisEnum,vx,P1Enum);
    885                                         element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    886 
    887                                         xDelete<IssmDouble>(debristhickness);
    888                                         xDelete<IssmDouble>(icethickness);
    889                                         xDelete<IssmDouble>(vx);
    890                                         xDelete<IssmDouble>(vy);
    891                                         xDelete<IssmDouble>(slopex);
    892                                         xDelete<IssmDouble>(slopey);
    893                                         xDelete<IssmDouble>(surface);
    894                                         break;
    895                                                  }
    896                                 case 2:
    897                                         // Do nothing
    898 
    899                                 default: _error_("Debris displacement model "<<EnumToStringx(displacementmodel)<<" not implemented yet");
    900                         }
    901                 }
    902         }
    903         //}
    904 }/*}}}*/       
     717        for(Object* & object : femmodel->elements->objects){
     718                element=xDynamicCast<Element*>(object);
     719
     720                int numvertices = element->GetNumberOfVertices();
     721
     722                IssmDouble* vx = xNew<IssmDouble>(numvertices);
     723                IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
     724                IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
     725                IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
     726                IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
     727
     728                element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
     729                element->GetInputListOnVertices(&vx[0],VxDebrisEnum);
     730                element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
     731                element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
     732                element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
     733
     734                IssmDouble slope,rad2deg=180./M_PI; //=57.2958
     735                IssmDouble vslipx,rhod=1900.;
     736                IssmDouble gravity=element->FindParam(ConstantsGEnum);
     737                IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
     738                IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
     739
     740                int step;
     741                IssmDouble dt, maxv;
     742                IssmDouble yts=31536000.;
     743                femmodel->parameters->FindParam(&step,StepEnum);
     744                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     745
     746                bool isminthicknessinelement=false;
     747                for(int i=0;i<numvertices;i++){
     748                        if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
     749                }
     750                if(isminthicknessinelement){
     751                        //do nothing
     752                        for(int i=0;i<numvertices;i++){
     753                                if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;
     754                        }
     755                }else{
     756                        for(int i=0;i<numvertices;i++){
     757                                //if(onsurface[i]>.5){
     758                                slope=fabs(slopex[i]);
     759                                //if((atan(slope)*rad2deg)>25.){
     760                                //if(debristhickness[i]>0.01){
     761                                vslipx=1.0/yts;
     762                                //maxv=10.0/2./dt;
     763                                //vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts;
     764                                vx[i]=vx[i]+vslipx;
     765                                //debristhickness[i]=debristhickness[i];
     766                                //if(vx[i]>maxv) vx[i]=maxv;
     767                                //}
     768                                //}
     769                                //}
     770                        }
     771                }
     772                //if(step%100==0)   
     773                element->AddInput(VxDebrisEnum,vx,P1Enum);
     774
     775                /* Free resources */
     776                xDelete<IssmDouble>(debristhickness);
     777                xDelete<IssmDouble>(icethickness);
     778                xDelete<IssmDouble>(vx);
     779                xDelete<IssmDouble>(slopex);
     780                xDelete<IssmDouble>(onsurface);
     781        }
     782
     783}/*}}}*/

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

    r27468 r27491  
    23602360        Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);
    23612361        return (input->GetInputMax()<=0.);
     2362}
     2363/*}}}*/
     2364bool       Element::IsAllMinThicknessInElement(){/*{{{*/
     2365
     2366        IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum);
     2367
     2368        Input* input=this->GetInput(ThicknessEnum); _assert_(input);
     2369        if(input->GetInputMax()<=(minthickness+0.00000001)){
     2370                return true;
     2371        }
     2372        else{
     2373                return false;
     2374        }
    23622375}
    23632376/*}}}*/

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

    r27439 r27491  
    155155                bool               IsOceanInElement();
    156156                bool               IsOceanOnlyInElement();
     157                bool               IsAllMinThicknessInElement();
    157158                bool               IsLandInElement();
    158159                void               Ismip6FloatingiceMeltingRate();

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r27462 r27491  
    921921                        }
    922922                        if(isdebris){
     923                                analyses_temp[numanalyses++]=L2ProjectionBaseAnalysisEnum;
     924                                analyses_temp[numanalyses++]=SmbAnalysisEnum;
     925                                analyses_temp[numanalyses++]=ExtrudeFromTopAnalysisEnum;
    923926                                analyses_temp[numanalyses++]=DebrisAnalysisEnum;
    924927                        }

  • issm/trunk-jpl/src/c/cores/debris_core.cpp

    r27298 r27491  
    3434        if(VerboseSolution()) _printf0_("   computing debris transport\n");
    3535
    36         // We need surface and bed slopes for removal model
     36        // We need surface slopes for removal model
    3737        surfaceslope_core(femmodel);
    38         bedslope_core(femmodel);
    3938
    4039        /*Transport Debris*/
    4140        if(VerboseSolution()) _printf0_("   call computational core\n");
    42         //InputDuplicatex(femmodel,VxEnum,VxDebrisEnum);
    4341        femmodel->inputs->DuplicateInput(VxEnum,VxDebrisEnum);
    44         //InputDuplicatex(femmodel,VyEnum,VyDebrisEnum);
     42        if(domaintype!=Domain2DverticalEnum){
     43                femmodel->inputs->DuplicateInput(VyEnum,VyDebrisEnum); 
     44        }
     45        femmodel->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
     46        femmodel->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
     47        depthaverage_core(femmodel);
     48        if(domaintype!=Domain2DverticalEnum){
     49                femmodel->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
     50                femmodel->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
     51                depthaverage_core(femmodel);
     52        }
     53
    4554        debris_analysis = new DebrisAnalysis();
    46         //debris_analysis->PreCore(femmodel);
    47         //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
    48         //extrudefromtop_core(femmodel);       
    49 
    5055        debris_analysis->Core(femmodel);
    51         delete debris_analysis;
     56        delete debris_analysis; 
    5257
    5358        femmodel->parameters->SetParam(DebrisThicknessEnum,InputToExtrudeEnum);
    5459        extrudefromtop_core(femmodel); 
    55         //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
    56         //extrudefromtop_core(femmodel);
    5760
    5861        if(save_results) femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    59 
    6062        if(solution_type==DebrisSolutionEnum)femmodel->RequestedDependentsx();
    6163

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

    r27283 r27491  
    1010#include "../modules.h"
    1111
    12 void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
     12void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
    1313        /* activate/deactivate nodes for levelset method according to IceMaskNodeActivation */
    1414
    … …  
    1616        int nodeactivationmask = IceMaskNodeActivationEnum;
    1717        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
     18        if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
    1819
    1920
    … …  
    5960}/*}}}*/
    6061
    61 void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology){/*{{{*/
     62void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
    6263
    6364        /*Determine which node activation to construct*/
    6465        int nodeactivationmask = IceMaskNodeActivationEnum;
    6566        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
     67        if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
    6668
    6769        /*Initialize vector with number of vertices*/
    … …  
    8385                        }
    8486                }
    85         }
    86         else{
     87        }else if(isdebris){
     88                for(Object* & object : femmodel->elements->objects){
     89                        Element* element=xDynamicCast<Element*>(object);
     90                        if(element->IsIceInElement() && !element->IsAllMinThicknessInElement()){
     91                                int nbv = element->GetNumberOfVertices();
     92                                for(int iv=0;iv<nbv;iv++){
     93                                        vec_mask_ice->SetValue(element->vertices[iv]->Pid(),1.,INS_VAL);
     94                                }
     95                        }
     96                }
     97        }else{
    8798                for(Object* & object : femmodel->elements->objects){
    8899                        Element* element=xDynamicCast<Element*>(object);
    … …  
    101112        delete vec_mask_ice;
    102113}/*}}}*/
    103 void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
     114void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology,bool isdebris){/*{{{*/
    104115
    105116        /*Set configuration to levelset*/
    … …  
    107118                /*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/
    108119                femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
    109         }
    110         else{
     120        }else if(isdebris){
     121                femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);
     122        }else{
    111123                femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
    112124        }
    … …  
    115127        int nodeactivationmask = IceMaskNodeActivationEnum;
    116128        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
     129        if(isdebris) nodeactivationmask = DebrisMaskNodeActivationEnum;
    117130
    118131        /*Create vector on gset*/
    … …  
    137150                                xDelete<int>(glist_local);
    138151                        }
    139                 }
    140                 else{
     152                }else if(isdebris){
     153                        if(element->IsIceInElement() && !element->IsAllMinThicknessInElement()){
     154                                int numnodes = element->GetNumberOfNodes();
     155                                int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
     156                                int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
     157                                IssmDouble* ones = xNew<IssmDouble>(gsize_local);
     158                                for(int n=0;n<gsize_local;n++) ones[n] = 1.;
     159                                vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
     160                                xDelete<IssmDouble>(ones);
     161                                xDelete<int>(glist_local);
     162                        }
     163                }else{
    141164                        if(element->IsIceInElement()){
    142165                                int numnodes = element->GetNumberOfNodes();

  • issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h

    r27282 r27491  
    88#include "../../classes/classes.h"
    99
    10 void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false);
    11 void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false);
    12 void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false);
     10void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
     11void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
     12void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false,bool isdebris=false);
    1313#endif  /* _SETACTIVENODESLSMX_H*/

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

    r27466 r27491  
    522522        IssmDouble AccG=0.1;                    // m w.e. /100m
    523523        IssmDouble AccMax=1.;                    // m w.e.
    524         IssmDouble ReferenceElevation=2200.;     // m M&L
     524        IssmDouble ReferenceElevation;
    525525        IssmDouble AblationDays=120.;            //
    526526
    … …  
    537537        IssmDouble Alphad=0.07;            //              debris albedo                 0.07
    538538        IssmDouble Alphai=0.4;             //              ice ablbedo
     539        IssmDouble Alphaeff;
    539540        IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
    540541        IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
    … …  
    552553        IssmDouble smb,yts,z,debris;
    553554        IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
     555        bool isdebris;
     556        int domaintype;
     557        femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
    554558
    555559        /*Get material parameters and constants */
    … …  
    558562        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
    559563        PhiD=0.;
    560         //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
     564        if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
    561565
    562566        /* Loop over all the elements of this partition */
    … …  
    573577                /* Get inputs */
    574578                element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
     579                element->FindParam(&domaintype,DomainTypeEnum);         
    575580
    576581                /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
    577582                for(int v=0;v<numvertices;v++){
    578583
    579                         /*Get vertex elevation */
    580                         z=surfacelist[v];
    581 
    582                         /* Get debris cover */
    583                         debris=debriscover[v];
    584 
    585                         /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
    586                         IssmDouble n=debris/dk;
    587                         IssmDouble nmax=1000;
    588                         IssmDouble Alphaeff;
    589                         if(n>nmax){
    590                                 Alphaeff=Alphad;
    591                         } else {
    592                                 Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
    593                         }*/
    594 
    595                         // M&L
    596                        IssmDouble Alphaeff=Alphad;
    597 
    598                         /* compute smb */
    599                         for (int ismb=0;ismb<2;ismb++){
    600                                 if(ismb==0){
    601                                         // calc a reference smb to identify accum and melt region; debris only develops in ablation area
    602                                         debris=0.;
    603                                 }else{
    604                                         // only in the meltregime debris develops
    605                                         if(-MassBalanceCmDayDebris<0) debris=debriscover[v];
    606                                 }
    607                                 MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    608                                     (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+
    609                                     (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
    610                                     WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z-
    611                                     ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+
    612                                     ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
    613                                     WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    614                                     K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z-
    615                                     ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    616                                     Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.)
     584                        /*Get vertex elevation */
     585                        z=surfacelist[v];
     586
     587                        /*Get top element*/
     588                        //if(domaintype==Domain3DEnum){
     589
     590                        //}else{
     591                        //      Alphaeff=Alphad;
     592                        //      ReferenceElevation=2200.;     // m M&L                         
     593                        //}
     594
     595                        /* compute smb */
     596                        for (int ismb=0;ismb<2;ismb++){
     597                                if(ismb==0){
     598                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
     599                                        debris=0.;
     600                                        PhiD=0.;
     601                                }else{
     602                                        // only in the meltregime debris develops
     603                                        if(-MassBalanceCmDayDebris<1e-14) debris=debriscover[v];
     604                                }
     605                                if(debris<=0.) debris=0.;
     606                                IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
     607                                IssmDouble n=debris/dk;
     608                                IssmDouble nmax=1000;
     609                                IssmDouble Alphaeff;
     610                                if(n>nmax){
     611                                        Alphaeff=Alphad;
     612                                } else {
     613                                        Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
     614                                }
     615                                ReferenceElevation=3200.;     // m HEF
     616
     617
     618                                Alphaeff=Alphad;
     619                                ReferenceElevation=2200.;     // m M&L 
     620
     621                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
     622                                                        (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+
     623                                                        (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
     624                                                                        WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z-
     625                                                                                ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+
     626                                                                        ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
     627                                                                                        WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
     628                                                                        K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z-
     629                                                                                                                ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)*
     630                                                                                        Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.)
    617631                                    -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
    618632                        }

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

    r27469 r27491  
    211211syn keyword cConstant FrictionLinearizeEnum
    212212syn keyword cConstant FrictionPseudoplasticityExponentEnum
     213syn keyword cConstant FrictionU0Enum
    213214syn keyword cConstant FrictionThresholdSpeedEnum
    214215syn keyword cConstant FrictionVoidRatioEnum
    … …  
    868869syn keyword cConstant HydrologyWaterVyEnum
    869870syn keyword cConstant HydrologyMaskNodeActivationEnum
     871syn keyword cConstant DebrisMaskNodeActivationEnum
    870872syn keyword cConstant IceEnum
    871873syn keyword cConstant IceMaskNodeActivationEnum
    … …  
    17241726syn keyword cType Cfsurfacesquare
    17251727syn keyword cType Channel
     1728syn keyword cType classes
    17261729syn keyword cType Constraint
    17271730syn keyword cType Constraints
    … …  
    17301733syn keyword cType ControlInput
    17311734syn keyword cType Covertree
     1735syn keyword cType DatasetInput
    17321736syn keyword cType DataSetParam
    1733 syn keyword cType DatasetInput
    17341737syn keyword cType Definition
    17351738syn keyword cType DependentObject
    … …  
    17441747syn keyword cType ElementInput
    17451748syn keyword cType ElementMatrix
     1749syn keyword cType Elements
    17461750syn keyword cType ElementVector
    1747 syn keyword cType Elements
    17481751syn keyword cType ExponentialVariogram
    17491752syn keyword cType ExternalResult
    … …  
    17521755syn keyword cType Friction
    17531756syn keyword cType Gauss
     1757syn keyword cType GaussianVariogram
     1758syn keyword cType gaussobjects
    17541759syn keyword cType GaussPenta
    17551760syn keyword cType GaussSeg
    17561761syn keyword cType GaussTetra
    17571762syn keyword cType GaussTria
    1758 syn keyword cType GaussianVariogram
    17591763syn keyword cType GenericExternalResult
    17601764syn keyword cType GenericOption
    … …  
    17731777syn keyword cType IssmDirectApplicInterface
    17741778syn keyword cType IssmParallelDirectApplicInterface
     1779syn keyword cType krigingobjects
    17751780syn keyword cType Load
    17761781syn keyword cType Loads
    … …  
    17831788syn keyword cType Matice
    17841789syn keyword cType Matlitho
     1790syn keyword cType matrixobjects
    17851791syn keyword cType MatrixParam
    17861792syn keyword cType Misfit
    … …  
    17951801syn keyword cType Observations
    17961802syn keyword cType Option
     1803syn keyword cType Options
    17971804syn keyword cType OptionUtilities
    1798 syn keyword cType Options
    17991805syn keyword cType Param
    18001806syn keyword cType Parameters
    … …  
    18101816syn keyword cType Regionaloutput
    18111817syn keyword cType Results
     1818syn keyword cType Riftfront
    18121819syn keyword cType RiftStruct
    1813 syn keyword cType Riftfront
    18141820syn keyword cType SealevelGeometry
    18151821syn keyword cType Seg
    18161822syn keyword cType SegInput
     1823syn keyword cType Segment
    18171824syn keyword cType SegRef
    1818 syn keyword cType Segment
    18191825syn keyword cType SpcDynamic
    18201826syn keyword cType SpcStatic
    … …  
    18351841syn keyword cType Vertex
    18361842syn keyword cType Vertices
    1837 syn keyword cType classes
    1838 syn keyword cType gaussobjects
    1839 syn keyword cType krigingobjects
    1840 syn keyword cType matrixobjects
    18411843syn keyword cType AdjointBalancethickness2Analysis
    18421844syn keyword cType AdjointBalancethicknessAnalysis

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

    r27476 r27491  
    865865        HydrologyWaterVyEnum,
    866866        HydrologyMaskNodeActivationEnum,
     867        DebrisMaskNodeActivationEnum,
    867868        IceEnum,
    868869        IceMaskNodeActivationEnum,

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

    r27469 r27491  
    213213                case FrictionLinearizeEnum : return "FrictionLinearize";
    214214                case FrictionPseudoplasticityExponentEnum : return "FrictionPseudoplasticityExponent";
     215                case FrictionU0Enum : return "FrictionU0";
    215216                case FrictionThresholdSpeedEnum : return "FrictionThresholdSpeed";
    216217                case FrictionVoidRatioEnum : return "FrictionVoidRatio";
    … …  
    870871                case HydrologyWaterVyEnum : return "HydrologyWaterVy";
    871872                case HydrologyMaskNodeActivationEnum : return "HydrologyMaskNodeActivation";
     873                case DebrisMaskNodeActivationEnum : return "DebrisMaskNodeActivation";
    872874                case IceEnum : return "Ice";
    873875                case IceMaskNodeActivationEnum : return "IceMaskNodeActivation";

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

    r27469 r27491  
    204204syn keyword juliaConstC FrictionLinearizeEnum
    205205syn keyword juliaConstC FrictionPseudoplasticityExponentEnum
     206syn keyword juliaConstC FrictionU0Enum
    206207syn keyword juliaConstC FrictionThresholdSpeedEnum
    207208syn keyword juliaConstC FrictionVoidRatioEnum
    … …  
    861862syn keyword juliaConstC HydrologyWaterVyEnum
    862863syn keyword juliaConstC HydrologyMaskNodeActivationEnum
     864syn keyword juliaConstC DebrisMaskNodeActivationEnum
    863865syn keyword juliaConstC IceEnum
    864866syn keyword juliaConstC IceMaskNodeActivationEnum

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

    r27469 r27491  
    216216              else if (strcmp(name,"FrictionLinearize")==0) return FrictionLinearizeEnum;
    217217              else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum;
     218              else if (strcmp(name,"FrictionU0")==0) return FrictionU0Enum;
    218219              else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum;
    219220              else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum;
    … …  
    259260              else if (strcmp(name,"HydrologyarmaNumBreaks")==0) return HydrologyarmaNumBreaksEnum;
    260261              else if (strcmp(name,"HydrologyarmaNumParams")==0) return HydrologyarmaNumParamsEnum;
    261               else if (strcmp(name,"Hydrologyarmapolyparams")==0) return HydrologyarmapolyparamsEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum;
     265              if (strcmp(name,"Hydrologyarmapolyparams")==0) return HydrologyarmapolyparamsEnum;
     266              else if (strcmp(name,"HydrologyarmaTimestep")==0) return HydrologyarmaTimestepEnum;
    266267              else if (strcmp(name,"HydrologyAveraging")==0) return HydrologyAveragingEnum;
    267268              else if (strcmp(name,"HydrologyCavitySpacing")==0) return HydrologyCavitySpacingEnum;
    … …  
    382383              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    383384              else if (strcmp(name,"MaterialsEffectiveconductivityAveraging")==0) return MaterialsEffectiveconductivityAveragingEnum;
    384               else if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
     388              if (strcmp(name,"MaterialsHeatcapacity")==0) return MaterialsHeatcapacityEnum;
     389              else if (strcmp(name,"MaterialsLatentheat")==0) return MaterialsLatentheatEnum;
    389390              else if (strcmp(name,"MaterialsMeltingpoint")==0) return MaterialsMeltingpointEnum;
    390391              else if (strcmp(name,"MaterialsMixedLayerCapacity")==0) return MaterialsMixedLayerCapacityEnum;
    … …  
    505506              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
    506507              else if (strcmp(name,"StochasticForcingCovariance")==0) return StochasticForcingCovarianceEnum;
    507               else if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
     511              if (strcmp(name,"StochasticForcingDefaultDimension")==0) return StochasticForcingDefaultDimensionEnum;
     512              else if (strcmp(name,"StochasticForcingDimensions")==0) return StochasticForcingDimensionsEnum;
    512513              else if (strcmp(name,"StochasticForcingFields")==0) return StochasticForcingFieldsEnum;
    513514              else if (strcmp(name,"StochasticForcingIsEffectivePressure")==0) return StochasticForcingIsEffectivePressureEnum;
    … …  
    628629              else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
    629630              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    630               else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
     634              if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     635              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    635636              else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
    636637              else if (strcmp(name,"ThermalRequestedOutputs")==0) return ThermalRequestedOutputsEnum;
    … …  
    751752              else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
    752753              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    753               else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"CalvingAblationrate")==0) return CalvingAblationrateEnum;
     757              if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     758              else if (strcmp(name,"CalvingAblationrate")==0) return CalvingAblationrateEnum;
    758759              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    759760              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    … …  
    874875              else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
    875876              else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
    876               else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
     880              if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
     881              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
    881882              else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
    882883              else if (strcmp(name,"HydrologyReynolds")==0) return HydrologyReynoldsEnum;
    … …  
    891892              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    892893              else if (strcmp(name,"HydrologyMaskNodeActivation")==0) return HydrologyMaskNodeActivationEnum;
     894              else if (strcmp(name,"DebrisMaskNodeActivation")==0) return DebrisMaskNodeActivationEnum;
    893895              else if (strcmp(name,"Ice")==0) return IceEnum;
    894896              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    … …  
    996998              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    997999              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    998               else if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
    999               else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
     1003              if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
     1004              else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
     1005              else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
    10041006              else if (strcmp(name,"SealevelchangeGrot")==0) return SealevelchangeGrotEnum;
    10051007              else if (strcmp(name,"SealevelchangeGSatGravirot")==0) return SealevelchangeGSatGravirotEnum;
    … …  
    11191121              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    11201122              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    1121               else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    1122               else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
     1126              if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
     1127              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
     1128              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
    11271129              else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum;
    11281130              else if (strcmp(name,"SmbTemperaturesPresentday")==0) return SmbTemperaturesPresentdayEnum;
    … …  
    12421244              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    12431245              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    1244               else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    1245               else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
     1249              if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
     1250              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
     1251              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    12501252              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    12511253              else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    … …  
    13651367              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    13661368              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    1367               else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    1368               else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    13691369         else stage=12;
    13701370   }
    13711371   if(stage==12){
    1372               if (strcmp(name,"IntInput")==0) return IntInputEnum;
     1372              if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     1373              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     1374              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    13731375              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
    13741376              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    … …  
    14881490              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    14891491              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    1490               else if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
    1491               else if (strcmp(name,"HydrologySubTime")==0) return HydrologySubTimeEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
     1495              if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
     1496              else if (strcmp(name,"HydrologySubTime")==0) return HydrologySubTimeEnum;
     1497              else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
    14961498              else if (strcmp(name,"Hydrologypism")==0) return HydrologypismEnum;
    14971499              else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum;
    … …  
    16111613              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    16121614              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    1613               else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    1614               else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
     1618              if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
     1619              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
     1620              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    16191621              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    16201622              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    … …  
    17341736              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    17351737              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
    1736               else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    1737               else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    17381738         else stage=15;
    17391739   }
    17401740   if(stage==15){
    1741               if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
     1741              if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
     1742              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
     1743              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    17421744              else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
    17431745              else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
Note: See TracChangeset for help on using the changeset viewer.