Changeset 22160


Ignore:
Timestamp:
10/10/17 11:14:30 (7 years ago)
Author:
Eric.Larour
Message:

CHG: diverse bug fixes

Location:
issm/branches/trunk-larour-NatGeoScience2016/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/analyses/SealevelriseAnalysis.cpp

    r22159 r22160  
    5858                        Element* element=(Element*)elements->GetObjectByOffset(counter);
    5959                        element->InputUpdateFromConstant(0.0,SealevelriseCumDeltathicknessEnum);
    60                         element->InputUpdateFromConstant(0.0,SealevelRateEnum);
    61                         element->InputUpdateFromConstant(0.0,SealevelUmotionRateEnum);
     60                        element->InputUpdateFromConstant(0.0,SealevelNEsaRateEnum);
     61                        element->InputUpdateFromConstant(0.0,SealevelNGiaRateEnum);
     62                        element->InputUpdateFromConstant(0.0,SealevelUEsaRateEnum);
     63                        element->InputUpdateFromConstant(0.0,SealevelUGiaRateEnum);
    6264                        counter++;
    6365                }

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/Elements/Element.cpp

    r22155 r22160  
    16411641                                name==MaterialsRheologyEsbarEnum ||
    16421642                                name==SealevelEnum ||
    1643                                 name==SealevelUmotionEnum ||
    1644                                 name==SealevelUmotionRateEnum ||
    1645                                 name==SealevelNmotionEnum ||
    1646                                 name==SealevelEmotionEnum ||
    1647                                 name==SealevelAbsoluteEnum ||
    1648                                 name==SealevelEustaticEnum ||
     1643                                name==SealevelUEsaEnum ||
     1644                                name==SealevelUEsaRateEnum ||
     1645                                name==SealevelNEsaEnum ||
     1646                                name==SealevelNEsaRateEnum ||
     1647                                name==SealevelUNorthEsaEnum ||
     1648                                name==SealevelUEastEsaEnum ||
     1649                                name==SealevelRSLEustaticEnum ||
     1650                                name==SealevelRSLEustaticRateEnum ||
     1651                                name==SealevelRSLEnum ||
     1652                                name==SealevelRSLRateEnum ||
     1653                                name==SealevelUGiaEnum ||
     1654                                name==SealevelUGiaRateEnum ||
     1655                                name==SealevelNGiaEnum ||
     1656                                name==SealevelNGiaRateEnum ||
    16491657                                name==SealevelriseDeltathicknessEnum ||
    16501658                                name==SealevelriseCumDeltathicknessEnum ||
    1651                                 name==SealevelRateEnum ||
    16521659                                name==EsaUmotionEnum ||
    16531660                                name==EsaNmotionEnum ||

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.cpp

    r22159 r22160  
    24222422#endif
    24232423#ifdef _HAVE_SEALEVELRISE_
    2424 void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
     2424void FemModel::SealevelriseEustatic(Vector<IssmDouble>* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/
    24252425
    24262426        /*serialized vectors:*/
    … …  
    24612461               
    24622462                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2463                         element->SealevelriseEustatic(pSgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
     2463                        element->SealevelriseEustatic(pRSLgi,&eustatic_cpu_e,latitude,longitude,radius,oceanarea,eartharea);
    24642464                        eustatic_cpu+=eustatic_cpu_e;
    24652465                }
    2466                 if(i%loop==0)pSgi->Assemble();
     2466                if(i%loop==0)pRSLgi->Assemble();
    24672467        }
    24682468        if(VerboseConvergence())_printf0_("\n");
    24692469               
    24702470        /*One last time: */
    2471         pSgi->Assemble();
     2471        pRSLgi->Assemble();
    24722472
    24732473        /*Sum all eustatic components from all cpus:*/
    … …  
    24812481}
    24822482/*}}}*/
    2483 void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
     2483void FemModel::SealevelriseNonEustatic(Vector<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/
    24842484
    24852485        /*serialized vectors:*/
    2486         IssmDouble* Sg_old=NULL;
     2486        IssmDouble* RSLg_old=NULL;
    24872487       
    24882488        IssmDouble  eartharea=0;
    … …  
    24922492       
    24932493        /*Serialize vectors from previous iteration:*/
    2494         Sg_old=pSg_old->ToMPISerial();
     2494        RSLg_old=pRSLg_old->ToMPISerial();
    24952495
    24962496        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    … …  
    25152515                        if(verboseconvolution)if(VerboseConvergence())if(i%100==0)_printf_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%   ");
    25162516                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2517                         element->SealevelriseNonEustatic(pSgo,Sg_old,latitude,longitude,radius,eartharea);
    2518                 }
    2519                 if(i%loop==0)pSgo->Assemble();
     2517                        element->SealevelriseNonEustatic(pRSLgo,RSLg_old,latitude,longitude,radius,eartharea);
     2518                }
     2519                if(i%loop==0)pRSLgo->Assemble();
    25202520        }
    25212521        if(verboseconvolution)if(VerboseConvergence())_printf_("\n");
    25222522       
    25232523        /*Free ressources:*/
    2524         xDelete<IssmDouble>(Sg_old);
    2525 }
    2526 /*}}}*/
    2527 void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
     2524        xDelete<IssmDouble>(RSLg_old);
     2525}
     2526/*}}}*/
     2527void FemModel::SealevelriseRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/
    25282528
    25292529        /*serialized vectors:*/
    2530         IssmDouble* Sg_old=NULL;
     2530        IssmDouble* RSLg_old=NULL;
    25312531        IssmDouble  eartharea=0;
    25322532        IssmDouble  eartharea_cpu=0;
    … …  
    25372537
    25382538        /*Serialize vectors from previous iteration:*/
    2539         Sg_old=pSg_old->ToMPISerial();
     2539        RSLg_old=pRSLg_old->ToMPISerial();
    25402540
    25412541        /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */
    … …  
    25512551        for(int i=0;i<elements->Size();i++){
    25522552                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2553                 element->SealevelriseMomentOfInertia(&moi_list[0],Sg_old,eartharea);
     2553                element->SealevelriseMomentOfInertia(&moi_list[0],RSLg_old,eartharea);
    25542554                moi_list_cpu[0] += moi_list[0];
    25552555                moi_list_cpu[1] += moi_list[1];
    … …  
    25952595                                                (-m3/6.0 + 0.5*m3*cos(2.0*lati) - 0.5*sin(2.*lati)*(m1*cos(longi)+m2*sin(longi)));
    25962596       
    2597                 pSgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
     2597                pRSLgo_rot->SetValue(sid,value,INS_VAL); //INS_VAL ensures that you don't add several times
    25982598        }
    25992599
    26002600        /*Assemble mesh velocity*/
    2601         pSgo_rot->Assemble();
     2601        pRSLgo_rot->Assemble();
    26022602       
    26032603        /*Free ressources:*/
    2604         xDelete<IssmDouble>(Sg_old);
    2605        
    2606 }
    2607 /*}}}*/
    2608 void FemModel::SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
     2604        xDelete<IssmDouble>(RSLg_old);
     2605       
     2606}
     2607/*}}}*/
     2608void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
    26092609
    26102610        /*serialized vectors:*/
    2611         IssmDouble* Sg=NULL;
     2611        IssmDouble* RSLg=NULL;
    26122612       
    26132613        IssmDouble  eartharea=0;
    … …  
    26172617       
    26182618        /*Serialize vectors from previous iteration:*/
    2619         Sg=pSg->ToMPISerial();
     2619        RSLg=pRSLg->ToMPISerial();
    26202620
    26212621        /*Go through elements, and add contribution from each element to the deflection vector wg:*/
    … …  
    26392639                        if(VerboseConvergence())if(i%100==0)_printf0_("\r" << "      convolution progress: " << (double)i/(double)ns*100 << "%  ");
    26402640                        Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2641                         element->SealevelriseGeodetic(pUp,pNorth,pEast,Sg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
     2641                        element->SealevelriseGeodetic(pUp,pNorth,pEast,RSLg,latitude,longitude,radius,xx,yy,zz,eartharea,horiz);
    26422642                }
    26432643                if(i%loop==0){
    … …  
    26592659
    26602660        /*Free ressources:*/
    2661         xDelete<IssmDouble>(Sg);
    2662 }
    2663 /*}}}*/
    2664 IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* Sg) { /*{{{*/
    2665 
    2666         IssmDouble* Sg_serial=NULL;
     2661        xDelete<IssmDouble>(RSLg);
     2662}
     2663/*}}}*/
     2664IssmDouble FemModel::SealevelriseOceanAverage(Vector<IssmDouble>* RSLg) { /*{{{*/
     2665
     2666        IssmDouble* RSLg_serial=NULL;
    26672667        IssmDouble  oceanvalue,oceanvalue_cpu;
    26682668        IssmDouble  oceanarea,oceanarea_cpu;
    26692669
    26702670        /*Serialize vectors from previous iteration:*/
    2671         Sg_serial=Sg->ToMPISerial();
     2671        RSLg_serial=RSLg->ToMPISerial();
    26722672
    26732673        /*Initialize:*/
    … …  
    26792679                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    26802680                oceanarea_cpu += element->OceanArea();
    2681                 oceanvalue_cpu += element->OceanAverage(Sg_serial);
     2681                oceanvalue_cpu += element->OceanAverage(RSLg_serial);
    26822682        }
    26832683        ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
    … …  
    26882688
    26892689        /*Free ressources:*/
    2690         xDelete<IssmDouble>(Sg_serial);
     2690        xDelete<IssmDouble>(RSLg_serial);
    26912691       
    26922692        return oceanvalue/oceanarea;

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/classes/FemModel.h

    r22114 r22160  
    129129                void SealevelriseNonEustatic(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop);
    130130                void SealevelriseRotationalFeedback(Vector<IssmDouble>* pSgo_rot, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius);
    131                 void SealevelriseGeodetic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
     131                void SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz);
    132132                IssmDouble SealevelriseOceanAverage(Vector<IssmDouble>* Sg);
    133133                #endif

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/cores.h

    r22155 r22160  
    5454void steric_core(FemModel* femmodel);
    5555Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
    56 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
    57 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* Sg);
     56Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic);
     57void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* RSLg);
     58void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia,Vector<IssmDouble>** pN_gia,FemModel*  femmodel,Vector<IssmDouble>* RSLg);
     59void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg);
    5860IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
    5961
    … …  
    6971void TransferSealevel(FemModel* femmodel,int forcingenum);
    7072void EarthMassTransport(FemModel* femmodel);
    71 void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs);
     73void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    7274
    7375//solution configuration

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp

    r22159 r22160  
    6161
    6262        /*variables:*/
    63         Vector<IssmDouble> *Sg    = NULL;
    64         Vector<IssmDouble> *Sg_rate    = NULL;
    65         Vector<IssmDouble> *Sg_eustatic  = NULL;
    66         Vector<IssmDouble> *U_radial  = NULL;
    67         Vector<IssmDouble> *U_radial_rate  = NULL;
    68         Vector<IssmDouble> *U_north   = NULL;
    69         Vector<IssmDouble> *U_east    = NULL;
    70         Vector<IssmDouble>* cumdeltathickness=NULL;
     63        Vector<IssmDouble> *RSLg    = NULL;
     64        Vector<IssmDouble> *RSLg_rate    = NULL;
     65        Vector<IssmDouble> *RSLg_eustatic  = NULL;
     66        Vector<IssmDouble> *RSLg_eustatic_rate  = NULL;
     67        Vector<IssmDouble> *U_esa  = NULL;
     68        Vector<IssmDouble> *U_esa_rate  = NULL;
     69        Vector<IssmDouble> *N_esa  = NULL;
     70        Vector<IssmDouble> *N_esa_rate  = NULL;
     71        Vector<IssmDouble> *U_north_esa   = NULL;
     72        Vector<IssmDouble> *U_east_esa    = NULL;
     73        Vector<IssmDouble> *N_gia= NULL;
     74        Vector<IssmDouble> *U_gia= NULL;
     75        Vector<IssmDouble> *N_gia_rate= NULL;
     76        Vector<IssmDouble> *U_gia_rate= NULL;
    7177
    7278        /*parameters:*/
    7379        bool iscoupler;
     80        int  solution_type;
    7481        int  configuration_type;
    7582        int  modelid,earthid;
    … …  
    93100        femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
    94101        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     102        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    95103
    96104        if(iscoupler){
    … …  
    132140
    133141                /*call eustatic core  (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
    134                 Sg_eustatic=sealevelrise_core_eustatic(femmodel);
     142                RSLg_eustatic=sealevelrise_core_eustatic(femmodel);
    135143
    136144                /*call non-eustatic core (ocean loading tems  - 2nd and 5th terms on the RHS of Farrel and Clark) */
    137                 Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic);
    138                
     145                RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic);
    139146
    140147                /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
    141                 sealevelrise_core_elastic(&U_radial,&U_north,&U_east,femmodel,Sg);
     148                sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg);
     149
     150                /*compute viscosus (GIA) geodetic signatures:*/
     151                sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
     152
     153                /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
     154                sealevelrise_diagnostics(femmodel,RSLg);
     155
     156                /*recover N_esa  = U_esa + RSLg:*/
     157                N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
    142158               
    143159                /*transform these values into rates (as we only run this once each frequency turn:*/
    144                 Sg_rate=Sg->Duplicate(); Sg->Copy(Sg_rate); Sg_rate->Scale(1/(dt*frequency));
    145                 U_radial_rate=U_radial->Duplicate(); U_radial->Copy(U_radial_rate); U_radial_rate->Scale(1/(dt*frequency));
    146                
    147 
    148                 /*get some results into elements:*/
    149                 InputUpdateFromVectorx(femmodel,Sg_rate,SealevelRateEnum,VertexSIdEnum);
    150                 InputUpdateFromVectorx(femmodel,U_radial_rate,SealevelUmotionRateEnum,VertexSIdEnum);
    151                 if (horiz){
    152                         InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion
    153                         InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion
    154                 }
     160                N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
     161                U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
     162                N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
     163                U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
     164                RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
     165                RSLg_eustatic_rate=RSLg_eustatic->Duplicate(); RSLg_eustatic->Copy(RSLg_eustatic_rate); RSLg_eustatic_rate->Scale(1/(dt*frequency));
     166               
     167                /*get some results into elements:{{{*/
     168                if(solution_type==SealevelriseSolutionEnum){
     169                        InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum);
     170                        InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum);
     171                        InputUpdateFromVectorx(femmodel,U_gia,SealevelUGiaEnum,VertexSIdEnum);
     172                        InputUpdateFromVectorx(femmodel,N_gia,SealevelNGiaEnum,VertexSIdEnum);
     173                        InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum);
     174                        InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum);
     175                        if (horiz){
     176                                InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum);       // north motion
     177                                InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEsaEnum,VertexSIdEnum);             // east motion
     178                        }
     179                }
     180                else{
     181                        InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
     182                        InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum);
     183                        InputUpdateFromVectorx(femmodel,U_gia_rate,SealevelUGiaRateEnum,VertexSIdEnum);
     184                        InputUpdateFromVectorx(femmodel,N_gia_rate,SealevelNGiaRateEnum,VertexSIdEnum);
     185                        InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum);
     186                        InputUpdateFromVectorx(femmodel,RSLg_eustatic_rate,SealevelRSLEustaticRateEnum,VertexSIdEnum);
     187                } /*}}}*/
    155188        }
    156189
    … …  
    158191        if(iscoupler){
    159192                /*transfer sea level back to ice caps:*/
    160                 TransferSealevel(femmodel,SealevelRateEnum);
    161                 TransferSealevel(femmodel,SealevelUmotionRateEnum);
     193                TransferSealevel(femmodel,SealevelNEsaRateEnum);
     194                TransferSealevel(femmodel,SealevelNGiaRateEnum);
     195                TransferSealevel(femmodel,SealevelUEsaRateEnum);
     196                TransferSealevel(femmodel,SealevelUGiaRateEnum);
    162197
    163198                //reset cumdeltathickness  to 0:
    … …  
    168203        femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
    169204
    170         /*Free ressources:*/   
    171         delete Sg_eustatic;
    172         delete Sg;
    173         delete Sg_rate;
    174         delete U_radial;
    175         delete U_radial_rate;
     205        /*free ressources:{{{*/
     206        delete RSLg;
     207        delete RSLg_rate;
     208        delete RSLg_eustatic;
     209        delete RSLg_eustatic_rate;
     210        delete U_esa;
     211        delete U_esa_rate;
     212        delete N_esa;
     213        delete N_esa_rate;
     214
    176215        if(horiz){
    177                 delete U_north;
    178                 delete U_east;
    179         }
     216                delete U_north_esa;
     217                delete U_east_esa;
     218        }
     219        delete N_gia;
     220        delete U_gia;
     221        delete N_gia_rate;
     222        delete U_gia_rate;
     223        /*}}}*/
     224
    180225}
    181226/*}}}*/
    … …  
    184229        /*variables:*/
    185230        Vector<IssmDouble> *bedrock  = NULL;
    186         Vector<IssmDouble> *Sg_absolute  = NULL;
     231        Vector<IssmDouble> *SL  = NULL;
    187232        Vector<IssmDouble> *steric_rate_g  = NULL;
    188         Vector<IssmDouble> *rsl_g    = NULL;
    189         Vector<IssmDouble> *bedrock_rate_g  = NULL;
     233        Vector<IssmDouble> *U_esa_rate= NULL;
     234        Vector<IssmDouble> *N_esa_rate= NULL;
     235        Vector<IssmDouble> *U_gia_rate= NULL;
     236        Vector<IssmDouble> *N_gia_rate= NULL;
    190237               
    191238        /*parameters: */
    … …  
    210257        if(VerboseSolution()) _printf0_("         computing steric sea level rise\n");
    211258
    212         /*Retrieve absolute sea level, RSL rate, bedrock uplift rate and steric rate, as vectors:*/
     259        /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
    213260        GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
    214         GetVectorFromInputsx(&Sg_absolute,femmodel,SealevelEnum,VertexSIdEnum);
     261        GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
    215262        GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum);
    216263        if(geodetic){
    217                 GetVectorFromInputsx(&rsl_g,femmodel,SealevelRateEnum,VertexSIdEnum);
    218                 GetVectorFromInputsx(&bedrock_rate_g,femmodel,SealevelUmotionRateEnum,VertexSIdEnum);
    219         }
    220 
    221         /*compute: absolute sea level change = initial absolute sea level + relative sea level change rate * dt + vertical motion rate * dt + steric_rate * dt*/
    222         if(rsl_g) Sg_absolute->AXPY(rsl_g,dt);
    223         if(bedrock_rate_g)Sg_absolute->AXPY(bedrock_rate_g,dt);
    224         Sg_absolute->AXPY(steric_rate_g,dt);
     264                GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
     265                GetVectorFromInputsx(&U_gia_rate,femmodel,SealevelUGiaRateEnum,VertexSIdEnum);
     266                GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
     267                GetVectorFromInputsx(&N_gia_rate,femmodel,SealevelNGiaRateEnum,VertexSIdEnum);
     268        }
     269
     270        /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate - U_gia_rate - U_esa_rate )  * dt + steric_rate * dt*/
     271        if(geodetic){
     272                SL->AXPY(N_gia_rate,dt);
     273                SL->AXPY(N_esa_rate,dt);
     274                SL->AXPY(U_gia_rate,-dt); //notice -dt!
     275                SL->AXPY(U_esa_rate,-dt); //same here
     276        }
     277        SL->AXPY(steric_rate_g,dt);
    225278
    226279        /*compute new bedrock position: */
    227         if(bedrock_rate_g)bedrock->AXPY(bedrock_rate_g,dt);
     280        if(geodetic){
     281                bedrock->AXPY(U_esa_rate,dt);
     282                bedrock->AXPY(U_gia_rate,dt);
     283        }
    228284
    229285        /*update element inputs:*/
    230286        InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
    231         InputUpdateFromVectorx(femmodel,bedrock,SealevelUmotionEnum,VertexSIdEnum);     
    232         InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelEnum,VertexSIdEnum);       
     287        InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);
    233288
    234289        /*Free ressources:*/   
    235290        delete bedrock;
    236         delete Sg_absolute;
     291        delete SL;
    237292        delete steric_rate_g;
    238         if(rsl_g)delete rsl_g;
    239         if(bedrock_rate_g)delete bedrock_rate_g;
    240 
     293        if(geodetic){
     294                delete U_esa_rate;
     295                delete U_gia_rate;
     296                delete N_esa_rate;
     297                delete N_gia_rate;
     298        }
    241299}
    242300/*}}}*/
    … …  
    245303        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    246304
    247         Vector<IssmDouble> *Sgi    = NULL;
    248         IssmDouble          Sgi_oceanaverage   = 0;
     305        Vector<IssmDouble> *RSLgi    = NULL;
     306        IssmDouble          RSLgi_oceanaverage   = 0;
    249307
    250308        /*parameters: */
    … …  
    271329       
    272330        /*Initialize:*/
    273         Sgi = new Vector<IssmDouble>(gsize);
     331        RSLgi = new Vector<IssmDouble>(gsize);
    274332
    275333        /*call the eustatic main module: */
    276         femmodel->SealevelriseEustatic(Sgi,&eustatic, latitude, longitude, radius,loop); //this computes
    277 
    278         /*we need to average Sgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    279         Sgi_oceanaverage=femmodel->SealevelriseOceanAverage(Sgi);
    280 
    281         /*Sg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
     334        femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes
     335
     336        /*we need to average RSLgi over the ocean: RHS term  4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
     337        RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi);
     338
     339        /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
    282340         * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
    283         Sgi->Shift(-eustatic-Sgi_oceanaverage);
    284 
    285         /*save eustatic value for results: */
    286         femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelEustaticEnum,-eustatic));
     341        RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
    287342
    288343        /*clean up and return:*/
    … …  
    290345        xDelete<IssmDouble>(longitude);
    291346        xDelete<IssmDouble>(radius);
    292         return Sgi;
     347        return RSLgi;
    293348}/*}}}*/
    294 Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic){ /*{{{*/
     349Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* RSLg_eustatic){ /*{{{*/
    295350
    296351        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    … …  
    298353
    299354
    300         Vector<IssmDouble> *Sg    = NULL;
    301         Vector<IssmDouble> *Sg_old    = NULL;
    302 
    303         Vector<IssmDouble> *Sgo    = NULL; //ocean convolution of the perturbation to gravity potential.
    304         Vector<IssmDouble> *Sgo_rot= NULL; // rotational feedback
    305         IssmDouble          Sgo_oceanaverage = 0;  //average of Sgo over the ocean.
     355        Vector<IssmDouble> *RSLg    = NULL;
     356        Vector<IssmDouble> *RSLg_old    = NULL;
     357
     358        Vector<IssmDouble> *RSLgo    = NULL; //ocean convolution of the perturbation to gravity potential.
     359        Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback
     360        IssmDouble          RSLgo_oceanaverage = 0;  //average of RSLgo over the ocean.
    306361
    307362        /*parameters: */
    … …  
    340395       
    341396        /*Initialize:*/
    342         Sg = new Vector<IssmDouble>(gsize);
    343         Sg->Assemble();
    344         Sg_eustatic->Copy(Sg);  //first initialize Sg with the eustatic component computed in sealevelrise_core_eustatic.
    345 
    346         Sg_old = new Vector<IssmDouble>(gsize);
    347         Sg_old->Assemble();
     397        RSLg = new Vector<IssmDouble>(gsize);
     398        RSLg->Assemble();
     399        RSLg_eustatic->Copy(RSLg);  //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.
     400
     401        RSLg_old = new Vector<IssmDouble>(gsize);
     402        RSLg_old->Assemble();
    348403
    349404        count=1;
    … …  
    354409
    355410                //save pointer to old sea level rise
    356                 delete Sg_old; Sg_old=Sg;
     411                delete RSLg_old; RSLg_old=RSLg;
    357412
    358413                /*Initialize solution vector: */
    359                 Sg  = new Vector<IssmDouble>(gsize); Sg->Assemble();
    360                 Sgo = new Vector<IssmDouble>(gsize); Sgo->Assemble();
     414                RSLg  = new Vector<IssmDouble>(gsize); RSLg->Assemble();
     415                RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
    361416
    362417                /*call the non eustatic module: */
    363                 femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution,loop);
     418                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
    364419       
    365420                /*assemble solution vector: */
    366                 Sgo->Assemble();
     421                RSLgo->Assemble();
    367422
    368423                if(rotation){
    369424                        /*call rotational feedback  module: */
    370                         Sgo_rot = new Vector<IssmDouble>(gsize); Sgo_rot->Assemble();
    371                         femmodel->SealevelriseRotationalFeedback(Sgo_rot,Sg_old,latitude,longitude,radius);
    372                         Sgo_rot->Assemble();
    373 
    374                         Sgo->AXPY(Sgo_rot,1);
    375                 }
    376                
    377                 /*we need to average Sgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    378                 Sgo_oceanaverage=femmodel->SealevelriseOceanAverage(Sgo);
    379        
    380                 /*Sg is the sum of the eustatic term, and the ocean terms: */
    381                 Sg_eustatic->Copy(Sg); Sg->AXPY(Sgo,1);
    382                 Sg->Shift(-Sgo_oceanaverage);
     425                        RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
     426                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,latitude,longitude,radius);
     427                        RSLgo_rot->Assemble();
     428
     429                        RSLgo->AXPY(RSLgo_rot,1);
     430                }
     431               
     432                /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
     433                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
     434       
     435                /*RSLg is the sum of the eustatic term, and the ocean terms: */
     436                RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
     437                RSLg->Shift(-RSLgo_oceanaverage);
    383438
    384439                /*convergence criterion:*/
    385                 slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs);
     440                slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
    386441
    387442                /*free ressources: */
    388                 delete Sgo;
     443                delete RSLgo;
    389444
    390445                /*Increase count: */
    … …  
    408463        xDelete<IssmDouble>(longitude);
    409464        xDelete<IssmDouble>(radius);
    410         delete Sg_old;
    411 
    412         return Sg;
     465        delete RSLg_old;
     466
     467        return RSLg;
    413468} /*}}}*/
    414 void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* Sg){ /*{{{*/
    415 
    416         Vector<IssmDouble> *U_radial  = NULL;
    417         Vector<IssmDouble> *U_north   = NULL;
    418         Vector<IssmDouble> *U_east    = NULL;
     469void sealevelrise_core_elastic(Vector<IssmDouble>** pU_esa, Vector<IssmDouble>** pU_north_esa,Vector<IssmDouble>** pU_east_esa,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
     470
     471        Vector<IssmDouble> *U_esa  = NULL;
     472        Vector<IssmDouble> *U_north_esa   = NULL;
     473        Vector<IssmDouble> *U_east_esa    = NULL;
    419474               
    420475        /*parameters: */
    … …  
    441496
    442497        /*intialize vectors:*/
    443         U_radial = new Vector<IssmDouble>(gsize);
     498        U_esa = new Vector<IssmDouble>(gsize);
    444499        if (horiz){
    445                 U_north = new Vector<IssmDouble>(gsize);
    446                 U_east = new Vector<IssmDouble>(gsize);
     500                U_north_esa = new Vector<IssmDouble>(gsize);
     501                U_east_esa = new Vector<IssmDouble>(gsize);
    447502        }
    448503
    … …  
    452507       
    453508        /*call the elastic main modlule:*/
    454         femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
     509        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
    455510
    456511        /*Assign output pointers:*/
    457         *pU_radial=U_radial;
     512        *pU_esa=U_esa;
    458513        if(horiz){
    459                 *pU_east=U_east;
    460                 *pU_north=U_north;
     514                *pU_east_esa=U_east_esa;
     515                *pU_north_esa=U_north_esa;
    461516        }
    462517
    … …  
    468523        delete yy;
    469524        delete zz;
     525}
     526/*}}}*/
     527void sealevelrise_core_viscous(Vector<IssmDouble>** pU_gia, Vector<IssmDouble>** pN_gia,FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
     528
     529        Vector<IssmDouble> *U_gia  = NULL;
     530        Vector<IssmDouble> *N_gia  = NULL;
     531               
     532        int configuration_type;
     533        int  gsize;
     534
     535        /*retrieve some parameters:*/
     536        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     537
     538        /*find size of vectors:*/
     539        gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
     540
     541        /*intialize vectors:*/
     542        U_gia = new Vector<IssmDouble>(gsize); U_gia->Set(0);
     543        N_gia = new Vector<IssmDouble>(gsize); N_gia->Set(0);
     544
     545        /*Assign output pointers:*/
     546        *pU_gia=U_gia;
     547        *pN_gia=N_gia;
     548
     549}
     550/*}}}*/
     551void sealevelrise_diagnostics(FemModel* femmodel,Vector<IssmDouble>* RSLg){ /*{{{*/
     552
     553        /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/
     554
    470555}
    471556/*}}}*/
    … …  
    748833
    749834} /*}}}*/
    750 void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
     835void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    751836       
    752837        bool converged=true;
    753838        IssmDouble ndS,nS;
    754         Vector<IssmDouble> *dSg    = NULL;
     839        Vector<IssmDouble> *dRSLg    = NULL;
    755840
    756841        //compute norm(du) and norm(u) if requested
    757         dSg=Sg_old->Duplicate(); Sg_old->Copy(dSg); dSg->AYPX(Sg,-1.0);
    758         ndS=dSg->Norm(NORM_TWO);
     842        dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
     843        ndS=dRSLg->Norm(NORM_TWO);
    759844       
    760845        if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
    761846       
    762847        if(!xIsNan<IssmDouble>(eps_rel)){
    763                 nS=Sg_old->Norm(NORM_TWO);
     848                nS=RSLg_old->Norm(NORM_TWO);
    764849                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
    765850        }
    … …  
    767852
    768853        //clean up
    769         delete dSg;
     854        delete dRSLg;
    770855
    771856        //print

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

    r22159 r22160  
    783783        /*Sea Level Rise{{{*/
    784784        SealevelEnum,
    785         SealevelRateEnum,
    786         SealevelUmotionEnum,
    787         SealevelUmotionRateEnum,
    788         SealevelNmotionEnum,
    789         SealevelEmotionEnum,
    790         SealevelAbsoluteEnum,
    791         SealevelEustaticEnum,
     785        SealevelUEsaEnum,
     786        SealevelUEsaRateEnum,
     787        SealevelNEsaEnum,
     788        SealevelNEsaRateEnum,
     789        SealevelUNorthEsaEnum,
     790        SealevelUEastEsaEnum,
     791        SealevelUGiaEnum,
     792        SealevelUGiaRateEnum,
     793        SealevelNGiaEnum,
     794        SealevelNGiaRateEnum,
     795        SealevelRSLEustaticEnum,
     796        SealevelRSLEustaticRateEnum,
     797        SealevelRSLEnum,
     798        SealevelRSLRateEnum,
    792799        SealevelObsEnum,
    793800        SealevelWeightsEnum,

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

    r22159 r22160  
    763763                case LevelsetReinitFrequencyEnum : return "LevelsetReinitFrequency";
    764764                case SealevelEnum : return "Sealevel";
    765                 case SealevelRateEnum : return "SealevelRate";
    766                 case SealevelUmotionEnum : return "SealevelUmotion";
    767                 case SealevelUmotionRateEnum : return "SealevelUmotionRate";
    768                 case SealevelNmotionEnum : return "SealevelNmotion";
    769                 case SealevelEmotionEnum : return "SealevelEmotion";
    770                 case SealevelAbsoluteEnum : return "SealevelAbsolute";
    771                 case SealevelEustaticEnum : return "SealevelEustatic";
     765                case SealevelUEsaEnum : return "SealevelUEsa";
     766                case SealevelUEsaRateEnum : return "SealevelUEsaRate";
     767                case SealevelNEsaEnum : return "SealevelNEsa";
     768                case SealevelNEsaRateEnum : return "SealevelNEsaRate";
     769                case SealevelUNorthEsaEnum : return "SealevelUNorthEsa";
     770                case SealevelUEastEsaEnum : return "SealevelUEastEsa";
     771                case SealevelUGiaEnum : return "SealevelUGia";
     772                case SealevelUGiaRateEnum : return "SealevelUGiaRate";
     773                case SealevelNGiaEnum : return "SealevelNGia";
     774                case SealevelNGiaRateEnum : return "SealevelNGiaRate";
     775                case SealevelRSLEustaticEnum : return "SealevelRSLEustatic";
     776                case SealevelRSLEustaticRateEnum : return "SealevelRSLEustaticRate";
     777                case SealevelRSLEnum : return "SealevelRSL";
     778                case SealevelRSLRateEnum : return "SealevelRSLRate";
    772779                case SealevelObsEnum : return "SealevelObs";
    773780                case SealevelWeightsEnum : return "SealevelWeights";

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

    r22159 r22160  
    781781              else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum;
    782782              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    783               else if (strcmp(name,"SealevelRate")==0) return SealevelRateEnum;
    784               else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
    785               else if (strcmp(name,"SealevelUmotionRate")==0) return SealevelUmotionRateEnum;
    786               else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
    787               else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    788               else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum;
    789               else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum;
     783              else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum;
     784              else if (strcmp(name,"SealevelUEsaRate")==0) return SealevelUEsaRateEnum;
     785              else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
     786              else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum;
     787              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
     788              else if (strcmp(name,"SealevelUEastEsa")==0) return SealevelUEastEsaEnum;
     789              else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum;
     790              else if (strcmp(name,"SealevelUGiaRate")==0) return SealevelUGiaRateEnum;
     791              else if (strcmp(name,"SealevelNGia")==0) return SealevelNGiaEnum;
     792              else if (strcmp(name,"SealevelNGiaRate")==0) return SealevelNGiaRateEnum;
     793              else if (strcmp(name,"SealevelRSLEustatic")==0) return SealevelRSLEustaticEnum;
     794              else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
     795              else if (strcmp(name,"SealevelRSL")==0) return SealevelRSLEnum;
     796              else if (strcmp(name,"SealevelRSLRate")==0) return SealevelRSLRateEnum;
    790797              else if (strcmp(name,"SealevelObs")==0) return SealevelObsEnum;
    791798              else if (strcmp(name,"SealevelWeights")==0) return SealevelWeightsEnum;
    … …  
    868875              else if (strcmp(name,"Node")==0) return NodeEnum;
    869876              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
    870               else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
    871881              else if (strcmp(name,"Param")==0) return ParamEnum;
    872882              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    … …  
    875885              else if (strcmp(name,"Profiler")==0) return ProfilerEnum;
    876886              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"Masscon")==0) return MassconEnum;
     887              else if (strcmp(name,"Masscon")==0) return MassconEnum;
    881888              else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
    882889              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    … …  
    991998              else if (strcmp(name,"IceMass")==0) return IceMassEnum;
    992999              else if (strcmp(name,"IceVolume")==0) return IceVolumeEnum;
    993               else if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"IceVolumeAboveFloatation")==0) return IceVolumeAboveFloatationEnum;
    9941004              else if (strcmp(name,"TotalFloatingBmb")==0) return TotalFloatingBmbEnum;
    9951005              else if (strcmp(name,"TotalGroundedBmb")==0) return TotalGroundedBmbEnum;
    … …  
    9981008              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
    9991009              else if (strcmp(name,"P1")==0) return P1Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"P1DG")==0) return P1DGEnum;
     1010              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    10041011              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    10051012              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
Note: See TracChangeset for help on using the changeset viewer.