Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19566) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 19567) @@ -2137,5 +2137,4 @@ int swIdx=0; IssmDouble cldFrac,t0wet, t0dry, K; - IssmDouble comp1,comp2; IssmDouble ulw; IssmDouble netSW; @@ -2146,4 +2145,5 @@ IssmDouble sumMass,dMass; bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux; + IssmDouble init_scaling=1.0; /*}}}*/ /*Output variables:{{{ */ @@ -2163,8 +2163,10 @@ IssmDouble mAdd; int m; + IssmDouble SmbMassBalance=0; /*}}}*/ /*only compute SMB at the surface: */ if (!IsOnSurface()) return; + /*Retrieve material properties and parameters:{{{ */ @@ -2190,4 +2192,5 @@ parameters->FindParam(&isdensification,SmbIsdensificationEnum); parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum); + parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum); /*}}}*/ /*Retrieve inputs: {{{*/ @@ -2225,8 +2228,9 @@ Vz_input->GetInputValue(&Vz,gauss); /*}}}*/ + /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/ - if(!isinitialized_input){ - - /*Initialize grid:*/ + if(!isinitialized_input){ + + if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n"); GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY); //if(this->Sid()==1) for(int i=0;iSid()); /*Snow, firn and ice albedo:*/ - if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m); + if(isalbedo)albedo(a,aIdx,re,d,cldFrac,aIce, aSnow,T,W,P,EC,t0wet,t0dry,K,smb_dt,m,this->Sid()); + /*Distribution of absorbed short wave radation with depth:*/ - if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m); + if(isshortwave)shortwave(swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid()); /*Thermal profile computation:*/ - if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, Vz, Tz,m); + if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid()); /*Change in thickness of top cell due to evaporation/condensation assuming same density as top cell. @@ -2311,16 +2325,12 @@ /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/ - if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow); + if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, Ta, P, dzMin, aSnow,this->Sid()); /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/ - comp2 = cellsum(dz,m); - if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin); - comp2 = (comp2 - cellsum(dz,m)); + if(ismelt)melt(&M, &R, &mAdd, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin,this->Sid()); /*Allow non-melt densification and determine compaction [m]*/ - comp1 = cellsum(dz,m); - if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,m); - comp1 = (comp1 - cellsum(dz,m)); + if(isdensification)densification(d,dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid()); /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every @@ -2333,5 +2343,5 @@ /*Calculate turbulent heat fluxes [W m-2]*/ - if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz); + if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,this->Sid()); /*Sum component mass changes [kg m-2]*/ @@ -2347,11 +2357,15 @@ dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd; dMass = round(dMass * 100.0)/100.0; - + /*Check mass conservation:*/ - if (dMass != 0.0) _error_("total system mass not conserved in MB function"); + if(!VerboseSmb())if (dMass != 0.0) _printf_("total system mass not conserved in MB function"); /*Check bottom grid cell T is unchanged:*/ - if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom"); - } + if(!VerboseSmb())if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n"); + + SmbMassBalance += sumP + sumEC - sumR; //increment SMB for the entire time span of ice-flow dynamics. + + + } //for (t=time;t<=time+dt;t=t+smb_dt) /*Save generated inputs: */ @@ -2361,8 +2375,11 @@ this->AddInput(new DoubleArrayInput(SmbGdnEnum,gdn,m)); this->AddInput(new DoubleArrayInput(SmbGspEnum,gsp,m)); + this->AddInput(new DoubleArrayInput(SmbTEnum,T,m)); this->AddInput(new DoubleInput(SmbECEnum,EC)); this->AddInput(new DoubleArrayInput(SmbWEnum,W,m)); this->AddInput(new DoubleArrayInput(SmbAEnum,a,m)); - this->AddInput(new DoubleArrayInput(SmbTEnum,T,m)); + this->AddInput(new DoubleArrayInput(SmbSwfEnum,swf,m)); + this->AddInput(new DoubleInput(SmbMassBalanceEnum,time)); + /*Free allocations:{{{*/