Ignore:
Timestamp:
05/23/11 18:46:37 (14 years ago)
Author:
schlegel
Message:

change accumulation rate to surface mass balance; introduce surface ablation and surface accumulation

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r8392 r8399  
    16611661        int        i,j,ig;
    16621662        double     xyz_list[NUMVERTICES][3];
    1663         double     dhdt_g,melting_g,accumulation_g,Jdettria;
     1663        double     dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
    16641664        double     L[NUMVERTICES];
    16651665        GaussTria* gauss=NULL;
    … …  
    16701670        /*Retrieve all inputs and parameters*/
    16711671        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1672         Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); _assert_(accumulation_input);
    1673         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
     1672        Input* surface_mass_balance_input=inputs->GetInput(SurfaceMassBalanceEnum); _assert_(surface_mass_balance_input);
     1673        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
    16741674        Input* dhdt_input=inputs->GetInput(DhDtEnum);                     _assert_(dhdt_input);
    16751675       
    … …  
    16801680                gauss->GaussPoint(ig);
    16811681
    1682                 accumulation_input->GetParameterValue(&accumulation_g,gauss);
    1683                 melting_input->GetParameterValue(&melting_g,gauss);
     1682                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     1683                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    16841684                dhdt_input->GetParameterValue(&dhdt_g,gauss);
    16851685
    … …  
    16871687                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    16881688
    1689                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
     1689                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
    16901690        }
    16911691
    … …  
    17041704        int        i,j,ig;
    17051705        double     xyz_list[NUMVERTICES][3];
    1706         double     melting_g,accumulation_g,dhdt_g,Jdettria;
     1706        double     basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
    17071707        double     L[NUMVERTICES];
    17081708        GaussTria* gauss=NULL;
    … …  
    17131713        /*Retrieve all inputs and parameters*/
    17141714        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1715         Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); _assert_(accumulation_input);
    1716         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
     1715        Input* surface_mass_balance_input=inputs->GetInput(SurfaceMassBalanceEnum); _assert_(surface_mass_balance_input);
     1716        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
    17171717        Input* dhdt_input=inputs->GetInput(DhDtEnum);                     _assert_(dhdt_input);
    17181718
    … …  
    17231723                gauss->GaussPoint(ig);
    17241724
    1725                 accumulation_input->GetParameterValue(&accumulation_g,gauss);
    1726                 melting_input->GetParameterValue(&melting_g,gauss);
     1725                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     1726                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    17271727                dhdt_input->GetParameterValue(&dhdt_g,gauss);
    17281728
    … …  
    17301730                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    17311731
    1732                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g-dhdt_g)*L[i];
     1732                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
    17331733        }
    17341734
    … …  
    17471747        int        i,j,ig;
    17481748        double     xyz_list[NUMVERTICES][3];
    1749         double     Jdettria,accumulation_g,melting_g;
     1749        double     Jdettria,surface_mass_balance_g,basal_melting_g;
    17501750        double     L[NUMVERTICES];
    17511751        GaussTria* gauss=NULL;
    … …  
    17561756        /*Retrieve all inputs and parameters*/
    17571757        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    1758         Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); _assert_(accumulation_input);
    1759         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
     1758        Input* surface_mass_balance_input=inputs->GetInput(SurfaceMassBalanceEnum); _assert_(surface_mass_balance_input);
     1759        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
    17601760
    17611761        /* Start  looping on the number of gaussian points: */
    … …  
    17681768                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    17691769
    1770                 accumulation_input->GetParameterValue(&accumulation_g,gauss);
    1771                 melting_input->GetParameterValue(&melting_g,gauss);
    1772 
    1773                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(accumulation_g-melting_g)*L[i];
     1770                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     1771                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
     1772
     1773                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g)*L[i];
    17741774        }
    17751775
    … …  
    17901790        double     xyz_list[NUMVERTICES][3];
    17911791        double     Jdet;
    1792         double     vx,vy,vz,dbdx,dbdy,meltingvalue;
     1792        double     vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
    17931793        double     slope[2];
    17941794        double     L[NUMVERTICES];
    … …  
    18021802        inputs->GetParameterValue(&approximation,ApproximationEnum);
    18031803        Input* bed_input=inputs->GetInput(BedEnum);             _assert_(bed_input);
    1804         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum); _assert_(melting_input);
     1804        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum); _assert_(basal_melting_input);
    18051805        Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    18061806        Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    … …  
    18161816                gauss->GaussPoint(ig);
    18171817
    1818                 melting_input->GetParameterValue(&meltingvalue, gauss);
     1818                basal_melting_input->GetParameterValue(&basalmeltingvalue, gauss);
    18191819                bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    18201820                vx_input->GetParameterValue(&vx, gauss);
    … …  
    18311831                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    18321832
    1833                 for(i=0;i<numdof;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i];
     1833                for(i=0;i<numdof;i++) pe->values[i]+=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*L[i];
    18341834        }
    18351835
    … …  
    23892389        int        i,j,ig;
    23902390        double     Jdettria,dt;
    2391         double     melting_g;
     2391        double     basal_melting_g;
    23922392        double     old_watercolumn_g;
    23932393        double     xyz_list[NUMVERTICES][3];
    … …  
    24012401        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    24022402        this->parameters->FindParam(&dt,DtEnum);
    2403         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
     2403        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
    24042404        Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);_assert_(old_watercolumn_input);
    24052405
    2406         /*Initialize melting_correction_g to 0, do not forget!:*/
     2406        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    24072407        /* Start  looping on the number of gaussian points: */
    24082408        gauss=new GaussTria(2);
    … …  
    24142414                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    24152415
    2416                 melting_input->GetParameterValue(&melting_g,gauss);
     2416                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    24172417                old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss);
    24182418       
    2419                 if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*melting_g)*L[i];
    2420                 else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*melting_g*L[i];
     2419                if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*L[i];
     2420                else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*L[i];
    24212421        }
    24222422               
    … …  
    24352435        int        i,j,ig;
    24362436        double     Jdettria,dt;
    2437         double     accumulation_g,melting_g,melting_correction_g,thickness_g;
     2437        double     surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
    24382438        double     xyz_list[NUMVERTICES][3];
    24392439        double     L[NUMVERTICES];
    … …  
    24462446        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    24472447        this->parameters->FindParam(&dt,DtEnum);
    2448         Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); _assert_(accumulation_input);
    2449         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
    2450         Input* melting_correction_input=inputs->GetInput(BasalMeltingRateCorrectionEnum);
     2448        Input* surface_mass_balance_input=inputs->GetInput(SurfaceMassBalanceEnum); _assert_(surface_mass_balance_input);
     2449        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
     2450        Input* basal_melting_correction_input=inputs->GetInput(BasalMeltingRateCorrectionEnum);
    24512451        Input* thickness_input=inputs->GetInput(ThicknessEnum);           _assert_(thickness_input);
    24522452
    2453         /*Initialize melting_correction_g to 0, do not forget!:*/
     2453        /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    24542454        /* Start  looping on the number of gaussian points: */
    24552455        gauss=new GaussTria(2);
    … …  
    24612461                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    24622462
    2463                 accumulation_input->GetParameterValue(&accumulation_g,gauss);
    2464                 melting_input->GetParameterValue(&melting_g,gauss);
     2463                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     2464                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    24652465                thickness_input->GetParameterValue(&thickness_g,gauss);
    2466                 if(melting_correction_input) melting_correction_input->GetParameterValue(&melting_correction_g,gauss);
    2467 
    2468                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g-melting_correction_g))*L[i];
     2466                if(basal_melting_correction_input) basal_melting_correction_input->GetParameterValue(&basal_melting_correction_g,gauss);
     2467
     2468                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*L[i];
    24692469        }
    24702470
    … …  
    24832483        int        i,j,ig;
    24842484        double     Jdettria,dt;
    2485         double     accumulation_g,melting_g,thickness_g;
     2485        double     surface_mass_balance_g,basal_melting_g,thickness_g;
    24862486        double     xyz_list[NUMVERTICES][3];
    24872487        double     L[NUMVERTICES];
    … …  
    24942494        this->parameters->FindParam(&dt,DtEnum);
    24952495        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    2496         Input* accumulation_input=inputs->GetInput(AccumulationRateEnum); _assert_(accumulation_input);
    2497         Input* melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(melting_input);
     2496        Input* surface_mass_balance_input=inputs->GetInput(SurfaceMassBalanceEnum); _assert_(surface_mass_balance_input);
     2497        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
    24982498        Input* thickness_input=inputs->GetInput(ThicknessEnum);           _assert_(thickness_input);
    24992499
    … …  
    25072507                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
    25082508
    2509                 accumulation_input->GetParameterValue(&accumulation_g,gauss);
    2510                 melting_input->GetParameterValue(&melting_g,gauss);
     2509                surface_mass_balance_input->GetParameterValue(&surface_mass_balance_g,gauss);
     2510                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    25112511                thickness_input->GetParameterValue(&thickness_g,gauss);
    25122512
    2513                 for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
     2513                for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*L[i];
    25142514        }
    25152515
    … …  
    38853885                this->inputs->AddInput(new TriaVertexInput(WaterColumnOldEnum,nodeinputs));
    38863886        }
    3887 
    3888         if (iomodel->accumulation_rate) {
    3889                 for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
    3890                 this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs));
     3887        if (iomodel->surface_accumulation_rate) {
     3888                for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface_accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
     3889                this->inputs->AddInput(new TriaVertexInput(SurfaceAccumulationRateEnum,nodeinputs));
     3890        }
     3891        if (iomodel->surface_ablation_rate) {
     3892                for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface_ablation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
     3893                this->inputs->AddInput(new TriaVertexInput(SurfaceAblationRateEnum,nodeinputs));
     3894        }
     3895        if (iomodel->surface_mass_balance) {
     3896                for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface_mass_balance[tria_vertex_ids[i]-1]/iomodel->yts;
     3897                this->inputs->AddInput(new TriaVertexInput(SurfaceMassBalanceEnum,nodeinputs));
    38913898        }
    38923899        if (iomodel->geothermalflux) {
    … …  
    45304537                                name==BasalMeltingRateEnum ||
    45314538                                name==WaterColumnEnum ||
    4532                                 name==AccumulationRateEnum ||
     4539                                name==SurfaceMassBalanceEnum ||
    45334540                                name==SurfaceAreaEnum||
    45344541                                name==VxEnum ||
    … …  
    51465153    this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf));
    51475154
    5148         /*If this element just  became ungrounded, set its melting rate at 50 m/yr:*/
     5155        /*If this element just  became ungrounded, set its basal melting rate at 50 m/yr:*/
    51495156        if(swap){
    51505157                Input* basal_melting_rate_input     =inputs->GetInput(BasalMeltingRateEnum);     _assert_(basal_melting_rate_input);
    … …  
    58865893
    58875894        if(step==0){
    5888                 /*This is the first time we are trying to replace the accumulation tria vertex input by an accumulation triav vertex
    5889                  *forcing, which is essentially a tria vertex input that can vary with time. So firt, if there is already
     5895                /*This is the first time we are trying to replace the smb tria vertex input by a smb triav vertex
     5896                 *forcing, which is essentially a tria vertex input that can vary with time. So first, if there is already
    58905897                 an input with enum FieldEnum, squash it: */
    58915898                this->inputs->AddInput(new TriaVertexForcing(FieldEnum,nodeinputs,time,iomodel->forcing_numtimesteps,this->parameters));
Note: See TracChangeset for help on using the changeset viewer.