Changeset 19237


Ignore:
Timestamp:
04/01/15 18:54:42 (10 years ago)
Author:
schlegel
Message:

CHG: Allow pdd to run on temporal step finer than 1 year and move Pdd functions to Element

Location:
issm/trunk-jpl/src/c/classes/Elements
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r19215 r19237  
    503503        *pdmudB=dmudB;
    504504
     505}
     506/*}}}*/
     507void       Element::Delta18oParameterization(void){/*{{{*/
     508
     509        /*Are we on the base? If not, return*/
     510        if(!IsOnBase()) return;
     511
     512        int        numvertices = this->GetNumberOfVertices();
     513
     514        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];
     519        IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
     520        IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
     521        IssmDouble time,yts,finaltime,time_yr;
     522
     523        /*Recover parameters*/
     524        this->parameters->FindParam(&time,TimeEnum);
     525        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     526        this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
     527        time_yr=floor(time/yts)*yts;
     528
     529        /*Recover present day temperature and precipitation*/
     530        Input* input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
     531        Input* input2=this->inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
     532        Input* input3=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
     533        /*loop over vertices: */
     534        Gauss* gauss=this->NewGauss();
     535        for(int month=0;month<12;month++){
     536                for(int iv=0;iv<numvertices;iv++){
     537                        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        }
     543
     544        /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/
     545        this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime);
     546        this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-(21000*yts)));
     547        this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
     548        this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime);
     549        this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-(21000*yts)));
     550        this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time);
     551
     552        /*Compute the temperature and precipitation*/
     553        for(int iv=0;iv<numvertices;iv++){
     554                ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
     555                                        Delta18oPresent, Delta18oLgm, Delta18oTime,
     556                                        &PrecipitationsPresentday[iv][0],
     557                                        &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
     558                                        &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     559        }
     560
     561        /*Update inputs*/
     562        TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
     563        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
     564        for (int imonth=0;imonth<12;imonth++) {
     565                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     566                switch(this->ObjectEnum()){
     567                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     568                        case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     569                        case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     570                        default: _error_("Not implemented yet");
     571                }
     572                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth];
     573                switch(this->ObjectEnum()){
     574                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     575                        case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     576                        case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     577                        default: _error_("Not implemented yet");
     578                }
     579        }
     580        NewTemperatureInput->Configure(this->parameters);
     581        NewPrecipitationInput->Configure(this->parameters);
     582
     583        this->inputs->AddInput(NewTemperatureInput);
     584        this->inputs->AddInput(NewPrecipitationInput);
     585
     586        switch(this->ObjectEnum()){
     587                case TriaEnum: break;
     588                case PentaEnum:
     589                case TetraEnum:
     590                                                        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     591                                                        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     592                                                        break;
     593                default: _error_("Not implemented yet");
     594        }
     595
     596        /*clean-up*/
     597        delete gauss;
     598}
     599/*}}}*/
     600void       Element::MungsmtpParameterization(void){/*{{{*/
     601        /*Are we on the base? If not, return*/
     602        if(!IsOnBase()) return;
     603
     604        int        numvertices = this->GetNumberOfVertices();
     605
     606        int        i;
     607        IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
     608        IssmDouble TemperaturesPresentday[numvertices][12],TemperaturesLgm[numvertices][12];
     609        IssmDouble PrecipitationsPresentday[numvertices][12],PrecipitationsLgm[numvertices][12];
     610        IssmDouble tmp[numvertices];
     611        IssmDouble TdiffTime,PfacTime;
     612        IssmDouble time,yts,time_yr;
     613        this->parameters->FindParam(&time,TimeEnum);
     614        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     615        time_yr=floor(time/yts)*yts;
     616
     617        /*Recover present day temperature and precipitation*/
     618        Input*     input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
     619        Input*     input2=this->inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
     620        Input*     input3=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
     621        Input*     input4=this->inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);        _assert_(input4);
     622        /*loop over vertices: */
     623        Gauss* gauss=this->NewGauss();
     624        for(int month=0;month<12;month++) {
     625                for(int iv=0;iv<numvertices;iv++) {
     626                        gauss->GaussVertex(iv);
     627                        input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
     628                        input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
     629                        input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
     630                        input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
     631                }
     632        }
     633
     634        /*Recover interpolation parameters at time t*/
     635        this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
     636        this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
     637
     638        /*Compute the temperature and precipitation*/
     639        for(int iv=0;iv<numvertices;iv++){
     640                ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
     641                                        &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
     642                                        &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
     643                                        &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     644        }
     645
     646        /*Update inputs*/
     647        TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
     648        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
     649        for (int imonth=0;imonth<12;imonth++) {
     650                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     651                switch(this->ObjectEnum()){
     652                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     653                        case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     654                        case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     655                        default: _error_("Not implemented yet");
     656                }
     657                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth];
     658                switch(this->ObjectEnum()){
     659                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     660                        case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     661                        case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     662                        default: _error_("Not implemented yet");
     663                }
     664        }
     665        NewTemperatureInput->Configure(this->parameters);
     666        NewPrecipitationInput->Configure(this->parameters);
     667
     668        this->inputs->AddInput(NewTemperatureInput);
     669        this->inputs->AddInput(NewPrecipitationInput);
     670
     671        switch(this->ObjectEnum()){
     672                case TriaEnum: break;
     673                case PentaEnum:
     674                case TetraEnum:
     675                                                        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     676                                                        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     677                                                        break;
     678                default: _error_("Not implemented yet");
     679        }
     680
     681        /*clean-up*/
     682        delete gauss;
     683}
     684/*}}}*/
     685void       Element::Delta18opdParameterization(void){/*{{{*/
     686        /*Are we on the base? If not, return*/
     687        if(!IsOnBase()) return;
     688
     689        int        numvertices = this->GetNumberOfVertices();
     690
     691        int        i;
     692        IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
     693        IssmDouble TemperaturesPresentday[numvertices][12];
     694        IssmDouble PrecipitationsPresentday[numvertices][12];
     695        IssmDouble tmp[numvertices];
     696        IssmDouble Delta18oTime;
     697        IssmDouble dpermil;
     698        IssmDouble time,yts,time_yr,month;
     699        this->parameters->FindParam(&time,TimeEnum);
     700        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     701        time_yr=floor(time/yts)*yts;
     702
     703        /*Get some pdd parameters*/
     704        dpermil=this->matpar->GetMaterialParameter(SurfaceforcingsDpermilEnum);
     705
     706        /*Recover present day temperature and precipitation*/
     707        Input*     input=this->inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
     708        Input*     input2=this->inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
     709        /*loop over vertices: */
     710        Gauss* gauss=this->NewGauss();
     711        for(int month=0;month<12;month++) {
     712                for(int iv=0;iv<numvertices;iv++) {
     713                        gauss->GaussVertex(iv);
     714                        input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
     715                        input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
     716                }
     717        }
     718
     719        /*Recover interpolation parameters at time t*/
     720        this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
     721
     722        /*Compute the temperature and precipitation*/
     723        for(int iv=0;iv<numvertices;iv++){
     724                ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
     725                                        &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0],
     726                                        &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
     727        }
     728
     729        /*Update inputs*/
     730        TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
     731        TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
     732        for (int imonth=0;imonth<12;imonth++) {
     733                for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     734                switch(this->ObjectEnum()){
     735                        case TriaEnum:  NewTemperatureInput->AddTimeInput(new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     736                        case PentaEnum: NewTemperatureInput->AddTimeInput(new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     737                        case TetraEnum: NewTemperatureInput->AddTimeInput(new TetraInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     738                        default: _error_("Not implemented yet");
     739                }
     740                for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth];
     741                switch(this->ObjectEnum()){
     742                        case TriaEnum:  NewPrecipitationInput->AddTimeInput(new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     743                        case PentaEnum: NewPrecipitationInput->AddTimeInput(new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     744                        case TetraEnum: NewPrecipitationInput->AddTimeInput(new TetraInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum),time_yr+imonth/12.*yts); break;
     745                        default: _error_("Not implemented yet");
     746                }
     747        }
     748        NewTemperatureInput->Configure(this->parameters);
     749        NewPrecipitationInput->Configure(this->parameters);
     750
     751        this->inputs->AddInput(NewTemperatureInput);
     752        this->inputs->AddInput(NewPrecipitationInput);
     753
     754        switch(this->ObjectEnum()){
     755                case TriaEnum: break;
     756                case PentaEnum:
     757                case TetraEnum:
     758                                                        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     759                                                        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     760                                                        break;
     761                default: _error_("Not implemented yet");
     762        }
     763
     764        /*clean-up*/
     765        delete gauss;
    505766}
    506767/*}}}*/
    … …  
    13781639ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
    13791640        return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
     1641}
     1642/*}}}*/
     1643void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
     1644
     1645        int  numvertices = this->GetNumberOfVertices();
     1646
     1647        int        i;
     1648        IssmDouble agd[numvertices];             // surface mass balance
     1649        IssmDouble monthlytemperatures[numvertices][12],monthlyprec[numvertices][12];
     1650        IssmDouble yearlytemperatures[numvertices]; memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
     1651        IssmDouble tmp[numvertices];
     1652        IssmDouble h[numvertices],s[numvertices];
     1653        IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
     1654        IssmDouble PfacTime,TdiffTime,sealevTime;
     1655        IssmDouble mavg=1./12.; //factor for monthly average
     1656
     1657        /*Get material parameters :*/
     1658        rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     1659        rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     1660
     1661        /*Get some pdd parameters*/
     1662        desfac=this->matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
     1663        s0p=this->matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
     1664        s0t=this->matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
     1665        rlaps=this->matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
     1666        rlapslgm=this->matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
     1667
     1668        /*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/
     1669        Input*     input=this->inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
     1670        Input*     input2=this->inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
     1671        IssmDouble time,yts,time_yr;
     1672        this->parameters->FindParam(&time,TimeEnum);
     1673        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     1674        time_yr=floor(time/yts)*yts;
     1675
     1676        /*loop over vertices: */
     1677        Gauss* gauss=this->NewGauss();
     1678        for(int month=0;month<12;month++) {
     1679                for(int iv=0;iv<numvertices;iv++) {
     1680                        gauss->GaussVertex(iv);
     1681                        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time_yr+month/12.*yts);
     1682                        yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin
     1683                        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
     1684                        input2->GetInputValue(&monthlyprec[iv][month],gauss,time_yr+month/12.*yts);
     1685                }
     1686        }
     1687
     1688        /*Recover Pfac, Tdiff and sealev at time t:
     1689         *     This parameters are used to interpolate the temperature
     1690         *         and precipitaton between PD and LGM when ismungsm==1 */
     1691        if (ismungsm==1){
     1692                this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
     1693                this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
     1694        }
     1695        else {
     1696                TdiffTime=0;
     1697                sealevTime=0;
     1698        }
     1699
     1700        /*Recover info at the vertices: */
     1701        GetInputListOnVertices(&h[0],ThicknessEnum);
     1702        GetInputListOnVertices(&s[0],SurfaceEnum);
     1703
     1704        /*measure the surface mass balance*/
     1705        for (int iv = 0; iv<numvertices; iv++){
     1706                agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
     1707                                        pdds, pds, signorm, yts, h[iv], s[iv],
     1708                                        desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
     1709                                        rho_water,rho_ice);
     1710        }
     1711
     1712        /*Update inputs*/
     1713        // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
     1714        // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
     1715        // for (int imonth=0;imonth<12;imonth++) {
     1716        //   for(i=0;i<numvertices;i++) tmp[i]=monthlytemperatures[i][imonth];
     1717        //   TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
     1718        //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
     1719        //
     1720        //   for(i=0;i<numvertices;i++) tmp[i]=monthlyprec[i][imonth];
     1721        //   TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
     1722        //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
     1723        // }
     1724        // NewTemperatureInput->Configure(this->parameters);
     1725        // NewPrecipitationInput->Configure(this->parameters);
     1726
     1727        switch(this->ObjectEnum()){
     1728                case TriaEnum: 
     1729                        this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     1730                        this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
     1731                        break;
     1732                case PentaEnum:
     1733                        this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     1734                        this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
     1735                        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     1736                        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     1737                        break;
     1738                case TetraEnum:
     1739                        this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     1740                        this->inputs->AddInput(new TetraInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
     1741                        this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     1742                        this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     1743                        break;
     1744                default: _error_("Not implemented yet");
     1745        }
     1746        // this->inputs->AddInput(NewTemperatureInput);
     1747        // this->inputs->AddInput(NewPrecipitationInput);
     1748        // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
     1749
     1750        //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
     1751        // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
     1752        // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
     1753
     1754        /*clean-up*/
     1755        delete gauss;
    13801756}
    13811757/*}}}*/

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

    r19215 r19237  
    6666                void               DeleteInput(int input_enum);
    6767                void               DeleteMaterials(void);
     68                void               Delta18oParameterization(void);
     69                void               MungsmtpParameterization(void);
     70                void               Delta18opdParameterization(void);
    6871                IssmDouble         Divergence(void);
    6972                void               dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    … …  
    119122                void               LinearFloatingiceMeltingRate();
    120123                void               MigrateGroundingLine(IssmDouble* sheet_ungrounding);
    121                 void Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses);
    122 
     124                void               Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction,int numanalyses);
    123125                ElementMatrix*     NewElementMatrix(int approximation_enum=NoneApproximationEnum);
    124126                ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
    125127                ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
     128                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
    126129                IssmDouble         PureIceEnthalpy(IssmDouble pressure);
    127130                void               ResetIceLevelset(void);
    … …  
    183186                virtual void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index)=0;
    184187                virtual void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum)=0;
    185                 virtual void       Delta18oParameterization(void)=0;
    186                 virtual void       MungsmtpParameterization(void)=0;
    187                 virtual void       Delta18opdParameterization(void)=0;
    188188                virtual void       ElementResponse(IssmDouble* presponse,int response_enum)=0;
    189189                virtual void       ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0;
    … …  
    258258                virtual int        NumberofNodesPressure(void)=0;
    259259                virtual int        NumberofNodesVelocity(void)=0;
    260                 virtual void       PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm)=0;
    261260                virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
    262261                virtual int        PressureInterpolation()=0;

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

    r19215 r19237  
    7777
    7878        int i;
    79         bool* hnodes_null=NULL; /*intermediary needed*/
     79        bool* hnodesi_null=NULL; /*intermediary needed*/
     80        bool  hnodes_null=true; /*this could be NULL on empty constructor*/
    8081        bool  hneighbors_null=true; /*don't deal with hneighbors, unless explicitely asked to*/
    8182
    … …  
    8586        if(marshall_direction==MARSHALLING_FORWARD || marshall_direction==MARSHALLING_SIZE){
    8687                if(this->hneighbors)hneighbors_null=false;
    87                 hnodes_null=xNew<bool>(numanalyses);
    88                 for(i=0;i<numanalyses;i++){
    89                         if(this->hnodes[i])hnodes_null[i]=false;
    90                         else hnodes_null[i]=true;
     88                if(this->hnodes){
     89                        hnodes_null=false;
     90                        hnodesi_null=xNew<bool>(numanalyses);
     91                        for(i=0;i<numanalyses;i++){
     92                                if(this->hnodes[i])hnodesi_null[i]=false;
     93                                else hnodesi_null[i]=true;
     94                        }
    9195                }
    9296        }
    … …  
    96100        MARSHALLING(numanalyses);
    97101        MARSHALLING(hneighbors_null);
    98         MARSHALLING_DYNAMIC(hnodes_null,bool,numanalyses);
     102        MARSHALLING(hnodes_null);
     103        MARSHALLING_DYNAMIC(hnodesi_null,bool,numanalyses);
    99104
    100105        if(marshall_direction==MARSHALLING_BACKWARD){
    101106               
    102                 this->hnodes      = new Hook*[numanalyses];
     107                if (!hnodes_null)this->hnodes = new Hook*[numanalyses];
     108                else this->hnodes=NULL;
    103109                this->hvertices   = new Hook();
    104110                this->hmaterial   = new Hook();
    105111                this->hmatpar     = new Hook();
    106112                if(!hneighbors_null)this->hneighbors  = new Hook();
     113                else this->hneighbors=NULL;
    107114
    108115                /*Initialize hnodes: */
    109                 for(int i=0;i<this->numanalyses;i++){
    110                         if(!hnodes_null[i])this->hnodes[i]=new Hook();
    111                         else this->hnodes[i]=NULL;
     116                if (this->hnodes){
     117                        for(int i=0;i<this->numanalyses;i++){
     118                                if(!hnodesi_null[i])this->hnodes[i]=new Hook();
     119                                else this->hnodes[i]=NULL;
     120                        }
    112121                }
    113122        }
    114123
    115         for (i=0;i<numanalyses;i++) if(this->hnodes[i])this->hnodes[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     124        if (this->hnodes){
     125                for (i=0;i<numanalyses;i++) if(this->hnodes[i])this->hnodes[i]->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     126        }
    116127        this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    117128        this->hmatpar->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    … …  
    120131
    121132        /*Free ressources: */
    122         xDelete<bool>(hnodes_null);
     133        xDelete<bool>(hnodesi_null);
    123134
    124135}

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

    r19172 r19237  
    580580
    581581}/*}}}*/
    582 void       Penta::Delta18oParameterization(void){/*{{{*/
    583 
    584         /*Are we on the base? If not, return*/
    585         if(!IsOnBase()) return;
    586 
    587         int        i;
    588         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    589         IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
    590         IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
    591         IssmDouble tmp[NUMVERTICES];
    592         IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
    593         IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
    594         IssmDouble time,yts,finaltime;
    595         this->parameters->FindParam(&time,TimeEnum);
    596         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    597         this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    598 
    599         /*Recover present day temperature and precipitation*/
    600         Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    601         Input*     input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
    602         Input*     input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
    603         GaussPenta* gauss=new GaussPenta();
    604         for(int month=0;month<12;month++) {
    605                 for(int iv=0;iv<NUMVERTICES;iv++) {
    606                         gauss->GaussVertex(iv);
    607                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    608                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    609                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    610                 }
    611         }
    612 
    613         /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/
    614         this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime);
    615         this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,finaltime-(21000*yts));
    616         this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
    617         this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime);
    618         this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,finaltime-(21000*yts));
    619         this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time);
    620 
    621         /*Compute the temperature and precipitation*/
    622         for(int iv=0;iv<NUMVERTICES;iv++){
    623                 ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
    624                                         Delta18oPresent, Delta18oLgm, Delta18oTime,
    625                                         &PrecipitationsPresentday[iv][0],
    626                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    627                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    628         }
    629 
    630         /*Update inputs*/
    631         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    632         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    633         for (int imonth=0;imonth<12;imonth++) {
    634                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    635                 PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    636                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    637 
    638                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    639                 PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    640                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    641         }
    642         NewTemperatureInput->Configure(this->parameters);
    643         NewPrecipitationInput->Configure(this->parameters);
    644 
    645         this->inputs->AddInput(NewTemperatureInput);
    646         this->inputs->AddInput(NewPrecipitationInput);
    647 
    648         this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    649         this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
    650 
    651         /*clean-up*/
    652         delete gauss;
    653 }
    654 /*}}}*/
    655 void       Penta::MungsmtpParameterization(void){/*{{{*/
    656         /*Are we on the base? If not, return*/
    657         if(!IsOnBase()) return;
    658 
    659         int        i;
    660         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    661         IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
    662         IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
    663         IssmDouble tmp[NUMVERTICES];
    664         IssmDouble TdiffTime,PfacTime;
    665         IssmDouble time,yts;
    666         this->parameters->FindParam(&time,TimeEnum);
    667         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    668 
    669         /*Recover present day temperature and precipitation*/
    670         Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    671         Input*     input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
    672         Input*     input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
    673         Input*     input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);        _assert_(input4);
    674         GaussPenta* gauss=new GaussPenta();
    675         for(int month=0;month<12;month++) {
    676                 for(int iv=0;iv<NUMVERTICES;iv++) {
    677                         gauss->GaussVertex(iv);
    678                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    679                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    680                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    681                         input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
    682                 }
    683         }
    684 
    685         /*Recover interpolation parameters at time t*/
    686         this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
    687         this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
    688 
    689         /*Compute the temperature and precipitation*/
    690         for(int iv=0;iv<NUMVERTICES;iv++){
    691           ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
    692                                         &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
    693                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    694                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    695         }
    696 
    697         /*Update inputs*/
    698         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    699         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    700         for (int imonth=0;imonth<12;imonth++) {
    701                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    702                 PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    703                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    704        
    705                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    706                 PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    707                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    708         }
    709         NewTemperatureInput->Configure(this->parameters);
    710         NewPrecipitationInput->Configure(this->parameters);
    711 
    712         this->inputs->AddInput(NewTemperatureInput);
    713         this->inputs->AddInput(NewPrecipitationInput);
    714        
    715         this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    716         this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
    717 
    718         /*clean-up*/
    719         delete gauss;
    720 }
    721 /*}}}*/
    722 void       Penta::Delta18opdParameterization(void){/*{{{*/
    723         /*Are we on the base? If not, return*/
    724         if(!IsOnBase()) return;
    725 
    726         int        i;
    727         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    728         IssmDouble TemperaturesPresentday[NUMVERTICES][12];
    729         IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
    730         IssmDouble tmp[NUMVERTICES];
    731         IssmDouble Delta18oTime;
    732         IssmDouble dpermil;
    733         IssmDouble time,yts;
    734         this->parameters->FindParam(&time,TimeEnum);
    735         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    736 
    737         /*Get some pdd parameters*/
    738         dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum);
    739 
    740         /*Recover present day temperature and precipitation*/
    741         Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    742         Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
    743         GaussPenta* gauss=new GaussPenta();
    744         for(int month=0;month<12;month++) {
    745                 for(int iv=0;iv<NUMVERTICES;iv++) {
    746                         gauss->GaussVertex(iv);
    747                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    748                         input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    749                 }
    750         }
    751 
    752         /*Recover interpolation parameters at time t*/
    753         this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
    754 
    755         /*Compute the temperature and precipitation*/
    756         for(int iv=0;iv<NUMVERTICES;iv++){
    757           ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
    758                                         &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0],
    759                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    760         }
    761 
    762         /*Update inputs*/
    763         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    764         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    765         for (int imonth=0;imonth<12;imonth++) {
    766                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    767                 PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    768                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    769 
    770                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    771                 PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    772                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    773         }
    774         NewTemperatureInput->Configure(this->parameters);
    775         NewPrecipitationInput->Configure(this->parameters);
    776 
    777         this->inputs->AddInput(NewTemperatureInput);
    778         this->inputs->AddInput(NewPrecipitationInput);
    779 
    780         this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    781         this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
    782 
    783         /*clean-up*/
    784         delete gauss;
    785 }
    786 /*}}}*/
    787582void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
    788583
    … …  
    22242019        return PentaEnum;
    22252020
    2226 }
    2227 /*}}}*/
    2228 void       Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
    2229 
    2230    int        i;
    2231    IssmDouble agd[NUMVERTICES];             // surface mass balance
    2232    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    2233    IssmDouble yearlytemperatures[NUMVERTICES]={0.};
    2234    IssmDouble tmp[NUMVERTICES];
    2235    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    2236    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    2237    IssmDouble PfacTime,TdiffTime,sealevTime;
    2238    IssmDouble mavg=1./12.; //factor for monthly average
    2239 
    2240    /*Get material parameters :*/
    2241    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    2242    rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    2243 
    2244   /*Get some pdd parameters*/
    2245   desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
    2246   s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
    2247   s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
    2248   rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
    2249   rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
    2250 
    2251    /*Recover monthly temperatures and precipitation*/
    2252    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    2253    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
    2254    GaussPenta* gauss=new GaussPenta();
    2255    IssmDouble time,yts;
    2256    this->parameters->FindParam(&time,TimeEnum);
    2257    this->parameters->FindParam(&yts,ConstantsYtsEnum);
    2258    
    2259    for(int month=0;month<12;month++) {
    2260      for(int iv=0;iv<NUMVERTICES;iv++) {
    2261        gauss->GaussVertex(iv);
    2262        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
    2263        yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin
    2264        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
    2265        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
    2266      }
    2267    }
    2268 
    2269   /*Recover Pfac, Tdiff and sealev at time t:
    2270     This parameters are used to interpolate the temperature
    2271     and precipitaton between PD and LGM when ismungsm==1 */
    2272   if (ismungsm==1){ 
    2273     this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
    2274     this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
    2275   }
    2276   else {
    2277     TdiffTime=0;
    2278     sealevTime=0; 
    2279   }
    2280 
    2281   /*Recover info at the vertices: */
    2282   GetInputListOnVertices(&h[0],ThicknessEnum);
    2283   GetInputListOnVertices(&s[0],SurfaceEnum);
    2284 
    2285 
    2286    /*measure the surface mass balance*/
    2287    for (int iv = 0; iv < NUMVERTICES; iv++){
    2288      agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], 
    2289                                   pdds,pds, signorm, yts, h[iv], s[iv],
    2290                                   desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
    2291                                   rho_water,rho_ice);
    2292    }
    2293 
    2294    /*Update inputs*/   
    2295    // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    2296    // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    2297    // for (int imonth=0;imonth<12;imonth++) {
    2298    //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    2299    //   PentaInput* newmonthinput1 = new PentaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    2300    //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    2301    //
    2302    //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    2303    //   PentaInput* newmonthinput2 = new PentaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    2304    //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    2305    // }
    2306    // NewTemperatureInput->Configure(this->parameters);
    2307    // NewPrecipitationInput->Configure(this->parameters);
    2308 
    2309 
    2310 
    2311    this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
    2312    this->inputs->AddInput(new PentaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    2313    // this->inputs->AddInput(NewTemperatureInput);
    2314    // this->inputs->AddInput(NewPrecipitationInput);
    2315    // //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2316     this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
    2317     this->InputExtrude(TemperatureEnum,-1);
    2318    // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    2319    // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
    2320 
    2321    /*clean-up*/
    2322    delete gauss;
    23232021}
    23242022/*}}}*/

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

    r19215 r19237  
    5959                void           ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    6060                ElementMatrix* CreateBasalMassMatrix(void);
    61                 void           Delta18oParameterization(void);
    62                 void           MungsmtpParameterization(void);
    63                 void           Delta18opdParameterization(void);
    6461                void           ElementResponse(IssmDouble* presponse,int response_enum);
    6562                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
    … …  
    141138                int            NumberofNodesPressure(void);
    142139                int            NumberofNodesVelocity(void);
    143                 void           PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
    144140                void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    145141                int            PressureInterpolation();

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

    r19215 r19237  
    5454                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    5555                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    56                 void        Delta18oParameterization(void){_error_("not implemented yet");};
    57                 void        MungsmtpParameterization(void){_error_("not implemented yet");};
    58                 void        Delta18opdParameterization(void){_error_("not implemented yet");};
    5956                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    6057                void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");};
    … …  
    133130                int         NumberofNodesPressure(void){_error_("not implemented yet");};
    134131                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    135                 void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
    136132                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    137133                int         PressureInterpolation(void){_error_("not implemented yet");};

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

    r19215 r19237  
    5454                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){_error_("not implemented yet");};
    5555                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){_error_("not implemented yet");};
    56                 void        Delta18oParameterization(void){_error_("not implemented yet");};
    57                 void        MungsmtpParameterization(void){_error_("not implemented yet");};
    58                 void        Delta18opdParameterization(void){_error_("not implemented yet");};
    5956                IssmDouble  DragCoefficientAbsGradient(void){_error_("not implemented yet");};
    6057                void        ElementResponse(IssmDouble* presponse,int response_enum){_error_("not implemented yet");};
    … …  
    139136                int         NumberofNodesPressure(void);
    140137                int         NumberofNodesVelocity(void);
    141                 void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){_error_("not implemented yet");};
    142138                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    143139                int         PressureInterpolation(void);

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

    r19215 r19237  
    123123        Element::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses);
    124124        TriaRef::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
     125
     126        vertices = (Vertex**)this->hvertices->deliverp();
     127        material = (Material*)this->hmaterial->delivers();
     128        matpar   = (Matpar*)this->hmatpar->delivers();
    125129
    126130}
    … …  
    654658
    655659}/*}}}*/
    656 void       Tria::Delta18oParameterization(void){/*{{{*/
    657 
    658         int        i;
    659         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    660         IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
    661         IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
    662         IssmDouble tmp[NUMVERTICES];
    663         IssmDouble Delta18oPresent,Delta18oLgm,Delta18oTime;
    664         IssmDouble Delta18oSurfacePresent,Delta18oSurfaceLgm,Delta18oSurfaceTime;
    665         IssmDouble time,yts,finaltime;
    666 
    667         /*Recover parameters*/
    668         this->parameters->FindParam(&time,TimeEnum);
    669         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    670         this->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    671 
    672         /*Recover present day temperature and precipitation*/
    673         Input* input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    674         Input* input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
    675         Input* input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
    676         GaussTria* gauss=new GaussTria();
    677         for(int month=0;month<12;month++){
    678                 for(int iv=0;iv<NUMVERTICES;iv++){
    679                         gauss->GaussVertex(iv);
    680                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    681                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    682                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    683                 }
    684         }
    685 
    686         /*Recover delta18o and Delta18oSurface at present day, lgm and at time t*/
    687         this->parameters->FindParam(&Delta18oPresent,SurfaceforcingsDelta18oEnum,finaltime);
    688         this->parameters->FindParam(&Delta18oLgm,SurfaceforcingsDelta18oEnum,(finaltime-(21000*yts)));
    689         this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
    690         this->parameters->FindParam(&Delta18oSurfacePresent,SurfaceforcingsDelta18oSurfaceEnum,finaltime);
    691         this->parameters->FindParam(&Delta18oSurfaceLgm,SurfaceforcingsDelta18oSurfaceEnum,(finaltime-(21000*yts)));
    692         this->parameters->FindParam(&Delta18oSurfaceTime,SurfaceforcingsDelta18oSurfaceEnum,time);
    693 
    694         /*Compute the temperature and precipitation*/
    695         for(int iv=0;iv<NUMVERTICES;iv++){
    696                 ComputeDelta18oTemperaturePrecipitation(Delta18oSurfacePresent, Delta18oSurfaceLgm, Delta18oSurfaceTime,
    697                                         Delta18oPresent, Delta18oLgm, Delta18oTime,
    698                                         &PrecipitationsPresentday[iv][0],
    699                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    700                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    701         }
    702 
    703         /*Update inputs*/
    704         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    705         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    706         for (int imonth=0;imonth<12;imonth++) {
    707                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    708                 TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    709                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    710 
    711                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    712                 TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    713                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    714         }
    715         NewTemperatureInput->Configure(this->parameters);
    716         NewPrecipitationInput->Configure(this->parameters);
    717 
    718         this->inputs->AddInput(NewTemperatureInput);
    719         this->inputs->AddInput(NewPrecipitationInput);
    720 
    721         /*clean-up*/
    722         delete gauss;
    723 }
    724 /*}}}*/
    725 void       Tria::MungsmtpParameterization(void){/*{{{*/
    726         /*Are we on the base? If not, return*/
    727         if(!IsOnBase()) return;
    728 
    729         int        i;
    730         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    731         IssmDouble TemperaturesPresentday[NUMVERTICES][12],TemperaturesLgm[NUMVERTICES][12];
    732         IssmDouble PrecipitationsPresentday[NUMVERTICES][12],PrecipitationsLgm[NUMVERTICES][12];
    733         IssmDouble tmp[NUMVERTICES];
    734         IssmDouble TdiffTime,PfacTime;
    735         IssmDouble time,yts;
    736         this->parameters->FindParam(&time,TimeEnum);
    737         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    738 
    739         /*Recover present day temperature and precipitation*/
    740         Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    741         Input*     input2=inputs->GetInput(SurfaceforcingsTemperaturesLgmEnum);          _assert_(input2);
    742         Input*     input3=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input3);
    743         Input*     input4=inputs->GetInput(SurfaceforcingsPrecipitationsLgmEnum);        _assert_(input4);
    744         GaussTria* gauss=new GaussTria();
    745         for(int month=0;month<12;month++) {
    746                 for(int iv=0;iv<NUMVERTICES;iv++) {
    747                         gauss->GaussVertex(iv);
    748                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    749                         input2->GetInputValue(&TemperaturesLgm[iv][month],gauss,month/12.*yts);
    750                         input3->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    751                         input4->GetInputValue(&PrecipitationsLgm[iv][month],gauss,month/12.*yts);
    752                 }
    753         }
    754 
    755         /*Recover interpolation parameters at time t*/
    756         this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
    757         this->parameters->FindParam(&PfacTime,SurfaceforcingsPfacEnum,time);
    758 
    759         /*Compute the temperature and precipitation*/
    760         for(int iv=0;iv<NUMVERTICES;iv++){
    761           ComputeMungsmTemperaturePrecipitation(TdiffTime,PfacTime,
    762                                         &PrecipitationsLgm[iv][0],&PrecipitationsPresentday[iv][0],
    763                                         &TemperaturesLgm[iv][0], &TemperaturesPresentday[iv][0],
    764                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    765         }
    766 
    767         /*Update inputs*/
    768         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    769         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    770         for (int imonth=0;imonth<12;imonth++) {
    771                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    772                 TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    773                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    774 
    775                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    776                 TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    777                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    778         }
    779         NewTemperatureInput->Configure(this->parameters);
    780         NewPrecipitationInput->Configure(this->parameters);
    781 
    782         this->inputs->AddInput(NewTemperatureInput);
    783         this->inputs->AddInput(NewPrecipitationInput);
    784 
    785         /*clean-up*/
    786         delete gauss;
    787 }
    788 /*}}}*/
    789 void       Tria::Delta18opdParameterization(void){/*{{{*/
    790         /*Are we on the base? If not, return*/
    791         if(!IsOnBase()) return;
    792 
    793         int        i;
    794         IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    795         IssmDouble TemperaturesPresentday[NUMVERTICES][12];
    796         IssmDouble PrecipitationsPresentday[NUMVERTICES][12];
    797         IssmDouble tmp[NUMVERTICES];
    798         IssmDouble Delta18oTime;
    799         IssmDouble dpermil;
    800         IssmDouble time,yts;
    801         this->parameters->FindParam(&time,TimeEnum);
    802         this->parameters->FindParam(&yts,ConstantsYtsEnum);
    803 
    804         /*Get some pdd parameters*/
    805         dpermil=matpar-> GetMaterialParameter(SurfaceforcingsDpermilEnum);
    806                        
    807         /*Recover present day temperature and precipitation*/
    808         Input*     input=inputs->GetInput(SurfaceforcingsTemperaturesPresentdayEnum);    _assert_(input);
    809         Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationsPresentdayEnum); _assert_(input2);
    810         GaussTria* gauss=new GaussTria();
    811         for(int month=0;month<12;month++) {
    812                 for(int iv=0;iv<NUMVERTICES;iv++) {
    813                         gauss->GaussVertex(iv);
    814                         input->GetInputValue(&TemperaturesPresentday[iv][month],gauss,month/12.*yts);
    815                         input2->GetInputValue(&PrecipitationsPresentday[iv][month],gauss,month/12.*yts);
    816                 }
    817         }
    818 
    819         /*Recover interpolation parameters at time t*/
    820         this->parameters->FindParam(&Delta18oTime,SurfaceforcingsDelta18oEnum,time);
    821 
    822         /*Compute the temperature and precipitation*/
    823         for(int iv=0;iv<NUMVERTICES;iv++){
    824           ComputeD18OTemperaturePrecipitationFromPD(Delta18oTime,dpermil,
    825                                         &PrecipitationsPresentday[iv][0], &TemperaturesPresentday[iv][0],
    826                                         &monthlytemperatures[iv][0], &monthlyprec[iv][0]);
    827         }
    828 
    829         /*Update inputs*/
    830         TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    831         TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    832         for (int imonth=0;imonth<12;imonth++) {
    833                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    834                 TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    835                 NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    836 
    837                 for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    838                 TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    839                 NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    840         }
    841         NewTemperatureInput->Configure(this->parameters);
    842         NewPrecipitationInput->Configure(this->parameters);
    843 
    844         this->inputs->AddInput(NewTemperatureInput);
    845         this->inputs->AddInput(NewPrecipitationInput);
    846 
    847         /*clean-up*/
    848         delete gauss;
    849 }
    850 /*}}}*/
    851660int        Tria::EdgeOnBaseIndex(void){/*{{{*/
    852661
    … …  
    26382447}
    26392448/*}}}*/
    2640 int        Tria::PressureInterpolation(void){/*{{{*/
    2641         return TriaRef::PressureInterpolation(this->element_type);
    2642 }
    2643 /*}}}*/
    26442449int        Tria::NumberofNodesPressure(void){/*{{{*/
    26452450        return TriaRef::NumberofNodes(this->PressureInterpolation());
    … …  
    26482453int        Tria::NumberofNodesVelocity(void){/*{{{*/
    26492454        return TriaRef::NumberofNodes(this->VelocityInterpolation());
    2650 }
    2651 /*}}}*/
    2652 void       Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm){/*{{{*/
    2653  
    2654    int        i;
    2655    IssmDouble agd[NUMVERTICES];             // surface mass balance
    2656    IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
    2657    IssmDouble yearlytemperatures[NUMVERTICES]={0.};
    2658    IssmDouble tmp[NUMVERTICES];
    2659    IssmDouble h[NUMVERTICES],s[NUMVERTICES];
    2660    IssmDouble rho_water,rho_ice,desfac,s0p,s0t,rlaps,rlapslgm;
    2661    IssmDouble PfacTime,TdiffTime,sealevTime;
    2662    IssmDouble mavg=1./12.; //factor for monthly average
    2663    
    2664    /*Get material parameters :*/
    2665    rho_water=matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    2666    rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    2667 
    2668   /*Get some pdd parameters*/
    2669   desfac=matpar->GetMaterialParameter(SurfaceforcingsDesfacEnum);
    2670   s0p=matpar->GetMaterialParameter(SurfaceforcingsS0pEnum);
    2671   s0t=matpar->GetMaterialParameter(SurfaceforcingsS0tEnum);
    2672   rlaps=matpar->GetMaterialParameter(SurfaceforcingsRlapsEnum);
    2673   rlapslgm=matpar->GetMaterialParameter(SurfaceforcingsRlapslgmEnum);
    2674 
    2675    /*Recover monthly temperatures and precipitation and compute the yearly mean temperatures*/
    2676    Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    2677    Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
    2678    GaussTria* gauss=new GaussTria();
    2679    IssmDouble time,yts;
    2680    this->parameters->FindParam(&time,TimeEnum);
    2681    this->parameters->FindParam(&yts,ConstantsYtsEnum);
    2682 
    2683    for(int month=0;month<12;month++) {
    2684      for(int iv=0;iv<NUMVERTICES;iv++) {
    2685        gauss->GaussVertex(iv);
    2686        input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
    2687        yearlytemperatures[iv]=yearlytemperatures[iv]+monthlytemperatures[iv][month]*mavg; // Has to be in Kelvin
    2688        monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius for PDD module
    2689        input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
    2690      }
    2691    }     
    2692 
    2693   /*Recover Pfac, Tdiff and sealev at time t:
    2694     This parameters are used to interpolate the temperature
    2695     and precipitaton between PD and LGM when ismungsm==1 */
    2696   if (ismungsm==1){ 
    2697     this->parameters->FindParam(&TdiffTime,SurfaceforcingsTdiffEnum,time);
    2698     this->parameters->FindParam(&sealevTime,SurfaceforcingsSealevEnum,time);
    2699   }
    2700   else {
    2701     TdiffTime=0;
    2702     sealevTime=0; 
    2703   }
    2704 
    2705   /*Recover info at the vertices: */
    2706   GetInputListOnVertices(&h[0],ThicknessEnum);
    2707   GetInputListOnVertices(&s[0],SurfaceEnum);
    2708      
    2709    /*measure the surface mass balance*/
    2710   for (int iv = 0; iv<NUMVERTICES; iv++){
    2711     agd[iv]=PddSurfaceMassBalance(&monthlytemperatures[iv][0], &monthlyprec[iv][0],
    2712                                   pdds, pds, signorm, yts, h[iv], s[iv],
    2713                                   desfac, s0t, s0p,rlaps,rlapslgm,TdiffTime,sealevTime,
    2714                                   rho_water,rho_ice);
    2715   }
    2716 
    2717    /*Update inputs*/   
    2718    // TransientInput* NewTemperatureInput = new TransientInput(SurfaceforcingsMonthlytemperaturesEnum);
    2719    // TransientInput* NewPrecipitationInput = new TransientInput(SurfaceforcingsPrecipitationEnum);
    2720    // for (int imonth=0;imonth<12;imonth++) {
    2721    //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlytemperatures[i][imonth];
    2722    //   TriaInput* newmonthinput1 = new TriaInput(SurfaceforcingsMonthlytemperaturesEnum,&tmp[0],P1Enum);
    2723    //   NewTemperatureInput->AddTimeInput(newmonthinput1,time+imonth/12.*yts);
    2724    //
    2725    //   for(i=0;i<NUMVERTICES;i++) tmp[i]=monthlyprec[i][imonth];
    2726    //   TriaInput* newmonthinput2 = new TriaInput(SurfaceforcingsPrecipitationEnum,&tmp[0],P1Enum);
    2727    //   NewPrecipitationInput->AddTimeInput(newmonthinput2,time+imonth/12.*yts);
    2728    // }
    2729    // NewTemperatureInput->Configure(this->parameters);
    2730    // NewPrecipitationInput->Configure(this->parameters);
    2731 
    2732 
    2733   this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
    2734   this->inputs->AddInput(new TriaInput(SurfaceforcingsMassBalanceEnum,&agd[0],P1Enum));
    2735    // this->inputs->AddInput(NewTemperatureInput);
    2736    // this->inputs->AddInput(NewPrecipitationInput);
    2737    // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2738 
    2739    //this->InputExtrude(SurfaceforcingsMassBalanceEnum,-1);
    2740    // this->InputExtrude(SurfaceforcingsMonthlytemperaturesEnum,-1);
    2741    // this->InputExtrude(SurfaceforcingsPrecipitationEnum,-1);
    2742 
    2743         /*clean-up*/
    2744         delete gauss;
    27452455}
    27462456/*}}}*/
    … …  
    27702480                }
    27712481        }
     2482}
     2483/*}}}*/
     2484int        Tria::PressureInterpolation(void){/*{{{*/
     2485        return TriaRef::PressureInterpolation(this->element_type);
    27722486}
    27732487/*}}}*/

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

    r19215 r19237  
    6464                void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
    6565                void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
    66                 void        Delta18oParameterization(void);
    67                 void        MungsmtpParameterization(void);
    68                 void        Delta18opdParameterization(void);
    6966                int         EdgeOnBaseIndex();
    7067                void        EdgeOnBaseIndices(int* pindex1,int* pindex);
    … …  
    114111                int         NumberofNodesPressure(void);
    115112                int         NumberofNodesVelocity(void);
    116                 void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm);
    117113                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    118114                int         PressureInterpolation();
Note: See TracChangeset for help on using the changeset viewer.