Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 600)
+++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 622)
@@ -316,5 +316,6 @@
       real dEtot, dEtots, masse, AtmToSurf_TurbFlux
       real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
-      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdif(ngridmx,nlayermx)
+      real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
+      real dEdiffs(ngridmx),dEdiff(ngridmx)
       real madjdE(ngridmx), lscaledE(ngridmx)
 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
@@ -559,5 +560,6 @@
 
             do ig=1,ngrid
-               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
+!already done above (JL12)
+!               qsurf_hist(ig,igcm_h2o_vap) = qsurf(ig,igcm_h2o_vap)
                if(ice_update)then
                   ice_initial(ig)=qsurf(ig,igcm_h2o_ice)
@@ -877,5 +879,4 @@
            do l=1,nlayer
              do ig=1,ngrid 
-               zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
                zdh(ig,l)=pdt(ig,l)/zpopsk(ig,l)
              enddo
@@ -931,16 +932,20 @@
             dEtot=0.0
             dEtots=0.0
-	    dEzdif(:,:)=0.
+	    dEzdiff(:,:)=0.
             AtmToSurf_TurbFlux=0.
-	    
+            dEdiffs(:)=0.
+            dEdiff(:)=0.
+	    	    
             do ig = 1, ngrid
                do l = 1, nlayer
                   masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                  dEtot = dEtot + cpp*masse*zdtdif(ig,l)*area(ig)
-		  dEzdif (ig,l)= cpp*masse*(zdtdif(ig,l))
-               enddo
-	       dEzdif(ig,1)= dEzdif(ig,1)+ sensibFlux(ig)! subtract flux to the ground
-               dEtot = dEtot + sensibFlux(ig)*area(ig)! subtract flux to the ground
-               dEtots = dEtots + capcal(ig)*zdtsdif(ig)*area(ig)-zflubid(ig)*area(ig)-sensibFlux(ig)*area(ig)
+		  dEzdiff (ig,l)= cpp*masse*(zdtdif(ig,l))
+		  dEdiff(ig)=dEdiff(ig)+dEzdiff (ig,l)
+               enddo
+	       dEzdiff(ig,1)= dEzdiff(ig,1)+ sensibFlux(ig)! subtract flux to the ground
+               dEdiff(ig)=dEdiff(ig)+ sensibFlux(ig)! subtract flux to the ground
+	       dEtot = dEtot + dEdiff(ig)*area(ig)
+	       dEdiffs(ig)=capcal(ig)*zdtsdif(ig)-zflubid(ig)-sensibFlux(ig)
+               dEtots = dEtots + dEdiffs(ig)*area(ig)
 	       AtmToSurf_TurbFlux=AtmToSurf_TurbFlux+sensibFlux(ig)*area(ig)
             enddo
@@ -1342,19 +1347,4 @@
                  rainout(ig)             = zdqsrain(ig)+zdqssnow(ig) ! diagnostic   
               enddo
-
-              !-------------------------
-              ! test energy conservation
-              if(enertest)then
-                 dEtot=0.0
-                 do ig = 1, ngrid
-                    do l = 1, nlayer
-                       masse = (pplev(ig,l) - pplev(ig,l+1))/g
-                       dEtot = dEtot + cpp*masse*zdtrain(ig,l)*area(ig)
-                    enddo
-                 enddo
-                 dEtot=dEtot/totarea
-                 print*,'In rain atmospheric energy change       =',dEtot,' W m-2'
-              endif
-              !-------------------------
 
 
@@ -1572,15 +1562,5 @@
          ! content of ocean gridpoints back to zero, in order
          ! to avoid rounding errors in vdifc, rain
-         do ig = 1, ngrid
-            do iq = 1, nq
-               if(iq.eq.igcm_h2o_vap .and. rnat(ig).eq.0)then ! if liquid water and terrain = ocean 
-                  qsurf_hist(ig,iq) = qsurf(ig,iq)
-                  !qsurf(ig,iq)      = qcol(ig,iq) 
-                  ! the value of qsurf we choose here makes NO DIFFERENCE TO ANYTHING AT ALL
-               else
-                  qsurf_hist(ig,iq) = qsurf(ig,iq)
-               endif
-            enddo
-         enddo
+         qsurf_hist(:,:) = qsurf(:,:)
 
          if(ice_update)then
@@ -2016,5 +1996,11 @@
 	
 	if(enertest) then
-	  if (calldifv) call writediagfi(ngridmx,"dEzdif","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdif)
+	  if (calldifv) then
+	     call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
+	     call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
+	     call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
+	     call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
+	     call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
+	  endif
 	  if (corrk) then
 	     call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
@@ -2023,5 +2009,7 @@
           if(watercond) then
 	     call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE) 
-	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE) 
+	     call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
+	     call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
+	     call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
           endif 
         endif
@@ -2054,10 +2042,8 @@
                           'kg m^-2',2,qcol(1,iq) )
 
-          if(water)then
-             call writediagfi(ngridmx,"H2Omaxcol","max. poss. H2O column","kg m^-2",2,H2Omaxcol)
-          endif
-
           if(watercond.or.CLFvarying)then
-             !call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
+!             call writediagfi(ngrid,"rneb_man","H2O cloud fraction (conv)"," ",3,rneb_man)
+!             call writediagfi(ngrid,"rneb_lsc","H2O cloud fraction (large scale)"," ",3,rneb_lsc)
+!             call writediagfi(ngrid,"CLF","H2O cloud fraction"," ",3,cloudfrac)
              call writediagfi(ngrid,"CLFt","H2O column cloud fraction"," ",2,totcloudfrac)
           endif
