Index: /codes/icosagcm/trunk/src/diagnostics/diagflux.F90
===================================================================
--- /codes/icosagcm/trunk/src/diagnostics/diagflux.F90	(revision 603)
+++ /codes/icosagcm/trunk/src/diagnostics/diagflux.F90	(revision 604)
@@ -10,6 +10,6 @@
        f_massfluxt(:), f_qfluxt(:), & ! time-integrated mass flux and tracer flux
        f_qfluxt_lon(:), f_qfluxt_lat(:), & ! scalar flux reconstructed at cell centers
-       f_epot(:), f_ekin(:), f_enthalpy(:), & ! time-averaged potential E, kinetic E and enthalpy
-       f_epotfluxt(:), f_ekinfluxt(:), f_enthalpyfluxt(:) ! time averaged 'fluxes' of epot, ekin and enthalpy
+       f_ulont(:), f_thetat(:), f_epot(:), f_ekin(:), f_enthalpy(:), & ! time-averaged potential E, kinetic E and enthalpy
+       f_ulonfluxt(:), f_thetafluxt(:), f_epotfluxt(:), f_ekinfluxt(:), f_enthalpyfluxt(:) ! time averaged 'fluxes' of epot, ekin and enthalpy
   LOGICAL :: diagflux_on
   !$OMP THREADPRIVATE(diagflux_on)
@@ -26,4 +26,6 @@
     ll = MERGE(llm,1,diagflux_on)
     CALL allocate_field(f_masst,         field_t,type_real,ll,       name="masst")
+    CALL allocate_field(f_ulont,         field_t,type_real,ll,       name="ulont")
+    CALL allocate_field(f_thetat,        field_t,type_real,ll,       name="thetat")
     CALL allocate_field(f_epot,          field_t,type_real,ll,       name="epot")
     CALL allocate_field(f_ekin,          field_t,type_real,ll,       name="ekin")
@@ -31,4 +33,6 @@
     CALL allocate_field(f_qmasst,        field_t,type_real,ll,nqtot, name="qmasst")
     CALL allocate_field(f_massfluxt,     field_u,type_real,ll,       name="massfluxt")
+    CALL allocate_field(f_ulonfluxt,     field_u,type_real,ll,       name="ulonfluxt")
+    CALL allocate_field(f_thetafluxt,    field_u,type_real,ll,       name="thetafluxt")
     CALL allocate_field(f_epotfluxt,     field_u,type_real,ll,       name="epotfluxt")
     CALL allocate_field(f_ekinfluxt,     field_u,type_real,ll,       name="ekinfluxt")
@@ -50,4 +54,6 @@
        CALL swap_dimensions(ind)
        ZERO2(f_masst)
+       ZERO2(f_ulont)
+       ZERO2(f_thetat)
        ZERO2(f_epot)
        ZERO2(f_ekin)
@@ -55,4 +61,6 @@
        ZERO3(f_qmasst)
        ZERO2(f_massfluxt)
+       ZERO2(f_ulonfluxt)
+       ZERO2(f_thetafluxt)
        ZERO2(f_epotfluxt)
        ZERO2(f_ekinfluxt)
@@ -110,11 +118,11 @@
 !------------------------------------ Compute energy fluxes ---------------------------------------
 
-  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
+  SUBROUTINE diagflux_energy(frac, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_pk, f_hfluxt)
     REAL(rstd), INTENT(IN) :: frac
-    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_hfluxt(:)
+    TYPE(t_field),POINTER :: f_phis(:),f_rhodz(:),f_theta_rhodz(:),f_u(:), f_geopot(:), f_theta(:), f_pk(:), f_hfluxt(:)
     REAL(rstd), POINTER :: phis(:), rhodz(:,:), theta_rhodz(:,:,:), u(:,:), &
