Changeset 21392


Ignore:
Timestamp:
11/17/16 21:09:42 (8 years ago)
Author:
Eric.Larour
Message:

CHG: rotational feedback.

Location:
issm/branches/trunk-larour-NatClimateChange2016/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Element.h

    r20999 r21392  
    301301                virtual IssmDouble    OceanAverage(IssmDouble* Sg)=0;
    302302                virtual IssmDouble    OceanArea(void)=0;
     303                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea)=0;
    303304                virtual void          SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
    304305                virtual void          SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Penta.h

    r20999 r21392  
    186186                IssmDouble    OceanArea(void){_error_("not implemented yet!");};
    187187                IssmDouble    OceanAverage(IssmDouble* Sg){_error_("not implemented yet!");};
     188                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    188189                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    189190                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Seg.h

    r20999 r21392  
    170170
    171171#ifdef _HAVE_SEALEVELRISE_
     172                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    172173                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    173174                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tetra.h

    r20999 r21392  
    176176#endif
    177177#ifdef _HAVE_SEALEVELRISE_
     178                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){_error_("not implemented yet!");};
    178179                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    179180                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.cpp

    r21206 r21392  
    36113611}
    36123612/*}}}*/
     3613void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/
     3614
     3615        /*early return if we are not on an ice cap OR ocean:*/
     3616        if(!(this->inputs->Max(MaskIceLevelsetEnum)<0) && !IsWaterInElement()) return;
     3617
     3618        /*Compute area of element:*/
     3619        IssmDouble area;
     3620        area=GetAreaSpherical();
     3621
     3622        /*Compute lat,long,radius of elemental centroid: */
     3623        bool spherical=true;
     3624        IssmDouble llr_list[NUMVERTICES][3];
     3625        IssmDouble late,longe,re;
     3626        /* Where is the centroid of this element?:{{{*/
     3627        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
     3628
     3629        IssmDouble minlong=400;
     3630        IssmDouble maxlong=-20;
     3631        for (int i=0;i<NUMVERTICES;i++){
     3632                llr_list[i][0]=(90-llr_list[i][0]);
     3633                if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
     3634                if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
     3635                if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
     3636        }
     3637        if(minlong==0 && maxlong>180){
     3638                if (llr_list[0][1]==0)llr_list[0][1]=360;
     3639                if (llr_list[1][1]==0)llr_list[1][1]=360;
     3640                if (llr_list[2][1]==0)llr_list[2][1]=360;
     3641        }
     3642
     3643        // correction at the north pole
     3644        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3645        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3646        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3647
     3648        //correction at the south pole
     3649        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
     3650        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
     3651        if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
     3652
     3653        late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
     3654        longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;
     3655
     3656        late=90-late;
     3657        if(longe>180)longe=(longe-180)-180;
     3658
     3659        late=late/180*PI;
     3660        longe=longe/180*PI;
     3661        /*}}}*/
     3662        re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;
     3663
     3664        if(IsWaterInElement()){
     3665                IssmDouble rho_water, S;
     3666               
     3667                /*recover material parameters: */
     3668                rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
     3669       
     3670                /*From Sg_old, recover water sea level rise:*/
     3671                S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
     3672               
     3673                /* Perturbation terms for moment of inertia (moi_list):
     3674                 * computed analytically (see Wu & Peltier, eqs 10 & 32)
     3675                 * also consistent with my GMD formulation!
     3676                 * ALL in geographic coordinates
     3677                 * */
     3678                dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
     3679                dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
     3680                dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
     3681        }
     3682        else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
     3683                IssmDouble rho_ice, I;
     3684               
     3685                /*recover material parameters: */
     3686                rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     3687       
     3688                /*Compute ice thickness change: */
     3689                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
     3690                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
     3691                deltathickness_input->GetInputAverage(&I);
     3692               
     3693                dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
     3694                dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
     3695                dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
     3696        }
     3697       
     3698        return;
     3699}/*}}}*/
    36133700void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    36143701

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/classes/Elements/Tria.h

    r20999 r21392  
    147147                IssmDouble OceanArea(void);
    148148                IssmDouble OceanAverage(IssmDouble* Sg);
     149                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea);
    149150                void    SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    150151                void    SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumDefinitions.h

    r21232 r21392  
    772772        SealevelriseRigidEnum,
    773773        SealevelriseElasticEnum,
     774        SealevelriseRotationEnum,
     775        SealevelriseTidalLoveHEnum,
     776        SealevelriseTidalLoveKEnum,
     777        SealevelriseFluidLoveEnum,
     778        SealevelriseEquatorialMoiEnum,
     779        SealevelrisePolarMoiEnum,
     780        SealevelriseAngularVelocityEnum,
    774781        SealevelriseGElasticEnum,
    775782        SealevelriseUElasticEnum,

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/EnumToStringx.cpp

    r21232 r21392  
    753753                case SealevelriseRigidEnum : return "SealevelriseRigid";
    754754                case SealevelriseElasticEnum : return "SealevelriseElastic";
     755                case SealevelriseRotationEnum : return "SealevelriseRotation";
     756                case SealevelriseTidalLoveHEnum : return "SealevelriseTidalLoveH";
     757                case SealevelriseTidalLoveKEnum : return "SealevelriseTidalLoveK";
     758                case SealevelriseFluidLoveEnum : return "SealevelriseFluidLove";
     759                case SealevelriseEquatorialMoiEnum : return "SealevelriseEquatorialMoi";
     760                case SealevelrisePolarMoiEnum : return "SealevelrisePolarMoi";
     761                case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";
    755762                case SealevelriseGElasticEnum : return "SealevelriseGElastic";
    756763                case SealevelriseUElasticEnum : return "SealevelriseUElastic";

  • issm/branches/trunk-larour-NatClimateChange2016/src/c/shared/Enum/StringToEnumx.cpp

    r21232 r21392  
    771771              else if (strcmp(name,"SealevelriseRigid")==0) return SealevelriseRigidEnum;
    772772              else if (strcmp(name,"SealevelriseElastic")==0) return SealevelriseElasticEnum;
     773              else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum;
     774              else if (strcmp(name,"SealevelriseTidalLoveH")==0) return SealevelriseTidalLoveHEnum;
     775              else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum;
     776              else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum;
     777              else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum;
     778              else if (strcmp(name,"SealevelrisePolarMoi")==0) return SealevelrisePolarMoiEnum;
     779              else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
    773780              else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum;
    774781              else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum;
    … …  
    868875              else if (strcmp(name,"HydrologyDCInefficientAnalysis")==0) return HydrologyDCInefficientAnalysisEnum;
    869876              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    870               else if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"HydrologySommersAnalysis")==0) return HydrologySommersAnalysisEnum;
    871881              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    872882              else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
    … …  
    875885              else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
    876886              else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
     887              else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum;
    881888              else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum;
    882889              else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum;

  • issm/branches/trunk-larour-NatClimateChange2016/src/m/classes/slr.m

    r21049 r21392  
    1616                tide_love_k    = 0; %ideam
    1717                tide_love_h    = 0; %ideam
     18                fluid_love     = 0;
     19                equatorial_moi = 0;
     20                polar_moi               = 0;
     21                angular_velocity = 0;
    1822                rigid          = 0;
    1923                elastic        = 0;
    … …  
    4448                self.rigid=1;
    4549                self.elastic=1;
    46                 self.rotation=1;
     50                self.rotation=0;
    4751
    4852                %tidal love numbers:
    4953                self.tide_love_h=0.6149; %degree 2
    5054                self.tide_love_k=0.3055; % degree 2
     55       
     56                %secular fluid love number:
     57                self.fluid_love=0.942;
     58               
     59                %moment of inertia:
     60                self.equatorial_moi=8.0077*10^37; % [kg m^2]
     61                self.polar_moi           =8.0345*10^37; % [kg m^2]
     62
     63                % mean rotational velocity of earth
     64                self.angular_velocity=7.2921*10^-5; % [s^-1]
    5165
    5266                %numerical discretization accuracy
    … …  
    7084                        md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1);
    7185                        md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1);
     86                        md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1);
     87                        md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1);
     88                        md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1);
     89                        md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1);
    7290                        md = checkfield(md,'fieldname','slr.reltol','size',[1 1]);
    7391                        md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
    … …  
    106124                        fielddisplay(self,'tide_love_k','tidal load Love number (deg 2)');
    107125                        fielddisplay(self,'tide_love_h','tidal load Love number (deg 2)');
    108                         fielddisplay(self,'rotation','earth rotational potential perturbation');
     126                        fielddisplay(self,'fluid_love','secular fluid Love number');
     127                        fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]');
     128                        fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]');
     129                        fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]');
    109130                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
    110131                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
     132                        fielddisplay(self,'rotation','earth rotational potential perturbation');
    111133                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
    112134                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    … …  
    125147                        WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double');
    126148                        WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double');
     149                        WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double');
     150                        WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double');
     151                        WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double');
     152                        WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double');
    127153                        WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean');
    128154                        WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean');
    … …  
    153179                        writejsdouble(fid,[modelname '.slr.tide_love_k'],self.tide_love_k);
    154180                        writejsdouble(fid,[modelname '.slr.tide_love_h'],self.tide_love_h);
     181                        writejsdouble(fid,[modelname '.slr.fluid_love'],self.fluid_love);
     182                        writejsdouble(fid,[modelname '.slr.equatorial_moi'],self.equatorial_moi);
     183                        writejsdouble(fid,[modelname '.slr.polar_moi'],self.polar_moi);
     184                        writejsdouble(fid,[modelname '.slr.angular_velocity'],self.angular_velocity);
    155185                        writejsdouble(fid,[modelname '.slr.rigid'],self.rigid);
     186                        writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
    156187                        writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
    157                         writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
    158188                        writejsdouble(fid,[modelname '.slr.degacc'],self.degacc);
    159189                        writejscellstring(fid,[modelname '.slr.requested_outputs'],self.requested_outputs);
Note: See TracChangeset for help on using the changeset viewer.