Changeset 27852


Ignore:
Timestamp:
07/24/23 23:49:48 (20 months ago)
Author:
rueckamp
Message:

CHG: revert to revision 27846

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

Legend:

Unmodified
Added
Removed

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

    r27847 r27852  
    1010
    1111#define FINITEELEMENT P1Enum
    12 #define EPS 1e-8
    13 #define RHO 2
    14 //#define KAPPA 5
     12#define EPS 1e-14
    1513
    1614/*Model processing*/
    … …  
    117115        parameters->AddObject(iomodel->CopyConstantObject("md.debris.min_thickness",DebrisMinThicknessEnum));
    118116        parameters->AddObject(iomodel->CopyConstantObject("md.debris.packingfraction",DebrisPackingFractionEnum));
    119         parameters->AddObject(iomodel->CopyConstantObject("md.debris.max_displacementvelocity",DebrisMaxdisplacementvelocityEnum));
    120117        parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_slope_threshold",DebrisRemovalSlopeThresholdEnum));
    121118        parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_stress_threshold",DebrisRemovalStressThresholdEnum));
    … …  
    131128void           DebrisAnalysis::Core(FemModel* femmodel){/*{{{*/
    132129
    133         PreProcessing(femmodel);
     130        //PreProcessing(femmodel);
     131        //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
     132        //extrudefromtop_core(femmodel);
    134133        femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);       
    135134        solutionsequence_linear(femmodel);
    136         femmodel->inputs->DuplicateInput(DebrisThicknessEnum,DebrisThicknessOldEnum);
    137135        PostProcessing(femmodel);
    138        
     136
    139137}/*}}}*/
    140138void           DebrisAnalysis::PreCore(FemModel* femmodel){/*{{{*/
    … …  
    199197        Input* vy_input=NULL;
    200198        if(dim>1){vy_input = topelement->GetInput(VyDebrisEnum); _assert_(vy_input);}
    201         Input* slopex_input  = topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    202         Input* slopey_input  = topelement->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
    203         Input* thickness_input  = topelement->GetInput(DebrisThicknessEnum); _assert_(thickness_input);
    204199        h=topelement->CharacteristicLength();
    205         //IssmDouble hx,hy,hz;
    206         //topelement->ElementSizes(&hx,&hy,&hz);
    207         //h=sqrt(hx*hx+hy*hy);
    208200
    209201        /* Start  looping on the number of gaussian points: */
    … …  
    217209                vx_input->GetInputValue(&vx,gauss);
    218210                if(dim==2) vy_input->GetInputValue(&vy,gauss);
    219        
     211
    220212                D_scalar=gauss->weight*Jdet;
    221213                for(int i=0;i<numnodes;i++) for(int j=0;j<numnodes;j++) Ke->values[i*numnodes+j] += D_scalar*basis[i]*basis[j];
    … …  
    223215                /*Advection terms*/
    224216                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    225                 dvxdx=dvx[0];           
    226217                D_scalar=dt*gauss->weight*Jdet;
    227218                if(dim==2){
    228219                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
     220                        dvxdx=dvx[0];
    229221                        dvydy=dvy[1];
    230222                        for(int i=0;i<numnodes;i++){
    … …  
    237229                        }
    238230                }else{
     231                        dvxdx=dvx[0];
    239232                        for(int i=0;i<numnodes;i++){
    240233                                for(int j=0;j<numnodes;j++){
    … …  
    245238                }
    246239
    247                 /*******************************************************************/
    248                 /* Diffusion */
    249                 bool isdisplacement=false;
    250                 int step;
    251                 topelement->FindParam(&step,StepEnum);
    252                 IssmDouble slope_threshold;
    253                 topelement->FindParam(&slope_threshold,DebrisRemovalSlopeThresholdEnum);
    254                 IssmDouble rad2deg=180./M_PI;
    255                 IssmDouble kappa,f,smb,debristhickness;//,slopex;
    256                 IssmDouble M=1,C=-0.2,m,alpha_t_deg=25.0;
    257                 IssmDouble alpha_t_rad=tan(alpha_t_deg/rad2deg);
    258                 IssmDouble Diff,fraction;
    259                 Diff=3.2/3e7;
    260                 //Input* slopex_input=topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    261                 Input* smb_input=topelement->GetInput(SmbMassBalanceEnum); _assert_(smb_input);
    262                 Input* debristhickness_input=topelement->GetInput(DebrisThicknessEnum); _assert_(debristhickness_input);
    263 
    264                 /*if(isdisplacement){
    265 
    266                         slopex_input->GetInputValue(&slopex, gauss);
    267                         smb_input->GetInputValue(&smb, gauss);
    268                         debristhickness_input->GetInputValue(&debristhickness, gauss);
    269                         if((atan(fabs(slopex))*rad2deg)<alpha_t_deg){
    270                                 f=0.;
    271                         }else if((atan(fabs(slopex))*rad2deg)>=alpha_t_deg & fabs(slopex)<=((1-C)/M)){
    272                                 m=1/(((1-C)/M)-alpha_t_rad);
    273                                 f=m*fabs(slopex)-m*alpha_t_rad;
    274                         }else{
    275                                 f=1;
    276                         }
    277                         //f=1;
    278                         //kappa=-5.6e16*smb*debristhickness*f;
    279                         //kappa=debristhickness/h*4e9*f;
    280                         //kappa=14.2809e8*f; // 25°
    281                         kappa=247*f; // 35°
    282                         if(dim==2){
    283                                 for(int i=0;i<numnodes;i++){
    284                                         for(int j=0;j<numnodes;j++){
    285                                                 Ke->values[i*numnodes+j] +=  D_scalar*kappa*(
    286                                                                 dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
    287                                                                 );
    288                                         }
    289                                 }
    290                         }else{
    291                                 for(int i=0;i<numnodes;i++){
    292                                         for(int j=0;j<numnodes;j++){
    293                                                 Ke->values[i*numnodes+j] += D_scalar*kappa*(dbasis[0*numnodes+j]*dbasis[0*numnodes+i]);
    294                                         }
    295                                 }
    296                         }               
    297                 }
    298                 /*******************************************************************/
    299 
    300                 IssmDouble rho,KAPPA,slopex,slopey, thickness;
     240                IssmDouble rho;
    301241                if(FINITEELEMENT==P1Enum){
    302                         rho=RHO;
     242                        rho=2.;
    303243                }else if(FINITEELEMENT==P2Enum){
    304                         rho=RHO*2.;
     244                        rho=4.;
    305245                }
    306246
    … …  
    309249                        /*SSA*/
    310250                        if(dim==1){
    311                                 //vx_input->GetInputValue(&vx,gauss);
    312                                 vx_input->GetInputAverage(&vx);
     251                                vx_input->GetInputValue(&vx,gauss);
    313252                                D[0]=h/rho*fabs(vx);
    314253                        }else{
    … …  
    323262                        if(dim==1){
    324263                                vx_input->GetInputValue(&vx,gauss);
    325                                 //vx_input->GetInputAverage(&vx);
    326264                                vel=fabs(vx)+EPS;
    327265                                D[0] = h/(rho*vel)*vx*vx;
    328266                        }else{
    329267                                vx_input->GetInputValue(&vx,gauss);
    330                                 //vx_input->GetInputAverage(&vx);
    331268                                vy_input->GetInputValue(&vy,gauss);
    332                                 //vy_input->GetInputAverage(&vy);
    333269                                vel=sqrt(vx*vx+vy*vy)+EPS;
    334270                                D[0*dim+0]=h/(rho*vel)*vx*vx;
    … …  
    341277                        if(dim==1){
    342278                                vx_input->GetInputValue(&vx,gauss);
    343                                 //vx_input->GetInputAverage(&vx);
    344                                 //tau=h/(rho*max(fabs(vx),EPS));
    345                                 tau=h/(rho*(fabs(vx)+EPS));
     279                                tau=h/(rho*max(fabs(vx),EPS));
    346280                        }else{
    347281                                vx_input->GetInputValue(&vx,gauss);
    348282                                vy_input->GetInputValue(&vy,gauss);
    349                                 vx_input->GetInputAverage(&vx);
    350                                 vy_input->GetInputAverage(&vy);
    351                                 vel=sqrt(vx*vx+vy*vy)+EPS;
    352                                
    353                                 KAPPA=vel;
    354                                 // adapt diffusion coefficient from here: http://websrv.cs.umt.edu/isis/index.php/Solving_the_equation_for_thickness_evolution
    355                                 slopex_input->GetInputValue(&slopex,gauss);
    356                                 slopey_input->GetInputValue(&slopey,gauss);
    357                                 thickness_input->GetInputValue(&thickness,gauss);
    358                                 //KAPPA=vel*thickness/(pow(pow(slopex,2)+pow(slopey,2),0.5)+EPS);
    359                                 tau=h/(rho*vel)*(cosh(vel*h/2/KAPPA)/sinh(vel*h/2/KAPPA)-1./(vel*h/2/KAPPA));
    360 
    361                                 /*kappa=KAPPA/yts;
    362                                 if(vel*h/(3*2*kappa)<1){
    363                                         tau=pow(h,2)/(3*2*2*kappa);
    364                                 }
    365                                 else tau=h/(2*vel);*/
    366                                 //tau=dt*3;
     283                                tau=h/(rho*sqrt(vx*vx+vy*vy)+EPS);
    367284                        }
    368285                }
    … …  
    410327                                }
    411328
     329
    412330                                /*Advection matrix - part 2, B*/
    413331                                for(int i=0;i<numnodes;i++){
    … …  
    420338                                for(int i=0;i<numnodes;i++){
    421339                                        for(int j=0;j<numnodes;j++){
    422                                                 Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);
     340                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);;
    423341                                        }
    424342                                }
    … …  
    524442        Input* smb_input        = topelement->GetInput(SmbMassBalanceEnum);  _assert_(smb_input);
    525443        Input* thickness_input  = topelement->GetInput(DebrisThicknessEnum); _assert_(thickness_input);
    526         Input* slopex_input  = topelement->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    527         Input* slopey_input  = topelement->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
    528444        Input* vx_input  = topelement->GetInput(VxDebrisEnum);  _assert_(vx_input);
    529445        Input* vy_input=NULL;
    … …  
    532448        }
    533449        h=topelement->CharacteristicLength();
    534         //IssmDouble hx,hy,hz;
    535         //topelement->ElementSizes(&hx,&hy,&hz);
    536         //h=sqrt(hx*hx+hy*hy);
    537 
    538         IssmDouble rho,KAPPA,slopex,slopey;
     450
     451        IssmDouble rho;
    539452        if(FINITEELEMENT==P1Enum){
    540                 rho=RHO;
     453                rho=2.;
    541454        }else if(FINITEELEMENT==P2Enum){
    542                 rho=RHO*2.;
     455                rho=4.;
    543456        }
    544457
    … …  
    553466                thickness_input->GetInputValue(&thickness,gauss);
    554467                if(smb>0.){
    555                         psi=thickness-0*dt*smb*pf;
     468                        psi=thickness-0.*dt*smb*pf;
    556469                }else{
    557470                        psi=thickness-dt*smb*pf;
    558471                }
    559                 //psi=thickness-dt*smb*pf;
     472
    560473                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*psi*basis[i];
    561474
    … …  
    567480                        if(dim==1){
    568481                                vx_input->GetInputValue(&vx,gauss);
    569                                 //vx_input->GetInputAverage(&vx);
    570                                 //tau=h/(rho*max(fabs(vx),EPS));
    571                                 tau=h/(rho*(fabs(vx)+EPS));
     482                                tau=h/(rho*max(fabs(vx),EPS));
    572483
    573484                                /*Force vector - part 2*/
    … …  
    583494                                vx_input->GetInputValue(&vx,gauss);
    584495                                vy_input->GetInputValue(&vy,gauss);
    585                                 vx_input->GetInputAverage(&vx);
    586                                 vy_input->GetInputAverage(&vy);
    587496                                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
    588                                 dvydy=dvy[1];                           
    589 
    590                                 vel=sqrt(vx*vx+vy*vy)+EPS;
    591                                 // tuned value
    592                                 KAPPA=vel;
    593                                 // adapt diffusion coefficient from here: http://websrv.cs.umt.edu/isis/index.php/Solving_the_equation_for_thickness_evolution
    594                                 slopex_input->GetInputValue(&slopex,gauss);
    595                                 slopey_input->GetInputValue(&slopey,gauss);
    596                                 //KAPPA=vel*thickness/(pow(pow(slopex,2)+pow(slopey,2),0.5)+EPS);
    597                                 tau=h/(rho*vel)*(cosh(vel*h/2/KAPPA)/sinh(vel*h/2/KAPPA)-1./(vel*h/2/KAPPA));
    598 
    599                                 /*kappa=KAPPA/yts;
    600                                 if(vel*h/(3*2*kappa)<1){
    601                                         tau=pow(h,2)/(3*2*2*kappa);
    602                                 }
    603                                 else tau=h/(2*vel);*/
    604                                 //tau=dt*3;
     497                                vel=sqrt(vx*vx+vy*vy);
     498                                dvydy=dvy[1];
     499                                tau=h/(rho*vel+EPS);
    605500
    606501                                /*Force vector - part 2*/
    … …  
    636531        int numnodes = element->GetNumberOfNodes();
    637532        IssmDouble* newthickness  = xNew<IssmDouble>(numnodes);
    638         IssmDouble* icethickness  = xNew<IssmDouble>(numnodes);
    639        
    640533
    641534        /*Use the dof list to index into the solution vector: */
    642535        IssmDouble minthickness = element->FindParam(DebrisMinThicknessEnum);
    643536        element->GetDofListLocal(&ddoflist,NoneApproximationEnum,GsetEnum);
    644         element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    645537
    646538        for(int i=0;i<numnodes;i++){
    … …  
    658550        // Free resources
    659551        xDelete<IssmDouble>(newthickness);
    660         xDelete<IssmDouble>(icethickness);
    661552        xDelete<int>(ddoflist);
    662553        //*/
    … …  
    708599        femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum);
    709600        Element* element= NULL;
     601        Element*    topelement = NULL;
    710602
    711603        if(removalmodel==0){
    … …  
    724616                        IssmDouble* slopey         = xNew<IssmDouble>(numnodes);
    725617                        IssmDouble* onsurface      = xNew<IssmDouble>(numnodes);
    726                         //IssmDouble* ls_active      = xNew<IssmDouble>(numnodes);
    727                         element->GetInputListOnNodes(debristhickness,DebrisThicknessOldEnum);
     618                        IssmDouble* ls_active      = xNew<IssmDouble>(numnodes);
     619                        element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
    728620                        element->GetInputListOnNodes(icethickness,ThicknessEnum);
    729621                        element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
    730                         //element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
     622                        element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
    731623
    732624                        dim=1;
    … …  
    751643                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
    752644                                               }
    753                                                //isminthicknessinelement=true;
    754                                                //if(kk<numnodes && isminthicknessinelement){
    755                                                if(isminthicknessinelement){                                           
     645                                               isminthicknessinelement=true;
     646                                               if(kk<numnodes && isminthicknessinelement){
    756647                                                       for(k=0; k<numnodes;k++){
    757648                                                               slope=fabs(slopex[k]);
    … …  
    768659                                                       }
    769660                                               }
    770                                                element->AddInput(DebrisThicknessEnum,debristhickness,P1DGEnum);
     661                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    771662
    772663                                               xDelete<IssmDouble>(debristhickness);
    773                                                xDelete<IssmDouble>(onsurface);
    774664                                               xDelete<IssmDouble>(icethickness);
    775665                                               xDelete<IssmDouble>(slopex);
    776666                                               xDelete<IssmDouble>(slopey);
    777                                                //xDelete<IssmDouble>(ls_active);
    778667                                               break;
    779668                                       }
    780669                                case 2:{
    781                                                IssmDouble stress_threshold=element->FindParam(DebrisRemovalStressThresholdEnum);
    782                                                IssmDouble gravity=element->FindParam(ConstantsGEnum);
     670                                               IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
     671                                               IssmDouble gravity = element->FindParam(ConstantsGEnum);
    783672                                               IssmDouble stress,rhod=1900.;
    784673                                               int kk=0;
    … …  
    787676                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
    788677                                               }
    789                                                //isminthicknessinelement=true;
    790                                                //if(isminthicknessinelement){
     678                                               isminthicknessinelement=true;
    791679                                               if(kk<numnodes && isminthicknessinelement){
     680                                                       //stress=0;
     681                                                       IssmDouble stress_sum=0.;
    792682                                                       for(k=0; k<numnodes;k++){
    793683                                                               slope=fabs(slopex[k]);
    794684                                                               if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
    795                                                                stress=rhod*gravity*debristhickness[k]*slope;
    796                                                                //if(stress>stress_threshold) debristhickness[k]=0;
     685                                                               stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
     686                                                               //stress_sum=stress_sum+stress;
    797687                                                               if(stress>stress_threshold) remove_debris=true;
    798                                                         }
    799                                                         if(remove_debris){
    800                                                                for(k=0; k<numnodes;k++){
    801                                                                         debristhickness[k]=0.;
    802                                                                }
    803                                                         }
     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                                                       }
    804696                                               }
    805                                                element->AddInput(DebrisThicknessEnum,debristhickness,P1DGEnum);
     697                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    806698
    807699                                               xDelete<IssmDouble>(debristhickness);
    … …  
    809701                                               xDelete<IssmDouble>(slopex);
    810702                                               xDelete<IssmDouble>(slopey);
     703                                               xDelete<IssmDouble>(ls_active);
    811704                                               break;
    812705                                       }
    … …  
    821714        if(VerboseSolution()) _printf0_("   Debris preprocessing\n");
    822715
    823         int removalmodel, displacementmodel;
    824         femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum);
    825         femmodel->parameters->FindParam(&displacementmodel,DebrisDisplacementmodelEnum);
    826 
    827716        Element* element= NULL;
    828717        for(Object* & object : femmodel->elements->objects){
    829718                element=xDynamicCast<Element*>(object);
    830719
    831                 int numvertices=element->GetNumberOfVertices();
    832                 int numnodes=element->GetNumberOfNodes();
    833                 numvertices=numnodes;
     720                int numvertices = element->GetNumberOfVertices();
    834721
    835722                IssmDouble* vx = xNew<IssmDouble>(numvertices);
    836                 IssmDouble* vy = xNew<IssmDouble>(numvertices);
    837                 IssmDouble* vx_new = xNew<IssmDouble>(numvertices);
    838                 IssmDouble* vy_new = xNew<IssmDouble>(numvertices);             
    839723                IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
    840724                IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
    841                 IssmDouble* slopey         = xNew<IssmDouble>(numvertices);
    842                 IssmDouble* icethickness   = xNew<IssmDouble>(numvertices);
    843                 int dim=1,domaintype;
    844                 femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    845                 if(domaintype!=Domain2DverticalEnum){
    846                         dim=2;
    847                 }
     725                IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
     726                IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
     727
    848728                element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
    849                 element->GetInputListOnVertices(&vx[0],VxEnum);
     729                element->GetInputListOnVertices(&vx[0],VxDebrisEnum);
    850730                element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
     731                element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
    851732                element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    852                 if(dim==2){
    853                         element->GetInputListOnVertices(&vy[0],VyEnum);
    854                         element->GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum);                 
    855                 }
    856 
    857                 IssmDouble rad2deg=180./M_PI; //=57.2958
    858                 IssmDouble vslipx,vslipy,stress,rhod=1900.;
    859                 IssmDouble taud_plus, taud_minus, taudx, taudy, taud;
     733
     734                IssmDouble slope,rad2deg=180./M_PI; //=57.2958
     735                IssmDouble vslipx,rhod=1900.;
    860736                IssmDouble gravity=element->FindParam(ConstantsGEnum);
    861737                IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
    862                 IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
    863                 IssmDouble maxv = element->FindParam(DebrisMaxdisplacementvelocityEnum);
    864738                IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
    865                 IssmDouble hx,hy,hz;
    866                 element->ElementSizes(&hx,&hy,&hz);
    867739
    868740                int step;
    869                 IssmDouble dt;
     741                IssmDouble dt, maxv;
    870742                IssmDouble yts=31536000.;
    871743                femmodel->parameters->FindParam(&step,StepEnum);
    872744                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     745
    873746                bool isminthicknessinelement=false;
    874                 int kk=0;               
    875 
    876                 if(displacementmodel==1){
    877 
     747                for(int i=0;i<numvertices;i++){
     748                        if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
     749                }
     750                if(isminthicknessinelement){
     751                        //do nothing
    878752                        for(int i=0;i<numvertices;i++){
    879                                 if(icethickness[i]<=(iceminthickness+0.00001)){
    880                                         isminthicknessinelement=true;
    881                                         kk++;
    882                                 }
    883                         }
    884                         //isminthicknessinelement=true;
     753                                if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;
     754                        }
     755                }else{
    885756                        for(int i=0;i<numvertices;i++){
    886                         //isminthicknessinelement=true;
    887                         //if(icethickness[i]<=(iceminthickness+0.01)){
    888                         //if(kk<numnodes && isminthicknessinelement){
    889                        
    890                         /*if(isminthicknessinelement && kk==numvertices){
    891                                 vx[i]=0.;
    892                                 if(dim==2) vy[i]=0.;
    893                         }else if(isminthicknessinelement && kk<numvertices){
    894                                 vx[i]=0.;
    895                                 if(dim==2) vy[i]=0.;
    896                         }else*/
    897                         /*if(isminthicknessinelement){
    898                                 if(kk==numvertices){
    899                                         vx[i]=0.;
    900                                         if(dim==2) vy[i]=0.;
    901                                 }
    902                         }else{*/
    903                                 if(~isminthicknessinelement){
    904                                                 if(removalmodel==1){
    905                                                         taud_plus=slope_threshold; taud_minus=slope_threshold/2;
    906                                                 }else if(removalmodel==2){
    907                                                         taud_plus=stress_threshold; taud_minus=stress_threshold/2;
    908                                                 }
    909                                                 taudx=rhod*gravity*debristhickness[i]*slopex[i]; // 0.1
    910                                                 if(dim==2) taudy=rhod*gravity*debristhickness[i]*slopey[i];
    911                                                 vslipx=-0.1*taudx/yts;//+0./yts*fabs(vx[i])/(vx[i]+1e-10);
    912                                                 if(dim==2) vslipy=-0.1*taudy/yts;//+0./yts*fabs(vy[i])/(vy[i]+1e-10);
    913 
    914                                                 if(removalmodel==1){
    915                                                         taud=atan(slopex[i]*slopex[i]+slopey[i]*slopey[i])*rad2deg;
    916                                                 }else if(removalmodel==2){
    917                                                         taud=sqrt(taudx*taudx+taudy*taudy);
    918                                                 }
    919                                                
    920                                                 if(taud>=taud_plus){
    921                                                         vx[i]=vx[i]+vslipx;
    922                                                         if(dim==2) vy[i]=vy[i]+vslipy;
    923                                                 }else if(taud>=taud_minus & taud<taud_plus){
    924                                                         //vx[i]=vx[i]+vslipx*(1-(taud_plus-fabs(taudx))/(taud_plus-taud_minus));
    925                                                         //if(dim==2) vy[i]=vy[i]+vslipy*(1-(taud_plus-fabs(taudy))/(taud_plus-taud_minus));
    926                                                         vx[i]=vx[i]+vslipx*(1+cos(PI*fabs(taudx)/(taud_plus-taud_minus)))/2.;
    927                                                         if(dim==2) vy[i]=vy[i]+vslipy*(1+cos(PI*fabs(taudy)/(taud_plus-taud_minus)))/2.;
    928                                                 }else if(taud<taud_minus){
    929                                                         vx[i]=vx[i];
    930                                                         if(dim==2) vy[i]=vy[i];
    931                                                 }
    932                                                 if(fabs(vx[i])>maxv/yts) vx[i]=maxv/yts*fabs(vx[i])/(vx[i]+1e-10);
    933                                                 if(fabs(vy[i])>maxv/yts) vy[i]=maxv/yts*fabs(vy[i])/(vy[i]+1e-10);
    934                                 }
    935                         }
    936                         element->AddInput(VxDebrisEnum,vx,P1Enum);
    937                         if(dim==2) element->AddInput(VyDebrisEnum,vy,P1Enum);
    938                 }
     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
    939775                /* Free resources */
    940776                xDelete<IssmDouble>(debristhickness);
    941777                xDelete<IssmDouble>(icethickness);
    942778                xDelete<IssmDouble>(vx);
    943                 xDelete<IssmDouble>(vx_new);
    944779                xDelete<IssmDouble>(slopex);
    945                 xDelete<IssmDouble>(vy);
    946                 xDelete<IssmDouble>(vy_new);
    947                 xDelete<IssmDouble>(slopey);
    948         }
    949 }/*}}}*/
     780                xDelete<IssmDouble>(onsurface);
     781        }
     782
     783}/*}}}*/

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

    r27847 r27852  
    242242                        }
    243243                        break;
    244                 case SMBdebrisEvattEnum:
    245                         iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
    246                         iomodel->FetchDataToInput(inputs,elements,"md.smb.s0t",SmbS0tEnum);
    247                         iomodel->FetchDataToInput(inputs,elements,"md.smb.snowheight",SmbSnowheightEnum);
    248                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
    249                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyprecipitation",SmbPrecipitationEnum);
    250                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydsradiation",SmbMonthlydsradiationEnum);
    251                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlydlradiation",SmbMonthlydlradiationEnum);
    252                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlywindspeed",SmbMonthlywindspeedEnum);
    253                         iomodel->FetchDataToDatasetInput(inputs,elements,"md.smb.monthlyairhumidity",SmbMonthlyairhumidityEnum);
    254                         iomodel->FetchDataToInput(inputs,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
    255                         iomodel->FetchDataToInput(inputs,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
    256                         iomodel->FetchDataToInput(inputs,elements,"md.smb.dsradiation_anomaly",SmbDsradiationAnomalyEnum);
    257                         iomodel->FetchDataToInput(inputs,elements,"md.smb.dlradiation_anomaly",SmbDlradiationAnomalyEnum);
    258                         iomodel->FetchDataToInput(inputs,elements,"md.smb.windspeed_anomaly",SmbWindspeedAnomalyEnum);
    259                         iomodel->FetchDataToInput(inputs,elements,"md.smb.airhumidity_anomaly",SmbAirhumidityAnomalyEnum);
    260                         break;
     244                case SMBdebrisMLEnum:
     245                        iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
     246                        break;
    261247                default:
    262248                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    … …  
    487473                        /*Nothing to add to parameters*/
    488474                        break;
    489                 case SMBdebrisEvattEnum:
    490                         /*Nothing to add to parameters*/
    491                         break;
     475                case SMBdebrisMLEnum:
     476                        /*Nothing to add to parameters*/
     477                        break;
    492478                default:
    493479                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    … …  
    587573                        #endif //_HAVE_SEMIC_
    588574                        break;
    589                 case SMBdebrisEvattEnum:
    590                         if(VerboseSolution())_printf0_("        call smb Evatt debris module\n");
    591                         SmbDebrisEvattx(femmodel);
    592                         break;
     575                case SMBdebrisMLEnum:
     576                        if(VerboseSolution())_printf0_("        call smb debris Mayer & Liculli module\n");
     577                        SmbDebrisMLx(femmodel);
     578                        break;
    593579                default:
    594580                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");

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

    r27847 r27852  
    2323#include "../Inputs/ArrayInput.h"
    2424#include "../Inputs/IntArrayInput.h"
    25 #include <cstdlib>
    2625/*}}}*/
    2726#define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/
    … …  
    23512350
    23522351        IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum);
    2353        
     2352
    23542353        Input* input=this->GetInput(ThicknessEnum); _assert_(input);
    2355         if(input->GetInputMin()<=(minthickness+0.00000001)){
     2354        if(input->GetInputMax()<=(minthickness+0.00000001)){
    23562355                return true;
    23572356        }
    … …  
    34013400                        monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
    34023401                        dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month);
    3403                         monthlyprec[iv*12+month]=monthlyprec[iv*12+month];//*yts;
     3402                        monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
    34043403                }
    34053404        }
    … …  
    35953594        Gauss* gauss=this->NewGauss();
    35963595        for(int month=0;month<12;month++){
     3596
    35973597                for(int iv=0;iv<NUM_VERTICES;iv++){
    35983598                        gauss->GaussVertex(iv);
    … …  
    37273727        xDelete<IssmDouble>(smbcorr);
    37283728        xDelete<IssmDouble>(melt_star);
    3729 }
    3730 /*}}}*/
    3731 void       Element::SmbDebrisEvatt(){/*{{{*/
    3732 
    3733         const int NUM_VERTICES          = this->GetNumberOfVertices();
    3734         const int NUM_VERTICES_MONTHS_PER_YEAR  = NUM_VERTICES * 365;
    3735 
    3736         int             i,vertexlids[MAXVERTICES];;
    3737         IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);         // surface mass balance
    3738         IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);                // melting comp. of surface mass balance
    3739         IssmDouble* summermelt=xNew<IssmDouble>(NUM_VERTICES);
    3740         IssmDouble* albedo=xNew<IssmDouble>(NUM_VERTICES);              // melting comp. of surface mass balance
    3741         IssmDouble* summeralbedo=xNew<IssmDouble>(NUM_VERTICES);
    3742         IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES);                // accuumulation comp. of surface mass balance
    3743         IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3744         IssmDouble* monthlyprecip=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3745         IssmDouble* monthlylw=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3746         IssmDouble* monthlysw=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3747         IssmDouble* monthlywind=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3748         IssmDouble* monthlyhumidity=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
    3749         IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
    3750         IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);      // precip anomaly
    3751         IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);      // temperature anomaly
    3752         IssmDouble* lw_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3753         IssmDouble* sw_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3754         IssmDouble* wind_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3755         IssmDouble* humidity_ampl=xNew<IssmDouble>(NUM_VERTICES);
    3756        
    3757         IssmDouble* surface=xNew<IssmDouble>(NUM_VERTICES);             // actual surface height
    3758         IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);         // reference elevation
    3759         IssmDouble* snowheight=xNew<IssmDouble>(NUM_VERTICES);
    3760         IssmDouble* debriscover=xNew<IssmDouble>(NUM_VERTICES);
    3761         IssmDouble rho_water,rho_ice,Tf,debris,debris_here;
    3762         IssmDouble qlaps,rlaps,dslaps,dllaps,windspeedlaps,humiditylaps,Tm;
    3763         IssmDouble inv_twelve=1./365.;                          //factor for monthly average
    3764         IssmDouble time,yts,time_yr,lambda;
    3765         IssmDouble MonthlyMelt,CleanIceMonthlyMelt, CumMonthlyMelt=0,CleanIceMelt,CumMonthlySummerMelt=0;
    3766         IssmDouble MeanAlbedo=0, MeanSummerAlbedo=0;
    3767         bool isdebris,ismonthly,isAnderson,iscryokarst;
    3768         this->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
    3769         this->parameters->FindParam(&ismonthly,SmbDebrisIsMonthlyEnum);
    3770         this->parameters->FindParam(&isAnderson,SmbDebrisIsAndersonEnum);
    3771         this->parameters->FindParam(&iscryokarst,SmbDebrisIsCryokarstEnum);
    3772         IssmDouble PhiD=0.,p;
    3773         IssmDouble icealbedo=this->FindParam(SmbIcealbedoEnum);
    3774         IssmDouble snowalbedo=this->FindParam(SmbSnowalbedoEnum);
    3775         IssmDouble debrisalbedo=this->FindParam(SmbDebrisalbedoEnum);
    3776         IssmDouble Lm=this->FindParam(MaterialsLatentheatEnum);
    3777         IssmDouble D0=this->FindParam(SmbDebrisAndersonD0Enum);
    3778         int step;
    3779         this->FindParam(&step,StepEnum);
    3780        
    3781         // cryokarst
    3782         int dim=1,domaintype;
    3783         this->parameters->FindParam(&domaintype,DomainTypeEnum);
    3784         if(domaintype!=Domain2DverticalEnum) dim=2;
    3785         IssmDouble taud_plus=110e3, taud_minus=60e3;
    3786         IssmDouble taud, slope, gravity, taudx, taudy;
    3787         this->parameters->FindParam(&gravity,ConstantsGEnum);           
    3788         IssmDouble* slopex         = xNew<IssmDouble>(NUM_VERTICES);
    3789         IssmDouble* slopey         = xNew<IssmDouble>(NUM_VERTICES);
    3790         IssmDouble* icethickness   = xNew<IssmDouble>(NUM_VERTICES);
    3791 
    3792         /*Get material parameters :*/
    3793         rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
    3794         rho_ice=this->FindParam(MaterialsRhoIceEnum);
    3795         IssmDouble sconv=(rho_water/rho_ice);
    3796         Tf=this->FindParam(MaterialsMeltingpointEnum);
    3797 
    3798         /*Get parameters for height corrections*/
    3799         qlaps=this->FindParam(SmbDesfacEnum); // comment MR; on alpine galciers we dont have the dertification effect
    3800         rlaps=this->FindParam(SmbRlapsEnum);
    3801         dslaps=this->FindParam(SmbSWlapsEnum);
    3802         dllaps=this->FindParam(SmbLWlapsEnum);
    3803         windspeedlaps=this->FindParam(SmbWindspeedlapsEnum);
    3804         humiditylaps=this->FindParam(SmbHumiditylapsEnum);
    3805 
    3806         /* Get time */
    3807         this->parameters->FindParam(&time,TimeEnum);
    3808         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    3809         time_yr=floor(time/yts)*yts;
    3810 
    3811         /*Get inputs*/
    3812         DatasetInput* tempmonthly     =this->GetDatasetInput(SmbMonthlytemperaturesEnum); _assert_(tempmonthly);
    3813         DatasetInput* precipmonthly   =this->GetDatasetInput(SmbPrecipitationEnum);       _assert_(precipmonthly);
    3814         DatasetInput* lwmonthly       =this->GetDatasetInput(SmbMonthlydlradiationEnum); _assert_(lwmonthly);
    3815         DatasetInput* swmonthly       =this->GetDatasetInput(SmbMonthlydsradiationEnum);       _assert_(swmonthly);     
    3816         DatasetInput* windmonthly     =this->GetDatasetInput(SmbMonthlywindspeedEnum); _assert_(windmonthly);
    3817         DatasetInput* humiditymonthly =this->GetDatasetInput(SmbMonthlyairhumidityEnum);       _assert_(humiditymonthly);
    3818 
    3819         /*loop over vertices: */
    3820         Gauss* gauss=this->NewGauss();
    3821         for(int month=0;month<365;month++){
    3822                 for(int iv=0;iv<NUM_VERTICES;iv++){
    3823                         gauss->GaussVertex(iv);
    3824                         tempmonthly->GetInputValue(&monthlytemperatures[iv*365+month],gauss,month);
    3825                         monthlytemperatures[iv*365+month]=monthlytemperatures[iv*365+month]-Tf; // conversion from Kelvin to celcius for PDD module
    3826                         precipmonthly->GetInputValue(&monthlyprecip[iv*365+month],gauss,month);
    3827                         monthlyprecip[iv*365+month]=monthlyprecip[iv*365+month]*yts; // from m/s to m/a
    3828                         lwmonthly->GetInputValue(&monthlylw[iv*365+month],gauss,month);
    3829                         swmonthly->GetInputValue(&monthlysw[iv*365+month],gauss,month);
    3830                         windmonthly->GetInputValue(&monthlywind[iv*365+month],gauss,month);
    3831                         humiditymonthly->GetInputValue(&monthlyhumidity[iv*365+month],gauss,month);
    3832                 }
    3833         }
    3834 
    3835         /*Recover info at the vertices: */
    3836         GetInputListOnVertices(&surface[0],SurfaceEnum);
    3837         GetInputListOnVertices(&s0t[0],SmbS0tEnum);
    3838         GetInputListOnVertices(&snowheight[0],SmbSnowheightEnum);
    3839         GetInputListOnVertices(&debriscover[0],DebrisThicknessEnum);
    3840         GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum);
    3841         GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum);
    3842         GetInputListOnVertices(&lw_ampl[0],SmbDsradiationAnomalyEnum);
    3843         GetInputListOnVertices(&sw_ampl[0],SmbDlradiationAnomalyEnum);
    3844         GetInputListOnVertices(&wind_ampl[0],SmbWindspeedAnomalyEnum);
    3845         GetInputListOnVertices(&humidity_ampl[0],SmbAirhumidityAnomalyEnum);
    3846         if(iscryokarst){
    3847                 GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
    3848                 GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    3849                 if(dim==2){
    3850                         GetInputListOnVertices(&slopey[0],SurfaceSlopeYEnum);                 
    3851                 }
    3852                 taudx=rho_ice*gravity*icethickness[i]*slopex[i];
    3853                 if(dim==2) taudy=rho_ice*gravity*icethickness[i]*slopey[i];
    3854                 taud=sqrt(taudx*taudx+taudy*taudy);
    3855         }
    3856         IssmDouble Alphaeff,Alphaeff_cleanice;
    3857 
    3858         /*measure the surface mass balance*/
    3859         for (int iv = 0; iv<NUM_VERTICES; iv++){
    3860 
    3861                 IssmDouble st=(surface[iv]-s0t[iv])/1000.;
    3862                
    3863                 int ismb_end=1;
    3864                 if(isdebris & !isAnderson) ismb_end=2;
    3865                 for (int ismb=0;ismb<ismb_end;ismb++){
    3866                         if(ismb==0){
    3867                                 // calc a reference smb to identify accum and melt region; debris only develops in ablation area
    3868                                 debris=0.;
    3869                                 PhiD=0.;
    3870                                 if(isAnderson) debris_here=debriscover[iv]; // store debris for later
    3871                         }else{
    3872                                 // debris only develops in ablation area
    3873                                 /*if((accu[iv]/yts-CleanIceMelt)<(-1e-2)/yts){
    3874                                         debris=debriscover[iv];
    3875                                 }else{
    3876                                         debris=0.;
    3877                                 }*/
    3878                                 debris=0.;
    3879                                 if(debris<=0.) debris=0.;
    3880                                 if(isdebris) PhiD=FindParam(DebrisPackingFractionEnum);
    3881                                 CumMonthlyMelt=0;
    3882                                 CumMonthlySummerMelt=0;
    3883                                 debris_here=debriscover[iv];
    3884                         }
    3885 
    3886                         /* Now run the debris part */
    3887                        
    3888                         // Climate inputs
    3889                         IssmDouble Tm;          // C air temperature
    3890                         IssmDouble In;          // Wm^-2 incoming long wave
    3891                         IssmDouble Q;           // Wm^-2 incoming short wave
    3892                         IssmDouble Um;          // ms^-1 measured wind speed
    3893                         IssmDouble Humidity;    // relative humidity
    3894                         IssmDouble P;           // precip
    3895 
    3896                         // other parameters
    3897                         IssmDouble Qh=0.006;   // kg m^-3      saturated humidity level // not used
    3898                         IssmDouble Qm=0.8*Qh;  // kg m^-3      measured humiditiy level // not used
    3899                         IssmDouble Rhoaa=1.22; // kgm^-3       air densitiy
    3900                         IssmDouble K=0.585;    // Wm^-1K^-1    thermal conductivity          0.585
    3901                         IssmDouble Xr=0.01;    // ms^-1        surface roughness             0.01
    3902                         IssmDouble Ustar=0.16; // ms^-1        friction velocity             0.16
    3903                         IssmDouble Ca=1000;    // jkg^-1K^-1   specific heat capacity of air
    3904                         IssmDouble Lv=2.50E+06;// jkg^-1K^-1   latent heat of evaporation
    3905                         IssmDouble Eps=0.95;   //              thermal emissivity
    3906                         IssmDouble Sigma=5.67E-08;// Wm^-2K^-4    Stefan Boltzmann constant
    3907                         IssmDouble Gamma=180.;    // m^-1         wind speed attenuation        234
    3908                
    3909                         // Calculate effective albedo
    3910                         IssmDouble dk=1e-2;
    3911                         IssmDouble n=debris/dk;
    3912                         IssmDouble nmax=(25*25)/2 * 0.02/dk;
    3913                         IssmDouble Alphaeff,Alphaeff_cleanice;
    3914                         IssmDouble mean_ela,delta=2000;
    3915                        
    3916                         // compute cleanice albedo based on previous SMB distribution
    3917                         if(step==1){
    3918                                 mean_ela=3000; //FIXME
    3919                         }else{
    3920                                 mean_ela=FindParam(SmbMeanElaEnum);
    3921                         }
    3922                         Alphaeff_cleanice=icealbedo+(snowalbedo-icealbedo)*(1+tanh(PI*(surface[iv]-mean_ela)/delta))/2;
    3923                         Alphaeff=Alphaeff_cleanice; // will be updated below
    3924                        
    3925                         yearlytemperatures[iv]=0; // tmp
    3926                         accu[iv]=0.;
    3927                         for (int imonth=0;imonth<365;imonth++) {
    3928 
    3929                                 Tm=monthlytemperatures[iv*365+imonth]-st*rlaps;//+t_ampl[iv];//+(rand()%10-5)/5;
    3930                                 In=monthlylw[iv*365+imonth]-st*dllaps+lw_ampl[iv];
    3931                                 Q=monthlysw[iv*365+imonth]+st*dslaps+sw_ampl[iv];
    3932                                 Humidity=monthlyhumidity[iv*365+imonth]-st*humiditylaps+humidity_ampl[iv];
    3933                                 Um=monthlywind[iv*365+imonth]-st*windspeedlaps+wind_ampl[iv];
    3934                                 P=(qlaps*st*monthlyprecip[iv*365+imonth]+monthlyprecip[iv*365+imonth]+p_ampl[iv])*sconv/365.; // convert precip from w.e. -> i.e
    3935                                 /*Partition of precip in solid and liquid parts */
    3936                                 IssmDouble temp_plus=1;
    3937                                 IssmDouble temp_minus=-1.;
    3938                                 IssmDouble frac_solid;
    3939                                 if(Tm>=temp_plus){
    3940                                         frac_solid=0;
    3941                                 }else if(Tm<=temp_minus){
    3942                                         frac_solid=1;
    3943                                 }else{
    3944                                         frac_solid=1*(1-cos(PI*(temp_plus-Tm)/(temp_plus-temp_minus)))/2;
    3945                                 }
    3946                
    3947                                 /*Get yearly temperatures and accumulation */
    3948                                 yearlytemperatures[iv]=yearlytemperatures[iv]+((monthlytemperatures[iv*365+imonth]-rlaps*st+Tf+t_ampl[iv]))/365; // Has to be in Kelvin
    3949                                 accu[iv]=accu[iv]+P*frac_solid;
    3950                                 if(yearlytemperatures[iv]>Tf) yearlytemperatures[iv]=Tf;
    3951                                
    3952                                 CleanIceMonthlyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    3953                                         Q*(1.-Alphaeff)+
    3954                                         (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+
    3955                                         ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    3956                                         K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    3957                                         rho_ice*Lm*Ustar)/(((Um
    3958                                         -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)));
    3959                                 if(CleanIceMonthlyMelt<0) CleanIceMonthlyMelt=0.;
    3960                                 MonthlyMelt=CleanIceMonthlyMelt;
    3961                                
    3962                                 if(ismb==1){
    3963                                         snowheight[iv]=snowheight[iv]+(P-CleanIceMonthlyMelt*yts/365);
    3964                                         if(snowheight[iv]<=0) snowheight[iv]=0.;
    3965                                         if(snowheight[iv]<=0.0001){
    3966                                                 p=debris_here*PhiD/(2*0.2*0.01); //Eq. 51 from Evatt et al 2015 without source term g*t
    3967                                                 if(p>1.) p=1.;
    3968                                                 if(p>=0.999){
    3969                                                         Alphaeff=debrisalbedo;
    3970                                                 } else {
    3971                                                         Alphaeff=Alphaeff_cleanice+p*(debrisalbedo-Alphaeff_cleanice);
    3972                                                 }
    3973                                                 debris=debris_here;
    3974                                                 MonthlyMelt=((In-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    3975                                                         Q*(1.-Alphaeff)+
    3976                                                         (Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))*Tm)/((1-PhiD)*rho_ice*Lm)/(1.+
    3977                                                         ((Rhoaa*Ca*Ustar*Ustar)/(Um-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    3978                                                         K*debris)-(Lv*Ustar*Ustar*((Humidity))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    3979                                                         rho_ice*Lm*Ustar)/(((Um-2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)));
    3980                                                 if(MonthlyMelt<0) MonthlyMelt=0.;
    3981                                                 MeanSummerAlbedo=MeanSummerAlbedo+Alphaeff;
    3982                                                 CumMonthlySummerMelt=CumMonthlySummerMelt+MonthlyMelt/365;
    3983                                         }
    3984                                 }
    3985                                 CumMonthlyMelt=CumMonthlyMelt+MonthlyMelt/365;
    3986                         }
    3987                         MeanAlbedo=MeanAlbedo+Alphaeff;
    3988                         if(ismb==0) CleanIceMelt=CumMonthlyMelt;
    3989                 }
    3990 
    3991                 if(iscryokarst){
    3992                         if(taud>=taud_plus){
    3993                                 lambda=0;
    3994                         }else if(taud>=taud_minus & taud<taud_plus){
    3995                                 lambda=0.1*(1-cos(PI*(taud_plus-taud)/(taud_plus-taud_minus)))/2;
    3996                         }else if(taud<taud_minus){
    3997                                 lambda=0.1;
    3998                         }                       
    3999                 }
    4000                
    4001                 // update values
    4002                 melt[iv]=CumMonthlyMelt; // is already in m/s
    4003                 accu[iv]=accu[iv]/yts;
    4004                 if(isAnderson){
    4005                         smb[iv]=(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4006                         if(iscryokarst){
    4007                                 smb[iv]=lambda*(accu[iv]-melt[iv])+(1-lambda)*(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4008                         }else{
    4009                                 smb[iv]=(accu[iv]-melt[iv])*D0/(D0+debris_here);
    4010                         }
    4011                 }else{
    4012                         if(iscryokarst){
    4013                                 smb[iv]=lambda*(accu[iv]-CleanIceMelt)+(1-lambda)*(accu[iv]-melt[iv]);
    4014                         }else{
    4015                                 smb[iv]=(accu[iv]-melt[iv]);
    4016                         }
    4017                 }
    4018                 albedo[iv]=MeanAlbedo;
    4019                 summeralbedo[iv]=MeanSummerAlbedo;
    4020                 summermelt[iv]=CumMonthlySummerMelt;
    4021         }
    4022        
    4023         this->AddInput(SmbMassBalanceEnum,smb,P1Enum);
    4024         this->AddInput(SmbAccumulationEnum,accu,P1Enum);
    4025         this->AddInput(SmbMeltEnum,melt,P1Enum);
    4026         this->AddInput(SmbSummerMeltEnum,summermelt,P1Enum);
    4027         this->AddInput(SmbSnowheightEnum,snowheight,P1Enum);
    4028         this->AddInput(SmbAlbedoEnum,albedo,P1Enum);
    4029         this->AddInput(SmbSummerAlbedoEnum,summeralbedo,P1Enum);
    4030         this->AddInput(TemperaturePDDEnum,yearlytemperatures,P1Enum); // TemperaturePDD is wrong here, but don't want to create new Enum ...
    4031        
    4032         /*clean-up*/
    4033         xDelete<IssmDouble>(monthlytemperatures);
    4034         xDelete<IssmDouble>(monthlyprecip);
    4035         xDelete<IssmDouble>(monthlylw);
    4036         xDelete<IssmDouble>(monthlysw);
    4037         xDelete<IssmDouble>(monthlywind);
    4038         xDelete<IssmDouble>(monthlyhumidity);
    4039         xDelete<IssmDouble>(smb);
    4040         xDelete<IssmDouble>(surface);
    4041         xDelete<IssmDouble>(melt);
    4042         xDelete<IssmDouble>(summermelt);
    4043         xDelete<IssmDouble>(albedo);
    4044         xDelete<IssmDouble>(summeralbedo);
    4045         xDelete<IssmDouble>(accu);
    4046         xDelete<IssmDouble>(yearlytemperatures);
    4047         xDelete<IssmDouble>(s0t);
    4048         xDelete<IssmDouble>(snowheight);
    4049         xDelete<IssmDouble>(debriscover);
    4050         xDelete<IssmDouble>(t_ampl);
    4051         xDelete<IssmDouble>(p_ampl);
    4052         xDelete<IssmDouble>(lw_ampl);
    4053         xDelete<IssmDouble>(sw_ampl);
    4054         xDelete<IssmDouble>(humidity_ampl);
    4055         xDelete<IssmDouble>(wind_ampl);
    4056         xDelete<IssmDouble>(slopex);
    4057         xDelete<IssmDouble>(slopey);
    4058         xDelete<IssmDouble>(icethickness);
    40593729}
    40603730/*}}}*/

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

    r27850 r27852  
    177177                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    178178                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
    179                 void               SmbDebrisEvatt();
    180179                void               RignotMeltParameterization();
    181180                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);
    … …  
    257256                virtual void       ComputeStressTensor(void)=0;
    258257                virtual void       ComputeEsaStrainAndVorticity(void)=0;
    259                 //virtual void       ComputeMeanEla(IssmDouble* paltitude,int* pcounter)=0;
    260258                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs* inputsin)=0;
    261259                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp)=0;

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

    r27850 r27852  
    52835283/*}}}*/
    52845284#endif
    5285 //void       Penta::ComputeMeanEla(IssmDouble* paltitude, int* pcounter){/*{{{*/
    5286 //
    5287 //      if(!this->IsOnSurface() & !this->IsIceInElement()) return;
    5288 //        int     ipos=0,ineg=0, counter=0;
    5289 //        IssmDouble mean_element_surface=0;
    5290        
    5291         /* input */
    5292 //        IssmDouble smb[NUMVERTICES],surface[NUMVERTICES];
    5293 //      Element::GetInputListOnVertices(&smb[0],SmbMassBalanceEnum);
    5294 //      Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
    5295 
    5296 //      for(int iv=0;iv<NUMVERTICES;iv++){
    5297 //              if(smb[iv]>0) ipos=ipos+1;
    5298 //              if(smb[iv]<0) ineg=ineg+1;
    5299 //      } 
    5300 
    5301         /* we define ELA if an element has pos and neg smb */
    5302 //      if(ipos>0 & ineg>0){
    5303 //              for(int iv=0;iv<NUMVERTICES;iv++) mean_element_surface=mean_element_surface+surface[iv]/double(NUMVERTICES);
    5304 //              *paltitude+=mean_element_surface;
    5305 //              *pcounter+=counter+1;
    5306 //      }
    5307 //}
    5308 /*}}}*/

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

    r27850 r27852  
    5050                void        ComputeSigmaNN(){_error_("not implemented yet");};
    5151                void        ComputeStressTensor(){_error_("not implemented yet");};
    52                 //void        ComputeMeanEla(IssmDouble* paltitude, int* pcounter){_error_("not implemented yet");};
    5352                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs* inputsin){_error_("not implemented yet");};
    5453                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp){_error_("not implemented yet");};

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

    r27850 r27852  
    4646                void        ComputeSigmaNN(){_error_("not implemented yet");};
    4747                void        ComputeStressTensor(){_error_("not implemented yet");};
    48                 //void        ComputeMeanEla(IssmDouble* paltitude, int* pcounter){_error_("not implemented yet");};
    4948                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    5049                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};

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

    r27850 r27852  
    7272                void        ComputeStressTensor();
    7373                void        ComputeSurfaceNormalVelocity();
    74                 //void        ComputeMeanEla(IssmDouble* paltitude, int* pcounter){_error_("not implemented yet");};
    7574                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs* inputsin);
    7675                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int M,int N,int interp);

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

    r27848 r27852  
    2727
    2828        /*recover parameters: */
    29         femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    3029        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    3130        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    4443                femmodel->inputs->DuplicateInput(VyEnum,VyDebrisEnum); 
    4544        }
     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        }
    4653
    4754        debris_analysis = new DebrisAnalysis();
    … …  
    5259        extrudefromtop_core(femmodel); 
    5360
    54         IssmDouble mean_ela;
    55         //ComputeMeanElax(&mean_ela,femmodel);
    56         //femmodel->parameters->SetParam(mean_ela,SmbMeanElaEnum);
    57         _printf_("core ELA: " << mean_ela << "\n");
    58        
    5961        if(save_results) femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    6062        if(solution_type==DebrisSolutionEnum)femmodel->RequestedDependentsx();

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r27847 r27852  
    489489                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismethod",SmbSemicMethodEnum));
    490490                        break;
    491                 case SMBdebrisEvattEnum:
    492                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.qlaps",SmbDesfacEnum));
    493                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
    494                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.dslaps",SmbSWlapsEnum));
    495                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.dllaps",SmbLWlapsEnum));
    496                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.windspeedlaps",SmbWindspeedlapsEnum));
    497                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.humiditylaps",SmbHumiditylapsEnum));
    498                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.icealbedo",SmbIcealbedoEnum));
    499                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.snowalbedo",SmbSnowalbedoEnum));
    500                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.debrisalbedo",SmbDebrisalbedoEnum));
    501                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismonthly",SmbDebrisIsMonthlyEnum));
    502                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.isAnderson",SmbDebrisIsAndersonEnum));
    503                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.iscryokarst",SmbDebrisIsCryokarstEnum));
    504                         parameters->AddObject(iomodel->CopyConstantObject("md.smb.AndersonD0",SmbDebrisAndersonD0Enum));
     491                case SMBdebrisMLEnum:
    505492                        break;
    506493                default:

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

    r27847 r27852  
    504504
    505505}/*}}}*/
    506 void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
    507         for(Object* & object : femmodel->elements->objects){
    508                 Element* element=xDynamicCast<Element*>(object);
    509                 element->SmbDebrisEvatt();
    510         }
     506void SmbDebrisMLx(FemModel* femmodel){/*{{{*/
     507
     508        //      The function is based on:
     509        //      Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2
     510        //      Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7
     511        //      function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90
     512
     513        /*Intermediaries*/
     514        // altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002)
     515        IssmDouble LW=2.9;          // W/m^2 /100m                       2.9
     516        IssmDouble SW=1.3;          // W/m^2 /100m                       1.3
     517        IssmDouble HumidityG=0;     // % /100m         rough estimate
     518        IssmDouble AirTemp=0.7;     // C /100m
     519        IssmDouble WindSpeed=0.02;  // m/s /100m       rough estimate    0.2
     520
     521        // accumulation follows a linear increase above the ELA up to a plateau
     522        IssmDouble AccG=0.1;                    // m w.e. /100m
     523        IssmDouble AccMax=1.;                    // m w.e.
     524        IssmDouble ReferenceElevation;
     525        IssmDouble AblationDays=120.;            //
     526
     527        IssmDouble In=100.;                 // Wm^-2        incoming long wave
     528        IssmDouble Q=500.;                  // Wm^-2        incoming short wave
     529        IssmDouble K=0.585;                // Wm^-1K^-1    thermal conductivity          0.585
     530        IssmDouble Qm=0.0012;              // kg m^-3      measured humiditiy level
     531        IssmDouble Qh=0.006 ;              // kg m^-3      saturated humidity level
     532        IssmDouble Tm=2.;                   // C            air temperature
     533        IssmDouble Rhoaa=1.22;             // kgm^-3       air densitiy
     534        IssmDouble Um=1.5;                 // ms^-1        measured wind speed
     535        IssmDouble Xm=1.5;                 // ms^-1        measurement height
     536        IssmDouble Xr=0.01;                // ms^-1        surface roughness             0.01
     537        IssmDouble Alphad=0.07;            //              debris albedo                 0.07
     538        IssmDouble Alphai=0.4;             //              ice ablbedo
     539        IssmDouble Alphaeff;
     540        IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
     541        IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
     542        IssmDouble Lm;//=3.34E+05;            // jkg^-1K^-1   latent heat of ice melt
     543        IssmDouble Lv=2.50E+06;            // jkg^-1K^-1   latent heat of evaporation
     544        IssmDouble Tf=273.;                 // K            water freeezing temperature
     545        IssmDouble Eps=0.95;               //              thermal emissivity
     546        IssmDouble Rhoi=900.;               // kgm^-3       ice density
     547        IssmDouble Sigma=5.67E-08;         // Wm^-2K^-4    Stefan Boltzmann constant
     548        IssmDouble Kstar=0.4;              //              von kármán constant
     549        IssmDouble Gamma=180.;              // m^-1         wind speed attenuation        234
     550        IssmDouble PhiD;//=0.005;              //              debris packing fraction       0.01
     551        IssmDouble Humidity=0.2;           //              relative humidity
     552
     553        IssmDouble smb,yts,z,debris;
     554        IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
     555        bool isdebris;
     556        int domaintype;
     557        femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
     558
     559        /*Get material parameters and constants */
     560        //femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as  benchmark was run with different densities for debris and FS
     561        femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum);
     562        femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     563        PhiD=0.;
     564        if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
     565
     566        /* Loop over all the elements of this partition */
     567        for(Object* & object : femmodel->elements->objects){
     568                Element* element=xDynamicCast<Element*>(object);
     569
     570                /* Allocate all arrays */
     571                int         numvertices=element->GetNumberOfVertices();
     572                IssmDouble* surfacelist=xNew<IssmDouble>(numvertices);
     573                IssmDouble* smb=xNew<IssmDouble>(numvertices);
     574                IssmDouble* debriscover=xNew<IssmDouble>(numvertices);
     575                element->GetInputListOnVertices(surfacelist,SurfaceEnum);
     576
     577                /* Get inputs */
     578                element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
     579                element->FindParam(&domaintype,DomainTypeEnum);         
     580
     581                /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
     582                for(int v=0;v<numvertices;v++){
     583
     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.)
     631                                    -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
     632                        }
     633
     634                        /* account form ablation days, and convert to m/s */
     635                        MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts;
     636
     637                        /*Update array accordingly*/
     638                        smb[v]=MassBalanceMYearDebris;
     639                }
     640
     641                /*Add input to element and Free memory*/
     642                element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
     643                xDelete<IssmDouble>(surfacelist);
     644                xDelete<IssmDouble>(smb);
     645                xDelete<IssmDouble>(debriscover);
     646        }
    511647}/*}}}*/
    512648void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/

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

    r27847 r27852  
    2323void SmbMeltComponentsx(FemModel* femmodel);
    2424void SmbGradientsComponentsx(FemModel* femmodel);
    25 void SmbDebrisEvattx(FemModel* femmodel);
     25void SmbDebrisMLx(FemModel* femmodel);
    2626/* SEMIC: */
    2727void SmbSemicx(FemModel* femmodel, int ismethod);

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

    r27849 r27852  
    8686#include "./ResetConstraintsx/ResetConstraintsx.h"
    8787#include "./ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.h"
    88 //#include "./ComputeMeanElax/ComputeMeanElax.h"
    8988#include "./RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
    9089#include "./RheologyBAbsGradientx/RheologyBAbsGradientx.h"

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

    r27847 r27852  
    176176syn keyword cConstant DebrisRemovalStressThresholdEnum
    177177syn keyword cConstant DebrisPackingFractionEnum
    178 syn keyword cConstant DebrisMaxdisplacementvelocityEnum
    179178syn keyword cConstant DebugProfilingEnum
    180179syn keyword cConstant DomainDimensionEnum
    … …  
    540539syn keyword cConstant SmbARMAmaOrderEnum
    541540syn keyword cConstant SmbAveragingEnum
    542 syn keyword cConstant SmbDebrisalbedoEnum
    543 syn keyword cConstant SmbDebrisIsMonthlyEnum
    544 syn keyword cConstant SmbDebrisIsAndersonEnum
    545 syn keyword cConstant SmbDebrisIsCryokarstEnum
    546 syn keyword cConstant SmbDebrisAndersonD0Enum
    547541syn keyword cConstant SmbDesfacEnum
    548542syn keyword cConstant SmbDesfacElevEnum
    … …  
    558552syn keyword cConstant SmbEIdxEnum
    559553syn keyword cConstant SmbFEnum
    560 syn keyword cConstant SmbHumiditylapsEnum
    561554syn keyword cConstant SmbInitDensityScalingEnum
    562 syn keyword cConstant SmbIcealbedoEnum
    563555syn keyword cConstant SmbIsaccumulationEnum
    564556syn keyword cConstant SmbIsalbedoEnum
    … …  
    580572syn keyword cConstant SmbKEnum
    581573syn keyword cConstant SmbLapseRatesEnum
    582 syn keyword cConstant SmbLWlapsEnum
    583 syn keyword cConstant SmbSWlapsEnum
    584 syn keyword cConstant SmbSnowalbedoEnum
    585 syn keyword cConstant SmbMeanElaEnum
    586574syn keyword cConstant SmbNumBasinsEnum
    587575syn keyword cConstant SmbNumBreaksEnum
    … …  
    624612syn keyword cConstant SmbTemperaturesReconstructedYearsEnum
    625613syn keyword cConstant SmbPrecipitationsReconstructedYearsEnum
    626 syn keyword cConstant SmbWindspeedlapsEnum
    627614syn keyword cConstant SmoothThicknessMultiplierEnum
    628615syn keyword cConstant SolutionTypeEnum
    … …  
    792779syn keyword cConstant DamageFEnum
    793780syn keyword cConstant DebrisThicknessEnum
    794 syn keyword cConstant DebrisThicknessOldEnum
    795781syn keyword cConstant DegreeOfChannelizationEnum
    796782syn keyword cConstant DepthBelowSurfaceEnum
    … …  
    11211107syn keyword cConstant SmbMeltEnum
    11221108syn keyword cConstant SmbMonthlytemperaturesEnum
    1123 syn keyword cConstant SmbMonthlydsradiationEnum
    1124 syn keyword cConstant SmbMonthlydlradiationEnum
    1125 syn keyword cConstant SmbMonthlywindspeedEnum
    1126 syn keyword cConstant SmbMonthlyairhumidityEnum
    11271109syn keyword cConstant SmbMSurfEnum
    11281110syn keyword cConstant SmbNetLWEnum
    … …  
    11341116syn keyword cConstant SmbPrecipitationEnum
    11351117syn keyword cConstant SmbPrecipitationsAnomalyEnum
    1136 syn keyword cConstant SmbDsradiationAnomalyEnum
    1137 syn keyword cConstant SmbDlradiationAnomalyEnum
    1138 syn keyword cConstant SmbWindspeedAnomalyEnum
    1139 syn keyword cConstant SmbAirhumidityAnomalyEnum
    11401118syn keyword cConstant SmbPrecipitationsLgmEnum
    11411119syn keyword cConstant SmbPrecipitationsPresentdayEnum
    … …  
    11571135syn keyword cConstant SmbSmbrefEnum
    11581136syn keyword cConstant SmbSzaValueEnum
    1159 syn keyword cConstant SmbSummerMeltEnum
    1160 syn keyword cConstant SmbSummerAlbedoEnum
    1161 syn keyword cConstant SmbSnowheightEnum
    11621137syn keyword cConstant SmbTEnum
    11631138syn keyword cConstant SmbTaEnum
    … …  
    16791654syn keyword cConstant SMBarmaEnum
    16801655syn keyword cConstant SMBcomponentsEnum
    1681 syn keyword cConstant SMBdebrisEvattEnum
     1656syn keyword cConstant SMBdebrisMLEnum
    16821657syn keyword cConstant SMBd18opddEnum
    16831658syn keyword cConstant SMBforcingEnum
    … …  
    18051780syn keyword cType Cfsurfacesquaretransient
    18061781syn keyword cType Channel
    1807 syn keyword cType classes
    18081782syn keyword cType Constraint
    18091783syn keyword cType Constraints
    … …  
    18131787syn keyword cType ControlParam
    18141788syn keyword cType Covertree
     1789syn keyword cType DataSetParam
    18151790syn keyword cType DatasetInput
    1816 syn keyword cType DataSetParam
    18171791syn keyword cType Definition
    18181792syn keyword cType DependentObject
    … …  
    18271801syn keyword cType ElementInput
    18281802syn keyword cType ElementMatrix
     1803syn keyword cType ElementVector
    18291804syn keyword cType Elements
    1830 syn keyword cType ElementVector
    18311805syn keyword cType ExponentialVariogram
    18321806syn keyword cType ExternalResult
    … …  
    18351809syn keyword cType Friction
    18361810syn keyword cType Gauss
    1837 syn keyword cType GaussianVariogram
    1838 syn keyword cType gaussobjects
    18391811syn keyword cType GaussPenta
    18401812syn keyword cType GaussSeg
    18411813syn keyword cType GaussTetra
    18421814syn keyword cType GaussTria
     1815syn keyword cType GaussianVariogram
    18431816syn keyword cType GenericExternalResult
    18441817syn keyword cType GenericOption
    … …  
    18571830syn keyword cType IssmDirectApplicInterface
    18581831syn keyword cType IssmParallelDirectApplicInterface
    1859 syn keyword cType krigingobjects
    18601832syn keyword cType Load
    18611833syn keyword cType Loads
    … …  
    18681840syn keyword cType Matice
    18691841syn keyword cType Matlitho
    1870 syn keyword cType matrixobjects
    18711842syn keyword cType MatrixParam
    18721843syn keyword cType Misfit
    … …  
    18811852syn keyword cType Observations
    18821853syn keyword cType Option
     1854syn keyword cType OptionUtilities
    18831855syn keyword cType Options
    1884 syn keyword cType OptionUtilities
    18851856syn keyword cType Param
    18861857syn keyword cType Parameters
    … …  
    18961867syn keyword cType Regionaloutput
    18971868syn keyword cType Results
     1869syn keyword cType RiftStruct
    18981870syn keyword cType Riftfront
    1899 syn keyword cType RiftStruct
    19001871syn keyword cType SealevelGeometry
    19011872syn keyword cType Seg
    19021873syn keyword cType SegInput
     1874syn keyword cType SegRef
    19031875syn keyword cType Segment
    1904 syn keyword cType SegRef
    19051876syn keyword cType SpcDynamic
    19061877syn keyword cType SpcStatic
    … …  
    19211892syn keyword cType Vertex
    19221893syn keyword cType Vertices
     1894syn keyword cType classes
     1895syn keyword cType gaussobjects
     1896syn keyword cType krigingobjects
     1897syn keyword cType matrixobjects
    19231898syn keyword cType AdjointBalancethickness2Analysis
    19241899syn keyword cType AdjointBalancethicknessAnalysis

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

    r27847 r27852  
    170170        DebrisRemovalStressThresholdEnum,
    171171        DebrisPackingFractionEnum,
    172         DebrisMaxdisplacementvelocityEnum,
    173172        DebugProfilingEnum,
    174173        DomainDimensionEnum,
    … …  
    531530        SmbAlbedoLandEnum,
    532531        SmbARMATimestepEnum,
    533         SmbARMAarOrderEnum,
    534         SmbARMAmaOrderEnum,
     532   SmbARMAarOrderEnum,
     533   SmbARMAmaOrderEnum,
    535534        SmbAveragingEnum,
    536         SmbDebrisalbedoEnum,
    537         SmbDebrisIsMonthlyEnum,
    538         SmbDebrisIsAndersonEnum,
    539         SmbDebrisIsCryokarstEnum,
    540         SmbDebrisAndersonD0Enum,
    541535        SmbDesfacEnum,
    542536        SmbDesfacElevEnum,
    543537        SmbDpermilEnum,
    544538        SmbDsnowIdxEnum,
    545         SmbElevationBinsEnum,
     539   SmbElevationBinsEnum,
    546540        SmbCldFracEnum,
    547541        SmbDelta18oEnum,
    … …  
    552546        SmbEIdxEnum,
    553547        SmbFEnum,
    554         SmbHumiditylapsEnum,
    555548        SmbInitDensityScalingEnum,
    556         SmbIcealbedoEnum,
    557549        SmbIsaccumulationEnum,
    558550        SmbIsalbedoEnum,
    … …  
    573565        SmbIsturbulentfluxEnum,
    574566        SmbKEnum,
    575         SmbLapseRatesEnum,
    576         SmbLWlapsEnum,
    577         SmbSWlapsEnum,
    578         SmbSnowalbedoEnum,
    579         SmbMeanElaEnum,
     567   SmbLapseRatesEnum,
    580568        SmbNumBasinsEnum,
    581569        SmbNumBreaksEnum,
    … …  
    618606        SmbTemperaturesReconstructedYearsEnum,
    619607        SmbPrecipitationsReconstructedYearsEnum,
    620         SmbWindspeedlapsEnum,
    621608        SmoothThicknessMultiplierEnum,
    622609        SolutionTypeEnum,
    … …  
    788775        DamageFEnum,
    789776        DebrisThicknessEnum,
    790         DebrisThicknessOldEnum,
    791777        DegreeOfChannelizationEnum,
    792778        DepthBelowSurfaceEnum,
    … …  
    10571043        SmbAdiffiniEnum,
    10581044        SmbAiniEnum,
    1059         SmbARMANoiseEnum,
     1045   SmbARMANoiseEnum,
    10601046        SmbBasinsIdEnum,
    10611047        SmbBMaxEnum,
    … …  
    11181104        SmbMInitnum,
    11191105        SmbMonthlytemperaturesEnum,
    1120         SmbMonthlydsradiationEnum,
    1121         SmbMonthlydlradiationEnum,
    1122         SmbMonthlywindspeedEnum,
    1123         SmbMonthlyairhumidityEnum,
    11241106        SmbMSurfEnum,
    11251107        SmbNetLWEnum,
    … …  
    11311113        SmbPrecipitationEnum,
    11321114        SmbPrecipitationsAnomalyEnum,
    1133         SmbDsradiationAnomalyEnum,
    1134         SmbDlradiationAnomalyEnum,
    1135         SmbWindspeedAnomalyEnum,
    1136         SmbAirhumidityAnomalyEnum,
    11371115        SmbPrecipitationsLgmEnum,
    11381116        SmbPrecipitationsPresentdayEnum,
    … …  
    11541132        SmbSmbrefEnum,
    11551133        SmbSzaValueEnum,
    1156         SmbSummerMeltEnum,
    1157         SmbSummerAlbedoEnum,
    1158         SmbSnowheightEnum,
    11591134        SmbTEnum,
    11601135        SmbTaEnum,
    … …  
    16781653        SMBarmaEnum,
    16791654        SMBcomponentsEnum,
    1680         SMBdebrisEvattEnum,
     1655        SMBdebrisMLEnum,
    16811656        SMBd18opddEnum,
    16821657        SMBforcingEnum,

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

    r27847 r27852  
    178178                case DebrisRemovalStressThresholdEnum : return "DebrisRemovalStressThreshold";
    179179                case DebrisPackingFractionEnum : return "DebrisPackingFraction";
    180                 case DebrisMaxdisplacementvelocityEnum : return "DebrisMaxdisplacementvelocity";
    181180                case DebugProfilingEnum : return "DebugProfiling";
    182181                case DomainDimensionEnum : return "DomainDimension";
    … …  
    542541                case SmbARMAmaOrderEnum : return "SmbARMAmaOrder";
    543542                case SmbAveragingEnum : return "SmbAveraging";
    544                 case SmbDebrisalbedoEnum : return "SmbDebrisalbedo";
    545                 case SmbDebrisIsMonthlyEnum : return "SmbDebrisIsMonthly";
    546                 case SmbDebrisIsAndersonEnum : return "SmbDebrisIsAnderson";
    547                 case SmbDebrisIsCryokarstEnum : return "SmbDebrisIsCryokarst";
    548                 case SmbDebrisAndersonD0Enum : return "SmbDebrisAndersonD0";
    549543                case SmbDesfacEnum : return "SmbDesfac";
    550544                case SmbDesfacElevEnum : return "SmbDesfacElev";
    … …  
    560554                case SmbEIdxEnum : return "SmbEIdx";
    561555                case SmbFEnum : return "SmbF";
    562                 case SmbHumiditylapsEnum : return "SmbHumiditylaps";
    563556                case SmbInitDensityScalingEnum : return "SmbInitDensityScaling";
    564                 case SmbIcealbedoEnum : return "SmbIcealbedo";
    565557                case SmbIsaccumulationEnum : return "SmbIsaccumulation";
    566558                case SmbIsalbedoEnum : return "SmbIsalbedo";
    … …  
    582574                case SmbKEnum : return "SmbK";
    583575                case SmbLapseRatesEnum : return "SmbLapseRates";
    584                 case SmbLWlapsEnum : return "SmbLWlaps";
    585                 case SmbSWlapsEnum : return "SmbSWlaps";
    586                 case SmbSnowalbedoEnum : return "SmbSnowalbedo";
    587                 case SmbMeanElaEnum : return "SmbMeanEla";
    588576                case SmbNumBasinsEnum : return "SmbNumBasins";
    589577                case SmbNumBreaksEnum : return "SmbNumBreaks";
    … …  
    626614                case SmbTemperaturesReconstructedYearsEnum : return "SmbTemperaturesReconstructedYears";
    627615                case SmbPrecipitationsReconstructedYearsEnum : return "SmbPrecipitationsReconstructedYears";
    628                 case SmbWindspeedlapsEnum : return "SmbWindspeedlaps";
    629616                case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier";
    630617                case SolutionTypeEnum : return "SolutionType";
    … …  
    794781                case DamageFEnum : return "DamageF";
    795782                case DebrisThicknessEnum : return "DebrisThickness";
    796                 case DebrisThicknessOldEnum : return "DebrisThicknessOld";
    797783                case DegreeOfChannelizationEnum : return "DegreeOfChannelization";
    798784                case DepthBelowSurfaceEnum : return "DepthBelowSurface";
    … …  
    11231109                case SmbMeltEnum : return "SmbMelt";
    11241110                case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures";
    1125                 case SmbMonthlydsradiationEnum : return "SmbMonthlydsradiation";
    1126                 case SmbMonthlydlradiationEnum : return "SmbMonthlydlradiation";
    1127                 case SmbMonthlywindspeedEnum : return "SmbMonthlywindspeed";
    1128                 case SmbMonthlyairhumidityEnum : return "SmbMonthlyairhumidity";
    11291111                case SmbMSurfEnum : return "SmbMSurf";
    11301112                case SmbNetLWEnum : return "SmbNetLW";
    … …  
    11361118                case SmbPrecipitationEnum : return "SmbPrecipitation";
    11371119                case SmbPrecipitationsAnomalyEnum : return "SmbPrecipitationsAnomaly";
    1138                 case SmbDsradiationAnomalyEnum : return "SmbDsradiationAnomaly";
    1139                 case SmbDlradiationAnomalyEnum : return "SmbDlradiationAnomaly";
    1140                 case SmbWindspeedAnomalyEnum : return "SmbWindspeedAnomaly";
    1141                 case SmbAirhumidityAnomalyEnum : return "SmbAirhumidityAnomaly";
    11421120                case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm";
    11431121                case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday";
    … …  
    11591137                case SmbSmbrefEnum : return "SmbSmbref";
    11601138                case SmbSzaValueEnum : return "SmbSzaValue";
    1161                 case SmbSummerMeltEnum : return "SmbSummerMelt";
    1162                 case SmbSummerAlbedoEnum : return "SmbSummerAlbedo";
    1163                 case SmbSnowheightEnum : return "SmbSnowheight";
    11641139                case SmbTEnum : return "SmbT";
    11651140                case SmbTaEnum : return "SmbTa";
    … …  
    16811656                case SMBarmaEnum : return "SMBarma";
    16821657                case SMBcomponentsEnum : return "SMBcomponents";
    1683                 case SMBdebrisEvattEnum : return "SMBdebrisEvatt";
     1658                case SMBdebrisMLEnum : return "SMBdebrisML";
    16841659                case SMBd18opddEnum : return "SMBd18opdd";
    16851660                case SMBforcingEnum : return "SMBforcing";

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

    r27847 r27852  
    169169syn keyword juliaConstC DebrisRemovalStressThresholdEnum
    170170syn keyword juliaConstC DebrisPackingFractionEnum
    171 syn keyword juliaConstC DebrisMaxdisplacementvelocityEnum
    172171syn keyword juliaConstC DebugProfilingEnum
    173172syn keyword juliaConstC DomainDimensionEnum
    … …  
    533532syn keyword juliaConstC SmbARMAmaOrderEnum
    534533syn keyword juliaConstC SmbAveragingEnum
    535 syn keyword juliaConstC SmbDebrisalbedoEnum
    536 syn keyword juliaConstC SmbDebrisIsMonthlyEnum
    537 syn keyword juliaConstC SmbDebrisIsAndersonEnum
    538 syn keyword juliaConstC SmbDebrisIsCryokarstEnum
    539 syn keyword juliaConstC SmbDebrisAndersonD0Enum
    540534syn keyword juliaConstC SmbDesfacEnum
    541535syn keyword juliaConstC SmbDesfacElevEnum
    … …  
    551545syn keyword juliaConstC SmbEIdxEnum
    552546syn keyword juliaConstC SmbFEnum
    553 syn keyword juliaConstC SmbHumiditylapsEnum
    554547syn keyword juliaConstC SmbInitDensityScalingEnum
    555 syn keyword juliaConstC SmbIcealbedoEnum
    556548syn keyword juliaConstC SmbIsaccumulationEnum
    557549syn keyword juliaConstC SmbIsalbedoEnum
    … …  
    573565syn keyword juliaConstC SmbKEnum
    574566syn keyword juliaConstC SmbLapseRatesEnum
    575 syn keyword juliaConstC SmbLWlapsEnum
    576 syn keyword juliaConstC SmbSWlapsEnum
    577 syn keyword juliaConstC SmbSnowalbedoEnum
    578 syn keyword juliaConstC SmbMeanElaEnum
    579567syn keyword juliaConstC SmbNumBasinsEnum
    580568syn keyword juliaConstC SmbNumBreaksEnum
    … …  
    617605syn keyword juliaConstC SmbTemperaturesReconstructedYearsEnum
    618606syn keyword juliaConstC SmbPrecipitationsReconstructedYearsEnum
    619 syn keyword juliaConstC SmbWindspeedlapsEnum
    620607syn keyword juliaConstC SmoothThicknessMultiplierEnum
    621608syn keyword juliaConstC SolutionTypeEnum
    … …  
    785772syn keyword juliaConstC DamageFEnum
    786773syn keyword juliaConstC DebrisThicknessEnum
    787 syn keyword juliaConstC DebrisThicknessOldEnum
    788774syn keyword juliaConstC DegreeOfChannelizationEnum
    789775syn keyword juliaConstC DepthBelowSurfaceEnum
    … …  
    11141100syn keyword juliaConstC SmbMeltEnum
    11151101syn keyword juliaConstC SmbMonthlytemperaturesEnum
    1116 syn keyword juliaConstC SmbMonthlydsradiationEnum
    1117 syn keyword juliaConstC SmbMonthlydlradiationEnum
    1118 syn keyword juliaConstC SmbMonthlywindspeedEnum
    1119 syn keyword juliaConstC SmbMonthlyairhumidityEnum
    11201102syn keyword juliaConstC SmbMSurfEnum
    11211103syn keyword juliaConstC SmbNetLWEnum
    … …  
    11271109syn keyword juliaConstC SmbPrecipitationEnum
    11281110syn keyword juliaConstC SmbPrecipitationsAnomalyEnum
    1129 syn keyword juliaConstC SmbDsradiationAnomalyEnum
    1130 syn keyword juliaConstC SmbDlradiationAnomalyEnum
    1131 syn keyword juliaConstC SmbWindspeedAnomalyEnum
    1132 syn keyword juliaConstC SmbAirhumidityAnomalyEnum
    11331111syn keyword juliaConstC SmbPrecipitationsLgmEnum
    11341112syn keyword juliaConstC SmbPrecipitationsPresentdayEnum
    … …  
    11501128syn keyword juliaConstC SmbSmbrefEnum
    11511129syn keyword juliaConstC SmbSzaValueEnum
    1152 syn keyword juliaConstC SmbSummerMeltEnum
    1153 syn keyword juliaConstC SmbSummerAlbedoEnum
    1154 syn keyword juliaConstC SmbSnowheightEnum
    11551130syn keyword juliaConstC SmbTEnum
    11561131syn keyword juliaConstC SmbTaEnum
    … …  
    16721647syn keyword juliaConstC SMBarmaEnum
    16731648syn keyword juliaConstC SMBcomponentsEnum
    1674 syn keyword juliaConstC SMBdebrisEvattEnum
     1649syn keyword juliaConstC SMBdebrisMLEnum
    16751650syn keyword juliaConstC SMBd18opddEnum
    16761651syn keyword juliaConstC SMBforcingEnum

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

    r27847 r27852  
    181181              else if (strcmp(name,"DebrisRemovalStressThreshold")==0) return DebrisRemovalStressThresholdEnum;
    182182              else if (strcmp(name,"DebrisPackingFraction")==0) return DebrisPackingFractionEnum;
    183               else if (strcmp(name,"DebrisMaxdisplacementvelocity")==0) return DebrisMaxdisplacementvelocityEnum;
    184183              else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
    185184              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    … …  
    260259              else if (strcmp(name,"Hydrologyarmadatebreaks")==0) return HydrologyarmadatebreaksEnum;
    261260              else if (strcmp(name,"Hydrologyarmamalagcoefs")==0) return HydrologyarmamalagcoefsEnum;
     261              else if (strcmp(name,"HydrologyarmamaOrder")==0) return HydrologyarmamaOrderEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"HydrologyarmamaOrder")==0) return HydrologyarmamaOrderEnum;
    266               else if (strcmp(name,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum;
     265              if (strcmp(name,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum;
    267266              else if (strcmp(name,"HydrologyarmaNumBreaks")==0) return HydrologyarmaNumBreaksEnum;
    268267              else if (strcmp(name,"HydrologyarmaNumParams")==0) return HydrologyarmaNumParamsEnum;
    … …  
    383382              else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
    384383              else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum;
     384              else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    389               else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
     388              if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    390389              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    391390              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    … …  
    506505              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    507506              else if (strcmp(name,"SolidearthSettingsTimeAcc")==0) return SolidearthSettingsTimeAccEnum;
     507              else if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
    512               else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
     511              if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
    513512              else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
    514513              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
    … …  
    554553              else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum;
    555554              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    556               else if (strcmp(name,"SmbDebrisalbedo")==0) return SmbDebrisalbedoEnum;
    557               else if (strcmp(name,"SmbDebrisIsMonthly")==0) return SmbDebrisIsMonthlyEnum;
    558               else if (strcmp(name,"SmbDebrisIsAnderson")==0) return SmbDebrisIsAndersonEnum;
    559               else if (strcmp(name,"SmbDebrisIsCryokarst")==0) return SmbDebrisIsCryokarstEnum;
    560               else if (strcmp(name,"SmbDebrisAndersonD0")==0) return SmbDebrisAndersonD0Enum;
    561555              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    562556              else if (strcmp(name,"SmbDesfacElev")==0) return SmbDesfacElevEnum;
    … …  
    572566              else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum;
    573567              else if (strcmp(name,"SmbF")==0) return SmbFEnum;
    574               else if (strcmp(name,"SmbHumiditylaps")==0) return SmbHumiditylapsEnum;
    575568              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
    576               else if (strcmp(name,"SmbIcealbedo")==0) return SmbIcealbedoEnum;
    577569              else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum;
    578570              else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
    … …  
    594586              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    595587              else if (strcmp(name,"SmbLapseRates")==0) return SmbLapseRatesEnum;
    596               else if (strcmp(name,"SmbLWlaps")==0) return SmbLWlapsEnum;
    597               else if (strcmp(name,"SmbSWlaps")==0) return SmbSWlapsEnum;
    598               else if (strcmp(name,"SmbSnowalbedo")==0) return SmbSnowalbedoEnum;
    599               else if (strcmp(name,"SmbMeanEla")==0) return SmbMeanElaEnum;
    600588              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    601589              else if (strcmp(name,"SmbNumBreaks")==0) return SmbNumBreaksEnum;
    … …  
    629617              else if (strcmp(name,"SmbSemicTmax")==0) return SmbSemicTmaxEnum;
    630618              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
     619              else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
    635620              else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
    636621              else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
    … …  
    641626              else if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum;
    642627              else if (strcmp(name,"SmbPrecipitationsReconstructedYears")==0) return SmbPrecipitationsReconstructedYearsEnum;
    643               else if (strcmp(name,"SmbWindspeedlaps")==0) return SmbWindspeedlapsEnum;
    644628              else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
    645629              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    646630              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
    647               else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    648635              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    649636              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    … …  
    752739              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
    753740              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
     741              else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
    758742              else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
    759743              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    … …  
    768752              else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
    769753              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    770               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    771758              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    772759              else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
    … …  
    812799              else if (strcmp(name,"DamageF")==0) return DamageFEnum;
    813800              else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum;
    814               else if (strcmp(name,"DebrisThicknessOld")==0) return DebrisThicknessOldEnum;
    815801              else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
    816802              else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
    … …  
    875861              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    876862              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
     863              else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    881864              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    882865              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    … …  
    892875              else if (strcmp(name,"UGia")==0) return UGiaEnum;
    893876              else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum;
    894               else if (strcmp(name,"Gradient")==0) return GradientEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Gradient")==0) return GradientEnum;
    895881              else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum;
    896882              else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum;
    … …  
    998984              else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    999985              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
     986              else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    1004987              else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum;
    1005988              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    … …  
    1015998              else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
    1016999              else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
    1017               else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
    10181004              else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
    10191005              else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
    … …  
    11211107              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    11221108              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
     1109              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    11271110              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    11281111              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    … …  
    11381121              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    11391122              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    1140               else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    11411127              else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum;
    11421128              else if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum;
    … …  
    11501136              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    11511137              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
    1152               else if (strcmp(name,"SmbMonthlydsradiation")==0) return SmbMonthlydsradiationEnum;
    1153               else if (strcmp(name,"SmbMonthlydlradiation")==0) return SmbMonthlydlradiationEnum;
    1154               else if (strcmp(name,"SmbMonthlywindspeed")==0) return SmbMonthlywindspeedEnum;
    1155               else if (strcmp(name,"SmbMonthlyairhumidity")==0) return SmbMonthlyairhumidityEnum;
    11561138              else if (strcmp(name,"SmbMSurf")==0) return SmbMSurfEnum;
    11571139              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
    … …  
    11631145              else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
    11641146              else if (strcmp(name,"SmbPrecipitationsAnomaly")==0) return SmbPrecipitationsAnomalyEnum;
    1165               else if (strcmp(name,"SmbDsradiationAnomaly")==0) return SmbDsradiationAnomalyEnum;
    1166               else if (strcmp(name,"SmbDlradiationAnomaly")==0) return SmbDlradiationAnomalyEnum;
    1167               else if (strcmp(name,"SmbWindspeedAnomaly")==0) return SmbWindspeedAnomalyEnum;
    1168               else if (strcmp(name,"SmbAirhumidityAnomaly")==0) return SmbAirhumidityAnomalyEnum;
    11691147              else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
    11701148              else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
    … …  
    11861164              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    11871165              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    1188               else if (strcmp(name,"SmbSummerMelt")==0) return SmbSummerMeltEnum;
    1189               else if (strcmp(name,"SmbSummerAlbedo")==0) return SmbSummerAlbedoEnum;
    1190               else if (strcmp(name,"SmbSnowheight")==0) return SmbSnowheightEnum;
    11911166              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    11921167              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    … …  
    12441219              else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
    12451220              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
     1221              else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    12501222              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    12511223              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    … …  
    12721244              else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
    12731245              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    1274               else if (strcmp(name,"Vel")==0) return VelEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Vel")==0) return VelEnum;
    12751250              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    12761251              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    … …  
    13671342              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    13681343              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
     1344              else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    13731345              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    13741346              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    … …  
    13951367              else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum;
    13961368              else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum;
    1397               else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
    13981373              else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
    13991374              else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
    … …  
    14901465              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    14911466              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"Dense")==0) return DenseEnum;
     1467              else if (strcmp(name,"Dense")==0) return DenseEnum;
    14961468              else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
    14971469              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    … …  
    15181490              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    15191491              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    1520               else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    15211496              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    15221497              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    … …  
    16131588              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    16141589              else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    1615          else stage=14;
    1616    }
    1617    if(stage==14){
    1618               if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
     1590              else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
    16191591              else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
    16201592              else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
    … …  
    16411613              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    16421614              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    1643               else if (strcmp(name,"Matice")==0) return MaticeEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"Matice")==0) return MaticeEnum;
    16441619              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    16451620              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    … …  
    17201695              else if (strcmp(name,"SMBarma")==0) return SMBarmaEnum;
    17211696              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    1722               else if (strcmp(name,"SMBdebrisEvatt")==0) return SMBdebrisEvattEnum;
     1697              else if (strcmp(name,"SMBdebrisML")==0) return SMBdebrisMLEnum;
    17231698              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    17241699              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    … …  
    17361711              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    17371712              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    1738          else stage=15;
    1739    }
    1740    if(stage==15){
    1741               if (strcmp(name,"Scaled")==0) return ScaledEnum;
     1713              else if (strcmp(name,"Scaled")==0) return ScaledEnum;
    17421714              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    17431715              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    … …  
    17641736              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    17651737              else if (strcmp(name,"Sset")==0) return SsetEnum;
    1766               else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
     1738         else stage=15;
     1739   }
     1740   if(stage==15){
     1741              if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
    17671742              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    17681743              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r27847 r27852  
    262262                case 12: return SMBsemicEnum;   
    263263                case 13: return SMBarmaEnum;
    264                 case 14: return SMBdebrisEvattEnum;
     264                case 14: return SMBdebrisMLEnum;
    265265                default: _error_("Marshalled SMB code \""<<enum_in<<"\" not supported yet");
    266266        }

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

    r27847 r27852  
    1212                removalmodel             = 0;
    1313                displacementmodel        = 0;
    14                 max_displacementvelocity = 0;
    1514                removal_slope_threshold  = 0;
    1615                removal_stress_threshold = 0;
    … …  
    6766                function list = defaultoutputs(self,md) % {{{
    6867
    69                         list = {'DebrisThickness','DebrisMaskNodeActivation','VxDebris','VyDebris'};
     68                        list = {'DebrisThickness'};
    7069
    7170                end % }}}
    … …  
    9392                        self.removal_stress_threshold=0;
    9493
    95                         %Max velocity for displacementmodel (1)
    96                         self.max_displacementvelocity=0;
    97 
    9894                        %default output
    9995                        self.requested_outputs={'default'};
    … …  
    105101
    106102                        md = checkfield(md,'fieldname','debris.spcthickness');
    107                         md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3 4 5]);
     103                        md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3]);
    108104                        md = checkfield(md,'fieldname','debris.min_thickness','>=',0);
    109105                        md = checkfield(md,'fieldname','debris.packingfraction','>=',0);
    … …  
    112108                        md = checkfield(md,'fieldname','debris.removal_slope_threshold','>=',0);
    113109                        md = checkfield(md,'fieldname','debris.removal_stress_threshold','>=',0);
    114                         md = checkfield(md,'fieldname','debris.max_displacementvelocity','>=',0);
    115110                        md = checkfield(md,'fieldname','debris.requested_outputs','stringrow',1);
    116111                        if ~any(isnan(md.stressbalance.vertex_pairing)),
    … …  
    125120                        fielddisplay(self,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: streamline upwind Petrov-Galerkin (SUPG)');
    126121                        fielddisplay(self,'removalmodel','frontal removal of debris. 0: no removal, 1: Slope-triggered debris removal, 2: driving-stress triggered debris removal');
    127                         fielddisplay(self,'max_displacementvelocity','maximum velocity of debris transport (v_ice + v_displacement) (m/a)');
    128                         fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: additional debris velocity above the critical slope/stress threshold');
     122                        fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: ...');
    129123                        fielddisplay(self,'removal_slope_threshold','critical slope (degrees) for removalmodel (1)');
    130124                        fielddisplay(self,'removal_stress_threshold','critical stress (Pa) for removalmodel (2)');
    … …  
    141135                        WriteData(fid,prefix,'object',self,'fieldname','removalmodel','format','Integer');
    142136                        WriteData(fid,prefix,'object',self,'fieldname','displacementmodel','format','Integer');
    143                         WriteData(fid,prefix,'object',self,'fieldname','max_displacementvelocity','format','Double');
    144137                        WriteData(fid,prefix,'object',self,'fieldname','removal_slope_threshold','format','Double');
    145138                        WriteData(fid,prefix,'object',self,'fieldname','removal_stress_threshold','format','Double');
    … …  
    163156                        writejsdouble(fid,[modelname '.debris.removalmodel'],self.removalmodel);
    164157                        writejsdouble(fid,[modelname '.debris.displacementmodel'],self.displacementmodel);
    165                         writejsdouble(fid,[modelname '.debris.max_displacementvelocity'],self.displacementmodel);
    166158                        writejsdouble(fid,[modelname '.debris.removal_slope_threshold'],self.removal_slope_threshold);
    167159                        writejsdouble(fid,[modelname '.debris.removal_stress_threshold'],self.removal_stress_threshold);

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

    r27847 r27852  
    12891289                        if isfield(structmd,'cfl_coefficient'), md.timestepping.cfl_coefficient=structmd.cfl_coefficient; end
    12901290                        if isfield(structmd,'spcthickness'), md.masstransport.spcthickness=structmd.spcthickness; end
    1291                         if isfield(structmd,'spcthickness'), md.debris.spcthickness=structmd.spcthickness; end
    12921291                        if isfield(structmd,'artificial_diffusivity'), md.masstransport.stabilization=structmd.artificial_diffusivity; end
    12931292                        if isfield(structmd,'hydrostatic_adjustment'), md.masstransport.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end

  • issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m

    r27847 r27852  
    240240                field = field*yts;
    241241        elseif strcmp(fieldname,'VyAverage'),
    242                 field = field*yts;
    243         elseif strcmp(fieldname,'VxDebris'),
    244                 field = field*yts;
    245         elseif strcmp(fieldname,'VyDebris'),
    246242                field = field*yts;
    247243        elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),
Note: See TracChangeset for help on using the changeset viewer.