-         geopot(:,:), pk(:,:,:), hfluxt(:,:), &
-         epot(:,:), ekin(:,:), enthalpy(:,:), &
-         epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
+         geopot(:,:), theta(:,:,:), pk(:,:), hfluxt(:,:), &
+         ulont(:,:), thetat(:,:), epot(:,:), ekin(:,:), enthalpy(:,:), &
+         thetaflux(:,:), ulonflux(:,:), epotflux(:,:), ekinflux(:,:), enthalpyflux(:,:)
     INTEGER :: ind
     DO ind=1,ndomain
@@ -128,31 +136,38 @@
        u = f_u(ind)
        geopot = f_geopot(ind)
-       pk = f_theta(ind) ! buffer
+       theta = f_theta(ind) ! buffer
+       pk = f_pk(ind) ! buffer
+       ulont = f_ulont(ind)
+       thetat = f_thetat(ind)
        epot = f_epot(ind)
        ekin = f_ekin(ind)
        enthalpy = f_enthalpy(ind)
+       ulonflux = f_ulonfluxt(ind)
+       thetaflux = f_thetafluxt(ind)
        epotflux = f_epotfluxt(ind)
        ekinflux = f_ekinfluxt(ind)
        enthalpyflux = f_enthalpyfluxt(ind)
-       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epotflux, ekinflux, enthalpyflux)
+       CALL compute_diagflux_energy(frac,hfluxt, phis,rhodz,theta_rhodz,u, geopot,theta,pk, &
+            ulont, thetat, epot, ekin, enthalpy, &
+            ulonflux, thetaflux, epotflux, ekinflux, enthalpyflux)
     END DO
   END SUBROUTINE diagflux_energy
 
-  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,u, geopot,pk, epot,ekin,enthalpy, epot_flux, ekin_flux, enthalpy_flux)
+  SUBROUTINE compute_diagflux_energy(frac, massflux, phis,rhodz,theta_rhodz,ue, geopot,theta,pk, &
+       ulon, thetat, epot, ekin, enthalpy, &
+       ulon_flux, thetat_flux, epot_flux, ekin_flux, enthalpy_flux)
     USE disvert_mod, ONLY : ptop
     REAL(rstd), INTENT(IN) :: frac
-    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), u(3*iim*jjm,llm),&
+    REAL(rstd), INTENT(IN) :: massflux(3*iim*jjm,llm), ue(3*iim*jjm,llm),&
                               phis(iim*jjm), rhodz(iim*jjm,llm), theta_rhodz(iim*jjm,llm,nqtot)
-    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), pk(iim*jjm,llm) ! pk = buffer
-    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  epot, ekin, enthalpy
-    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  epot_flux, ekin_flux, enthalpy_flux    
-    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd 
+    REAL(rstd), INTENT(INOUT) :: geopot(iim*jjm,llm+1), theta(iim*jjm,llm), pk(iim*jjm,llm) ! theta,pk = buffers
+    REAL(rstd), INTENT(INOUT), DIMENSION(iim*jjm, llm)   ::  ulon, thetat, epot, ekin, enthalpy
+    REAL(rstd), INTENT(INOUT), DIMENSION(3*iim*jjm, llm) ::  ulon_flux, thetat_flux, epot_flux, ekin_flux, enthalpy_flux    
+    REAL(rstd) :: energy, p_ik, theta_ik, temp_ik, gv, Rd, cx,cy,cz, ux,uy,uz, ue_le,ulon_i
     INTEGER :: ij, l, ij_omp_begin_ext, ij_omp_end_ext
     Rd = kappa*cpp
     ! even if loops are of the _ext variant, we still need halo exchanges before reconstructing fluxes at cell centers
     ! => loop over interior region
