Changeset 21740


Ignore:
Timestamp:
05/23/17 00:22:34 (8 years ago)
Author:
sjohnsen
Message:

NEW: adding a friction effective pressure flag

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

Legend:

Unmodified
Added
Removed

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

    r21722 r21740  
    174174        switch(frictionlaw){
    175175                case 1:
     176                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
    176177                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
    177178                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
    178179                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
     180                        if (FrictionCoupling==1){
     181                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     182                        }
    179183                        break;
    180184                case 2:
    … …  
    187191                        iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
    188192                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
    189                         if (FrictionCoupling==0){
     193                        if (FrictionCoupling==1){
    190194                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
    191195                        }
    … …  
    245249        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    246250        if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    247         if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     251        if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
    248252        if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    249253}/*}}}*/

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

    r21600 r21740  
    781781        switch(frictionlaw){
    782782                case 1:
     783                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
    783784                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
    784785                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
    785786                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
     787                        if(FrictionCoupling==1){
     788                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     789                        }
    786790                        break;
    787791                case 2:
    … …  
    794798                        iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
    795799                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
    796                         if(FrictionCoupling==0){
     800                        if(FrictionCoupling==1){
    797801                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
    798802                        }
    … …  
    893897        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    894898        if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    895         if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     899        if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
    896900        if(frictionlaw==5) parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum));
    897901        if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));

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

    r21542 r21740  
    157157        switch(frictionlaw){
    158158                case 1:
     159                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
    159160                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
    160161                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
    161162                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
     163                        if (FrictionCoupling==1){
     164                          iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     165                        }
    162166                        break;
    163167                case 2:
    … …  
    170174                        iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
    171175                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
    172                         if (FrictionCoupling==0){
    173                                 iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     176                        if (FrictionCoupling==1){
     177                          iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
    174178                        }
    175179                        break;
    … …  
    220224        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    221225        if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    222         if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     226        if(frictionlaw==3 || frictionlaw==1) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     227
    223228}/*}}}*/
    224229

  • issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

    r21600 r21740  
    7070        IssmDouble  vx,vy,vz,vmag;
    7171        IssmDouble  alpha_complement;
     72        IssmDouble  base,sealevel,thickness;
    7273
    7374        /*Recover parameters: */
    … …  
    7778        element->GetInputValue(&As,gauss,FrictionAsEnum);
    7879        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
     80        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     81        element->GetInputValue(&base, gauss,BaseEnum);
     82        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     83
     84        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     85        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     86        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    7987        element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    8088
    81         if (CoupledFlag==1){
    82                 element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
    83         }
    84         else{
    85                 element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     89        switch(CoupledFlag){
     90                case 0:
     91                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     92                        break;
     93                case 1:
     94                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     95                        break;
     96                case 2:
     97                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     98                        break;
     99                default:
     100                        _error_("not supported");
    86101        }
    87102
    … …  
    159174void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
    160175
    161         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     176        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p.
    162177         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    163178         * alpha_complement= Neff ^r * vel ^s*/
    164179
    165180        /*diverse: */
     181        int         CoupledFlag;
    166182        IssmDouble  r,s;
    167183        IssmDouble  vx,vy,vz,vmag;
    … …  
    169185        IssmDouble  Neff;
    170186        IssmDouble  drag_coefficient;
    171         IssmDouble  bed,thickness,sealevel;
     187        IssmDouble  base,thickness,sealevel;
    172188        IssmDouble  alpha_complement;
    173 
    174         /*Recover parameters: */
    175         element->GetInputValue(&drag_p,FrictionPEnum);
    176         element->GetInputValue(&drag_q,FrictionQEnum);
    177         element->GetInputValue(&thickness, gauss,ThicknessEnum);
    178         element->GetInputValue(&bed, gauss,BaseEnum);
    179         element->GetInputValue(&sealevel, gauss,SealevelEnum);
    180         element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    181         IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    182         IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    183         IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    184 
    185         //compute r and q coefficients: */
    186         r=drag_q/drag_p;
    187         s=1./drag_p;
    188 
    189         //From bed and thickness, compute effective pressure when drag is viscous:
    190         Neff=gravity*(rho_ice*thickness+rho_water*(bed-sealevel));
    191         if(Neff<0)Neff=0;
    192 
    193         //We need the velocity magnitude to evaluate the basal stress:
    194         switch(dim){
    195                 case 1:
    196                         element->GetInputValue(&vx,gauss,VxEnum);
    197                         vmag=sqrt(vx*vx);
    198                         break;
    199                 case 2:
    200                         element->GetInputValue(&vx,gauss,VxEnum);
    201                         element->GetInputValue(&vy,gauss,VyEnum);
    202                         vmag=sqrt(vx*vx+vy*vy);
    203                         break;
    204                 case 3:
    205                         element->GetInputValue(&vx,gauss,VxEnum);
    206                         element->GetInputValue(&vy,gauss,VyEnum);
    207                         element->GetInputValue(&vz,gauss,VzEnum);
    208                         vmag=sqrt(vx*vx+vy*vy+vz*vz);
    209                         break;
    210                 default:
    211                         _error_("not supported");
    212         }
    213 
    214         /*Check to prevent dividing by zero if vmag==0*/
    215         if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
    216         else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
    217 
    218         /*Assign output pointers:*/
    219         *palpha_complement=alpha_complement;
    220 }
    221 /*}}}*/
    222 void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    223 
    224         switch(this->law){
    225                 case 1:
    226                         GetAlpha2Viscous(palpha2,gauss);
    227                         break;
    228                 case 2:
    229                         GetAlpha2Weertman(palpha2,gauss);
    230                         break;
    231                 case 3:
    232                         GetAlpha2Hydro(palpha2,gauss);
    233                         break;
    234                 case 4:
    235                         GetAlpha2Temp(palpha2,gauss);
    236                         break;
    237                 case 5:
    238                         GetAlpha2WaterLayer(palpha2,gauss);
    239                         break;
    240                 case 6:
    241                         GetAlpha2WeertmanTemp(palpha2,gauss);
    242                         break;
    243                 case 7:
    244                         GetAlpha2Coulomb(palpha2,gauss);
    245                         break;
    246                 case 8:
    247                         GetAlpha2Sommers(palpha2,gauss);
    248                         break;
    249                 case 9:
    250                         GetAlpha2Josh(palpha2,gauss);
    251                         break;
    252           default:
    253                         _error_("Friction law "<< this->law <<" not supported");
    254         }
    255 
    256 }/*}}}*/
    257 void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    258 
    259         /*This routine calculates the basal friction coefficient
    260           alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    261 
    262         /*diverse: */
    263         IssmDouble  r,s;
    264         IssmDouble  drag_p, drag_q;
    265         IssmDouble  Neff;
    266         IssmDouble  thickness,base,bed,floatation_thickness,sealevel;
    267         IssmDouble  vx,vy,vz,vmag;
    268         IssmDouble  drag_coefficient,drag_coefficient_coulomb;
    269         IssmDouble  alpha2,alpha2_coulomb;
    270189
    271190        /*Recover parameters: */
    … …  
    275194        element->GetInputValue(&base, gauss,BaseEnum);
    276195        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    277         element->GetInputValue(&bed, gauss,BedEnum);
     196        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
     197        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     198        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     199        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
     200        element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
     201        //compute r and q coefficients: */
     202        r=drag_q/drag_p;
     203        s=1./drag_p;
     204
     205        //From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing or coupled to hydrologymodel:
     206        switch(CoupledFlag){
     207                case 0:
     208                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     209                        break; 
     210                case 1:
     211                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     212                        break;
     213                case 2:
     214                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     215                        break;
     216                default:
     217                        _error_("not supported");
     218        }
     219        if(Neff<0)Neff=0;
     220
     221        //We need the velocity magnitude to evaluate the basal stress:
     222        switch(dim){
     223                case 1:
     224                        element->GetInputValue(&vx,gauss,VxEnum);
     225                        vmag=sqrt(vx*vx);
     226                        break;
     227                case 2:
     228                        element->GetInputValue(&vx,gauss,VxEnum);
     229                        element->GetInputValue(&vy,gauss,VyEnum);
     230                        vmag=sqrt(vx*vx+vy*vy);
     231                        break;
     232                case 3:
     233                        element->GetInputValue(&vx,gauss,VxEnum);
     234                        element->GetInputValue(&vy,gauss,VyEnum);
     235                        element->GetInputValue(&vz,gauss,VzEnum);
     236                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
     237                        break;
     238                default:
     239                        _error_("not supported");
     240        }
     241
     242        /*Check to prevent dividing by zero if vmag==0*/
     243        if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
     244        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
     245
     246        /*Assign output pointers:*/
     247        *palpha_complement=alpha_complement;
     248}
     249/*}}}*/
     250void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     251
     252        switch(this->law){
     253                case 1:
     254                        GetAlpha2Viscous(palpha2,gauss);
     255                        break;
     256                case 2:
     257                        GetAlpha2Weertman(palpha2,gauss);
     258                        break;
     259                case 3:
     260                        GetAlpha2Hydro(palpha2,gauss);
     261                        break;
     262                case 4:
     263                        GetAlpha2Temp(palpha2,gauss);
     264                        break;
     265                case 5:
     266                        GetAlpha2WaterLayer(palpha2,gauss);
     267                        break;
     268                case 6:
     269                        GetAlpha2WeertmanTemp(palpha2,gauss);
     270                        break;
     271                case 7:
     272                        GetAlpha2Coulomb(palpha2,gauss);
     273                        break;
     274                case 8:
     275                        GetAlpha2Sommers(palpha2,gauss);
     276                        break;
     277                case 9:
     278                        GetAlpha2Josh(palpha2,gauss);
     279                        break;
     280          default:
     281                        _error_("Friction law "<< this->law <<" not supported");
     282        }
     283
     284}/*}}}*/
     285void Friction::GetAlpha2Coulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
     286
     287        /*This routine calculates the basal friction coefficient
     288          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
     289
     290        /*diverse: */
     291        IssmDouble  r,s;
     292        IssmDouble  drag_p, drag_q;
     293        IssmDouble  Neff;
     294        IssmDouble  thickness,base,floatation_thickness,sealevel;
     295        IssmDouble  vx,vy,vz,vmag;
     296        IssmDouble  drag_coefficient,drag_coefficient_coulomb;
     297        IssmDouble  alpha2,alpha2_coulomb;
     298
     299        /*Recover parameters: */
     300        element->GetInputValue(&drag_p,FrictionPEnum);
     301        element->GetInputValue(&drag_q,FrictionQEnum);
     302        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     303        element->GetInputValue(&base, gauss,BaseEnum);
     304        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    278305        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    279306        element->GetInputValue(&drag_coefficient_coulomb, gauss,FrictionCoefficientcoulombEnum);
    … …  
    315342
    316343        floatation_thickness=0;
    317         if(bed<0) floatation_thickness=-rho_water/rho_ice*bed;
     344        if(base<0) floatation_thickness=-rho_water/rho_ice*base;
    318345        if(vmag==0.) alpha2_coulomb=0.;
    319346        else alpha2_coulomb=drag_coefficient_coulomb*drag_coefficient_coulomb*rho_water*gravity*(thickness-floatation_thickness)/vmag;
    … …  
    349376        IssmDouble  vx,vy,vz,vmag;
    350377        IssmDouble  alpha2;
     378        IssmDouble  base,thickness,sealevel;
    351379
    352380        /*Recover parameters: */
    … …  
    355383        element->GetInputValue(&As,gauss,FrictionAsEnum);
    356384        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
     385        element->GetInputValue(&thickness, gauss,ThicknessEnum);
     386        element->GetInputValue(&base, gauss,BaseEnum);
     387        element->GetInputValue(&sealevel, gauss,SealevelEnum);
     388        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     389        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
     390        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
     391        element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    357392       
    358         element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    359         if (CoupledFlag==1){
    360                 element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
    361         }
    362         else{
    363                 element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
    364         }
    365                
     393        switch(CoupledFlag){
     394                case 0:
     395                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     396                        break; 
     397                case 1:
     398                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     399                        break;
     400                case 2:
     401                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     402                        break;
     403                default:
     404                        _error_("not supported");
     405        }
     406
    366407        if(Neff<0)Neff=0;
    367408
    … …  
    407448void Friction::GetAlpha2Sommers(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
    408449
    409         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-bed)*/
     450        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff, with Neff=rho_ice*g*thickness+rho_ice*g*(head-base)*/
    410451
    411452        /*diverse: */
    … …  
    413454        IssmDouble  Neff;
    414455        IssmDouble  drag_coefficient;
    415         IssmDouble  bed,thickness,head,sealevel;
     456        IssmDouble  base,thickness,head,sealevel;
    416457        IssmDouble  alpha2;
    417458
    418459        /*Recover parameters: */
    419460        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    420         element->GetInputValue(&bed, gauss,BaseEnum);
     461        element->GetInputValue(&base, gauss,BaseEnum);
    421462        element->GetInputValue(&head, gauss,HydrologyHeadEnum);
    422463        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    … …  
    426467        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    427468
    428         //From bed and thickness, compute effective pressure when drag is viscous:
     469        //From base and thickness, compute effective pressure when drag is viscous:
    429470        pressure_ice   = rho_ice*gravity*thickness;
    430         pressure_water = rho_water*gravity*(head-bed+sealevel);
     471        pressure_water = rho_water*gravity*(head-base+sealevel);
    431472        Neff=pressure_ice-pressure_water;
    432473        if(Neff<0.) Neff=0.;
    … …  
    511552
    512553        /*This routine calculates the basal friction coefficient
    513           alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     554          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    514555
    515556        /*diverse: */
     557        int         CoupledFlag;
    516558        IssmDouble  r,s;
    517559        IssmDouble  drag_p, drag_q;
    … …  
    532574        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
    533575        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
    534 
     576        element->parameters->FindParam(&CoupledFlag,FrictionCouplingEnum);
    535577        //compute r and q coefficients: */
    536578        r=drag_q/drag_p;
    537579        s=1./drag_p;
    538580
    539         //From base and thickness, compute effective pressure when drag is viscous:
    540         Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     581        //From base and thickness, compute effective pressure when drag is viscous, or get Neff from forcing:
     582        switch(CoupledFlag){
     583                case 0:
     584                        Neff=gravity*(rho_ice*thickness+rho_water*(base-sealevel));
     585                        break;
     586                case 1:
     587                        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
     588                        break;
     589                case 2:
     590                        element->GetInputValue(&Neff,gauss,EffectivePressureEnum);
     591                        break;
     592                default:
     593                        _error_("not supported");
     594        }
    541595        if(Neff<0)Neff=0;
    542596
    … …  
    572626
    573627        /*This routine calculates the basal friction coefficient
    574           alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
     628          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*base, r=q/p and s=1/p**/
    575629
    576630        /*diverse: */
    … …  
    578632        IssmDouble  drag_p, drag_q;
    579633        IssmDouble  Neff,F;
    580         IssmDouble  thickness,bed,sealevel;
     634        IssmDouble  thickness,base,sealevel;
    581635        IssmDouble  vx,vy,vz,vmag;
    582636        IssmDouble  drag_coefficient,water_layer;
    … …  
    588642        element->GetInputValue(&drag_q,FrictionQEnum);
    589643        element->GetInputValue(&thickness, gauss,ThicknessEnum);
    590         element->GetInputValue(&bed, gauss,BaseEnum);
     644        element->GetInputValue(&base, gauss,BaseEnum);
    591645        element->GetInputValue(&sealevel, gauss,SealevelEnum);
    592646        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
    … …  
    600654        s=1./drag_p;
    601655
    602         //From bed and thickness, compute effective pressure when drag is viscous:
    603         if(bed>0) bed=0;
    604         if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(bed-sealevel);
     656        //From base and thickness, compute effective pressure when drag is viscous:
     657        if(base>0) base=0;
     658        if(water_layer==0) Neff=gravity*rho_ice*thickness+gravity*rho_water*(base-sealevel);
    605659        else if(water_layer>0) Neff=gravity*rho_ice*thickness*F;
    606660        else _error_("negative water layer thickness");
Note: See TracChangeset for help on using the changeset viewer.