Changeset 19567

Timestamp:

09/19/15 14:55:36 (10 years ago)

Author:

Eric.Larour

Message:

CHG: adding verbosity + fixing some serious bugs

File:

• issm/trunk-jpl/src/c/classes/Elements/Element.cpp (modified) (15 diffs)

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

@@ -2137 +2137 @@
         int        swIdx=0;
         IssmDouble cldFrac,t0wet, t0dry, K;
-        IssmDouble comp1,comp2;
         IssmDouble ulw;
         IssmDouble netSW;
@@ -2146 +2145 @@
     IssmDouble sumMass,dMass;
         bool isgraingrowth,isalbedo,isshortwave,isthermal,isaccumulation,ismelt,isdensification,isturbulentflux;
+        IssmDouble init_scaling=1.0;
         /*}}}*/
         /*Output variables:{{{ */
@@ -2163 +2163 @@
         IssmDouble  mAdd;
         int         m;
+        IssmDouble  SmbMassBalance=0;
         /*}}}*/
 
         /*only compute SMB at the surface: */
         if (!IsOnSurface()) return;
+
 
         /*Retrieve material properties and parameters:{{{ */
@@ -2190 +2192 @@
         parameters->FindParam(&isdensification,SmbIsdensificationEnum);
         parameters->FindParam(&isturbulentflux,SmbIsturbulentfluxEnum);
+        parameters->FindParam(&init_scaling,SmbInitDensityScalingEnum);
         /*}}}*/
         /*Retrieve inputs: {{{*/
@@ -2225 +2228 @@
         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;i<m;i++)_printf_("z[" << i << "]=" <<
@@ -2234 +2238 @@
 
                 /*initialize profile variables:*/
-                d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice; //ice density 
+                d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=rho_ice*init_scaling; //ice density scaled by a factor
                 re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=2.5;         //set grain size to old snow [mm] 
                 gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=0;         //set grain dentricity to old snow 
@@ -2242 +2246 @@
                 a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aSnow;         //set albedo equal to fresh snow [fraction] 
                 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]
+                swf = xNewZeroInit<IssmDouble>(m); 
 
                 //fixed lower temperatuer bounday condition - T is fixed
@@ -2251 +2256 @@
         else{ 
                 /*Recover inputs: */
-                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum));
-                DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));
-                DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));
-                DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));
-                DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));
-                DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));
-                DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));
-                DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));
-                DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));
+                DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
+                DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
+                DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
+                DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
+                DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
+                DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
+                DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
+                DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
+                DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
+                DoubleArrayInput* swf_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbSwfEnum));_assert_(swf_input);
                 
                 /*Recover arrays: */
@@ -2271 +2277 @@
                 a_input->GetValues(&a,&m);
                 T_input->GetValues(&T,&m);
+                swf_input->GetValues(&swf,&m);
 
         } /*}}}*/
@@ -2280 +2287 @@
     sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
 
+        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT. 
+        //go back to time - deltaT: 
+        time-=dt;
+
         /*Start loop: */
         for (t=time;t<=time+dt;t=t+smb_dt){
+
+                if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << "\n");
 
                 /*extract daily data:{{{*/
@@ -2295 +2308 @@
 
                 /*Snow grain metamorphism:*/
-                if(isgraingrowth)grainGrowth(re, gdn, gsp, T, dz, d, W, smb_dt, m, aIdx);
+                if(isgraingrowth)grainGrowth(re, gdn, gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
 
                 /*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 +2325 @@
                 
                 /*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 +2343 @@
                 
                 /*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 +2357 @@
         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 +2375 @@
         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:{{{*/