Changeset 12704


Ignore:
Timestamp:
07/24/12 10:25:08 (13 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk and trunk-jpl

Location:
issm/trunk-jpl
Files:
26 edited
1 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl

  • issm/trunk-jpl/externalpackages/vim/addons/vim/syntax/c.vim

    r12672 r12704  
    706706syn keyword cConstant SurfaceforcingsPrecipitationEnum
    707707syn keyword cConstant SurfaceforcingsMassBalanceEnum
    708 syn keyword cConstant SurfaceforcingsIspddEnum
    709 syn keyword cConstant SurfaceforcingsIssmbgradientsEnum
    710 syn keyword cConstant SurfaceforcingsMonthlytemperaturesEnum
    711 syn keyword cConstant SurfaceforcingsHcEnum
    712 syn keyword cConstant SurfaceforcingsSmbPosMaxEnum
    713 syn keyword cConstant SurfaceforcingsSmbPosMinEnum
    714 syn keyword cConstant SurfaceforcingsAPosEnum
    715 syn keyword cConstant SurfaceforcingsBPosEnum
    716 syn keyword cConstant SurfaceforcingsANegEnum
    717 syn keyword cConstant SurfaceforcingsBNegEnum
    718708syn keyword cConstant ThermalMaxiterEnum
    719709syn keyword cConstant ThermalPenaltyFactorEnum

  • issm/trunk-jpl/src/c/Container/Elements.cpp

    r12493 r12704  
    237237                        for(int j=0;j<this->Size();j++){
    238238                                Element* element=(Element*)this->GetObjectByOffset(j);
    239                                 element->GetVectorFromResults(vector,i,resultssizes[i]);
     239                                element->GetVectorFromResults(vector,i,resultsenums[i],resultssizes[i]);
    240240                        }
    241241                        vector->Assemble();

  • issm/trunk-jpl/src/c/Makefile.am

    r12696 r12704  
    205205                                        ./shared/Elements/GetGlobalDofList.cpp\
    206206                                        ./shared/Elements/GetNumberOfDofs.cpp\
     207                                        ./shared/Elements/PddSurfaceMassBalance.cpp\
    207208                                        ./shared/String/sharedstring.h\
    208209                                        ./shared/Wrapper/wrappershared.h\

  • issm/trunk-jpl/src/c/matlab/io/OptionParse.cpp

    r12519 r12704  
    88#endif
    99
     10#include <cstring>
     11#include <mex.h>
    1012#include "../../shared/shared.h"
    1113#include "../../io/io.h"
    1214#include "../../include/include.h"
    1315#include "./matlabio.h"
    14 
    15 #include <mex.h>
    1616
    1717/*FUNCTION OptionDoubleParse {{{*/

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

    r12572 r12704  
    6464    tstar = it*DT-siglim;
    6565    tint = tsint;
    66     pddt = 0;
     66    pddt = 0.;
    6767    for ( jj = 0; jj < 600; jj++){
    6868      if (tint > (tstar+siglim)){break;}
    … …  
    9090    //    tstar = REAL(it)*DT-siglimc;
    9191    tint = tsint;          // start against upper bound
    92     pd = 0;
     92    pd = 0.;
    9393    for (jj = 0; jj < 600; jj++){
    9494      if (tint<(tstar-siglimc)) {break;}
    … …  
    9898    pds[it] = pd*snormfac;  // gaussian integral lookup table for snow fraction
    9999  }
    100   pds[itm+1] = 0;
     100  pds[itm+1] = 0.;
    101101  //     *******END initialize PDD
    102102 

  • issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.cpp

    r12561 r12704  
    126126}
    127127/*}}}*/
     128/*FUNCTION DoubleElementResult::GetVectorFromResults{{{1*/
     129void DoubleElementResult::GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){
     130
     131        _error_("cannot return vector on vertices");
     132} /*}}}*/
     133/*FUNCTION DoubleElementResult::GetElementVectorFromResults{{{1*/
     134void DoubleElementResult::GetElementVectorFromResults(Vector* vector,int dof){
     135
     136        vector->SetValue(dof,value,INS_VAL);
     137} /*}}}*/

  • issm/trunk-jpl/src/c/objects/ElementResults/DoubleElementResult.h

    r12561 r12704  
    4848                /*DoubleElementResult management: {{{*/
    4949                int   InstanceEnum();
    50                 void GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs){_error2_("not implemented");};
    51                 void GetElementVectorFromResults(Vector* vector,int dof){_error2_("not implemented");};
     50                void GetVectorFromResults(Vector* vector,int* doflist,int* connectivitylist,int numdofs);
     51                void GetElementVectorFromResults(Vector* vector,int dof);
    5252                /*}}}*/
    5353};

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

    r12677 r12704  
    6363                virtual void   InputScale(int enum_type,IssmDouble scale_factor)=0;
    6464                virtual void   GetVectorFromInputs(Vector* vector, int name_enum)=0;
    65                 virtual void   GetVectorFromResults(Vector* vector,int id,int interp)=0;
     65                virtual void   GetVectorFromResults(Vector* vector,int id,int enum_in,int interp)=0;
    6666                virtual void   InputArtificialNoise(int enum_type,IssmDouble min,IssmDouble max)=0;
    6767                virtual bool   InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums)=0;

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

    r12677 r12704  
    10801080/*}}}*/
    10811081/*FUNCTION Penta::GetVectorFromResults{{{*/
    1082 void  Penta::GetVectorFromResults(Vector* vector,int offset,int interp){
     1082void  Penta::GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp){
    10831083
    10841084        /*Get result*/
    10851085        ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset);
     1086        if(elementresult->InstanceEnum()!=enum_in){
     1087                _error_("Results of offset %i is %s, when %s was expected",offset,EnumToStringx(elementresult->InstanceEnum()),EnumToStringx(enum_in));
     1088        } 
    10861089        if(interp==P1Enum){
    10871090                int doflist1[NUMVERTICES];
    … …  
    22552258void  Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){
    22562259
    2257 
    2258    int    i,iqj,imonth;
    2259    IssmDouble agd[NUMVERTICES];  // surface and basal
    2260    IssmDouble saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    2261    IssmDouble smelt[NUMVERTICES] = {0};     // yearly melt
    2262    IssmDouble precrunoff[NUMVERTICES];      // yearly runoff
    2263    IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
    2264    IssmDouble water; //water=rain + snowmelt
    2265    IssmDouble runoff; //meltwater only, does not include rain
    2266    IssmDouble sconv; //rhow_rain/rhoi / 12 months
    2267 
    2268    IssmDouble  rho_water,rho_ice,density;
    2269    IssmDouble lapser=6.5/1000, sealev=0;    // lapse rate. degrees per meter.
    2270    IssmDouble desfac = 0.5;                 //desert elevation factor
    2271    IssmDouble s0p[NUMVERTICES]={0};         //should be set to elevation from precip source
    2272    IssmDouble s0t[NUMVERTICES]={0};         //should be set to elevation from temperature source
    2273    IssmDouble st;             // elevation between altitude of the temp record and current altitude
    2274    IssmDouble sp;             // elevation between altitude of the prec record and current altitude
    2275 
    2276 
    2277    // PDD and PD constants and variables
    2278    IssmDouble siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
    2279    IssmDouble signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
    2280    IssmDouble siglimc, siglim0, siglim0c;
    2281    IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
    2282    IssmDouble DT = 0.02;
    2283    IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
    2284    
    2285    IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
    2286    IssmDouble qm[NUMVERTICES] = {0};        // snow part of the precipitation
    2287    IssmDouble qmt[NUMVERTICES] = {0};       // precipitation without desertification effect adjustment
    2288    IssmDouble qmp[NUMVERTICES] = {0};       // desertification taken into account
    2289    IssmDouble pdd[NUMVERTICES] = {0};     
    2290    IssmDouble frzndd[NUMVERTICES] = {0}; 
    2291 
    2292    IssmDouble tstar;                        // monthly mean surface temp
    2293    IssmDouble Tsum[NUMVERTICES]= {0};       // average summer (JJA) temperature
    2294    IssmDouble Tsurf[NUMVERTICES] = {0};     // average annual temperature   
    2295    
    2296    IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
    2297    IssmDouble t[NUMVERTICES][12],prec[NUMVERTICES][12];
    2298    IssmDouble deltm=1/12;
    2299    int    ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
    2300 
    2301    IssmDouble snwm;  // snow that could have been melted in a year.
    2302    IssmDouble snwmf; //  ablation factor for snow per positive degree day.
    2303    IssmDouble smf;   //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
    2304 
    2305    IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
    2306    IssmDouble supice,supcap,diffndd;
    2307    IssmDouble fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
    2308    IssmDouble pddtj[NUMVERTICES], hmx2;
    2309 
    2310    /*Recover info at the vertices: */
    2311    GetInputListOnVertices(&h[0],ThicknessEnum);
    2312    GetInputListOnVertices(&s[0],SurfaceEnum);
    2313    GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
    2314    GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
    2315 
    2316    for(i=0;i<NUMVERTICES;i++) ttmp[i]=ttmp[i]-273.15; // convertion from Kelvin to celcius
    2317 
    2318    for(i=0;i<NUMVERTICES;i++)
    2319      for(imonth=0;imonth<12;imonth++){
    2320        t[i][imonth]=ttmp[i];
    2321        prec[i][imonth]=prectmp[i];
     2260   IssmDouble agd[NUMVERTICES];             // surface mass balance
     2261   IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
     2262   IssmDouble h[NUMVERTICES],s[NUMVERTICES]; // ,b
     2263   IssmDouble rho_water,rho_ice;
     2264
     2265   /*Recover monthly temperatures and precipitation*/
     2266   Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
     2267   Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
     2268   GaussPenta* gauss=new GaussPenta();
     2269   IssmDouble time,yts;
     2270   this->parameters->FindParam(&time,TimeEnum);
     2271   this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2272   for(int month=0;month<12;month++) {
     2273     for(int iv=0;iv<NUMVERTICES;iv++) {
     2274       gauss->GaussVertex(iv);
     2275       input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
     2276       monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
     2277       input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
     2278       monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
     2279     }
     2280   }
     2281
     2282  /*Recover info at the vertices: */
     2283  GetInputListOnVertices(&h[0],ThicknessEnum);
     2284  GetInputListOnVertices(&s[0],SurfaceEnum);
     2285
     2286  /*Get material parameters :*/
     2287  rho_ice=matpar->GetRhoIce();
     2288  rho_water=matpar->GetRhoFreshwater();
     2289
     2290   /*measure the surface mass balance*/
     2291   for (int iv = 0; iv < NUMVERTICES; iv++){
     2292     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water);
    23222293   }
    23232294
    2324    /*Get material parameters :*/
    2325    rho_ice=matpar->GetRhoIce();
    2326    //rho_freshwater=matpar->GetRhoWater();
    2327    
    2328    sconv=(1000/rho_ice)/12.; //rhow_rain/rhoi / 12 months
    2329      
    2330      /*PDD constant*/
    2331    siglim = 2.5*signorm;
    2332    siglimc = 2.5*signormc;
    2333    siglim0 =  siglim/DT + 0.5;
    2334    siglim0c =  siglimc/DT + 0.5;
    2335    PDup = siglimc+PDCUT;
    2336    
    2337    // seasonal loop
    2338    for (iqj = 0; iqj < 12; iqj++){
    2339      imonth =  ismon[iqj];
    2340      for (i = 0; i < NUMVERTICES; i++){
    2341        st=(s[i]-s0t[i])/1000;
    2342        tstar = t[i][imonth] - lapser *max(st,sealev);
    2343        Tsurf[i] = tstar*deltm+Tsurf[i];       
    2344        
    2345        /*********compute PD ****************/
    2346        if (tstar < PDup){
    2347          pd = 1;
    2348          if (tstar >= -siglimc){ pd = pds[int(tstar/DT + siglim0c)];}}
    2349        else {
    2350          pd = 0;}
    2351        
    2352        /******exp des/elev precip reduction*******/
    2353        sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
    2354        if (sp>0.0){q = exp(-desfac*sp);}
    2355        else {q = 1.0;}
    2356        
    2357        qmt[i]= qmt[i] + prec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent
    2358        qmpt= q*prec[i][imonth]*sconv;           
    2359        qmp[i]= qmp[i] + qmpt;
    2360        qm[i]= qm[i] + qmpt*pd;
    2361        
    2362        /*********compute PDD************/
    2363        // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
    2364        // gaussian=T_m, so ndd=-(Tsurf-pdd)
    2365        if (iqj>6 &&  iqj<10){ Tsum[i]=Tsum[i]+tstar;}
    2366        if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
    2367        else if (tstar> -siglim){
    2368          pddsig=pdds[int(tstar/DT + siglim0)];
    2369          pdd[i] = pdd[i] + pddsig*deltm;
    2370          frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
    2371        else{frzndd[i] = frzndd[i] - tstar*deltm; }
    2372      }
    2373    } // end of seasonal loop
    2374    
    2375    //******************************************************************
    2376    for(i=0;i<NUMVERTICES;i++){
    2377      saccu[i] = qm[i];
    2378      prect = qmp[i];     // total precipitation during 1 year taking into account des. ef.
    2379      Tsum[i]=Tsum[i]/3;
    2380      
    2381      /***** determine PDD factors *****/
    2382      if(Tsum[i]< -1.) {
    2383        snwmf=2.65*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
    2384        smf=17.22*0.001;    //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
    2385      }
    2386      else if(Tsum[i]< 10){
    2387        snwmf = (0.15*Tsum[i] + 2.8)*0.001;
    2388        smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
    2389      }
    2390      else{
    2391        snwmf=4.3*0.001;
    2392        smf=8.3*0.001;
    2393      }
    2394      snwmf=0.95*snwmf;
    2395      smf=0.95*smf;
    2396      
    2397      /*****  compute PDD ablation and refreezing *****/
    2398      pddt = pdd[i] *365;
    2399      snwm = snwmf*pddt;       // snow that could have been melted in a year
    2400      hmx2 = min(h[i],dfrz);   // refreeze active layer max depth: dfrz
    2401      
    2402      if(snwm < saccu[i]) {
    2403        water=prect-saccu[i] + snwm; //water=rain + snowmelt
    2404        //     l 2.2= capillary factor
    2405        //     Should refreezing be controlled by frzndd or by mean annual Tsurf?
    2406        //     dCovrLm concept is of warming of active layer (thickness =d@=1-
    2407        //     >2m)
    2408        //     problem with water seepage into ice: should be sealed after
    2409        //     refreezing
    2410        //     so everything needs to be predicated on 1 year scale, except for
    2411        //     thermal
    2412        //     conductivity through ice
    2413        //     also, need to account that melt season has low accum, so what's
    2414        //     going to
    2415        //     hold the meltwater around for refreezing? And melt-time will have
    2416        //     low seasonal frzndd
    2417        
    2418        //      Superimposed ice :  Pfeffer et al. 1991, Tarasov 2002
    2419        
    2420        supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
    2421        supcap=min(2.2*(saccu[i]-snwm),water);
    2422        runoff=snwm - supice;  //meltwater only, does not include rain
    2423      }
    2424      else {  //all snow melted
    2425        supice= min(hmx2*CovrLm*frzndd[i], prect );
    2426        runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
    2427        supcap=0;
    2428      }
    2429      //     pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
    2430      //     except pdd melt heat source is atmosphere, while refreeze is
    2431      //     ground/ice stored interim
    2432      //     assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
    2433      //     assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
    2434      //     <0
    2435      //     assume ndd>pdd, little melt => little supice
    2436      //     bottom line: compare for Tsurf<0 : supice and no supice case,
    2437      //     expect Tsurf difference
    2438      //     except some of cooling flux comes from atmosphere//
    2439      //     1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
    2440      //     does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
    2441      //     < 0.1
    2442      
    2443      //     make more sense to just use residual pdd-ndd except that pdd
    2444      //     residual not clear yet
    2445      //     frzndd should not be used up by refreezing in snow, so stick in
    2446      //     supcap.
    2447      diffndd=0;
    2448      if (frzndd[i]>0) {
    2449        diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
    2450        frzndd[i]=frzndd[i]-diffndd;
    2451      }
    2452      if(runoff<0){
    2453        saccu[i]= saccu[i] -runoff;
    2454        smelt[i] = 0;
    2455        precrunoff[i]=prect-saccu[i];
    2456        //here assume pdd residual is 0, =>
    2457        Tsurf[i]= max(Tsurf[i],-frzndd[i]);
    2458      }
    2459      else {
    2460        smelt[i] = runoff;
    2461        precrunoff[i]=prect-max(0.,supice)-saccu[i];}
    2462      //here really need pdd balance, try 0.5 fudge factor?
    2463      //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
    2464      //yet from site plots, can be ice free with Tsurf=-5.5C
    2465      if(Tsurf[i]<0) {
    2466        Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
    2467      
    2468      agd[i] = -smelt[i]+saccu[i];
    2469      pddtj[i]=pddt;
    2470      
    2471      /*Update inputs*/   
    2472      this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
    2473      // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2474      this->inputs->AddInput(new PentaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
    2475      
    2476    }       //end of the for loop over the vertices
     2295   /*Update inputs*/   
     2296   this->inputs->AddInput(new PentaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
     2297   //this->inputs->AddInput(new PentaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
     2298   this->InputExtrude(SurfaceforcingsMassBalanceEnum,ElementEnum);
    24772299}
    24782300/*}}}*/

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

    r12677 r12704  
    8989                IssmDouble GetZcoord(GaussPenta* gauss);
    9090                void   GetVectorFromInputs(Vector* vector,int name_enum);
    91                 void   GetVectorFromResults(Vector* vector,int offset,int interp);
     91                void   GetVectorFromResults(Vector* vector,int offset,int name_enum,int interp);
    9292               
    9393                int    Sid();

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

    r12677 r12704  
    11921192/*}}}*/
    11931193/*FUNCTION Tria::GetVectorFromResults{{{*/
    1194 void  Tria::GetVectorFromResults(Vector* vector,int offset,int interp){
     1194void  Tria::GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp){
    11951195
    11961196        /*Get result*/
    11971197        ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset);
     1198        if(elementresult->InstanceEnum()!=enum_in){
     1199                _error_("Results of offset %i is %s, when %s was expected",offset,EnumToStringx(elementresult->InstanceEnum()),EnumToStringx(enum_in));
     1200        }
    11981201        if(interp==P1Enum){
    11991202                int doflist1[NUMVERTICES];
    … …  
    20862089void  Tria::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){
    20872090
    2088    int    i,iqj,imonth;
    20892091   IssmDouble agd[NUMVERTICES];             // surface mass balance
    2090    IssmDouble saccu[NUMVERTICES] = {0};     // yearly surface accumulation
    2091    IssmDouble smelt[NUMVERTICES] = {0};     // yearly melt
    2092    IssmDouble precrunoff[NUMVERTICES];      // yearly runoff
    2093    IssmDouble prect; // total precipitation during 1 year taking into account des. ef.
    2094    IssmDouble water; //water=rain + snowmelt
    2095    IssmDouble runoff; //meltwater only, does not include rain
    2096    IssmDouble sconv; //rhow_rain/rhoi / 12 months
    2097 
    2098    IssmDouble rho_water,rho_ice,density;
    2099    IssmDouble lapser=6.5/1000, sealev=0;    // lapse rate. degrees per meter. 7.5 lev's 99 paper, 9 Marshall 99 paper
    2100    IssmDouble desfac = 0.5;                 // desert elevation factor
    2101    IssmDouble s0p[NUMVERTICES]={0};         // should be set to elevation from precip source
    2102    IssmDouble s0t[NUMVERTICES]={0};         // should be set to elevation from temperature source
    2103    IssmDouble st;             // elevation between altitude of the temp record and current altitude
    2104    IssmDouble sp;             // elevation between altitude of the prec record and current altitude
    2105 
    2106    // PDD and PD constants and variables
    2107    IssmDouble siglim;          // sigma limit for the integration which is equal to 2.5 sigmanorm
    2108    IssmDouble signormc = signorm - 0.5;     // sigma of the temperature distribution for cloudy day
    2109    IssmDouble siglimc, siglim0, siglim0c;
    2110    IssmDouble PDup, pddsig, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C)
    2111    IssmDouble DT = 0.02;
    2112    IssmDouble pddt, pd; // pd: snow/precip fraction, precipitation falling as snow
    2113    
    2114    IssmDouble q, qmpt; // q is desert/elev. fact, hnpfac is huybrect fact, and pd is normal dist.
    2115    IssmDouble qm[NUMVERTICES] = {0};        // snow part of the precipitation
    2116    IssmDouble qmt[NUMVERTICES] = {0};       // precipitation without desertification effect adjustment
    2117    IssmDouble qmp[NUMVERTICES] = {0};       // desertification taken into account
    2118    IssmDouble pdd[NUMVERTICES] = {0};     
    2119    IssmDouble frzndd[NUMVERTICES] = {0}; 
    2120 
    2121    IssmDouble tstar;                        // monthly mean surface temp
    2122    IssmDouble Tsum[NUMVERTICES]= {0};       // average summer (JJA) temperature
    2123    IssmDouble Tsurf[NUMVERTICES] = {0};     // average annual temperature   
    2124    
    2125    IssmDouble h[NUMVERTICES],s[NUMVERTICES],ttmp[NUMVERTICES],prectmp[NUMVERTICES]; // ,b[NUMVERTICES]
    2126    IssmDouble t[NUMVERTICES+1][12],prec[NUMVERTICES+1][12];
    2127    IssmDouble deltm=1/12;
    2128    int    ismon[12]={12,1,2,3,4,5,6,7,8,9,10,11};
    2129 
    2130    IssmDouble snwm;  // snow that could have been melted in a year.
    2131    IssmDouble snwmf; //  ablation factor for snow per positive degree day.
    2132    IssmDouble smf;   //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002).
    2133 
    2134    IssmDouble dfrz=1.5, CovrLm=2009./3.35e+5, dCovrLm=dfrz*CovrLm; //m*J kg^-1 C^-1 /(J kg^-1)=m/C yr
    2135    IssmDouble supice,supcap,diffndd;
    2136    IssmDouble fsupT=0.5,  fsupndd=0.5;  // Tsurf mode factors for supice
    2137    IssmDouble pddtj[NUMVERTICES], hmx2;
    2138 
    2139    /*Recover info at the vertices: */
    2140    GetInputListOnVertices(&h[0],ThicknessEnum);
    2141    GetInputListOnVertices(&s[0],SurfaceEnum);
    2142    GetInputListOnVertices(&ttmp[0],ThermalSpctemperatureEnum);
    2143    GetInputListOnVertices(&prectmp[0],SurfaceforcingsPrecipitationEnum);
    2144 
    2145         /*Recover monthly temperatures*/
    2146         Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
    2147         GaussTria* gauss=new GaussTria();
    2148         IssmDouble 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         }
    2158 
    2159    for(i=0;i<NUMVERTICES;i++)
    2160      for(imonth=0;imonth<12;imonth++){
    2161        t[i][imonth]=ttmp[i];
    2162        prec[i][imonth]=prectmp[i];
     2092   IssmDouble monthlytemperatures[NUMVERTICES][12],monthlyprec[NUMVERTICES][12];
     2093   IssmDouble h[NUMVERTICES],s[NUMVERTICES]; // ,b
     2094   IssmDouble rho_water,rho_ice;
     2095
     2096   /*Recover monthly temperatures and precipitation*/
     2097   Input*     input=inputs->GetInput(SurfaceforcingsMonthlytemperaturesEnum); _assert_(input);
     2098   Input*     input2=inputs->GetInput(SurfaceforcingsPrecipitationEnum); _assert_(input2);
     2099   GaussTria* gauss=new GaussTria();
     2100   IssmDouble time,yts;
     2101   this->parameters->FindParam(&time,TimeEnum);
     2102   this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2103   for(int month=0;month<12;month++) {
     2104     for(int iv=0;iv<NUMVERTICES;iv++) {
     2105       gauss->GaussVertex(iv);
     2106       input->GetInputValue(&monthlytemperatures[iv][month],gauss,time+month/12.*yts);
     2107       monthlytemperatures[iv][month]=monthlytemperatures[iv][month]-273.15; // conversion from Kelvin to celcius
     2108       input2->GetInputValue(&monthlyprec[iv][month],gauss,time+month/12.*yts);
     2109       monthlyprec[iv][month]=monthlyprec[iv][month]*yts; // convertion to m/yr
     2110     }
    21632111   }
    21642112
    2165    /*Get material parameters :*/
    2166    rho_ice=matpar->GetRhoIce();
    2167    rho_water=matpar->GetRhoFreshwater();
    2168    
    2169    sconv=(rho_water/rho_ice)/12.; //rhow_rain/rhoi / 12 months
    2170      
    2171      /*PDD constant*/
    2172    siglim = 2.5*signorm;
    2173    siglimc = 2.5*signormc;
    2174    siglim0 =  siglim/DT + 0.5;
    2175    siglim0c =  siglimc/DT + 0.5;
    2176    PDup = siglimc+PDCUT;
    2177    
    2178    // seasonal loop
    2179    for (iqj = 0; iqj < 12; iqj++){
    2180      imonth =  ismon[iqj];
    2181      for (i = 0; i < NUMVERTICES; i++){
    2182        st=(s[i]-s0t[i])/1000;
    2183        tstar = t[i][imonth] - lapser *max(st,sealev);
    2184        Tsurf[i] = tstar*deltm+Tsurf[i];       
    2185        
    2186        /*********compute PD ****************/
    2187        if (tstar < PDup){
    2188          pd = 1;
    2189          if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
    2190        else {
    2191          pd = 0;}
    2192        
    2193        /******exp des/elev precip reduction*******/
    2194        sp=(s[i]-s0p[i])/1000; // deselev effect is wrt chng in topo
    2195        if (sp>0.0){q = exp(-desfac*sp);}
    2196        else {q = 1.0;}
    2197        
    2198        qmt[i]= qmt[i] + prec[i][imonth]*sconv;  //*sconv to convert in m of ice equivalent per month
    2199        qmpt= q*prec[i][imonth]*sconv;           
    2200        qmp[i]= qmp[i] + qmpt;
    2201        qm[i]= qm[i] + qmpt*pd;
    2202        
    2203        /*********compute PDD************/
    2204        // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
    2205        // gaussian=T_m, so ndd=-(Tsurf-pdd)
    2206        if (iqj>6 &&  iqj<10){ Tsum[i]=Tsum[i]+tstar;}
    2207        if (tstar >= siglim) {pdd[i] = pdd[i] + tstar*deltm;}
    2208        else if (tstar> -siglim){
    2209          pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
    2210          pdd[i] = pdd[i] + pddsig*deltm;
    2211          frzndd[i] = frzndd[i] - (tstar-pddsig)*deltm;}
    2212        else{frzndd[i] = frzndd[i] - tstar*deltm; }
    2213      }
    2214    } // end of seasonal loop
    2215    
    2216    //******************************************************************
    2217    for(i=0;i<NUMVERTICES;i++){
    2218      saccu[i] = qm[i];
    2219      prect = qmp[i];     // total precipitation during 1 year taking into account des. ef.
    2220      Tsum[i]=Tsum[i]/3;
    2221      
    2222      /***** determine PDD factors *****/
    2223      if(Tsum[i]< -1.) {
    2224        snwmf=2.65*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
    2225        smf=17.22*0.001;    //  ablation factor for ice per pdd (Braithwaite 1995 from tarasov 2002)
    2226      }
    2227      else if(Tsum[i]< 10){
    2228        snwmf = (0.15*Tsum[i] + 2.8)*0.001;
    2229        smf = (0.0067*pow((10.-Tsum[i]),3) + 8.3)*0.001;
    2230      }
    2231      else{
    2232        snwmf=4.3*0.001;
    2233        smf=8.3*0.001;
    2234      }
    2235      snwmf=0.95*snwmf;
    2236      smf=0.95*smf;
    2237      
    2238      /*****  compute PDD ablation and refreezing *****/
    2239      pddt = pdd[i] *365;
    2240      snwm = snwmf*pddt;       // snow that could have been melted in a year
    2241      hmx2 = min(h[i],dfrz);   // refreeze active layer max depth: dfrz
    2242      
    2243      if(snwm < saccu[i]) {
    2244        water=prect-saccu[i] + snwm; //water=rain + snowmelt
    2245        //     l 2.2= capillary factor
    2246        //     Should refreezing be controlled by frzndd or by mean annual Tsurf?
    2247        //     dCovrLm concept is of warming of active layer (thickness =d@=1-
    2248        //     >2m)
    2249        //     problem with water seepage into ice: should be sealed after
    2250        //     refreezing
    2251        //     so everything needs to be predicated on 1 year scale, except for
    2252        //     thermal
    2253        //     conductivity through ice
    2254        //     also, need to account that melt season has low accum, so what's
    2255        //     going to
    2256        //     hold the meltwater around for refreezing? And melt-time will have
    2257        //     low seasonal frzndd
    2258        
    2259        //      Superimposed ice :  Pfeffer et al. 1991, Tarasov 2002
    2260        
    2261        supice= min(hmx2*CovrLm*frzndd[i]+2.2*(saccu[i]-snwm), water); // superimposed ice
    2262        supcap=min(2.2*(saccu[i]-snwm),water);
    2263        runoff=snwm - supice;  //meltwater only, does not include rain
    2264      }
    2265      else {  //all snow melted
    2266        supice= min(hmx2*CovrLm*frzndd[i], prect );
    2267        runoff= saccu[i] + smf*(pddt-saccu[i]/snwmf) - supice;
    2268        supcap=0;
    2269      }
    2270      //     pdd melting doesn't cool Tsurf, so ndd refreeze shouldn't warm it
    2271      //     except pdd melt heat source is atmosphere, while refreeze is
    2272      //     ground/ice stored interim
    2273      //     assume pdd=ndd, then melt should equal refreeze and Tsurf should=0
    2274      //     assume ndd=2*pdd, then all supice is refrozen, but Tsurf should be
    2275      //     <0
    2276      //     assume ndd>pdd, little melt => little supice
    2277      //     bottom line: compare for Tsurf<0 : supice and no supice case,
    2278      //     expect Tsurf difference
    2279      //     except some of cooling flux comes from atmosphere//
    2280      //     1 dm supice should not raise Tsurf by 1/dCovrLm = 16.675C
    2281      //     does supice make sense when H< 0.1m? then d=thermoactive ice layer ////
    2282      //     < 0.1
    2283      
    2284      //     make more sense to just use residual pdd-ndd except that pdd
    2285      //     residual not clear yet
    2286      //     frzndd should not be used up by refreezing in snow, so stick in
    2287      //     supcap.
    2288      diffndd=0;
    2289      if (frzndd[i]>0) {
    2290        diffndd=fsupndd*min((supice-supcap)/dCovrLm ,frzndd[i]);
    2291        frzndd[i]=frzndd[i]-diffndd;
    2292      }
    2293      if(runoff<0){
    2294        saccu[i]= saccu[i] -runoff;
    2295        smelt[i] = 0;
    2296        precrunoff[i]=prect-saccu[i];
    2297        //here assume pdd residual is 0, =>
    2298        Tsurf[i]= max(Tsurf[i],-frzndd[i]);
    2299      }
    2300      else {
    2301        smelt[i] = runoff;
    2302        precrunoff[i]=prect-max(0.,supice)-saccu[i];}
    2303      //here really need pdd balance, try 0.5 fudge factor?
    2304      //at least runoff>0 => it's fairly warm, so Tsurf is !<<0,
    2305      //yet from site plots, can be ice free with Tsurf=-5.5C
    2306      if(Tsurf[i]<0) {
    2307        Tsurf[i]= min(Tsurf[i]+fsupT*diffndd , 0.);}
    2308      
    2309      agd[i] = -smelt[i]+saccu[i];
    2310      pddtj[i]=pddt;
    2311      
    2312      /*Update inputs*/   
    2313      this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0])); ////////verifier le nom
    2314      // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    2315      this->inputs->AddInput(new TriaP1Input(ThermalSpctemperatureEnum,&Tsurf[0]));
    2316      
    2317    }       //end of the for loop over the vertices
     2113  /*Recover info at the vertices: */
     2114  GetInputListOnVertices(&h[0],ThicknessEnum);
     2115  GetInputListOnVertices(&s[0],SurfaceEnum);
     2116
     2117  /*Get material parameters :*/
     2118  rho_ice=matpar->GetRhoIce();
     2119  rho_water=matpar->GetRhoFreshwater();
     2120
     2121   /*measure the surface mass balance*/
     2122   for (int iv = 0; iv<NUMVERTICES; iv++){
     2123     agd[iv]=PddSurfaceMassBlance(&monthlytemperatures[iv][0], &monthlyprec[iv][0], pdds, pds, signorm, yts, h[iv], s[iv], rho_ice, rho_water);
     2124   }
     2125
     2126   /*Update inputs*/   
     2127   this->inputs->AddInput(new TriaP1Input(SurfaceforcingsMassBalanceEnum,&agd[0]));
     2128   // this->inputs->AddInput(new TriaVertexInput(ThermalSpcTemperatureEnum,&Tsurf[0]));
    23182129}
    23192130/*}}}*/

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

    r12677 r12704  
    8989                void   GetSolutionFromInputs(Vector* solution);
    9090                void   GetVectorFromInputs(Vector* vector, int name_enum);
    91                 void   GetVectorFromResults(Vector* vector,int offset,int interp);
     91                void   GetVectorFromResults(Vector* vector,int offset,int enum_in,int interp);
    9292                void   InputArtificialNoise(int enum_type,IssmDouble min, IssmDouble max);
    9393                bool   InputConvergence(IssmDouble* eps, int* enums,int num_enums,int* criterionenums,IssmDouble* criterionvalues,int num_criterionenums);

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

    r12554 r12704  
    4949                void GetInputValue(IssmDouble* pvalue);
    5050                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     51                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5455                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    5454                void GetInputValue(IssmDouble* pvalue);
    5555                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     56                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
    5657                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    57                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
     58                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5859                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5960                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    4444                void Configure(Parameters* parameters);
    4545                /*}}}*/
    46                 /*numeriics: {{{*/
     46                /*numerics: {{{*/
    4747                void GetInputValue(bool* pvalue){_error2_("not implemented yet");};
    4848                void GetInputValue(int* pvalue){_error2_("not implemented yet");};
    4949                void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");};
    5050                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error2_("not implemented yet");};
     51                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index);
    5455                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    4848                void GetInputValue(IssmDouble* pvalue);
    4949                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     50                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
    5051                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    51                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
     52                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5253                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    2121               
    2222                virtual        ~Input(){};
    23                 /*Virtual functions:{{{*/
     23
    2424                virtual int  InstanceEnum()=0;
    2525                virtual void GetInputValue(bool* pvalue)=0;
    … …  
    2727                virtual void GetInputValue(IssmDouble* pvalue)=0;
    2828                virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss)=0;
     29                virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0;
    2930                virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time)=0;
    30                 virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss)=0;
     31                virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time)=0;
    3132                virtual void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index)=0;
    3233                virtual void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index)=0;
    … …  
    6566                virtual Input* PointwiseMin(Input* inputmin)=0;
    6667                virtual ElementResult* SpawnResult(int step, IssmDouble time)=0;
    67 
    68                 /*}}}*/
    69 
    7068};
    7169#endif

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

    r12550 r12704  
    4949                void GetInputValue(IssmDouble* pvalue);
    5050                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     51                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5455                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    4949                void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");};
    5050                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss){_error2_("not implemented yet");};
     51                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss);
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5455                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    172172}
    173173/*}}}*/
     174/*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){{{*/
     175void TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){
     176        IssmDouble time;
     177
     178        /*First, recover current time from parameters: */
     179        this->parameters->FindParam(&time,TimeEnum);
     180
     181        /*Retrieve interpolated values for this time step: */
     182        Input* input=GetTimeInput(time);
     183
     184        /*Call input function*/
     185        input->GetInputValue(pvalue,gauss);
     186
     187        delete input;
     188}
     189/*}}}*/
    174190/*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){{{*/
    175191void TransientInput::GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){
     192
     193        /*Retrieve interpolated values for this time step: */
     194        Input* input=GetTimeInput(time);
     195
     196        /*Call input function*/
     197        input->GetInputValue(pvalue,gauss);
     198
     199        delete input;
     200}
     201/*}}}*/
     202/*FUNCTION TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){{{*/
     203void TransientInput::GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){
    176204
    177205        /*Retrieve interpolated values for this time step: */

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

    r12550 r12704  
    5151                void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");};
    5252                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time);
    54                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
     55                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5556                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5657                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss ,int index){_error2_("not implemented yet");};

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

    r12550 r12704  
    4949                void GetInputValue(IssmDouble* pvalue){_error2_("not implemented yet");}
    5050                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss);
     51                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
    5152                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss,IssmDouble time){_error2_("not implemented yet");};
    52                 void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss){_error2_("not implemented yet");};
     53                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,IssmDouble time){_error2_("not implemented yet");};
    5354                void GetInputValue(IssmDouble* pvalue,GaussTria* gauss ,int index){_error2_("not implemented yet");};
    5455                void GetInputValue(IssmDouble* pvalue,GaussPenta* gauss,int index){_error2_("not implemented yet");};

  • issm/trunk-jpl/src/c/shared/Elements/elements.h

    r12581 r12704  
    1313IssmDouble Paterson(IssmDouble temperature);
    1414IssmDouble Arrhenius(IssmDouble temperature,IssmDouble depth,IssmDouble n);
     15IssmDouble PddSurfaceMassBlance(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec, IssmDouble* pdds, IssmDouble* pds, IssmDouble signorm, IssmDouble yts, IssmDouble h, IssmDouble s, IssmDouble rho_ice, IssmDouble rho_water);
    1516void   GetVerticesCoordinates(IssmDouble* xyz,  Node** nodes, int numvertices);
    1617int    GetNumberOfDofs( Node** nodes,int numnodes,int setenum,int approximation_enum);

  • issm/trunk-jpl/src/c/shared/Exp/DomainOutlineRead.cpp

    r12494 r12704  
    77
    88#include <stdio.h>
     9#include <cstring>
    910#include "../Alloc/alloc.h"
    1011#include "../../include/include.h"

  • issm/trunk-jpl/src/m/model/extrude.m

    r11742 r12704  
    143143md.surfaceforcings.precipitation=project3d(md,'vector',md.surfaceforcings.precipitation,'type','node');
    144144md.balancethickness.thickening_rate=project3d(md,'vector',md.balancethickness.thickening_rate,'type','node');
     145md.surfaceforcings.monthlytemperatures=project3d(md,'vector',md.surfaceforcings.monthlytemperatures,'type','node');
    145146
    146147%results
Note: See TracChangeset for help on using the changeset viewer.