Changeset 12326


Ignore:
Timestamp:
06/01/12 17:06:57 (13 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk and trunk-jpl

Location:
issm/trunk-jpl
Files:
25 edited
4 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl

  • issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

    r12266 r12326  
    103103        MaterialsRhoIceEnum,
    104104        MaterialsRhoWaterEnum,
     105        MaterialsRhoFreshwaterEnum,
    105106        MaterialsMuWaterEnum,
    106107        MaterialsThermalExchangeVelocityEnum,
    … …  
    160161        SurfaceforcingsPrecipitationEnum,
    161162        SurfaceforcingsMassBalanceEnum,
     163        SurfaceforcingsIspddEnum,
     164        SurfaceforcingsMonthlytemperaturesEnum,
    162165        ThermalMaxiterEnum,
    163166        ThermalPenaltyFactorEnum,

  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r12319 r12326  
    109109                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    110110                case MaterialsRhoWaterEnum : return "MaterialsRhoWater";
     111                case MaterialsRhoFreshwaterEnum : return "MaterialsRhoFreshwater";
    111112                case MaterialsMuWaterEnum : return "MaterialsMuWater";
    112113                case MaterialsThermalExchangeVelocityEnum : return "MaterialsThermalExchangeVelocity";
    … …  
    166167                case SurfaceforcingsPrecipitationEnum : return "SurfaceforcingsPrecipitation";
    167168                case SurfaceforcingsMassBalanceEnum : return "SurfaceforcingsMassBalance";
     169                case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd";
     170                case SurfaceforcingsMonthlytemperaturesEnum : return "SurfaceforcingsMonthlytemperatures";
    168171                case ThermalMaxiterEnum : return "ThermalMaxiter";
    169172                case ThermalPenaltyFactorEnum : return "ThermalPenaltyFactor";

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r11964 r12326  
    8989        parameters->AddObject(iomodel->CopyConstantObject(InversionIscontrolEnum));
    9090        parameters->AddObject(iomodel->CopyConstantObject(InversionTaoEnum));
     91        parameters->AddObject(iomodel->CopyConstantObject(SurfaceforcingsIspddEnum));
    9192
    9293        /*some parameters that did not come with the iomodel: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/UpdateElementsPrognostic.cpp

    r11719 r12326  
    2020        int    stabilization;
    2121        bool   dakota_analysis;
     22        bool   ispdd;
    2223
    2324        /*Fetch data needed: */
    … …  
    2728        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    2829        iomodel->FetchData(1,MeshElementsEnum);
     30        iomodel->Constant(&ispdd,SurfaceforcingsIspddEnum);
    2931
    3032        /*Update elements: */
    … …  
    4446        iomodel->FetchDataToInput(elements,MaskElementonfloatingiceEnum);
    4547        iomodel->FetchDataToInput(elements,MaskElementonwaterEnum);
    46         iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
    4748        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
    4849        iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateCorrectionEnum);
    … …  
    6667                iomodel->FetchDataToInput(elements,TemperatureEnum);
    6768        }
     69        if(ispdd){
     70          iomodel->FetchDataToInput(elements,VyEnum);
     71          iomodel->FetchDataToInput(elements,ThermalSpctemperatureEnum);
     72          iomodel->FetchDataToInput(elements,SurfaceforcingsPrecipitationEnum);
     73          iomodel->FetchDataToInput(elements,SurfaceforcingsMonthlytemperaturesEnum);
     74        }
     75        else{
     76                iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum);
     77        }
    6878
    6979        /*Free data: */

  • issm/trunk-jpl/src/c/modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp

    r11994 r12326  
    2626  double signorm = 5.5;      // signorm : sigma of the temperature distribution for a normal day
    2727  double siglim;       // sigma limit for the integration which is equal to 2.5 sigmanorm
    28   double signormc;     // sigma of the temperature distribution for cloudy day
    29   double siglimc=0, siglim0, siglim0c;
     28  double signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
     29  double siglimc, siglim0, siglim0c;
    3030  double tstep, tsint, tint, tstepc;
    3131  int    NPDMAX = 1504, NPDCMAX = 1454;
    … …  
    4343  pds=(double*)xmalloc(NPDCMAX*sizeof(double)+1);
    4444 
    45  
    46   // PDD constant
    47   siglim = 2.5*signorm;
    48  
    4945  // initialize PDD (creation of a lookup table)
    5046  tstep = 0.1;
    … …  
    5349  snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0)));
    5450  siglim = 2.5*signorm;
     51  siglimc = 2.5*signormc;
     52  siglim0 =  siglim/DT + 0.5;
     53  siglim0c =  siglimc/DT + 0.5;
     54  PDup = siglimc+PDCUT;
     55
    5556  itm = (int)(2*siglim/DT + 1.5);
    5657 

  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r12319 r12326  
    110110              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    111111              else if (strcmp(name,"MaterialsRhoWater")==0) return MaterialsRhoWaterEnum;
     112              else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
    112113              else if (strcmp(name,"MaterialsMuWater")==0) return MaterialsMuWaterEnum;
    113114              else if (strcmp(name,"MaterialsThermalExchangeVelocity")==0) return MaterialsThermalExchangeVelocityEnum;
    … …  
    138139              else if (strcmp(name,"PrognosticMinThickness")==0) return PrognosticMinThicknessEnum;
    139140              else if (strcmp(name,"PrognosticPenaltyFactor")==0) return PrognosticPenaltyFactorEnum;
    140               else if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
    141141         else stage=2;
    142142   }
    143143   if(stage==2){
    144               if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
     144              if (strcmp(name,"PrognosticSpcthickness")==0) return PrognosticSpcthicknessEnum;
     145              else if (strcmp(name,"PrognosticStabilization")==0) return PrognosticStabilizationEnum;
    145146              else if (strcmp(name,"PrognosticVertexPairing")==0) return PrognosticVertexPairingEnum;
    146147              else if (strcmp(name,"QmuIsdakota")==0) return QmuIsdakotaEnum;
    … …  
    170171              else if (strcmp(name,"SurfaceforcingsPrecipitation")==0) return SurfaceforcingsPrecipitationEnum;
    171172              else if (strcmp(name,"SurfaceforcingsMassBalance")==0) return SurfaceforcingsMassBalanceEnum;
     173              else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum;
     174              else if (strcmp(name,"SurfaceforcingsMonthlytemperatures")==0) return SurfaceforcingsMonthlytemperaturesEnum;
    172175              else if (strcmp(name,"ThermalMaxiter")==0) return ThermalMaxiterEnum;
    173176              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    … …  
    265268   }
    266269   if(stage==3){
    267               if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
     270              if (strcmp(name,"IntVecParam")==0) return IntVecParamEnum;
     271              else if (strcmp(name,"MacAyeal2dIceFront")==0) return MacAyeal2dIceFrontEnum;
     272              else if (strcmp(name,"MacAyeal3dIceFront")==0) return MacAyeal3dIceFrontEnum;
    268273              else if (strcmp(name,"Matice")==0) return MaticeEnum;
    269274              else if (strcmp(name,"Matpar")==0) return MatparEnum;

  • issm/trunk-jpl/src/c/modules/modules.h

    r12164 r12326  
    9494#include "./ConstraintsStatex/ConstraintsStatex.h"
    9595#include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
     96#include "./PositiveDegreeDayx/PositiveDegreeDayx.h"
    9697#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
    9798#include "./Dakotax/Dakotax.h"
    … …  
    120121#include "./VerticesDofx/VerticesDofx.h"
    121122#include "./VecMergex/VecMergex.h"
    122 
    123123#endif

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

    r12014 r12326  
    22772277   // PDD and PD constants and variables
    22782278   double siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
    2279    double siglimc=0, siglim0, siglim0c;
     2279   double signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
     2280   double siglimc, siglim0, siglim0c;
    22802281   double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
    22812282   double DT = 0.02;
    … …  
    23232324   /*Get material parameters :*/
    23242325   rho_ice=matpar->GetRhoIce();
    2325    rho_water=matpar->GetRhoWater();
    2326    density=rho_ice/rho_water;
     2326   //rho_freshwater=matpar->GetRhoWater();
    23272327   
    2328    sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
     2328   sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
    23292329     
    23302330     /*PDD constant*/
    23312331   siglim = 2.5*signorm;
     2332   siglimc = 2.5*signormc;
    23322333   siglim0 =  siglim/DT + 0.5;
    23332334   siglim0c =  siglimc/DT + 0.5;

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

    r12272 r12326  
    20872087
    20882088   int    i,iqj,imonth;
    2089    double agd[NUMVERTICES];  // surface and basal
     2089   double agd[NUMVERTICES];             // surface mass balance
    20902090   double saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    20912091   double smelt[NUMVERTICES] = {0};     // yearly melt
    … …  
    20962096   double sconv; //rhow_rain/rhoi / 12 months
    20972097
    2098    double  rho_water,rho_ice,density;
    2099    double lapser=6.5/1000, sealev=0;    // lapse rate. degrees per meter.
    2100    double desfac = 0.5;                 //desert elevation factor
    2101    double s0p[NUMVERTICES]={0};         //should be set to elevation from precip source
    2102    double s0t[NUMVERTICES]={0};         //should be set to elevation from temperature source
     2098   double rho_water,rho_ice,density;
     2099   double lapser=6.5/1000, sealev=0;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
     2100   double desfac = 0.5;                 // desert elevation factor
     2101   double s0p[NUMVERTICES]={0};         // should be set to elevation from precip source
     2102   double s0t[NUMVERTICES]={0};         // should be set to elevation from temperature source
    21032103   double st;             // elevation between altitude of the temp record and current altitude
    21042104   double sp;             // elevation between altitude of the prec record and current altitude
    21052105
    2106 
    21072106   // PDD and PD constants and variables
    21082107   double siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
    2109    double siglimc=0, siglim0, siglim0c;
     2108   double signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
     2109   double siglimc, siglim0, siglim0c;
    21102110   double PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
    21112111   double DT = 0.02;
    … …  
    21242124   
    21252125   double h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
    2126    double t[NUMVERTICES][12],prec[NUMVERTICES][12];
     2126   double t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12];
    21272127   double deltm=1/12;
    21282128   int    ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
    … …  
    21432143   GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
    21442144
    2145    for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion from Kelvin to celcius
     2145        /*Recover monthly temperatures*/
     2146        Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
     2147        GaussTria* gauss=new GaussTria();
     2148        double time,yts;
     2149        this->parameters->FindParam(&time,TimeEnum);
     2150        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2151        for(int month=0;month<12;month++){
     2152                for(int iv=0;iv<NUMVERTICES;iv++){
     2153                        gauss->GaussVertex(iv);
     2154                        input->GetInputValue(&t[iv][month],gauss,time+month/12*yts);
     2155                        t[iv][month]=t[iv][month]-273.15; // conversion from Kelvin to celcius
     2156                }
     2157        }
    21462158
    21472159   for(i=0;i<NUMVERTICES;i++)
    … …  
    21532165   /*Get material parameters :*/
    21542166   rho_ice=matpar->GetRhoIce();
    2155    rho_water=matpar->GetRhoWater();
    2156    density=rho_ice/rho_water;
     2167   rho_water=matpar->GetRhoFreshwater();
    21572168   
    21582169   sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
    … …  
    21602171     /*PDD constant*/
    21612172   siglim = 2.5*signorm;
     2173   siglimc = 2.5*signormc;
    21622174   siglim0 =  siglim/DT + 0.5;
    21632175   siglim0c =  siglimc/DT + 0.5;
    … …  
    21842196       else {q = 1.0;}
    21852197       
    2186        qmt[i]= qmt[i] + prec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent
     2198       qmt[i]= qmt[i] + prec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent per month
    21872199       qmpt= q*prec[i][imonth]*sconv;           
    21882200       qmp[i]= qmp[i] + qmpt;

  • issm/trunk-jpl/src/c/objects/Inputs/BoolInput.h

    r12014 r12326  
    4949                void GetInputValue(double* pvalue);
    5050                void GetInputValue(double* pvalue,GaussTria* gauss);
     51                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5152                void GetInputValue(double* pvalue,GaussPenta* gauss);
    5253                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/ControlInput.h

    r12014 r12326  
    5454                void GetInputValue(double* pvalue);
    5555                void GetInputValue(double* pvalue,GaussTria* gauss);
     56                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5657                void GetInputValue(double* pvalue,GaussPenta* gauss);
    5758                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/DatasetInput.h

    r12014 r12326  
    4949                void GetInputValue(double* pvalue){_error_("not implemented yet");};
    5050                void GetInputValue(double* pvalue,GaussTria* gauss){_error_("not implemented yet");};
     51                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5152                void GetInputValue(double* pvalue,GaussPenta* gauss){_error_("not implemented yet");};
    5253                void GetInputValue(double* pvalue,GaussTria* gauss ,int index);

  • issm/trunk-jpl/src/c/objects/Inputs/DoubleInput.h

    r12014 r12326  
    4848                void GetInputValue(double* pvalue);
    4949                void GetInputValue(double* pvalue,GaussTria* gauss);
     50                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5051                void GetInputValue(double* pvalue,GaussPenta* gauss);
    5152                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/Input.h

    r11695 r12326  
    2727                virtual void GetInputValue(double* pvalue)=0;
    2828                virtual void GetInputValue(double* pvalue,GaussTria* gauss)=0;
     29                virtual void GetInputValue(double* pvalue,GaussTria* gauss,double time)=0;
    2930                virtual void GetInputValue(double* pvalue,GaussPenta* gauss)=0;
    3031                virtual void GetInputValue(double* pvalue,GaussTria* gauss ,int index)=0;

  • issm/trunk-jpl/src/c/objects/Inputs/IntInput.h

    r12014 r12326  
    4949                void GetInputValue(double* pvalue);
    5050                void GetInputValue(double* pvalue,GaussTria* gauss);
     51                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5152                void GetInputValue(double* pvalue,GaussPenta* gauss);
    5253                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/PentaP1Input.h

    r12014 r12326  
    4949                void GetInputValue(double* pvalue){_error_("not implemented yet");};
    5050                void GetInputValue(double* pvalue,GaussTria* gauss){_error_("not implemented yet");};
     51                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5152                void GetInputValue(double* pvalue,GaussPenta* gauss);
    5253                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/TransientInput.cpp

    r12014 r12326  
    22 * \brief: implementation of the TransientInput object
    33 */
    4 /*Headers{{{1*/
     4/*Headers{{{*/
    55#ifdef HAVE_CONFIG_H
    66        #include <config.h>
    … …  
    1919
    2020/*TransientInput constructors and destructor*/
    21 /*FUNCTION TransientInput::TransientInput(){{{1*/
     21/*FUNCTION TransientInput::TransientInput(){{{*/
    2222TransientInput::TransientInput(){
    2323
    … …  
    3030}
    3131/*}}}*/
    32 /*FUNCTION TransientInput::TransientInput(int in_enum_type){{{1*/
     32/*FUNCTION TransientInput::TransientInput(int in_enum_type){{{*/
    3333TransientInput::TransientInput(int in_enum_type)
    3434{
    … …  
    4444}
    4545/*}}}*/
    46 /*FUNCTION TransientInput::~TransientInput{{{1*/
     46/*FUNCTION TransientInput::~TransientInput{{{*/
    4747TransientInput::~TransientInput(){
    4848        xfree((void**)&this->timesteps);
    … …  
    5454}
    5555/*}}}*/
    56 /*FUNCTION void TransientInput::AddTimeInput(Input* input,double time){{{1*/
    57 void TransientInput::AddTimeInput(Input* input,double time){
    58 
    59         /*insert values at time step: */
    60         if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially");
    61 
    62         //copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
    63    double* old_timesteps=NULL;
    64 
    65         if (this->numtimesteps > 0){
    66                 old_timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
    67                 memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(double));
    68                 xfree((void**)&this->timesteps);
    69         }
    70 
    71         this->numtimesteps=this->numtimesteps+1;
    72    this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
    73 
    74         if (this->numtimesteps > 1){
    75                 memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(double));
    76                 xfree((void**)&old_timesteps);
    77         }
    78 
    79         /*go ahead and plug: */
    80         this->timesteps[this->numtimesteps-1]=time;
    81         inputs->AddObject(input);
    82 
    83 }
    84 /*}}}*/
    8556
    8657/*Object virtual functions definitions:*/
    87 /*FUNCTION TransientInput::Echo {{{1*/
     58/*FUNCTION TransientInput::Echo {{{*/
    8859void TransientInput::Echo(void){
    8960        this->DeepEcho();
    9061}
    9162/*}}}*/
    92 /*FUNCTION TransientInput::DeepEcho{{{1*/
     63/*FUNCTION TransientInput::DeepEcho{{{*/
    9364void TransientInput::DeepEcho(void){
    9465
    … …  
    10576}
    10677/*}}}*/
    107 /*FUNCTION TransientInput::Id{{{1*/
     78/*FUNCTION TransientInput::Id{{{*/
    10879int    TransientInput::Id(void){ return -1; }
    10980/*}}}*/
    110 /*FUNCTION TransientInput::MyRank{{{1*/
     81/*FUNCTION TransientInput::MyRank{{{*/
    11182int    TransientInput::MyRank(void){
    11283        extern int my_rank;
    … …  
    11485}
    11586/*}}}*/
    116 /*FUNCTION TransientInput::ObjectEnum{{{1*/
     87/*FUNCTION TransientInput::ObjectEnum{{{*/
    11788int TransientInput::ObjectEnum(void){
    11889
    … …  
    12192}
    12293/*}}}*/
    123 /*FUNCTION TransientInput::copy{{{1*/
     94/*FUNCTION TransientInput::copy{{{*/
    12495Object* TransientInput::copy() {
    12596
    … …  
    140111       
    141112/*TransientInput management*/
    142 /*FUNCTION TransientInput::InstanceEnum{{{1*/
     113/*FUNCTION TransientInput::InstanceEnum{{{*/
    143114int TransientInput::InstanceEnum(void){
    144115
    … …  
    147118}
    148119/*}}}*/
    149 /*FUNCTION TransientInput::SpawnTriaInput{{{1*/
     120/*FUNCTION TransientInput::SpawnTriaInput{{{*/
    150121Input* TransientInput::SpawnTriaInput(int* indices){
    151122
    … …  
    167138}
    168139/*}}}*/
    169 /*FUNCTION TransientInput::SpawnResult{{{1*/
     140/*FUNCTION TransientInput::SpawnResult{{{*/
    170141ElementResult* TransientInput::SpawnResult(int step, double time){
    171142
    … …  
    185156
    186157/*Object functions*/
    187 /*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){{{1*/
     158/*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){{{*/
    188159void TransientInput::GetInputValue(double* pvalue,GaussTria* gauss){
    189        
    190160        double time;
    191161
    … …  
    200170
    201171        delete input;
    202 
    203 }
    204 /*}}}*/
    205 /*FUNCTION TransientInput::GetInputDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{1*/
     172}
     173/*}}}*/
     174/*FUNCTION TransientInput::GetInputValue(double* pvalue,GaussTria* gauss,double time){{{*/
     175void TransientInput::GetInputValue(double* pvalue,GaussTria* gauss,double time){
     176
     177        /*Retrieve interpolated values for this time step: */
     178        Input* input=GetTimeInput(time);
     179
     180        /*Call input function*/
     181        input->GetInputValue(pvalue,gauss);
     182
     183        delete input;
     184}
     185/*}}}*/
     186/*FUNCTION TransientInput::GetInputDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){{{*/
    206187void TransientInput::GetInputDerivativeValue(double* p, double* xyz_list, GaussTria* gauss){
    207188
    … …  
    221202}
    222203/*}}}*/
    223 /*FUNCTION TransientInput::ChangeEnum{{{1*/
     204/*FUNCTION TransientInput::ChangeEnum{{{*/
    224205void TransientInput::ChangeEnum(int newenumtype){
    225206        this->enum_type=newenumtype;
    226207}
    227208/*}}}*/
    228 /*FUNCTION TransientInput::GetInputAverage{{{1*/
     209/*FUNCTION TransientInput::GetInputAverage{{{*/
    229210void TransientInput::GetInputAverage(double* pvalue){
    230211       
    … …  
    246227
    247228/*Intermediary*/
    248 /*FUNCTION TransientInput::SquareMin{{{1*/
     229/*FUNCTION TransientInput::AddTimeInput{{{*/
     230void TransientInput::AddTimeInput(Input* input,double time){
     231
     232        /*insert values at time step: */
     233        if (this->numtimesteps>0 && time<=this->timesteps[this->numtimesteps-1]) _assert_("timestep values must increase sequentially");
     234
     235        //copy timesteps, add the new time, delete previous timesteps, and add the new input: inputs->AddObject(input);
     236        double* old_timesteps=NULL;
     237
     238        if (this->numtimesteps > 0){
     239                old_timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
     240                memcpy(old_timesteps,this->timesteps,this->numtimesteps*sizeof(double));
     241                xfree((void**)&this->timesteps);
     242        }
     243
     244        this->numtimesteps=this->numtimesteps+1;
     245        this->timesteps=(double*)xmalloc(this->numtimesteps*sizeof(double));
     246
     247        if (this->numtimesteps > 1){
     248                memcpy(this->timesteps,old_timesteps,(this->numtimesteps-1)*sizeof(double));
     249                xfree((void**)&old_timesteps);
     250        }
     251
     252        /*go ahead and plug: */
     253        this->timesteps[this->numtimesteps-1]=time;
     254        inputs->AddObject(input);
     255
     256}
     257/*}}}*/
     258/*FUNCTION TransientInput::SquareMin{{{*/
    249259void TransientInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
    250260
    … …  
    264274}
    265275/*}}}*/
    266 /*FUNCTION TransientInput::InfinityNorm{{{1*/
     276/*FUNCTION TransientInput::InfinityNorm{{{*/
    267277double TransientInput::InfinityNorm(void){
    268278
    … …  
    284294}
    285295/*}}}*/
    286 /*FUNCTION TransientInput::Max{{{1*/
     296/*FUNCTION TransientInput::Max{{{*/
    287297double TransientInput::Max(void){
    288298
    … …  
    304314}
    305315/*}}}*/
    306 /*FUNCTION TransientInput::MaxAbs{{{1*/
     316/*FUNCTION TransientInput::MaxAbs{{{*/
    307317double TransientInput::MaxAbs(void){
    308318
    … …  
    325335}
    326336/*}}}*/
    327 /*FUNCTION TransientInput::Min{{{1*/
     337/*FUNCTION TransientInput::Min{{{*/
    328338double TransientInput::Min(void){
    329339
    … …  
    346356}
    347357/*}}}*/
    348 /*FUNCTION TransientInput::MinAbs{{{1*/
     358/*FUNCTION TransientInput::MinAbs{{{*/
    349359double TransientInput::MinAbs(void){
    350360
    … …  
    366376}
    367377/*}}}*/
    368 /*FUNCTION TransientInput::GetVectorFromInputs{{{1*/
     378/*FUNCTION TransientInput::GetVectorFromInputs{{{*/
    369379void TransientInput::GetVectorFromInputs(Vector* vector,int* doflist){
    370380
    … …  
    383393
    384394} /*}}}*/
    385 /*FUNCTION TransientInput::GetTimeInput{{{1*/
     395/*FUNCTION TransientInput::GetTimeInput{{{*/
    386396Input* TransientInput::GetTimeInput(double intime){
    387397
    … …  
    442452}
    443453/*}}}*/
    444 /*FUNCTION TransientInput::Configure{{{1*/
     454/*FUNCTION TransientInput::Configure{{{*/
    445455void TransientInput::Configure(Parameters* parameters){
    446456        this->parameters=parameters;

  • issm/trunk-jpl/src/c/objects/Inputs/TransientInput.h

    r12014 r12326  
    5151                void GetInputValue(double* pvalue){_error_("not implemented yet");};
    5252                void GetInputValue(double* pvalue,GaussTria* gauss);
     53                void GetInputValue(double* pvalue,GaussTria* gauss,double time);
    5354                void GetInputValue(double* pvalue,GaussPenta* gauss){_error_("not implemented yet");};
    5455                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Inputs/TriaP1Input.h

    r12014 r12326  
    4949                void GetInputValue(double* pvalue){_error_("not implemented yet");}
    5050                void GetInputValue(double* pvalue,GaussTria* gauss);
     51                void GetInputValue(double* pvalue,GaussTria* gauss,double time){_error_("not implemented yet");};
    5152                void GetInputValue(double* pvalue,GaussPenta* gauss){_error_("not implemented yet");};
    5253                void GetInputValue(double* pvalue,GaussTria* gauss ,int index){_error_("not implemented yet");};

  • issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp

    r12014 r12326  
    2525Matpar::Matpar(int matpar_mid, IoModel* iomodel){
    2626
    27         this->mid                       = matpar_mid;
     27        this->mid = matpar_mid;
    2828        iomodel->Constant(&this->rho_ice,MaterialsRhoIceEnum);
    2929        iomodel->Constant(&this->rho_water,MaterialsRhoWaterEnum);
     30        iomodel->Constant(&this->rho_freshwater,MaterialsRhoFreshwaterEnum);
    3031        iomodel->Constant(&this->mu_water,MaterialsMuWaterEnum);
    3132        iomodel->Constant(&this->heatcapacity,MaterialsHeatcapacityEnum);
    … …  
    6061        printf("   rho_ice: %g\n",rho_ice);
    6162        printf("   rho_water: %g\n",rho_water);
     63        printf("   rho_freshwater: %g\n",rho_freshwater);
    6264        printf("   mu_water: %g\n",mu_water);
    6365        printf("   heatcapacity: %g\n",heatcapacity);
    … …  
    7678void Matpar::DeepEcho(void){
    7779
    78         printf("Matpar:\n");
    79         printf("   mid: %i\n",mid);
    80         printf("   rho_ice: %g\n",rho_ice);
    81         printf("   rho_water: %g\n",rho_water);
    82         printf("   mu_water: %g\n",mu_water);
    83         printf("   heatcapacity: %g\n",heatcapacity);
    84         printf("   thermalconductivity: %g\n",thermalconductivity);
    85         printf("   latentheat: %g\n",latentheat);
    86         printf("   beta: %g\n",beta);
    87         printf("   meltingpoint: %g\n",meltingpoint);
    88         printf("   referencetemperature: %g\n",referencetemperature);
    89         printf("   mixed_layer_capacity: %g\n",mixed_layer_capacity);
    90         printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
    91         printf("   g: %g\n",g);
    92         return;
     80        this->Echo();
    9381}               
    9482/*}}}1*/
    … …  
    158146                        this->rho_ice=constant;
    159147                        break;
    160 
    161148                case MaterialsRhoWaterEnum:
    162149                        this->rho_water=constant;
    163150                        break;
    164 
     151                case MaterialsRhoFreshwaterEnum:
     152                        this->rho_freshwater=constant;
     153                        break;
    165154                case MaterialsMuWaterEnum:
    166155                        this->mu_water=constant;
    167156                        break;
    168 
    169157                case MaterialsHeatcapacityEnum:
    170158                        this->heatcapacity=constant;
    171159                        break;
    172 
    173160                case MaterialsThermalconductivityEnum:
    174161                        this->thermalconductivity=constant;
    175162                        break;
    176 
    177163                case  MaterialsLatentheatEnum:
    178164                        this->latentheat=constant;
    179165                        break;
    180 
    181166                case  MaterialsBetaEnum:
    182167                        this->beta=constant;
    183168                        break;
    184 
    185169                case  MaterialsMeltingpointEnum:
    186170                        this->meltingpoint=constant;
    187171                        break;
    188 
    189172                case  ConstantsReferencetemperatureEnum:
    190173                        this->referencetemperature=constant;
    191174                        break;
    192 
    193175                case  MaterialsMixedLayerCapacityEnum:
    194176                        this->mixed_layer_capacity=constant;
    195177                        break;
    196 
    197178                case  MaterialsThermalExchangeVelocityEnum:
    198179                        this->thermalconductivity=constant;
    199180                        break;
    200 
    201181                case  ConstantsGEnum:
    202182                        this->g=constant;
    203183                        break;
    204 
    205184                default:
    206185                        break;
    … …  
    277256double Matpar::GetRhoWater(){
    278257        return rho_water;
     258}
     259/*}}}1*/
     260/*FUNCTION Matpar::GetRhoFreshwater {{{1*/
     261double Matpar::GetRhoFreshwater(){
     262        return rho_freshwater;
    279263}
    280264/*}}}1*/

  • issm/trunk-jpl/src/c/objects/Materials/Matpar.h

    r12014 r12326  
    1515
    1616        private:
    17                 int     mid;
    18                 double  rho_ice;
    19                 double  rho_water;
     17                int       mid;
     18                double  rho_ice;
     19                double  rho_water;
     20                double  rho_freshwater;
    2021                double  mu_water;
    2122                double  heatcapacity;
    … …  
    7273                double GetRhoIce();
    7374                double GetRhoWater();
     75                double GetRhoFreshwater();
    7476                double GetMuWater();
    7577                double GetMixedLayerCapacity();

  • issm/trunk-jpl/src/c/solutions/prognostic_core.cpp

    r11832 r12326  
    1616        /*parameters: */
    1717        bool save_results;
     18        bool ispdd;
    1819
    1920        /*activate formulation: */
    … …  
    2223        /*recover parameters: */
    2324        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
     25        femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum);
     26
     27        if(ispdd){
     28          _printf_(VerboseSolution(),"   call positive degree day module\n");
     29          PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
     30        }       
    2431
    2532        _printf_(VerboseSolution(),"   call computational core\n");
    2633        solver_linear(femmodel);
    27                
     34       
    2835        if(save_results){
    2936                _printf_(VerboseSolution(),"   saving results\n");

  • issm/trunk-jpl/src/m/classes/materials.m

    r11869 r12326  
    88                rho_ice                    = 0;
    99                rho_water                  = 0;
     10                rho_freshwater             = 0;
    1011                mu_water                   = 0;
    1112                heatcapacity               = 0;
    … …  
    3435                        obj.rho_ice=917;
    3536
    36                         %water density (kg/m^3)
     37                        %ocean water density (kg/m^3)
    3738                        obj.rho_water=1023;
     39
     40                        %fresh water density (kg/m^3)
     41                        obj.rho_freshwater=1000;
    3842
    3943                        %water viscosity (N.s/m^2)
    … …  
    6872                        checkfield(md,'materials.rho_ice','>',0);
    6973                        checkfield(md,'materials.rho_water','>',0);
     74                        checkfield(md,'materials.rho_freshwater','>',0);
    7075                        checkfield(md,'materials.mu_water','>',0);
    7176                        checkfield(md,'materials.rheology_B','>',0,'size',[md.mesh.numberofvertices 1]);
    … …  
    7782
    7883                        fielddisplay(obj,'rho_ice','ice density [kg/m^3]');
    79                         fielddisplay(obj,'rho_water','water density [kg/m^3]');
     84                        fielddisplay(obj,'rho_water','ocean water density [kg/m^3]');
     85                        fielddisplay(obj,'freshrho_water','fresh water density [kg/m^3]');
    8086                        fielddisplay(obj,'mu_water','water viscosity [N s/m^2]');
    8187                        fielddisplay(obj,'heatcapacity','heat capacity [J/kg/K]');
    … …  
    9399                        WriteData(fid,'object',obj,'fieldname','rho_ice','format','Double');
    94100                        WriteData(fid,'object',obj,'fieldname','rho_water','format','Double');
     101                        WriteData(fid,'object',obj,'fieldname','rho_freshwater','format','Double');
    95102                        WriteData(fid,'object',obj,'fieldname','mu_water','format','Double');
    96103                        WriteData(fid,'object',obj,'fieldname','heatcapacity','format','Double');

  • issm/trunk-jpl/src/m/classes/surfaceforcings.m

    r11869 r12326  
    88                precipitation = NaN;
    99                mass_balance  = NaN;
     10                ispdd = 0;
     11                monthlytemperatures = NaN;
    1012        end
    1113        methods
    … …  
    1921                end % }}}
    2022                function obj = setdefaultparameters(obj) % {{{
     23                 
     24                  %pdd method not used in default mode
     25                  obj.ispdd=0;
    2126
    2227                end % }}}
    … …  
    2429
    2530                        if ismember(PrognosticAnalysisEnum,analyses),
    26                                 checkfield(md,'surfaceforcings.mass_balance','forcing',1,'NaN',1);
     31                                checkfield(md,'surfaceforcings.ispdd','numel',1,'values',[0 1]);
     32                                if(obj.ispdd)
     33                                        checkfield(md,'surfaceforcings.monthlytemperatures','forcing',1,'NaN',1);
     34                                else
     35                                        checkfield(md,'surfaceforcings.mass_balance','forcing',1,'NaN',1);
     36                                end
    2737                        end
    2838                        if ismember(BalancethicknessAnalysisEnum,analyses),
    … …  
    3545                        fielddisplay(obj,'precipitation','surface precipitation [m/yr water eq]');
    3646                        fielddisplay(obj,'mass_balance','surface mass balance [m/yr ice eq]');
     47                        fielddisplay(obj,'ispdd','is pdd activated (0 or 1, default is 0)');
     48                        fielddisplay(obj,'monthlytemperatures','monthly surface temperatures required if pdd is activated');
    3749
    3850                end % }}}
    … …  
    4052                        WriteData(fid,'object',obj,'fieldname','precipitation','format','DoubleMat','mattype',1);
    4153                        WriteData(fid,'object',obj,'fieldname','mass_balance','format','DoubleMat','mattype',1);
     54                        WriteData(fid,'object',obj,'fieldname','ispdd','format','Boolean');
     55                        if obj.ispdd,
     56                                WriteData(fid,'object',obj,'fieldname','monthlytemperatures','format','DoubleMat','mattype',1);
     57                        end
     58
    4259                end % }}}
    4360        end
Note: See TracChangeset for help on using the changeset viewer.