Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 596)
+++ /trunk/LMDZ.GENERIC/README	(revision 597)
@@ -639,2 +639,3 @@
 - Cleaned continuum part in optc*
 - Added .def files for a typical 1d earth case in deftank (dry case for the moment)
+- Corrects a bug on potential temperature calculation in physic 
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 596)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 597)
@@ -1,3 +1,3 @@
-     subroutine physiq(ngrid,nlayer,nq,   &
+      subroutine physiq(ngrid,nlayer,nq,   &
                   firstcall,lastcall,      &
                   pday,ptime,ptimestep,    &
@@ -202,7 +202,7 @@
       real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
  
-      real sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
+      real,save :: sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
  
-      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,sensibFlux
+      save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
 
 
@@ -226,4 +226,5 @@
       real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
       real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
+      real zdtdif(ngridmx,nlayermx)			      ! (K/s)
       real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
       real zdhadj(ngridmx,nlayermx)                           ! (K/s)
@@ -267,5 +268,5 @@
       character*2 str2
       character*5 str5
-      real zdtdif(ngridmx,nlayermx), zdtadj(ngridmx,nlayermx)
+      real zdtadj(ngridmx,nlayermx)
       real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
       save ztprevious 
@@ -285,5 +286,5 @@
       real tau_col(ngridmx)
       save tau_col
-
+      
 !     included by RW to reduce insanity of code
       real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS
@@ -311,4 +312,5 @@
 !     included by RW to account for surface cooling by evaporation
       real dtsurfh2olat(ngridmx)
+
 
 !     to test energy conservation (RW+JL)
@@ -388,5 +390,4 @@
       real reffH2O(ngridmx,nlayermx)
       real reffcol(ngridmx,2)
-
 
 !     included by RW for sourceevol
@@ -550,5 +551,5 @@
 
             if(ngrid.eq.1)then
-               qzero1D               = 20.0
+               qzero1D               = 10000.0
                qsurf(1,igcm_h2o_vap) = qzero1D
                do l=1, nlayer
@@ -633,8 +634,8 @@
 
 
-         
-      do l=1,nlayer
+      do l=1,nlayer         
          do ig=1,ngrid 
             zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
+            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
          enddo
       enddo
@@ -938,6 +939,7 @@
                   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))*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)
@@ -947,7 +949,13 @@
             dEtots = dEtots/totarea
 	    AtmToSurf_TurbFlux=AtmToSurf_TurbFlux/totarea
-            print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
-            print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
-            print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+            if (UseTurbDiff) then
+	       print*,'In TurbDiff sensible flux (atm=>surf) =',AtmToSurf_TurbFlux,' W m-2'
+               print*,'In TurbDiff non-cons atm nrj change   =',dEtot,' W m-2'
+               print*,'In TurbDiff (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	    else
+	       print*,'In vdifc sensible flux (atm=>surf)    =',AtmToSurf_TurbFlux,' W m-2'
+               print*,'In vdifc non-cons atm nrj change      =',dEtot,' W m-2'
+               print*,'In vdifc (correc rad+latent heat) surf nrj change =',dEtots,' W m-2'
+	    end if
 ! JL12 note that the black body radiative flux emitted by the surface has been updated by the implicit scheme
 !    but not given back elsewhere
@@ -1609,5 +1617,5 @@
          enddo
          dEtots=dEtots/totarea
-         print*,'Surface energy change=',dEtots,' W m-2'
+         print*,'Surface energy change                 =',dEtots,' W m-2'
       endif
 !-------------------------
@@ -1992,4 +2000,5 @@
         call writediagfi(ngrid,"ps","Surface pressure","Pa",2,ps)
         call writediagfi(ngrid,"temp","temperature","K",3,zt)
+        call writediagfi(ngrid,"teta","potential temperature","K",3,zh)
         call writediagfi(ngrid,"u","Zonal wind","m.s-1",3,zu)
         call writediagfi(ngrid,"v","Meridional wind","m.s-1",3,zv)
