Changeset 19583

Timestamp:

09/25/15 09:12:29 (9 years ago)

Author:

Eric.Larour

Message:

CHG: finalized SMBgemb method and validated against Alex Gardner's model.

File:

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

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

@@ -2164 +2164 @@
         int         m;
         IssmDouble  SmbMassBalance=0;
+        int         count=0;
         /*}}}*/
 
@@ -2246 +2247 @@
                 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
@@ -2265 +2265 @@
                 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: */
@@ -2277 +2276 @@
                 a_input->GetValues(&a,&m);
                 T_input->GetValues(&T,&m);
-                swf_input->GetValues(&swf,&m);
+                
+                //fixed lower temperatuer bounday condition - T is fixed
+                T_bottom=T[m-1];
 
         } /*}}}*/
@@ -2292 +2293 @@
 
         /*Start loop: */
+        count=1;
         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");
+                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" << setprecision(3) << " Step: " << count << "\n");
 
                 /*extract daily data:{{{*/
@@ -2315 +2317 @@
                                         
                 /*Distribution of absorbed short wave radation with depth:*/
-        if(isshortwave)shortwave(swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
+        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,m,this->Sid());
                 
+                /*Calculate net shortwave [W m-2]*/
+        netSW = cellsum(swf,m);
+
                 /*Thermal profile computation:*/
         if(isthermal)thermo(&EC, T, dz, d, swf, dlw, Ta, V, eAir, pAir, W[0], smb_dt, m, Vz, Tz,this->Sid());     
@@ -2338 +2343 @@
         ulw = 5.67E-8 * pow(T[0],4.0);
 
-                /*Calculate net shortwave and longwave [W m-2]*/
-        netSW = cellsum(swf,m);
+                /*Calculate net longwave [W m-2]*/
         netLW = dlw - ulw;
                 
                 /*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,this->Sid());
+                
+                /*Verbose some resuls in debug mode: {{{*/
+                if(VerboseSmb() && 0){ 
+                        _printf_("smb log: count[" << count << "] m[" << m << "] " 
+                                << setprecision(16)   << "T[" << cellsum(T,m)  << "] " 
+                                                          << "d[" << cellsum(d,m)  << "] "
+                                                          << "dz[" << cellsum(dz,m)  << "] "
+                                                          << "a[" << cellsum(a,m)  << "] "
+                                                          << "W[" << cellsum(W,m)  << "] "
+                                                          << "re[" << cellsum(re,m)  << "] "
+                                                          << "gdn[" << cellsum(gdn,m)  << "] "
+                                                          << "gsp[" << cellsum(gsp,m)  << "] "
+                                                          << "swf[" << netSW << "] "
+                                                                          << "\n");
+                } /*}}}*/
                 
                 /*Sum component mass changes [kg m-2]*/
@@ -2359 +2378 @@
 
                 /*Check mass conservation:*/
-        if(!VerboseSmb())if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
+        if (dMass != 0.0) _printf_("total system mass not conserved in MB function");
                 
                 /*Check bottom grid cell T is unchanged:*/
-        if(!VerboseSmb())if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
+        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.
 
+                count++;
 
         } //for (t=time;t<=time+dt;t=t+smb_dt)
@@ -2392 +2412 @@
         xDelete<IssmDouble>(a);
         xDelete<IssmDouble>(T);
+        xDelete<IssmDouble>(swf);
+        delete gauss;
         /*}}}*/
 }