-    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin_ext,ij_omp_end_ext)
-    ij_omp_begin_ext = ij_omp_begin_ext+ij_begin-1
-    ij_omp_end_ext = ij_omp_end_ext+ij_begin-1
+    CALL distrib_level(ij_begin_ext, ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
 #include "../kernels/energy_fluxes.k90"
   END SUBROUTINE compute_diagflux_energy
Index: /codes/icosagcm/trunk/src/diagnostics/observable.f90
===================================================================
--- /codes/icosagcm/trunk/src/diagnostics/observable.f90	(revision 603)
+++ /codes/icosagcm/trunk/src/diagnostics/observable.f90	(revision 604)
@@ -1,4 +1,6 @@
 MODULE observable_mod
   USE icosa
+  USE diagflux_mod
+  USE output_field_mod
   IMPLICIT NONE
   PRIVATE
@@ -38,5 +40,4 @@
     USE disvert_mod
     USE wind_mod
-    USE output_field_mod
     USE omp_para
     USE time_mod
@@ -47,5 +48,4 @@
     USE theta2theta_rhodz_mod
     USE omega_mod
-    USE diagflux_mod
     LOGICAL, INTENT(IN) :: init
     INTEGER :: l
@@ -166,21 +166,9 @@
           CALL output_field("massflux_lat",f_buf_ulat)
 
-          CALL transfert_request(f_epotfluxt,req_e1_vect) 
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_epotfluxt, f_buf_ulon, f_buf_ulat)
-          CALL output_field("epot_t", f_epot)
-          CALL output_field("epotflux_lon",f_buf_ulon)
-          CALL output_field("epotflux_lat",f_buf_ulat)
-
-          CALL transfert_request(f_ekinfluxt,req_e1_vect) 
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_ekinfluxt, f_buf_ulon, f_buf_ulat)
-          CALL output_field("ekin_t", f_ekin)
-          CALL output_field("ekinflux_lon",f_buf_ulon)
-          CALL output_field("ekinflux_lat",f_buf_ulat)
-
-          CALL transfert_request(f_enthalpyfluxt,req_e1_vect) 
-          CALL flux_centered_lonlat(1./(itau_out*dt) , f_enthalpyfluxt, f_buf_ulon, f_buf_ulat)
-          CALL output_field("enthalpy_t", f_enthalpy)
-          CALL output_field("enthalpyflux_lon",f_buf_ulon)
-          CALL output_field("enthalpyflux_lat",f_buf_ulat)
+          CALL output_energyflux(f_ulont, f_ulonfluxt, "ulon_t", "ulonflux_lon", "ulonflux_lat")
+          CALL output_energyflux(f_thetat, f_thetafluxt, "theta_t", "thetaflux_lon", "thetaflux_lat")
+          CALL output_energyflux(f_epot, f_epotfluxt, "epot_t", "epotflux_lon", "epotflux_lat")
+          CALL output_energyflux(f_ekin, f_ekinfluxt, "ekin_t", "ekinflux_lon", "ekinflux_lat")
+          CALL output_energyflux(f_enthalpy, f_enthalpyfluxt, "enthalpy_t", "enthalpyflux_lon", "enthalpyflux_lat")
 
           CALL qflux_centered_lonlat(1./(itau_out*dt) , f_qfluxt, f_qfluxt_lon, f_qfluxt_lat)
@@ -192,4 +180,14 @@
     END IF
   END SUBROUTINE write_output_fields_basic
+
+  SUBROUTINE output_energyflux(f_energy, f_flux, name_energy, name_fluxlon, name_fluxlat)
+    TYPE(t_field), POINTER :: f_energy(:), f_flux(:)
+    CHARACTER(*), INTENT(IN) :: name_energy, name_fluxlon, name_fluxlat
+    CALL transfert_request(f_flux,req_e1_vect)
+    CALL flux_centered_lonlat(1./(itau_out*dt) , f_flux, f_buf_ulon, f_buf_ulat)
+    CALL output_field(name_energy,  f_energy)
+    CALL output_field(name_fluxlon, f_buf_ulon)
+    CALL output_field(name_fluxlat, f_buf_ulat)
+  END SUBROUTINE output_energyflux
   
  !------------------- Conversion from prognostic to observable variables ------------------
Index: /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90
===================================================================
--- /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 603)
+++ /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_base.F90	(revision 604)
@@ -43,7 +43,5 @@
 !$OMP BARRIER
 
-    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
-    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
-    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
+    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
 
     Rd = kappa*cpp
@@ -190,7 +188,5 @@
     CALL trace_start("compute_caldyn_vert")
 
