Changeset 19274


Ignore:
Timestamp:
04/08/15 19:15:23 (10 years ago)
Author:
schlegel
Message:

CHG: allocate pdd matrices as xNew<IssmDouble>

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r19264 r19274  
    513513
    514514        int        i;
    515         IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
    516         IssmDouble TemperaturesPresentday[numvertices][12],TemperaturesLgm[numvertices][12];
    517         IssmDouble PrecipitationsPresentday[numvertices][12];
    518         IssmDouble tmp[numvertices];
     515        IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
     516        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
     517        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
     518        IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
     519        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
     520        IssmDouble* tmp=xNew<IssmDouble>(numvertices);
    519521        IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
    520522        IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
    … …  
    536538                for(int iv=0;iv<numvertices;iv++){
    537539                        gauss->GaussVertex(iv);
    538                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    539                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    540                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    541 
    542                         PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]*yts;
     540                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
     541                        input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts);
     542                        input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts);
     543
     544                        PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
    543545                }
    544546        }
    … …  
    556558                ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
    557559                                        Delta18oPresent, Delta18oLgm, Delta18oTime,
    558                                         &PrecipitationsPresentday[iv][0],
    559                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    560                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     560                                        &PrecipitationsPresentday[iv*12],
     561                                        &TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12],
     562                                        &monthlytemperatures[iv*12], &monthlyprec[iv*12]);
    561563        }
    562564
    … …  
    565567        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    566568        for (int imonth=0;imonth<12;imonth++) {
    567                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     569                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    568570                switch(this->ObjectEnum()){
    569571                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    572574                        default: _error_("Not implemented yet");
    573575                }
    574                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth]/yts;
     576                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    575577                switch(this->ObjectEnum()){
    576578                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    607609
    608610        int        i;
    609         IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
    610         IssmDouble TemperaturesPresentday[numvertices][12],TemperaturesLgm[numvertices][12];
    611         IssmDouble PrecipitationsPresentday[numvertices][12],PrecipitationsLgm[numvertices][12];
    612         IssmDouble tmp[numvertices];
     611        IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
     612        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
     613        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
     614        IssmDouble* TemperaturesLgm=xNew<IssmDouble>(12*numvertices);
     615        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
     616        IssmDouble* PrecipitationsLgm=xNew<IssmDouble>(12*numvertices);
     617        IssmDouble* tmp=xNew<IssmDouble>(numvertices);
    613618        IssmDouble TdiffTime,PfacTime;
    614619        IssmDouble time,yts,time_yr;
    … …  
    627632                for(int iv=0;iv<numvertices;iv++) {
    628633                        gauss->GaussVertex(iv);
    629                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    630                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    631                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    632                         input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
    633 
    634                         PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]*yts;
    635                         PrecipitationsLgm[iv][month]=PrecipitationsLgm[iv][month]*yts;
     634                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
     635                        input2->GetInputValue(&TemperaturesLgm[iv*12+month],gauss,month/12.*yts);
     636                        input3->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts);
     637                        input4->GetInputValue(&PrecipitationsLgm[iv*12+month],gauss,month/12.*yts);
     638
     639                        PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
     640                        PrecipitationsLgm[iv*12+month]=PrecipitationsLgm[iv*12+month]*yts;
    636641                }
    637642        }
    … …  
    644649        for(int iv=0;iv<numvertices;iv++){
    645650                ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
    646                                         &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
    647                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    648                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     651                                        &PrecipitationsLgm[iv*12],&PrecipitationsPresentday[iv*12],
     652                                        &TemperaturesLgm[iv*12], &TemperaturesPresentday[iv*12],
     653                                        &monthlytemperatures[iv*12], &monthlyprec[iv*12]);
    649654        }
    650655
    … …  
    653658        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    654659        for (int imonth=0;imonth<12;imonth++) {
    655                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     660                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    656661                switch(this->ObjectEnum()){
    657662                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    660665                        default: _error_("Not implemented yet");
    661666                }
    662                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth]/yts;
     667                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    663668                switch(this->ObjectEnum()){
    664669                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    695700
    696701        int        i;
    697         IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
    698         IssmDouble TemperaturesPresentday[numvertices][12];
    699         IssmDouble PrecipitationsPresentday[numvertices][12];
    700         IssmDouble tmp[numvertices];
     702        IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
     703        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
     704        IssmDouble* TemperaturesPresentday=xNew<IssmDouble>(12*numvertices);
     705        IssmDouble* PrecipitationsPresentday=xNew<IssmDouble>(12*numvertices);
     706        IssmDouble* tmp=xNew<IssmDouble>(numvertices);
    701707        IssmDouble Delta18oTime;
    702708        IssmDouble dpermil;
    … …  
    717723                for(int iv=0;iv<numvertices;iv++) {
    718724                        gauss->GaussVertex(iv);
    719                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    720                         input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    721 
    722                         PrecipitationsPresentday[iv][month]=PrecipitationsPresentday[iv][month]*yts;
     725                        input->GetInputValue(&TemperaturesPresentday[iv*12+month],gauss,month/12.*yts);
     726                        input2->GetInputValue(&PrecipitationsPresentday[iv*12+month],gauss,month/12.*yts);
     727
     728                        PrecipitationsPresentday[iv*12+month]=PrecipitationsPresentday[iv*12+month]*yts;
    723729                }
    724730        }
    … …  
    730736        for(int iv=0;iv<numvertices;iv++){
    731737                ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
    732                                         &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0],
    733                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     738                                        &PrecipitationsPresentday[iv*12], &TemperaturesPresentday[iv*12],
     739                                        &monthlytemperatures[iv*12], &monthlyprec[iv*12]);
    734740        }
    735741
    … …  
    738744        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    739745        for (int imonth=0;imonth<12;imonth++) {
    740                 for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     746                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    741747                switch(this->ObjectEnum()){
    742748                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    745751                        default: _error_("Not implemented yet");
    746752                }
    747                 for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth]/yts;
     753                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    748754                switch(this->ObjectEnum()){
    749755                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
    … …  
    16571663
    16581664        int        i;
    1659         IssmDouble agd[numvertices];             // surface mass balance
    1660         IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
    1661         IssmDouble yearlytemperatures[numvertices]; memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
    1662         IssmDouble tmp[numvertices];
    1663         IssmDouble h[numvertices],s[numvertices];
     1665        IssmDouble* agd=xNew<IssmDouble>(numvertices); // surface mass balance
     1666        IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
     1667        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
     1668        IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
     1669        IssmDouble* tmp=xNew<IssmDouble>(numvertices);
     1670        IssmDouble* h=xNew<IssmDouble>(numvertices);
     1671        IssmDouble* s=xNew<IssmDouble>(numvertices);
    16641672        IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    16651673        IssmDouble PfacTime,TdiffTime,sealevTime;
    … …  
    16901698                for(int iv=0;iv<numvertices;iv++) {
    16911699                        gauss->GaussVertex(iv);
    1692                         input->GetInputValue(&monthlytemperatures[iv][month],gauss,time_yr+month/12.*yts);
    1693                         yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin
    1694                         monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
    1695                         input2->GetInputValue(&monthlyprec[iv][month],gauss,time_yr+month/12.*yts);
    1696                         monthlyprec[iv][month]=monthlyprec[iv][month]*yts;
     1700                        input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,time_yr+month/12.*yts);
     1701                        yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv*12+month]*mavg; // Has to be in Kelvin
     1702                        monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
     1703                        input2->GetInputValue(&monthlyprec[iv*12+month],gauss,time_yr+month/12.*yts);
     1704                        monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
    16971705                }
    16981706        }
    … …  
    17161724        /*measure the surface mass balance*/
    17171725        for (int iv = 0; iv<numvertices; iv++){
    1718                 agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
     1726                agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
    17191727                                        pdds, pds, signorm, yts, h[iv], s[iv],
    17201728                                        desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
    … …  
    17261734        // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    17271735        // for (int imonth=0;imonth<12;imonth++) {
    1728         //   for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     1736        //   for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i*12+imonth];
    17291737        //   TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    17301738        //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    17311739        //
    1732         //   for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth];
     1740        //   for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i*12+imonth]/yts;
    17331741        //   TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    17341742        //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
Note: See TracChangeset for help on using the changeset viewer.