Changeset 27829


Ignore:
Timestamp:
07/16/23 18:22:25 (20 months ago)
Author:
inwoo
Message:

CHG: add calculation of TotalSmbMelt and TotalSmbRefreeze for SEMIC.

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

Legend:

Unmodified
Added
Removed

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

    r27789 r27829  
    43184318                _printf0_("smb core: dt       : " << dt <<"\n");
    43194319        }
    4320         /* loop over vertices and days */ //FIXME account for leap years (365 -> 366)
     4320        /* loop over vertices and days */
    43214321        Gauss* gauss=this->NewGauss();
    43224322        /* Retrieve inputs: */
    … …  
    43474347                enum_qmr        =SmbSemicQmrInitEnum;
    43484348        }
    4349         if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n");
     4349        //if(isverbose && this->Sid()==0)_printf0_("smb core: assign temp.\n");
    43504350        Input* tsurf_input       = this->GetInput(enum_temp); _assert_(tsurf_in);
    4351         if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n");
     4351        //if(isverbose && this->Sid()==0)_printf0_("smb core: assign mask.\n");
    43524352        Input* mask_input        = this->GetInput(SmbMaskEnum); _assert_(mask_input);
    4353         if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n");
     4353        //if(isverbose && this->Sid()==0)_printf0_("smb core: assign Tamp.\n");
    43544354        Input* Tamp_input        = this->GetInput(SmbTampEnum); _assert_(Tamp_input);
    4355         if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n");
     4355        //if(isverbose && this->Sid()==0)_printf0_("smb core: assign albedo.\n");
    43564356        Input* albedo_input      = this->GetInput(enum_albedo); _assert_(albedo_input);
    43574357        Input* albedo_snow_input = this->GetInput(enum_albedo_snow); _assert_(albedo_snow_input);
    43584358        Input* hice_input        = this->GetInput(enum_hice); _assert_(hice_input);
    43594359        Input* hsnow_input       = this->GetInput(enum_hsnow); _assert_(hsnow_input);
    4360         if(isverbose && this->Sid()==0)_printf0_("smb core: assign qmr.\n");
    43614360        Input* qmr_input         = this->GetInput(enum_qmr); _assert_(qmr_input);
    43624361
    … …  
    44384437                 */
    44394438                smb_out[iv]  = smb_out[iv]*yts;  // w.e. m/sec -> m/yr
    4440                 smbi_out[iv] = smbi_out[iv]*rho_water/rho_ice;
     4439                smbi_out[iv] = smbi_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
    44414440                smbs_out[iv] = smbs_out[iv]*yts; // w.e. m/sec -> m/yr
    44424441                saccu_out[iv] = saccu_out[iv]*yts; // w.e. m/sec -> m/yr
    4443         }
    4444         /*
    4445          * unit conversion is not required for smelt and refr_out.
    4446          */
     4442                smelt_out[iv] = smelt_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
     4443                refr_out[iv]  = refr_out[iv]*rho_water/rho_ice; // w.e. m/sec -> ice m/yr
     4444        }
    44474445
    44484446        if(isverbose && this->Sid()==0){
    … …  
    54065404}
    54075405/*}}}*/
     5406IssmDouble Element::TotalSmbMelt(IssmDouble* mask, bool scaled){/*{{{*/
     5407
     5408        /*Retrieve values of the mask defining the element: */
     5409        for(int i=0;i<this->GetNumberOfVertices();i++){
     5410                if(mask[this->vertices[i]->Sid()]<=0.){
     5411                        return 0.;
     5412                }
     5413        }
     5414
     5415        /*Return: */
     5416        return this->TotalSmbMelt(scaled);
     5417}
     5418/*}}}*/
     5419IssmDouble Element::TotalSmbRefreeze(IssmDouble* mask, bool scaled){/*{{{*/
     5420
     5421        /*Retrieve values of the mask defining the element: */
     5422        for(int i=0;i<this->GetNumberOfVertices();i++){
     5423                if(mask[this->vertices[i]->Sid()]<=0.){
     5424                        return 0.;
     5425                }
     5426        }
     5427
     5428        /*Return: */
     5429        return this->TotalSmbRefreeze(scaled);
     5430}
     5431/*}}}*/
    54085432void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
    54095433

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

    r27789 r27829  
    202202                IssmDouble         TotalGroundedBmb(IssmDouble* mask, bool scaled);
    203203                IssmDouble         TotalSmb(IssmDouble* mask, bool scaled);
     204                IssmDouble         TotalSmbMelt(IssmDouble* mask, bool scaled);
     205                IssmDouble         TotalSmbRefreeze(IssmDouble* mask, bool scaled);
    204206                void               TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int cs_enum);
    205207                void               TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int cs_enum);
    … …  
    387389                virtual IssmDouble TotalGroundedBmb(bool scaled)=0;
    388390                virtual IssmDouble TotalSmb(bool scaled)=0;
     391                virtual IssmDouble TotalSmbMelt(bool scaled)=0;
     392                virtual IssmDouble TotalSmbRefreeze(bool scaled)=0;
    389393                virtual void       Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finite_element)=0;
    390394                virtual void       UpdateConstraintsExtrudeFromBase(void)=0;

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

    r27724 r27829  
    43584358        /*Return: */
    43594359        return Total_Smb;
     4360}
     4361/*}}}*/
     4362IssmDouble Penta::TotalSmbMelt(bool scaled){/*{{{*/
     4363
     4364        /*The smbmelt[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
     4365        IssmDouble base,smbmelt,rho_ice,scalefactor;
     4366        IssmDouble Total_SmbMelt=0;
     4367        IssmDouble lsf[NUMVERTICES];
     4368        IssmDouble xyz_list[NUMVERTICES][3];
     4369
     4370        /*Get material parameters :*/
     4371        rho_ice=FindParam(MaterialsRhoIceEnum);
     4372
     4373        if(!IsIceInElement() || !IsOnSurface()) return 0.;
     4374
     4375        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     4376
     4377        /*First calculate the area of the base (cross section triangle)
     4378         * http://en.wikipedia.org/wiki/Triangle
     4379         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     4380        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     4381
     4382        /*Now get the average SMB over the element*/
     4383        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     4384   if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
     4385                /*Partially ice-covered element*/
     4386                bool mainlyice;
     4387      int point;
     4388      IssmDouble* smbmelt_vertices = xNew<IssmDouble>(NUMVERTICES);
     4389                IssmDouble  weights[NUMVERTICES2D];
     4390                IssmDouble  lsf2d[NUMVERTICES2D];
     4391      IssmDouble f1,f2,phi;
     4392      Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum);
     4393                for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
     4394                GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
     4395                smbmelt = 0.0;
     4396                for(int i=0;i<NUMVERTICES2D;i++) smbmelt += weights[i]*smbmelt_vertices[i];
     4397
     4398                if(scaled==true){
     4399         IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
     4400         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     4401         /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
     4402         scalefactor = 0.0;
     4403         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     4404         xDelete<IssmDouble>(scalefactor_vertices);
     4405                }
     4406                else scalefactor = 1.0;
     4407
     4408                /*Cleanup*/
     4409      xDelete<IssmDouble>(smbmelt_vertices);
     4410        }
     4411
     4412        else{
     4413                /*Fully ice-covered element*/
     4414                Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input);
     4415                smbmelt_input->GetInputAverage(&smbmelt);
     4416
     4417                if(scaled==true){
     4418                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     4419                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     4420                }
     4421                else scalefactor=1.0;
     4422        }
     4423
     4424        Total_SmbMelt=rho_ice*base*smbmelt*scalefactor;// smbmelt on element in kg s-1
     4425
     4426        /*Return: */
     4427        return Total_SmbMelt;
     4428}
     4429/*}}}*/
     4430IssmDouble Penta::TotalSmbRefreeze(bool scaled){/*{{{*/
     4431
     4432        /*The smbrefreeze[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
     4433        IssmDouble base,smbrefreeze,rho_ice,scalefactor;
     4434        IssmDouble Total_SmbRefreeze=0;
     4435        IssmDouble lsf[NUMVERTICES];
     4436        IssmDouble xyz_list[NUMVERTICES][3];
     4437
     4438        /*Get material parameters :*/
     4439        rho_ice=FindParam(MaterialsRhoIceEnum);
     4440
     4441        if(!IsIceInElement() || !IsOnSurface()) return 0.;
     4442
     4443        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     4444
     4445        /*First calculate the area of the base (cross section triangle)
     4446         * http://en.wikipedia.org/wiki/Triangle
     4447         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     4448        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
     4449
     4450        /*Now get the average SMB over the element*/
     4451        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     4452   if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
     4453                /*Partially ice-covered element*/
     4454                bool mainlyice;
     4455      int point;
     4456      IssmDouble* smbrefreeze_vertices = xNew<IssmDouble>(NUMVERTICES);
     4457                IssmDouble  weights[NUMVERTICES2D];
     4458                IssmDouble  lsf2d[NUMVERTICES2D];
     4459      IssmDouble f1,f2,phi;
     4460      Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum);
     4461                for(int i=0;i<NUMVERTICES2D;i++) lsf2d[i] = lsf[i];
     4462                GetFractionGeometry2D(weights,&phi,&point,&f1,&f2,&mainlyice,lsf2d);
     4463                smbrefreeze = 0.0;
     4464                for(int i=0;i<NUMVERTICES2D;i++) smbrefreeze += weights[i]*smbrefreeze_vertices[i];
     4465
     4466                if(scaled==true){
     4467         IssmDouble* scalefactor_vertices   = xNew<IssmDouble>(NUMVERTICES);
     4468         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     4469         /*Compute loop only over lower vertices: i<NUMVERTICES2D*/
     4470         scalefactor = 0.0;
     4471         for(int i=0;i<NUMVERTICES2D;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     4472         xDelete<IssmDouble>(scalefactor_vertices);
     4473                }
     4474                else scalefactor = 1.0;
     4475
     4476                /*Cleanup*/
     4477      xDelete<IssmDouble>(smbrefreeze_vertices);
     4478        }
     4479
     4480        else{
     4481                /*Fully ice-covered element*/
     4482                Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input);
     4483                smbrefreeze_input->GetInputAverage(&smbrefreeze);
     4484
     4485                if(scaled==true){
     4486                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     4487                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     4488                }
     4489                else scalefactor=1.0;
     4490        }
     4491
     4492        Total_SmbRefreeze=rho_ice*base*smbrefreeze*scalefactor;// smbrefreeze on element in kg s-1
     4493
     4494        /*Return: */
     4495        return Total_SmbRefreeze;
    43604496}
    43614497/*}}}*/

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

    r27308 r27829  
    200200                IssmDouble     TotalGroundedBmb(bool scaled);
    201201                IssmDouble     TotalSmb(bool scaled);
     202                IssmDouble     TotalSmbMelt(bool scaled);
     203                IssmDouble     TotalSmbRefreeze(bool scaled);
    202204                void           Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    203205                void           UpdateConstraintsExtrudeFromBase(void);

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

    r27308 r27829  
    158158                IssmDouble  TotalGroundedBmb(bool scaled){_error_("not implemented yet");};
    159159                IssmDouble  TotalSmb(bool scaled){_error_("not implemented yet");};
     160                IssmDouble  TotalSmbMelt(bool scaled){_error_("not implemented yet");};
     161                IssmDouble  TotalSmbRefreeze(bool scaled){_error_("not implemented yet");};
    160162                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement){_error_("not implemented yet");};
    161163                void        UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};

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

    r27308 r27829  
    167167                IssmDouble  TotalGroundedBmb(bool scaled){_error_("not implemented yet");};
    168168                IssmDouble  TotalSmb(bool scaled){_error_("not implemented yet");};
     169                IssmDouble  TotalSmbMelt(bool scaled){_error_("not implemented yet");};
     170                IssmDouble  TotalSmbRefreeze(bool scaled){_error_("not implemented yet");};
    169171                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    170172                void        UpdateConstraintsExtrudeFromBase(){_error_("not implemented");};

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

    r27828 r27829  
    58665866        /*Return: */
    58675867        return Total_Smb;
     5868}
     5869/*}}}*/
     5870IssmDouble Tria::TotalSmbMelt(bool scaled){/*{{{*/
     5871
     5872        /*The smbmelt[kg yr-1] of one element is area[m2] * smbmelt [kg m^-2 yr^-1]*/
     5873        IssmDouble base,smbmelt,rho_ice,scalefactor;
     5874        IssmDouble Total_Melt=0;
     5875        IssmDouble lsf[NUMVERTICES];
     5876        IssmDouble xyz_list[NUMVERTICES][3];
     5877
     5878        /*Get material parameters :*/
     5879        rho_ice=FindParam(MaterialsRhoIceEnum);
     5880
     5881   if(!IsIceInElement())return 0;
     5882
     5883        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5884
     5885        /*First calculate the area of the base (cross section triangle)
     5886         * http://en.wikipedia.org/wiki/Triangle
     5887         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     5888        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
     5889       
     5890        /*Now get the average SMB over the element*/
     5891        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     5892        if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
     5893                /*Partially ice-covered element*/
     5894                bool mainlyice;
     5895      int point;
     5896      IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
     5897      IssmDouble* smbmelt_vertices  = xNew<IssmDouble>(NUMVERTICES);
     5898      IssmDouble f1,f2,phi;
     5899               
     5900                Element::GetInputListOnVertices(&smbmelt_vertices[0],SmbMeltEnum);
     5901                GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
     5902                smbmelt = 0.0;
     5903                for(int i=0;i<NUMVERTICES;i++) smbmelt += weights[i]*smbmelt_vertices[i];
     5904       
     5905                if(scaled==true){
     5906         IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
     5907         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     5908         scalefactor = 0.0;
     5909         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     5910         xDelete<IssmDouble>(scalefactor_vertices);
     5911      }
     5912                else scalefactor = 1.0;
     5913
     5914                /*Cleanup*/
     5915      xDelete<IssmDouble>(weights);
     5916      xDelete<IssmDouble>(smbmelt_vertices);
     5917        }
     5918        else{
     5919                /*Fully ice-covered element*/
     5920                Input* smbmelt_input = this->GetInput(SmbMeltEnum); _assert_(smbmelt_input);
     5921                smbmelt_input->GetInputAverage(&smbmelt);   // average smbmelt on element in m ice s-1
     5922
     5923                if(scaled==true){
     5924                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     5925                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5926                }
     5927                else scalefactor=1.0;
     5928        }
     5929       
     5930   Total_Melt=rho_ice*base*smbmelt*scalefactor; // smbmelt on element in kg s-1
     5931
     5932        /*Return: */
     5933        return Total_Melt;
     5934}
     5935/*}}}*/
     5936IssmDouble Tria::TotalSmbRefreeze(bool scaled){/*{{{*/
     5937
     5938        /*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
     5939        IssmDouble base,smbrefreeze,rho_ice,scalefactor;
     5940        IssmDouble Total_Refreeze=0;
     5941        IssmDouble lsf[NUMVERTICES];
     5942        IssmDouble xyz_list[NUMVERTICES][3];
     5943
     5944        /*Get material parameters :*/
     5945        rho_ice=FindParam(MaterialsRhoIceEnum);
     5946
     5947   if(!IsIceInElement())return 0;
     5948
     5949        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
     5950
     5951        /*First calculate the area of the base (cross section triangle)
     5952         * http://en.wikipedia.org/wiki/Triangle
     5953         * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
     5954        base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1])); // area of element in m2
     5955       
     5956        /*Now get the average SMB over the element*/
     5957        Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
     5958        if(lsf[0]*lsf[1]<=0 || lsf[0]*lsf[2]<=0 || lsf[1]*lsf[2]<=0){
     5959                /*Partially ice-covered element*/
     5960                bool mainlyice;
     5961      int point;
     5962      IssmDouble* weights       = xNew<IssmDouble>(NUMVERTICES);
     5963      IssmDouble* smbrefreeze_vertices  = xNew<IssmDouble>(NUMVERTICES);
     5964      IssmDouble f1,f2,phi;
     5965               
     5966                Element::GetInputListOnVertices(&smbrefreeze_vertices[0],SmbRefreezeEnum);
     5967                GetFractionGeometry(weights,&phi,&point,&f1,&f2,&mainlyice,lsf);
     5968                smbrefreeze = 0.0;
     5969                for(int i=0;i<NUMVERTICES;i++) smbrefreeze += weights[i]*smbrefreeze_vertices[i];
     5970       
     5971                if(scaled==true){
     5972         IssmDouble* scalefactor_vertices = xNew<IssmDouble>(NUMVERTICES);
     5973         Element::GetInputListOnVertices(&scalefactor_vertices[0],MeshScaleFactorEnum);
     5974         scalefactor = 0.0;
     5975         for(int i=0;i<NUMVERTICES;i++) scalefactor += weights[i]/phi*scalefactor_vertices[i];
     5976         xDelete<IssmDouble>(scalefactor_vertices);
     5977      }
     5978                else scalefactor = 1.0;
     5979
     5980                /*Cleanup*/
     5981      xDelete<IssmDouble>(weights);
     5982      xDelete<IssmDouble>(smbrefreeze_vertices);
     5983        }
     5984        else{
     5985                /*Fully ice-covered element*/
     5986                Input* smbrefreeze_input = this->GetInput(SmbRefreezeEnum); _assert_(smbrefreeze_input);
     5987                smbrefreeze_input->GetInputAverage(&smbrefreeze);   // average smbrefreeze on element in m ice s-1
     5988
     5989                if(scaled==true){
     5990                        Input* scalefactor_input = this->GetInput(MeshScaleFactorEnum); _assert_(scalefactor_input);
     5991                        scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
     5992                }
     5993                else scalefactor=1.0;
     5994        }
     5995       
     5996   Total_Refreeze=rho_ice*base*smbrefreeze*scalefactor; // smbrefreeze on element in kg s-1
     5997
     5998        /*Return: */
     5999        return Total_Refreeze;
    58686000}
    58696001/*}}}*/

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

    r27789 r27829  
    158158                IssmDouble  TotalGroundedBmb(bool scaled);
    159159                IssmDouble  TotalSmb(bool scaled);
     160                IssmDouble  TotalSmbMelt(bool scaled);
     161                IssmDouble  TotalSmbRefreeze(bool scaled);
    160162                void        Update(Inputs* inputs,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
    161163                int         UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r27728 r27829  
    24722472                                        case TotalGroundedBmbScaledEnum:         this->TotalGroundedBmbx(&double_result,true);          break;
    24732473                                        case TotalSmbEnum:                       this->TotalSmbx(&double_result,false);                 break;
     2474                                        case TotalSmbMeltEnum:                   this->TotalSmbMeltx(&double_result,false);             break;
     2475                                        case TotalSmbRefreezeEnum:               this->TotalSmbRefreezex(&double_result,false);         break;
    24742476                                        case TotalSmbScaledEnum:                 this->TotalSmbx(&double_result,true);                  break;
    24752477
    … …  
    27372739                case TotalGroundedBmbScaledEnum:                          this->TotalGroundedBmbx(responses, true); break;
    27382740                case TotalSmbEnum:                                              this->TotalSmbx(responses, false); break;
     2741                case TotalSmbMeltEnum:                                       this->TotalSmbMeltx(responses, false); break;
     2742                case TotalSmbRefreezeEnum:                                        this->TotalSmbRefreezex(responses, false); break;
    27392743                case TotalSmbScaledEnum:                                          this->TotalSmbx(responses, true); break;
    27402744                case MaterialsRheologyBbarEnum:          this->ElementResponsex(responses,MaterialsRheologyBbarEnum); break;
    … …  
    31273131        /*Assign output pointers: */
    31283132        *pSmb=total_smb;
     3133
     3134}/*}}}*/
     3135void FemModel::TotalSmbMeltx(IssmDouble* pSmbMelt, bool scaled){/*{{{*/
     3136
     3137        IssmDouble local_smbmelt = 0;
     3138        IssmDouble total_smbmelt;
     3139
     3140        for(Object* & object : this->elements->objects){
     3141                Element* element = xDynamicCast<Element*>(object);
     3142                local_smbmelt+=element->TotalSmbMelt(scaled);
     3143        }
     3144        ISSM_MPI_Reduce(&local_smbmelt,&total_smbmelt,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     3145        ISSM_MPI_Bcast(&total_smbmelt,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     3146
     3147        /*Assign output pointers: */
     3148        *pSmbMelt=total_smbmelt;
     3149
     3150}/*}}}*/
     3151void FemModel::TotalSmbRefreezex(IssmDouble* pSmbRefreeze, bool scaled){/*{{{*/
     3152
     3153        IssmDouble local_smbrefreeze = 0;
     3154        IssmDouble total_smbrefreeze;
     3155
     3156        for(Object* & object : this->elements->objects){
     3157                Element* element = xDynamicCast<Element*>(object);
     3158                local_smbrefreeze+=element->TotalSmbRefreeze(scaled);
     3159        }
     3160        ISSM_MPI_Reduce(&local_smbrefreeze,&total_smbrefreeze,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
     3161        ISSM_MPI_Bcast(&total_smbrefreeze,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     3162
     3163        /*Assign output pointers: */
     3164        *pSmbRefreeze=total_smbrefreeze;
    31293165
    31303166}/*}}}*/

  • issm/trunk-jpl/src/c/classes/FemModel.h

    r26890 r27829  
    143143                void TotalGroundedBmbx(IssmDouble* pGbmb, bool scaled);
    144144                void TotalSmbx(IssmDouble* pSmb, bool scaled);
     145                void TotalSmbMeltx(IssmDouble* pSmbMelt, bool scaled);
     146                void TotalSmbRefreezex(IssmDouble* pSmbRefreeze, bool scaled);
    145147                #ifdef  _HAVE_DAKOTA_
    146148                void DakotaResponsesx(double* d_responses,char** responses_descriptors,int numresponsedescriptors,int d_numresponses);

  • issm/trunk-jpl/src/c/classes/Regionaloutput.cpp

    r25539 r27829  
    143143                                val_t+=element->TotalSmb(this->mask,true);
    144144                                break;
     145                        case TotalSmbMeltEnum:
     146                                val_t+=element->TotalSmbMelt(this->mask,true);
     147                                break;
     148                        case TotalSmbRefreezeEnum:
     149                                val_t+=element->TotalSmbRefreeze(this->mask,true);
     150                                break;
    145151                        default:
    146152                                _error_("Regional output type " << this->outputname << " not supported yet!");

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

    r27780 r27829  
    13921392syn keyword cConstant CalvingTestEnum
    13931393syn keyword cConstant CalvingParameterizationEnum
     1394syn keyword cConstant CalvingCalvingMIPEnum
    13941395syn keyword cConstant CalvingVonmisesEnum
    13951396syn keyword cConstant CalvingVonmisesADEnum
    … …  
    17251726syn keyword cConstant TotalSmbEnum
    17261727syn keyword cConstant TotalSmbScaledEnum
     1728syn keyword cConstant TotalSmbRefreezeEnum
     1729syn keyword cConstant TotalSmbMeltEnum
    17271730syn keyword cConstant TransientArrayParamEnum
    17281731syn keyword cConstant TransientInputEnum
    … …  
    17771780syn keyword cType Cfsurfacesquaretransient
    17781781syn keyword cType Channel
    1779 syn keyword cType classes
    17801782syn keyword cType Constraint
    17811783syn keyword cType Constraints
    … …  
    17851787syn keyword cType ControlParam
    17861788syn keyword cType Covertree
     1789syn keyword cType DataSetParam
    17871790syn keyword cType DatasetInput
    1788 syn keyword cType DataSetParam
    17891791syn keyword cType Definition
    17901792syn keyword cType DependentObject
    … …  
    17991801syn keyword cType ElementInput
    18001802syn keyword cType ElementMatrix
     1803syn keyword cType ElementVector
    18011804syn keyword cType Elements
    1802 syn keyword cType ElementVector
    18031805syn keyword cType ExponentialVariogram
    18041806syn keyword cType ExternalResult
    … …  
    18071809syn keyword cType Friction
    18081810syn keyword cType Gauss
    1809 syn keyword cType GaussianVariogram
    1810 syn keyword cType gaussobjects
    18111811syn keyword cType GaussPenta
    18121812syn keyword cType GaussSeg
    18131813syn keyword cType GaussTetra
    18141814syn keyword cType GaussTria
     1815syn keyword cType GaussianVariogram
    18151816syn keyword cType GenericExternalResult
    18161817syn keyword cType GenericOption
    … …  
    18291830syn keyword cType IssmDirectApplicInterface
    18301831syn keyword cType IssmParallelDirectApplicInterface
    1831 syn keyword cType krigingobjects
    18321832syn keyword cType Load
    18331833syn keyword cType Loads
    … …  
    18401840syn keyword cType Matice
    18411841syn keyword cType Matlitho
    1842 syn keyword cType matrixobjects
    18431842syn keyword cType MatrixParam
    18441843syn keyword cType Misfit
    … …  
    18531852syn keyword cType Observations
    18541853syn keyword cType Option
     1854syn keyword cType OptionUtilities
    18551855syn keyword cType Options
    1856 syn keyword cType OptionUtilities
    18571856syn keyword cType Param
    18581857syn keyword cType Parameters
    … …  
    18681867syn keyword cType Regionaloutput
    18691868syn keyword cType Results
     1869syn keyword cType RiftStruct
    18701870syn keyword cType Riftfront
    1871 syn keyword cType RiftStruct
    18721871syn keyword cType SealevelGeometry
    18731872syn keyword cType Seg
    18741873syn keyword cType SegInput
     1874syn keyword cType SegRef
    18751875syn keyword cType Segment
    1876 syn keyword cType SegRef
    18771876syn keyword cType SpcDynamic
    18781877syn keyword cType SpcStatic
    … …  
    18931892syn keyword cType Vertex
    18941893syn keyword cType Vertices
     1894syn keyword cType classes
     1895syn keyword cType gaussobjects
     1896syn keyword cType krigingobjects
     1897syn keyword cType matrixobjects
    18951898syn keyword cType AdjointBalancethickness2Analysis
    18961899syn keyword cType AdjointBalancethicknessAnalysis

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

    r27789 r27829  
    17251725        TotalSmbEnum,
    17261726        TotalSmbScaledEnum,
     1727        TotalSmbRefreezeEnum,
     1728        TotalSmbMeltEnum,
    17271729        TransientArrayParamEnum,
    17281730        TransientInputEnum,

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

    r27780 r27829  
    13941394                case CalvingTestEnum : return "CalvingTest";
    13951395                case CalvingParameterizationEnum : return "CalvingParameterization";
     1396                case CalvingCalvingMIPEnum : return "CalvingCalvingMIP";
    13961397                case CalvingVonmisesEnum : return "CalvingVonmises";
    13971398                case CalvingVonmisesADEnum : return "CalvingVonmisesAD";
    … …  
    17271728                case TotalSmbEnum : return "TotalSmb";
    17281729                case TotalSmbScaledEnum : return "TotalSmbScaled";
     1730                case TotalSmbRefreezeEnum : return "TotalSmbRefreeze";
     1731                case TotalSmbMeltEnum : return "TotalSmbMelt";
    17291732                case TransientArrayParamEnum : return "TransientArrayParam";
    17301733                case TransientInputEnum : return "TransientInput";

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

    r27780 r27829  
    13851385syn keyword juliaConstC CalvingTestEnum
    13861386syn keyword juliaConstC CalvingParameterizationEnum
     1387syn keyword juliaConstC CalvingCalvingMIPEnum
    13871388syn keyword juliaConstC CalvingVonmisesEnum
    13881389syn keyword juliaConstC CalvingVonmisesADEnum
    … …  
    17181719syn keyword juliaConstC TotalSmbEnum
    17191720syn keyword juliaConstC TotalSmbScaledEnum
     1721syn keyword juliaConstC TotalSmbRefreezeEnum
     1722syn keyword juliaConstC TotalSmbMeltEnum
    17201723syn keyword juliaConstC TransientArrayParamEnum
    17211724syn keyword juliaConstC TransientInputEnum

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

    r27780 r27829  
    14271427              else if (strcmp(name,"CalvingTest")==0) return CalvingTestEnum;
    14281428              else if (strcmp(name,"CalvingParameterization")==0) return CalvingParameterizationEnum;
     1429              else if (strcmp(name,"CalvingCalvingMIP")==0) return CalvingCalvingMIPEnum;
    14291430              else if (strcmp(name,"CalvingVonmises")==0) return CalvingVonmisesEnum;
    14301431              else if (strcmp(name,"CalvingVonmisesAD")==0) return CalvingVonmisesADEnum;
    … …  
    14891490              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    14901491              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    1491               else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    14921492         else stage=13;
    14931493   }
    14941494   if(stage==13){
    1495               if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
     1495              if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
     1496              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
    14961497              else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
    14971498              else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
    … …  
    16121613              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    16131614              else if (strcmp(name,"Matestar")==0) return MatestarEnum;
    1614               else if (strcmp(name,"Matice")==0) return MaticeEnum;
    16151615         else stage=14;
    16161616   }
    16171617   if(stage==14){
    1618               if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
     1618              if (strcmp(name,"Matice")==0) return MaticeEnum;
     1619              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
    16191620              else if (strcmp(name,"Mathydro")==0) return MathydroEnum;
    16201621              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    … …  
    17351736              else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
    17361737              else if (strcmp(name,"Sset")==0) return SsetEnum;
    1737               else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
    17381738         else stage=15;
    17391739   }
    17401740   if(stage==15){
    1741               if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
     1741              if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum;
     1742              else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum;
    17421743              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    17431744              else if (strcmp(name,"StressbalanceAnalysis")==0) return StressbalanceAnalysisEnum;
    … …  
    17691770              else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum;
    17701771              else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum;
     1772              else if (strcmp(name,"TotalSmbRefreeze")==0) return TotalSmbRefreezeEnum;
     1773              else if (strcmp(name,"TotalSmbMelt")==0) return TotalSmbMeltEnum;
    17711774              else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum;
    17721775              else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;

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

    r27826 r27829  
    9191                                        'SmbMassBalanceSemic','SmbMelt','SmbRefreeze','SmbAccumulation',...
    9292                                        'SmbHIce','SmbHSnow','SmbAlbedo','SmbAlbedoSnow','TemperatureSEMIC',...
    93                                         'SmbSemicQmr'};
     93                                        'SmbSemicQmr','TotalSmb','TotalSmbMelt','TotalSmbRefreeze'};
    9494                        else
    9595                                list = {'default','SmbMassBalance'};

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

    r21063 r27829  
    133133        elseif strcmp(fieldname,'TotalSmb'),
    134134                field = field/10.^12*yts; %(GigaTon/year)
     135        elseif strcmp(fieldname,'TotalMelt'),
     136                field = field/10.^12*yts; %(GigaTon/year)
     137        elseif strcmp(fieldname,'TotalRefreeze'),
     138                field = field/10.^12*yts; %(GigaTon/year)
    135139        elseif strcmp(fieldname,'SmbMassBalance'),
    136140                field = field*yts;

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

    r26790 r27829  
    260260                field = field/10.^12*yts; %(GigaTon/year)
    261261        elseif strcmp(fieldname,'TotalSmbScaled'),
     262                field = field/10.^12*yts; %(GigaTon/year)
     263        elseif strcmp(fieldname,'TotalSmbMelt'),
     264                field = field/10.^12*yts; %(GigaTon/year)
     265        elseif strcmp(fieldname,'TotalSmbRefreeze'),
    262266                field = field/10.^12*yts; %(GigaTon/year)
    263267        elseif strcmp(fieldname,'GroundinglineMassFlux'),
Note: See TracChangeset for help on using the changeset viewer.