-    CALL distrib_level(ij_end-ij_begin+1,ij_omp_begin,ij_omp_end)
-    ij_omp_begin=ij_omp_begin+ij_begin-1
-    ij_omp_end=ij_omp_end+ij_begin-1
+    CALL distrib_level(ij_begin,ij_end, ij_omp_begin,ij_omp_end)
 
     !    REAL(rstd) :: wwuu(iim*3*jjm,llm+1) ! tmp var, don't know why but gain 30% on the whole code in opemp
Index: /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 603)
+++ /codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90	(revision 604)
@@ -122,7 +122,5 @@
     INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
 
-    CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext)
-    ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1
-    ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1
+    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
 
     IF(dysl) THEN
Index: /codes/icosagcm/trunk/src/kernels/energy_fluxes.k90
===================================================================
--- /codes/icosagcm/trunk/src/kernels/energy_fluxes.k90	(revision 603)
+++ /codes/icosagcm/trunk/src/kernels/energy_fluxes.k90	(revision 604)
@@ -11,4 +11,6 @@
       END DO
    END DO
+   ! NB : at this point pressure is stored in array pk
+   ! pk then serves as buffer to store temperature
    SELECT CASE(caldyn_thermo)
    CASE(thermo_theta)
@@ -17,4 +19,5 @@
             p_ik = pk(ij,l)
             theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
+            theta(ij,l) = theta_ik
             temp_ik = theta_ik*(p_ik/preff)**kappa
             gv = (g*Rd)*temp_ik/p_ik
@@ -29,4 +32,5 @@
             theta_ik = theta_rhodz(ij,l,1)/rhodz(ij,l)
             temp_ik = Treff*exp((theta_ik + Rd*log(p_ik/preff))/cpp)
+            theta(ij,l) = Treff*exp(theta_ik/cpp)
             gv = (g*Rd)*temp_ik/p_ik ! specific volume v = Rd*T/p
             pk(ij,l) = temp_ik
@@ -37,7 +41,7 @@
    !$OMP BARRIER
    ! Now accumulate energies and energy fluxes
-   ! enthalpy
    ! NB : at this point temperature is stored in array pk
    ! pk then serves as buffer to store energy
+   ! enthalpy
    DO l = ll_begin, ll_end
       !DIR$ SIMD
@@ -73,4 +77,21 @@
       END DO
    END DO
