Changeset 16816


Ignore:
Timestamp:
11/18/13 00:55:32 (11 years ago)
Author:
jbondzio
Message:

CHG: Drainage of waterfraction debugged. Drainage happens now in Penta::ComputeBasalMeltingrates()

Location:
issm/trunk-jpl
Files:
14 edited

Legend:

Unmodified
Added
Removed

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

    r16813 r16816  
    192192                virtual void UpdateBasalConstraintsEnthalpy(void)=0;
    193193                virtual void ComputeBasalMeltingrate(void)=0;
    194                 virtual void DrainWaterfraction(void)=0;
     194                virtual void DrainWaterfraction(IssmDouble* drainrate_element)=0;
    195195                #endif
    196196

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

    r16813 r16816  
    47504750        /* Intermediaries */
    47514751        bool        isenthalpy, checkpositivethickness, istemperatelayer;
    4752         int         i,j,analysis_type;
     4752        int         i,j,iv,analysis_type;
    47534753        IssmDouble  xyz_list[NUMVERTICES][3];
    47544754        IssmDouble  xyz_list_tria[NUMVERTICES2D][3];
    … …  
    47644764        IssmDouble  geothermalflux[NUMVERTICES];
    47654765        IssmDouble  dt, yts;
    4766         IssmDouble  melting_overshoot,meltingrate_enthalpy;
    4767         IssmDouble  lambda,heating;
     4766        IssmDouble  melting_overshoot,meltingrate_enthalpy[NUMVERTICES2D];
     4767        IssmDouble  drainrate_element[NUMVERTICES2D],drainrate_column[NUMVERTICES2D];
     4768        IssmDouble  lambda,heating[NUMVERTICES2D];
    47684769        Friction   *friction  = NULL;
     4770        Penta      *penta = NULL;
    47694771
    47704772        /* Only compute melt rates at the base of grounded ice*/
    … …  
    47954797        friction=new Friction("3d",inputs,matpar,analysis_type);
    47964798
    4797         /*Ok, get meltingrates now from basal conditions*/
     4799        /******** MELTING RATES  ************************************/
    47984800        GaussPenta* gauss=new GaussPenta();
    4799         for(int iv=0;iv<3;iv++){
     4801        for(iv=0;iv<NUMVERTICES2D;iv++){
    48004802
    48014803                gauss->GaussVertex(iv);
    … …  
    48054807
    48064808                /*Calculate basal meltingrate after Fig.5 of A.Aschwanden 2012*/
    4807                 meltingrate_enthalpy=0.;
    4808                 heating=0.;
     4809                meltingrate_enthalpy[iv]=0.;
     4810                heating[iv]=0.;
    48094811                if((watercolumn[iv]>0.) && (enthalpy[iv]<matpar->PureIceEnthalpy(pressure[iv]))){
    48104812                        /*ensure that no ice is at T<Tm(p), if water layer present*/
    … …  
    48204822                        /*ok, from here on all basal ice is temperate. Check for positive thickness of layer of temperate ice. */
    48214823                        istemperatelayer=false;
    4822                         if(enthalpy[iv+3]>=matpar->PureIceEnthalpy(pressure[iv+3])) istemperatelayer=true;
     4824                        if(enthalpy[iv+NUMVERTICES2D]>=matpar->PureIceEnthalpy(pressure[iv+NUMVERTICES2D])) istemperatelayer=true;
    48234825                        if(istemperatelayer) for(i=0;i<3;i++) vec_heatflux[i]=0.; // TODO: add -k*nabla T_pmp
    48244826                        else{
    … …  
    48394841                        matpar->EnthalpyToThermal(&temperature, &waterfraction, enthalpy[iv],pressure[iv]);
    48404842                        // -Mb= Fb-(q-q_geo)/((1-w)*L), cf Aschwanden 2012, eq.66
    4841                         heating=(heatflux+basalfriction+geothermalflux[iv]);
    4842                         meltingrate_enthalpy=heating/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent
    4843                 }
    4844 
    4845                 /*Update water column, basal meltingrate*/
    4846                 this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4843                        heating[iv]=(heatflux+basalfriction+geothermalflux[iv]);
     4844                        meltingrate_enthalpy[iv]=heating[iv]/((1-waterfraction)*latentheat*rho_ice); // m/s water equivalent
     4845                }               
     4846        }
     4847        // enthalpy might have been changed, update
     4848        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
     4849
     4850        /******** DRAINAGE *****************************************/
     4851        for(iv=0; iv<NUMVERTICES2D; iv++)       
     4852                drainrate_column[iv]=0.;
     4853        penta=this;
     4854
     4855        for(;;){
     4856                for(iv=0; iv<NUMVERTICES2D; iv++)       drainrate_element[iv]=0.;
     4857                penta->DrainWaterfraction(&drainrate_element[0]); // TODO: make sure every vertex is only drained once
     4858                for(iv=0; iv<NUMVERTICES2D; iv++)       drainrate_column[iv]+=drainrate_element[iv];
     4859
     4860                if(penta->IsOnSurface()) break;
     4861                penta=penta->GetUpperElement();                 
     4862        }
     4863        // add drained water to melting rate
     4864        for(iv=0; iv<NUMVERTICES2D;iv++)
     4865                meltingrate_enthalpy[iv]+=drainrate_column[iv];
     4866
     4867        /******** UPDATE MELTINGRATES AND WATERCOLUMN **************/
     4868        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4869        for(iv=0;iv<NUMVERTICES2D;iv++){
    48474870                if(reCast<bool,IssmDouble>(dt)){
    4848                         if(watercolumn[iv]+meltingrate_enthalpy*dt<0.){                         
    4849                                 melting_overshoot=watercolumn[iv]+meltingrate_enthalpy*dt;
    4850                                 lambda=melting_overshoot/(meltingrate_enthalpy*dt); _assert_(lambda>0); _assert_(lambda<1);
    4851                                 basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy;
     4871                        if(watercolumn[iv]+meltingrate_enthalpy[iv]*dt<0.){                             
     4872                                melting_overshoot=watercolumn[iv]+meltingrate_enthalpy[iv]*dt;
     4873                                lambda=melting_overshoot/(meltingrate_enthalpy[iv]*dt); _assert_(lambda>0); _assert_(lambda<1);
     4874                                basalmeltingrate[iv]=(1.-lambda)*meltingrate_enthalpy[iv];
    48524875                                watercolumn[iv]=0.;
    48534876                                yts=365*24*60*60;
    4854                                 enthalpy[iv]+=dt/yts*lambda*heating;
     4877                                enthalpy[iv]+=dt/yts*lambda*heating[iv];
    48554878                        }
    48564879                        else{
    4857                                 basalmeltingrate[iv]=meltingrate_enthalpy;
    4858                                 watercolumn[iv]+=dt*meltingrate_enthalpy;
     4880                                basalmeltingrate[iv]=meltingrate_enthalpy[iv];
     4881                                watercolumn[iv]+=dt*meltingrate_enthalpy[iv];
    48594882                        }
    48604883                }
    48614884                else{
    4862                         basalmeltingrate[iv]=meltingrate_enthalpy;
    4863                         watercolumn[iv]+=meltingrate_enthalpy;
     4885                        basalmeltingrate[iv]=meltingrate_enthalpy[iv];
     4886                        watercolumn[iv]+=meltingrate_enthalpy[iv];
    48644887                }         
    4865         } 
     4888               
     4889                _assert_(watercolumn[iv]>=0.);
     4890        }
     4891
    48664892        /*feed updated variables back into model*/
    48674893        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    … …  
    48754901/*}}}*/
    48764902/*FUNCTION Penta::DrainWaterfraction{{{*/
    4877 void Penta::DrainWaterfraction(void){
     4903void Penta::DrainWaterfraction(IssmDouble* drainrate_element){
    48784904
    48794905    /*Intermediaries*/
    48804906        bool isenthalpy;
     4907        int iv, index0;
    48814908        IssmDouble waterfraction[NUMVERTICES], temperature[NUMVERTICES];
    4882         IssmDouble watercolumnbase[NUMVERTICES];
     4909        IssmDouble dw[NUMVERTICES];
    48834910        IssmDouble enthalpy[NUMVERTICES], pressure[NUMVERTICES];
    4884         IssmDouble latentheat, dt;
    4885         IssmDouble dw, dwc;
    4886         Penta *pentabase = NULL;
     4911        IssmDouble dt, height_element;
     4912        IssmDouble xyz_list[NUMVERTICES][3];
     4913        IssmDouble rho_water, rho_ice;
    48874914
    48884915        /*Check wether enthalpy is activated*/
    … …  
    48904917        if(!isenthalpy) return;       
    48914918       
    4892         /*get basal element, needed for basal watercolumn*/
    4893         pentabase=(Penta*)this->GetBasalElement();
    4894        
     4919        rho_ice=matpar->GetRhoIce();
     4920        rho_water=matpar->GetRhoFreshwater();
     4921
     4922        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    48954923        this->GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
    48964924        this->GetInputListOnVertices(&pressure[0],PressureEnum);
    4897         pentabase->GetInputListOnVertices(&watercolumnbase[0], WatercolumnEnum);
    48984925
    48994926        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4900         latentheat=matpar->GetLatentHeat();
    4901 
    4902         for(int iv=0;iv<NUMVERTICES;iv++){
     4927        for(iv=0;iv<NUMVERTICES;iv++){
    49034928                matpar->EnthalpyToThermal(&temperature[iv],&waterfraction[iv], enthalpy[iv],pressure[iv]);
    4904                 dw=DrainageFunctionWaterfraction(waterfraction[iv], dt);
    4905                 /*drain water fraction & update enthalpy*/
    4906                 waterfraction[iv]-=dw;
    4907                 matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);       
    4908                 /*add drained water to watercolumn at base*/
    4909                 dwc=dw*this->IceVolume();
    4910                 watercolumnbase[iv%NUMVERTICES2D]+=dwc;
    4911         }
    4912 
    4913         /*feed updated results back into model*/
     4929                dw[iv]=DrainageFunctionWaterfraction(waterfraction[iv], dt);
     4930        }
     4931       
     4932        /*drain waterfraction, feed updated variables back into model*/
     4933        for(iv=0;iv<NUMVERTICES;iv++){
     4934                if(reCast<bool,IssmDouble>(dt))
     4935                        waterfraction[iv]-=dw[iv]*dt;
     4936                else
     4937                        waterfraction[iv]-=dw[iv];
     4938                matpar->ThermalToEnthalpy(&enthalpy[iv], temperature[iv], waterfraction[iv], pressure[iv]);
     4939        }
    49144940        this->inputs->AddInput(new PentaInput(EnthalpyEnum,enthalpy,P1Enum));
    49154941        this->inputs->AddInput(new PentaInput(WaterfractionEnum,waterfraction,P1Enum));
    4916         pentabase->inputs->AddInput(new PentaInput(WatercolumnEnum, watercolumnbase,P1Enum));
    4917 
     4942
     4943        /*return meltwater column equivalent to drained water*/
     4944        for(iv=0;iv<NUMVERTICES2D;iv++){
     4945                index0=(iv+NUMVERTICES2D);
     4946                height_element=fabs(xyz_list[index0][2]-xyz_list[iv][2]);
     4947                drainrate_element[iv]=(dw[iv]+dw[index0])/2.*rho_water/rho_ice*height_element;
     4948        }
    49184949}
    49194950/*}}}*/

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

    r16813 r16816  
    378378                void           UpdateBasalConstraintsEnthalpy(void);
    379379                void           ComputeBasalMeltingrate(void);
    380                 void           DrainWaterfraction(void);
     380                void           DrainWaterfraction(IssmDouble* drainrate_element);
    381381                #endif
    382382

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

    r16813 r16816  
    142142                void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
    143143                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
    144                 void DrainWaterfraction(void){_error_("not implemented yet");};
     144                void DrainWaterfraction(IssmDouble* drainrate_element){_error_("not implemented yet");};
    145145                #endif
    146146                #ifdef _HAVE_HYDROLOGY_

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

    r16813 r16816  
    327327                void UpdateBasalConstraintsEnthalpy(void){_error_("not implemented yet");};
    328328                void ComputeBasalMeltingrate(void){_error_("not implemented yet");};
    329                 void DrainWaterfraction(void){_error_("not implemented yet");};
     329                void DrainWaterfraction(IssmDouble* drainrate_element){_error_("not implemented yet");};
    330330                #endif
    331331

  • issm/trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp

    r16360 r16816  
    1414
    1515        IssmDouble w0=0.01, w1=0.02, w2=0.03;
    16         IssmDouble Dret, D0=0, D1=0.005, D2=0.05;
    17         IssmDouble yts=365*24*60*60;
    18         dt/=yts;
     16        IssmDouble yts=365.*24.*60.*60.;
     17        IssmDouble Dret, D0=0., D1=0.005/yts, D2=0.05/yts;
    1918
    2019        /*get drainage function value*/
    2120        if((w0==w1)||(w1==w2)||(w0==w2))
    2221                _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
     22       
    2323        if(waterfraction<=w0)
    2424                Dret=D0;
    25         if((waterfraction>w0) && (waterfraction<=w1))
     25        else if((waterfraction>w0) && (waterfraction<=w1))
    2626                Dret=(D1-D0)/(w1-w0)*(waterfraction-w0)+D0;
    27         if((waterfraction>w1) && (waterfraction<=w2))
     27        else if((waterfraction>w1) && (waterfraction<=w2))
    2828                Dret=(D2-D1)/(w2-w1)*(waterfraction-w1)+D1;
    2929        else
    3030                Dret=D2;
    3131
    32         /*check if dt*Dret>waterfraction. If so, drain whole waterfraction*/
     32        /*drain only up to w0*/
    3333        if(dt==0.){
    34                 if(Dret>waterfraction)
    35                         return waterfraction;
     34                if((waterfraction>w0) && (waterfraction-Dret*yts<w0))
     35                        return waterfraction-w0;
     36                else
     37                        return Dret*yts;
     38        }
     39        else{
     40                if((waterfraction>w0) && (waterfraction-dt*Dret<w0))
     41                        return (waterfraction-w0)/dt;
    3642                else
    3743                        return Dret;
    3844        }
    39         else{
    40                 if(dt*Dret>waterfraction)
    41                         return waterfraction;
    42                 else
    43                         return dt*Dret;
    44         }
    4545}

  • issm/trunk-jpl/src/c/shared/Elements/LliboutryDuval.cpp

    r16322 r16816  
    3636
    3737        /*check feasibility*/
    38   _assert_(pressure>0);
     38  _assert_(pressure>=0);
    3939  _assert_(n>0);
    40   _assert_(betaCC>0);
    41   _assert_(referencetemperature>0);
     40  _assert_(betaCC>=0);
     41  _assert_(referencetemperature>=0);
    4242  _assert_(heatcapacity>0);
    4343  _assert_(latentheat>0);
    … …  
    5858    Tstar=Tpmp;
    5959    waterfraction=(enthalpy - H_sp)/latentheat;
     60                if (waterfraction > 0.01)
     61                        waterfraction = 0.01;
    6062  }
    6163
    … …  
    7577  return B;
    7678}
    77 
    78 // /*Get stiffness from temperature, waterfraction and depth*/
    79 // IssmDouble LliboutryDuval(IssmDouble temperature, IssmDouble waterfraction, IssmDouble depth,IssmDouble n){
    80 //      /*Use parameterization for the rheology: Greve and Blatter 2009
    81 //       * get enthalpy from temperature and water fraction,
    82 //       * and use LliboutryDuval(IssmDouble enthalpy, IssmDouble pressure,IssmDouble n) */
    83 //
    84 //      /*TODO: update params from model*/
    85 //   IssmDouble rho_ice=910; // kg/m^3
    86 //   IssmDouble g=9.81; //kg*m/s^2
    87 //   IssmDouble heatcapacity=2009; // J/kg/K
    88 //   IssmDouble referencetemperature=253.15;
    89 //   IssmDouble betaCC=7.9*pow(10.,-8.);
    90 //   IssmDouble latentheat=3.34*pow(10,5.); // from Aschwanden 2012
    91 //
    92 //   IssmDouble Tstar, enthalpy, pressure, B;
    93 //   _assert_(temperature>0);
    94 //   _assert_(waterfraction>0);
    95 //   _assert_(depth>0);
    96 //
    97 //   /*get pressure*/
    98 //   pressure= rho_ice*g*depth;
    99 //   Tstar=temperature-betaCC*pressure; // TODO: check whether plus or minus
    100 //   /*get enthalpy*/
    101 //   if (Tstar < 273.15){
    102 //     enthalpy=heatcapacity*(Tstar - referencetemperature);
    103 //   }
    104 //   else{
    105 //     enthalpy=heatcapacity*(273.15 - referencetemperature) + waterfraction*latentheat;
    106 //   }
    107 //
    108 //   B=LliboutryDuval(enthalpy, pressure, n, betaCC, referencetemperature, heatcapacity, latentheat);
    109 //
    110 //   return B;
    111 // }
Note: See TracChangeset for help on using the changeset viewer.