Changeset 27847


Ignore:
Timestamp:
07/24/23 05:33:26 (20 months ago)
Author:
rueckamp
Message:

CHG: Several imprvements to Debris analysis

Location:
issm/trunk-jpl/src
Files:
1 added
23 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r27735 r27847  
    253253        ./modules/StochasticForcingx/StochasticForcingx.cpp \
    254254        ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp \
     255        ./modules/ComputeMeanElax/ComputeMeanElax.cpp \
    255256        ./cores/ProcessArguments.cpp \
    256257        ./cores/ResetBoundaryConditions.cpp \

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

    r27795 r27847  
    1010
    1111#define FINITEELEMENT P1Enum
    12 #define EPS 1e-14
     12#define EPS 1e-8
     13#define RHO 2
     14//#define KAPPA 5
    1315
    1416/*Model processing*/
    … …  
    115117        parameters->AddObject(iomodel->CopyConstantObject("md.debris.min_thickness",DebrisMinThicknessEnum));
    116118        parameters->AddObject(iomodel->CopyConstantObject("md.debris.packingfraction",DebrisPackingFractionEnum));
     119        parameters->AddObject(iomodel->CopyConstantObject("md.debris.max_displacementvelocity",DebrisMaxdisplacementvelocityEnum));
    117120        parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_slope_threshold",DebrisRemovalSlopeThresholdEnum));
    118121        parameters->AddObject(iomodel->CopyConstantObject("md.debris.removal_stress_threshold",DebrisRemovalStressThresholdEnum));
    … …  
    128131void           DebrisAnalysis::Core(FemModel* femmodel){/*{{{*/
    129132
    130         //PreProcessing(femmodel);
    131         //femmodel->parameters->SetParam(VxDebrisEnum,InputToExtrudeEnum);
    132         //extrudefromtop_core(femmodel);
     133        PreProcessing(femmodel);
    133134        femmodel->SetCurrentConfiguration(DebrisAnalysisEnum);       
    134135        solutionsequence_linear(femmodel);
     136        femmodel->inputs->DuplicateInput(DebrisThicknessEnum,DebrisThicknessOldEnum);
    135137        PostProcessing(femmodel);
    136 
     138       
    137139}/*}}}*/
    138140void           DebrisAnalysis::PreCore(FemModel* femmodel){/*{{{*/
    … …  
    197199        Input* vy_input=NULL;
    198200        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);
    199204        h=topelement->CharacteristicLength();
     205        //IssmDouble hx,hy,hz;
     206        //topelement->ElementSizes(&hx,&hy,&hz);
     207        //h=sqrt(hx*hx+hy*hy);
    200208
    201209        /* Start  looping on the number of gaussian points: */
    … …  
    209217                vx_input->GetInputValue(&vx,gauss);
    210218                if(dim==2) vy_input->GetInputValue(&vy,gauss);
    211 
     219       
    212220                D_scalar=gauss->weight*Jdet;
    213221                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];
    … …  
    215223                /*Advection terms*/
    216224                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
     225                dvxdx=dvx[0];           
    217226                D_scalar=dt*gauss->weight*Jdet;
    218227                if(dim==2){
    219228                        vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
    220                         dvxdx=dvx[0];
    221229                        dvydy=dvy[1];
    222230                        for(int i=0;i<numnodes;i++){
    … …  
    229237                        }
    230238                }else{
    231                         dvxdx=dvx[0];
    232239                        for(int i=0;i<numnodes;i++){
    233240                                for(int j=0;j<numnodes;j++){
    … …  
    238245                }
    239246
    240                 IssmDouble rho;
     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;
    241301                if(FINITEELEMENT==P1Enum){
    242                         rho=2.;
     302                        rho=RHO;
    243303                }else if(FINITEELEMENT==P2Enum){
    244                         rho=4.;
     304                        rho=RHO*2.;
    245305                }
    246306
    … …  
    249309                        /*SSA*/
    250310                        if(dim==1){
    251                                 vx_input->GetInputValue(&vx,gauss);
     311                                //vx_input->GetInputValue(&vx,gauss);
     312                                vx_input->GetInputAverage(&vx);
    252313                                D[0]=h/rho*fabs(vx);
    253314                        }else{
    … …  
    262323                        if(dim==1){
    263324                                vx_input->GetInputValue(&vx,gauss);
     325                                //vx_input->GetInputAverage(&vx);
    264326                                vel=fabs(vx)+EPS;
    265327                                D[0] = h/(rho*vel)*vx*vx;
    266328                        }else{
    267329                                vx_input->GetInputValue(&vx,gauss);
     330                                //vx_input->GetInputAverage(&vx);
    268331                                vy_input->GetInputValue(&vy,gauss);
     332                                //vy_input->GetInputAverage(&vy);
    269333                                vel=sqrt(vx*vx+vy*vy)+EPS;
    270334                                D[0*dim+0]=h/(rho*vel)*vx*vx;
    … …  
    277341                        if(dim==1){
    278342                                vx_input->GetInputValue(&vx,gauss);
    279                                 tau=h/(rho*max(fabs(vx),EPS));
     343                                //vx_input->GetInputAverage(&vx);
     344                                //tau=h/(rho*max(fabs(vx),EPS));
     345                                tau=h/(rho*(fabs(vx)+EPS));
    280346                        }else{
    281347                                vx_input->GetInputValue(&vx,gauss);
    282348                                vy_input->GetInputValue(&vy,gauss);
    283                                 tau=h/(rho*sqrt(vx*vx+vy*vy)+EPS);
     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;
    284367                        }
    285368                }
    … …  
    327410                                }
    328411
    329 
    330412                                /*Advection matrix - part 2, B*/
    331413                                for(int i=0;i<numnodes;i++){
    … …  
    338420                                for(int i=0;i<numnodes;i++){
    339421                                        for(int j=0;j<numnodes;j++){
    340                                                 Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);;
     422                                                Ke->values[i*numnodes+j]+=dt*gauss->weight*Jdet*tau*(basis[j]*dvxdx)*(basis[i]*dvxdx);
    341423                                        }
    342424                                }
    … …  
    442524        Input* smb_input        = topelement->GetInput(SmbMassBalanceEnum);  _assert_(smb_input);
    443525        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);
    444528        Input* vx_input  = topelement->GetInput(VxDebrisEnum);  _assert_(vx_input);
    445529        Input* vy_input=NULL;
    … …  
    448532        }
    449533        h=topelement->CharacteristicLength();
    450 
    451         IssmDouble rho;
     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;
    452539        if(FINITEELEMENT==P1Enum){
    453                 rho=2.;
     540                rho=RHO;
    454541        }else if(FINITEELEMENT==P2Enum){
    455                 rho=4.;
     542                rho=RHO*2.;
    456543        }
    457544
    … …  
    466553                thickness_input->GetInputValue(&thickness,gauss);
    467554                if(smb>0.){
    468                         psi=thickness-0.*dt*smb*pf;
     555                        psi=thickness-0*dt*smb*pf;
    469556                }else{
    470557                        psi=thickness-dt*smb*pf;
    471558                }
    472 
     559                //psi=thickness-dt*smb*pf;
    473560                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*psi*basis[i];
    474561
    … …  
    480567                        if(dim==1){
    481568                                vx_input->GetInputValue(&vx,gauss);
    482                                 tau=h/(rho*max(fabs(vx),EPS));
     569                                //vx_input->GetInputAverage(&vx);
     570                                //tau=h/(rho*max(fabs(vx),EPS));
     571                                tau=h/(rho*(fabs(vx)+EPS));
    483572
    484573                                /*Force vector - part 2*/
    … …  
    494583                                vx_input->GetInputValue(&vx,gauss);
    495584                                vy_input->GetInputValue(&vy,gauss);
     585                                vx_input->GetInputAverage(&vx);
     586                                vy_input->GetInputAverage(&vy);
    496587                                vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
    497                                 vel=sqrt(vx*vx+vy*vy);
    498                                 dvydy=dvy[1];
    499                                 tau=h/(rho*vel+EPS);
     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;
    500605
    501606                                /*Force vector - part 2*/
    … …  
    531636        int numnodes = element->GetNumberOfNodes();
    532637        IssmDouble* newthickness  = xNew<IssmDouble>(numnodes);
     638        IssmDouble* icethickness  = xNew<IssmDouble>(numnodes);
     639       
    533640
    534641        /*Use the dof list to index into the solution vector: */
    535642        IssmDouble minthickness = element->FindParam(DebrisMinThicknessEnum);
    536643        element->GetDofListLocal(&ddoflist,NoneApproximationEnum,GsetEnum);
     644        element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    537645
    538646        for(int i=0;i<numnodes;i++){
    … …  
    550658        // Free resources
    551659        xDelete<IssmDouble>(newthickness);
     660        xDelete<IssmDouble>(icethickness);
    552661        xDelete<int>(ddoflist);
    553662        //*/
    … …  
    599708        femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum);
    600709        Element* element= NULL;
    601         Element*    topelement = NULL;
    602710
    603711        if(removalmodel==0){
    … …  
    616724                        IssmDouble* slopey         = xNew<IssmDouble>(numnodes);
    617725                        IssmDouble* onsurface      = xNew<IssmDouble>(numnodes);
    618                         IssmDouble* ls_active      = xNew<IssmDouble>(numnodes);
    619                         element->GetInputListOnNodes(debristhickness,DebrisThicknessEnum);
     726                        //IssmDouble* ls_active      = xNew<IssmDouble>(numnodes);
     727                        element->GetInputListOnNodes(debristhickness,DebrisThicknessOldEnum);
    620728                        element->GetInputListOnNodes(icethickness,ThicknessEnum);
    621729                        element->GetInputListOnNodes(onsurface,MeshVertexonsurfaceEnum);
    622                         element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
     730                        //element->GetInputListOnNodes(ls_active,DebrisMaskNodeActivationEnum);
    623731
    624732                        dim=1;
    … …  
    643751                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
    644752                                               }
    645                                                isminthicknessinelement=true;
    646                                                if(kk<numnodes && isminthicknessinelement){
     753                                               //isminthicknessinelement=true;
     754                                               //if(kk<numnodes && isminthicknessinelement){
     755                                               if(isminthicknessinelement){                                           
    647756                                                       for(k=0; k<numnodes;k++){
    648757                                                               slope=fabs(slopex[k]);
    … …  
    659768                                                       }
    660769                                               }
    661                                                element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
     770                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1DGEnum);
     771
     772                                               xDelete<IssmDouble>(debristhickness);
     773                                               xDelete<IssmDouble>(onsurface);
     774                                               xDelete<IssmDouble>(icethickness);
     775                                               xDelete<IssmDouble>(slopex);
     776                                               xDelete<IssmDouble>(slopey);
     777                                               //xDelete<IssmDouble>(ls_active);
     778                                               break;
     779                                       }
     780                                case 2:{
     781                                               IssmDouble stress_threshold=element->FindParam(DebrisRemovalStressThresholdEnum);
     782                                               IssmDouble gravity=element->FindParam(ConstantsGEnum);
     783                                               IssmDouble stress,rhod=1900.;
     784                                               int kk=0;
     785                                               for(k=0; k<numnodes;k++){
     786                                                       if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
     787                                                       if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
     788                                               }
     789                                               //isminthicknessinelement=true;
     790                                               //if(isminthicknessinelement){
     791                                               if(kk<numnodes && isminthicknessinelement){
     792                                                       for(k=0; k<numnodes;k++){
     793                                                               slope=fabs(slopex[k]);
     794                                                               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;
     797                                                               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                                                        }
     804                                               }
     805                                               element->AddInput(DebrisThicknessEnum,debristhickness,P1DGEnum);
    662806
    663807                                               xDelete<IssmDouble>(debristhickness);
    … …  
    667811                                               break;
    668812                                       }
    669                                 case 2:{
    670                                                IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
    671                                                IssmDouble gravity = element->FindParam(ConstantsGEnum);
    672                                                IssmDouble stress,rhod=1900.;
    673                                                int kk=0;
    674                                                for(k=0; k<numnodes;k++){
    675                                                        if(icethickness[k]<=(iceminthickness+0.00001)) isminthicknessinelement=true;
    676                                                        if(icethickness[k]<=(iceminthickness+0.00001)) kk++;
    677                                                }
    678                                                isminthicknessinelement=true;
    679                                                if(kk<numnodes && isminthicknessinelement){
    680                                                        //stress=0;
    681                                                        IssmDouble stress_sum=0.;
    682                                                        for(k=0; k<numnodes;k++){
    683                                                                slope=fabs(slopex[k]);
    684                                                                if(dim==2) slope=pow(pow(slopex[k],2)+pow(slopey[k],2),0.5);
    685                                                                stress=rhod*gravity*debristhickness[k]*slope;//pow(slope*slope/(slope*slope+1),0.5);//sin(slope/rad2deg);
    686                                                                //stress_sum=stress_sum+stress;
    687                                                                if(stress>stress_threshold) remove_debris=true;
    688                                                                //if(stress>stress_threshold) debristhickness[k]=0.;
    689                                                        }
    690                                                        //if((stress_sum/double(kk))>stress_threshold) remove_debris=true;
    691                                                        if(remove_debris){
    692                                                                for(k=0; k<numnodes;k++){
    693                                                                        debristhickness[k]=0.;
    694                                                                }
    695                                                        }
    696                                                }
    697                                                element->AddInput(DebrisThicknessEnum,debristhickness,P1Enum);
    698 
    699                                                xDelete<IssmDouble>(debristhickness);
    700                                                xDelete<IssmDouble>(icethickness);
    701                                                xDelete<IssmDouble>(slopex);
    702                                                xDelete<IssmDouble>(slopey);
    703                                                xDelete<IssmDouble>(ls_active);
    704                                                break;
    705                                        }
    706813                                default: _error_("removalmodel "<<EnumToStringx(removalmodel)<<" not implemented yet");
    707814                        }
    … …  
    713820
    714821        if(VerboseSolution()) _printf0_("   Debris preprocessing\n");
     822
     823        int removalmodel, displacementmodel;
     824        femmodel->parameters->FindParam(&removalmodel,DebrisRemovalmodelEnum);
     825        femmodel->parameters->FindParam(&displacementmodel,DebrisDisplacementmodelEnum);
    715826
    716827        Element* element= NULL;
    … …  
    718829                element=xDynamicCast<Element*>(object);
    719830
    720                 int numvertices = element->GetNumberOfVertices();
     831                int numvertices=element->GetNumberOfVertices();
     832                int numnodes=element->GetNumberOfNodes();
     833                numvertices=numnodes;
    721834
    722835                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);             
    723839                IssmDouble* debristhickness = xNew<IssmDouble>(numvertices);
    724840                IssmDouble* slopex         = xNew<IssmDouble>(numvertices);
    725                 IssmDouble* onsurface      = xNew<IssmDouble>(numvertices);
    726                 IssmDouble* icethickness      = xNew<IssmDouble>(numvertices);
    727 
     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                }
    728848                element->GetInputListOnVertices(&debristhickness[0],DebrisThicknessEnum);
    729                 element->GetInputListOnVertices(&vx[0],VxDebrisEnum);
     849                element->GetInputListOnVertices(&vx[0],VxEnum);
    730850                element->GetInputListOnVertices(&slopex[0],SurfaceSlopeXEnum);
    731                 element->GetInputListOnVertices(&onsurface[0],MeshVertexonsurfaceEnum);
    732851                element->GetInputListOnVertices(&icethickness[0],ThicknessEnum);
    733 
    734                 IssmDouble slope,rad2deg=180./M_PI; //=57.2958
    735                 IssmDouble vslipx,rhod=1900.;
     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;
    736860                IssmDouble gravity=element->FindParam(ConstantsGEnum);
    737861                IssmDouble slope_threshold=element->FindParam(DebrisRemovalSlopeThresholdEnum);
     862                IssmDouble stress_threshold = element->FindParam(DebrisRemovalStressThresholdEnum);
     863                IssmDouble maxv = element->FindParam(DebrisMaxdisplacementvelocityEnum);
    738864                IssmDouble iceminthickness=element->FindParam(MasstransportMinThicknessEnum);
     865                IssmDouble hx,hy,hz;
     866                element->ElementSizes(&hx,&hy,&hz);
    739867
    740868                int step;
    741                 IssmDouble dt, maxv;
     869                IssmDouble dt;
    742870                IssmDouble yts=31536000.;
    743871                femmodel->parameters->FindParam(&step,StepEnum);
    744872                femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    745 
    746873                bool isminthicknessinelement=false;
    747                 for(int i=0;i<numvertices;i++){
    748                         if(icethickness[i]<=(iceminthickness+0.01)) isminthicknessinelement=true;
    749                 }
    750                 if(isminthicknessinelement){
    751                         //do nothing
     874                int kk=0;               
     875
     876                if(displacementmodel==1){
     877
    752878                        for(int i=0;i<numvertices;i++){
    753                                 if(icethickness[i]<=(iceminthickness+0.01)) vx[i]=0.;
    754                         }
    755                 }else{
     879                                if(icethickness[i]<=(iceminthickness+0.00001)){
     880                                        isminthicknessinelement=true;
     881                                        kk++;
     882                                }
     883                        }
     884                        //isminthicknessinelement=true;
    756885                        for(int i=0;i<numvertices;i++){
    757                                 //if(onsurface[i]>.5){
    758                                 slope=fabs(slopex[i]);
    759                                 //if((atan(slope)*rad2deg)>25.){
    760                                 //if(debristhickness[i]>0.01){
    761                                 vslipx=1.0/yts;
    762                                 //maxv=10.0/2./dt;
    763                                 //vslipx=-slope_threshold*rhod*gravity*debristhickness[i]*slopex[i]/yts;
    764                                 vx[i]=vx[i]+vslipx;
    765                                 //debristhickness[i]=debristhickness[i];
    766                                 //if(vx[i]>maxv) vx[i]=maxv;
    767                                 //}
    768                                 //}
    769                                 //}
    770                         }
    771                 }
    772                 //if(step%100==0)   
    773                 element->AddInput(VxDebrisEnum,vx,P1Enum);
    774 
     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                }
    775939                /* Free resources */
    776940                xDelete<IssmDouble>(debristhickness);
    777941                xDelete<IssmDouble>(icethickness);
    778942                xDelete<IssmDouble>(vx);
     943                xDelete<IssmDouble>(vx_new);
    779944                xDelete<IssmDouble>(slopex);
    780                 xDelete<IssmDouble>(onsurface);
    781         }
    782 
    783 }/*}}}*/
     945                xDelete<IssmDouble>(vy);
     946                xDelete<IssmDouble>(vy_new);
     947                xDelete<IssmDouble>(slopey);
     948        }
     949}/*}}}*/

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

    r27646 r27847  
    242242                        }
    243243                        break;
    244                 case SMBdebrisMLEnum:
    245                         iomodel->FetchDataToInput(inputs,elements,"md.initialization.debris",DebrisThicknessEnum);
    246                         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;
    247261                default:
    248262                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    … …  
    473487                        /*Nothing to add to parameters*/
    474488                        break;
    475                 case SMBdebrisMLEnum:
    476                         /*Nothing to add to parameters*/
    477                         break;
     489                case SMBdebrisEvattEnum:
     490                        /*Nothing to add to parameters*/
     491                        break;
    478492                default:
    479493                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    … …  
    573587                        #endif //_HAVE_SEMIC_
    574588                        break;
    575                 case SMBdebrisMLEnum:
    576                         if(VerboseSolution())_printf0_("        call smb debris Mayer & Liculli module\n");
    577                         SmbDebrisMLx(femmodel);
    578                         break;
     589                case SMBdebrisEvattEnum:
     590                        if(VerboseSolution())_printf0_("        call smb Evatt debris module\n");
     591                        SmbDebrisEvattx(femmodel);
     592                        break;
    579593                default:
    580594                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");

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

    r27844 r27847  
    2323#include "../Inputs/ArrayInput.h"
    2424#include "../Inputs/IntArrayInput.h"
     25#include <cstdlib>
    2526/*}}}*/
    2627#define MAXVERTICES 6 /*Maximum number of vertices per element, currently Penta, to avoid dynamic mem allocation*/
    … …  
    23502351
    23512352        IssmDouble minthickness = this->FindParam(MasstransportMinThicknessEnum);
    2352 
     2353       
    23532354        Input* input=this->GetInput(ThicknessEnum); _assert_(input);
    2354         if(input->GetInputMax()<=(minthickness+0.00000001)){
     2355        if(input->GetInputMin()<=(minthickness+0.00000001)){
    23552356                return true;
    23562357        }
    … …  
    34003401                        monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
    34013402                        dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month);
    3402                         monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
     3403                        monthlyprec[iv*12+month]=monthlyprec[iv*12+month];//*yts;
    34033404                }
    34043405        }
    … …  
    35943595        Gauss* gauss=this->NewGauss();
    35953596        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/*}}}*/
     3731void       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);
    37294059}
    37304060/*}}}*/

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

    r27829 r27847  
    177177                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    178178                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
     179                void               SmbDebrisEvatt();
    179180                void               RignotMeltParameterization();
    180181                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);
    … …  
    256257                virtual void       ComputeStressTensor(void)=0;
    257258                virtual void       ComputeEsaStrainAndVorticity(void)=0;
     259                virtual void       ComputeMeanEla(IssmDouble* paltitude,int* pcounter)=0;
    258260                virtual void       Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters,Inputs* inputsin)=0;
    259261                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

    r27829 r27847  
    52835283/*}}}*/
    52845284#endif
     5285void       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

    r27829 r27847  
    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");};
    5253                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs* inputsin){_error_("not implemented yet");};
    5354                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

    r27829 r27847  
    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");};
    4849                void        ComputeDeviatoricStressTensor(){_error_("not implemented yet");};
    4950                void        ComputeEsaStrainAndVorticity(){_error_("not implemented yet!");};

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

    r27829 r27847  
    7272                void        ComputeStressTensor();
    7373                void        ComputeSurfaceNormalVelocity();
     74                void        ComputeMeanEla(IssmDouble* paltitude, int* pcounter){_error_("not implemented yet");};
    7475                void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters,Inputs* inputsin);
    7576                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

    r27795 r27847  
    2727
    2828        /*recover parameters: */
     29        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    2930        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    3031        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    … …  
    4344                femmodel->inputs->DuplicateInput(VyEnum,VyDebrisEnum); 
    4445        }
    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         }
    5346
    5447        debris_analysis = new DebrisAnalysis();
    … …  
    5952        extrudefromtop_core(femmodel); 
    6053
     54        IssmDouble mean_ela;
     55        ComputeMeanElax(&mean_ela,femmodel);
     56        femmodel->parameters->SetParam(mean_ela,SmbMeanElaEnum);
     57        _printf_("core ELA: " << mean_ela << "\n");
     58       
    6159        if(save_results) femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    6260        if(solution_type==DebrisSolutionEnum)femmodel->RequestedDependentsx();

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

    r27780 r27847  
    489489                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismethod",SmbSemicMethodEnum));
    490490                        break;
    491                 case SMBdebrisMLEnum:
     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));
    492505                        break;
    493506                default:

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

    r27498 r27847  
    504504
    505505}/*}}}*/
    506 void 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         }
     506void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
     507        for(Object* & object : femmodel->elements->objects){
     508                Element* element=xDynamicCast<Element*>(object);
     509                element->SmbDebrisEvatt();
     510        }
    647511}/*}}}*/
    648512void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/

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

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

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

    r27146 r27847  
    8686#include "./ResetConstraintsx/ResetConstraintsx.h"
    8787#include "./ResetFSBasalBoundaryConditionx/ResetFSBasalBoundaryConditionx.h"
     88#include "./ComputeMeanElax/ComputeMeanElax.h"
    8889#include "./RheologyBbarAbsGradientx/RheologyBbarAbsGradientx.h"
    8990#include "./RheologyBAbsGradientx/RheologyBAbsGradientx.h"

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

    r27829 r27847  
    176176syn keyword cConstant DebrisRemovalStressThresholdEnum
    177177syn keyword cConstant DebrisPackingFractionEnum
     178syn keyword cConstant DebrisMaxdisplacementvelocityEnum
    178179syn keyword cConstant DebugProfilingEnum
    179180syn keyword cConstant DomainDimensionEnum
    … …  
    539540syn keyword cConstant SmbARMAmaOrderEnum
    540541syn keyword cConstant SmbAveragingEnum
     542syn keyword cConstant SmbDebrisalbedoEnum
     543syn keyword cConstant SmbDebrisIsMonthlyEnum
     544syn keyword cConstant SmbDebrisIsAndersonEnum
     545syn keyword cConstant SmbDebrisIsCryokarstEnum
     546syn keyword cConstant SmbDebrisAndersonD0Enum
    541547syn keyword cConstant SmbDesfacEnum
    542548syn keyword cConstant SmbDesfacElevEnum
    … …  
    552558syn keyword cConstant SmbEIdxEnum
    553559syn keyword cConstant SmbFEnum
     560syn keyword cConstant SmbHumiditylapsEnum
    554561syn keyword cConstant SmbInitDensityScalingEnum
     562syn keyword cConstant SmbIcealbedoEnum
    555563syn keyword cConstant SmbIsaccumulationEnum
    556564syn keyword cConstant SmbIsalbedoEnum
    … …  
    572580syn keyword cConstant SmbKEnum
    573581syn keyword cConstant SmbLapseRatesEnum
     582syn keyword cConstant SmbLWlapsEnum
     583syn keyword cConstant SmbSWlapsEnum
     584syn keyword cConstant SmbSnowalbedoEnum
     585syn keyword cConstant SmbMeanElaEnum
    574586syn keyword cConstant SmbNumBasinsEnum
    575587syn keyword cConstant SmbNumBreaksEnum
    … …  
    612624syn keyword cConstant SmbTemperaturesReconstructedYearsEnum
    613625syn keyword cConstant SmbPrecipitationsReconstructedYearsEnum
     626syn keyword cConstant SmbWindspeedlapsEnum
    614627syn keyword cConstant SmoothThicknessMultiplierEnum
    615628syn keyword cConstant SolutionTypeEnum
    … …  
    779792syn keyword cConstant DamageFEnum
    780793syn keyword cConstant DebrisThicknessEnum
     794syn keyword cConstant DebrisThicknessOldEnum
    781795syn keyword cConstant DegreeOfChannelizationEnum
    782796syn keyword cConstant DepthBelowSurfaceEnum
    … …  
    11071121syn keyword cConstant SmbMeltEnum
    11081122syn keyword cConstant SmbMonthlytemperaturesEnum
     1123syn keyword cConstant SmbMonthlydsradiationEnum
     1124syn keyword cConstant SmbMonthlydlradiationEnum
     1125syn keyword cConstant SmbMonthlywindspeedEnum
     1126syn keyword cConstant SmbMonthlyairhumidityEnum
    11091127syn keyword cConstant SmbMSurfEnum
    11101128syn keyword cConstant SmbNetLWEnum
    … …  
    11161134syn keyword cConstant SmbPrecipitationEnum
    11171135syn keyword cConstant SmbPrecipitationsAnomalyEnum
     1136syn keyword cConstant SmbDsradiationAnomalyEnum
     1137syn keyword cConstant SmbDlradiationAnomalyEnum
     1138syn keyword cConstant SmbWindspeedAnomalyEnum
     1139syn keyword cConstant SmbAirhumidityAnomalyEnum
    11181140syn keyword cConstant SmbPrecipitationsLgmEnum
    11191141syn keyword cConstant SmbPrecipitationsPresentdayEnum
    … …  
    11351157syn keyword cConstant SmbSmbrefEnum
    11361158syn keyword cConstant SmbSzaValueEnum
     1159syn keyword cConstant SmbSummerMeltEnum
     1160syn keyword cConstant SmbSummerAlbedoEnum
     1161syn keyword cConstant SmbSnowheightEnum
    11371162syn keyword cConstant SmbTEnum
    11381163syn keyword cConstant SmbTaEnum
    … …  
    16541679syn keyword cConstant SMBarmaEnum
    16551680syn keyword cConstant SMBcomponentsEnum
    1656 syn keyword cConstant SMBdebrisMLEnum
     1681syn keyword cConstant SMBdebrisEvattEnum
    16571682syn keyword cConstant SMBd18opddEnum
    16581683syn keyword cConstant SMBforcingEnum
    … …  
    17801805syn keyword cType Cfsurfacesquaretransient
    17811806syn keyword cType Channel
     1807syn keyword cType classes
    17821808syn keyword cType Constraint
    17831809syn keyword cType Constraints
    … …  
    17871813syn keyword cType ControlParam
    17881814syn keyword cType Covertree
     1815syn keyword cType DatasetInput
    17891816syn keyword cType DataSetParam
    1790 syn keyword cType DatasetInput
    17911817syn keyword cType Definition
    17921818syn keyword cType DependentObject
    … …  
    18011827syn keyword cType ElementInput
    18021828syn keyword cType ElementMatrix
     1829syn keyword cType Elements
    18031830syn keyword cType ElementVector
    1804 syn keyword cType Elements
    18051831syn keyword cType ExponentialVariogram
    18061832syn keyword cType ExternalResult
    … …  
    18091835syn keyword cType Friction
    18101836syn keyword cType Gauss
     1837syn keyword cType GaussianVariogram
     1838syn keyword cType gaussobjects
    18111839syn keyword cType GaussPenta
    18121840syn keyword cType GaussSeg
    18131841syn keyword cType GaussTetra
    18141842syn keyword cType GaussTria
    1815 syn keyword cType GaussianVariogram
    18161843syn keyword cType GenericExternalResult
    18171844syn keyword cType GenericOption
    … …  
    18301857syn keyword cType IssmDirectApplicInterface
    18311858syn keyword cType IssmParallelDirectApplicInterface
     1859syn keyword cType krigingobjects
    18321860syn keyword cType Load
    18331861syn keyword cType Loads
    … …  
    18401868syn keyword cType Matice
    18411869syn keyword cType Matlitho
     1870syn keyword cType matrixobjects
    18421871syn keyword cType MatrixParam
    18431872syn keyword cType Misfit
    … …  
    18521881syn keyword cType Observations
    18531882syn keyword cType Option
     1883syn keyword cType Options
    18541884syn keyword cType OptionUtilities
    1855 syn keyword cType Options
    18561885syn keyword cType Param
    18571886syn keyword cType Parameters
    … …  
    18671896syn keyword cType Regionaloutput
    18681897syn keyword cType Results
     1898syn keyword cType Riftfront
    18691899syn keyword cType RiftStruct
    1870 syn keyword cType Riftfront
    18711900syn keyword cType SealevelGeometry
    18721901syn keyword cType Seg
    18731902syn keyword cType SegInput
     1903syn keyword cType Segment
    18741904syn keyword cType SegRef
    1875 syn keyword cType Segment
    18761905syn keyword cType SpcDynamic
    18771906syn keyword cType SpcStatic
    … …  
    18921921syn keyword cType Vertex
    18931922syn keyword cType Vertices
    1894 syn keyword cType classes
    1895 syn keyword cType gaussobjects
    1896 syn keyword cType krigingobjects
    1897 syn keyword cType matrixobjects
    18981923syn keyword cType AdjointBalancethickness2Analysis
    18991924syn keyword cType AdjointBalancethicknessAnalysis

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

    r27829 r27847  
    170170        DebrisRemovalStressThresholdEnum,
    171171        DebrisPackingFractionEnum,
     172        DebrisMaxdisplacementvelocityEnum,
    172173        DebugProfilingEnum,
    173174        DomainDimensionEnum,
    … …  
    530531        SmbAlbedoLandEnum,
    531532        SmbARMATimestepEnum,
    532    SmbARMAarOrderEnum,
    533    SmbARMAmaOrderEnum,
     533        SmbARMAarOrderEnum,
     534        SmbARMAmaOrderEnum,
    534535        SmbAveragingEnum,
     536        SmbDebrisalbedoEnum,
     537        SmbDebrisIsMonthlyEnum,
     538        SmbDebrisIsAndersonEnum,
     539        SmbDebrisIsCryokarstEnum,
     540        SmbDebrisAndersonD0Enum,
    535541        SmbDesfacEnum,
    536542        SmbDesfacElevEnum,
    537543        SmbDpermilEnum,
    538544        SmbDsnowIdxEnum,
    539    SmbElevationBinsEnum,
     545        SmbElevationBinsEnum,
    540546        SmbCldFracEnum,
    541547        SmbDelta18oEnum,
    … …  
    546552        SmbEIdxEnum,
    547553        SmbFEnum,
     554        SmbHumiditylapsEnum,
    548555        SmbInitDensityScalingEnum,
     556        SmbIcealbedoEnum,
    549557        SmbIsaccumulationEnum,
    550558        SmbIsalbedoEnum,
    … …  
    565573        SmbIsturbulentfluxEnum,
    566574        SmbKEnum,
    567    SmbLapseRatesEnum,
     575        SmbLapseRatesEnum,
     576        SmbLWlapsEnum,
     577        SmbSWlapsEnum,
     578        SmbSnowalbedoEnum,
     579        SmbMeanElaEnum,
    568580        SmbNumBasinsEnum,
    569581        SmbNumBreaksEnum,
    … …  
    606618        SmbTemperaturesReconstructedYearsEnum,
    607619        SmbPrecipitationsReconstructedYearsEnum,
     620        SmbWindspeedlapsEnum,
    608621        SmoothThicknessMultiplierEnum,
    609622        SolutionTypeEnum,
    … …  
    775788        DamageFEnum,
    776789        DebrisThicknessEnum,
     790        DebrisThicknessOldEnum,
    777791        DegreeOfChannelizationEnum,
    778792        DepthBelowSurfaceEnum,
    … …  
    10431057        SmbAdiffiniEnum,
    10441058        SmbAiniEnum,
    1045    SmbARMANoiseEnum,
     1059        SmbARMANoiseEnum,
    10461060        SmbBasinsIdEnum,
    10471061        SmbBMaxEnum,
    … …  
    11041118        SmbMInitnum,
    11051119        SmbMonthlytemperaturesEnum,
     1120        SmbMonthlydsradiationEnum,
     1121        SmbMonthlydlradiationEnum,
     1122        SmbMonthlywindspeedEnum,
     1123        SmbMonthlyairhumidityEnum,
    11061124        SmbMSurfEnum,
    11071125        SmbNetLWEnum,
    … …  
    11131131        SmbPrecipitationEnum,
    11141132        SmbPrecipitationsAnomalyEnum,
     1133        SmbDsradiationAnomalyEnum,
     1134        SmbDlradiationAnomalyEnum,
     1135        SmbWindspeedAnomalyEnum,
     1136        SmbAirhumidityAnomalyEnum,
    11151137        SmbPrecipitationsLgmEnum,
    11161138        SmbPrecipitationsPresentdayEnum,
    … …  
    11321154        SmbSmbrefEnum,
    11331155        SmbSzaValueEnum,
     1156        SmbSummerMeltEnum,
     1157        SmbSummerAlbedoEnum,
     1158        SmbSnowheightEnum,
    11341159        SmbTEnum,
    11351160        SmbTaEnum,
    … …  
    16531678        SMBarmaEnum,
    16541679        SMBcomponentsEnum,
    1655         SMBdebrisMLEnum,
     1680        SMBdebrisEvattEnum,
    16561681        SMBd18opddEnum,
    16571682        SMBforcingEnum,

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

    r27829 r27847  
    178178                case DebrisRemovalStressThresholdEnum : return "DebrisRemovalStressThreshold";
    179179                case DebrisPackingFractionEnum : return "DebrisPackingFraction";
     180                case DebrisMaxdisplacementvelocityEnum : return "DebrisMaxdisplacementvelocity";
    180181                case DebugProfilingEnum : return "DebugProfiling";
    181182                case DomainDimensionEnum : return "DomainDimension";
    … …  
    541542                case SmbARMAmaOrderEnum : return "SmbARMAmaOrder";
    542543                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";
    543549                case SmbDesfacEnum : return "SmbDesfac";
    544550                case SmbDesfacElevEnum : return "SmbDesfacElev";
    … …  
    554560                case SmbEIdxEnum : return "SmbEIdx";
    555561                case SmbFEnum : return "SmbF";
     562                case SmbHumiditylapsEnum : return "SmbHumiditylaps";
    556563                case SmbInitDensityScalingEnum : return "SmbInitDensityScaling";
     564                case SmbIcealbedoEnum : return "SmbIcealbedo";
    557565                case SmbIsaccumulationEnum : return "SmbIsaccumulation";
    558566                case SmbIsalbedoEnum : return "SmbIsalbedo";
    … …  
    574582                case SmbKEnum : return "SmbK";
    575583                case SmbLapseRatesEnum : return "SmbLapseRates";
     584                case SmbLWlapsEnum : return "SmbLWlaps";
     585                case SmbSWlapsEnum : return "SmbSWlaps";
     586                case SmbSnowalbedoEnum : return "SmbSnowalbedo";
     587                case SmbMeanElaEnum : return "SmbMeanEla";
    576588                case SmbNumBasinsEnum : return "SmbNumBasins";
    577589                case SmbNumBreaksEnum : return "SmbNumBreaks";
    … …  
    614626                case SmbTemperaturesReconstructedYearsEnum : return "SmbTemperaturesReconstructedYears";
    615627                case SmbPrecipitationsReconstructedYearsEnum : return "SmbPrecipitationsReconstructedYears";
     628                case SmbWindspeedlapsEnum : return "SmbWindspeedlaps";
    616629                case SmoothThicknessMultiplierEnum : return "SmoothThicknessMultiplier";
    617630                case SolutionTypeEnum : return "SolutionType";
    … …  
    781794                case DamageFEnum : return "DamageF";
    782795                case DebrisThicknessEnum : return "DebrisThickness";
     796                case DebrisThicknessOldEnum : return "DebrisThicknessOld";
    783797                case DegreeOfChannelizationEnum : return "DegreeOfChannelization";
    784798                case DepthBelowSurfaceEnum : return "DepthBelowSurface";
    … …  
    11091123                case SmbMeltEnum : return "SmbMelt";
    11101124                case SmbMonthlytemperaturesEnum : return "SmbMonthlytemperatures";
     1125                case SmbMonthlydsradiationEnum : return "SmbMonthlydsradiation";
     1126                case SmbMonthlydlradiationEnum : return "SmbMonthlydlradiation";
     1127                case SmbMonthlywindspeedEnum : return "SmbMonthlywindspeed";
     1128                case SmbMonthlyairhumidityEnum : return "SmbMonthlyairhumidity";
    11111129                case SmbMSurfEnum : return "SmbMSurf";
    11121130                case SmbNetLWEnum : return "SmbNetLW";
    … …  
    11181136                case SmbPrecipitationEnum : return "SmbPrecipitation";
    11191137                case SmbPrecipitationsAnomalyEnum : return "SmbPrecipitationsAnomaly";
     1138                case SmbDsradiationAnomalyEnum : return "SmbDsradiationAnomaly";
     1139                case SmbDlradiationAnomalyEnum : return "SmbDlradiationAnomaly";
     1140                case SmbWindspeedAnomalyEnum : return "SmbWindspeedAnomaly";
     1141                case SmbAirhumidityAnomalyEnum : return "SmbAirhumidityAnomaly";
    11201142                case SmbPrecipitationsLgmEnum : return "SmbPrecipitationsLgm";
    11211143                case SmbPrecipitationsPresentdayEnum : return "SmbPrecipitationsPresentday";
    … …  
    11371159                case SmbSmbrefEnum : return "SmbSmbref";
    11381160                case SmbSzaValueEnum : return "SmbSzaValue";
     1161                case SmbSummerMeltEnum : return "SmbSummerMelt";
     1162                case SmbSummerAlbedoEnum : return "SmbSummerAlbedo";
     1163                case SmbSnowheightEnum : return "SmbSnowheight";
    11391164                case SmbTEnum : return "SmbT";
    11401165                case SmbTaEnum : return "SmbTa";
    … …  
    16561681                case SMBarmaEnum : return "SMBarma";
    16571682                case SMBcomponentsEnum : return "SMBcomponents";
    1658                 case SMBdebrisMLEnum : return "SMBdebrisML";
     1683                case SMBdebrisEvattEnum : return "SMBdebrisEvatt";
    16591684                case SMBd18opddEnum : return "SMBd18opdd";
    16601685                case SMBforcingEnum : return "SMBforcing";

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

    r27829 r27847  
    169169syn keyword juliaConstC DebrisRemovalStressThresholdEnum
    170170syn keyword juliaConstC DebrisPackingFractionEnum
     171syn keyword juliaConstC DebrisMaxdisplacementvelocityEnum
    171172syn keyword juliaConstC DebugProfilingEnum
    172173syn keyword juliaConstC DomainDimensionEnum
    … …  
    532533syn keyword juliaConstC SmbARMAmaOrderEnum
    533534syn keyword juliaConstC SmbAveragingEnum
     535syn keyword juliaConstC SmbDebrisalbedoEnum
     536syn keyword juliaConstC SmbDebrisIsMonthlyEnum
     537syn keyword juliaConstC SmbDebrisIsAndersonEnum
     538syn keyword juliaConstC SmbDebrisIsCryokarstEnum
     539syn keyword juliaConstC SmbDebrisAndersonD0Enum
    534540syn keyword juliaConstC SmbDesfacEnum
    535541syn keyword juliaConstC SmbDesfacElevEnum
    … …  
    545551syn keyword juliaConstC SmbEIdxEnum
    546552syn keyword juliaConstC SmbFEnum
     553syn keyword juliaConstC SmbHumiditylapsEnum
    547554syn keyword juliaConstC SmbInitDensityScalingEnum
     555syn keyword juliaConstC SmbIcealbedoEnum
    548556syn keyword juliaConstC SmbIsaccumulationEnum
    549557syn keyword juliaConstC SmbIsalbedoEnum
    … …  
    565573syn keyword juliaConstC SmbKEnum
    566574syn keyword juliaConstC SmbLapseRatesEnum
     575syn keyword juliaConstC SmbLWlapsEnum
     576syn keyword juliaConstC SmbSWlapsEnum
     577syn keyword juliaConstC SmbSnowalbedoEnum
     578syn keyword juliaConstC SmbMeanElaEnum
    567579syn keyword juliaConstC SmbNumBasinsEnum
    568580syn keyword juliaConstC SmbNumBreaksEnum
    … …  
    605617syn keyword juliaConstC SmbTemperaturesReconstructedYearsEnum
    606618syn keyword juliaConstC SmbPrecipitationsReconstructedYearsEnum
     619syn keyword juliaConstC SmbWindspeedlapsEnum
    607620syn keyword juliaConstC SmoothThicknessMultiplierEnum
    608621syn keyword juliaConstC SolutionTypeEnum
    … …  
    772785syn keyword juliaConstC DamageFEnum
    773786syn keyword juliaConstC DebrisThicknessEnum
     787syn keyword juliaConstC DebrisThicknessOldEnum
    774788syn keyword juliaConstC DegreeOfChannelizationEnum
    775789syn keyword juliaConstC DepthBelowSurfaceEnum
    … …  
    11001114syn keyword juliaConstC SmbMeltEnum
    11011115syn keyword juliaConstC SmbMonthlytemperaturesEnum
     1116syn keyword juliaConstC SmbMonthlydsradiationEnum
     1117syn keyword juliaConstC SmbMonthlydlradiationEnum
     1118syn keyword juliaConstC SmbMonthlywindspeedEnum
     1119syn keyword juliaConstC SmbMonthlyairhumidityEnum
    11021120syn keyword juliaConstC SmbMSurfEnum
    11031121syn keyword juliaConstC SmbNetLWEnum
    … …  
    11091127syn keyword juliaConstC SmbPrecipitationEnum
    11101128syn keyword juliaConstC SmbPrecipitationsAnomalyEnum
     1129syn keyword juliaConstC SmbDsradiationAnomalyEnum
     1130syn keyword juliaConstC SmbDlradiationAnomalyEnum
     1131syn keyword juliaConstC SmbWindspeedAnomalyEnum
     1132syn keyword juliaConstC SmbAirhumidityAnomalyEnum
    11111133syn keyword juliaConstC SmbPrecipitationsLgmEnum
    11121134syn keyword juliaConstC SmbPrecipitationsPresentdayEnum
    … …  
    11281150syn keyword juliaConstC SmbSmbrefEnum
    11291151syn keyword juliaConstC SmbSzaValueEnum
     1152syn keyword juliaConstC SmbSummerMeltEnum
     1153syn keyword juliaConstC SmbSummerAlbedoEnum
     1154syn keyword juliaConstC SmbSnowheightEnum
    11301155syn keyword juliaConstC SmbTEnum
    11311156syn keyword juliaConstC SmbTaEnum
    … …  
    16471672syn keyword juliaConstC SMBarmaEnum
    16481673syn keyword juliaConstC SMBcomponentsEnum
    1649 syn keyword juliaConstC SMBdebrisMLEnum
     1674syn keyword juliaConstC SMBdebrisEvattEnum
    16501675syn keyword juliaConstC SMBd18opddEnum
    16511676syn keyword juliaConstC SMBforcingEnum

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

    r27829 r27847  
    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;
    183184              else if (strcmp(name,"DebugProfiling")==0) return DebugProfilingEnum;
    184185              else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum;
    … …  
    259260              else if (strcmp(name,"Hydrologyarmadatebreaks")==0) return HydrologyarmadatebreaksEnum;
    260261              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,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum;
     265              if (strcmp(name,"HydrologyarmamaOrder")==0) return HydrologyarmamaOrderEnum;
     266              else if (strcmp(name,"HydrologyarmaMonthlyFactors")==0) return HydrologyarmaMonthlyFactorsEnum;
    266267              else if (strcmp(name,"HydrologyarmaNumBreaks")==0) return HydrologyarmaNumBreaksEnum;
    267268              else if (strcmp(name,"HydrologyarmaNumParams")==0) return HydrologyarmaNumParamsEnum;
    … …  
    382383              else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum;
    383384              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,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
     388              if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
     389              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
    389390              else if (strcmp(name,"MaterialsBeta")==0) return MaterialsBetaEnum;
    390391              else if (strcmp(name,"MaterialsEarthDensity")==0) return MaterialsEarthDensityEnum;
    … …  
    505506              else if (strcmp(name,"SolidearthSettingsRunFrequency")==0) return SolidearthSettingsRunFrequencyEnum;
    506507              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,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
     511              if (strcmp(name,"SolidearthSettingsHoriz")==0) return SolidearthSettingsHorizEnum;
     512              else if (strcmp(name,"SolidearthSettingsMaxiter")==0) return SolidearthSettingsMaxiterEnum;
    512513              else if (strcmp(name,"SolidearthSettingsGrdOcean")==0) return SolidearthSettingsGrdOceanEnum;
    513514              else if (strcmp(name,"SolidearthSettingsOceanAreaScaling")==0) return SolidearthSettingsOceanAreaScalingEnum;
    … …  
    553554              else if (strcmp(name,"SmbARMAmaOrder")==0) return SmbARMAmaOrderEnum;
    554555              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;
    555561              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
    556562              else if (strcmp(name,"SmbDesfacElev")==0) return SmbDesfacElevEnum;
    … …  
    566572              else if (strcmp(name,"SmbEIdx")==0) return SmbEIdxEnum;
    567573              else if (strcmp(name,"SmbF")==0) return SmbFEnum;
     574              else if (strcmp(name,"SmbHumiditylaps")==0) return SmbHumiditylapsEnum;
    568575              else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum;
     576              else if (strcmp(name,"SmbIcealbedo")==0) return SmbIcealbedoEnum;
    569577              else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum;
    570578              else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
    … …  
    586594              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    587595              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;
    588600              else if (strcmp(name,"SmbNumBasins")==0) return SmbNumBasinsEnum;
    589601              else if (strcmp(name,"SmbNumBreaks")==0) return SmbNumBreaksEnum;
    … …  
    617629              else if (strcmp(name,"SmbSemicTmax")==0) return SmbSemicTmaxEnum;
    618630              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
    619               else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
    620635              else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
    621636              else if (strcmp(name,"SmbT0wet")==0) return SmbT0wetEnum;
    … …  
    626641              else if (strcmp(name,"SmbTemperaturesReconstructedYears")==0) return SmbTemperaturesReconstructedYearsEnum;
    627642              else if (strcmp(name,"SmbPrecipitationsReconstructedYears")==0) return SmbPrecipitationsReconstructedYearsEnum;
     643              else if (strcmp(name,"SmbWindspeedlaps")==0) return SmbWindspeedlapsEnum;
    628644              else if (strcmp(name,"SmoothThicknessMultiplier")==0) return SmoothThicknessMultiplierEnum;
    629645              else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
    630646              else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
     647              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
    635648              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
    636649              else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
    … …  
    739752              else if (strcmp(name,"BasalforcingsSpatialDeepwaterElevation")==0) return BasalforcingsSpatialDeepwaterElevationEnum;
    740753              else if (strcmp(name,"BasalforcingsSpatialDeepwaterMeltingRate")==0) return BasalforcingsSpatialDeepwaterMeltingRateEnum;
    741               else if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"BasalforcingsSpatialUpperwaterElevation")==0) return BasalforcingsSpatialUpperwaterElevationEnum;
    742758              else if (strcmp(name,"BasalforcingsSpatialUpperwaterMeltingRate")==0) return BasalforcingsSpatialUpperwaterMeltingRateEnum;
    743759              else if (strcmp(name,"BasalforcingsIsmip6BasinId")==0) return BasalforcingsIsmip6BasinIdEnum;
    … …  
    752768              else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
    753769              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     770              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    758771              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    759772              else if (strcmp(name,"BasalStressx")==0) return BasalStressxEnum;
    … …  
    799812              else if (strcmp(name,"DamageF")==0) return DamageFEnum;
    800813              else if (strcmp(name,"DebrisThickness")==0) return DebrisThicknessEnum;
     814              else if (strcmp(name,"DebrisThicknessOld")==0) return DebrisThicknessOldEnum;
    801815              else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
    802816              else if (strcmp(name,"DepthBelowSurface")==0) return DepthBelowSurfaceEnum;
    … …  
    861875              else if (strcmp(name,"FrictionPressureAdjustedTemperature")==0) return FrictionPressureAdjustedTemperatureEnum;
    862876              else if (strcmp(name,"FrictionQ")==0) return FrictionQEnum;
    863               else if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"FrictionSedimentCompressibilityCoefficient")==0) return FrictionSedimentCompressibilityCoefficientEnum;
    864881              else if (strcmp(name,"FrictionTillFrictionAngle")==0) return FrictionTillFrictionAngleEnum;
    865882              else if (strcmp(name,"FrictionWaterLayer")==0) return FrictionWaterLayerEnum;
    … …  
    875892              else if (strcmp(name,"UGia")==0) return UGiaEnum;
    876893              else if (strcmp(name,"UGiaRate")==0) return UGiaRateEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"Gradient")==0) return GradientEnum;
     894              else if (strcmp(name,"Gradient")==0) return GradientEnum;
    881895              else if (strcmp(name,"GroundinglineHeight")==0) return GroundinglineHeightEnum;
    882896              else if (strcmp(name,"HydraulicPotential")==0) return HydraulicPotentialEnum;
    … …  
    984998              else if (strcmp(name,"SamplingTau")==0) return SamplingTauEnum;
    985999              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    986               else if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"SealevelGRD")==0) return SealevelGRDEnum;
    9871004              else if (strcmp(name,"SatGraviGRD")==0) return SatGraviGRDEnum;
    9881005              else if (strcmp(name,"SealevelBarystaticMask")==0) return SealevelBarystaticMaskEnum;
    … …  
    9981015              else if (strcmp(name,"SealevelBarystaticHydroLatbar")==0) return SealevelBarystaticHydroLatbarEnum;
    9991016              else if (strcmp(name,"SealevelBarystaticHydroLongbar")==0) return SealevelBarystaticHydroLongbarEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
     1017              else if (strcmp(name,"SealevelBarystaticHydroLoad")==0) return SealevelBarystaticHydroLoadEnum;
    10041018              else if (strcmp(name,"SealevelBarystaticBpMask")==0) return SealevelBarystaticBpMaskEnum;
    10051019              else if (strcmp(name,"SealevelBarystaticBpWeights")==0) return SealevelBarystaticBpWeightsEnum;
    … …  
    11071121              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    11081122              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    1109               else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    11101127              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    11111128              else if (strcmp(name,"SmbFAC")==0) return SmbFACEnum;
    … …  
    11211138              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    11221139              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
     1140              else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    11271141              else if (strcmp(name,"SmbMassBalanceSnow")==0) return SmbMassBalanceSnowEnum;
    11281142              else if (strcmp(name,"SmbMassBalanceIce")==0) return SmbMassBalanceIceEnum;
    … …  
    11361150              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    11371151              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;
    11381156              else if (strcmp(name,"SmbMSurf")==0) return SmbMSurfEnum;
    11391157              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
    … …  
    11451163              else if (strcmp(name,"SmbPrecipitation")==0) return SmbPrecipitationEnum;
    11461164              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;
    11471169              else if (strcmp(name,"SmbPrecipitationsLgm")==0) return SmbPrecipitationsLgmEnum;
    11481170              else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
    … …  
    11641186              else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum;
    11651187              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;
    11661191              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    11671192              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
    … …  
    12191244              else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
    12201245              else if (strcmp(name,"Surface")==0) return SurfaceEnum;
    1221               else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
    12221250              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
    12231251              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
    … …  
    12441272              else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
    12451273              else if (strcmp(name,"TransientAccumulatedDeltaIceThickness")==0) return TransientAccumulatedDeltaIceThicknessEnum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Vel")==0) return VelEnum;
     1274              else if (strcmp(name,"Vel")==0) return VelEnum;
    12501275              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    12511276              else if (strcmp(name,"VxBase")==0) return VxBaseEnum;
    … …  
    13421367              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    13431368              else if (strcmp(name,"Outputdefinition60")==0) return Outputdefinition60Enum;
    1344               else if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"Outputdefinition61")==0) return Outputdefinition61Enum;
    13451373              else if (strcmp(name,"Outputdefinition62")==0) return Outputdefinition62Enum;
    13461374              else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum;
    … …  
    13671395              else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum;
    13681396              else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
     1397              else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum;
    13731398              else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
    13741399              else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
    … …  
    14651490              else if (strcmp(name,"DefaultAnalysis")==0) return DefaultAnalysisEnum;
    14661491              else if (strcmp(name,"DefaultCalving")==0) return DefaultCalvingEnum;
    1467               else if (strcmp(name,"Dense")==0) return DenseEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"Dense")==0) return DenseEnum;
    14681496              else if (strcmp(name,"DependentObject")==0) return DependentObjectEnum;
    14691497              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;
    … …  
    14901518              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    14911519              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
     1520              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    14961521              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    14971522              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    … …  
    15881613              else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum;
    15891614              else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum;
    1590               else if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"LoveKf")==0) return LoveKfEnum;
    15911619              else if (strcmp(name,"LoveKt")==0) return LoveKtEnum;
    15921620              else if (strcmp(name,"LoveLf")==0) return LoveLfEnum;
    … …  
    16131641              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    16141642              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    1615          else stage=14;
    1616    }
    1617    if(stage==14){
    1618               if (strcmp(name,"Matice")==0) return MaticeEnum;
     1643              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    16191644              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    16201645              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    … …  
    16951720              else if (strcmp(name,"SMBarma")==0) return SMBarmaEnum;
    16961721              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    1697               else if (strcmp(name,"SMBdebrisML")==0) return SMBdebrisMLEnum;
     1722              else if (strcmp(name,"SMBdebrisEvatt")==0) return SMBdebrisEvattEnum;
    16981723              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    16991724              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    … …  
    17111736              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    17121737              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    1713               else if (strcmp(name,"Scaled")==0) return ScaledEnum;
     1738         else stage=15;
     1739   }
     1740   if(stage==15){
     1741              if (strcmp(name,"Scaled")==0) return ScaledEnum;
    17141742              else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    17151743              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    … …  
    17361764              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    17371765              else if (strcmp(name,"Sset")==0) return SsetEnum;
    1738          else stage=15;
    1739    }
    1740    if(stage==15){
    1741               if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
     1766              else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
    17421767              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    17431768              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;

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

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

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

    r27417 r27847  
    1212                removalmodel             = 0;
    1313                displacementmodel        = 0;
     14                max_displacementvelocity = 0;
    1415                removal_slope_threshold  = 0;
    1516                removal_stress_threshold = 0;
    … …  
    6667                function list = defaultoutputs(self,md) % {{{
    6768
    68                         list = {'DebrisThickness'};
     69                        list = {'DebrisThickness','DebrisMaskNodeActivation','VxDebris','VyDebris'};
    6970
    7071                end % }}}
    … …  
    9293                        self.removal_stress_threshold=0;
    9394
     95                        %Max velocity for displacementmodel (1)
     96                        self.max_displacementvelocity=0;
     97
    9498                        %default output
    9599                        self.requested_outputs={'default'};
    … …  
    101105
    102106                        md = checkfield(md,'fieldname','debris.spcthickness');
    103                         md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3]);
     107                        md = checkfield(md,'fieldname','debris.stabilization','values',[0 1 2 3 4 5]);
    104108                        md = checkfield(md,'fieldname','debris.min_thickness','>=',0);
    105109                        md = checkfield(md,'fieldname','debris.packingfraction','>=',0);
    … …  
    108112                        md = checkfield(md,'fieldname','debris.removal_slope_threshold','>=',0);
    109113                        md = checkfield(md,'fieldname','debris.removal_stress_threshold','>=',0);
     114                        md = checkfield(md,'fieldname','debris.max_displacementvelocity','>=',0);
    110115                        md = checkfield(md,'fieldname','debris.requested_outputs','stringrow',1);
    111116                        if ~any(isnan(md.stressbalance.vertex_pairing)),
    … …  
    120125                        fielddisplay(self,'stabilization','0: no stabilization, 1: artificial diffusion, 2: streamline upwinding, 3: streamline upwind Petrov-Galerkin (SUPG)');
    121126                        fielddisplay(self,'removalmodel','frontal removal of debris. 0: no removal, 1: Slope-triggered debris removal, 2: driving-stress triggered debris removal');
    122                         fielddisplay(self,'displacementmodel','debris displacement. 0: no displacement, 1: ...');
     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');
    123129                        fielddisplay(self,'removal_slope_threshold','critical slope (degrees) for removalmodel (1)');
    124130                        fielddisplay(self,'removal_stress_threshold','critical stress (Pa) for removalmodel (2)');
    … …  
    135141                        WriteData(fid,prefix,'object',self,'fieldname','removalmodel','format','Integer');
    136142                        WriteData(fid,prefix,'object',self,'fieldname','displacementmodel','format','Integer');
     143                        WriteData(fid,prefix,'object',self,'fieldname','max_displacementvelocity','format','Double');
    137144                        WriteData(fid,prefix,'object',self,'fieldname','removal_slope_threshold','format','Double');
    138145                        WriteData(fid,prefix,'object',self,'fieldname','removal_stress_threshold','format','Double');
    … …  
    156163                        writejsdouble(fid,[modelname '.debris.removalmodel'],self.removalmodel);
    157164                        writejsdouble(fid,[modelname '.debris.displacementmodel'],self.displacementmodel);
     165                        writejsdouble(fid,[modelname '.debris.max_displacementvelocity'],self.displacementmodel);
    158166                        writejsdouble(fid,[modelname '.debris.removal_slope_threshold'],self.removal_slope_threshold);
    159167                        writejsdouble(fid,[modelname '.debris.removal_stress_threshold'],self.removal_stress_threshold);

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

    r27648 r27847  
    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
    12911292                        if isfield(structmd,'artificial_diffusivity'), md.masstransport.stabilization=structmd.artificial_diffusivity; end
    12921293                        if isfield(structmd,'hydrostatic_adjustment'), md.masstransport.hydrostatic_adjustment=structmd.hydrostatic_adjustment; end

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

    r27829 r27847  
    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'),
    242246                field = field*yts;
    243247        elseif strcmp(fieldname,'BasalforcingsGroundediceMeltingRate'),
Note: See TracChangeset for help on using the changeset viewer.