+   ! theta
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         energy = theta(ij,l)
+         thetat(ij,l) = thetat(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         thetat_flux(ij+u_right,l) = thetat_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         thetat_flux(ij+u_lup,l) = thetat_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         thetat_flux(ij+u_ldown,l) = thetat_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
    ! kinetic energy
    DO l = ll_begin, ll_end
@@ -78,10 +99,10 @@
       DO ij=ij_begin_ext, ij_end_ext
          energy=0.d0
-         energy = energy + le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2
-         energy = energy + le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2
-         energy = energy + le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2
-         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2
-         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2
-         energy = energy + le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2
+         energy = energy + le(ij+u_rup)*de(ij+u_rup)*ue(ij+u_rup,l)**2
+         energy = energy + le(ij+u_lup)*de(ij+u_lup)*ue(ij+u_lup,l)**2
+         energy = energy + le(ij+u_left)*de(ij+u_left)*ue(ij+u_left,l)**2
+         energy = energy + le(ij+u_ldown)*de(ij+u_ldown)*ue(ij+u_ldown,l)**2
+         energy = energy + le(ij+u_rdown)*de(ij+u_rdown)*ue(ij+u_rdown,l)**2
+         energy = energy + le(ij+u_right)*de(ij+u_right)*ue(ij+u_right,l)**2
          energy = energy * (.25/Ai(ij))
          ekin(ij,l) = ekin(ij,l) + frac*rhodz(ij,l)*energy
@@ -97,4 +118,50 @@
       END DO
    END DO
+   ! ulon
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         cx=centroid(ij,1)
+         cy=centroid(ij,2)
+         cz=centroid(ij,3)
+         ux=0. ; uy=0. ; uz=0.
+         ue_le = ne_rup*ue(ij+u_rup,l)*le(ij+u_rup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_lup*ue(ij+u_lup,l)*le(ij+u_lup)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_up,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_up,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_up,3))-cz)
+         ue_le = ne_left*ue(ij+u_left,l)*le(ij+u_left)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_lup,1)+xyz_v(ij+z_ldown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_lup,2)+xyz_v(ij+z_ldown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_lup,3)+xyz_v(ij+z_ldown,3))-cz)
+         ue_le = ne_ldown*ue(ij+u_ldown,l)*le(ij+u_ldown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_ldown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_ldown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_ldown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_rdown*ue(ij+u_rdown,l)*le(ij+u_rdown)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rdown,1)+xyz_v(ij+z_down,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rdown,2)+xyz_v(ij+z_down,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rdown,3)+xyz_v(ij+z_down,3))-cz)
+         ue_le = ne_right*ue(ij+u_right,l)*le(ij+u_right)
+         ux = ux + ue_le*(.5*(xyz_v(ij+z_rup,1)+xyz_v(ij+z_rdown,1))-cx)
+         uy = uy + ue_le*(.5*(xyz_v(ij+z_rup,2)+xyz_v(ij+z_rdown,2))-cy)
+         uz = uz + ue_le*(.5*(xyz_v(ij+z_rup,3)+xyz_v(ij+z_rdown,3))-cz)
+         ulon_i = ux*elon_i(ij,1) + uy*elon_i(ij,2) + uz*elon_i(ij,3)
+         energy = ulon_i*(1./Ai(ij))
+         ulon(ij,l) = ulon(ij,l) + frac*rhodz(ij,l)*energy
+         pk(ij,l) = energy
+      END DO
+   END DO
+   DO l = ll_begin, ll_end
+      !DIR$ SIMD
+      DO ij=ij_begin_ext, ij_end_ext
+         ulon_flux(ij+u_right,l) = ulon_flux(ij+u_right,l) + .5*massflux(ij+u_right,l)*(pk(ij,l)+pk(ij+t_right,l))
+         ulon_flux(ij+u_lup,l) = ulon_flux(ij+u_lup,l) + .5*massflux(ij+u_lup,l)*(pk(ij,l)+pk(ij+t_lup,l))
+         ulon_flux(ij+u_ldown,l) = ulon_flux(ij+u_ldown,l) + .5*massflux(ij+u_ldown,l)*(pk(ij,l)+pk(ij+t_ldown,l))
+      END DO
+   END DO
    !---------------------------- energy_fluxes ----------------------------------
    !--------------------------------------------------------------------------
Index: /codes/icosagcm/trunk/src/parallel/omp_para.F90
===================================================================
--- /codes/icosagcm/trunk/src/parallel/omp_para.F90	(revision 603)
+++ /codes/icosagcm/trunk/src/parallel/omp_para.F90	(revision 604)
@@ -204,18 +204,18 @@
   END SUBROUTINE init_omp_para
 
-  SUBROUTINE distrib_level(size,lbegin,lend)
-  IMPLICIT NONE
-    INTEGER,INTENT(IN)  :: size  
+  SUBROUTINE distrib_level(ibegin,iend, lbegin,lend)
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)  :: ibegin,iend  
     INTEGER,INTENT(OUT) :: lbegin  
     INTEGER,INTENT(OUT) :: lend  
-    INTEGER :: div,rest
-    
+    INTEGER :: size,div,rest
+    size=iend-ibegin+1
     div=size/omp_level_size
     rest=MOD(size,omp_level_size)
     IF (omp_level_rank<rest) THEN
-      lbegin=(div+1)*omp_level_rank+1
+      lbegin=(div+1)*omp_level_rank + ibegin
       lend=lbegin+div
     ELSE
-      lbegin=(div+1)*rest + (omp_level_rank-rest)*div+1
+      lbegin=(div+1)*rest + (omp_level_rank-rest)*div + ibegin
       lend=lbegin+div-1
     ENDIF
