Changeset 22155


Ignore:
Timestamp:
10/09/17 15:30:06 (7 years ago)
Author:
Eric.Larour
Message:

CHG: rearranging the entire slr core

Location:
issm/branches/trunk-larour-NatGeoScience2016/src
Files:
2 deleted
14 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/Makefile.am

    r21759 r22155  
    476476if SEALEVELRISE
    477477issm_sources +=  ./cores/sealevelrise_core.cpp\
    478                                  ./cores/sealevelrise_core_eustatic.cpp\
    479                                  ./cores/sealevelrise_core_noneustatic.cpp\
    480478                                 ./analyses/SealevelriseAnalysis.cpp
    481479endif

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

    r22143 r22155  
    2020void SealevelriseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    2121
     22        int geodetic=0;
    2223
    2324        /*Update elements: */
    … …  
    3435        iomodel->FetchDataToInput(elements,"md.mask.groundedice_levelset",MaskGroundediceLevelsetEnum);
    3536        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    36         iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
    37         iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
     37        //those only if we have requested geodetic computations:
     38        iomodel->FetchData(&geodetic,"md.slr.geodetic");
     39        if (geodetic){
     40                iomodel->FetchDataToInput(elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum);
     41                iomodel->FetchDataToInput(elements,"md.mask.land_levelset",MaskLandLevelsetEnum);
     42        }
    3843        iomodel->FetchDataToInput(elements,"md.slr.deltathickness",SealevelriseDeltathicknessEnum);
    3944        iomodel->FetchDataToInput(elements,"md.slr.spcthickness",SealevelriseSpcthicknessEnum);
    4045        iomodel->FetchDataToInput(elements,"md.slr.sealevel",SealevelEnum,0);
     46        iomodel->FetchDataToInput(elements,"md.geometry.bed",BedEnum);
    4147
    4248        /*Initialize cumdeltalthickness input: unfortunately, we don't have femmodel, so we need to iterate on
    … …  
    96102        parameters->AddObject(iomodel->CopyConstantObject("md.slr.angular_velocity",SealevelriseAngularVelocityEnum));
    97103        parameters->AddObject(iomodel->CopyConstantObject("md.slr.ocean_area_scaling",SealevelriseOceanAreaScalingEnum));
     104        parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic",SealevelriseGeodeticEnum));
    98105
    99106        iomodel->FetchData(&elastic,"md.slr.elastic");

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

    r21759 r22155  
    16421642                                name==SealevelEnum ||
    16431643                                name==SealevelUmotionEnum ||
     1644                                name==SealevelUmotionRateEnum ||
    16441645                                name==SealevelNmotionEnum ||
    16451646                                name==SealevelEmotionEnum ||
    … …  
    16481649                                name==SealevelriseDeltathicknessEnum ||
    16491650                                name==SealevelriseCumDeltathicknessEnum ||
     1651                                name==SealevelRateEnum ||
    16501652                                name==EsaUmotionEnum ||
    16511653                                name==EsaNmotionEnum ||

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

    r22143 r22155  
    5151void damage_core(FemModel* femmodel);
    5252void sealevelrise_core(FemModel* femmodel);
     53void geodetic_core(FemModel* femmodel);
     54void steric_core(FemModel* femmodel);
    5355Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
    5456Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
     57void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* Sg);
    5558IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
    5659
    … …  
    6669void TransferSealevel(FemModel* femmodel,int forcingenum);
    6770void EarthMassTransport(FemModel* femmodel);
     71void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs);
    6872
    6973//solution configuration

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

    r22144 r22155  
    1010#include "../solutionsequences/solutionsequences.h"
    1111
     12/*cores:*/
    1213void sealevelrise_core(FemModel* femmodel){ /*{{{*/
    1314
     15
     16        /*Parameters, variables:*/
     17        bool save_results;
     18        bool isslr=0;
     19        int solution_type;
     20               
     21        /*Retrieve parameters:*/
     22        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
     23        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     24        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     25       
     26        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
     27        if(solution_type==SealevelriseSolutionEnum)isslr=1;
     28
     29        /*Should we be here?:*/
     30        if(!isslr)return;
     31
     32        /*Verbose: */
     33        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
     34
     35        /*set SLR configuration: */
     36        femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
     37
     38        /*Run geodetic:*/
     39        geodetic_core(femmodel);
     40
     41        /*Run steric core for sure:*/
     42        steric_core(femmodel);
     43       
     44        /*Save results: */
     45        if(save_results){
     46                int        numoutputs        = 0;
     47                char     **requested_outputs = NULL;
     48
     49                if(VerboseSolution()) _printf0_("   saving results\n");
     50                femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
     51                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
     52                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
     53        }
     54
     55        /*requested dependents: */
     56        if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
     57}
     58/*}}}*/
     59void geodetic_core(FemModel* femmodel){ /*{{{*/
     60
     61
     62        /*variables:*/
    1463        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;
     71
     72        /*parameters:*/
     73        bool iscoupler;
     74        int  configuration_type;
     75        int  modelid,earthid;
     76        bool istransientmasstransport;
     77        int  frequency,count;
     78        int  horiz;
     79        int  geodetic=0;
     80        IssmDouble          dt;
     81
     82        /*Should we even be here?:*/
     83        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
     84       
     85        /*Verbose: */
     86        if(VerboseSolution()) _printf0_("         computing geodetic sea level rise\n");
     87
     88        /*retrieve more parameters:*/
     89        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
     90        femmodel->parameters->FindParam(&frequency,SealevelriseGeodeticRunFrequencyEnum);
     91        femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
     92        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     93        femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
     94        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     95
     96        if(iscoupler){
     97                femmodel->parameters->FindParam(&modelid,ModelIdEnum);
     98                femmodel->parameters->FindParam(&earthid,EarthIdEnum);
     99                femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
     100        }
     101        else{
     102                /* we are here, we are not running in a coupler, so we will indeed compute SLR,
     103                 * so make sure we are identified as being the Earth.:*/
     104                modelid=1; earthid=1;
     105                /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we
     106                 * run, irresepective of the time settings:*/
     107                count=frequency;
     108        }
     109
     110        /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if
     111         * not already done by the mass trasnport module. For ice caps, they rely on the transient mass
     112         * transport module exclusively:*/
     113        if(iscoupler) if(modelid==earthid) if(istransientmasstransport) EarthMassTransport(femmodel);
     114
     115        /*increment counter, or call solution core if count==frequency:*/
     116        if (count<frequency){
     117                count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum);
     118                return;
     119        }
     120
     121        /*call sea-level rise sub cores:*/
     122        if(iscoupler){
     123                /*transfer cumulated deltathickness forcing from ice caps to earth model: */
     124                TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
     125
     126                /*we have accumulated thicknesses, dump them in deltathcikness: */
     127                if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SealevelriseDeltathicknessEnum);
     128        }
     129
     130        /*run cores:*/
     131        if(modelid==earthid){
     132
     133                /*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);
     135
     136                /*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               
     139
     140                /*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);
     142               
     143                /*transform these values into rates (as we only run this once each frequency turns:*/
     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                }
     155        }
     156
     157
     158        if(iscoupler){
     159                /*transfer sea level back to ice caps:*/
     160                TransferSealevel(femmodel,SealevelRateEnum);
     161                TransferSealevel(femmodel,SealevelUmotionRateEnum);
     162
     163                //reset cumdeltathickness  to 0:
     164                InputUpdateFromConstantx(femmodel,0,SealevelriseCumDeltathicknessEnum);
     165        }
     166
     167        /*reset counter to 1:*/
     168        femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
     169
     170        /*Free ressources:*/   
     171        delete Sg_eustatic;
     172        delete Sg;
     173        delete Sg_rate;
     174        delete U_radial;
     175        delete U_radial_rate;
     176        if(horiz){
     177                delete U_north;
     178                delete U_east;
     179        }
     180}
     181/*}}}*/
     182void steric_core(FemModel* femmodel){ /*{{{*/
     183
     184        /*variables:*/
     185        Vector<IssmDouble> *bedrock  = NULL;
    15186        Vector<IssmDouble> *Sg_absolute  = NULL;
    16         Vector<IssmDouble> *Sg_eustatic  = NULL;
    17187        Vector<IssmDouble> *steric_rate_g  = NULL;
     188        Vector<IssmDouble> *rsl_g    = NULL;
     189        Vector<IssmDouble> *bedrock_rate_g  = NULL;
     190               
     191        /*parameters: */
     192        bool isslr=0;
     193        int  solution_type;
     194        IssmDouble          dt;
     195        int  geodetic=0;
     196
     197        /*Retrieve parameters:*/
     198        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
     199        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     200        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
     201        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum);
     202
     203        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
     204        if(solution_type==SealevelriseSolutionEnum)isslr=1;
     205       
     206        /*Should we be here?:*/
     207        if(!isslr)return;
     208
     209        /*Verbose: */
     210        if(VerboseSolution()) _printf0_("         computing steric sea level rise\n");
     211
     212        /*Retrieve absolute sea level, RSL rate, bedrock uplift rate and steric rate, as vectors:*/
     213        GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
     214        GetVectorFromInputsx(&Sg_absolute,femmodel,SealevelEnum,VertexSIdEnum);
     215        GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum);
     216        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);
     225
     226        /*compute new bedrock position: */
     227        if(bedrock_rate_g)bedrock->AXPY(bedrock_rate_g,dt);
     228
     229        /*update element inputs:*/
     230        InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
     231        InputUpdateFromVectorx(femmodel,bedrock,SealevelUmotionEnum,VertexSIdEnum);     
     232        InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelEnum,VertexSIdEnum);       
     233
     234        /*Free ressources:*/   
     235        delete bedrock;
     236        delete Sg_absolute;
     237        delete steric_rate_g;
     238        if(rsl_g)delete rsl_g;
     239        if(bedrock_rate_g)delete bedrock_rate_g;
     240
     241}
     242/*}}}*/
     243Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
     244 
     245        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
     246
     247        Vector<IssmDouble> *Sgi    = NULL;
     248        IssmDouble          Sgi_oceanaverage   = 0;
     249
     250        /*parameters: */
     251        int  configuration_type;
     252        int  gsize;
     253        bool spherical=true;
     254        IssmDouble *latitude  = NULL;
     255        IssmDouble *longitude = NULL;
     256        IssmDouble *radius    = NULL;
     257        int         loop;
     258
     259        /*outputs:*/
     260        IssmDouble eustatic;
     261
     262        /*recover parameters:*/
     263        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     264        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
     265
     266        /*first, recover lat,long and radius vectors from vertices: */
     267        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
     268
     269        /*Figure out size of g-set deflection vector and allocate solution vector: */
     270        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
     271       
     272        /*Initialize:*/
     273        Sgi = new Vector<IssmDouble>(gsize);
     274
     275        /*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
     282         * 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));
     287
     288        /*clean up and return:*/
     289        xDelete<IssmDouble>(latitude);
     290        xDelete<IssmDouble>(longitude);
     291        xDelete<IssmDouble>(radius);
     292        return Sgi;
     293}/*}}}*/
     294Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic){ /*{{{*/
     295
     296        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
     297          non eustatic core of the SLR solution */
     298
     299
     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.
     306
     307        /*parameters: */
     308        int count;
     309        bool save_results;
     310        int  gsize;
     311        int  configuration_type;
     312        bool spherical=true;
     313        bool converged=true;
     314        bool rotation=true;
     315        bool verboseconvolution=true;
     316        int max_nonlinear_iterations;
     317        IssmDouble           eps_rel;
     318        IssmDouble           eps_abs;
     319        IssmDouble          *latitude    = NULL;
     320        IssmDouble          *longitude    = NULL;
     321        IssmDouble          *radius    = NULL;
     322        IssmDouble           eustatic;
     323        int                      loop;
     324
     325        /*Recover some parameters: */
     326        femmodel->parameters->FindParam(&max_nonlinear_iterations,SealevelriseMaxiterEnum);
     327        femmodel->parameters->FindParam(&eps_rel,SealevelriseReltolEnum);
     328        femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
     329        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     330        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
     331       
     332        /*computational flag: */
     333        femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
     334
     335        /*first, recover lat,long and radius vectors from vertices: */
     336        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
     337
     338        /*Figure out size of g-set deflection vector and allocate solution vector: */
     339        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
     340       
     341        /*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();
     348
     349        count=1;
     350        converged=false;
     351       
     352        /*Start loop: */
     353        for(;;){
     354
     355                //save pointer to old sea level rise
     356                delete Sg_old; Sg_old=Sg;
     357
     358                /*Initialize solution vector: */
     359                Sg  = new Vector<IssmDouble>(gsize); Sg->Assemble();
     360                Sgo = new Vector<IssmDouble>(gsize); Sgo->Assemble();
     361
     362                /*call the non eustatic module: */
     363                femmodel->SealevelriseNonEustatic(Sgo, Sg_old, latitude, longitude, radius,verboseconvolution,loop);
     364       
     365                /*assemble solution vector: */
     366                Sgo->Assemble();
     367
     368                if(rotation){
     369                        /*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);
     383
     384                /*convergence criterion:*/
     385                slrconvergence(&converged,Sg,Sg_old,eps_rel,eps_abs);
     386
     387                /*free ressources: */
     388                delete Sgo;
     389
     390                /*Increase count: */
     391                count++;
     392                if(converged==true){
     393                        break;
     394                }
     395                if(count>=max_nonlinear_iterations){
     396                        _printf0_("   maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
     397                        converged=true;
     398                        break;
     399                }       
     400               
     401                /*some minor verbosing adjustment:*/
     402                if(count>1)verboseconvolution=false;
     403               
     404        }
     405        if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count-1 << "\n");
     406
     407        xDelete<IssmDouble>(latitude);
     408        xDelete<IssmDouble>(longitude);
     409        xDelete<IssmDouble>(radius);
     410        delete Sg_old;
     411
     412        return Sg;
     413} /*}}}*/
     414void sealevelrise_core_elastic(Vector<IssmDouble>** pU_radial, Vector<IssmDouble>** pU_north,Vector<IssmDouble>** pU_east,FemModel* femmodel,Vector<IssmDouble>* Sg){ /*{{{*/
     415
    18416        Vector<IssmDouble> *U_radial  = NULL;
    19417        Vector<IssmDouble> *U_north   = NULL;
    20418        Vector<IssmDouble> *U_east    = NULL;
    21         Vector<IssmDouble>* cumdeltathickness=NULL;
    22         bool save_results,isslr,iscoupler;
     419               
     420        /*parameters: */
    23421        int configuration_type;
    24         int solution_type;
    25         int        numoutputs        = 0;
    26         char     **requested_outputs = NULL;
    27        
    28         /*additional parameters: */
    29422        int  gsize;
    30423        bool spherical=true;
     424
    31425        IssmDouble          *latitude   = NULL;
    32426        IssmDouble          *longitude  = NULL;
    … …  
    35429        IssmDouble          *yy     = NULL;
    36430        IssmDouble          *zz     = NULL;
    37         IssmDouble          dt,steric_rate;
    38         int                 frequency,count;
    39431        int  loop;
    40432        int  horiz;
    41433
    42         /*Recover some parameters: */
     434        /*retrieve some parameters:*/
    43435        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    44         femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    45         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    46         femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    47         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
    48         femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum);
    49         femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
    50 
    51        
    52         /*first, recover lat,long and radius vectors from vertices: */
     436        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
     437        femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
     438
     439        /*find size of vectors:*/
     440        gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
     441
     442        /*intialize vectors:*/
     443        U_radial = new Vector<IssmDouble>(gsize);
     444        if (horiz){
     445                U_north = new Vector<IssmDouble>(gsize);
     446                U_east = new Vector<IssmDouble>(gsize);
     447        }
     448
     449        /*retrieve geometric information: */
    53450        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    54 
    55         /*recover x,y,z vectors from vertices: */
    56451        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    57 
    58         /*Figure out size of g-set deflection vector and allocate solution vector: */
    59         gsize      = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    60 
    61         /*several cases here, depending on value of iscoupler and isslr:
    62         solution_type == SealevelriseSolutionEnum)       we are running sea level rise core (no coupler)
    63         ( !iscoupler & !isslr)       we are not interested in being here :)
    64         ( !iscoupler & isslr)        we are running in uncoupled mode
    65         ( iscoupler & isslr)         we are running in coupled mode, and better be earth
    66         ( iscoupler & !isslr)        we are running in coupled mode, and better be an ice cap
    67         */
    68 
    69         if(solution_type==SealevelriseSolutionEnum){
    70                 isslr=1;
    71                 iscoupler=0;
    72                 count=frequency;
    73         }
    74 
    75         /*early return: */
    76         if( !iscoupler & !isslr) return;  //we are not interested in being here :)
    77 
    78         /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/
    79         if(VerboseSolution()) _printf0_("   computing sea level rise\n");
    80 
    81         /*set configuration: */
    82         if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
    83 
    84         /*figure out whether deltathickness was computed on the earth (mass transport module should
    85          * have taken care of it for ice caps:  */
    86         if(iscoupler){
    87                 int modelid,earthid;
    88                 bool ismasstransport;
    89                 femmodel->parameters->FindParam(&modelid,ModelIdEnum);
    90                 femmodel->parameters->FindParam(&earthid,EarthIdEnum);
    91                 femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
    92                 if((modelid==earthid) && (ismasstransport==0)) EarthMassTransport(femmodel);
    93         }
    94 
    95         /*transfer cum deltathickness forcing from ice caps to earth model: */
    96         if(iscoupler & (count==frequency)) TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
    97 
    98         /*call sea-level rise sub cores:*/
    99         if(isslr & (count==frequency)){
    100                        
    101                 if(iscoupler){
    102                         /*we have accumulated thicknesses, dump them in deltathcikness: */
    103                         GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
    104                         InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
    105                         delete cumdeltathickness;
    106                 }
    107 
    108                 femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
    109                 femmodel->parameters->FindParam(&horiz,SealevelriseHorizEnum);
    110 
    111                 Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
    112 
    113                 Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //ocean loading tems  (2nd and 5th terms on the RHS of Farrel and Clark)
    114 
    115                 /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
    116                 /*Initialize:*/
    117                 U_radial = new Vector<IssmDouble>(gsize);
    118                 if (horiz){
    119                         U_north = new Vector<IssmDouble>(gsize);
    120                         U_east = new Vector<IssmDouble>(gsize);
    121                 }
    122                 Sg_absolute = new Vector<IssmDouble>(gsize);
    123                
    124                 /*call the geodetic main modlule:*/
    125                 femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
    126                                
    127                 /*Now deal with steric ocean expansion by just shifting Sg by a spatial rate pattern : */
    128                 femmodel->parameters->FindParam(&frequency,SealevelriseRunFrequencyEnum);
    129                 femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    130                 GetVectorFromInputsx(&steric_rate_g,femmodel,SealevelriseStericRateEnum,VertexSIdEnum);
    131                 Sg->AXPY(steric_rate_g,dt*frequency);
    132                
    133                 /*compute: absolute sea level change = relative sea level change + vertical motion*/
    134                 Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1);
    135                
    136                 /*get results into elements:*/
    137                 InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum);    // radial displacement
    138                 if (horiz){
    139                         InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion
    140                         InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion
    141                 }
    142                 InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum); //absolute sea level
    143                 InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum); //relative sea level
    144        
    145                 if(save_results){
    146                         if(VerboseSolution()) _printf0_("   saving results\n");
    147                         femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
    148                         femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    149                 }
    150 
    151                 /*requested dependents: */
    152                 if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
    153 
    154                 if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    155                
    156                 /*Free ressources:*/   
    157                 delete U_radial;
    158                 if(horiz){
    159                         delete U_north;
    160                         delete U_east;
    161                 }
    162                 delete Sg_absolute;
    163                 delete Sg_eustatic;
    164                 delete Sg;
    165                 delete steric_rate_g;
    166 
    167         }
    168                
    169         /*transfer sea-level back to ice caps, reset cum thicknesses to 0 and reset counter: */
    170         if(iscoupler){
    171                 if (count==frequency) {
    172                         //transfer sea level back to ice caps:
    173                         TransferSealevel(femmodel,SealevelEnum);
    174 
    175                         //reset cumdeltathickness  to 0:
    176                         GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
    177                         cumdeltathickness->Set(0);
    178                         InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
    179                         delete cumdeltathickness;
    180 
    181                         //reset counter to 1:
    182                         femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
    183                 }
    184                 else{
    185                         //increment counter:
    186                         count++;
    187                         femmodel->parameters->SetParam(count,SealevelriseRunCountEnum);
    188                 }
    189         }
    190        
     452       
     453        /*call the elastic main modlule:*/
     454        femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
     455
     456        /*Assign output pointers:*/
     457        *pU_radial=U_radial;
     458        if(horiz){
     459                *pU_east=U_east;
     460                *pU_north=U_north;
     461        }
     462
    191463        /*Free ressources: */
    192464        delete longitude;
    … …  
    196468        delete yy;
    197469        delete zz;
    198 
    199 }
     470}
    200471/*}}}*/
     472
     473/*support routines:*/
    201474void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
    202475
    … …  
    475748
    476749} /*}}}*/
     750void slrconvergence(bool* pconverged, Vector<IssmDouble>* Sg,Vector<IssmDouble>* Sg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
     751       
     752        bool converged=true;
     753        IssmDouble ndS,nS;
     754        Vector<IssmDouble> *dSg    = NULL;
     755
     756        //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);
     759       
     760        if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
     761       
     762        if(!xIsNan<IssmDouble>(eps_rel)){
     763                nS=Sg_old->Norm(NORM_TWO);
     764                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
     765        }
     766
     767
     768        //clean up
     769        delete dSg;
     770
     771        //print
     772        if(!xIsNan<IssmDouble>(eps_rel)){
     773                if((ndS/nS)<eps_rel){
     774                        if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
     775                }
     776                else{
     777                        if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
     778                        converged=false;
     779                }
     780        }
     781        if(!xIsNan<IssmDouble>(eps_abs)){
     782                if(ndS<eps_abs){
     783                        if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
     784                }
     785                else{
     786                        if(VerboseConvergence()) _printf0_(setw(50) << left << "      convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
     787                        converged=false;
     788                }
     789        }
     790
     791        /*assign output*/
     792        *pconverged=converged;
     793
     794} /*}}}*/

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

    r21759 r22155  
    2121        /*parameters: */
    2222        IssmDouble finaltime,dt,yts;
    23         bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology;
     23        bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,ismovingfront,isdamageevolution,ishydrology;
    2424        bool       save_results,dakota_analysis;
    2525        bool       time_adapt;
    … …  
    5555        femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
    5656        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    57         femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
    5857        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
    5958        femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
    … …  
    154153
    155154                /*Sea level rise: */
    156                 if(isslr || iscoupler) sealevelrise_core(femmodel);
     155                if(isslr) sealevelrise_core(femmodel);
    157156
    158157                /*unload results*/

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r22143 r22155  
    6666        parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
    6767        parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
    68         parameters->AddObject(iomodel->CopyConstantObject("md.slr.run_frequency",SealevelriseRunFrequencyEnum));
     68        parameters->AddObject(iomodel->CopyConstantObject("md.slr.geodetic_run_frequency",SealevelriseGeodeticRunFrequencyEnum));
    6969        parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); 
    7070        {/*This is specific to ice...*/

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

    r22114 r22155  
    782782        /*Sea Level Rise{{{*/
    783783        SealevelEnum,
     784        SealevelRateEnum,
    784785        SealevelUmotionEnum,
     786        SealevelUmotionRateEnum,
    785787        SealevelNmotionEnum,
    786788        SealevelEmotionEnum,
    … …  
    801803        SealevelriseOceanAreaScalingEnum,
    802804        SealevelriseStericRateEnum,
    803         SealevelriseRunFrequencyEnum,
     805        SealevelriseGeodeticRunFrequencyEnum,
     806        SealevelriseGeodeticEnum,
    804807        SealevelriseRunCountEnum,
    805808        SealevelriseRotationEnum,

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

    r22114 r22155  
    762762                case LevelsetReinitFrequencyEnum : return "LevelsetReinitFrequency";
    763763                case SealevelEnum : return "Sealevel";
     764                case SealevelRateEnum : return "SealevelRate";
    764765                case SealevelUmotionEnum : return "SealevelUmotion";
     766                case SealevelUmotionRateEnum : return "SealevelUmotionRate";
    765767                case SealevelNmotionEnum : return "SealevelNmotion";
    766768                case SealevelEmotionEnum : return "SealevelEmotion";
    … …  
    781783                case SealevelriseOceanAreaScalingEnum : return "SealevelriseOceanAreaScaling";
    782784                case SealevelriseStericRateEnum : return "SealevelriseStericRate";
    783                 case SealevelriseRunFrequencyEnum : return "SealevelriseRunFrequency";
     785                case SealevelriseGeodeticRunFrequencyEnum : return "SealevelriseGeodeticRunFrequency";
     786                case SealevelriseGeodeticEnum : return "SealevelriseGeodetic";
    784787                case SealevelriseRunCountEnum : return "SealevelriseRunCount";
    785788                case SealevelriseRotationEnum : return "SealevelriseRotation";

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

    r22114 r22155  
    780780              else if (strcmp(name,"LevelsetReinitFrequency")==0) return LevelsetReinitFrequencyEnum;
    781781              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
     782              else if (strcmp(name,"SealevelRate")==0) return SealevelRateEnum;
    782783              else if (strcmp(name,"SealevelUmotion")==0) return SealevelUmotionEnum;
     784              else if (strcmp(name,"SealevelUmotionRate")==0) return SealevelUmotionRateEnum;
    783785              else if (strcmp(name,"SealevelNmotion")==0) return SealevelNmotionEnum;
    784786              else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum;
    … …  
    799801              else if (strcmp(name,"SealevelriseOceanAreaScaling")==0) return SealevelriseOceanAreaScalingEnum;
    800802              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
    801               else if (strcmp(name,"SealevelriseRunFrequency")==0) return SealevelriseRunFrequencyEnum;
     803              else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum;
     804              else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum;
    802805              else if (strcmp(name,"SealevelriseRunCount")==0) return SealevelriseRunCountEnum;
    803806              else if (strcmp(name,"SealevelriseRotation")==0) return SealevelriseRotationEnum;
    … …  
    872875              else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
    873876              else if (strcmp(name,"Masscon")==0) return MassconEnum;
    874               else if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
    875               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    876               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
     880              if (strcmp(name,"Massconaxpby")==0) return MassconaxpbyEnum;
     881              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
     882              else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
     883              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    881884              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    882885              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    … …  
    995998              else if (strcmp(name,"P1")==0) return P1Enum;
    996999              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    997               else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    998               else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    999               else if (strcmp(name,"P2")==0) return P2Enum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
     1003              if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
     1004              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
     1005              else if (strcmp(name,"P2")==0) return P2Enum;
     1006              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    10041007              else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    10051008              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/pairoptions.m

    r20840 r22155  
    127127                        for i=1:numoptions,
    128128                                if ~self.list{i,3},
    129                                         disp(['WARNING: option ' self.list{i,1} ' was not used'])
     129                                        %disp(['WARNING: option ' self.list{i,1} ' was not used'])
    130130                                end
    131131                        end

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/classes/sealevelmodel.m

    r22143 r22155  
    6161                        %check that run_frequency is the same everywhere:
    6262                        for i=1:length(slm.icecaps),
    63                                 if slm.icecaps{i}.slr.run_frequency~=slm.earth.slr.run_frequency,
     63                                if slm.icecaps{i}.slr.geodetic_run_frequency~=slm.earth.slr.geodetic_run_frequency,
    6464                                        error(sprintf('sealevelmodel checkconsistenty error:  icecap model %s should have the same run frequency as earth!',slm.icecaps{i}.miscellaneous.name));
    6565                                end

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

    r22114 r22155  
    2626                ocean_area_scaling = 0;
    2727                steric_rate    = 0; %rate of ocean expansion from steric effects.
    28                 run_frequency  = 1; %how many time steps we skip before we run SLR solver during transient
     28                geodetic_run_frequency  = 1; %how many time steps we skip before we run the geodetic part of the solver during transient
     29                geodetic       = 0; %compute geodetic SLR? (in addition to steric?)
    2930                degacc         = 0;
    3031                loop_increment = 0;
    … …  
    5354
    5455                %computational flags:
     56                self.geodetic=0;
    5557                self.rigid=1;
    5658                self.elastic=1;
    … …  
    8082       
    8183                %how many time steps we skip before we run SLR solver during transient
    82                 self.run_frequency=1;
     84                self.geodetic_run_frequency=1;
    8385       
    8486                %output default:
    … …  
    110112                        md = checkfield(md,'fieldname','slr.abstol','size',[1 1]);
    111113                        md = checkfield(md,'fieldname','slr.maxiter','size',[1 1],'>=',1);
    112                         md = checkfield(md,'fieldname','slr.run_frequency','size',[1 1],'>=',1);
     114                        md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1 1],'>=',1);
    113115                        md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1]);
    114116                        md = checkfield(md,'fieldname','slr.degacc','size',[1 1],'>=',1e-10);
    … …  
    129131                                error('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!');
    130132                        end
     133
     134                        %check that  if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
     135                        %a coupler to a planet model is provided.
     136                        if self.geodetic,
     137                                if md.transient.iscoupler,
     138                                        %we are good;
     139                                else
     140                                        if strcmpi(class(md.mesh),'mesh3dsurface'),
     141                                                %we are good
     142                                        else
     143                                                error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!');
     144                                        end
     145                                end
     146                        end
     147
    131148
    132149                end % }}}
    … …  
    154171                        fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
    155172                        fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)');
    156                         fielddisplay(self,'run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
     173                        fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)');
    157174                        fielddisplay(self,'rigid','rigid earth graviational potential perturbation');
    158175                        fielddisplay(self,'elastic','elastic earth graviational potential perturbation');
    159176                        fielddisplay(self,'rotation','earth rotational potential perturbation');
    160                         fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]');
    161                         fielddisplay(self,'steric_rate','rate of steric ocean expansion (in mm/yr)');
    162177                        fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions');
    163178                        fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
    … …  
    186201                        WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean');
    187202                        WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean');
    188                         WriteData(fid,prefix,'object',self,'fieldname','run_frequency','format','Integer');
     203                        WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer');
    189204                        WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts);
    190205                        WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double');
    … …  
    192207                        WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer');
    193208                        WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer');
     209                        WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer');
    194210                       
    195211                        %process requested outputs
    … …  
    224240                        writejsdouble(fid,[modelname '.slr.rotation'],self.rotation);
    225241                        writejsdouble(fid,[modelname '.slr.ocean_area_scaling'],self.ocean_area_scaling);
    226                         writejsdouble(fid,[modelname '.slr.run_frequency'],self.run_frequency);
     242                        writejsdouble(fid,[modelname '.slr.geodetic_run_frequency'],self.geodetic_run_frequency);
    227243                        writejsdouble(fid,[modelname '.slr.elastic'],self.elastic);
    228244                        writejs1Darray(fid,[modelname '.slr.steric_rate'],self.steric_rate);

  • issm/branches/trunk-larour-NatGeoScience2016/src/m/solve/solveslm.m

    r22112 r22155  
    3434options=pairoptions(varargin{:},'solutionstring',solutionstring);
    3535
    36 %figure out if the sum of cluster processors requested sums up correctly:
     36%make sure we request sum of cluster processors
    3737totalnp=0;
    3838for i=1:length(slm.icecaps), totalnp=totalnp+slm.icecaps{i}.cluster.np; end
Note: See TracChangeset for help on using the changeset viewer.