Changeset 19554


Ignore:
Timestamp:
09/16/15 17:15:08 (10 years ago)
Author:
Eric.Larour
Message:

CHG: finished first prototype implementation of the GEMB model inside the Gembx module of the SurfaceMassBalancex module,
and inside the SMBgemb.m class, as well as the SMBGemb method of the Element.cpp class.
Had to introduce the concept also of a DoubleArrayInput to keep track of the vertical gridded datasets within each element.

Location:
issm/trunk-jpl
Files:
24 added
30 edited

Legend:

Unmodified
Added
Removed

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

    r19528 r19554  
    5858                                        ./classes/Inputs/IntInput.cpp\
    5959                                        ./classes/Inputs/DoubleInput.cpp\
     60                                        ./classes/Inputs/DoubleArrayInput.cpp\
    6061                                        ./classes/Inputs/DatasetInput.cpp\
    6162                                        ./classes/Materials/Materials.cpp\
    … …  
    171172                                        ./modules/SpcNodesx/SpcNodesx.cpp\
    172173                                        ./modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp\
     174                                        ./modules/SurfaceMassBalancex/Gembx.cpp\
    173175                                        ./modules/Reducevectorgtofx/Reducevectorgtofx.cpp\
    174176                                        ./modules/Reduceloadx/Reduceloadx.cpp\

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r19528 r19554  
    2323        bool   isdelta18o,ismungsm,isd18opd;
    2424       
    25         /*Figure out smb model: */
    26         iomodel->Constant(&smb_model,SmbEnum);
    27 
    2825        /*Update elements: */
    2926        int counter=0;
    … …  
    3633        }
    3734       
     35        /*Figure out smb model: */
     36        iomodel->Constant(&smb_model,SmbEnum);
     37
     38                       
     39        Input* bof=NULL;
     40        Element* element=NULL;
    3841        switch(smb_model){
    3942                case SMBforcingEnum:
    4043                        iomodel->FetchDataToInput(elements,SmbMassBalanceEnum,0.);
     44                        break;
     45                case SMBgembEnum:
     46                        iomodel->FetchDataToInput(elements,SmbTaEnum);
     47                        iomodel->FetchDataToInput(elements,SmbVEnum);
     48                        iomodel->FetchDataToInput(elements,SmbDswrfEnum);
     49                        iomodel->FetchDataToInput(elements,SmbDlwrfEnum);
     50                        iomodel->FetchDataToInput(elements,SmbPEnum);
     51                        iomodel->FetchDataToInput(elements,SmbEAirEnum);
     52                        iomodel->FetchDataToInput(elements,SmbPAirEnum);
     53                        iomodel->FetchDataToInput(elements,SmbZTopEnum);
     54                        iomodel->FetchDataToInput(elements,SmbDzTopEnum);
     55                        iomodel->FetchDataToInput(elements,SmbDzMinEnum);
     56                        iomodel->FetchDataToInput(elements,SmbZYEnum);
     57                        iomodel->FetchDataToInput(elements,SmbZMaxEnum);
     58                        iomodel->FetchDataToInput(elements,SmbZMinEnum);
     59                        iomodel->FetchDataToInput(elements,SmbTmeanEnum);
     60                        iomodel->FetchDataToInput(elements,SmbCEnum);
     61                        iomodel->FetchDataToInput(elements,SmbTzEnum);
     62                        iomodel->FetchDataToInput(elements,SmbVzEnum);
    4163                        break;
    4264                case SMBpddEnum:
    … …  
    92114                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    93115        }
     116       
     117       
    94118
    95119}/*}}}*/
    … …  
    111135                case SMBforcingEnum:
    112136                        /*Nothing to add to parameters*/
     137                        break;
     138                case SMBgembEnum:
     139                        parameters->AddObject(iomodel->CopyConstantObject(SmbAIdxEnum));
     140                        parameters->AddObject(iomodel->CopyConstantObject(SmbSwIdxEnum));
     141                        parameters->AddObject(iomodel->CopyConstantObject(SmbDenIdxEnum));
     142                        parameters->AddObject(iomodel->CopyConstantObject(SmbOutputFreqEnum));
     143                        parameters->AddObject(iomodel->CopyConstantObject(SmbCldFracEnum));
     144                        parameters->AddObject(iomodel->CopyConstantObject(SmbT0wetEnum));
     145                        parameters->AddObject(iomodel->CopyConstantObject(SmbT0dryEnum));
     146                        parameters->AddObject(iomodel->CopyConstantObject(SmbKEnum));
     147                        parameters->AddObject(iomodel->CopyConstantObject(SmbASnowEnum));
     148                        parameters->AddObject(iomodel->CopyConstantObject(SmbAIceEnum));
     149                        parameters->AddObject(iomodel->CopyConstantObject(SmbDtEnum));
     150                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsgraingrowthEnum));
     151                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsalbedoEnum));
     152                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsshortwaveEnum));
     153                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsthermalEnum));
     154                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsaccumulationEnum));
     155                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsmeltEnum));
     156                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsdensificationEnum));
     157                        parameters->AddObject(iomodel->CopyConstantObject(SmbIsturbulentfluxEnum));
    113158                        break;
    114159                case SMBpddEnum:
    … …  
    196241                        /*Nothing to be done*/
    197242                        break;
     243                case SMBgembEnum:
     244                        Gembx(femmodel);
     245                        break;
    198246                case SMBpddEnum:
    199247                        bool isdelta18o,ismungsm;

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

    r19527 r19554  
    1414#include "../classes.h"
    1515#include "../../shared/shared.h"
     16#include "../../modules/SurfaceMassBalancex/SurfaceMassBalancex.h"
    1617/*}}}*/
    1718
    … …  
    13661367        }
    13671368        else if(vector_type==2){ //element vector
     1369
     1370                IssmDouble value;
     1371
    13681372                /*Are we in transient or static? */
    13691373                if(M==iomodel->numberofelements){
    … …  
    13791383                        else _error_("could not recognize nature of vector from code " << code);
    13801384                }
    1381                 else {
    1382                         _error_("transient element inputs not supported yet!");
    1383                 }
    1384         }
    1385         else{
    1386                 _error_("Cannot add input for vector type " << vector_type << " (not supported)");
    1387         }
    1388 }/*}}}*/
     1385                else if(M==iomodel->numberofelements+1){
     1386                        /*create transient input: */
     1387                        IssmDouble* times = xNew<IssmDouble>(N);
     1388                        for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
     1389                        TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1390                        TriaInput* bof=NULL;
     1391                        for(t=0;t<N;t++){
     1392                                value=vector[N*this->Sid()+t];
     1393                                switch(this->ObjectEnum()){
     1394                                        case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
     1395                                        case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
     1396                                        case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
     1397                                        default: _error_("Not implemented yet");
     1398                                }
     1399                        }
     1400                        this->inputs->AddInput(transientinput);
     1401                        xDelete<IssmDouble>(times);
     1402                }
     1403                else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     1404        }
     1405        else _error_("Cannot add input for vector type " << vector_type << " (not supported)");
     1406}
     1407/*}}}*/
    13891408void       Element::InputDuplicate(int original_enum,int new_enum){/*{{{*/
    13901409
    … …  
    18661885
    18671886}/*}}}*/
    1868 void       Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/
     1887void       Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int* parray_size, int output_enum){/*{{{*/
    18691888
    18701889        Input* input=this->inputs->GetInput(output_enum);
    … …  
    19601979        *pinterpolation   = input->GetResultInterpolation();
    19611980        *pnodesperelement = input->GetResultNumberOfNodes();
     1981        *parray_size      = input->GetResultArraySize();
    19621982}/*}}}*/
    19631983void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/
    … …  
    19671987
    19681988        input->ResultToPatch(values,nodesperelement,this->Sid());
     1989
     1990} /*}}}*/
     1991void       Element::ResultToMatrix(IssmDouble* values,int ncols,int output_enum){/*{{{*/
     1992
     1993        Input* input=this->inputs->GetInput(output_enum);
     1994        if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
     1995
     1996        input->ResultToMatrix(values,ncols,this->Sid());
    19691997
    19701998} /*}}}*/
    … …  
    20932121}
    20942122/*}}}*/
     2123void       Element::SmbGemb(){/*{{{*/
     2124
     2125        /*Intermediary variables: {{{*/
     2126        bool       isinitialized=false;
     2127        IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin;
     2128        IssmDouble Tmean;
     2129        IssmDouble C;
     2130        IssmDouble Tz,Vz;
     2131        IssmDouble rho_ice, aSnow,aIce;
     2132        IssmDouble time,dt;
     2133        IssmDouble t,smb_dt;
     2134        IssmDouble Ta,V,dlw,dsw,P,eAir,pAir;
     2135        int        aIdx=0;
     2136        int        denIdx=0;
     2137        int        swIdx=0;
     2138        IssmDouble cldFrac,t0wet, t0dry, K;
     2139        IssmDouble comp1,comp2;
     2140        IssmDouble ulw;
     2141        IssmDouble netSW;
     2142        IssmDouble netLW;
     2143        IssmDouble lhf, shf, dayEC;
     2144        IssmDouble initMass;
     2145    IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd;
     2146    IssmDouble sumMass,dMass;
     2147        bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
     2148        /*}}}*/
     2149        /*Output variables:{{{ */
     2150        IssmDouble* dz=NULL;
     2151        IssmDouble* d = NULL;
     2152        IssmDouble* re = NULL;
     2153        IssmDouble* gdn = NULL;
     2154        IssmDouble* gsp = NULL;
     2155        IssmDouble  EC = 0;
     2156        IssmDouble* W = NULL;
     2157        IssmDouble* a = NULL;
     2158        IssmDouble* swf=NULL;
     2159        IssmDouble* T = NULL;
     2160        IssmDouble  T_bottom;
     2161        IssmDouble  M;
     2162        IssmDouble  R;
     2163        IssmDouble  mAdd;
     2164        int         m;
     2165        /*}}}*/
     2166
     2167        /*only compute SMB at the surface: */
     2168        if (!IsOnSurface()) return;
     2169
     2170        /*Retrieve material properties and parameters:{{{ */
     2171        rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     2172        parameters->FindParam(&aSnow,SmbASnowEnum);
     2173        parameters->FindParam(&aIce,SmbAIceEnum);
     2174        parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
     2175        parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
     2176        parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
     2177        parameters->FindParam(&aIdx,SmbAIdxEnum);
     2178        parameters->FindParam(&denIdx,SmbDenIdxEnum);
     2179        parameters->FindParam(&swIdx,SmbSwIdxEnum);
     2180        parameters->FindParam(&cldFrac,SmbCldFracEnum);
     2181        parameters->FindParam(&t0wet,SmbT0wetEnum);
     2182        parameters->FindParam(&t0dry,SmbT0dryEnum);
     2183        parameters->FindParam(&K,SmbKEnum);
     2184        parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum);
     2185        parameters->FindParam(&isalbedo,SmbIsalbedoEnum);
     2186        parameters->FindParam(&isshortwave,SmbIsshortwaveEnum);
     2187        parameters->FindParam(&isthermal,SmbIsthermalEnum);
     2188        parameters->FindParam(&isaccumulation,SmbIsaccumulationEnum);
     2189        parameters->FindParam(&ismelt,SmbIsmeltEnum);
     2190        parameters->FindParam(&isdensification,SmbIsdensificationEnum);
     2191        parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
     2192        /*}}}*/
     2193        /*Retrieve inputs: {{{*/
     2194        Input* zTop_input=this->GetInput(SmbZTopEnum); _assert_(zTop_input);
     2195        Input* dzTop_input=this->GetInput(SmbDzTopEnum); _assert_(dzTop_input);
     2196        Input* dzMin_input=this->GetInput(SmbDzMinEnum); _assert_(dzMin_input);
     2197        Input* zMax_input=this->GetInput(SmbZMaxEnum); _assert_(zMax_input);
     2198        Input* zMin_input=this->GetInput(SmbZMinEnum); _assert_(zMin_input);
     2199        Input* zY_input=this->GetInput(SmbZYEnum); _assert_(zY_input);
     2200        Input* Tmean_input=this->GetInput(SmbTmeanEnum); _assert_(Tmean_input);
     2201        Input* C_input=this->GetInput(SmbCEnum); _assert_(C_input);
     2202        Input* Tz_input=this->GetInput(SmbTzEnum); _assert_(Tz_input);
     2203        Input* Vz_input=this->GetInput(SmbVzEnum); _assert_(Vz_input);
     2204        Input* Ta_input=this->GetInput(SmbTaEnum); _assert_(Ta_input);
     2205        Input* V_input=this->GetInput(SmbVEnum); _assert_(V_input);
     2206        Input* Dlwr_input=this->GetInput(SmbDlwrfEnum); _assert_(Dlwr_input);
     2207        Input* Dswr_input=this->GetInput(SmbDswrfEnum); _assert_(Dswr_input);
     2208        Input* P_input=this->GetInput(SmbPEnum); _assert_(P_input);
     2209        Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input);
     2210        Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input);
     2211        Input* isinitialized_input=this->inputs->GetInput(SmbIsInitializedEnum);
     2212       
     2213        /*Retrieve input values:*/
     2214        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
     2215
     2216        zTop_input->GetInputValue(&zTop,gauss);
     2217        dzTop_input->GetInputValue(&dzTop,gauss);
     2218        dzMin_input->GetInputValue(&dzMin,gauss);
     2219        zMax_input->GetInputValue(&zMax,gauss);
     2220        zMin_input->GetInputValue(&zMin,gauss);
     2221        zY_input->GetInputValue(&zY,gauss);
     2222        Tmean_input->GetInputValue(&Tmean,gauss);
     2223        C_input->GetInputValue(&C,gauss);
     2224        Tz_input->GetInputValue(&Tz,gauss);
     2225        Vz_input->GetInputValue(&Vz,gauss);
     2226        /*}}}*/
     2227        /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
     2228        if(!isinitialized_input){
     2229
     2230                /*Initialize grid:*/
     2231                GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
     2232                //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
     2233                //dz[i] << "\n");
     2234
     2235                /*initialize profile variables:*/
     2236                d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice; //ice density
     2237                re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5;         //set grain size to old snow [mm]
     2238                gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0;         //set grain dentricity to old snow
     2239                gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0;         //set grain sphericity to old snow
     2240                EC = 0;                                                                //surface evaporation (-) condensation (+) [kg m-2]
     2241                W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0;             //set water content to zero [kg m-2]
     2242                a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //set albedo equal to fresh snow [fraction]
     2243                T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
     2244
     2245                //fixed lower temperatuer bounday condition - T is fixed
     2246                T_bottom=T[m-1];
     2247
     2248                /*Flag the initialization:*/
     2249                this->AddInput(new BoolInput(SmbIsInitializedEnum,true));
     2250        }
     2251        else{
     2252                /*Recover inputs: */
     2253                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum));
     2254                DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));
     2255                DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));
     2256                DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));
     2257                DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));
     2258                DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));
     2259                DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));
     2260                DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));
     2261                DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));
     2262               
     2263                /*Recover arrays: */
     2264                dz_input->GetValues(&dz,&m);
     2265                d_input->GetValues(&d,&m);
     2266                re_input->GetValues(&re,&m);
     2267                gdn_input->GetValues(&gdn,&m);
     2268                gsp_input->GetValues(&gsp,&m);
     2269                EC_input->GetInputValue(&EC);
     2270                W_input->GetValues(&W,&m);
     2271                a_input->GetValues(&a,&m);
     2272                T_input->GetValues(&T,&m);
     2273
     2274        } /*}}}*/
     2275
     2276        // determine initial mass [kg]
     2277        initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i];
     2278   
     2279    // initialize cumulative variables
     2280    sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
     2281
     2282        /*Start loop: */
     2283        for (t=time;t<=time+dt;t=t+smb_dt){
     2284
     2285                /*extract daily data:{{{*/
     2286                Ta_input->GetInputValue(&Ta,gauss,t);//screen level air temperature [K]
     2287                V_input->GetInputValue(&V,gauss,t);  //wind speed [m s-1]
     2288                Dlwr_input->GetInputValue(&dlw,gauss,t);   //downward longwave radiation flux [W m-2]
     2289                Dswr_input->GetInputValue(&dsw,gauss,t);   //downward shortwave radiation flux [W m-2]
     2290                P_input->GetInputValue(&P,gauss,t);        //precipitation [kg m-2]
     2291                eAir_input->GetInputValue(&eAir,gauss,t);  //screen level vapor pressure [Pa]
     2292                pAir_input->GetInputValue(&pAir,gauss,t);  // screen level air pressure [Pa]
     2293                //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
     2294                /*}}}*/
     2295
     2296                /*Snow grain metamorphism:*/
     2297                if(isgraingrowth)grainGrowth(re, gdn, gsp, T, dz, d, W, smb_dt, m, aIdx);
     2298
     2299                /*Snow, firn and ice albedo:*/
     2300                if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m);
     2301               
     2302                /*Distribution of absorbed short wave radation with depth:*/
     2303        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m);
     2304               
     2305                /*Thermal profile computation:*/
     2306        if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, Vz, Tz,m);     
     2307
     2308                /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
     2309                 * need to fix this in case all or more of cell evaporates */
     2310        dz[0] = dz[0] + EC / d[0];
     2311               
     2312                /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
     2313        if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow);
     2314
     2315                /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
     2316                 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
     2317                comp2 = cellsum(dz,m);
     2318                if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin);
     2319        comp2 = (comp2 - cellsum(dz,m));
     2320
     2321                /*Allow non-melt densification and determine compaction [m]*/
     2322        comp1 = cellsum(dz,m);
     2323        if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,m);
     2324        comp1 = (comp1 - cellsum(dz,m));
     2325               
     2326                /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
     2327                 * sub-time step in thermo equations*/
     2328        ulw = 5.67E-8 * pow(T[0],4.0);
     2329
     2330                /*Calculate net shortwave and longwave [W m-2]*/
     2331        netSW = cellsum(swf,m);
     2332        netLW = dlw - ulw;
     2333               
     2334                /*Calculate turbulent heat fluxes [W m-2]*/
     2335        if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz);
     2336               
     2337                /*Sum component mass changes [kg m-2]*/
     2338        sumMassAdd = mAdd + sumMassAdd;
     2339        sumM = M + sumM;
     2340        sumR = R + sumR;
     2341        sumW = cellsum(W,m);
     2342        sumP = P +  sumP;
     2343        sumEC = sumEC + EC;  // evap (-)/cond(+)
     2344
     2345                /*Calculate total system mass:*/
     2346        sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i];
     2347        dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
     2348        dMass = round(dMass * 100.0)/100.0;
     2349               
     2350                /*Check mass conservation:*/
     2351        if (dMass != 0.0) _error_("total system mass not conserved in MB function");
     2352               
     2353                /*Check bottom grid cell T is unchanged:*/
     2354        if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom");
     2355        }
     2356
     2357        /*Save generated inputs: */
     2358        this->AddInput(new DoubleArrayInput(SmbDzEnum,dz,m));
     2359        this->AddInput(new DoubleArrayInput(SmbDEnum,d,m));
     2360        this->AddInput(new DoubleArrayInput(SmbReEnum,re,m));
     2361        this->AddInput(new DoubleArrayInput(SmbGdnEnum,gdn,m));
     2362        this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m));
     2363        this->AddInput(new DoubleInput(SmbECEnum,EC));
     2364        this->AddInput(new DoubleArrayInput(SmbWEnum,W,m));
     2365        this->AddInput(new DoubleArrayInput(SmbAEnum,a,m));
     2366        this->AddInput(new DoubleArrayInput(SmbTEnum,T,m));
     2367
     2368        /*Free allocations:{{{*/
     2369        xDelete<IssmDouble>(dz);
     2370        xDelete<IssmDouble>(d);
     2371        xDelete<IssmDouble>(re);
     2372        xDelete<IssmDouble>(gdn);
     2373        xDelete<IssmDouble>(gsp);
     2374        xDelete<IssmDouble>(W);
     2375        xDelete<IssmDouble>(a);
     2376        xDelete<IssmDouble>(T);
     2377        /*}}}*/
     2378}
     2379/*}}}*/
    20952380void       Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
    20962381        /*Compute the 3d Strain Rate (6 components):

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

    r19518 r19554  
    130130                IssmDouble         PureIceEnthalpy(IssmDouble pressure);
    131131                void               ResetIceLevelset(void);
    132                 void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum);
     132                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);
    133133                void               ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum);
     134                void               ResultToMatrix(IssmDouble* values,int ncols,int output_enum);
    134135                void               ResultToVector(Vector<IssmDouble>* vector,int output_enum);
    135136                void               SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);
    136137                int                Sid();
     138                void               SmbGemb();
    137139                void               StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    138140                void               StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);

  • issm/trunk-jpl/src/c/classes/FemModel.cpp

    r19528 r19554  
    10401040
    10411041                                                /*Vector layout*/
    1042                                                 int interpolation,nodesperelement,size;
    1043                                                 int rank_interpolation=-1,rank_nodesperelement=-1;
     1042                                                int interpolation,nodesperelement,size,nlines,ncols,array_size;
     1043                                                int rank_interpolation=-1,rank_nodesperelement=-1,rank_arraysize=-1,max_rank_arraysize=0;
     1044                                                bool isarray=false;
    10441045
    10451046                                                /*Get interpolation (and compute input if necessary)*/
    10461047                                                for(int j=0;j<elements->Size();j++){
    10471048                                                        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
    1048                                                         element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,output_enum);
     1049                                                        element->ResultInterpolation(&rank_interpolation,&rank_nodesperelement,&rank_arraysize,output_enum);
     1050                                                        if (rank_arraysize>max_rank_arraysize)max_rank_arraysize=rank_arraysize;
    10491051                                                }
     1052                                                rank_arraysize=max_rank_arraysize;
    10501053
    10511054                                                /*Broadcast for cpus that do not have any elements*/
    10521055                                                ISSM_MPI_Reduce(&rank_interpolation,&interpolation,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    10531056                                                ISSM_MPI_Reduce(&rank_nodesperelement,&nodesperelement,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
     1057                                                ISSM_MPI_Reduce(&rank_arraysize,&array_size,1,ISSM_MPI_INT,ISSM_MPI_MAX,0,IssmComm::GetComm());
    10541058                                                ISSM_MPI_Bcast(&interpolation,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    10551059                                                ISSM_MPI_Bcast(&nodesperelement,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    1056 
     1060                                                ISSM_MPI_Bcast(&array_size,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1061                                               
    10571062                                                if(results_on_nodes){
    10581063
    … …  
    10801085                                                        /*Allocate vector depending on interpolation*/
    10811086                                                        switch(interpolation){
    1082                                                                 case P0Enum: size = this->elements->NumberOfElements(); break;
    1083                                                                 case P1Enum: size = this->vertices->NumberOfVertices(); break;
     1087                                                                case P0Enum: isarray = false; size = this->elements->NumberOfElements(); break;
     1088                                                                case P1Enum: isarray = false; size = this->vertices->NumberOfVertices(); break;
     1089                                                                case P0ArrayEnum: isarray = true; nlines = this->elements->NumberOfElements(); ncols= array_size; break;
    10841090                                                                default:     _error_("Interpolation "<<EnumToStringx(interpolation)<<" not supported yet");
    10851091
    10861092                                                        }
    1087                                                         Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
    1088 
    1089                                                         /*Fill in vector*/
    1090                                                         for(int j=0;j<elements->Size();j++){
    1091                                                                 Element* element=(Element*)elements->GetObjectByOffset(j);
    1092                                                                 element->ResultToVector(vector_result,output_enum);
     1093                                                        if (!isarray){
     1094                                                                Vector<IssmDouble> *vector_result = new Vector<IssmDouble>(size);
     1095
     1096                                                                /*Fill in vector*/
     1097                                                                for(int j=0;j<elements->Size();j++){
     1098                                                                        Element* element=(Element*)elements->GetObjectByOffset(j);
     1099                                                                        element->ResultToVector(vector_result,output_enum);
     1100                                                                }
     1101                                                                vector_result->Assemble();
     1102
     1103                                                                if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
    10931104                                                        }
    1094                                                         vector_result->Assemble();
    1095 
    1096                                                         if(save_results)results->AddResult(new GenericExternalResult<Vector<IssmDouble>*>(results->Size()+1,output_enum,vector_result,step,time));
     1105                                                        else{
     1106                                                                IssmDouble* values    = xNewZeroInit<IssmDouble>(nlines*ncols);
     1107                                                                IssmDouble* allvalues = xNew<IssmDouble>(nlines*ncols);
     1108                                                               
     1109                                                                /*Fill-in matrix*/
     1110                                                                for(int j=0;j<elements->Size();j++){
     1111                                                                        Element* element=xDynamicCast<Element*>(this->elements->GetObjectByOffset(j));
     1112                                                                        element->ResultToMatrix(values,ncols, output_enum);
     1113                                                                }
     1114                                                                /*Gather from all cpus*/
     1115                                                                ISSM_MPI_Allreduce((void*)values,(void*)allvalues,ncols*nlines,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
     1116                                                                xDelete<IssmDouble>(values);
     1117                                                               
     1118                                                                if(save_results)results->AddResult(new GenericExternalResult<IssmDouble*>(results->Size()+1,output_enum,allvalues,nlines,ncols,step,time));
     1119                                                                xDelete<IssmDouble>(allvalues);
     1120                                                        }
    10971121                                                }
    10981122                                                isvec = true;

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

    r19254 r19554  
    3838                int  GetResultInterpolation(void){return P0Enum;};
    3939                int  GetResultNumberOfNodes(void){return 1;};
     40                int  GetResultArraySize(void){return 1;};
    4041                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    4142                void Configure(Parameters* parameters);

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

    r19254 r19554  
    7979                int  GetResultInterpolation(void);
    8080                int  GetResultNumberOfNodes(void);
     81                int  GetResultArraySize(void){return 1;};
    8182                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    8283                void GetGradient(Vector<IssmDouble>* gradient_vec,int* doflist);

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

    r19254 r19554  
    7474                int GetResultInterpolation(void){_error_("not implemented yet");};
    7575                int GetResultNumberOfNodes(void){_error_("not implemented yet");};
     76                int GetResultArraySize(void){_error_("not implemented yet");};
    7677                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    7778                void GetGradient(Vector<IssmDouble>* gradient_vec,int* doflist){_error_("not implemented yet");};

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

    r19254 r19554  
    4141                int  GetResultInterpolation(void){return P0Enum;};
    4242                int  GetResultNumberOfNodes(void){return 1;};
     43                int  GetResultArraySize(void){return 1;};
    4344                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    4445                void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

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

    r18450 r19554  
    6161                virtual int  GetResultInterpolation(void)=0;
    6262                virtual int  GetResultNumberOfNodes(void)=0;
    63                 virtual void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
     63                virtual int  GetResultArraySize(void)=0;
     64                virtual void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
     65                virtual void ResultToMatrix(IssmDouble* values,int ncols,int sid){_error_("not supported yet");};
    6466};
    6567#endif

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

    r19254 r19554  
    4242                int  GetResultInterpolation(void){return P0Enum;};
    4343                int  GetResultNumberOfNodes(void){return 1;};
     44                int  GetResultArraySize(void){return 1;};
    4445                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    4546                void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/PentaInput.h

    r19254 r19554  
    4343                int  GetResultInterpolation(void);
    4444                int  GetResultNumberOfNodes(void);
     45                int  GetResultArraySize(void){return 1;};
    4546                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid);
    4647                void AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/SegInput.h

    r19254 r19554  
    4343                int  GetResultInterpolation(void){_error_("not implemented");};
    4444                int  GetResultNumberOfNodes(void){_error_("not implemented");};
     45                int  GetResultArraySize(void){_error_("not implemented");};
    4546                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    4647                void   AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

  • issm/trunk-jpl/src/c/classes/Inputs/TetraInput.h

    r19254 r19554  
    4343                int    GetResultInterpolation(void);
    4444                int    GetResultNumberOfNodes(void);
     45                int    GetResultArraySize(void){return 1;};
    4546                void   ResultToPatch(IssmDouble* values,int nodesperelement,int sid);
    4647                void   AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

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

    r19254 r19554  
    357357}
    358358/*}}}*/
     359int  TransientInput::GetResultArraySize(void){/*{{{*/
     360
     361        return 1;
     362}
     363/*}}}*/
    359364void TransientInput::Extrude(int start){/*{{{*/
    360365

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

    r19254 r19554  
    4848                int  GetResultInterpolation(void);
    4949                int  GetResultNumberOfNodes(void);
     50                int  GetResultArraySize(void);
    5051                void ResultToPatch(IssmDouble* values,int nodesperelement,int sid){_error_("not supported yet");};
    5152                void Configure(Parameters* parameters);

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.cpp

    r19254 r19554  
    141141}
    142142/*}}}*/
     143int  TriaInput::GetResultArraySize(void){/*{{{*/
     144
     145        return 1;
     146
     147}
     148/*}}}*/
    143149void TriaInput::ResultToPatch(IssmDouble* values,int nodesperelement,int sid){/*{{{*/
    144150

  • issm/trunk-jpl/src/c/classes/Inputs/TriaInput.h

    r19254 r19554  
    4343                int    GetResultInterpolation(void);
    4444                int    GetResultNumberOfNodes(void);
     45                int    GetResultArraySize(void);
    4546                void   ResultToPatch(IssmDouble* values,int nodesperelement,int sid);
    4647                void   AddTimeValues(IssmDouble* values,int step,IssmDouble time){_error_("not supported yet");};

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

    r19527 r19554  
    3939        dpermil=0;
    4040
     41        albedo_snow=0;
     42        albedo_ice=0;
     43
    4144        sediment_compressibility=0;
    4245        sediment_porosity=0;
    … …  
    9598                                case SMBforcingEnum:
    9699                                        /*Nothing to add*/
     100                                        break;
     101                                case SMBgembEnum:
     102                                        iomodel->Constant(&this->albedo_ice,SmbAIceEnum);
     103                                        iomodel->Constant(&this->albedo_snow,SmbASnowEnum);
    97104                                        break;
    98105                                case SMBpddEnum:

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

    r19309 r19554  
    3535                IssmDouble  rlapslgm;
    3636                IssmDouble  dpermil;
     37
     38                /*albedo: */
     39                IssmDouble albedo_ice;
     40                IssmDouble albedo_snow;
    3741
    3842                /*hydrology Dual Porous Continuum: */   

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

    r19087 r19554  
    6262#include "./Inputs/BoolInput.h"
    6363#include "./Inputs/DoubleInput.h"
     64#include "./Inputs/DoubleArrayInput.h"
    6465#include "./Inputs/IntInput.h"
    6566#include "./Inputs/TetraInput.h"

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

    r19172 r19554  
    1717void SmbHenningx(FemModel* femmodel);
    1818void SmbComponentsx(FemModel* femmodel);
    19 void SmbMeltComponentsx(FemModel* femmodel);
     19void SmbMeltComponentsx(FemModel* femmodel); 
    2020
     21/*GEMB: */
     22void       Gembx(FemModel* femmodel);
     23void       GembgridInitialize(IssmDouble** pdz, int* psize, IssmDouble zTop, IssmDouble dzTop, IssmDouble zMax, IssmDouble zY);
     24IssmDouble Marbouty(IssmDouble T, IssmDouble d, IssmDouble dT);
     25void grainGrowth(IssmDouble* pre, IssmDouble* pgdn, IssmDouble* pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx);
     26void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m);
     27void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m);
     28void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz);
     29void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow);
     30void melt(IssmDouble* pM, IssmDouble* pR, IssmDouble* pmAdd, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin);
     31void densification(IssmDouble* d,IssmDouble* dz, IssmDouble* T, IssmDouble* re, int denIdx, IssmDouble C, IssmDouble dt, IssmDouble Tmean,int m);
     32void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz);
    2133#endif  /* _SurfaceMassBalancex_H*/

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

    r19528 r19554  
    354354        SmbNumRequestedOutputsEnum,
    355355        SmbRequestedOutputsEnum,
     356        SmbIsInitializedEnum,
    356357        /*SMBforcing*/
    357358        SMBforcingEnum,
    … …  
    369370        SmbCEnum,
    370371        SmbTzEnum,
    371         SmbVzEnum,
    372         SmbSpinUpEnum,
     372        SmbVzEnum,
     373        SmbDtEnum,
     374        SmbDzEnum,
    373375        SmbAIdxEnum,
    374376        SmbSwIdxEnum,
    … …  
    387389        SmbT0dryEnum,
    388390        SmbKEnum,
     391        SmbDEnum,
     392        SmbReEnum,
     393        SmbGdnEnum,
     394        SmbGspEnum,
     395        SmbECEnum,
     396        SmbWEnum,
     397        SmbAEnum,
     398        SmbTEnum,
     399        SmbIsgraingrowthEnum,
     400        SmbIsalbedoEnum,
     401        SmbIsshortwaveEnum,
     402        SmbIsthermalEnum,
     403        SmbIsaccumulationEnum,
     404        SmbIsmeltEnum,
     405        SmbIsdensificationEnum,
     406        SmbIsturbulentfluxEnum,
    389407        /*SMBpdd*/
    390408        SMBpddEnum,     
    … …  
    520538        DatasetInputEnum,
    521539        DoubleInputEnum,
     540        DoubleArrayInputEnum,
    522541        DataSetParamEnum,
    523542        DoubleMatArrayParamEnum,
    … …  
    696715        /*Element Interpolations{{{*/
    697716        P0Enum,
     717        P0ArrayEnum,
    698718        P1Enum,
    699719        P1DGEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r19528 r19554  
    358358                case SmbNumRequestedOutputsEnum : return "SmbNumRequestedOutputs";
    359359                case SmbRequestedOutputsEnum : return "SmbRequestedOutputs";
     360                case SmbIsInitializedEnum : return "SmbIsInitialized";
    360361                case SMBforcingEnum : return "SMBforcing";
    361362                case SmbMassBalanceEnum : return "SmbMassBalance";
    … …  
    372373                case SmbTzEnum : return "SmbTz";
    373374                case SmbVzEnum : return "SmbVz";
    374                 case SmbSpinUpEnum : return "SmbSpinUp";
     375                case SmbDtEnum : return "SmbDt";
     376                case SmbDzEnum : return "SmbDz";
    375377                case SmbAIdxEnum : return "SmbAIdx";
    376378                case SmbSwIdxEnum : return "SmbSwIdx";
    … …  
    389391                case SmbT0dryEnum : return "SmbT0dry";
    390392                case SmbKEnum : return "SmbK";
     393                case SmbDEnum : return "SmbD";
     394                case SmbReEnum : return "SmbRe";
     395                case SmbGdnEnum : return "SmbGdn";
     396                case SmbGspEnum : return "SmbGsp";
     397                case SmbECEnum : return "SmbEC";
     398                case SmbWEnum : return "SmbW";
     399                case SmbAEnum : return "SmbA";
     400                case SmbTEnum : return "SmbT";
     401                case SmbIsgraingrowthEnum : return "SmbIsgraingrowth";
     402                case SmbIsalbedoEnum : return "SmbIsalbedo";
     403                case SmbIsshortwaveEnum : return "SmbIsshortwave";
     404                case SmbIsthermalEnum : return "SmbIsthermal";
     405                case SmbIsaccumulationEnum : return "SmbIsaccumulation";
     406                case SmbIsmeltEnum : return "SmbIsmelt";
     407                case SmbIsdensificationEnum : return "SmbIsdensification";
     408                case SmbIsturbulentfluxEnum : return "SmbIsturbulentflux";
    391409                case SMBpddEnum : return "SMBpdd";
    392410                case SmbDelta18oEnum : return "SmbDelta18o";
    … …  
    512530                case DatasetInputEnum : return "DatasetInput";
    513531                case DoubleInputEnum : return "DoubleInput";
     532                case DoubleArrayInputEnum : return "DoubleArrayInput";
    514533                case DataSetParamEnum : return "DataSetParam";
    515534                case DoubleMatArrayParamEnum : return "DoubleMatArrayParam";
    … …  
    680699                case GiaWEnum : return "GiaW";
    681700                case P0Enum : return "P0";
     701                case P0ArrayEnum : return "P0Array";
    682702                case P1Enum : return "P1";
    683703                case P1DGEnum : return "P1DG";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r19528 r19554  
    364364              else if (strcmp(name,"SmbNumRequestedOutputs")==0) return SmbNumRequestedOutputsEnum;
    365365              else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum;
     366              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
    366367              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;
    367368              else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum;
    … …  
    378379              else if (strcmp(name,"SmbTz")==0) return SmbTzEnum;
    379380              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    380               else if (strcmp(name,"SmbSpinUp")==0) return SmbSpinUpEnum;
     381              else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
     382              else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
    381383              else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
    382384              else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
    383               else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
    384               else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
     388              if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
     389              else if (strcmp(name,"SmbZTop")==0) return SmbZTopEnum;
     390              else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
    389391              else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
    390392              else if (strcmp(name,"SmbZY")==0) return SmbZYEnum;
    … …  
    398400              else if (strcmp(name,"SmbT0dry")==0) return SmbT0dryEnum;
    399401              else if (strcmp(name,"SmbK")==0) return SmbKEnum;
     402              else if (strcmp(name,"SmbD")==0) return SmbDEnum;
     403              else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
     404              else if (strcmp(name,"SmbGdn")==0) return SmbGdnEnum;
     405              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
     406              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
     407              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
     408              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
     409              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
     410              else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum;
     411              else if (strcmp(name,"SmbIsalbedo")==0) return SmbIsalbedoEnum;
     412              else if (strcmp(name,"SmbIsshortwave")==0) return SmbIsshortwaveEnum;
     413              else if (strcmp(name,"SmbIsthermal")==0) return SmbIsthermalEnum;
     414              else if (strcmp(name,"SmbIsaccumulation")==0) return SmbIsaccumulationEnum;
     415              else if (strcmp(name,"SmbIsmelt")==0) return SmbIsmeltEnum;
     416              else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum;
     417              else if (strcmp(name,"SmbIsturbulentflux")==0) return SmbIsturbulentfluxEnum;
    400418              else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum;
    401419              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    … …  
    488506              else if (strcmp(name,"MeshdeformationSolution")==0) return MeshdeformationSolutionEnum;
    489507              else if (strcmp(name,"MeshdeformationAnalysis")==0) return MeshdeformationAnalysisEnum;
    490               else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    491512              else if (strcmp(name,"LevelsetStabilization")==0) return LevelsetStabilizationEnum;
    492513              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    … …  
    506527              else if (strcmp(name,"DataSet")==0) return DataSetEnum;
    507528              else if (strcmp(name,"Constraints")==0) return ConstraintsEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"Loads")==0) return LoadsEnum;
     529              else if (strcmp(name,"Loads")==0) return LoadsEnum;
    512530              else if (strcmp(name,"Materials")==0) return MaterialsEnum;
    513531              else if (strcmp(name,"Nodes")==0) return NodesEnum;
    … …  
    524542              else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum;
    525543              else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
     544              else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum;
    526545              else if (strcmp(name,"DataSetParam")==0) return DataSetParamEnum;
    527546              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
    … …  
    610629              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    611630              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    612               else if (strcmp(name,"Misfit")==0) return MisfitEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Misfit")==0) return MisfitEnum;
    613635              else if (strcmp(name,"Pressure")==0) return PressureEnum;
    614636              else if (strcmp(name,"PressurePicard")==0) return PressurePicardEnum;
    … …  
    629651              else if (strcmp(name,"ThicknessAbsMisfit")==0) return ThicknessAbsMisfitEnum;
    630652              else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Vel")==0) return VelEnum;
     653              else if (strcmp(name,"Vel")==0) return VelEnum;
    635654              else if (strcmp(name,"Velocity")==0) return VelocityEnum;
    636655              else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
    … …  
    695714              else if (strcmp(name,"GiaW")==0) return GiaWEnum;
    696715              else if (strcmp(name,"P0")==0) return P0Enum;
     716              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
    697717              else if (strcmp(name,"P1")==0) return P1Enum;
    698718              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    … …  
    732752              else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum;
    733753              else if (strcmp(name,"Outputdefinition5")==0) return Outputdefinition5Enum;
    734               else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
    735758              else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
    736759              else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
    … …  
    752775              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    753776              else if (strcmp(name,"Outputdefinition25")==0) return Outputdefinition25Enum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
     777              else if (strcmp(name,"Outputdefinition26")==0) return Outputdefinition26Enum;
    758778              else if (strcmp(name,"Outputdefinition27")==0) return Outputdefinition27Enum;
    759779              else if (strcmp(name,"Outputdefinition28")==0) return Outputdefinition28Enum;
    … …  
    855875              else if (strcmp(name,"MaxVx")==0) return MaxVxEnum;
    856876              else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
    857               else if (strcmp(name,"MinVy")==0) return MinVyEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"MinVy")==0) return MinVyEnum;
    858881              else if (strcmp(name,"MaxVy")==0) return MaxVyEnum;
    859882              else if (strcmp(name,"MaxAbsVy")==0) return MaxAbsVyEnum;
    … …  
    875898              else if (strcmp(name,"None")==0) return NoneEnum;
    876899              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
     900              else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
    881901              else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
    882902              else if (strcmp(name,"SubelementMigration2")==0) return SubelementMigration2Enum;

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r19030 r19554  
    557557        for(int i=0;i<4;i++) X[i]=Ainv[i][0]*B[0] + Ainv[i][1]*B[1] + Ainv[i][2]*B[2] + Ainv[i][3]*B[3];
    558558}/*}}}*/
     559
     560void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){  /*{{{*/
     561
     562        IssmDouble* cell=NULL;
     563        IssmDouble* dummy=NULL;
     564
     565        /*recover pointer: */
     566        cell=*pcell;
     567
     568        /*reallocate:*/
     569        dummy=xNew<IssmDouble>(m+1);
     570
     571        /*copy data:*/
     572        for(int i=0;i<m;i++)dummy[i+1]=cell[i];
     573
     574        /*top or bottom layer:*/
     575        if(top) dummy[0]=newvalue;
     576        else dummy[m]=newvalue;
     577       
     578        /*delete and reassign: */
     579        xDelete<IssmDouble>(cell); cell=dummy;
     580
     581        /*assign output pointer:*/
     582        *pcell=cell;
     583} /*}}}*/
     584IssmDouble  cellsum(IssmDouble* cell, int m){ /*{{{*/
     585
     586        IssmDouble sum=0;
     587
     588        for(int i=0;i<m;i++)sum+=cell[i];
     589
     590        return sum;
     591} /*}}}*/
     592void celldelete(IssmDouble** pcell, int m, int* indices, int nind){ /*{{{*/
     593
     594        /*input: */
     595        IssmDouble* cell=*pcell;
     596       
     597        /*output: */
     598        IssmDouble* newcell=xNew<IssmDouble>(nind);
     599
     600        for(int i=0;i<nind;i++){
     601                newcell[i]=cell[indices[i]];
     602        }
     603       
     604        /*free allocation:*/
     605        xDelete<IssmDouble>(cell);
     606
     607        /*assign output pointers: */
     608        *pcell=newcell;
     609} /*}}}*/
     610void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale) { /*{{{*/
     611
     612        /*input: */
     613        IssmDouble* cell=*pcell;
     614       
     615        /*output: */
     616        IssmDouble* newcell=xNew<IssmDouble>(m+1);
     617
     618        for(int j=0;j<i;j++)newcell[j]=cell[j];
     619        newcell[i]=scale*cell[i];
     620        newcell[i+1]=scale* cell[i]/2;
     621        for(int j=i+2;j<m+1;j++)newcell[j]=cell[j-1];
     622       
     623        /*free allocation:*/
     624        xDelete<IssmDouble>(cell);
     625
     626        /*assign output pointers: */
     627        *pcell=newcell;
     628} /*}}}*/

  • issm/trunk-jpl/src/c/shared/Matrix/matrix.h

    r19030 r19554  
    2525void Matrix4x4Determinant(IssmDouble* Adet,IssmDouble* A);
    2626void Matrix4x4Solve(IssmDouble* X,IssmDouble* A,IssmDouble* B);
    27 
     27 
     28void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m);
     29IssmDouble  cellsum(IssmDouble* cell, int m);
     30void celldelete(IssmDouble** pcell, int m, int* indices, int nind);
     31void cellsplit(IssmDouble** pcell, int m, int i,IssmDouble scale) ;
    2832#endif //ifndef _MATRIXUTILS_H_

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

    r19528 r19554  
    1212                %from an automatic weather stations (AWS). Each property is therefore a matrix, of size (numberofvertices x number
    1313                %of time steps. )
     14
     15                %solution choices
     16                isgraingrowth;
     17                isalbedo;
     18                isshortwave;
     19                isthermal;
     20                isaccumulation;
     21                ismelt;
     22                isdensification;
     23                isturbulentflux;
    1424
    1525                %inputs:
    … …  
    2838
    2939                %settings:
    30                 spinUp = NaN; %number of cycles of met data run before output is calculated (default is 0)
    3140                aIdx   = NaN; %method for calculating albedo and subsurface absorption (default is 1)
    3241                              % 1: effective grain radius [Gardner & Sharp, 2009]
    … …  
    97106                function self = setdefaultparameters(self,mesh,geometry) % {{{
    98107
    99                
    100                 self.spinUp=0;
     108                self.isgraingrowth=1;
     109                self.isalbedo=1;
     110                self.isshortwave=1;
     111                self.isthermal=1;
     112                self.isaccumulation=1;
     113                self.ismelt=1;
     114                self.isdensification=1;
     115                self.isturbulentflux=1;
     116       
    101117                self.aIdx = 1;
    102118                self.swIdx = 1;
    … …  
    123139function md = checkconsistency(self,md,solution,analyses) % {{{
    124140
     141
     142                md = checkfield(md,'fieldname','smb.isgraingrowth','values',[0 1]);
     143                md = checkfield(md,'fieldname','smb.isalbedo','values',[0 1]);
     144                md = checkfield(md,'fieldname','smb.isshortwave','values',[0 1]);
     145                md = checkfield(md,'fieldname','smb.isthermal','values',[0 1]);
     146                md = checkfield(md,'fieldname','smb.isaccumulation','values',[0 1]);
     147                md = checkfield(md,'fieldname','smb.ismelt','values',[0 1]);
     148                md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]);
     149                md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]);
     150               
    125151                md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'>',273-60,'<',273+60); %60 celsius max value
    126152                md = checkfield(md,'fieldname','smb.V','timeseries',1,'NaN',1,'>=',0,'<',45); %max 500 km/h
    … …  
    130156                md = checkfield(md,'fieldname','smb.eAir','timeseries',1,'NaN',1);
    131157               
    132                 md = checkfield(md,'fieldname','smb.Tmean','NaN',1,'>',273-60,'<',273+60); %60 celsius max value
    133                 md = checkfield(md,'fieldname','smb.C','NaN',1,'>=',0);
    134                 md = checkfield(md,'fieldname','smb.Tz','NaN',1,'>=',0,'<=',5000);
    135                 md = checkfield(md,'fieldname','smb.Vz','NaN',1,'>=',0,'<=',5000);
    136                
    137                 md = checkfield(md,'fieldname','smb.spinUp','NaN',1,'>=',0);
     158                md = checkfield(md,'fieldname','smb.Tmean','size',[md.mesh.numberofelements 1],'NaN',1,'>',273-60,'<',273+60); %60 celsius max value
     159                md = checkfield(md,'fieldname','smb.C','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0);
     160                md = checkfield(md,'fieldname','smb.Tz','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0,'<=',5000);
     161                md = checkfield(md,'fieldname','smb.Vz','size',[md.mesh.numberofelements 1],'NaN',1,'>=',0,'<=',5000);
     162               
    138163                md = checkfield(md,'fieldname','smb.aIdx','NaN',1,'values',[1,2,3,4]);
    139164                md = checkfield(md,'fieldname','smb.swIdx','NaN',1,'values',[0,1]);
    … …  
    170195                        disp(sprintf('   surface forcings for SMB GEMB model :'));
    171196                       
     197                        fielddisplay(self,'isgraingrowth','run grain growth module (default true)');
     198                        fielddisplay(self,'isalbedo','run albedo module (default true)');
     199                        fielddisplay(self,'isshortwave','run short wave module (default true)');
     200                        fielddisplay(self,'isthermal','run thermal module (default true)');
     201                        fielddisplay(self,'isaccumulation','run accumulation module (default true)');
     202                        fielddisplay(self,'ismelt','run melting  module (default true)');
     203                        fielddisplay(self,'isdensification','run densification module (default true)');
     204                        fielddisplay(self,'isturbulentflux','run turbulant heat fluxes module (default true)');
    172205                        fielddisplay(self,'Ta','2 m air temperature, in Kelvin');
    173206                        fielddisplay(self,'V','wind speed (m/s-1)');
    … …  
    188221                        fielddisplay(self,'zY','strech grid cells bellow top_z by a [top_dz * y ^ (cells bellow top_z)]');
    189222                        fielddisplay(self,'outputFreq','output frequency in days (default is monthly, 30)');
    190                         fielddisplay(self,'spinUp','number of cycles of met data run before output is calcualted (default is 0)');
    191223                        fielddisplay(self,'aIdx',{'method for calculating albedo and subsurface absorption (default is 1)',...
    192224                                                                        '1: effective grain radius [Gardner & Sharp, 2009]',...
    … …  
    198230                        case {1 2}
    199231                                fielddisplay(self,'aSnow','new snow albedo (0.64 - 0.89)');
    200                                 fielddisplay(self,'aIce','range 0.27-0.58 for old snow');
     232                                fielddisplay(self,'aIce','albedo of ice (0.27-0.58)');
    201233                        case 3
    202234                                fielddisplay(self,'cldFrac','average cloud amount');
    … …  
    224256                        WriteData(fid,'enum',SmbEnum(),'data',SMBgembEnum(),'format','Integer');
    225257                       
    226                         WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    227                         WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    228                         WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    229                         WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    230                         WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    231                         WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    232                         WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',1,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
    233                        
    234                         WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','Double','scale',1);
    235                         WriteData(fid,'object',self,'class','smb','fieldname','C','format','Double','scale',1);
    236                         WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','Double','scale',1);
    237                         WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','Double','scale',1);
    238                         WriteData(fid,'object',self,'class','smb','fieldname','spinUp','format','Integer','scale',1);
     258                        WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean');
     259                        WriteData(fid,'object',self,'class','smb','fieldname','isalbedo','format','Boolean');
     260                        WriteData(fid,'object',self,'class','smb','fieldname','isshortwave','format','Boolean');
     261                        WriteData(fid,'object',self,'class','smb','fieldname','isthermal','format','Boolean');
     262                        WriteData(fid,'object',self,'class','smb','fieldname','isaccumulation','format','Boolean');
     263                        WriteData(fid,'object',self,'class','smb','fieldname','ismelt','format','Boolean');
     264                        WriteData(fid,'object',self,'class','smb','fieldname','isdensification','format','Boolean');
     265                        WriteData(fid,'object',self,'class','smb','fieldname','isturbulentflux','format','Boolean');
     266                        WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean');
     267                        WriteData(fid,'object',self,'class','smb','fieldname','isgraingrowth','format','Boolean');
     268                       
     269                        WriteData(fid,'object',self,'class','smb','fieldname','Ta','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     270                        WriteData(fid,'object',self,'class','smb','fieldname','V','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     271                        WriteData(fid,'object',self,'class','smb','fieldname','dswrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     272                        WriteData(fid,'object',self,'class','smb','fieldname','dlwrf','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     273                        WriteData(fid,'object',self,'class','smb','fieldname','P','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     274                        WriteData(fid,'object',self,'class','smb','fieldname','eAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     275                        WriteData(fid,'object',self,'class','smb','fieldname','pAir','format','DoubleMat','mattype',2,'scale',1,'timeserieslength',md.mesh.numberofelements+1);
     276                       
     277                        WriteData(fid,'object',self,'class','smb','fieldname','Tmean','format','DoubleMat','mattype',2,'scale',1);
     278                        WriteData(fid,'object',self,'class','smb','fieldname','C','format','DoubleMat','mattype',2,'scale',1);
     279                        WriteData(fid,'object',self,'class','smb','fieldname','Tz','format','DoubleMat','mattype',2,'scale',1);
     280                        WriteData(fid,'object',self,'class','smb','fieldname','Vz','format','DoubleMat','mattype',2,'scale',1);
     281                        WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',2,'scale',1);
     282                        WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',2,'scale',1);
     283                        WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',2,'scale',1);
     284                        WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',2,'scale',1);
     285                        WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',2,'scale',1);
     286                        WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',2,'scale',1);
     287               
    239288                        WriteData(fid,'object',self,'class','smb','fieldname','aIdx','format','Integer','scale',1);
    240289                        WriteData(fid,'object',self,'class','smb','fieldname','swIdx','format','Integer','scale',1);
    241290                        WriteData(fid,'object',self,'class','smb','fieldname','denIdx','format','Integer','scale',1);
    242291
    243                         WriteData(fid,'object',self,'class','smb','fieldname','zTop','format','DoubleMat','mattype',1,'scale',1);
    244                         WriteData(fid,'object',self,'class','smb','fieldname','dzTop','format','DoubleMat','mattype',1,'scale',1);
    245                         WriteData(fid,'object',self,'class','smb','fieldname','dzMin','format','DoubleMat','mattype',1,'scale',1);
    246                         WriteData(fid,'object',self,'class','smb','fieldname','zY','format','DoubleMat','mattype',1,'scale',1);
    247                         WriteData(fid,'object',self,'class','smb','fieldname','zMax','format','DoubleMat','mattype',1,'scale',1);
    248                         WriteData(fid,'object',self,'class','smb','fieldname','zMin','format','DoubleMat','mattype',1,'scale',1);
    249292
    250293                        WriteData(fid,'object',self,'class','smb','fieldname','outputFreq','format','Double','scale',1);
    … …  
    255298                        WriteData(fid,'object',self,'class','smb','fieldname','t0dry','format','Double','scale',1);
    256299                        WriteData(fid,'object',self,'class','smb','fieldname','K','format','Double','scale',1);
     300
     301                        %figure out dt from forcings:
     302                        time=self.Ta(end,:); %assume all forcings are one the same time step
     303                        dtime=diff(time,1);
     304                        dt=min(dtime);
     305                        WriteData(fid,'data',dt,'enum',SmbDtEnum,'format','Double','scale',yts);
    257306                       
    258307                        %process requested outputs
    … …  
    264313                        end
    265314                        WriteData(fid,'data',outputs,'enum',SmbRequestedOutputsEnum,'format','StringArray');
    266 
    267315                end % }}}
    268316        end

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r19528 r19554  
    350350def SmbNumRequestedOutputsEnum(): return StringToEnum("SmbNumRequestedOutputs")[0]
    351351def SmbRequestedOutputsEnum(): return StringToEnum("SmbRequestedOutputs")[0]
     352def SmbIsInitializedEnum(): return StringToEnum("SmbIsInitialized")[0]
    352353def SMBforcingEnum(): return StringToEnum("SMBforcing")[0]
    353354def SmbMassBalanceEnum(): return StringToEnum("SmbMassBalance")[0]
    … …  
    364365def SmbTzEnum(): return StringToEnum("SmbTz")[0]
    365366def SmbVzEnum(): return StringToEnum("SmbVz")[0]
    366 def SmbSpinUpEnum(): return StringToEnum("SmbSpinUp")[0]
     367def SmbDtEnum(): return StringToEnum("SmbDt")[0]
     368def SmbDzEnum(): return StringToEnum("SmbDz")[0]
    367369def SmbAIdxEnum(): return StringToEnum("SmbAIdx")[0]
    368370def SmbSwIdxEnum(): return StringToEnum("SmbSwIdx")[0]
    … …  
    381383def SmbT0dryEnum(): return StringToEnum("SmbT0dry")[0]
    382384def SmbKEnum(): return StringToEnum("SmbK")[0]
     385def SmbDEnum(): return StringToEnum("SmbD")[0]
     386def SmbReEnum(): return StringToEnum("SmbRe")[0]
     387def SmbGdnEnum(): return StringToEnum("SmbGdn")[0]
     388def SmbGspEnum(): return StringToEnum("SmbGsp")[0]
     389def SmbECEnum(): return StringToEnum("SmbEC")[0]
     390def SmbWEnum(): return StringToEnum("SmbW")[0]
     391def SmbAEnum(): return StringToEnum("SmbA")[0]
     392def SmbTEnum(): return StringToEnum("SmbT")[0]
     393def SmbIsgraingrowthEnum(): return StringToEnum("SmbIsgraingrowth")[0]
     394def SmbIsalbedoEnum(): return StringToEnum("SmbIsalbedo")[0]
     395def SmbIsshortwaveEnum(): return StringToEnum("SmbIsshortwave")[0]
     396def SmbIsthermalEnum(): return StringToEnum("SmbIsthermal")[0]
     397def SmbIsaccumulationEnum(): return StringToEnum("SmbIsaccumulation")[0]
     398def SmbIsmeltEnum(): return StringToEnum("SmbIsmelt")[0]
     399def SmbIsdensificationEnum(): return StringToEnum("SmbIsdensification")[0]
     400def SmbIsturbulentfluxEnum(): return StringToEnum("SmbIsturbulentflux")[0]
    383401def SMBpddEnum(): return StringToEnum("SMBpdd")[0]
    384402def SmbDelta18oEnum(): return StringToEnum("SmbDelta18o")[0]
    … …  
    504522def DatasetInputEnum(): return StringToEnum("DatasetInput")[0]
    505523def DoubleInputEnum(): return StringToEnum("DoubleInput")[0]
     524def DoubleArrayInputEnum(): return StringToEnum("DoubleArrayInput")[0]
    506525def DataSetParamEnum(): return StringToEnum("DataSetParam")[0]
    507526def DoubleMatArrayParamEnum(): return StringToEnum("DoubleMatArrayParam")[0]
    … …  
    672691def GiaWEnum(): return StringToEnum("GiaW")[0]
    673692def P0Enum(): return StringToEnum("P0")[0]
     693def P0ArrayEnum(): return StringToEnum("P0Array")[0]
    674694def P1Enum(): return StringToEnum("P1")[0]
    675695def P1DGEnum(): return StringToEnum("P1DG")[0]

  • issm/trunk-jpl/test/NightlyRun/gemb.m

    r19527 r19554  
    2323md.smb.pAir=[repmat(inputs.pAir0',md.mesh.numberofelements,1);dateN'];
    2424
    25 md.smb.Vz=inputs.LP.Vz;
    26 md.smb.Tz=inputs.LP.Tz;
    27 md.smb.Tmean=inputs.LP.Tmean;
    28 md.smb.C=inputs.LP.C;
     25md.smb.Vz=repmat(inputs.LP.Vz,md.mesh.numberofelements,1);
     26md.smb.Tz=repmat(inputs.LP.Tz,md.mesh.numberofelements,1);
     27md.smb.Tmean=repmat(inputs.LP.Tmean,md.mesh.numberofelements,1);
     28md.smb.C=repmat(inputs.LP.C,md.mesh.numberofelements,1);
    2929
    3030%settings
    31 md.smb.spinUp=2;
     31md.smb.ismelt=0;
     32md.smb.isthermal=0;
     33md.smb.isalbedo=0;
     34md.transient.isstressbalance=0;
     35md.transient.ismasstransport=0;
     36md.transient.isthermal=0;
     37md.timestepping.start_time=2000;
     38md.timestepping.final_time=2001;
     39md.timestepping.time_step=1/12; %monthly
     40
     41md.verbose=verbose('11111111');
    3242
    3343%Run transient
Note: See TracChangeset for help on using the changeset viewer.