Index: /codes/icosagcm/trunk/src/parallel/transfert_mpi.f90
===================================================================
--- /codes/icosagcm/trunk/src/parallel/transfert_mpi.f90	(revision 603)
+++ /codes/icosagcm/trunk/src/parallel/transfert_mpi.f90	(revision 604)
@@ -1237,5 +1237,5 @@
 
           dim3=size(field(ind)%rval3d,2)
-          CALL distrib_level(dim3,lbegin,lend)
+          CALL distrib_level(1,dim3, lbegin,lend)
 
           rval3d=>field(ind)%rval3d
@@ -1295,5 +1295,5 @@
           IF (.NOT. assigned_domain(ind) ) CYCLE
           dim3=size(field(ind)%rval3d,2)
-          CALL distrib_level(dim3,lbegin,lend)
+          CALL distrib_level(1,dim3, lbegin,lend)
           rval3d=>field(ind)%rval3d
           req=>message%request(ind)
@@ -1349,5 +1349,5 @@
 
           dim3=size(field(ind)%rval4d,2)
-          CALL distrib_level(dim3,lbegin,lend)
+          CALL distrib_level(1,dim3, lbegin,lend)
           dim4=size(field(ind)%rval4d,3)
           rval4d=>field(ind)%rval4d
@@ -1411,5 +1411,5 @@
           
           dim3=size(field(ind)%rval4d,2)
-          CALL distrib_level(dim3,lbegin,lend)
+          CALL distrib_level(1,dim3, lbegin,lend)
           dim4=size(field(ind)%rval4d,3)
           rval4d=>field(ind)%rval4d
@@ -1589,5 +1589,5 @@
               dim3=size(rval3d,2)
     
-              CALL distrib_level(dim3,lbegin,lend)
+              CALL distrib_level(1,dim3, lbegin,lend)
               offset=recv%offset*dim3 + (lbegin-1)*recv%size
               CALL trace_start("copy_from_buffer")
@@ -1637,5 +1637,5 @@
 
               dim3=size(rval4d,2)
-              CALL distrib_level(dim3,lbegin,lend)
+              CALL distrib_level(1,dim3, lbegin,lend)
               dim4=size(rval4d,3)
               CALL trace_start("copy_from_buffer")
Index: /codes/icosagcm/trunk/src/sphere/geometry.f90
===================================================================
--- /codes/icosagcm/trunk/src/sphere/geometry.f90	(revision 603)
+++ /codes/icosagcm/trunk/src/sphere/geometry.f90	(revision 604)
@@ -639,4 +639,5 @@
           IF (ii_glo==1 .AND. jj_glo==1) THEN
             le(n+u_ldown)=0
+            le_de(n+u_ldown)=0
             xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
                        
@@ -645,4 +646,5 @@
           IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
             le(n+u_right)=0
+            le_de(n+u_right)=0
             xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
           ENDIF
@@ -650,4 +652,5 @@
           IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
             le(n+u_rup)=0
+            le_de(n+u_rup)=0
             xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
           ENDIF
@@ -655,4 +658,5 @@
           IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
             le(n+u_lup)=0
+            le_de(n+u_lup)=0
             xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
           ENDIF
Index: /codes/icosagcm/trunk/src/time/timeloop_gcm.f90
===================================================================
--- /codes/icosagcm/trunk/src/time/timeloop_gcm.f90	(revision 603)
+++ /codes/icosagcm/trunk/src/time/timeloop_gcm.f90	(revision 604)
@@ -309,5 +309,5 @@
           ! At this point advect_tracer has obtained the halos of u and rhodz,
           ! needed for correct computation of kinetic energy
-          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta, f_hfluxt)
+          IF(diagflux_on) CALL diagflux_energy(adv_over_out, f_phis,f_rhodz,f_theta_rhodz,f_u, f_geopot,f_theta,f_buf_i, f_hfluxt)
 
           IF (check_rhodz) THEN ! check that rhodz is consistent with ps
