source: issm/oecreview/Archive/24684-25833/ISSM-24939-24940.diff@ 25834

Last change on this file since 25834 was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 52.0 KB

  • ../trunk-jpl/src/m/classes/slr.m

     
    3232                horiz                  = 0;
    3333                Ngia                   = NaN;
    3434                Ugia                   = NaN;
     35                eartharea              = 4*pi*planetradius('earth')^2;
    3536
    3637                requested_outputs      = {};
    3738                transitions            = {};
    … …  
    101102                                return;
    102103                        end
    103104
    104                         md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
     105                        %md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements 1]);
     106                        md = checkfield(md,'fieldname','slr.deltathickness','timeseries',1,'NaN',1,'Inf',1);
    105107                        md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    106108                        md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1);
    107109                        md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1);
    … …  
    131133                        end
    132134
    133135                        %cross check that whereever we have an ice load, the mask is <0 on each vertex:
    134                         pos=find(self.deltathickness);
    135                         maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
    136                         [els,vertices]=find(maskpos>0);
    137                         if length(els),
    138                                 warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
    139                         end
     136                        %pos=find(self.deltathickness);
     137                        %maskpos=md.mask.ice_levelset(md.mesh.elements(pos,:));
     138                        %[els,vertices]=find(maskpos>0);
     139                        %if length(els),
     140                        %       warning('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
     141                        %end
    140142
    141143                        %check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
    142144                        %a coupler to a planet model is provided.
    … …  
    220222                        WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
    221223                        WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
    222224                        WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
     225                        WriteData(fid,prefix,'object',self,'fieldname','eartharea','format','Double');
    223226                       
    224227                        %process requested outputs
    225228                        outputs = self.requested_outputs;

  • ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

     
    320320        parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
    321321        parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
    322322        parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
     323        parameters->AddObject(iomodel->CopyConstantObject("md.slr.eartharea",SealevelEarthAreaEnum));
    323324
    324325        /*Deal with dsl multi-model ensembles: {{{*/
    325326        iomodel->FetchData(&dslmodel,"md.dsl.model");

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

     
    301301syn keyword cConstant ResultsEnum
    302302syn keyword cConstant RootPathEnum
    303303syn keyword cConstant SaveResultsEnum
     304syn keyword cConstant SealevelEarthAreaEnum
    304305syn keyword cConstant SealevelEustaticEnum
    305306syn keyword cConstant SealevelriseAbstolEnum
    306307syn keyword cConstant SealevelriseAngularVelocityEnum
    … …  
    14301431syn keyword cType Results
    14311432syn keyword cType RiftStruct
    14321433syn keyword cType Riftfront
     1434syn keyword cType SealevelMasks
    14331435syn keyword cType Seg
    14341436syn keyword cType SegInput2
    14351437syn keyword cType SegRef

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

     
    295295        ResultsEnum,
    296296        RootPathEnum,
    297297        SaveResultsEnum,
     298        SealevelEarthAreaEnum,
    298299        SealevelEustaticEnum,
    299300        SealevelriseAbstolEnum,
    300301        SealevelriseAngularVelocityEnum,

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

     
    303303                case ResultsEnum : return "Results";
    304304                case RootPathEnum : return "RootPath";
    305305                case SaveResultsEnum : return "SaveResults";
     306                case SealevelEarthAreaEnum : return "SealevelEarthArea";
    306307                case SealevelEustaticEnum : return "SealevelEustatic";
    307308                case SealevelriseAbstolEnum : return "SealevelriseAbstol";
    308309                case SealevelriseAngularVelocityEnum : return "SealevelriseAngularVelocity";

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

     
    309309              else if (strcmp(name,"Results")==0) return ResultsEnum;
    310310              else if (strcmp(name,"RootPath")==0) return RootPathEnum;
    311311              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
     312              else if (strcmp(name,"SealevelEarthArea")==0) return SealevelEarthAreaEnum;
    312313              else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum;
    313314              else if (strcmp(name,"SealevelriseAbstol")==0) return SealevelriseAbstolEnum;
    314315              else if (strcmp(name,"SealevelriseAngularVelocity")==0) return SealevelriseAngularVelocityEnum;
    … …  
    381382              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
    382383              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    383384              else if (strcmp(name,"SmbPfac")==0) return SmbPfacEnum;
    384               else if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
     388              if (strcmp(name,"SmbRdl")==0) return SmbRdlEnum;
     389              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
    389390              else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum;
    390391              else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum;
    391392              else if (strcmp(name,"SmbRunoffalti")==0) return SmbRunoffaltiEnum;
    … …  
    504505              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    505506              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    506507              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    507               else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
     511              if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     512              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
    512513              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    513514              else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
    514515              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    … …  
    627628              else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    628629              else if (strcmp(name,"Ice")==0) return IceEnum;
    629630              else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
    630               else if (strcmp(name,"Input")==0) return InputEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
     634              if (strcmp(name,"Input")==0) return InputEnum;
     635              else if (strcmp(name,"InversionCostFunctionsCoefficients")==0) return InversionCostFunctionsCoefficientsEnum;
    635636              else if (strcmp(name,"InversionSurfaceObs")==0) return InversionSurfaceObsEnum;
    636637              else if (strcmp(name,"InversionThicknessObs")==0) return InversionThicknessObsEnum;
    637638              else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum;
    … …  
    750751              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    751752              else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum;
    752753              else if (strcmp(name,"SmbMassBalanceClimate")==0) return SmbMassBalanceClimateEnum;
    753               else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
     757              if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
     758              else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
    758759              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
    759760              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    760761              else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    … …  
    873874              else if (strcmp(name,"Outputdefinition10")==0) return Outputdefinition10Enum;
    874875              else if (strcmp(name,"Outputdefinition11")==0) return Outputdefinition11Enum;
    875876              else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
    876               else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
     880              if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
     881              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    881882              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    882883              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    883884              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    … …  
    996997              else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum;
    997998              else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum;
    998999              else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
    999               else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
     1003              if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
     1004              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    10041005              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    10051006              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    10061007              else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
    … …  
    11191120              else if (strcmp(name,"Hydrologyshakti")==0) return HydrologyshaktiEnum;
    11201121              else if (strcmp(name,"Hydrologyshreve")==0) return HydrologyshreveEnum;
    11211122              else if (strcmp(name,"IceMass")==0) return IceMassEnum;
    1122               else if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
     1126              if (strcmp(name,"IceMassScaled")==0) return IceMassScaledEnum;
     1127              else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    11271128              else if (strcmp(name,"IceVolumeAboveFloatationScaled")==0) return IceVolumeAboveFloatationScaledEnum;
    11281129              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    11291130              else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
    … …  
    12421243              else if (strcmp(name,"Paterson")==0) return PatersonEnum;
    12431244              else if (strcmp(name,"Pengrid")==0) return PengridEnum;
    12441245              else if (strcmp(name,"Penpair")==0) return PenpairEnum;
    1245               else if (strcmp(name,"Penta")==0) return PentaEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
     1249              if (strcmp(name,"Penta")==0) return PentaEnum;
     1250              else if (strcmp(name,"PentaInput")==0) return PentaInputEnum;
    12501251              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    12511252              else if (strcmp(name,"ProfilingCurrentFlops")==0) return ProfilingCurrentFlopsEnum;
    12521253              else if (strcmp(name,"ProfilingCurrentMem")==0) return ProfilingCurrentMemEnum;

  • ../trunk-jpl/src/c/cores/cores.h

     
    5757void steric_core(FemModel* femmodel);
    5858void sealevelrise_core_geometry(FemModel* femmodel);
    5959SealevelMasks* sealevelrise_core_masks(FemModel* femmodel);
    60 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* peartharea,IssmDouble* poceanarea);
    61 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea);
    62 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks);
     60Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* mask, IssmDouble* poceanarea);
     61Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea);
     62void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks);
    6363void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel*  femmodel,Vector<IssmDouble>* RSLg);
    6464void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);
    6565IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);

  • ../trunk-jpl/src/c/cores/sealevelrise_core.cpp

     
    9595        int  horiz;
    9696        int  geodetic=0;
    9797        IssmDouble          dt;
    98         IssmDouble          eartharea,oceanarea;
     98        IssmDouble          oceanarea;
    9999
    100100        /*Should we even be here?:*/
    101101        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
    … …  
    152152                masks=sealevelrise_core_masks(femmodel);
    153153
    154154                /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
    155                 RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&eartharea,&oceanarea);
     155                RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea);
    156156
    157157                /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
    158                 RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,eartharea,oceanarea);
     158                RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea);
    159159
    160160                /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
    161                 sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea,masks);
     161                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
    162162
    163163                /*compute viscosus (GIA) geodetic signatures:*/
    164164                sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
    … …  
    353353
    354354
    355355}/*}}}*/
    356 Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/
     356Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel,SealevelMasks* masks, IssmDouble* poceanarea){ /*{{{*/
    357357
    358358        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    359359
    … …  
    367367        IssmDouble *longitude = NULL;
    368368        IssmDouble *radius    = NULL;
    369369        int         loop;
    370         IssmDouble eartharea,oceanarea;
     370        IssmDouble oceanarea;
    371371
    372372        /*outputs:*/
    373373        IssmDouble eustatic;
    … …  
    387387        RSLgi = new Vector<IssmDouble>(gsize);
    388388
    389389        /*call the eustatic main module: */
    390         femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes
     390        femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks, latitude, longitude, radius,loop); //this computes
    391391
    392392        /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    393393        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
    … …  
    405405        xDelete<IssmDouble>(radius);
    406406
    407407        /*Assign output pointers and return: */
    408         *peartharea=eartharea;
    409408        *poceanarea=oceanarea;
    410409        return RSLgi;
    411410}/*}}}*/
    412 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/
     411Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel, SealevelMasks* masks, Vector<IssmDouble>* RSLg_eustatic,IssmDouble oceanarea){ /*{{{*/
    413412
    414413        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    415414          non eustatic core of the SLR solution */
    … …  
    478477                RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
    479478
    480479                /*call the non eustatic module: */
    481                 femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  eartharea, masks, latitude,longitude, radius,verboseconvolution,loop);
     480                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old,  masks, latitude,longitude, radius,verboseconvolution,loop);
    482481
    483482                /*assemble solution vector: */
    484483                RSLgo->Assemble();
    … …  
    487486
    488487                        /*call rotational feedback  module: */
    489488                        RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
    490                         femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, eartharea, masks, latitude,longitude,radius);
     489                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks, latitude,longitude,radius);
    491490                        RSLgo_rot->Assemble();
    492491
    493492                        /*save changes in inertia tensor as results: */
    … …  
    536535
    537536        return RSLg;
    538537} /*}}}*/
    539 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg,IssmDouble eartharea, SealevelMasks* masks){ /*{{{*/
     538void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg, SealevelMasks* masks){ /*{{{*/
    540539
    541540        Vector<IssmDouble> *U_esa  = NULL;
    542541        Vector<IssmDouble> *U_north_esa   = NULL;
    … …  
    576575        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    577576
    578577        /*call the elastic main modlule:*/
    579         femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz);
     578        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks, latitude,longitude,radius,xx,yy,zz,loop,horiz);
    580579
    581580        /*Assign output pointers:*/
    582581        *pU_esa=U_esa;

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

     
    208208                #endif
    209209                #ifdef _HAVE_ESA_
    210210                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    211                 void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     211                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    212212                #endif
    213213                #ifdef _HAVE_SEALEVELRISE_
    214214                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    215215                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    216                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
     216                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    217217                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    218                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    219                 void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    220                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     218                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
     219                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     220                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
    221221                #endif
    222222
    223223                /*}}}*/

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

     
    168168#endif
    169169#ifdef _HAVE_ESA_
    170170                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    171                 void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     171                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    172172#endif
    173173#ifdef _HAVE_SEALEVELRISE_
    174                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
     174                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    175175                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    176176                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    177                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    178                 void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    179                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     177                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
     178                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     179                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
    180180                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    181181#endif
    182182

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

     
    174174#endif
    175175#ifdef _HAVE_ESA_
    176176                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){_error_("not implemented yet!");};
    177                 void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){_error_("not implemented yet!");};
     177                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){_error_("not implemented yet!");};
    178178#endif
    179179#ifdef _HAVE_SEALEVELRISE_
    180180                void    SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");};
    181                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea){_error_("not implemented yet!");};
     181                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");};
    182182                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
    183                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){_error_("not implemented yet!");};
    184                 void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){_error_("not implemented yet!");};
    185                 void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){_error_("not implemented yet!");};
     183                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){_error_("not implemented yet!");};
     184                void    SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){_error_("not implemented yet!");};
     185                void    SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){_error_("not implemented yet!");};
    186186                IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");};
    187187#endif
    188188

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

     
    370370                #endif
    371371                #ifdef _HAVE_ESA_
    372372                virtual void          EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast, Vector<IssmDouble>* pX, Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy)=0;
    373                 virtual void          EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea)=0;
     373                virtual void          EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz)=0;
    374374                #endif
    375375                #ifdef _HAVE_SEALEVELRISE_
    376376                virtual void          SetSealevelMasks(SealevelMasks* masks)=0;
    … …  
    377377                virtual IssmDouble    GetArea3D(void)=0;
    378378                virtual IssmDouble    GetAreaSpherical(void)=0;
    379379                virtual IssmDouble    OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0;
    380                 virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble eartharea)=0;
    381                 virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea)=0;
     380                virtual void          SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0;
     381                virtual void          SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea)=0;
    382382                virtual void          SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
    383                 virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea)=0;
    384                 virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz)=0;
     383                virtual void          SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius)=0;
     384                virtual void          SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz)=0;
    385385                #endif
    386386
    387387};

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

     
    159159                #endif
    160160                #ifdef _HAVE_ESA_
    161161                void    EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,  Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy);
    162                 void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea);
     162                void    EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz);
    163163                #endif
    164164                #ifdef _HAVE_SEALEVELRISE_
    165165                IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks);
    166                 void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble eartharea);
     166                void    SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks);
    167167                void    SetSealevelMasks(SealevelMasks* masks);
    168168                void    SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
    169                 void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    170                 void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    171                 void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea);
    172                 void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea);
    173                 void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz);
     169                void    SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
     170                void    SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
     171                void    SealevelriseEustaticBottomPressure(IssmDouble* Sgi, IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea);
     172                void    SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius);
     173                void    SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz);
    174174                #endif
    175175                /*}}}*/
    176176                /*Tria specific routines:{{{*/

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

     
    162162                void EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz);
    163163                #endif
    164164                #ifdef _HAVE_SEALEVELRISE_
    165                 void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
     165                void SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop);
    166166                void SealevelriseGeometry(IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    167                 void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    168                 void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    169                 void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble eartharea,SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
     167                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
     168                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
     169                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
    170170                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg,SealevelMasks* masks, IssmDouble oceanarea);
    171171                #endif
    172172                void HydrologyEPLupdateDomainx(IssmDouble* pEplcount);

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

     
    46344634/*}}}*/
    46354635void FemModel::EsaGeodetic3D(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){/*{{{*/
    46364636
    4637         IssmDouble  eartharea=0;
    4638         IssmDouble  eartharea_cpu=0;
    4639 
    46404637        int         ns,nsmax;
    46414638
    46424639        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    … …  
    46454642        /*First, figure out the surface area of Earth: */
    46464643        for(int i=0;i<ns;i++){
    46474644                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4648                 eartharea_cpu += element->GetAreaSpherical();
    46494645        }
    4650         ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4651         ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    46524646
    46534647        /*Figure out max of ns: */
    46544648        ISSM_MPI_Reduce(&ns,&nsmax,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    … …  
    46584652        for(int i=0;i<nsmax;i++){
    46594653                if(i<ns){
    46604654                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4661                         element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz,eartharea);
     4655                        element->EsaGeodetic3D(pUp,pNorth,pEast,latitude,longitude,radius,xx,yy,zz);
    46624656                }
    46634657                if(i%100==0){
    46644658                        pUp->Assemble();
    … …  
    46934687
    46944688}
    46954689/*}}}*/
    4696 void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
     4690void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
    46974691
    46984692        /*serialized vectors:*/
    46994693        IssmDouble  eustatic       = 0.;
    … …  
    47014695        IssmDouble  eustatic_cpu_e = 0.;
    47024696        IssmDouble  oceanarea      = 0.;
    47034697        IssmDouble  oceanarea_cpu  = 0.;
    4704         IssmDouble  eartharea      = 0.;
    4705         IssmDouble  eartharea_cpu = 0.;
    47064698
    47074699   /*Initialize temporary vector that will be used to sum eustatic components
    47084700    * on all local elements, prior to assembly:*/
    … …  
    47154707        for(int i=0;i<elements->Size();i++){
    47164708                Element*   element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    47174709                IssmDouble area=element->GetAreaSpherical();
    4718                 eartharea_cpu += area;
    47194710                if (masks->isoceanin[i]) oceanarea_cpu += area;
    47204711        }
    47214712        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    47224713        ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    47234714        _assert_(oceanarea>0.);
    47244715
    4725         ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    4726         ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    4727 
    47284716        /*Call the sea level rise core: */
    47294717        for(int i=0;i<elements->Size();i++){
    47304718                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4731                 element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea,eartharea);
     4719                element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, latitude,longitude,radius,oceanarea);
    47324720                eustatic_cpu+=eustatic_cpu_e;
    47334721        }
    47344722
    … …  
    47464734        xDelete<IssmDouble>(RSLgi);
    47474735
    47484736        /*Assign output pointers:*/
    4749         *peartharea = eartharea;
    47504737        *poceanarea = oceanarea;
    47514738        *peustatic  = eustatic;
    47524739
    47534740}
    47544741/*}}}*/
    4755 void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
     4742void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
    47564743
    47574744        /*serialized vectors:*/
    47584745        IssmDouble* RSLg_old=NULL;
    … …  
    47734760        /*Call the sea level rise core: */
    47744761        for(int i=0;i<elements->Size();i++){
    47754762                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4776                 element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius,eartharea);
     4763                element->SealevelriseNonEustatic(RSLgo,RSLg_old,masks, latitude,longitude,radius);
    47774764        }
    47784765        pRSLgo->SetValues(gsize,indices,RSLgo,ADD_VAL);
    47794766        pRSLgo->Assemble();
    … …  
    47854772
    47864773}
    47874774/*}}}*/
    4788 void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
     4775void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
    47894776
    47904777        /*serialized vectors:*/
    47914778        IssmDouble* RSLg_old=NULL;
    … …  
    48014788        IssmDouble moi_list_cpu[3]={0,0,0};
    48024789        for(int i=0;i<elements->Size();i++){
    48034790                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4804                 element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks, eartharea);
     4791                element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,masks );
    48054792                moi_list_cpu[0] += moi_list[0];
    48064793                moi_list_cpu[1] += moi_list[1];
    48074794                moi_list_cpu[2] += moi_list[2];
    … …  
    48614848
    48624849}
    48634850/*}}}*/
    4864 void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble eartharea, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
     4851void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, SealevelMasks* masks, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
    48654852
    48664853        /*serialized vectors:*/
    48674854        IssmDouble* RSLg=NULL;
    … …  
    48874874        /*Call the sea level rise core: */
    48884875        for(int i=0;i<elements->Size();i++){
    48894876                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    4890                 element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
     4877                element->SealevelriseGeodetic(Up,North,East,RSLg,masks, latitude,longitude,radius,xx,yy,zz,horiz);
    48914878        }
    48924879
    48934880        pUp->SetValues(gsize,indices,Up,ADD_VAL);

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

     
    52885288        return;
    52895289}
    52905290/*}}}*/
    5291 void    Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/
     5291void    Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz){ /*{{{*/
    52925292
    52935293        /*diverse:*/
    52945294        int gsize;
    … …  
    52955295        bool spherical=true;
    52965296        IssmDouble llr_list[NUMVERTICES][3];
    52975297        IssmDouble xyz_list[NUMVERTICES][3];
    5298         IssmDouble area;
     5298        IssmDouble area,eartharea;
    52995299        IssmDouble I;           //ice/water loading
    53005300        IssmDouble late,longe,re;
    53015301        IssmDouble lati,longi,ri;
    … …  
    53275327        /*recover material parameters: */
    53285328        rho_ice=FindParam(MaterialsRhoIceEnum);
    53295329        rho_earth=FindParam(MaterialsEarthDensityEnum);
     5330       
     5331        /*recover earth area: */
     5332        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
    53305333
    53315334        /*how many dofs are we working with here? */
    53325335        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    … …  
    54645467
    54655468}
    54665469/*}}}*/
    5467 void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks, IssmDouble eartharea){/*{{{*/
     5470void    Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/
    54685471        /*early return if we are not on an ice cap OR ocean:*/
    54695472        if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){
    54705473                dI_list[0] = 0.0; // this is important!!!
    … …  
    54745477        }
    54755478
    54765479        /*Compute area of element:*/
    5477         IssmDouble area;
     5480        IssmDouble area,eartharea;
    54785481        area=GetAreaSpherical();
    54795482
     5483        /*recover earth area: */
     5484        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     5485
    54805486        /*Compute lat,long,radius of elemental centroid: */
    54815487        bool spherical=true;
    54825488        IssmDouble llr_list[NUMVERTICES][3];
    … …  
    56605666        return;
    56615667}
    56625668/*}}}*/
    5663 void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5669void    Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
    56645670
    56655671        /*Computational flags:*/
    56665672        int bp_compute_fingerprints= 0;
    … …  
    56715677        if(!masks->isoceanin[this->lid]){
    56725678                /*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
    56735679                 *bottom pressure fingerprints:*/
    5674                 if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
     5680                if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(Sgi,peustatic,latitude,longitude,radius,oceanarea);
    56755681        }
    56765682        //if(!IsIceInElement()){
    56775683                /*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
    5678                 this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea,eartharea);
     5684                this->SealevelriseEustaticIce(Sgi,peustatic,masks, latitude,longitude,radius,oceanarea);
    56795685        //}
    56805686
    56815687}
    56825688/*}}}*/
    5683 void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5689void    Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
    56845690
    56855691        /*diverse:*/
    56865692        int gsize,dummy;
    56875693        bool spherical=true;
    56885694        IssmDouble llr_list[NUMVERTICES][3];
    5689         IssmDouble area;
     5695        IssmDouble area,eartharea;
    56905696        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    56915697        IssmDouble rho;
    56925698        IssmDouble late,longe,re;
    … …  
    57445750        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    57455751        this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);
    57465752
     5753        /*recover earth area: */
     5754        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     5755
    57475756        /*recover precomputed green function kernels:*/
    57485757        DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter);
    57495758        parameter->GetParameterValueByPointer(&G_rigid_precomputed,&M);
    … …  
    58355844        return;
    58365845}
    58375846/*}}}*/
    5838 void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
     5847void    Tria::SealevelriseEustaticBottomPressure(IssmDouble* Sgi,IssmDouble* peustatic, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea){ /*{{{*/
    58395848
    58405849        /*diverse:*/
    58415850        int gsize;
    58425851        bool spherical=true;
    58435852        IssmDouble llr_list[NUMVERTICES][3];
    5844         IssmDouble area;
     5853        IssmDouble area,eartharea;
    58455854        IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
    58465855        IssmDouble rho;
    58475856        IssmDouble late,longe,re;
    … …  
    58755884        rho_water=FindParam(MaterialsRhoSeawaterEnum);
    58765885        rho_earth=FindParam(MaterialsEarthDensityEnum);
    58775886
     5887        /*recover earth area: */
     5888        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     5889
    58785890        /*recover love numbers and computational flags: */
    58795891        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    58805892        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    … …  
    59815993        return;
    59825994}
    59835995/*}}}*/
    5984 void    Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/
     5996void    Tria::SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius){ /*{{{*/
    59855997
    59865998        /*diverse:*/
    59875999        int gsize,dummy;
    59886000        bool spherical=true;
    59896001        IssmDouble llr_list[NUMVERTICES][3];
    5990         IssmDouble area;
     6002        IssmDouble area,eartharea;
    59916003        IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
    59926004        IssmDouble late,longe;
    59936005        IssmDouble lati,longi,ri;
    … …  
    60256037        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    60266038        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    60276039
     6040        /*recover earth area: */
     6041        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     6042
    60286043        /*early return if rigid or elastic not requested:*/
    60296044        if(!computerigid && !computeelastic) return;
    60306045
    … …  
    60806095        return;
    60816096}
    60826097/*}}}*/
    6083 void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/
     6098void    Tria::SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks, IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,int horiz){ /*{{{*/
    60846099
    60856100        /*diverse:*/
    60866101        int gsize,dummy;
    … …  
    60876102        bool spherical=true;
    60886103        IssmDouble llr_list[NUMVERTICES][3];
    60896104        IssmDouble xyz_list[NUMVERTICES][3];
    6090         IssmDouble area;
     6105        IssmDouble area,eartharea;
    60916106        IssmDouble I, S;                //change in relative ice thickness and sea level
    60926107        IssmDouble late,longe,re;
    60936108        IssmDouble lati,longi,ri;
    … …  
    61406155        /*how many dofs are we working with here? */
    61416156        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    61426157
     6158        /*recover earth area: */
     6159        this->parameters->FindParam(&eartharea,SealevelEarthAreaEnum);
     6160
    61436161        /*retrieve indices:*/
    61446162        if(computerigid){this->inputs2->GetArrayPtr(SealevelriseIndicesEnum,this->lid,&indices,&dummy); _assert_(dummy==gsize);}
    61456163
Note: See TracBrowser for help on using the repository browser.