![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 21232
- Timestamp:
- 09/27/16 16:49:37 (8 years ago)
- Location:
- issm/trunk-jpl/src
- Files:
-
- 7 edited
-
c/analyses/SmbAnalysis.cpp (modified) (1 diff)
-
c/classes/Elements/Element.cpp (modified) (13 diffs)
-
c/modules/SurfaceMassBalancex/Gembx.cpp (modified) (27 diffs)
-
c/shared/Enum/EnumDefinitions.h (modified) (3 diffs)
-
c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
-
c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
-
m/classes/SMBgemb.m (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp
r21216 r21232 58 58 iomodel->FetchDataToInput(elements,"md.smb.Tz",SmbTzEnum); 59 59 iomodel->FetchDataToInput(elements,"md.smb.Vz",SmbVzEnum); 60 iomodel->FetchDataToInput(elements,"md.smb.isInitialized",SmbIsInitializedEnum); 61 iomodel->FetchDataToInput(elements,"md.smb.isrestart",SmbIsrestartEnum); 62 iomodel->FetchDataToInput(elements,"md.smb.Dzrst",SmbDzrstEnum); 63 iomodel->FetchDataToInput(elements,"md.smb.Drst",SmbDrstEnum); 64 iomodel->FetchDataToInput(elements,"md.smb.Rerst",SmbRerstEnum); 65 iomodel->FetchDataToInput(elements,"md.smb.Gdnrst",SmbGdnrstEnum); 66 iomodel->FetchDataToInput(elements,"md.smb.Gsprst",SmbGsprstEnum); 67 iomodel->FetchDataToInput(elements,"md.smb.ECrst",SmbECrstEnum); 68 iomodel->FetchDataToInput(elements,"md.smb.Wrst",SmbWrstEnum); 69 iomodel->FetchDataToInput(elements,"md.smb.Arst",SmbArstEnum); 70 iomodel->FetchDataToInput(elements,"md.smb.Trst",SmbTrstEnum); 71 iomodel->FetchDataToInput(elements,"md.smb.Sizerst",SmbSizerstEnum); 60 InputUpdateFromConstantx(elements,0.,SmbIsInitializedEnum); 61 iomodel->FetchDataToInput(elements,"md.smb.Dzini",SmbDziniEnum); 62 iomodel->FetchDataToInput(elements,"md.smb.Dini",SmbDiniEnum); 63 iomodel->FetchDataToInput(elements,"md.smb.Reini",SmbReiniEnum); 64 iomodel->FetchDataToInput(elements,"md.smb.Gdnini",SmbGdniniEnum); 65 iomodel->FetchDataToInput(elements,"md.smb.Gspini",SmbGspiniEnum); 66 iomodel->FetchDataToInput(elements,"md.smb.ECini",SmbECiniEnum); 67 iomodel->FetchDataToInput(elements,"md.smb.Wini",SmbWiniEnum); 68 iomodel->FetchDataToInput(elements,"md.smb.Aini",SmbAiniEnum); 69 iomodel->FetchDataToInput(elements,"md.smb.Tini",SmbTiniEnum); 70 iomodel->FetchDataToInput(elements,"md.smb.Sizeini",SmbSizeiniEnum); 72 71 break; 73 72 case SMBpddEnum: -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r21224 r21232 1475 1475 name==SmbRefreezeEnum || 1476 1476 name==SmbEvaporationEnum || 1477 name==SmbIsInitializedEnum || 1477 1478 name==BasalforcingsGroundediceMeltingRateEnum || 1478 1479 name==BasalforcingsFloatingiceMeltingRateEnum || … … 2246 2247 2247 2248 /*Intermediary variables: {{{*/ 2248 bool isinitialized; 2249 bool isrestart; 2250 IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; 2249 IssmDouble isinitialized; 2250 IssmDouble zTop,dzTop,zMax,zMin,zY,dzMin; 2251 2251 IssmDouble Tmean; 2252 2252 IssmDouble C; … … 2266 2266 IssmDouble lhf, shf, dayEC; 2267 2267 IssmDouble initMass; 2268 2269 2270 2268 IssmDouble sumR, sumM, sumEC, sumP, sumW,sumMassAdd; 2269 IssmDouble sumdz_add; 2270 IssmDouble sumMass,dMass; 2271 2271 bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; 2272 2272 IssmDouble init_scaling; 2273 2273 2274 /*}}}*/ 2274 2275 /*Output variables:{{{ */ … … 2287 2288 IssmDouble R; 2288 2289 IssmDouble mAdd; 2289 2290 IssmDouble dz_add; 2290 2291 2291 IssmDouble* dzrst=NULL;2292 IssmDouble* drst= NULL;2293 IssmDouble* rerst= NULL;2294 IssmDouble* gdnrst= NULL;2295 IssmDouble* gsprst= NULL;2296 IssmDouble* Wrst= NULL;2297 IssmDouble* arst= NULL;2298 IssmDouble* Trst= NULL;2299 2292 IssmDouble* dzini=NULL; 2293 IssmDouble* dini = NULL; 2294 IssmDouble* reini = NULL; 2295 IssmDouble* gdnini = NULL; 2296 IssmDouble* gspini = NULL; 2297 IssmDouble* Wini = NULL; 2298 IssmDouble* aini = NULL; 2299 IssmDouble* Tini = NULL; 2300 2300 2301 int m; 2301 2302 int count=0; … … 2313 2314 parameters->FindParam(&time,TimeEnum); /*transient core time at which we run the smb core*/ 2314 2315 parameters->FindParam(&dt,TimesteppingTimeStepEnum); /*transient core time step*/ 2315 2316 parameters->FindParam(&yts,ConstantsYtsEnum); 2316 2317 parameters->FindParam(&smb_dt,SmbDtEnum); /*time period for the smb solution, usually smaller than the glaciological dt*/ 2317 2318 parameters->FindParam(&aIdx,SmbAIdxEnum); … … 2331 2332 parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); 2332 2333 parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); 2333 2334 2334 2335 /*}}}*/ 2335 2336 /*Retrieve inputs: {{{*/ … … 2351 2352 Input* eAir_input=this->GetInput(SmbEAirEnum); _assert_(eAir_input); 2352 2353 Input* pAir_input=this->GetInput(SmbPAirEnum); _assert_(pAir_input); 2353 Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input); 2354 Input* isrestart_input=this->GetInput(SmbIsrestartEnum); _assert_(isrestart_input); 2355 2354 Input* isinitialized_input=this->GetInput(SmbIsInitializedEnum); _assert_(isinitialized_input); 2356 2355 /*Retrieve input values:*/ 2357 2356 Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0); … … 2367 2366 Tz_input->GetInputValue(&Tz,gauss); 2368 2367 Vz_input->GetInputValue(&Vz,gauss); 2369 isinitialized_input->GetInputValue(&isinitialized); 2370 isrestart_input->GetInputValue(&isrestart); 2368 isinitialized_input->GetInputValue(&isinitialized); 2371 2369 /*}}}*/ 2372 2370 2373 2374 if(isinitialized){2371 /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ 2372 if(isinitialized==0.0){ 2375 2373 if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n"); 2376 if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n");2377 GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);2378 2374 //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" << 2379 2375 //dz[i] << "\n"); 2380 /*initialize profile variables:*/2381 d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor2382 re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5; //set grain size to old snow [mm]2383 gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0; //set grain dentricity to old snow2384 gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=0; //set grain sphericity to old snow2385 EC = 0; //surface evaporation (-) condensation (+) [kg m-2]2386 W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=0; //set water content to zero [kg m-2]2387 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow; //set albedo equal to fresh snow [fraction]2388 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]2389 2376 2377 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input); 2378 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input); 2379 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input); 2380 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input); 2381 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input); 2382 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input); 2383 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input); 2384 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input); 2385 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input); 2386 2387 dz_input->GetValues(&dzini,&m); 2388 d_input->GetValues(&dini,&m); 2389 re_input->GetValues(&reini,&m); 2390 gdn_input->GetValues(&gdnini,&m); 2391 gsp_input->GetValues(&gspini,&m); 2392 EC_input->GetInputValue(&EC); 2393 W_input->GetValues(&Wini,&m); 2394 a_input->GetValues(&aini,&m); 2395 T_input->GetValues(&Tini,&m); 2396 2397 /*Retrive the correct value of m (without the zeroes at the end)*/ 2398 Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input); 2399 Size_input->GetInputValue(&m); 2400 2401 if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too 2402 // if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w DEFAULT values\n"); 2403 2404 /*initialize profile variables:*/ 2405 GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); 2406 2407 d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3] 2408 re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0]; //set grain size to old snow [mm] 2409 gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0]; //set grain dentricity to old snow 2410 gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0]; //set grain sphericity to old snow 2411 W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0]; //set water content to zero [kg m-2] 2412 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0]; //set albedo equal to fresh snow [fraction] 2413 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] 2414 /* /!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties) because don't know Tmean yet when set default values. 2415 Default value of 0C given in SMBgemb.m is overwritten here with value of Tmean*/ 2416 2417 //fixed lower temperature bounday condition - T is fixed 2418 T_bottom=T[m-1]; 2419 } 2420 else{ //Retrieve snow properties from previous run. Need to provide values for all layers 2421 // if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n"); 2422 2423 dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i]; 2424 d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i]; 2425 re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i]; 2426 gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i]; 2427 gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i]; 2428 W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i]; 2429 a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i]; 2430 T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i]; 2431 2432 //fixed lower temperature bounday condition - T is fixed 2433 T_bottom=T[m-1]; 2434 } 2435 2436 /*Flag the initialization:*/ 2437 this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0)); 2438 } 2439 else{ 2440 /*Recover inputs: */ 2441 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input); 2442 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input); 2443 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input); 2444 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input); 2445 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input); 2446 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input); 2447 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input); 2448 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input); 2449 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input); 2450 2451 /*Recover arrays: */ 2452 dz_input->GetValues(&dz,&m); 2453 d_input->GetValues(&d,&m); 2454 re_input->GetValues(&re,&m); 2455 gdn_input->GetValues(&gdn,&m); 2456 gsp_input->GetValues(&gsp,&m); 2457 EC_input->GetInputValue(&EC); 2458 W_input->GetValues(&W,&m); 2459 a_input->GetValues(&a,&m); 2460 T_input->GetValues(&T,&m); 2461 2390 2462 //fixed lower temperature bounday condition - T is fixed 2391 2463 T_bottom=T[m-1]; 2392 2393 /*Flag the initialization:*/2394 this->AddInput(new BoolInput(SmbIsInitializedEnum,false));2395 }2396 else if(isrestart){ //Retrieve the snow properties from previous run2397 if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");2398 2399 /*Recover inputs: */2400 DoubleArrayInput* dzrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzrstEnum)); _assert_(dzrst_input);2401 DoubleArrayInput* drst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDrstEnum));_assert_(drst_input);2402 DoubleArrayInput* rerst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbRerstEnum));_assert_(rerst_input);2403 DoubleArrayInput* gdnrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnrstEnum));_assert_(gdnrst_input);2404 DoubleArrayInput* gsprst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGsprstEnum));_assert_(gsprst_input);2405 DoubleInput* ECrst_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECrstEnum));_assert_(ECrst_input);2406 DoubleArrayInput* Wrst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWrstEnum));_assert_(Wrst_input);2407 DoubleArrayInput* arst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbArstEnum));_assert_(arst_input);2408 DoubleArrayInput* Trst_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTrstEnum));_assert_(Trst_input);2409 2410 /*Recover arrays: */2411 dzrst_input->GetValues(&dzrst,&m);2412 drst_input->GetValues(&drst,&m);2413 rerst_input->GetValues(&rerst,&m);2414 gdnrst_input->GetValues(&gdnrst,&m);2415 gsprst_input->GetValues(&gsprst,&m);2416 ECrst_input->GetInputValue(&EC);2417 Wrst_input->GetValues(&Wrst,&m);2418 arst_input->GetValues(&arst,&m);2419 Trst_input->GetValues(&Trst,&m);2420 2421 /*Retrive the correct value of m (without the zeroes at the end)*/2422 Input* Sizerst_input=this->GetInput(SmbSizerstEnum); _assert_(Sizerst_input);2423 Sizerst_input->GetInputValue(&m);2424 2425 dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzrst[i];2426 d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=drst[i];2427 re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=rerst[i];2428 gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnrst[i];2429 gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gsprst[i];2430 W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wrst[i];2431 a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=arst[i];2432 T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Trst[i];2433 2434 //fixed lower temperatuer bounday condition - T is fixed2435 T_bottom=T[m-1];2436 2437 /*Flag the initialization:*/2438 this->AddInput(new BoolInput(SmbIsrestartEnum,false));2439 }2440 else{2441 /*Recover inputs: */2442 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);2443 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);2444 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);2445 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);2446 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);2447 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);2448 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);2449 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);2450 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);2451 2452 /*Recover arrays: */2453 dz_input->GetValues(&dz,&m);2454 d_input->GetValues(&d,&m);2455 re_input->GetValues(&re,&m);2456 gdn_input->GetValues(&gdn,&m);2457 gsp_input->GetValues(&gsp,&m);2458 EC_input->GetInputValue(&EC);2459 W_input->GetValues(&W,&m);2460 a_input->GetValues(&a,&m);2461 T_input->GetValues(&T,&m);2462 2463 //fixed lower temperatuer bounday condition - T is fixed2464 T_bottom=T[m-1];2465 2464 2466 2465 } /*}}}*/ … … 2469 2468 initMass=0; for(int i=0;i<m;i++) initMass += dz[i]*d[i] + W[i]; 2470 2469 2471 // initialize cumulative variables2472 2473 2470 // initialize cumulative variables 2471 sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0; 2472 sumdz_add=0; 2474 2473 2475 2474 //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. … … 2502 2501 2503 2502 /*Distribution of absorbed short wave radation with depth:*/ 2504 2503 if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid()); 2505 2504 2506 2505 /*Calculate net shortwave [W m-2]*/ 2507 2506 netSW = cellsum(swf,m); 2508 2507 2509 2508 /*Thermal profile computation:*/ 2510 if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid()); 2509 if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid()); 2511 2510 2512 2511 /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. 2513 2512 * need to fix this in case all or more of cell evaporates */ 2514 2513 dz[0] = dz[0] + EC / d[0]; 2515 2514 2516 2515 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ 2517 2518 2519 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K 2516 if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,this->Sid()); 2517 2518 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K 2520 2519 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ 2521 2520 if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,this->Sid()); 2522 2521 2523 2522 /*Allow non-melt densification and determine compaction [m]*/ 2524 2523 if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); 2525 2524 2526 2525 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every 2527 2526 * sub-time step in thermo equations*/ 2528 2527 ulw = 5.67E-8 * pow(T[0],4.0); 2529 2528 2530 2529 /*Calculate net longwave [W m-2]*/ 2531 2530 netLW = dlw - ulw; 2532 2531 2533 2532 /*Calculate turbulent heat fluxes [W m-2]*/ 2534 2533 if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid()); 2535 2534 2536 2535 /*Verbose some resuls in debug mode: {{{*/ … … 2550 2549 2551 2550 /*Sum component mass changes [kg m-2]*/ 2552 2553 2554 2555 2556 2557 2558 sumdz_add=dz_add+sumdz_add; //*CL* 2551 sumMassAdd = mAdd + sumMassAdd; 2552 sumM = M + sumM; 2553 sumR = R + sumR; 2554 sumW = cellsum(W,m); 2555 sumP = P + sumP; 2556 sumEC = sumEC + EC; // evap (-)/cond(+) 2557 sumdz_add=dz_add+sumdz_add; 2559 2558 2560 2559 /*Calculate total system mass:*/ 2561 2562 2563 2564 2560 sumMass=0; for(int i=0;i<m;i++) sumMass += dz[i]*d[i]; 2561 2562 #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable. 2563 dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd; 2565 2564 dMass = round(dMass * 100.0)/100.0; 2566 2565 2567 2566 /*Check mass conservation:*/ 2568 if (dMass != 0.0) _printf_("total system mass not conserved in MB function");2567 if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n"); 2569 2568 #endif 2570 2569 2571 2570 /*Check bottom grid cell T is unchanged:*/ 2572 2571 if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); 2573 2572 2574 2573 /*Free ressources: */ … … 2579 2578 } //for (t=time;t<time+dt;t=t+smb_dt) 2580 2579 2581 2580 /*Save generated inputs: */ 2582 2581 this->AddInput(new DoubleArrayInput(SmbDzEnum,dz,m)); 2583 2582 this->AddInput(new DoubleArrayInput(SmbDEnum,d,m)); … … 2589 2588 this->AddInput(new DoubleArrayInput(SmbWEnum,W,m)); 2590 2589 this->AddInput(new DoubleArrayInput(SmbAEnum,a,m)); 2591 2592 2593 2594 2595 2596 2597 2590 this->AddInput(new DoubleInput(SmbMassBalanceEnum,(sumP + sumEC -sumR)/yts)); 2591 this->AddInput(new DoubleInput(SmbRunoffEnum,sumR/yts)); 2592 this->AddInput(new DoubleInput(SmbPrecipitationEnum,sumP/yts)); 2593 this->AddInput(new DoubleInput(SmbDz_addEnum,sumdz_add/yts)); 2594 this->AddInput(new DoubleInput(SmbM_addEnum,sumMassAdd/yts)); 2595 2596 /*Free allocations:{{{*/ 2598 2597 xDelete<IssmDouble>(dz); 2599 2598 xDelete<IssmDouble>(d); -
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp
r21228 r21232 12 12 13 13 for(int i=0;i<femmodel->elements->Size();i++){ 14 15 14 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i)); 15 element->SmbGemb(); 16 16 } 17 17 … … 67 67 gpBottom++; 68 68 } 69 69 //initialize bottom vectors 70 70 dzB = xNewZeroInit<IssmDouble>(gpBottom); 71 71 gp0 = dzTop; … … 98 98 99 99 // calculates grain growth according to Fig. 9 of Marbouty, 1980 100 // ------NO GRAIN GROWTH FOR d > 400 kg m-3 because (H is set to zero)------100 // ------NO GRAIN GROWTH FOR d > 400 kg m-3 because H is set to zero------ 101 101 // ---------------this is a major limitation of the model------------------- 102 102 … … 354 354 355 355 if(aIdx==1){ 356 356 //function of effective grain radius 357 357 358 358 //convert effective radius to specific surface area [cm2 g-1] 359 359 IssmDouble S = 3.0 / (.091 * re[0]); 360 360 361 361 //determine broadband albedo … … 364 364 else if(aIdx==2){ 365 365 366 366 // Spectral fractions (Lefebre et al., 2003) 367 367 // [0.3-0.8um 0.8-1.5um 1.5-2.8um] 368 368 369 369 IssmDouble sF[3] = {0.606, 0.301, 0.093}; 370 370 371 371 // convert effective radius to grain size in meters … … 386 386 else if(aIdx==3){ 387 387 388 388 // a as a function of density 389 389 390 390 // calculate albedo … … 393 393 else if(aIdx==4){ 394 394 395 395 // exponential time decay & wetness 396 396 397 397 // change in albedo with time: … … 402 402 403 403 // initialize variables 404 405 406 407 404 IssmDouble* t0=xNew<IssmDouble>(m); 405 IssmDouble* T=xNew<IssmDouble>(m); 406 IssmDouble* t0warm=xNew<IssmDouble>(m); 407 IssmDouble* d_a=xNew<IssmDouble>(m); 408 408 409 409 // specify constants … … 418 418 419 419 // determine timescale for albedo decay 420 421 422 420 for(int i=0;i<m;i++)if(W[i]>0)t0[i]=t0wet; // wet snow timescale 421 for(int i=0;i<m;i++)T[i]=TK[i] - 273.15; // change T from K to degC 422 for(int i=0;i<m;i++) t0warm[i]= fabs(T[i]) * K + t0dry; //// 'warm' snow timescale 423 423 for(int i=0;i<m;i++)if(W[i]==0.0 && T[i]>=-10)t0[i]= t0warm[i]; 424 424 for(int i=0;i<m;i++)if(T[i]<-10) t0[i] = 10 * K + t0dry; // 'cold' snow timescale … … 440 440 // a_surf = a_wet - (a_wet - a_surf) * exp(-W_surf/W0); 441 441 442 443 444 445 446 442 /*Free ressources:*/ 443 xDelete<IssmDouble>(t0); 444 xDelete<IssmDouble>(T); 445 xDelete<IssmDouble>(t0warm); 446 xDelete<IssmDouble>(d_a); 447 447 448 448 } … … 1219 1219 IssmDouble Wi; 1220 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1221 IssmDouble Ztot; 1222 IssmDouble T_bot; 1223 IssmDouble m_bot; 1224 IssmDouble dz_bot; 1225 IssmDouble d_bot; 1226 IssmDouble W_bot; 1227 IssmDouble a_bot; 1228 IssmDouble re_bot; 1229 IssmDouble gdn_bot; 1230 IssmDouble gsp_bot; 1231 bool top=false; 1232 1232 1233 1234 1235 1236 1237 1238 1233 IssmDouble* Zcum=NULL; 1234 IssmDouble* dzMin2=NULL; 1235 IssmDouble zY2=1.025; 1236 bool lastCellFlag; 1237 int X1=0; 1238 int X2=0; 1239 1239 1240 1240 int D_size; … … 1242 1242 /*outputs:*/ 1243 1243 IssmDouble mAdd; 1244 1244 IssmDouble dz_add; 1245 1245 IssmDouble Rsum; 1246 1246 IssmDouble* T=*pT; … … 1284 1284 mAdd = 0; // mass added/removed to/from base of model [kg] 1285 1285 addE = 0; // energy added/removed to/from base of model [J] 1286 1286 dz_add=0; // thickness of the layer added/removed to/from base of model [m] 1287 1287 1288 1288 // calculate temperature excess above 0 deg C … … 1322 1322 } 1323 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i];//*CL* 1324 //// MELT, PERCOLATION AND REFREEZE 1325 1326 // run melt algorithm if there is melt water or excess pore water 1327 if ((cellsum(exsT,n) > 0) || (cellsum(exsW,n) > 0)){ 1328 // _printf_(""MELT OCCURS"); 1329 // check to see if thermal energy exceeds energy to melt entire cell 1330 // if so redistribute temperature to lower cells (temperature surplus) 1331 // (maximum T of snow before entire grid cell melts is a constant 1332 // LF/CI = 159.1342) 1333 surpT=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpT[i] = fmax(0, exsT [i]- 159.1342); 1334 1335 if (cellsum(surpT,n) > 0 ){ 1336 // _printf_("T Surplus"); 1337 // calculate surplus energy 1338 surpE=xNew<IssmDouble>(n); for(int i=0;i<n;i++)surpE[i] = surpT[i] * CI * m[i]; 1339 1339 1340 1341 1342 1343 T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];//*CL* 1340 int i = 0; 1341 while (cellsum(surpE,n) > 0){ 1342 // use surplus energy to increase the temperature of lower cell 1343 T[i+1] = surpE[i]/m[i+1]/CI + T[i+1]; 1344 1344 1345 1346 1345 exsT[i+1] = fmax(0, T[i+1] - CtoK) + exsT[i+1]; 1346 T[i+1] = fmin(CtoK, T[i+1]); 1347 1347 1348 1349 surpE[i+1] = surpT[i+1] * CI * m[i+1];//*CL* 1348 surpT[i+1] = fmax(0, exsT[i+1] - 159.1342); 1349 surpE[i+1] = surpT[i+1] * CI * m[i+1]; 1350 1350 1351 1352 1353 surpE[i] = 0; 1354 1355 1356 1351 // adjust current cell properties (again 159.1342 is the max T) 1352 exsT[i] = 159.1342; 1353 surpE[i] = 0; 1354 i = i + 1; 1355 } 1356 } 1357 1357 1358 1358 // convert temperature excess to melt [kg] … … 1364 1364 1365 1365 // initialize refreeze, runoff, flxDn and dW vectors [kg] 1366 IssmDouble* F = xNewZeroInit<IssmDouble>(n); 1367 IssmDouble* R=xNewZeroInit<IssmDouble>(n); 1366 IssmDouble* F = xNewZeroInit<IssmDouble>(n); 1367 IssmDouble* R=xNewZeroInit<IssmDouble>(n); 1368 1368 1369 1369 for(int i=0;i<n;i++)dW[i] = 0; … … 1390 1390 1391 1391 // if reaches impermeable ice layer all liquid water runs off (R) 1392 else if (d[i] >= dIce){ 1392 else if (d[i] >= dIce){ // dPHC = pore hole close off [kg m-3] 1393 1393 // _printf_("ICE LAYER"); 1394 1394 // no water freezes in this cell … … 1401 1401 R[i] = fmax(0, inM - dW[i]); // runoff 1402 1402 } 1403 // check if no energy to refreeze meltwater 1403 // check if no energy to refreeze meltwater 1404 1404 else if (maxF[i] == 0){ 1405 1405 // _printf_("REFREEZE == 0"); … … 1432 1432 IssmDouble F2 = 0; 1433 1433 1434 if (dW[i] < 0){ 1434 if (dW[i] < 0){ // excess pore water 1435 1435 dMax = (dIce - d[i])*dz_0; // maximum refreeze 1436 1436 IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze … … 1464 1464 1465 1465 // delete all cells with zero mass 1466 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; D=xNew<int>(D_size); 1466 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0)D_size++; D=xNew<int>(D_size); 1467 1467 D_size=0; for(int i=0;i<n;i++)if(m[i]!=0){ D[D_size] = i; D_size++;} 1468 1468 … … 1484 1484 for(int i=0;i<n;i++)dz[i] = m[i] / d[i]; 1485 1485 1486 //calculate Rsum: 1486 //calculate Rsum: 1487 1487 Rsum=cellsum(R,n); 1488 1488 … … 1492 1492 } 1493 1493 1494 1495 1496 1494 //Merging of cells as they are burried under snow. 1495 Zcum=xNew<IssmDouble>(n); 1496 dzMin2=xNew<IssmDouble>(n); 1497 1497 1498 1498 Zcum[0]=dz[0]; // Compute a cumulative depth vector 1499 1499 1500 1501 1502 1500 for (int i=1;i<n;i++){ 1501 Zcum[i]=Zcum[i-1]+dz[i]; 1502 } 1503 1503 1504 1505 1506 1507 1508 1509 1510 1511 1504 for (int i=0;i<n;i++){ 1505 if (Zcum[i]<=zTop){ 1506 dzMin2[i]=dzMin; 1507 } 1508 else{ 1509 dzMin2[i]=zY2*dzMin2[i-1]; 1510 } 1511 } 1512 1512 1513 1513 // check if depth is too small 1514 1514 X = 0; 1515 1515 for(int i=n-1;i>=0;i--){ 1516 1516 if(dz[i]<dzMin2[i]){ 1517 1517 X=i; 1518 1518 break; … … 1520 1520 } 1521 1521 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 if (dz [i] < dzMin2[i]){// merge top cells1533 //_printf_("dz > dzMin * 2')1534 1535 1536 1537 1538 1539 1540 1522 //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell) 1523 if(X==n-1){ 1524 lastCellFlag = true; 1525 X = X-1; 1526 } 1527 else{ 1528 lastCellFlag = false; 1529 } 1530 1531 for (int i = 0; i<=X;i++){ 1532 if (dz [i] < dzMin2[i]){ // merge top cells 1533 // _printf_("dz > dzMin * 2') 1534 // adjust variables as a linearly weighted function of mass 1535 IssmDouble m_new = m[i] + m[i+1]; 1536 T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new; 1537 a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new; 1538 re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new; 1539 gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new; 1540 gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new; 1541 1541 1542 1543 1544 1545 1546 1542 // merge with underlying grid cell and delete old cell 1543 dz [i+1] = dz[i] + dz[i+1]; // combine cell depths 1544 d[i+1] = m_new / dz[i+1]; // combine top densities 1545 W[i+1] = W[i+1] + W[i]; // combine liquid water 1546 m[i+1] = m_new; // combine top masses 1547 1547 1548 1549 1550 1551 1552 1553 1554 1555 //find closest cell to merge with1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1548 // set cell to 99999 for deletion 1549 m[i] = 99999; 1550 } 1551 } 1552 1553 //If last cell has to be merged 1554 if(lastCellFlag){ 1555 //find closest cell to merge with 1556 for(int i=n-2;i>=0;i--){ 1557 if(m[i]!=99999){ 1558 X2=i; 1559 X1=n-1; 1560 break; 1561 } 1562 } 1563 1564 // adjust variables as a linearly weighted function of mass 1565 IssmDouble m_new = m[X2] + m[X1]; 1566 T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new; 1567 a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new; 1568 re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new; 1569 gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new; 1570 gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new; 1571 1572 // merge with underlying grid cell and delete old cell 1573 dz [X1] = dz[X2] + dz[X1]; // combine cell depths 1574 d[X1] = m_new / dz[X1]; // combine top densities 1575 W[X1] = W[X1] + W[X2]; // combine liquid water 1576 m[X1] = m_new; // combine top masses 1577 1578 // set cell to 99999 for deletion 1579 m[X2] = 99999; 1580 } 1581 1581 1582 1582 // delete combined cells … … 1597 1597 n=D_size; 1598 1598 xDelete<int>(D); 1599 1599 1600 1600 // check if any of the top 10 cell depths are too large 1601 1601 X=0; … … 1611 1611 if (dz [i] > dzMin *2){ 1612 1612 1613 // _printf_("dz > dzMin * 2");1614 // split in two1615 cellsplit(&dz, n, i,.5);1616 cellsplit(&W, n, i,.5);1617 cellsplit(&m, n, i,.5);1618 cellsplit(&T, n, i,1.0);1619 cellsplit(&d, n, i,1.0);1620 cellsplit(&a, n, i,1.0);1621 cellsplit(&EI, n, i,1.0);1622 cellsplit(&EW, n, i,1.0);1623 cellsplit(&re, n, i,1.0);1624 cellsplit(&gdn, n, i,1.0);1625 cellsplit(&gsp, n, i,1.0);1626 n++;1627 X=X+1;1613 // _printf_("dz > dzMin * 2"); 1614 // split in two 1615 cellsplit(&dz, n, i,.5); 1616 cellsplit(&W, n, i,.5); 1617 cellsplit(&m, n, i,.5); 1618 cellsplit(&T, n, i,1.0); 1619 cellsplit(&d, n, i,1.0); 1620 cellsplit(&a, n, i,1.0); 1621 cellsplit(&EI, n, i,1.0); 1622 cellsplit(&EW, n, i,1.0); 1623 cellsplit(&re, n, i,1.0); 1624 cellsplit(&gdn, n, i,1.0); 1625 cellsplit(&gsp, n, i,1.0); 1626 n++; 1627 X=X+1; 1628 1628 } 1629 1629 else i++; 1630 1630 } 1631 1631 1632 1633 1634 1635 1636 1632 //// CORRECT FOR TOTAL MODEL DEPTH 1633 // WORKS FINE BUT HAS BEEN DISABLED FOR CONVIENCE OF MODEL OUTPUT 1634 // INTERPRETATION 1635 // // calculate total model depth 1636 Ztot = cellsum(dz,n); 1637 1637 1638 1639 //printf("Total depth < zMin %f \n", Ztot);1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 //printf("Total depth > zMax %f \n", Ztot);1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686 1687 1688 1689 1690 1691 1692 1638 if (Ztot < zMin){ 1639 // printf("Total depth < zMin %f \n", Ztot); 1640 // mass and energy to be added 1641 mAdd = m[n-1]+W[n-1]; 1642 addE = T[n-1]*m[n-1]*CI; 1643 1644 // add a grid cell of the same size and temperature to the bottom 1645 dz_bot=dz[n-1]; 1646 T_bot=T[n-1]; 1647 W_bot=W[n-1]; 1648 m_bot=m[n-1]; 1649 d_bot=d[n-1]; 1650 a_bot=a[n-1]; 1651 re_bot=re[n-1]; 1652 gdn_bot=gdn[n-1]; 1653 gsp_bot=gsp[n-1]; 1654 1655 dz_add=dz_bot; 1656 1657 newcell(&dz,dz_bot,top,n); 1658 newcell(&T,T_bot,top,n); 1659 newcell(&W,W_bot,top,n); 1660 newcell(&m,m_bot,top,n); 1661 newcell(&d,d_bot,top,n); 1662 newcell(&a,a_bot,top,n); 1663 newcell(&re,re_bot,top,n); 1664 newcell(&gdn,gdn_bot,top,n); 1665 newcell(&gsp,gsp_bot,top,n); 1666 n=n+1; 1667 } 1668 else if (Ztot > zMax){ 1669 // printf("Total depth > zMax %f \n", Ztot); 1670 // mass and energy loss 1671 mAdd = -(m[n-1]+W[n-1]); 1672 addE = -(T[n-1]*m[n-1]*CI); 1673 dz_add=-(dz[n-1]); 1674 1675 // add a grid cell of the same size and temperature to the bottom 1676 D_size=n-1; 1677 D=xNew<int>(D_size); 1678 1679 for(int i=0;i<n-1;i++) D[i]=i; 1680 celldelete(&dz,n,D,D_size); 1681 celldelete(&T,n,D,D_size); 1682 celldelete(&W,n,D,D_size); 1683 celldelete(&m,n,D,D_size); 1684 celldelete(&d,n,D,D_size); 1685 celldelete(&a,n,D,D_size); 1686 celldelete(&re,n,D,D_size); 1687 celldelete(&gdn,n,D,D_size); 1688 celldelete(&gsp,n,D,D_size); 1689 n=D_size; 1690 xDelete<int>(D); 1691 1692 } 1693 1693 1694 1694 //// CHECK FOR MASS AND ENERGY CONSERVATION … … 1709 1709 dm = round(mSum0 - mSum1 + mAdd); 1710 1710 dE = round(sumE0 - sumE1 - sumER + addE); 1711 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n" 1711 if (dm !=0 || dE !=0) _error_("mass or energy are not conserved in melt equations\n" 1712 1712 << "dm: " << dm << " dE: " << dE << "\n"); 1713 1713 #endif … … 1726 1726 if(D)xDelete<int>(D); 1727 1727 if(M)xDelete<IssmDouble>(M); 1728 1729 1728 xDelete<IssmDouble>(Zcum); 1729 xDelete<IssmDouble>(dzMin2); 1730 1730 1731 1731 /*Assign output pointers:*/ … … 1733 1733 *pR=Rsum; 1734 1734 *pmAdd=mAdd; 1735 1735 *pdz_add=dz_add; 1736 1736 1737 1737 *pT=T; -
issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h
r21216 r21232 166 166 HydrologysommersEnum, 167 167 HydrologyHeadEnum, 168 168 HydrologyHeadOldEnum, 169 169 HydrologyGapHeightEnum, 170 170 HydrologyBumpSpacingEnum, … … 361 361 SmbRequestedOutputsEnum, 362 362 SmbIsInitializedEnum, 363 SmbIsrestartEnum, 364 SmbDzrstEnum, 365 SmbDrstEnum, 366 SmbRerstEnum, 367 SmbGdnrstEnum, 368 SmbGsprstEnum, 369 SmbECrstEnum, 370 SmbWrstEnum, 371 SmbArstEnum, 372 SmbTrstEnum, 373 SmbSizerstEnum, 363 SmbDziniEnum, 364 SmbDiniEnum, 365 SmbReiniEnum, 366 SmbGdniniEnum, 367 SmbGspiniEnum, 368 SmbECiniEnum, 369 SmbWiniEnum, 370 SmbAiniEnum, 371 SmbTiniEnum, 372 SmbSizeiniEnum, 374 373 /*SMBforcing*/ 375 374 SMBforcingEnum, … … 423 422 SmbIsdensificationEnum, 424 423 SmbIsturbulentfluxEnum, 425 426 424 SmbDz_addEnum, 425 SmbM_addEnum, 427 426 /*SMBpdd*/ 428 427 SMBpddEnum, -
issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp
r21220 r21232 365 365 case SmbRequestedOutputsEnum : return "SmbRequestedOutputs"; 366 366 case SmbIsInitializedEnum : return "SmbIsInitialized"; 367 case SmbIsrestartEnum : return "SmbIsrestart"; 368 case SmbDzrstEnum : return "SmbDzrst"; 369 case SmbDrstEnum : return "SmbDrst"; 370 case SmbRerstEnum : return "SmbRerst"; 371 case SmbGdnrstEnum : return "SmbGdnrst"; 372 case SmbGsprstEnum : return "SmbGsprst"; 373 case SmbECrstEnum : return "SmbECrst"; 374 case SmbWrstEnum : return "SmbWrst"; 375 case SmbArstEnum : return "SmbArst"; 376 case SmbTrstEnum : return "SmbTrst"; 377 case SmbSizerstEnum : return "SmbSizerst"; 367 case SmbDziniEnum : return "SmbDzini"; 368 case SmbDiniEnum : return "SmbDini"; 369 case SmbReiniEnum : return "SmbReini"; 370 case SmbGdniniEnum : return "SmbGdnini"; 371 case SmbGspiniEnum : return "SmbGspini"; 372 case SmbECiniEnum : return "SmbECini"; 373 case SmbWiniEnum : return "SmbWini"; 374 case SmbAiniEnum : return "SmbAini"; 375 case SmbTiniEnum : return "SmbTini"; 376 case SmbSizeiniEnum : return "SmbSizeini"; 378 377 case SMBforcingEnum : return "SMBforcing"; 379 378 case SmbMassBalanceEnum : return "SmbMassBalance"; -
issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp
r21220 r21232 371 371 else if (strcmp(name,"SmbRequestedOutputs")==0) return SmbRequestedOutputsEnum; 372 372 else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum; 373 else if (strcmp(name,"SmbIsrestart")==0) return SmbIsrestartEnum; 374 else if (strcmp(name,"SmbDzrst")==0) return SmbDzrstEnum; 375 else if (strcmp(name,"SmbDrst")==0) return SmbDrstEnum; 376 else if (strcmp(name,"SmbRerst")==0) return SmbRerstEnum; 377 else if (strcmp(name,"SmbGdnrst")==0) return SmbGdnrstEnum; 378 else if (strcmp(name,"SmbGsprst")==0) return SmbGsprstEnum; 379 else if (strcmp(name,"SmbECrst")==0) return SmbECrstEnum; 380 else if (strcmp(name,"SmbWrst")==0) return SmbWrstEnum; 381 else if (strcmp(name,"SmbArst")==0) return SmbArstEnum; 382 else if (strcmp(name,"SmbTrst")==0) return SmbTrstEnum; 383 else if (strcmp(name,"SmbSizerst")==0) return SmbSizerstEnum; 373 else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum; 374 else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum; 375 else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum; 376 else if (strcmp(name,"SmbGdnini")==0) return SmbGdniniEnum; 377 else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum; 378 else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum; 379 else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum; 380 else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum; 381 else if (strcmp(name,"SmbTini")==0) return SmbTiniEnum; 382 else if (strcmp(name,"SmbSizeini")==0) return SmbSizeiniEnum; 384 383 else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum; 384 else if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; 385 385 else stage=4; 386 386 } 387 387 if(stage==4){ 388 if (strcmp(name,"SmbMassBalance")==0) return SmbMassBalanceEnum; 389 else if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; 388 if (strcmp(name,"SMBgemb")==0) return SMBgembEnum; 390 389 else if (strcmp(name,"SmbInitDensityScaling")==0) return SmbInitDensityScalingEnum; 391 390 else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; … … 506 505 else if (strcmp(name,"VxAverage")==0) return VxAverageEnum; 507 506 else if (strcmp(name,"Vx")==0) return VxEnum; 507 else if (strcmp(name,"VxPicard")==0) return VxPicardEnum; 508 508 else stage=5; 509 509 } 510 510 if(stage==5){ 511 if (strcmp(name,"VxPicard")==0) return VxPicardEnum; 512 else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; 511 if (strcmp(name,"VyAverage")==0) return VyAverageEnum; 513 512 else if (strcmp(name,"Vy")==0) return VyEnum; 514 513 else if (strcmp(name,"VyPicard")==0) return VyPicardEnum; … … 629 628 else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; 630 629 else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; 630 else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; 631 631 else stage=6; 632 632 } 633 633 if(stage==6){ 634 if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; 635 else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; 634 if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; 636 635 else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; 637 636 else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; … … 752 751 else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; 753 752 else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; 753 else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; 754 754 else stage=7; 755 755 } 756 756 if(stage==7){ 757 if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum; 758 else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; 757 if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; 759 758 else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; 760 759 else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; … … 875 874 else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum; 876 875 else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; 876 else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; 877 877 else stage=8; 878 878 } 879 879 if(stage==8){ 880 if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; 881 else if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; 880 if (strcmp(name,"ExtrudeFromBaseAnalysis")==0) return ExtrudeFromBaseAnalysisEnum; 882 881 else if (strcmp(name,"ExtrudeFromTopAnalysis")==0) return ExtrudeFromTopAnalysisEnum; 883 882 else if (strcmp(name,"DepthAverageAnalysis")==0) return DepthAverageAnalysisEnum; -
issm/trunk-jpl/src/m/classes/SMBgemb.m
r21224 r21232 21 21 ismelt; 22 22 isdensification; 23 isturbulentflux; 24 isInitialized; 25 isrestart; 23 isturbulentflux; 26 24 27 25 %inputs: … … 39 37 Vz = NaN; %height above ground at which wind (V) eas sampled [m] 40 38 41 % Initialization of snow properties w restart 42 Dzrst= NaN; %cell depth (m)43 Drst= NaN; %snow density (kg m-3)44 Rerst= NaN; %effective grain size (mm)45 Gdnrst= NaN; %grain dendricity (0-1)46 Gsprst= NaN; %grain sphericity (0-1)47 ECrst= NaN; %evaporation/condensation (kg m-2)48 Wrst= NaN; %Water content (kg m-2)49 Arst= NaN; %albedo (0-1)50 Trst= NaN; %snow temperature (K)51 Sizerst= NaN; %Number of layers39 % Initialization of snow properties 40 Dzini = NaN; %cell depth (m) 41 Dini = NaN; %snow density (kg m-3) 42 Reini = NaN; %effective grain size (mm) 43 Gdnini = NaN; %grain dendricity (0-1) 44 Gspini = NaN; %grain sphericity (0-1) 45 ECini = NaN; %evaporation/condensation (kg m-2) 46 Wini = NaN; %Water content (kg m-2) 47 Aini = NaN; %albedo (0-1) 48 Tini = NaN; %snow temperature (K) 49 Sizeini = NaN; %Number of layers 52 50 53 51 %settings: … … 131 129 self.isdensification=1; 132 130 self.isturbulentflux=1; 133 self.isInitialized = true*ones(mesh.numberofelements,1);134 self.isrestart = false*ones(mesh.numberofelements,1);135 131 136 132 self.aIdx = 1; … … 142 138 self.InitDensityScaling = 1.0; 143 139 144 140 self.zMax=250*ones(mesh.numberofelements,1); 145 141 self.zMin=130*ones(mesh.numberofelements,1); 146 142 self.zY = 1.10*ones(mesh.numberofelements,1); … … 155 151 self.K = 7; 156 152 157 self.Dzrst=0.05*ones(mesh.numberofelements,1); 158 self.Drst=910.0*ones(mesh.numberofelements,1); 159 self.Rerst=2.5*ones(mesh.numberofelements,1); 160 self.Gdnrst=0.0*ones(mesh.numberofelements,1); 161 self.Gsprst=0.0*ones(mesh.numberofelements,1); 162 self.ECrst=0.0*ones(mesh.numberofelements,1); 163 self.Wrst=0.0*ones(mesh.numberofelements,1); 164 self.Arst=self.aSnow*ones(mesh.numberofelements,1); 165 self.Trst=273.15*ones(mesh.numberofelements,1); %*CL* Default value must be Tmean 166 self.Sizerst=200*ones(mesh.numberofelements,1); 153 self.Dzini=0.05*ones(mesh.numberofelements,2); 154 self.Dini=910.0*ones(mesh.numberofelements,2); 155 self.Reini=2.5*ones(mesh.numberofelements,2); 156 self.Gdnini=0.0*ones(mesh.numberofelements,2); 157 self.Gspini=0.0*ones(mesh.numberofelements,2); 158 self.ECini=0.0*ones(mesh.numberofelements,1); 159 self.Wini=0.0*ones(mesh.numberofelements,2); 160 self.Aini=self.aSnow*ones(mesh.numberofelements,2); 161 self.Tini=273.15*ones(mesh.numberofelements,2); 162 % /!\ Default value of Tini must be equal to Tmean but don't know Tmean yet (computed when atmospheric forcings are interpolated on mesh). 163 % If initialization without restart, this value will be overwritten when snow parameters are retrieved in Element.cpp 164 self.Sizeini=2*ones(mesh.numberofelements,1); 167 165 168 166 end % }}} … … 178 176 md = checkfield(md,'fieldname','smb.isdensification','values',[0 1]); 179 177 md = checkfield(md,'fieldname','smb.isturbulentflux','values',[0 1]); 180 md = checkfield(md,'fieldname','smb.isInitialized','values',[false true]);181 md = checkfield(md,'fieldname','smb.isrestart','values',[false true]);182 178 183 179 md = checkfield(md,'fieldname','smb.Ta','timeseries',1,'NaN',1,'Inf',1,'>',273-100,'<',273+100); %-100/100 celsius min/max value … … 222 218 end 223 219 md = checkfield(md,'fieldname','smb.requested_outputs','stringrow',1); 224 225 % Check that isinitialization and isrestart have different values226 if any(self.isInitialized==self.isrestart),227 error('isInitialized and isrestart have the same value');228 end229 220 230 221 end % }}} … … 266 257 '4: exponential time decay & wetness [Bougamont & Bamber, 2005]'}); 267 258 268 269 fielddisplay(self,'Dzrst','Initial cel depth when restart [m]');270 fielddisplay(self,'Drst','Initial snow density when restart [kg m-3]');271 fielddisplay(self,'Rerst','Initial grain size when restart [mm]');272 fielddisplay(self,'Gdnrst','Initial grain dendricity when restart [-]');273 fielddisplay(self,'Gsprst','Initial grain sphericity when restart [-]');274 fielddisplay(self,'ECrst','Initial evaporation/condensation when restart [kg m-2]');275 fielddisplay(self,'Wrst','Initial snow water content when restart [kg m-2]');276 fielddisplay(self,'Arst','Initial albedo when restart [-]');277 fielddisplay(self,'Trst','Initial snow temperature when restart [K]');278 fielddisplay(self,'Sizerst','Initial number of layers when restart [K]');259 %snow properties init 260 fielddisplay(self,'Dzini','Initial cel depth when restart [m]'); 261 fielddisplay(self,'Dini','Initial snow density when restart [kg m-3]'); 262 fielddisplay(self,'Reini','Initial grain size when restart [mm]'); 263 fielddisplay(self,'Gdnini','Initial grain dendricity when restart [-]'); 264 fielddisplay(self,'Gspini','Initial grain sphericity when restart [-]'); 265 fielddisplay(self,'ECini','Initial evaporation/condensation when restart [kg m-2]'); 266 fielddisplay(self,'Wini','Initial snow water content when restart [kg m-2]'); 267 fielddisplay(self,'Aini','Initial albedo when restart [-]'); 268 fielddisplay(self,'Tini','Initial snow temperature when restart [K]'); 269 fielddisplay(self,'Sizeini','Initial number of layers when restart [K]'); 279 270 280 271 … … 341 332 WriteData(fid,prefix,'object',self,'class','smb','fieldname','denIdx','format','Integer'); 342 333 WriteData(fid,prefix,'object',self,'class','smb','fieldname','InitDensityScaling','format','Double'); 343 344 WriteData(fid,prefix,'object',self,'class','smb','fieldname','isInitialized','format','BooleanMat','mattype',2);345 WriteData(fid,prefix,'object',self,'class','smb','fieldname','isrestart','format','BooleanMat','mattype',2);346 347 334 WriteData(fid,prefix,'object',self,'class','smb','fieldname','outputFreq','format','Double'); 348 335 WriteData(fid,prefix,'object',self,'class','smb','fieldname','aSnow','format','Double'); … … 353 340 WriteData(fid,prefix,'object',self,'class','smb','fieldname','K','format','Double'); 354 341 355 356 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzrst','format','DoubleMat','mattype',3);357 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Drst','format','DoubleMat','mattype',3);358 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Rerst','format','DoubleMat','mattype',3);359 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnrst','format','DoubleMat','mattype',3);360 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gsprst','format','DoubleMat','mattype',3);361 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECrst','format','DoubleMat','mattype',2);362 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wrst','format','DoubleMat','mattype',3);363 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Arst','format','DoubleMat','mattype',3);364 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Trst','format','DoubleMat','mattype',3);365 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizerst','format','IntMat','mattype',2);342 %snow properties init 343 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dzini','format','DoubleMat','mattype',3); 344 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Dini','format','DoubleMat','mattype',3); 345 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Reini','format','DoubleMat','mattype',3); 346 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gdnini','format','DoubleMat','mattype',3); 347 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Gspini','format','DoubleMat','mattype',3); 348 WriteData(fid,prefix,'object',self,'class','smb','fieldname','ECini','format','DoubleMat','mattype',2); 349 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Wini','format','DoubleMat','mattype',3); 350 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Aini','format','DoubleMat','mattype',3); 351 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tini','format','DoubleMat','mattype',3); 352 WriteData(fid,prefix,'object',self,'class','smb','fieldname','Sizeini','format','IntMat','mattype',2); 366 353 367 354 %figure out dt from forcings: … … 372 359 WriteData(fid,prefix,'data',dt,'name','md.smb.dt','format','Double','scale',yts); 373 360 374 375 361 % Check if smb_dt goes evenly into transient core time step 362 if (mod(md.timestepping.time_step,dt) >= 1e-10) 376 363 error('smb_dt/dt = %f. The number of SMB time steps in one transient core time step has to be an an integer',md.timestepping.time_step/dt); 377 364 end 378 365 379 366 %process requested outputs
Note:
See TracChangeset
for help on using the changeset viewer.
