Changeset 25150


Ignore:
Timestamp:
06/24/20 22:26:13 (5 years ago)
Author:
Eric.Larour
Message:

CHG: materials hydro class was missing fresh water density for hydrological applications.
Also added hydrological fingerprints capabilities for sealevelchange core.

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

Legend:

Unmodified
Added
Removed

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

    r25118 r25150  
    56625662        this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea);
    56635663
     5664        /*Compute hydrological contribution to sea level rise: */
     5665        this->SealevelriseEustaticHydro(Sgi,peustatic,masks, oceanarea);
     5666
     5667
    56645668}
    56655669/*}}}*/
    … …  
    57775781        /*convolve:*/
    57785782        for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*I;
     5783
     5784        /*Assign output pointer:*/
     5785        _assert_(!xIsNan<IssmDouble>(eustatic));
     5786        *peustatic=eustatic;
     5787        return;
     5788}
     5789/*}}}*/
     5790void    Tria::SealevelriseEustaticHydro(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/
     5791
     5792        /*diverse:*/
     5793        int gsize;
     5794        IssmDouble area;
     5795        IssmDouble phi=1.0; //WARNING: do not touch this, default is entire elemnt contributes eustatic
     5796        IssmDouble W;  //change in water height thickness (Farrel and Clarke, Equ. 4)
     5797        bool notfullygrounded=false;
     5798        bool scaleoceanarea= false;
     5799
     5800        /*elastic green function:*/
     5801        IssmDouble* G=NULL;
     5802
     5803        /*ice properties: */
     5804        IssmDouble rho_water;
     5805        IssmDouble rho_freshwater;
     5806
     5807        /*constants:*/
     5808        IssmDouble constant=0;
     5809
     5810        /*Initialize eustatic component: do not skip this step :):*/
     5811        IssmDouble eustatic = 0.;
     5812
     5813        /*early return if we are on an ice cap:*/
     5814        if(masks->isiceonly[this->lid]){ *peustatic=0; return; }
     5815
     5816        /*early return if we are fully floating:*/
     5817        if (masks->isfullyfloating[this->lid]){ constant=0; *peustatic=0; return; }
     5818
     5819        /*If we are here, we are on earth, not on ice: */
     5820
     5821        /*recover material parameters: */
     5822        rho_water=FindParam(MaterialsRhoSeawaterEnum);
     5823        rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum);
     5824
     5825        /*recover ocean area scaling: */
     5826        this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum);
     5827
     5828        /*retrieve precomputed G:*/
     5829        this->inputs2->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize);
     5830
     5831        /*Get area of element: precomputed in the sealevelrise_core_geometry:*/
     5832        this->GetInput2Value(&area,AreaEnum);
     5833
     5834        /*Retrieve water height at vertices: */
     5835        Input2* deltathickness_input=this->GetInput2(SurfaceloadWaterHeightChangeEnum);
     5836        if (!deltathickness_input)_error_("SurfaceloadWaterHeightChangeEnum input needed to compute sea level rise!");
     5837        deltathickness_input->GetInputAverage(&W);
     5838
     5839        /*Compute eustatic component:*/
     5840        _assert_(oceanarea>0.);
     5841        if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
     5842        eustatic += rho_freshwater*area*phi*W/(oceanarea*rho_water);
     5843
     5844        /*convert from m to kg/m^2:*/
     5845        W=W*rho_freshwater*phi;
     5846
     5847        /*convolve:*/
     5848        for(int i=0;i<gsize;i++) Sgi[i]+=G[i]*W;
    57795849
    57805850        /*Assign output pointer:*/

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

    r25066 r25150  
    169169                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea);
    170170                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea);
     171                void    SealevelriseEustaticHydro(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea);
    171172                void    SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks);
    172173                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks);

  • issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp

    r25128 r25150  
    264264        Vector<IssmDouble> *steric_rate_g  = NULL;
    265265        Vector<IssmDouble> *dynamic_rate_g = NULL;
    266         Vector<IssmDouble> *hydro_rate_g  = NULL;
    267266        Vector<IssmDouble> *U_esa_rate= NULL;
    268267        Vector<IssmDouble> *N_esa_rate= NULL;
    … …  
    296295        GetStericRate(&steric_rate_g,femmodel);
    297296        GetDynamicRate(&dynamic_rate_g,femmodel);
    298         GetVectorFromInputsx(&hydro_rate_g,femmodel,SurfaceloadWaterHeightChangeEnum,VertexSIdEnum);
    299297        if(geodetic){
    300298                GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
    … …  
    304302        }
    305303
    306         /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate + hydro_rate* dt*/
     304        /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate)  * dt + steric_rate + dynamic_rate dt*/
    307305        if(geodetic){
    308306                SL->AXPY(N_gia_rate,dt);
    … …  
    311309        SL->AXPY(steric_rate_g,dt);
    312310        SL->AXPY(dynamic_rate_g,dt);
    313         SL->AXPY(hydro_rate_g,dt);
    314311
    315312        /*compute new bedrock position: */
    … …  
    328325        delete steric_rate_g;
    329326        delete dynamic_rate_g;
    330         delete hydro_rate_g;
    331327        if(geodetic){
    332328                delete U_esa_rate;

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

    r25066 r25150  
    365365                                                parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_water",MaterialsRhoSeawaterEnum));
    366366                                                parameters->AddObject(iomodel->CopyConstantObject("md.materials.earth_density",MaterialsEarthDensityEnum));
     367                                                parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_freshwater",MaterialsRhoFreshwaterEnum));
    367368                                                break;
    368369                                }

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

    r24756 r25150  
    5858                                        self.addprop('rho_ice');
    5959                                        self.addprop('rho_water');
     60                                        self.addprop('rho_freshwater');
    6061                                        self.addprop('earth_density');
    6162                                otherwise
    … …  
    138139                                        % average density of the Earth (kg/m^3)
    139140                                        self.earth_density=5512;
     141
     142                                        %fresh water density (kg/m^3)
     143                                        self.rho_freshwater=1000.;
    140144
    141145                                otherwise
    … …  
    184188                                        fielddisplay(self,'rho_water','ocean water density [kg/m^3]');
    185189                                        fielddisplay(self,'earth_density','mantle density [kg/m^3]');
     190                                        fielddisplay(self,'rho_freshwater','fresh water density [kg/m^3]');
    186191
    187192                                otherwise
    … …  
    233238                                        md = checkfield(md,'fieldname','materials.rho_water','>',0);
    234239                                        md = checkfield(md,'fieldname','materials.earth_density','>',0,'numel',1);
     240                                        md = checkfield(md,'fieldname','materials.rho_freshwater','>',0);
    235241
    236242                                otherwise
    … …  
    280286                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_water','format','Double');
    281287                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','earth_density','format','Double');
     288                                        WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_freshwater','format','Double');
    282289
    283290                                otherwise
    … …  
    334341                                        writejsdouble(fid,[modelname '.materials.rho_water'],self.rho_water);
    335342                                        writejsdouble(fid,[modelname '.materials.earth_density'],self.earth_density);
     343                                        writejsdouble(fid,[modelname '.materials.rho_freshwater'],self.rho_freshwater);
    336344
    337345                                otherwise
Note: See TracChangeset for help on using the changeset viewer.