Changeset 24013


Ignore:
Timestamp:
06/12/19 14:40:27 (6 years ago)
Author:
jdquinn
Message:

BUG: Attempting to correct Windows issue with non-constants.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r24010 r24013  
    6666
    6767        /*Allocate arrays*/
    68         int numvertices = this->GetNumberOfVertices();
    69         IssmDouble* lambdas = xNew<IssmDouble>(numvertices);
     68        const int NUM_VERTICES = this->GetNumberOfVertices();
     69
     70        IssmDouble* lambdas = xNew<IssmDouble>(NUM_VERTICES);
    7071
    7172        /* Start looping on the number of vertices: */
    7273        Gauss* gauss=this->NewGauss();
    73         for (int iv=0;iv<numvertices;iv++){
     74        for (int iv=0;iv<NUM_VERTICES;iv++){
    7475                gauss->GaussVertex(iv);
    7576
    … …  
    245246
    246247        /*Allocate arrays*/
    247         int numvertices = this->GetNumberOfVertices();
    248         IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
    249         IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
    250         IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
    251         IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
    252         IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
    253         IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);
    254         IssmDouble* eps_ef = xNew<IssmDouble>(numvertices);
     248        const int NUM_VERTICES = this->GetNumberOfVertices();
     249
     250        IssmDouble* eps_xx = xNew<IssmDouble>(NUM_VERTICES);
     251        IssmDouble* eps_yy = xNew<IssmDouble>(NUM_VERTICES);
     252        IssmDouble* eps_zz = xNew<IssmDouble>(NUM_VERTICES);
     253        IssmDouble* eps_xy = xNew<IssmDouble>(NUM_VERTICES);
     254        IssmDouble* eps_xz = xNew<IssmDouble>(NUM_VERTICES);
     255        IssmDouble* eps_yz = xNew<IssmDouble>(NUM_VERTICES);
     256        IssmDouble* eps_ef = xNew<IssmDouble>(NUM_VERTICES);
    255257
    256258        /* Start looping on the number of vertices: */
    257259        Gauss* gauss=this->NewGauss();
    258         for (int iv=0;iv<numvertices;iv++){
     260        for (int iv=0;iv<NUM_VERTICES;iv++){
    259261                gauss->GaussVertex(iv);
    260262
    … …  
    388390        _printf_("   sid: "<<this->sid<<"\n");
    389391        if(vertices){
    390                 int numvertices = this->GetNumberOfVertices();
    391                 for(int i=0;i<numvertices;i++) vertices[i]->Echo();
     392                const int NUM_VERTICES = this->GetNumberOfVertices();
     393                for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
    392394        }
    393395        else _printf_("vertices = NULL\n");
    … …  
    427429        if(!IsOnBase()) return;
    428430
    429         int        numvertices = this->GetNumberOfVertices();
     431        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
     432        const int NUM_VERTICES_MONTHS_PER_YEAR  = NUM_VERTICES * 12;
    430433
    431434        int        i;
    432         IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
    433         IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    434         IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
    435         IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
    436         IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
    437         IssmDouble* tmp=xNew<IssmDouble>(numvertices);
     435        IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     436        IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     437        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     438        IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     439        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     440        IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
    438441        IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
    439442        IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
    … …  
    453456        Gauss* gauss=this->NewGauss();
    454457        for(int month=0;month<12;month++){
    455                 for(int iv=0;iv<numvertices;iv++){
     458                for(int iv=0;iv<NUM_VERTICES;iv++){
    456459                        gauss->GaussVertex(iv);
    457460                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
    … …  
    472475
    473476        /*Compute the temperature and precipitation*/
    474         for(int iv=0;iv<numvertices;iv++){
     477        for(int iv=0;iv<NUM_VERTICES;iv++){
    475478                ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
    476479                                        Delta18oPresent, Delta18oLgm, Delta18oTime,
    … …  
    484487        TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
    485488        for (int imonth=0;imonth<12;imonth++) {
    486                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
     489                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    487490                switch(this->ObjectEnum()){
    488491                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    491494                        default: _error_("Not implemented yet");
    492495                }
    493                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
     496                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    494497                switch(this->ObjectEnum()){
    495498                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    530533        if(!IsOnBase()) return;
    531534
    532         int        numvertices = this->GetNumberOfVertices();
    533 
    534         int        i;
    535         IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
    536         IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    537         IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
    538         IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
    539         IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(12*numvertices);
    540         IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(12*numvertices);
    541         IssmDouble* tmp=xNew<IssmDouble>(numvertices);
     535        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
     536        const int NUM_VERTICES_MONTHS_PER_YEAR  = NUM_VERTICES * 12;
     537
     538        int             i;
     539        IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     540        IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     541        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     542        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     543        IssmDouble* TemperaturesReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     544        IssmDouble* PrecipitationsReconstructed=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     545        IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
    542546        IssmDouble Delta18oTime;
    543547        IssmDouble dpermil,f;
    … …  
    583587        Gauss* gauss=this->NewGauss();
    584588        for(int month=0;month<12;month++) {
    585                 for(int iv=0;iv<numvertices;iv++) {
     589                for(int iv=0;iv<NUM_VERTICES;iv++) {
    586590                        gauss->GaussVertex(iv);
    587591                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,time_clim+month/12.*yts);
    … …  
    604608
    605609        /*Compute the temperature and precipitation*/
    606         for(int iv=0;iv<numvertices;iv++){
     610        for(int iv=0;iv<NUM_VERTICES;iv++){
    607611                ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,isTemperatureScaled,isPrecipScaled,
    608612                                        f,&PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
    … …  
    615619        TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
    616620        for (int imonth=0;imonth<12;imonth++) {
    617                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
     621                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    618622                switch(this->ObjectEnum()){
    619623                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    622626                        default: _error_("Not implemented yet");
    623627                }
    624                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
     628                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    625629                switch(this->ObjectEnum()){
    626630                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    662666        /*Are we on the base? If not, return*/
    663667        if(!IsOnBase()) return;
    664         int        numvertices = this->GetNumberOfVertices();
     668        const int NUM_VERTICES = this->GetNumberOfVertices();
    665669
    666670        int        i;
    … …  
    671675        IssmDouble time;
    672676
    673         IssmDouble*             smb              = xNew<IssmDouble>(numvertices);
    674         IssmDouble*             surf     = xNew<IssmDouble>(numvertices);
    675         IssmDouble*             accu     = xNew<IssmDouble>(numvertices);
    676         IssmDouble*             runoff = xNew<IssmDouble>(numvertices);
     677        IssmDouble*             smb             = xNew<IssmDouble>(NUM_VERTICES);
     678        IssmDouble*             surf    = xNew<IssmDouble>(NUM_VERTICES);
     679        IssmDouble*             accu    = xNew<IssmDouble>(NUM_VERTICES);
     680        IssmDouble*             runoff  = xNew<IssmDouble>(NUM_VERTICES);
    677681
    678682        /*Get material parameters :*/
    … …  
    695699
    696700        /*Compute the temperature and precipitation*/
    697         for(int iv=0;iv<numvertices;iv++){
     701        for(int iv=0;iv<NUM_VERTICES;iv++){
    698702                accu[iv]=max(0.,(accuref+(surf[iv]-accualti)*accugrad));
    699703                runoff[iv]=max(0.,(runoffref+(surf[iv]-runoffalti)*runoffgrad));
    … …  
    914918        _printf_("   sid: "<<this->sid<<"\n");
    915919        if(vertices){
    916                 int numvertices = this->GetNumberOfVertices();
    917                 for(int i=0;i<numvertices;i++) vertices[i]->Echo();
     920                const int NUM_VERTICES = this->GetNumberOfVertices();
     921                for(int i=0;i<NUM_VERTICES;i++) vertices[i]->Echo();
    918922        }
    919923        else _printf_("vertices = NULL\n");
    … …  
    11731177
    11741178        /*Fetch number vertices for this element*/
    1175         int numvertices = this->GetNumberOfVertices();
     1179        const int NUM_VERTICES = this->GetNumberOfVertices();
    11761180
    11771181        /*Checks in debugging mode*/
    … …  
    11801184        /* Start looping on the number of vertices: */
    11811185        Gauss*gauss=this->NewGauss();
    1182         for(int iv=0;iv<numvertices;iv++){
     1186        for(int iv=0;iv<NUM_VERTICES;iv++){
    11831187                gauss->GaussVertex(iv);
    11841188                input->GetInputValue(&pvalue[iv],gauss);
    … …  
    11961200
    11971201        /*Fetch number vertices for this element*/
    1198         int numvertices = this->GetNumberOfVertices();
     1202        const int NUM_VERTICES = this->GetNumberOfVertices();
    11991203
    12001204        /*Checks in debugging mode*/
    … …  
    12031207        /* Start looping on the number of vertices: */
    12041208        Gauss*gauss=this->NewGauss();
    1205         for(int iv=0;iv<numvertices;iv++){
     1209        for(int iv=0;iv<NUM_VERTICES;iv++){
    12061210                gauss->GaussVertex(iv);
    12071211                input->GetInputValue(&pvalue[iv],gauss,time);
    … …  
    12211225
    12221226        /*Fetch number vertices for this element*/
    1223         int numvertices = this->GetNumberOfVertices();
     1227        const int NUM_VERTICES = this->GetNumberOfVertices();
    12241228
    12251229        /* Start looping on the number of vertices: */
    12261230        if (input){
    12271231                Gauss* gauss=this->NewGauss();
    1228                 for (int iv=0;iv<numvertices;iv++){
     1232                for (int iv=0;iv<NUM_VERTICES;iv++){
    12291233                        gauss->GaussVertex(iv);
    12301234                        input->GetInputValue(&pvalue[iv],gauss);
    … …  
    12331237        }
    12341238        else{
    1235                 for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
     1239                for(int iv=0;iv<NUM_VERTICES;iv++) pvalue[iv]=defaultvalue;
    12361240        }
    12371241}
    … …  
    14061410
    14071411/*      /\*Fetch number vertices for this element and allocate arrays*\/ */
    1408 /*      int numvertices = this->GetNumberOfVertices(); */
    1409 /*      int*        vertexpidlist = xNew<int>(numvertices); */
    1410 /*      IssmDouble* values        = xNew<IssmDouble>(numvertices); */
     1412/*      const int NUM_VERTICES = this->GetNumberOfVertices(); */
     1413/*      int*        vertexpidlist = xNew<int>(NUM_VERTICES); */
     1414/*      IssmDouble* values        = xNew<IssmDouble>(NUM_VERTICES); */
    14111415
    14121416/*      /\*Fill in values*\/ */
    14131417/*      this->GetVerticesPidList(vertexpidlist); */
    14141418/*      this->GetInputListOnVertices(values,input_enum); */
    1415 /*      vector->SetValues(numvertices,vertexpidlist,values,INS_VAL); */
     1419/*      vector->SetValues(NUM_VERTICES,vertexpidlist,values,INS_VAL); */
    14161420
    14171421/*      /\*Clean up*\/ */
    … …  
    14241428
    14251429        /*Fetch number vertices for this element and allocate arrays*/
    1426         int         numvertices = this->GetNumberOfVertices();
     1430        const int NUM_VERTICES = this->GetNumberOfVertices();
     1431
    14271432        int         numnodes    = this->GetNumberOfNodes();
    14281433        int*        doflist     = NULL;
    … …  
    14381443                        break;
    14391444                case VertexPIdEnum:
    1440                         doflist = xNew<int>(numvertices);
    1441                         values = xNew<IssmDouble>(numvertices);
     1445                        doflist = xNew<int>(NUM_VERTICES);
     1446                        values = xNew<IssmDouble>(NUM_VERTICES);
    14421447                        /*Fill in values*/
    14431448                        this->GetVerticesPidList(doflist);
    14441449                        this->GetInputListOnVertices(values,input_enum);
    1445                         vector->SetValues(numvertices,doflist,values,INS_VAL);
     1450                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    14461451                        break;
    14471452                case VertexSIdEnum:
    1448                         doflist = xNew<int>(numvertices);
    1449                         values = xNew<IssmDouble>(numvertices);
     1453                        doflist = xNew<int>(NUM_VERTICES);
     1454                        values = xNew<IssmDouble>(NUM_VERTICES);
    14501455                        /*Fill in values*/
    14511456                        this->GetVerticesSidList(doflist);
    14521457                        this->GetInputListOnVertices(values,input_enum);
    1453                         vector->SetValues(numvertices,doflist,values,INS_VAL);
     1458                        vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    14541459                        break;
    14551460                case NodesEnum:
    … …  
    14821487
    14831488        /*Fetch number vertices for this element and allocate arrays*/
    1484         int         numvertices = this->GetNumberOfVertices();
     1489        const int NUM_VERTICES = this->GetNumberOfVertices();
     1490
    14851491        int         numnodes    = this->GetNumberOfNodes();
    14861492        int*        doflist     = NULL;
    … …  
    14891495        switch(type){
    14901496        case VertexPIdEnum:
    1491                 doflist = xNew<int>(numvertices);
    1492                 values = xNew<IssmDouble>(numvertices);
     1497                doflist = xNew<int>(NUM_VERTICES);
     1498                values = xNew<IssmDouble>(NUM_VERTICES);
    14931499                /*Fill in values*/
    14941500                this->GetVerticesPidList(doflist);
    14951501                this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1496                 vector->SetValues(numvertices,doflist,values,INS_VAL);
     1502                vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    14971503                break;
    14981504        case VertexSIdEnum:
    1499                 doflist = xNew<int>(numvertices);
    1500                 values = xNew<IssmDouble>(numvertices);
     1505                doflist = xNew<int>(NUM_VERTICES);
     1506                values = xNew<IssmDouble>(NUM_VERTICES);
    15011507                /*Fill in values*/
    15021508                this->GetVerticesSidList(doflist);
    15031509                this->GetInputListOnVerticesAtTime(values,input_enum,time);
    1504                 vector->SetValues(numvertices,doflist,values,INS_VAL);
     1510                vector->SetValues(NUM_VERTICES,doflist,values,INS_VAL);
    15051511                break;
    15061512        default:
    … …  
    15161522void       Element::GetVerticesLidList(int* lidlist){/*{{{*/
    15171523
    1518         int numvertices = this->GetNumberOfVertices();
    1519         for(int i=0;i<numvertices;i++) lidlist[i]=vertices[i]->Lid();
     1524        const int NUM_VERTICES = this->GetNumberOfVertices();
     1525        for(int i=0;i<NUM_VERTICES;i++) lidlist[i]=vertices[i]->Lid();
    15201526
    15211527}
    … …  
    15231529void       Element::GetVerticesPidList(int* pidlist){/*{{{*/
    15241530
    1525         int numvertices = this->GetNumberOfVertices();
    1526         for(int i=0;i<numvertices;i++) pidlist[i]=vertices[i]->Pid();
     1531        const int NUM_VERTICES = this->GetNumberOfVertices();
     1532        for(int i=0;i<NUM_VERTICES;i++) pidlist[i]=vertices[i]->Pid();
    15271533
    15281534}
    … …  
    15301536void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
    15311537
    1532         int numvertices = this->GetNumberOfVertices();
    1533         for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
     1538        const int NUM_VERTICES = this->GetNumberOfVertices();
     1539        for(int i=0;i<NUM_VERTICES;i++) connectivity[i]=this->vertices[i]->Connectivity();
    15341540}
    15351541/*}}}*/
    15361542void       Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/
    15371543
    1538         int         numvertices = this->GetNumberOfVertices();
    1539         IssmDouble* xyz_list    = xNew<IssmDouble>(numvertices*3);
    1540         ::GetVerticesCoordinates(xyz_list,this->vertices,numvertices);
     1544        const int NUM_VERTICES = this->GetNumberOfVertices();
     1545
     1546        IssmDouble* xyz_list    = xNew<IssmDouble>(NUM_VERTICES*3);
     1547        ::GetVerticesCoordinates(xyz_list,this->vertices,NUM_VERTICES);
    15411548
    15421549        *pxyz_list = xyz_list;
    … …  
    15451552void       Element::GetVerticesSidList(int* sidlist){/*{{{*/
    15461553
    1547         int numvertices = this->GetNumberOfVertices();
    1548         for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid();
     1554        const int NUM_VERTICES = this->GetNumberOfVertices();
     1555        for(int i=0;i<NUM_VERTICES;i++) sidlist[i]=this->vertices[i]->Sid();
    15491556}
    15501557/*}}}*/
    … …  
    15551562
    15561563        /*Create list of x*/
    1557         int         numvertices = this->GetNumberOfVertices();
    1558         IssmDouble* x_list      = xNew<IssmDouble>(numvertices);
    1559 
    1560         for(int i=0;i<numvertices;i++) x_list[i]=xyz_list[i*3+0];
     1564        const int NUM_VERTICES = this->GetNumberOfVertices();
     1565
     1566        IssmDouble* x_list = xNew<IssmDouble>(NUM_VERTICES);
     1567
     1568        for(int i=0;i<NUM_VERTICES;i++) x_list[i]=xyz_list[i*3+0];
    15611569        ValueP1OnGauss(&x,x_list,gauss);
    15621570
    … …  
    15701578
    15711579        /*Create list of y*/
    1572         int         numvertices = this->GetNumberOfVertices();
    1573         IssmDouble* y_list      = xNew<IssmDouble>(numvertices);
    1574 
    1575         for(int i=0;i<numvertices;i++) y_list[i]=xyz_list[i*3+1];
     1580        const int NUM_VERTICES = this->GetNumberOfVertices();
     1581
     1582        IssmDouble* y_list      = xNew<IssmDouble>(NUM_VERTICES);
     1583
     1584        for(int i=0;i<NUM_VERTICES;i++) y_list[i]=xyz_list[i*3+1];
    15761585        ValueP1OnGauss(&y,y_list,gauss);
    15771586
    … …  
    15851594
    15861595        /*Create list of z*/
    1587         int         numvertices = this->GetNumberOfVertices();
    1588         IssmDouble* z_list      = xNew<IssmDouble>(numvertices);
    1589 
    1590         for(int i=0;i<numvertices;i++) z_list[i]=xyz_list[i*3+2];
     1596        const int NUM_VERTICES = this->GetNumberOfVertices();
     1597
     1598        IssmDouble* z_list      = xNew<IssmDouble>(NUM_VERTICES);
     1599
     1600        for(int i=0;i<NUM_VERTICES;i++) z_list[i]=xyz_list[i*3+2];
    15911601        ValueP1OnGauss(&z,z_list,gauss);
    15921602
    … …  
    16031613
    16041614        /*Get number of vertices*/
    1605         int numvertices = this->GetNumberOfVertices();
     1615        const int NUM_VERTICES = this->GetNumberOfVertices();
    16061616        if(isautodiff){
    16071617                int* N=NULL;
    … …  
    16171627
    16181628                for(int n=0;n<N[control_index];n++){
    1619                         for(int i=0;i<numvertices;i++){
    1620                                 indexing[i+n*numvertices]=this->vertices[i]->Sid() + start + n*M[control_index];
     1629                        for(int i=0;i<NUM_VERTICES;i++){
     1630                                indexing[i+n*NUM_VERTICES]=this->vertices[i]->Sid() + start + n*M[control_index];
    16211631                        }
    16221632                }
    … …  
    16261636                parameters->FindParam(&M,ControlInputSizeMEnum);
    16271637                /*get gradient indices*/
    1628                 for(int i=0;i<numvertices;i++){
     1638                for(int i=0;i<NUM_VERTICES;i++){
    16291639                        indexing[i]=this->vertices[i]->Sid() + (control_index)*M;
    16301640                }
    … …  
    17201730    if(vector_type==1){ //nodal vector
    17211731
    1722         int         numvertices = this->GetNumberOfVertices();
    1723         int        *vertexids   = xNew<int>(numvertices);
    1724         IssmDouble *values      = xNew<IssmDouble>(numvertices);
     1732        const int NUM_VERTICES = this->GetNumberOfVertices();
     1733
     1734        int        *vertexids   = xNew<int>(NUM_VERTICES);
     1735        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
    17251736
    17261737        /*Recover vertices ids needed to initialize inputs*/
    17271738        _assert_(iomodel->elements);
    1728         for(i=0;i<numvertices;i++){
    1729             vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1739        for(i=0;i<NUM_VERTICES;i++){
     1740            vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    17301741        }
    17311742
    … …  
    17361747                  }
    17371748                  else if(M==iomodel->numberofvertices){
    1738             for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
     1749            for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    17391750            this->AddInput(vector_enum,values,P1Enum);
    17401751        }
    … …  
    17451756            TransientInput* transientinput=new TransientInput(vector_enum,times,N);
    17461757            for(t=0;t<N;t++){
    1747                 for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1758                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    17481759                switch(this->ObjectEnum()){
    17491760                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
    … …  
    18331844
    18341845        /*Intermediaries*/
    1835         int         numvertices = this->GetNumberOfVertices();
    1836         int        *vertexids   = xNew<int>(numvertices);
    1837         IssmDouble *values      = xNew<IssmDouble>(numvertices);
    1838         IssmDouble *values_min  = xNew<IssmDouble>(numvertices);
    1839         IssmDouble *values_max  = xNew<IssmDouble>(numvertices);
     1846        const int NUM_VERTICES = this->GetNumberOfVertices();
     1847
     1848        int        *vertexids   = xNew<int>(NUM_VERTICES);
     1849        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
     1850        IssmDouble *values_min  = xNew<IssmDouble>(NUM_VERTICES);
     1851        IssmDouble *values_max  = xNew<IssmDouble>(NUM_VERTICES);
    18401852
    18411853        /*Some sanity checks*/
    … …  
    18501862        /*Recover vertices ids needed to initialize inputs*/
    18511863        _assert_(iomodel->elements);
    1852         for(int i=0;i<numvertices;i++){
    1853                 vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1864        for(int i=0;i<NUM_VERTICES;i++){
     1865                vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    18541866        }
    18551867
    18561868        /*Are we in transient or static? */
    18571869        if(M==iomodel->numberofvertices){
    1858                 for(int i=0;i<numvertices;i++){
     1870                for(int i=0;i<NUM_VERTICES;i++){
    18591871                        values[i]=vector[vertexids[i]-1];
    18601872                        values_min[i] = min_vector[vertexids[i]-1];
    … …  
    18741886                           TransientInput* grad_input = new TransientInput(ControlInputGradEnum);
    18751887                                for(int t=0;t<N;t++){
    1876                 for(int i=0;i<numvertices;i++){
     1888                for(int i=0;i<NUM_VERTICES;i++){
    18771889                                                values[i]=vector[N*(vertexids[i]-1)+t];
    18781890                                                values_min[i] = min_vector[N*(vertexids[i]-1)+t];
    … …  
    19351947    if(vector_type==1){ //nodal vector
    19361948
    1937         int         numvertices = this->GetNumberOfVertices();
    1938         int        *vertexids   = xNew<int>(numvertices);
    1939         IssmDouble *values      = xNew<IssmDouble>(numvertices);
     1949        const int NUM_VERTICES = this->GetNumberOfVertices();
     1950
     1951        int        *vertexids   = xNew<int>(NUM_VERTICES);
     1952        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
    19401953
    19411954        /*Recover vertices ids needed to initialize inputs*/
    19421955        _assert_(iomodel->elements);
    1943         for(i=0;i<numvertices;i++){
    1944             vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
     1956        for(i=0;i<NUM_VERTICES;i++){
     1957            vertexids[i]=reCast<int>(iomodel->elements[NUM_VERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
    19451958        }
    19461959
    … …  
    19561969                  }
    19571970                  else if(M==iomodel->numberofvertices){
    1958             for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
     1971            for(i=0;i<NUM_VERTICES;i++) values[i]=vector[vertexids[i]-1];
    19591972                                switch(this->ObjectEnum()){
    19601973                    case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
    … …  
    19691982            TransientInput* transientinput=new TransientInput(input_enum,times,N);
    19701983            for(t=0;t<N;t++){
    1971                 for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
     1984                for(i=0;i<NUM_VERTICES;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    19721985                switch(this->ObjectEnum()){
    19731986                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break;
    … …  
    21612174        if(!this->IsIceInElement() || !this->IsFloating()) return;
    21622175
    2163         int         numvertices = this->GetNumberOfVertices();
     2176        const int NUM_VERTICES = this->GetNumberOfVertices();
     2177
    21642178        int         basinid,num_basins,M,N;
    21652179        IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
    2166         IssmDouble  basalmeltrate[numvertices];
     2180        IssmDouble  basalmeltrate[NUM_VERTICES];
    21672181        bool        islocal;
    21682182        IssmDouble* delta_t = NULL;
    … …  
    21942208        /*Compute melt rate for Local and Nonlocal parameterizations*/
    21952209        Gauss* gauss=this->NewGauss();
    2196         for(int i=0;i<numvertices;i++){
     2210        for(int i=0;i<NUM_VERTICES;i++){
    21972211                gauss->GaussVertex(i);
    21982212                tf_input->GetInputValue(&tf,gauss);
    … …  
    22202234void       Element::LinearFloatingiceMeltingRate(){/*{{{*/
    22212235
    2222         int numvertices      = this->GetNumberOfVertices();
     2236        const int NUM_VERTICES = this->GetNumberOfVertices();
     2237
    22232238        IssmDouble  deepwaterel,upperwaterel,deepwatermelt,upperwatermelt;
    2224         IssmDouble* base     = xNew<IssmDouble>(numvertices);
    2225         IssmDouble* values   = xNew<IssmDouble>(numvertices);
    2226         IssmDouble time;
     2239        IssmDouble* base                = xNew<IssmDouble>(NUM_VERTICES);
     2240        IssmDouble* values              = xNew<IssmDouble>(NUM_VERTICES);
     2241        IssmDouble      time;
    22272242
    22282243        parameters->FindParam(&time,TimeEnum);
    … …  
    22312246        parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum,time);
    22322247        parameters->FindParam(&upperwatermelt,BasalforcingsUpperwaterMeltingRateEnum,time);
    2233         _assert_(upperwaterel>deepwaterel); 
     2248        _assert_(upperwaterel>deepwaterel);
    22342249
    22352250        this->GetInputListOnVertices(base,BaseEnum);
    2236         for(int i=0;i<numvertices;i++){
     2251        for(int i=0;i<NUM_VERTICES;i++){
    22372252                if(base[i]>=upperwaterel){
    22382253                        values[i]=upperwatermelt;
    … …  
    22532268void       Element::SpatialLinearFloatingiceMeltingRate(){/*{{{*/
    22542269
    2255         int numvertices      = this->GetNumberOfVertices();
    2256         IssmDouble* deepwatermelt     = xNew<IssmDouble>(numvertices);
    2257         IssmDouble* deepwaterel     = xNew<IssmDouble>(numvertices);
    2258         IssmDouble* upperwaterel     = xNew<IssmDouble>(numvertices);
    2259         IssmDouble* base     = xNew<IssmDouble>(numvertices);
    2260         IssmDouble* values   = xNew<IssmDouble>(numvertices);
     2270        const int NUM_VERTICES = this->GetNumberOfVertices();
     2271
     2272        IssmDouble* deepwatermelt       = xNew<IssmDouble>(NUM_VERTICES);
     2273        IssmDouble* deepwaterel     = xNew<IssmDouble>(NUM_VERTICES);
     2274        IssmDouble* upperwaterel        = xNew<IssmDouble>(NUM_VERTICES);
     2275        IssmDouble* base                        = xNew<IssmDouble>(NUM_VERTICES);
     2276        IssmDouble* values                      = xNew<IssmDouble>(NUM_VERTICES);
    22612277
    22622278        this->GetInputListOnVertices(base,BaseEnum);
    … …  
    22652281        this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
    22662282
    2267         for(int i=0;i<numvertices;i++){
     2283        for(int i=0;i<NUM_VERTICES;i++){
    22682284                if(base[i]>upperwaterel[i])      values[i]=0;
    22692285                else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
    … …  
    22812297void       Element::MantlePlumeGeothermalFlux(){/*{{{*/
    22822298
    2283         int numvertices      = this->GetNumberOfVertices();
     2299        const int NUM_VERTICES = this->GetNumberOfVertices();
    22842300        IssmDouble  mantleconductivity,nusselt,dtbg,plumeradius,topplumedepth,bottomplumedepth,plumex,plumey;
    22852301        IssmDouble  crustthickness,uppercrustthickness,uppercrustheat,lowercrustheat;
    22862302        IssmDouble  crustheat,plumeheat,dt,middleplumedepth,a,e,eprime,A0,lambda,Alambda,dAlambda;
    22872303        IssmDouble  x,y,z,c;
    2288         IssmDouble* values   = xNew<IssmDouble>(numvertices);
     2304        IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
    22892305        IssmDouble *xyz_list = NULL;
    22902306
    … …  
    23072323        e=pow(a*a-c*c,1./2.)/a;
    23082324        A0=(1-pow(e,2.))/pow(e,3.)*(1./2.*log((1+e)/(1-e))-e);
    2309         for(int i=0;i<numvertices;i++){
     2325        for(int i=0;i<NUM_VERTICES;i++){
    23102326                y=xyz_list[i*3+0]-plumex;
    23112327                z=xyz_list[i*3+1]-plumey;
    … …  
    23462362void       Element::MigrateGroundingLine(IssmDouble* phi_ungrounding){/*{{{*/
    23472363
    2348         int         numvertices = this->GetNumberOfVertices();
     2364        const int  NUM_VERTICES = this->GetNumberOfVertices();
    23492365        int        i,migration_style;
    23502366        IssmDouble bed_hydro,yts;
    23512367        IssmDouble rho_water,rho_ice,density;
    2352         IssmDouble* melting = xNew<IssmDouble>(numvertices);
    2353         IssmDouble* phi     = xNew<IssmDouble>(numvertices);
    2354         IssmDouble* h       = xNew<IssmDouble>(numvertices);
    2355         IssmDouble* s       = xNew<IssmDouble>(numvertices);
    2356         IssmDouble* b       = xNew<IssmDouble>(numvertices);
    2357         IssmDouble* r       = xNew<IssmDouble>(numvertices);
    2358         IssmDouble* sl      = xNew<IssmDouble>(numvertices);
     2368        IssmDouble* melting = xNew<IssmDouble>(NUM_VERTICES);
     2369        IssmDouble* phi     = xNew<IssmDouble>(NUM_VERTICES);
     2370        IssmDouble* h       = xNew<IssmDouble>(NUM_VERTICES);
     2371        IssmDouble* s       = xNew<IssmDouble>(NUM_VERTICES);
     2372        IssmDouble* b       = xNew<IssmDouble>(NUM_VERTICES);
     2373        IssmDouble* r       = xNew<IssmDouble>(NUM_VERTICES);
     2374        IssmDouble* sl      = xNew<IssmDouble>(NUM_VERTICES);
    23592375
    23602376        /*Recover info at the vertices: */
    … …  
    23722388
    23732389        /*go through vertices, and update inputs, considering them to be TriaVertex type: */
    2374         for(i=0;i<numvertices;i++){
     2390        for(i=0;i<NUM_VERTICES;i++){
    23752391                /* Contact FS*/
    23762392                if(migration_style == ContactEnum){
    … …  
    24082424
    24092425        /*Recalculate phi*/
    2410         for(i=0;i<numvertices;i++){
     2426        for(i=0;i<NUM_VERTICES;i++){
    24112427                if(migration_style==SoftMigrationEnum){
    24122428                        bed_hydro=-density*h[i]+sl[i];
    … …  
    24382454/*}}}*/
    24392455void       Element::MismipFloatingiceMeltingRate(){/*{{{*/
    2440 
    2441         int numvertices      = this->GetNumberOfVertices();
     2456        const int NUM_VERTICES = this->GetNumberOfVertices();
     2457
    24422458        IssmDouble  meltratefactor,thresholdthickness,upperdepthmelt;
    2443         IssmDouble* base     = xNew<IssmDouble>(numvertices);
    2444         IssmDouble* bed      = xNew<IssmDouble>(numvertices);
    2445         IssmDouble* values   = xNew<IssmDouble>(numvertices);
     2459        IssmDouble* base     = xNew<IssmDouble>(NUM_VERTICES);
     2460        IssmDouble* bed      = xNew<IssmDouble>(NUM_VERTICES);
     2461        IssmDouble* values   = xNew<IssmDouble>(NUM_VERTICES);
    24462462
    24472463        parameters->FindParam(&meltratefactor,BasalforcingsMeltrateFactorEnum);
    … …  
    24512467        this->GetInputListOnVertices(base,BaseEnum);
    24522468        this->GetInputListOnVertices(bed,BedEnum);
    2453         for(int i=0;i<numvertices;i++){
     2469        for(int i=0;i<NUM_VERTICES;i++){
    24542470                if(base[i]>upperdepthmelt){
    24552471                        values[i]=0;
    … …  
    24702486        if(!IsOnBase()) return;
    24712487
    2472         int        numvertices = this->GetNumberOfVertices();
    2473 
    2474         int        i;
    2475         IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
    2476         IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    2477         IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
    2478         IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
    2479         IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
    2480         IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(12*numvertices);
    2481         IssmDouble* tmp=xNew<IssmDouble>(numvertices);
     2488        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
     2489        const int NUM_VERTICES_MONTHS_PER_YEAR  = NUM_VERTICES * 12;
     2490
     2491        int i;
     2492        IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2493        IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2494        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2495        IssmDouble* TemperaturesLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2496        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2497        IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2498        IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
    24822499        IssmDouble TdiffTime,PfacTime;
    24832500        IssmDouble time,yts,time_yr;
    … …  
    24942511        Gauss* gauss=this->NewGauss();
    24952512        for(int month=0;month<12;month++) {
    2496                 for(int iv=0;iv<numvertices;iv++) {
     2513                for(int iv=0;iv<NUM_VERTICES;iv++) {
    24972514                        gauss->GaussVertex(iv);
    24982515                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
    … …  
    25112528
    25122529        /*Compute the temperature and precipitation*/
    2513         for(int iv=0;iv<numvertices;iv++){
     2530        for(int iv=0;iv<NUM_VERTICES;iv++){
    25142531                ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
    25152532                                        &PrecipitationsLgm[iv*12],&PrecipitationsPresentday[iv*12],
    … …  
    25222539        TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
    25232540        for (int imonth=0;imonth<12;imonth++) {
    2524                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
     2541                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    25252542                switch(this->ObjectEnum()){
    25262543                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    25292546                        default: _error_("Not implemented yet");
    25302547                }
    2531                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
     2548                for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    25322549                switch(this->ObjectEnum()){
    25332550                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    26162633
    26172634        this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));
    2618        
     2635
    26192636}/*}}}*/
    26202637void       Element::PicoUpdateBox(int loopboxid){/*{{{*/
    … …  
    26262643        if(loopboxid!=boxid) return;
    26272644
    2628         int        NUMVERTICES = this->GetNumberOfVertices();
     2645        const int NUM_VERTICES = this->GetNumberOfVertices();
     2646
    26292647        int        basinid, maxbox, num_basins, numnodes, M;
    26302648        IssmDouble gamma_T, overturning_coeff, thickness;
    … …  
    26642682        IssmDouble g1               = area_boxi*gamma_T;
    26652683
    2666         IssmDouble basalmeltrates_shelf[NUMVERTICES];
    2667         IssmDouble potential_pressure_melting_point[NUMVERTICES];
    2668         IssmDouble Tocs[NUMVERTICES];
    2669         IssmDouble Socs[NUMVERTICES];
     2684        IssmDouble basalmeltrates_shelf[NUM_VERTICES];
     2685        IssmDouble potential_pressure_melting_point[NUM_VERTICES];
     2686        IssmDouble Tocs[NUM_VERTICES];
     2687        IssmDouble Socs[NUM_VERTICES];
    26702688
    26712689        /* First box calculations */
    … …  
    26772695                this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
    26782696                IssmDouble s1 = soc_farocean/(nu*lambda);
    2679                 IssmDouble overturnings[NUMVERTICES];
     2697                IssmDouble overturnings[NUM_VERTICES];
    26802698
    26812699                /* Start looping on the number of verticies and calculate ocean vars */
    26822700                Gauss* gauss=this->NewGauss();
    2683                 for(int i=0;i<NUMVERTICES;i++){
     2701                for(int i=0;i<NUM_VERTICES;i++){
    26842702                        gauss->GaussVertex(i);
    26852703                        thickness_input->GetInputValue(&thickness,gauss);
    … …  
    27242742                /* Start looping on the number of verticies and calculate ocean vars */
    27252743                Gauss* gauss=this->NewGauss();
    2726                 for(int i=0;i<NUMVERTICES;i++){
     2744                for(int i=0;i<NUM_VERTICES;i++){
    27272745                        gauss->GaussVertex(i);
    27282746                        thickness_input->GetInputValue(&thickness,gauss);
    … …  
    27482766        /*Cleanup and return*/
    27492767        xDelete<IssmDouble>(boxareas);
    2750        
     2768
    27512769}/*}}}*/
    27522770void       Element::PicoComputeBasalMelt(){/*{{{*/
    … …  
    27542772        if(!this->IsIceInElement() || !this->IsFloating()) return;
    27552773
    2756         int        NUMVERTICES = this->GetNumberOfVertices();
     2774        const int NUM_VERTICES = this->GetNumberOfVertices();
     2775
    27572776        IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
    27582777        IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
    … …  
    27742793
    27752794        /*Define arrays*/
    2776         IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
     2795        IssmDouble basalmeltrates_shelf[NUM_VERTICES];  //Basal melt-rate
    27772796
    27782797        /*Polynomial coefficients*/
    … …  
    28022821        /*Loop over the number of vertices in this element*/
    28032822        Gauss* gauss=this->NewGauss();
    2804         for(int i=0;i<NUMVERTICES;i++){
     2823        for(int i=0;i<NUM_VERTICES;i++){
    28052824                gauss->GaussVertex(i);
    28062825
    … …  
    28512870        /*Cleanup and return*/
    28522871        delete gauss;
    2853                
     2872
    28542873}/*}}}*/
    28552874void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
    28562875
    2857         int  numvertices = this->GetNumberOfVertices();
    2858 
    2859         int        i;
    2860         IssmDouble* agd=xNew<IssmDouble>(numvertices); // surface mass balance
    2861         IssmDouble* melt=xNew<IssmDouble>(numvertices); // surface mass balance
    2862         IssmDouble* accu=xNew<IssmDouble>(numvertices); // surface mass balance
    2863         IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
    2864         IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    2865         IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
    2866         IssmDouble* tmp=xNew<IssmDouble>(numvertices);
    2867         IssmDouble* h=xNew<IssmDouble>(numvertices);
    2868         IssmDouble* s=xNew<IssmDouble>(numvertices);
    2869         IssmDouble* s0p=xNew<IssmDouble>(numvertices);
    2870         IssmDouble* s0t=xNew<IssmDouble>(numvertices);
     2876        const int NUM_VERTICES          = this->GetNumberOfVertices();
     2877        const int NUM_VERTICES_MONTHS_PER_YEAR = NUM_VERTICES * 12;
     2878
     2879        int             i;
     2880        IssmDouble* agd=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
     2881        IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
     2882        IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance
     2883        IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2884        IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     2885        IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
     2886        IssmDouble* tmp=xNew<IssmDouble>(NUM_VERTICES);
     2887        IssmDouble* h=xNew<IssmDouble>(NUM_VERTICES);
     2888        IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
     2889        IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);
     2890        IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);
    28712891        IssmDouble pddsnowfac = -1.;
    28722892        IssmDouble pddicefac = -1.;
    … …  
    28952915        Gauss* gauss=this->NewGauss();
    28962916        for(int month=0;month<12;month++) {
    2897                 for(int iv=0;iv<numvertices;iv++) {
     2917                for(int iv=0;iv<NUM_VERTICES;iv++) {
    28982918                        gauss->GaussVertex(iv);
    28992919                        input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,time_yr+month/12.*yts);
    … …  
    29322952
    29332953        /*measure the surface mass balance*/
    2934         for (int iv = 0; iv<numvertices; iv++){
     2954        for (int iv = 0; iv<NUM_VERTICES; iv++){
    29352955                gauss->GaussVertex(iv);
    29362956                pddsnowfac=-1.;
    … …  
    29522972        // TransientInput* NewPrecipitationInput = new TransientInput(SmbPrecipitationEnum);
    29532973        // for (int imonth=0;imonth<12;imonth++) {
    2954         //   for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
     2974        //   for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    29552975        //   TriaInput* newmonthinput1 = new TriaInput(SmbMonthlytemperaturesEnum,&tmp[0],P1Enum);
    29562976        //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    29572977        //
    2958         //   for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
     2978        //   for(i=0;i<NUM_VERTICES;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    29592979        //   TriaInput* newmonthinput2 = new TriaInput(SmbPrecipitationEnum,&tmp[0],P1Enum);
    29602980        //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    … …  
    30423062        /* General FIXMEs: get Tmelting point, pddicefactor, pddsnowfactor, sigma from parameters/user input */
    30433063
    3044         int  numvertices = this->GetNumberOfVertices();
    3045 
    3046         int        i;
    3047         IssmDouble* smb=xNew<IssmDouble>(numvertices);          // surface mass balance
    3048         IssmDouble* melt=xNew<IssmDouble>(numvertices);         // melting comp. of surface mass balance
    3049         IssmDouble* accu=xNew<IssmDouble>(numvertices);         // accuumulation comp. of surface mass balance
    3050         IssmDouble* melt_star=xNew<IssmDouble>(numvertices);
    3051         IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
    3052         IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
    3053         IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
    3054         IssmDouble* s=xNew<IssmDouble>(numvertices);                    // actual surface height
    3055         IssmDouble* s0p=xNew<IssmDouble>(numvertices);          // reference elevation for precip.
    3056         IssmDouble* s0t=xNew<IssmDouble>(numvertices);          // reference elevation for temperature
    3057         IssmDouble* smbcorr=xNew<IssmDouble>(numvertices); // surface mass balance correction; will be added after pdd call
    3058         IssmDouble* p_ampl=xNew<IssmDouble>(numvertices);       // precip anomaly
    3059         IssmDouble* t_ampl=xNew<IssmDouble>(numvertices);       // remperature anomaly
     3064        const int NUM_VERTICES          = this->GetNumberOfVertices();
     3065        const int NUM_VERTICES_MONTHS_PER_YEAR  = NUM_VERTICES * 12;
     3066
     3067        int             i;
     3068        IssmDouble* smb=xNew<IssmDouble>(NUM_VERTICES);         // surface mass balance
     3069        IssmDouble* melt=xNew<IssmDouble>(NUM_VERTICES);                // melting comp. of surface mass balance
     3070        IssmDouble* accu=xNew<IssmDouble>(NUM_VERTICES);                // accuumulation comp. of surface mass balance
     3071        IssmDouble* melt_star=xNew<IssmDouble>(NUM_VERTICES);
     3072        IssmDouble* monthlytemperatures=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     3073        IssmDouble* monthlyprec=xNew<IssmDouble>(NUM_VERTICES_MONTHS_PER_YEAR);
     3074        IssmDouble* yearlytemperatures=xNew<IssmDouble>(NUM_VERTICES); memset(yearlytemperatures, 0., NUM_VERTICES*sizeof(IssmDouble));
     3075        IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);                   // actual surface height
     3076        IssmDouble* s0p=xNew<IssmDouble>(NUM_VERTICES);         // reference elevation for precip.
     3077        IssmDouble* s0t=xNew<IssmDouble>(NUM_VERTICES);         // reference elevation for temperature
     3078        IssmDouble* smbcorr=xNew<IssmDouble>(NUM_VERTICES); // surface mass balance correction; will be added after pdd call
     3079        IssmDouble* p_ampl=xNew<IssmDouble>(NUM_VERTICES);      // precip anomaly
     3080        IssmDouble* t_ampl=xNew<IssmDouble>(NUM_VERTICES);      // remperature anomaly
    30603081        IssmDouble rho_water,rho_ice,desfac,rlaps;
    30613082        IssmDouble inv_twelve=1./12.;                                                           //factor for monthly average
    … …  
    30883109        Gauss* gauss=this->NewGauss();
    30893110        for(int month=0;month<12;month++){
    3090                 for(int iv=0;iv<numvertices;iv++){
     3111                for(int iv=0;iv<NUM_VERTICES;iv++){
    30913112                        gauss->GaussVertex(iv);
    30923113                        input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,(month+1)/12.*yts);
    … …  
    31063127
    31073128        /*measure the surface mass balance*/
    3108         for (int iv = 0; iv<numvertices; iv++){
     3129        for (int iv = 0; iv<NUM_VERTICES; iv++){
    31093130                smb[iv]=PddSurfaceMassBalanceSicopolis(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
    31103131                                        &melt[iv], &accu[iv], &melt_star[iv], &t_ampl[iv], &p_ampl[iv], yts, s[iv],
    … …  
    33243345                }
    33253346                case P1Enum:{
    3326                         int         numvertices = this->GetNumberOfVertices();
    3327                         IssmDouble *values      = xNew<IssmDouble>(numvertices);
    3328                         int        *connectivity= xNew<int>(numvertices);
    3329                         int        *sidlist     = xNew<int>(numvertices);
     3347                        const int NUM_VERTICES = this->GetNumberOfVertices();
     3348
     3349                        IssmDouble *values      = xNew<IssmDouble>(NUM_VERTICES);
     3350                        int        *connectivity= xNew<int>(NUM_VERTICES);
     3351                        int        *sidlist     = xNew<int>(NUM_VERTICES);
    33303352
    33313353                        this->GetVerticesSidList(sidlist);
    33323354                        this->GetVerticesConnectivityList(connectivity);
    33333355                        this->GetInputListOnVertices(values,output_enum);
    3334                         for(int i=0;i<numvertices;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
    3335 
    3336                         vector->SetValues(numvertices,sidlist,values,ADD_VAL);
     3356                        for(int i=0;i<NUM_VERTICES;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);
     3357
     3358                        vector->SetValues(NUM_VERTICES,sidlist,values,ADD_VAL);
    33373359
    33383360                        xDelete<IssmDouble>(values);
    … …  
    34213443        if (!IsOnSurface()) return;
    34223444
    3423         int  numvertices = this->GetNumberOfVertices();
    3424 
    3425         IssmDouble* s=xNew<IssmDouble>(numvertices);
    3426         IssmDouble* s0gcm=xNew<IssmDouble>(numvertices);
    3427         IssmDouble* st=xNew<IssmDouble>(numvertices);
     3445        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
     3446        const int NUM_VERTICES_DAYS_PER_YEAR    = NUM_VERTICES * 365;
     3447
     3448        IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
     3449        IssmDouble* s0gcm=xNew<IssmDouble>(NUM_VERTICES);
     3450        IssmDouble* st=xNew<IssmDouble>(NUM_VERTICES);
    34283451
    34293452        // daily forcing inputs
    3430         IssmDouble* dailyrainfall=xNew<IssmDouble>(365*numvertices);
    3431         IssmDouble* dailysnowfall=xNew<IssmDouble>(365*numvertices);
    3432         IssmDouble* dailydlradiation=xNew<IssmDouble>(365*numvertices);
    3433         IssmDouble* dailydsradiation=xNew<IssmDouble>(365*numvertices);
    3434         IssmDouble* dailywindspeed=xNew<IssmDouble>(365*numvertices);
    3435         IssmDouble* dailypressure=xNew<IssmDouble>(365*numvertices);
    3436         IssmDouble* dailyairdensity=xNew<IssmDouble>(365*numvertices);
    3437         IssmDouble* dailyairhumidity=xNew<IssmDouble>(365*numvertices);
    3438         IssmDouble* dailytemperature=xNew<IssmDouble>(365*numvertices);
     3453        IssmDouble* dailyrainfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3454        IssmDouble* dailysnowfall=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3455        IssmDouble* dailydlradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3456        IssmDouble* dailydsradiation=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3457        IssmDouble* dailywindspeed=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3458        IssmDouble* dailypressure=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3459        IssmDouble* dailyairdensity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3460        IssmDouble* dailyairhumidity=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
     3461        IssmDouble* dailytemperature=xNew<IssmDouble>(NUM_VERTICES_DAYS_PER_YEAR);
    34393462        // daily outputs
    3440         IssmDouble* tsurf_out=xNew<IssmDouble>(numvertices); memset(tsurf_out, 0., numvertices*sizeof(IssmDouble));
    3441         IssmDouble* smb_out=xNew<IssmDouble>(numvertices); memset(smb_out, 0., numvertices*sizeof(IssmDouble));
    3442         IssmDouble* saccu_out=xNew<IssmDouble>(numvertices); memset(saccu_out, 0., numvertices*sizeof(IssmDouble));
    3443         IssmDouble* smelt_out=xNew<IssmDouble>(numvertices); memset(smelt_out, 0., numvertices*sizeof(IssmDouble));
     3463        IssmDouble* tsurf_out=xNew<IssmDouble>(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     3464        IssmDouble* smb_out=xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     3465        IssmDouble* saccu_out=xNew<IssmDouble>(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     3466        IssmDouble* smelt_out=xNew<IssmDouble>(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble));
    34443467
    34453468        IssmDouble rho_water,rho_ice,desfac,rlaps,rdl;
    … …  
    34763499        Gauss* gauss=this->NewGauss();
    34773500        for (int iday = 0; iday < 365; iday++){
    3478                 for(int iv=0;iv<numvertices;iv++) {
     3501                for(int iv=0;iv<NUM_VERTICES;iv++) {
    34793502                        gauss->GaussVertex(iv);
    34803503                        /* get forcing */
    … …  
    35083531        }
    35093532
    3510         for (int iv = 0; iv<numvertices; iv++){
     3533        for (int iv = 0; iv<NUM_VERTICES; iv++){
    35113534                /* call semic */
    3512                 run_semic_(&dailysnowfall[iv*365], &dailyrainfall[iv*365], &dailydsradiation[iv*365], &dailydlradiation[iv*365], 
     3535                run_semic_(&dailysnowfall[iv*365], &dailyrainfall[iv*365], &dailydsradiation[iv*365], &dailydlradiation[iv*365],
    35133536                                        &dailywindspeed[iv*365], &dailypressure[iv*365], &dailyairdensity[iv*365], &dailyairhumidity[iv*365], &dailytemperature[iv*365],
    35143537                                        &tsurf_out[iv], &smb_out[iv], &saccu_out[iv], &smelt_out[iv]);
    … …  
    35163539
    35173540        switch(this->ObjectEnum()){
    3518                 case TriaEnum: 
     3541                case TriaEnum:
    35193542                        this->inputs->AddInput(new TriaInput(TemperatureSEMICEnum,&tsurf_out[0],P1Enum)); // TODO add TemperatureSEMICEnum to EnumDefinitions
    35203543                        this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb_out[0],P1Enum));
    … …  
    35253548                        // TODO
    35263549                        break;
    3527                 case TetraEnum: 
     3550                case TetraEnum:
    35283551                        // TODO
    35293552                        break;
    … …  
    35403563        xDelete<IssmDouble>(dailypressure);
    35413564        xDelete<IssmDouble>(dailyairdensity);
    3542         xDelete<IssmDouble>(dailyairhumidity); 
     3565        xDelete<IssmDouble>(dailyairhumidity);
    35433566        xDelete<IssmDouble>(dailypressure);
    35443567        xDelete<IssmDouble>(dailytemperature);
    … …  
    35483571        xDelete<IssmDouble>(tsurf_out);
    35493572        xDelete<IssmDouble>(s);
    3550         xDelete<IssmDouble>(st);       
     3573        xDelete<IssmDouble>(st);
    35513574        xDelete<IssmDouble>(s0gcm);
    35523575}
    … …  
    41614184
    41624185        /*Fetch number vertices and allocate memory*/
    4163         int         numvertices  = this->GetNumberOfVertices();
    4164         IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
     4186        const int NUM_VERTICES  = this->GetNumberOfVertices();
     4187
     4188        IssmDouble* maxprincipal = xNew<IssmDouble>(NUM_VERTICES);
    41654189
    41664190        /*Retrieve all inputs and parameters*/
    … …  
    41804204        /*loop over vertices: */
    41814205        Gauss* gauss=this->NewGauss();
    4182         for (int iv=0;iv<numvertices;iv++){
     4206        for (int iv=0;iv<NUM_VERTICES;iv++){
    41834207                gauss->GaussVertex(iv);
    41844208
    … …  
    45034527
    45044528        /*Fetch number vertices and allocate memory*/
    4505         int         numvertices    = this->GetNumberOfVertices();
    4506         IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
     4529        const int NUM_VERTICES = this->GetNumberOfVertices();
     4530
     4531        IssmDouble* viscousheating = xNew<IssmDouble>(NUM_VERTICES);
    45074532
    45084533        /*Retrieve all inputs and parameters*/
    … …  
    45144539        /*loop over vertices: */
    45154540        Gauss* gauss=this->NewGauss();
    4516         for (int iv=0;iv<numvertices;iv++){
     4541        for (int iv=0;iv<NUM_VERTICES;iv++){
    45174542                gauss->GaussVertex(iv);
    45184543
Note: See TracChangeset for help on using the changeset viewer.