         NiuE2normq = 0.5*dg::blas2::dot(Ni, w2d,uE2);
         
         
-        Edfnorm = uE2norm + dnenorm*dnenorm ;
+        Edfnorm = uE2norm + 0.5*dnenorm*dnenorm;
         Enorm = NiuE2norm + nlnnnorm ;
         Eqnorm = NiuE2normq + nlnnnormq;
-        unorm = sqrt(dnenorm*dnenorm+ phinorm);
+
+        unorm = sqrt(dnenorm*dnenorm + phinorm*phinorm);
+        usupnorm = sqrt(dnesupnorm*dnesupnorm + phisupnorm*phisupnorm);
         //normalize
         if (i==0) {
-            phisupnorm0=1.+phisupnorm; //is initialy 0
-            phinorm0=1.+phinorm;//is initialy 0
+            phisupnorm0=1.+phisupnorm; //is zero at t=0
+            phinorm0=1.+phinorm; //is zero at t=0
             dnesupnorm0=dnesupnorm;
             dnenorm0=dnenorm;
             Edfnorm0=Edfnorm;
             Enorm0=Enorm;
             Eqnorm0=Eqnorm;
             unorm0=unorm;
+            usupnorm0=usupnorm;
         }
         //write norm data
         std::cout << time << " " <<  (1.+phisupnorm)/phisupnorm0<< " " 

                                 <<  dnenorm/dnenorm0<< " "
                                 <<  Edfnorm/Edfnorm0<<" "
                                 <<  Enorm/Enorm0<<" "
                                 <<  Eqnorm/Eqnorm0<<" "
                                 <<  unorm/unorm0<<"\n";

        time += deltaT;
    }

    double deltaT = p.dt*p.itstp;     //define timestep

    

    double phisupnorm, phisupnorm0, phinorm, phinorm0,dnesupnorm, dnesupnorm0, dnenorm, dnenorm0, dNinorm, dNinorm0, Enorm, Enorm0, Eqnorm, Eqnorm0, Edfnorm, Edfnorm0, uE2norm, NiuE2norm, NiuE2norm0, nlnnnorm, nlnnnorm0,NiuE2normq, NiuE2normq0, nlnnnormq, nlnnnormq0, unorm, unorm0, usupnorm, usupnorm0=0.;

    for( unsigned i=imin; i<=imax; i++)//timestepping
    {
