Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 150)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 151)
@@ -53,4 +53,5 @@
   SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d)
     USE trace
+    USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
@@ -72,5 +73,5 @@
     ! Compute gradient at triangles solving a linear system
     ! arr = area of triangle joining centroids of hexagons
-    Do l = 1,llm 
+     DO l = ll_begin,ll_end 
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -90,5 +91,5 @@
       
     DO k=1,3
-      DO l =1,llm
+      DO l =ll_begin,ll_end
         DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -103,6 +104,5 @@
 
     !============================================================================================= LIMITING 
-    ! GO TO 120 
-    DO l =1,llm
+    DO l =ll_begin,ll_end
        DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -110,9 +110,4 @@
              maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:))
              maggrd = sqrt(maggrd) 
-!             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
-!                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
-!                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2))
-!             maxq_c = qi(n,l) + maggrd*sqrt(leng(n))
-!             minq_c = qi(n,l) - maggrd*sqrt(leng(n))
              maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n)
              minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n)
@@ -137,4 +132,5 @@
   SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
     USE trace
+    USE omp_para
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
@@ -152,5 +148,5 @@
     
     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -213,4 +209,5 @@
   SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -230,5 +227,5 @@
     ! q(cc)= q0 + gradq.(cc-xyz_i)  with xi centroid of hexagon
     ! ne*hfluxt>0 iff outgoing
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin-1,jj_end+1
           DO i=ii_begin-1,ii_end+1
@@ -265,5 +262,5 @@
 
     ! update q and, if update_mass, update mass
-    DO l = 1,llm
+    DO l = ll_begin,ll_end
        DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 150)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 151)
@@ -9,6 +9,14 @@
   TYPE(t_field),POINTER :: f_cc(:)  ! starting point of backward-trajectory (Miura approach)
   TYPE(t_field),POINTER :: f_one_over_sqrt_leng(:)
- 
+
+  TYPE(t_message) :: req_u, req_wfluxt, req_q, req_rhodz, req_gradq3d
+
   REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
+
+! temporary shared variable for vlz
+  TYPE(t_field),POINTER :: f_dzqw(:)   ! vertical finite difference of q 
+  TYPE(t_field),POINTER :: f_adzqw(:)  ! abs(dzqw)
+  TYPE(t_field),POINTER :: f_dzq(:)    ! limited slope of q
+  TYPE(t_field),POINTER :: f_wq(:)     ! time-integrated flux of q
 
   PUBLIC init_advect_tracer, advect_tracer
@@ -28,5 +36,9 @@
     CALL allocate_field(f_cc,field_u,type_real,llm,3, name='cc')
     CALL allocate_field(f_one_over_sqrt_leng,field_t,type_real, name='one_over_sqrt_leng')
-
+    CALL allocate_field(f_dzqw, field_t, type_real, llm, name='dzqw')
+    CALL allocate_field(f_adzqw, field_t, type_real, llm, name='adzqw')
+    CALL allocate_field(f_dzq, field_t, type_real, llm, name='dzq')
+    CALL allocate_field(f_wq, field_t, type_real, llm+1, name='wq')
+    
     DO ind=1,ndomain
        CALL swap_dimensions(ind)
@@ -56,47 +68,36 @@
     REAL(rstd),POINTER :: hfluxt(:,:), wfluxt(:,:)
     REAL(rstd),POINTER :: rhodz(:,:), u(:,:) 
-    INTEGER :: ind,k
+! temporary shared variable for vlz
+    REAL(rstd),POINTER ::  dzqw(:,:)         ! vertical finite difference of q 
+    REAL(rstd),POINTER ::  adzqw(:,:)        ! abs(dzqw)
+    REAL(rstd),POINTER ::  dzq(:,:)          ! limited slope of q
+    REAL(rstd),POINTER ::  wq(:,:)           ! time-integrated flux of q
+    
+     INTEGER :: ind,k
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_wfluxt,req_i1,req_wfluxt)
+      CALL init_message(f_q,req_i1,req_q)
+      CALL init_message(f_rhodz,req_i1,req_rhodz)
+      CALL init_message(f_gradq3d,req_i1,req_gradq3d)
+    ENDIF
+    
+!$OMP BARRIER
 
     CALL trace_start("advect_tracer") 
 
-    CALL transfert_request(f_u,req_e1_vect)
-!    CALL transfert_request(f_hfluxt,req_e1)  ! BUG : This (unnecessary) transfer makes the computation go wrong
-    CALL transfert_request(f_wfluxt,req_i1)
-    CALL transfert_request(f_q,req_i1)
-    CALL transfert_request(f_rhodz,req_i1)
-	
-!    IF (is_mpi_root) PRINT *, 'Advection scheme'
-
-!    DO ind=1,ndomain
-!       CALL swap_dimensions(ind)
-!       CALL swap_geometry(ind)
-!       normal  = f_normal(ind)
-!       tangent = f_tangent(ind)
-!       cc      = f_cc(ind)
-!       u       = f_u(ind)
-!       q       = f_q(ind)
-!       rhodz   = f_rhodz(ind)
-!       hfluxt  = f_hfluxt(ind) 
-!       wfluxt  = f_wfluxt(ind) 
-!       gradq3d = f_gradq3d(ind)
-!
-!       ! 1/2 vertical transport
-!       DO k = 1, nqtot
-!          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k))
-!       END DO
-!
-!       ! horizontal transport
-!       CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
-!       DO k = 1,nqtot
-!          CALL compute_gradq3d(q(:,:,k),gradq3d)
-!          CALL compute_advect_horiz(k==nqtot,hfluxt,cc,gradq3d, rhodz,q(:,:,k))
-!       END DO
-!
-!       ! 1/2 vertical transport
-!       DO k = 1,nqtot
-!          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k))
-!       END DO
-!    END DO
-
+    CALL send_message(f_u,req_u)
+    CALL send_message(f_wfluxt,req_wfluxt)
+    CALL send_message(f_q,req_q)
+    CALL send_message(f_rhodz,req_rhodz)
+    CALL wait_message(req_u)
+    CALL wait_message(req_wfluxt)
+    CALL wait_message(req_q)
+    CALL wait_message(req_rhodz)
+    
     ! 1/2 vertical transport + back-trajectories
     DO ind=1,ndomain
@@ -110,17 +111,20 @@
        rhodz   = f_rhodz(ind)
        wfluxt  = f_wfluxt(ind) 
+       dzqw    = f_dzqw(ind)
+       adzqw   = f_adzqw(ind)
+       dzq     = f_dzq(ind)
+       wq      = f_wq(ind) 
 
        DO k = 1, nqtot
-          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k),1)
+          CALL vlz(k==nqtot,0.5, wfluxt,rhodz,q(:,:,k),1,dzqw, adzqw, dzq, wq)
        END DO
 
        CALL compute_backward_traj(tangent,normal,u,0.5*dt*itau_adv, cc) 
-    END DO
-
-!    CALL transfert_request(f_q,req_i1)      ! necessary ?
-!    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
+
+    END DO
+
 
     ! horizontal transport - split in two to place transfer of gradq3d
-
+!$OMP BARRIER
     DO k = 1, nqtot
        DO ind=1,ndomain
@@ -133,6 +137,6 @@
        END DO
 
-       CALL transfert_request(f_gradq3d,req_i1)
-
+       CALL send_message(f_gradq3d,req_gradq3d)
+       CALL wait_message(req_gradq3d)
 
 
@@ -149,8 +153,7 @@
     END DO 
     
-!    CALL transfert_request(f_q,req_i1)      ! necessary ?
-!    CALL transfert_request(f_rhodz,req_i1)  ! necessary ?
-    
     ! 1/2 vertical transport
+!$OMP BARRIER
+
     DO ind=1,ndomain
        CALL swap_dimensions(ind)
@@ -159,14 +162,22 @@
        rhodz   = f_rhodz(ind)
        wfluxt  = f_wfluxt(ind) 
+       dzqw    = f_dzqw(ind)
+       adzqw   = f_adzqw(ind)
+       dzq     = f_dzq(ind)
+       wq      = f_wq(ind) 
+
        DO k = 1,nqtot
-          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k),0)
+          CALL vlz(k==nqtot, 0.5,wfluxt,rhodz, q(:,:,k),0, dzqw, adzqw, dzq, wq)
        END DO
+
     END DO
 
     CALL trace_end("advect_tracer")
 
+!$OMP BARRIER
+
   END SUBROUTINE advect_tracer
 
-  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q, halo)
+  SUBROUTINE vlz(update_mass, fac,wfluxt,mass, q, halo, dzqw, adzqw, dzq, wq)
     !
     !     Auteurs:   P.Le Van, F.Hourdin, F.Forget, T. Dubos
@@ -178,4 +189,5 @@
     !    ********************************************************************
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -185,9 +197,11 @@
     INTEGER, INTENT(IN) :: halo
 
-    REAL(rstd) :: dq(iim*jjm,llm), & ! increase of q
-         dzqw(iim*jjm,llm),        & ! vertical finite difference of q 
-         adzqw(iim*jjm,llm),       & ! abs(dzqw)
-         dzq(iim*jjm,llm),         & ! limited slope of q
-         wq(iim*jjm,llm+1)           ! time-integrated flux of q
+! temporary shared variable
+    REAL(rstd),INTENT(INOUT) :: dzqw(iim*jjm,llm),        & ! vertical finite difference of q 
+                                adzqw(iim*jjm,llm),       & ! abs(dzqw)
+                                dzq(iim*jjm,llm),         & ! limited slope of q
+                                wq(iim*jjm,llm+1)           ! time-integrated flux of q
+
+
     REAL(rstd) :: dzqmax, newmass, sigw, qq, w
     INTEGER :: i,ij,l,j
@@ -196,5 +210,6 @@
 
     ! finite difference of q
-    DO l=2,llm
+
+     DO l=ll_beginp1,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
@@ -206,7 +221,11 @@
     ENDDO
 
+!--> flush dzqw, adzqw
+!$OMP BARRIER
+
     ! minmod-limited slope of q
     ! dzq = slope*dz, i.e. the reconstructed q varies by dzq inside level l
-    DO l=2,llm-1
+
+     DO l=ll_beginp1,ll_endm1
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
@@ -223,42 +242,69 @@
     ENDDO
 
+
     ! 0 slope in top and bottom layers
-    DO j=jj_begin-halo,jj_end+halo
-       DO i=ii_begin-halo,ii_end+halo
-          ij=(j-1)*iim+i
-          dzq(ij,1)=0.
+    IF (omp_first) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+           ij=(j-1)*iim+i
+           dzq(ij,1)=0.
+         ENDDO
+      ENDDO
+    ENDIF
+      
+    IF (omp_last) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
           dzq(ij,llm)=0.
-       ENDDO
-    ENDDO
+         ENDDO
+      ENDDO
+    ENDIF
+
+!---> flush dzq
+!$OMP BARRIER  
 
     ! sigw = fraction of mass that leaves level l/l+1
     ! then amount of q leaving level l/l+1 = wq = w * qq
-    DO l = 1,llm-1
+     DO l=ll_beginp1,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
              ij=(j-1)*iim+i
-             w = fac*wfluxt(ij,l+1)
+             w = fac*wfluxt(ij,l)
              IF(w>0.) THEN  ! upward transport, upwind side is at level l 
+                sigw       = w/mass(ij,l-1)
+                qq         = q(ij,l-1)+0.5*(1.-sigw)*dzq(ij,l-1) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
+             ELSE           ! downward transport, upwind side is at level l+1
                 sigw       = w/mass(ij,l)
-                qq         = q(ij,l)+0.5*(1.-sigw)*dzq(ij,l) ! qq = q if sigw=1 , qq = q+dzq/2 if sigw=0
-             ELSE           ! downward transport, upwind side is at level l+1
-                sigw       = w/mass(ij,l+1)
-                qq         = q(ij,l+1)-0.5*(1.+sigw)*dzq(ij,l+1) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
+                qq         = q(ij,l)-0.5*(1.+sigw)*dzq(ij,l) ! qq = q if sigw=-1 , qq = q-dzq/2 if sigw=0                
              ENDIF
-             wq(ij,l+1) = w*qq
+             wq(ij,l) = w*qq
           ENDDO
        ENDDO
     END DO
     ! wq = 0 at top and bottom
-    DO j=jj_begin-halo,jj_end+halo
-       DO i=ii_begin-halo,ii_end+halo
-          ij=(j-1)*iim+i
-          wq(ij,llm+1)=0.
-          wq(ij,1)=0.
-       ENDDO
-    END DO
+    IF (omp_first) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+            ij=(j-1)*iim+i
+            wq(ij,1)=0.
+         ENDDO
+      END DO
+    ENDIF
+    
+    IF (omp_last) THEN
+      DO j=jj_begin-halo,jj_end+halo
+         DO i=ii_begin-halo,ii_end+halo
+            ij=(j-1)*iim+i
+            wq(ij,llm+1)=0.
+         ENDDO
+      END DO
+    ENDIF
+
+! --> flush wq
+!$OMP BARRIER
+
 
     ! update q, mass is updated only after all q's have been updated
-    DO l=1,llm
+    DO l=ll_begin,ll_end
        DO j=jj_begin-halo,jj_end+halo
           DO i=ii_begin-halo,ii_end+halo
Index: codes/icosagcm/trunk/src/caldyn.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn.f90	(revision 150)
+++ codes/icosagcm/trunk/src/caldyn.f90	(revision 151)
@@ -11,5 +11,4 @@
   USE icosa
   USE caldyn_gcm_mod, ONLY : init_caldyn_gcm=>init_caldyn
-  USE caldyn_gcm_opt_mod, ONLY : init_caldyn_gcm_opt=>init_caldyn
   USE caldyn_adv_mod, ONLY : init_caldyn_adv=>init_caldyn
   IMPLICIT NONE
@@ -21,6 +20,4 @@
       CASE('gcm')
         CALL init_caldyn_gcm
-      CASE('gcm_opt')
-        CALL init_caldyn_gcm_opt
       CASE('adv')
         CALL init_caldyn_adv
@@ -37,5 +34,4 @@
   USE icosa
   USE caldyn_gcm_mod, ONLY : caldyn_gcm=>caldyn
-  USE caldyn_gcm_opt_mod, ONLY : caldyn_gcm_opt=>caldyn
   USE caldyn_adv_mod, ONLY : caldyn_adv=>caldyn
   IMPLICIT NONE
@@ -56,7 +52,4 @@
         CALL caldyn_gcm(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
              f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
-      CASE('gcm_opt')
-        CALL caldyn_gcm_opt(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
-             f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
       CASE('adv')
         CALL caldyn_adv(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 150)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 151)
@@ -1,5 +1,5 @@
 MODULE caldyn_gcm_mod
   USE icosa
-
+  USE transfert_mod
   PRIVATE
 
@@ -11,7 +11,17 @@
   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
 
-  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
-
+! temporary shared variable for caldyn
+  TYPE(t_field),POINTER :: f_theta(:)
+  TYPE(t_field),POINTER :: f_pk(:)
+  TYPE(t_field),POINTER :: f_pks(:)
+  TYPE(t_field),POINTER :: f_phi(:)
+  TYPE(t_field),POINTER :: f_divm(:)
+  TYPE(t_field),POINTER :: f_wwuu(:)
+   
   INTEGER :: caldyn_hydrostat, caldyn_conserv
+  
+  TYPE(t_message) :: req_ps, req_theta_rhodz, req_u, req_qu
+
+  PUBLIC init_caldyn, caldyn, write_output_fields,req_ps
 
 CONTAINS
@@ -69,16 +79,23 @@
   IMPLICIT NONE
 
-  CALL allocate_field(f_out_u,field_u,type_real,llm) 
-  CALL allocate_field(f_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_rhodz,field_t,type_real,llm) 
-  CALL allocate_field(f_qu,field_u,type_real,llm) 
-  
-  CALL allocate_field(f_buf_i,field_t,type_real,llm)
-  CALL allocate_field(f_buf_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
-  CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
-  CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
-  CALL allocate_field(f_buf_v,field_z,type_real,llm)
-  CALL allocate_field(f_buf_s,field_t,type_real)
+    CALL allocate_field(f_out_u,field_u,type_real,llm) 
+    CALL allocate_field(f_p,field_t,type_real,llm+1)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm) 
+    CALL allocate_field(f_qu,field_u,type_real,llm) 
+  
+    CALL allocate_field(f_buf_i,field_t,type_real,llm)
+    CALL allocate_field(f_buf_p,field_t,type_real,llm+1) 
+    CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
+    CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
+    CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
+    CALL allocate_field(f_buf_v,field_z,type_real,llm)
+    CALL allocate_field(f_buf_s,field_t,type_real)
+
+    CALL allocate_field(f_theta,field_t,type_real,llm)  ! potential temperature
+    CALL allocate_field(f_pk,field_t,type_real,llm)
+    CALL allocate_field(f_pks,field_t,type_real) ! Exner function
+    CALL allocate_field(f_phi,field_t,type_real,llm)    ! geopotential
+    CALL allocate_field(f_divm,field_t,type_real,llm)  ! mass flux divergence
+    CALL allocate_field(f_wwuu,field_u,type_real,llm+1) 
 
   END SUBROUTINE allocate_caldyn
@@ -92,4 +109,5 @@
     USE mpipara
     USE trace
+    USE omp_para
     IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
@@ -108,11 +126,31 @@
     REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
     REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
+
+! temporary shared variable
+    REAL(rstd),POINTER  :: theta(:,:)  
+    REAL(rstd),POINTER  :: pk(:,:), pks(:)
+    REAL(rstd),POINTER  :: phi(:,:)    
+    REAL(rstd),POINTER  :: divm(:,:) 
+    REAL(rstd),POINTER  :: wwuu(:,:)
+    
+    
     INTEGER :: ind,ij
+    LOGICAL,SAVE :: first=.TRUE.
+!$OMP THREADPRIVATE(first)
+
+    
+    IF (first) THEN
+      first=.FALSE.
+      CALL init_message(f_ps,req_i1,req_ps)
+      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
+      CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_qu,req_e1_scal,req_qu)
+      CALL send_message(f_ps,req_ps) 
+    ENDIF
     
     CALL trace_start("caldyn")
-    CALL transfert_request(f_phis,req_i1) 
-    CALL transfert_request(f_ps,req_i1) 
-    CALL transfert_request(f_theta_rhodz,req_i1) 
-    CALL transfert_request(f_u,req_e1_vect)
+
+    CALL send_message(f_u,req_u)
+    CALL send_message(f_theta_rhodz,req_theta_rhodz) 
 
     SELECT CASE(caldyn_conserv)
@@ -126,10 +164,10 @@
           qu=f_qu(ind)
           u=f_u(ind)
-          !$OMP PARALLEL DEFAULT(SHARED)
+
           CALL compute_pvort(ps, u, p,rhodz,qu)
-          !$OMP END PARALLEL
+
        ENDDO
 
-       CALL transfert_request(f_qu,req_e1_scal)
+       CALL send_message(f_qu,req_qu)
 
        DO ind=1,ndomain
@@ -148,9 +186,14 @@
           u=f_u(ind)
           du=f_du(ind)
-          out_u=f_out_u(ind) 
-          !$OMP PARALLEL DEFAULT(SHARED)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-          !$OMP END PARALLEL
+          out_u=f_out_u(ind)
+          theta = f_theta(ind)
+          pk = f_pk(ind)
+          pks = f_pks(ind)
+          phi = f_phi(ind)  
+          divm = f_divm(ind)
+          wwuu=f_wwuu(ind)
+
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du, &
+                              theta,pk, pks, phi, divm, wwuu)
        ENDDO       
        
@@ -172,9 +215,8 @@
           du=f_du(ind)
           out_u=f_out_u(ind) 
-          !$OMP PARALLEL DEFAULT(SHARED)
+          wwuu=f_wwuu(ind)
           CALL compute_pvort(ps, u, p,rhodz,qu)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-          !$OMP END PARALLEL
+          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
+                              theta,pk, pks, phi, divm, wwuu)
        ENDDO
        
@@ -183,12 +225,22 @@
     END SELECT
 
+!$OMP BARRIER
     IF (write_out) THEN
+
        IF (is_mpi_root) PRINT *,'CALL write_output_fields'
-       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+
+! ---> for openMP test to fix later
+!       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
+!            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
+       CALL writefield("ps",f_ps)
+       CALL writefield("dps",f_dps)
+       CALL writefield("vort",f_qu)
+       CALL writefield("theta",f_theta)
+
     END IF
     
     !    CALL check_mass_conservation(f_ps,f_dps)
     CALL trace_end("caldyn")
+!$OMP BARRIER
     
 END SUBROUTINE caldyn
@@ -199,4 +251,5 @@
   USE exner_mod
   USE trace
+  USE omp_para
   IMPLICIT NONE
   REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
@@ -208,19 +261,14 @@
   INTEGER :: i,j,ij,l
   REAL(rstd) :: etav,hv
-  REAL(rstd),ALLOCATABLE,SAVE :: qv(:,:)     ! potential velocity  
-  
-  LOGICAL,SAVE :: first=.TRUE.
-  !$OMP THREADPRIVATE(first)
-
+  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
+  
 
   CALL trace_start("compute_pvort")  
-  !$OMP BARRIER      
-  !$OMP MASTER  
-  !    IF (first) THEN
-  ALLOCATE(qv(2*iim*jjm,llm))     ! potential velocity  
-  
+
+  CALL wait_message(req_ps)   
 !!! Compute pressure
-  DO    l    = 1, llm+1
-     !$OMP DO
+   DO    l    = ll_begin, ll_endp1
+    CALL test_message(req_u) 
+
      DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -230,8 +278,10 @@
      ENDDO
   ENDDO
-  
+
+!$OMP BARRIER  
+
 !!! Compute mass
-  DO l = 1, llm
-     !$OMP DO
+   DO l = ll_begin,ll_end
+     CALL test_message(req_u) 
      DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -241,15 +291,18 @@
      ENDDO
   ENDDO
+
+  CALL wait_message(req_u)   
   
 !!! Compute shallow-water potential vorticity
-  DO l = 1,llm
-!$OMP DO
+  DO l = ll_begin,ll_end
+    CALL test_message(req_theta_rhodz) 
+
     DO j=jj_begin-1,jj_end+1
        DO i=ii_begin-1,ii_end+1
           ij=(j-1)*iim+i
           
-          etav= 1./Av(ij+z_up)*(  ne(ij,rup)        * u(ij+u_rup,l)        * de(ij+u_rup)         &
-                                + ne(ij+t_rup,left) * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
-                                - ne(ij,lup)        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
+          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
+                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
+                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
 
           hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
@@ -259,7 +312,7 @@
           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
           
-          etav = 1./Av(ij+z_down)*(  ne(ij,ldown)         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
-                                   + ne(ij+t_ldown,right) * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
-                                   - ne(ij,rdown)         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
+          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
+                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
+                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
        
           hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
@@ -283,18 +336,16 @@
    ENDDO
 
-!!$OMP BARRIER
-!!$OMP MASTER  
-    DEALLOCATE(qv)       ! potential velocity  
-!!$OMP END MASTER
-!!$OMP BARRIER   
-    CALL trace_end("compute_pvort")
+   CALL trace_end("compute_pvort")
                                                        
   END SUBROUTINE compute_pvort
-  
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
+
+  
+  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du,   &
+                            theta,pk, pks, phi, divm, wwuu) 
   USE icosa
   USE disvert_mod
   USE exner_mod
   USE trace
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
@@ -310,36 +361,26 @@
     REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
     REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
-    
-    REAL(rstd),ALLOCATABLE,SAVE :: theta(:,:)  ! potential temperature
-    REAL(rstd),ALLOCATABLE,SAVE :: pk(:,:), pks(:) ! Exner function
-    REAL(rstd),ALLOCATABLE,SAVE :: alpha(:,:), beta(:,:)
-    REAL(rstd),ALLOCATABLE,SAVE :: phi(:,:)    ! geopotential
-    REAL(rstd),ALLOCATABLE,SAVE :: Ftheta(:,:) ! theta flux
-    REAL(rstd),ALLOCATABLE,SAVE :: divm(:,:)  ! mass flux divergence
-    REAL(rstd),ALLOCATABLE,SAVE :: berni(:,:)  ! Bernouilli function
+
+! temporary variable
+    REAL(rstd),INTENT(INOUT) :: theta(iim*jjm,llm)  ! potential temperature
+    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm), pks(iim*jjm) ! Exner function
+    REAL(rstd),INTENT(INOUT) :: phi(iim*jjm,llm)    ! geopotential
+    REAL(rstd),INTENT(INOUT) :: divm(iim*jjm,llm)  ! mass flux divergence
+    REAL(rstd),INTENT(INOUT) :: wwuu(iim*3*jjm,llm+1)
+
+    
+    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
+    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function
     
     INTEGER :: i,j,ij,l
     REAL(rstd) :: ww,uu
 
-    LOGICAL,SAVE :: first=.TRUE.
-!$OMP THREADPRIVATE(first)
-    
     CALL trace_start("compute_caldyn")
-!$OMP BARRIER      
-!$OMP MASTER  
-!    IF (first) THEN
-    ALLOCATE(theta(iim*jjm,llm))    ! potential temperature
-    ALLOCATE(pk(iim*jjm,llm))       ! Exner function
-    ALLOCATE(pks(iim*jjm))
-    ALLOCATE(alpha(iim*jjm,llm))
-    ALLOCATE(beta(iim*jjm,llm))
-    ALLOCATE(phi(iim*jjm,llm))      ! geopotential
-    ALLOCATE(Ftheta(3*iim*jjm,llm)) ! theta flux   
-    ALLOCATE(divm(iim*jjm,llm))    ! mass flux divvergence
-    ALLOCATE(berni(iim*jjm,llm))    ! bernouilli term
+
+    CALL wait_message(req_theta_rhodz) 
 
 !!! Compute theta
-   DO l = 1, llm
-!$OMP DO
+   DO l = ll_begin, ll_end
+      IF (caldyn_conserv==energy) CALL test_message(req_qu) 
       DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
@@ -350,6 +391,7 @@
    ENDDO
 
-  DO l = 1, llm
+  DO l = ll_begin, ll_end
 !!!  Compute mass and theta fluxes
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
     DO j=jj_begin-1,jj_end+1
       DO i=ii_begin-1,ii_end+1
@@ -370,26 +412,28 @@
         ij=(j-1)*iim+i  
         ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        divm(ij,l)= 1./Ai(ij)*(ne(ij,right)*hflux(ij+u_right,l)   +  &
-                                ne(ij,rup)*hflux(ij+u_rup,l)       +  &  
-                                ne(ij,lup)*hflux(ij+u_lup,l)       +  &  
-                                ne(ij,left)*hflux(ij+u_left,l)     +  &  
-                                ne(ij,ldown)*hflux(ij+u_ldown,l)   +  &  
-                                ne(ij,rdown)*hflux(ij+u_rdown,l))    
+        divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
+                                ne_rup*hflux(ij+u_rup,l)       +  &  
+                                ne_lup*hflux(ij+u_lup,l)       +  &  
+                                ne_left*hflux(ij+u_left,l)     +  &  
+                                ne_ldown*hflux(ij+u_ldown,l)   +  &  
+                                ne_rdown*hflux(ij+u_rdown,l))    
 
         ! signe ? attention d (rho theta dz)
         ! dtheta_rhodz =  -div(flux.theta)
-        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
-                                 ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
-                                 ne(ij,lup)*Ftheta(ij+u_lup,l)        +  &  
-                                 ne(ij,left)*Ftheta(ij+u_left,l)      +  &  
-                                 ne(ij,ldown)*Ftheta(ij+u_ldown,l)    +  &  
-                                 ne(ij,rdown)*Ftheta(ij+u_rdown,l))
+        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
+                                 ne_rup*Ftheta(ij+u_rup,l)        +  &  
+                                 ne_lup*Ftheta(ij+u_lup,l)        +  &  
+                                 ne_left*Ftheta(ij+u_left,l)      +  &  
+                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
+                                 ne_rdown*Ftheta(ij+u_rdown,l))
       ENDDO
     ENDDO
   ENDDO
-   
+
+!$OMP BARRIER   
 !!! cumulate mass flux divergence from top to bottom
   DO  l = llm-1, 1, -1
-!$OMP DO
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
+!$OMP DO SCHEDULE(STATIC) 
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -399,8 +443,12 @@
     ENDDO
   ENDDO
+
+! IMPLICIT FLUSH on divm
+!!!!!!!!!!!!!!!!!!!!!!!!!
   
 !!! Compute vertical mass flux
-  DO l = 1,llm-1
-!$OMP DO
+!  DO l = 2,llm
+  DO l=ll_beginp1,ll_end
+    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -408,5 +456,5 @@
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
-        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
+        wflux( ij, l ) = divm( ij, l ) - bp(l) * divm( ij, 1 )
       ENDDO
     ENDDO
@@ -414,15 +462,24 @@
 
 ! compute dps, set vertical mass flux at the surface to 0
-!$OMP DO
-  DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      wflux(ij,1)  = 0.
-      wflux(ij,llm+1)  = 0.
-      ! dps/dt = -int(div flux)dz
-      dps(ij)=-divm(ij,1) * g 
-    ENDDO
-  ENDDO
-
+  IF (omp_first) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wflux(ij,1)  = 0.
+        ! dps/dt = -int(div flux)dz
+        dps(ij)=-divm(ij,1) * g 
+      ENDDO
+    ENDDO
+  ENDIF
+
+
+  IF (omp_last) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wflux(ij,llm+1)  = 0.
+      ENDDO
+    ENDDO
+  ENDIF
 !!! Compute potential vorticity (Coriolis) contribution to du
 
@@ -430,6 +487,7 @@
     CASE(energy) ! energy-conserving TRiSK
 
-       DO l=1,llm
-          !$OMP DO
+       CALL wait_message(req_qu)
+
+        DO l=ll_begin,ll_end
           DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -479,6 +537,5 @@
     CASE(enstrophy) ! enstrophy-conserving TRiSK
   
-       DO l=1,llm
-          !$OMP DO
+        DO l=ll_begin,ll_end
           DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -531,19 +588,53 @@
   
 !!! Compute Exner function
-!    PRINT *, 'Computing Exner'
-    CALL compute_exner(ps,p,pks,pk,1) 
+!    CALL compute_exner(ps,p,pks,pk,1) 
+! replaced in source
+    IF (omp_first) THEN    
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             ij=(j-1)*iim+i
+             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+          ENDDO
+       ENDDO
+     ENDIF
+     
+       ! 3D : pk
+     DO l = ll_begin, ll_end
+        DO j=jj_begin-1,jj_end+1
+           DO i=ii_begin-1,ii_end+1
+              ij=(j-1)*iim+i
+              pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
+           ENDDO
+        ENDDO
+     ENDDO
+
+!---> flush pk,pks, theta
+!$OMP BARRIER
 
 !!! Compute geopotential
 
   ! for first layer
-!$OMP DO
-   DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-       ij=(j-1)*iim+i
-       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+  IF (omp_first) THEN
+    DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+        ij=(j-1)*iim+i
+        phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
+      ENDDO
+    ENDDO
+  ENDIF
+!!-> implicit flush on phi(:,1)
+          
+  ! for other layers
+   DO l = ll_beginp1, ll_end
+     DO j=jj_begin-1,jj_end+1
+       DO i=ii_begin-1,ii_end+1
+         ij=(j-1)*iim+i
+         phi(ij,l) =  0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
+                          * (  pk(ij,l-1) -  pk(ij,l)    )
+       ENDDO
      ENDDO
    ENDDO
-          
-  ! for other layers
+
+!$OMP BARRIER
    DO l = 2, llm
 !$OMP DO
@@ -551,14 +642,12 @@
        DO i=ii_begin-1,ii_end+1
          ij=(j-1)*iim+i
-         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                       * (  pk(ij,l-1) -  pk(ij,l)    )
+         phi(ij,l) = phi(ij,l)+ phi(ij,l-1)
        ENDDO
      ENDDO
    ENDDO
-
-   
+! --> IMPLICIT FLUSH on phi
+      
 !!! Compute bernouilli term = Kinetic Energy + geopotential
-  DO l=1,llm
-!$OMP DO
+   DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -576,84 +665,100 @@
     ENDDO
   ENDDO
-  
+
  
 !!! gradients of Bernoulli and Exner functions
-  DO l=1,llm
-!$OMP DO
+
+   DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
         
-        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
-                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
-                             *( ne(ij,right)*pk(ij,l)+ne(ij+t_right,left)*pk(ij+t_right,l)) &
-                        + ne(ij,right)*berni(ij,l)+ne(ij+t_right,left)*berni(ij+t_right,l) )
+        du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
+                          0.5*(theta(ij,l)+theta(ij+t_right,l))                &
+                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l))    &
+                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
         
-        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
         
-        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
-                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
-                             *( ne(ij,lup)*pk(ij,l)+ne(ij+t_lup,rdown)*pk(ij+t_lup,l))    &
-                        + ne(ij,lup)*berni(ij,l)+ne(ij+t_lup,rdown)*berni(ij+t_lup,l) )
-        
-        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
+        du(ij+u_lup,l) = du(ij+u_lup,l) +  1/de(ij+u_lup) * (                  &
+                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                  &
+                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))       &
+                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
                 
-        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
-                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
-                             *( ne(ij,ldown)*pk(ij,l)+ne(ij+t_ldown,rup)*pk(ij+t_ldown,l))    &
-                        + ne(ij,ldown)*berni(ij,l)+ne(ij+t_ldown,rup)*berni(ij+t_ldown,l) )
-        
-        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
+        du(ij+u_ldown,l) = du(ij+u_ldown,l) + 1/de(ij+u_ldown) * (             &
+                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                &
+                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))     &
+                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
 
       ENDDO
     ENDDO
   ENDDO
- 
-!!! contributions of vertical advection to du, dtheta
- 
-  DO l=1,llm-1
-!$OMP DO
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-         ! ww>0 <=> upward transport
-
-        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
-        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
-        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
-        uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
-        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
-        du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
-        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
-        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
-        du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
-        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
-        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
-        du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))
-
+
+   
+  DO l=ll_begin,ll_endm1
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  DO l=ll_beginp1,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
+      ENDDO
+    ENDDO
+  ENDDO
+  
+  IF (omp_first) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,1)=0   
+        wwuu(ij+u_lup,1)=0   
+        wwuu(ij+u_ldown,1)=0
+      ENDDO
+    ENDDO
+  ENDIF   
+
+  IF (omp_last) THEN
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,llm+1)=0   
+        wwuu(ij+u_lup,llm+1)=0   
+        wwuu(ij+u_ldown,llm+1)=0
+      ENDDO
+    ENDDO
+  ENDIF   
+  
+  DO l=ll_beginp1,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
+        wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
+        wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
+       ENDDO
+     ENDDO
+   ENDDO
+
+ !--> flush wwuu  
+ !$OMP BARRIER
+
+  DO l=ll_begin,ll_end
+    DO j=jj_begin,jj_end
+      DO i=ii_begin,ii_end
+        ij=(j-1)*iim+i
+        du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
+        du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
+        du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
       ENDDO
     ENDDO
   ENDDO      
-   
-!!$OMP BARRIER
-!!$OMP MASTER  
-    DEALLOCATE(theta)   ! potential temperature
-    DEALLOCATE(pk)      ! Exner function
-    DEALLOCATE(pks)
-    DEALLOCATE(alpha)
-    DEALLOCATE(beta)
-    DEALLOCATE(phi)     ! geopotential
-    DEALLOCATE(Ftheta)  ! theta flux   
-    DEALLOCATE(divm)    ! mass flux divergence
-    DEALLOCATE(berni)   ! bernouilli term
-!!$OMP END MASTER
-!!$OMP BARRIER                                                      
-    CALL trace_end("compute_caldyn")
+  
+  CALL trace_end("compute_caldyn")
 
   END SUBROUTINE compute_caldyn
@@ -709,4 +814,5 @@
     USE write_field
     USE vertical_interp_mod
+    USE wind_mod
     TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
          f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
@@ -722,15 +828,15 @@
     
     CALL writefield("ps",f_ps)
-!    CALL writefield("dps",f_dps)
-!    CALL writefield("phis",f_phis)
-!    CALL vorticity(f_u,f_buf_v)
-!    CALL writefield("vort",f_buf_v)
-
-!    CALL w_omega(f_ps, f_u, f_buf_i)
-!    CALL writefield('omega', f_buf_i)
-!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-!    ENDIF
+    CALL writefield("dps",f_dps)
+    CALL writefield("phis",f_phis)
+    CALL vorticity(f_u,f_buf_v)
+    CALL writefield("vort",f_buf_v)
+
+    CALL w_omega(f_ps, f_u, f_buf_i)
+    CALL writefield('omega', f_buf_i)
+    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
+      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
+      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
+    ENDIF
     
     ! Temperature
@@ -754,5 +860,5 @@
    
     ! velocity components
-    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
+    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
     CALL writefield("ulon",f_buf1_i)
     CALL writefield("ulat",f_buf2_i)
@@ -767,8 +873,8 @@
     ! geopotential
     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-!    CALL writefield("p",f_buf_p)
-!    CALL writefield("phi",f_buf_i)
-     CALL writefield("theta",f_buf1_i) ! potential temperature
-!    CALL writefield("pk",f_buf2_i) ! Exner pressure
+    CALL writefield("p",f_buf_p)
+    CALL writefield("phi",f_buf_i)
+    CALL writefield("theta",f_buf1_i) ! potential temperature
+    CALL writefield("pk",f_buf2_i) ! Exner pressure
   
   
@@ -806,23 +912,4 @@
   END SUBROUTINE thetarhodz2geopot
   
-  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
-    USE field_mod
-    USE wind_mod
-    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
-         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
-    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
-    INTEGER :: ind
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       u3d=f_u3d(ind)
-       CALL compute_wind_centered(u,u3d)
-       ulon=f_ulon(ind)
-       ulat=f_ulat(ind)
-       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
-    END DO
-  END SUBROUTINE un2ulonlat
-
   SUBROUTINE Tv2T(f_Tv, f_q, f_T)
   USE icosa
Index: codes/icosagcm/trunk/src/caldyn_gcm_opt.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm_opt.f90	(revision 150)
+++ 	(revision )
@@ -1,791 +1,0 @@
-MODULE caldyn_gcm_opt_mod
-  USE icosa
-
-  PRIVATE
-
-  INTEGER, PARAMETER :: energy=1, enstrophy=2
-  TYPE(t_field),POINTER :: f_out_u(:), f_p(:), f_rhodz(:), f_qu(:)
-  REAL(rstd),POINTER :: out_u(:,:), p(:,:), rhodz(:,:), qu(:,:)
-
-  TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
-  TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
-
-  PUBLIC init_caldyn, caldyn, write_output_fields
-
-  INTEGER :: caldyn_hydrostat, caldyn_conserv
-
-CONTAINS
-  
-  SUBROUTINE init_caldyn
-    USE icosa
-    USE exner_mod
-    USE mpipara
-    IMPLICIT NONE
-    CHARACTER(len=255) :: def
-  
-    def='enstrophy'
-    CALL getin('caldyn_conserv',def)
-    SELECT CASE(TRIM(def))
-    CASE('energy')
-       caldyn_conserv=energy
-    CASE('enstrophy')
-       caldyn_conserv=enstrophy
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', TRIM(def),'> options are <energy>, <enstrophy>'
-       STOP
-    END SELECT
-    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
-
-    def='direct'
-    CALL getin('caldyn_exner',def)
-    SELECT CASE(TRIM(def))
-    CASE('lmdz')
-       caldyn_exner=lmdz
-    CASE('direct')
-       caldyn_exner=direct
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_exner : <', TRIM(def),'> options are <lmdz>, <direct>'
-       STOP
-    END SELECT
-
-    def='direct'
-    CALL getin('caldyn_hydrostat',def)
-    SELECT CASE(TRIM(def))
-    CASE('lmdz')
-       caldyn_hydrostat=lmdz
-    CASE('direct')
-       caldyn_hydrostat=direct
-    CASE DEFAULT
-       IF (is_mpi_root) PRINT*,'Bad selector for variable caldyn_hydrostat : <', TRIM(def),'> options are <lmdz>, <direct>'
-       STOP
-    END SELECT
-    
-    CALL allocate_caldyn
-  
-  END SUBROUTINE init_caldyn
- 
-  SUBROUTINE allocate_caldyn
-  USE icosa
-  IMPLICIT NONE
-
-  CALL allocate_field(f_out_u,field_u,type_real,llm) 
-  CALL allocate_field(f_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_rhodz,field_t,type_real,llm) 
-  CALL allocate_field(f_qu,field_u,type_real,llm) 
-  
-  CALL allocate_field(f_buf_i,field_t,type_real,llm)
-  CALL allocate_field(f_buf_p,field_t,type_real,llm+1)
-  CALL allocate_field(f_buf_u3d,field_t,type_real,3,llm)  ! 3D vel at cell centers
-  CALL allocate_field(f_buf_ulon,field_t,type_real,llm)
-  CALL allocate_field(f_buf_ulat,field_t,type_real,llm)
-  CALL allocate_field(f_buf_v,field_z,type_real,llm)
-  CALL allocate_field(f_buf_s,field_t,type_real)
-
-  END SUBROUTINE allocate_caldyn
-   
-  SUBROUTINE caldyn(write_out,f_phis, f_ps, f_theta_rhodz, f_u, f_q, &
-       f_hflux, f_wflux, f_dps, f_dtheta_rhodz, f_du)
-    USE icosa
-    USE vorticity_mod
-    USE kinetic_mod
-    USE theta2theta_rhodz_mod
-    USE mpipara
-    USE trace
-    IMPLICIT NONE
-    LOGICAL,INTENT(IN)    :: write_out
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
-    TYPE(t_field),POINTER :: f_dps(:)
-    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
-    TYPE(t_field),POINTER :: f_du(:)
-    
-    REAL(rstd),POINTER :: phis(:), ps(:), dps(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:), dtheta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:), du(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: p(:,:), rhodz(:,:), qu(:,:)
-    INTEGER :: ind,ij
-    
-    CALL trace_start("caldyn")
-    CALL transfert_request(f_ps,req_i1) 
-    CALL transfert_request(f_theta_rhodz,req_i1) 
-    CALL transfert_request(f_u,req_e1_vect)
-
-    SELECT CASE(caldyn_conserv)
-    CASE(energy) ! energy-conserving
-       DO ind=1,ndomain
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          ps=f_ps(ind)
-          rhodz=f_rhodz(ind)
-          p=f_p(ind)
-          qu=f_qu(ind)
-          u=f_u(ind)
-          CALL compute_pvort(ps, u, p,rhodz,qu)
-       ENDDO
-
-       CALL transfert_request(f_qu,req_e1_scal)
-
-       DO ind=1,ndomain
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          phis=f_phis(ind)
-          hflux=f_hflux(ind)
-          wflux=f_wflux(ind)
-          ps=f_ps(ind)
-          dps=f_dps(ind)
-          theta_rhodz=f_theta_rhodz(ind)
-          dtheta_rhodz=f_dtheta_rhodz(ind)
-          rhodz=f_rhodz(ind)
-          p=f_p(ind)
-          qu=f_qu(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          out_u=f_out_u(ind) 
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-       ENDDO       
-       
-    CASE(enstrophy) ! enstrophy-conserving
-       DO ind=1,ndomain
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          phis=f_phis(ind)
-          ps=f_ps(ind)
-          dps=f_dps(ind)
-          hflux=f_hflux(ind)
-          wflux=f_wflux(ind)
-          theta_rhodz=f_theta_rhodz(ind)
-          dtheta_rhodz=f_dtheta_rhodz(ind)
-          rhodz=f_rhodz(ind)
-          p=f_p(ind)
-          qu=f_qu(ind)
-          u=f_u(ind)
-          du=f_du(ind)
-          out_u=f_out_u(ind) 
-          CALL compute_pvort(ps, u, p,rhodz,qu)
-          CALL compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, &
-               hflux, wflux, dps, dtheta_rhodz, du)
-       ENDDO
-       
-    CASE DEFAULT
-       STOP
-    END SELECT
-
-    IF (write_out) THEN
-       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
-       CALL write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-            f_buf_i, f_buf_v, f_buf_u3d, f_buf_ulon, f_buf_ulat, f_buf_s, f_buf_p)
-    END IF
-    
-    !    CALL check_mass_conservation(f_ps,f_dps)
-    CALL trace_end("caldyn")
-    
-END SUBROUTINE caldyn
-    
-SUBROUTINE compute_pvort(ps, u, p,rhodz,qu)
-  USE icosa
-  USE disvert_mod
-  USE exner_mod
-  USE trace
-  IMPLICIT NONE
-  REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-  REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-  REAL(rstd),INTENT(OUT) :: p(iim*jjm,llm+1)
-  REAL(rstd),INTENT(OUT) :: rhodz(iim*jjm,llm)
-  REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
-  
-  INTEGER :: i,j,ij,l
-  REAL(rstd) :: etav,hv
-  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
-  
-
-  CALL trace_start("compute_pvort")  
-  
-!!! Compute pressure
-  DO    l    = 1, llm+1
-     DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-           ij=(j-1)*iim+i
-           p(ij,l) = ap(l) + bp(l) * ps(ij)
-        ENDDO
-     ENDDO
-  ENDDO
-  
-!!! Compute mass
-  DO l = 1, llm
-     DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-           ij=(j-1)*iim+i
-           rhodz(ij,l) = ( p(ij,l) - p(ij,l+1) )/g
-        ENDDO
-     ENDDO
-  ENDDO
-  
-!!! Compute shallow-water potential vorticity
-  DO l = 1,llm
-    DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          
-          etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
-                                + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
-                                - ne_lup        * u(ij+u_lup,l)        * de(ij+u_lup) )                               
-
-          hv =  Riv2(ij,vup)          * rhodz(ij,l)            &
-              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-     
-          qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
-          
-          etav = 1./Av(ij+z_down)*(  ne_ldown         * u(ij+u_ldown,l)          * de(ij+u_ldown)          &
-                                   + ne_right * u(ij+t_ldown+u_right,l)  * de(ij+t_ldown+u_right)  &
-                                   - ne_rdown         * u(ij+u_rdown,l)          * de(ij+u_rdown) )
-       
-          hv =  Riv2(ij,vdown)        * rhodz(ij,l)          &
-              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-                       
-          qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
-
-        ENDDO
-      ENDDO
-
-      DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
-            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
-            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
-            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
-         END DO
-      END DO
-      
-   ENDDO
-
-    CALL trace_end("compute_pvort")
-                                                       
-  END SUBROUTINE compute_pvort
-  
-  SUBROUTINE compute_caldyn(ps, u, p,rhodz,qu, phis, theta_rhodz, hflux, wflux, dps, dtheta_rhodz, du)
-  USE icosa
-  USE disvert_mod
-  USE exner_mod
-  USE trace
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN)  :: phis(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: ps(iim*jjm)
-    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)  :: qu(iim*3*jjm,llm)
-
-    REAL(rstd),INTENT(OUT) :: du(iim*3*jjm,llm), hflux(iim*3*jjm,llm) ! hflux in kg/s
-    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: dps(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: wflux(iim*jjm,llm+1) ! vertical mass flux (kg/m2/s)
-    
-    REAL(rstd) :: theta(iim*jjm,llm)  ! potential temperature
-    REAL(rstd) :: pk(iim*jjm,llm), pks(iim*jjm) ! Exner function
-    REAL(rstd) :: alpha(iim*jjm,llm), beta(iim*jjm,llm)
-    REAL(rstd) :: phi(iim*jjm,llm)    ! geopotential
-    REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! theta flux
-    REAL(rstd) :: divm(iim*jjm,llm)  ! mass flux divergence
-    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernouilli function
-    
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: ww,uu
-
-    LOGICAL,SAVE :: first=.TRUE.
-    
-    CALL trace_start("compute_caldyn")
-
-
-!!! Compute theta
-   DO l = 1, llm
-      DO j=jj_begin-1,jj_end+1
-         DO i=ii_begin-1,ii_end+1
-            ij=(j-1)*iim+i
-            theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-         ENDDO
-      ENDDO
-   ENDDO
-
-  DO l = 1, llm
-!!!  Compute mass and theta fluxes
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
-        hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-        hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-        hflux(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
-
-        Ftheta(ij+u_right,l)=0.5*(theta(ij,l)+theta(ij+t_right,l))*hflux(ij+u_right,l)
-        Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
-        Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
-      ENDDO
-    ENDDO
-
-!!! compute horizontal divergence of fluxes
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i  
-        ! divm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-        divm(ij,l)= 1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
-                                ne_rup*hflux(ij+u_rup,l)       +  &  
-                                ne_lup*hflux(ij+u_lup,l)       +  &  
-                                ne_left*hflux(ij+u_left,l)     +  &  
-                                ne_ldown*hflux(ij+u_ldown,l)   +  &  
-                                ne_rdown*hflux(ij+u_rdown,l))    
-
-        ! signe ? attention d (rho theta dz)
-        ! dtheta_rhodz =  -div(flux.theta)
-        dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne_right*Ftheta(ij+u_right,l)    +  &
-                                 ne_rup*Ftheta(ij+u_rup,l)        +  &  
-                                 ne_lup*Ftheta(ij+u_lup,l)        +  &  
-                                 ne_left*Ftheta(ij+u_left,l)      +  &  
-                                 ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
-                                 ne_rdown*Ftheta(ij+u_rdown,l))
-      ENDDO
-    ENDDO
-  ENDDO
-   
-!!! cumulate mass flux divergence from top to bottom
-  DO  l = llm-1, 1, -1
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        divm(ij,l) = divm(ij,l) + divm(ij,l+1)
-      ENDDO
-    ENDDO
-  ENDDO
-  
-!!! Compute vertical mass flux
-  DO l = 1,llm-1
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
-        ! => w>0 for upward transport
-        wflux( ij, l+1 ) = divm( ij, l+1 ) - bp(l+1) * divm( ij, 1 )
-      ENDDO
-    ENDDO
-  ENDDO
-
-! compute dps, set vertical mass flux at the surface to 0
-  DO j=jj_begin,jj_end
-    DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      wflux(ij,1)  = 0.
-      wflux(ij,llm+1)  = 0.
-      ! dps/dt = -int(div flux)dz
-      dps(ij)=-divm(ij,1) * g 
-    ENDDO
-  ENDDO
-
-!!! Compute potential vorticity (Coriolis) contribution to du
-
-    SELECT CASE(caldyn_conserv)
-    CASE(energy) ! energy-conserving TRiSK
-
-       DO l=1,llm
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-
-                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
-                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
-                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
-                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
-                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
-                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
-                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
-                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
-                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
-                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
-                du(ij+u_right,l) = .5*uu/de(ij+u_right)
-                
-                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
-                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
-                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
-                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
-                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
-                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
-                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
-                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
-                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
-                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
-                du(ij+u_lup,l) = .5*uu/de(ij+u_lup)
-
-                
-                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
-                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
-                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
-                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
-                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
-                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
-                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
-                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
-                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
-                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
-                du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
-                
-             ENDDO
-          ENDDO
-       ENDDO
-       
-    CASE(enstrophy) ! enstrophy-conserving TRiSK
-  
-       DO l=1,llm
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-
-                uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
-                     wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
-                     wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
-                     wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
-                     wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
-                     wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
-                     wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
-                     wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
-                     wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
-                     wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
-                du(ij+u_right,l) = qu(ij+u_right,l)*uu/de(ij+u_right)
-                
-                
-                uu = wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
-                     wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
-                     wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
-                     wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
-                     wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
-                     wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
-                     wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
-                     wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
-                     wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
-                     wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
-                du(ij+u_lup,l) = qu(ij+u_lup,l)*uu/de(ij+u_lup)
-
-                uu = wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
-                     wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
-                     wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
-                     wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
-                     wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
-                     wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
-                     wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
-                     wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
-                     wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
-                     wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
-                du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
-     
-             ENDDO
-          ENDDO
-       ENDDO
-       
-    CASE DEFAULT
-       STOP
-    END SELECT
-  
-!!! Compute Exner function
-!    PRINT *, 'Computing Exner'
-    CALL compute_exner(ps,p,pks,pk,1) 
-
-!!! Compute geopotential
-
-  ! for first layer
-   DO j=jj_begin-1,jj_end+1
-     DO i=ii_begin-1,ii_end+1
-       ij=(j-1)*iim+i
-       phi( ij,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
-     ENDDO
-   ENDDO
-          
-  ! for other layers
-   DO l = 2, llm
-     DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-         ij=(j-1)*iim+i
-         phi(ij,l) = phi(ij,l-1) + 0.5 * ( theta(ij,l)  + theta(ij,l-1) )  & 
-                                       * (  pk(ij,l-1) -  pk(ij,l)    )
-       ENDDO
-     ENDDO
-   ENDDO
-
-   
-!!! Compute bernouilli term = Kinetic Energy + geopotential
-  DO l=1,llm
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-
-        berni(ij,l) =  phi(ij,l) &                                                          
-                     + 1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +    &
-                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +       &
-                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-       
-      ENDDO
-    ENDDO
-  ENDDO
-  
- 
-!!! gradients of Bernoulli and Exner functions
-  DO l=1,llm
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        
-        out_u(ij+u_right,l)=  1/de(ij+u_right) * (                              &
-                          0.5*(theta(ij,l)+theta(ij+t_right,l))                             &
-                             *( ne_right*pk(ij,l)+ne_left*pk(ij+t_right,l)) &
-                        + ne_right*berni(ij,l)+ne_left*berni(ij+t_right,l) )
-        
-        du(ij+u_right,l) = du(ij+u_right,l) + out_u(ij+u_right,l)
-        
-        out_u(ij+u_lup,l)=  1/de(ij+u_lup) * (                                  &
-                          0.5*(theta(ij,l)+theta(ij+t_lup,l))                             &
-                             *( ne_lup*pk(ij,l)+ne_rdown*pk(ij+t_lup,l))    &
-                        + ne_lup*berni(ij,l)+ne_rdown*berni(ij+t_lup,l) )
-        
-        du(ij+u_lup,l) = du(ij+u_lup,l) + out_u(ij+u_lup,l)
-                
-        out_u(ij+u_ldown,l)= 1/de(ij+u_ldown) * (                                &
-                          0.5*(theta(ij,l)+theta(ij+t_ldown,l))                               &
-                             *( ne_ldown*pk(ij,l)+ne_rup*pk(ij+t_ldown,l))    &
-                        + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
-        
-        du(ij+u_ldown,l) = du(ij+u_ldown,l) + out_u(ij+u_ldown,l)
-
-      ENDDO
-    ENDDO
-  ENDDO
- 
-!!! contributions of vertical advection to du, dtheta
- 
-  DO l=1,llm-1
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-         ! ww>0 <=> upward transport
-
-        ww            =   0.5 * wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
-        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
-        dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_right,l+1))
-        uu  = u(ij+u_right,l+1) - u(ij+u_right,l)  
-        du(ij+u_right, l )   = du(ij+u_right,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_right,l)))
-        du(ij+u_right, l+1 ) = du(ij+u_right,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_right,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_lup,l+1))
-        uu  = u(ij+u_lup,l+1) - u(ij+u_lup,l)  
-        du(ij+u_lup, l )   = du(ij+u_lup,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l)))
-        du(ij+u_lup, l+1 ) = du(ij+u_lup,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_lup,l+1)))
-
-        ww  = 0.5 * ( wflux(ij,l+1) + wflux(ij+t_ldown,l+1))
-        uu  = u(ij+u_ldown,l+1) - u(ij+u_ldown,l)  
-        du(ij+u_ldown, l )   = du(ij+u_ldown,l)   - 0.5 * ww * uu / (0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l)))
-        du(ij+u_ldown, l+1 ) = du(ij+u_ldown,l+1) - 0.5 * ww * uu / (0.5*(rhodz(ij,l+1)+rhodz(ij+t_ldown,l+1)))
-
-      ENDDO
-    ENDDO
-  ENDDO      
-   
-
-    CALL trace_end("compute_caldyn")
-
-  END SUBROUTINE compute_caldyn
-
-!-------------------------------- Diagnostics ----------------------------
-
-  SUBROUTINE check_mass_conservation(f_ps,f_dps)
-  USE icosa
-  USE mpipara
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_dps(:)
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: dps(:)
-    REAL(rstd) :: mass_tot,dmass_tot
-    INTEGER :: ind,i,j,ij
-    
-    mass_tot=0
-    dmass_tot=0
-    
-    CALL transfert_request(f_dps,req_i1)
-    CALL transfert_request(f_ps,req_i1)
-
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-
-      ps=f_ps(ind)
-      dps=f_dps(ind)
-
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-            mass_tot=mass_tot+ps(ij)*Ai(ij)/g
-            dmass_tot=dmass_tot+dps(ij)*Ai(ij)/g
-          ENDIF
-        ENDDO
-      ENDDO
-   
-    ENDDO
-    IF (is_mpi_root) PRINT*, "mass_tot ", mass_tot,"      dmass_tot ",dmass_tot        
-
-  END SUBROUTINE check_mass_conservation  
-  
-  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
-       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
-    USE icosa
-    USE vorticity_mod
-    USE theta2theta_rhodz_mod
-    USE pression_mod
-    USE omega_mod
-    USE write_field
-    USE vertical_interp_mod
-    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
-         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
-    
-    REAL(rstd) :: out_pression_lev
-    CHARACTER(LEN=255) :: str_pression
-    CHARACTER(LEN=255) :: physics_type
-    
-    out_pression_level=0
-    CALL getin("out_pression_level",out_pression_level) 
-    WRITE(str_pression,*) INT(out_pression_level/100)
-    str_pression=ADJUSTL(str_pression)
-    
-    CALL writefield("ps",f_ps)
-    CALL writefield("dps",f_dps)
-    CALL writefield("phis",f_phis)
-    CALL vorticity(f_u,f_buf_v)
-    CALL writefield("vort",f_buf_v)
-
-    CALL w_omega(f_ps, f_u, f_buf_i)
-    CALL writefield('omega', f_buf_i)
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! Temperature
-    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
-    
-    CALL getin('physics',physics_type)
-    IF (TRIM(physics_type)=='dcmip') THEN
-      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
-      CALL writefield("T",f_buf1_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ELSE
-      CALL writefield("T",f_buf_i)
-      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-        CALL writefield("T"//TRIM(str_pression),f_buf_s)
-      ENDIF
-    ENDIF
-   
-    ! velocity components
-    CALL un2ulonlat(f_u, f_buf_i3, f_buf1_i, f_buf2_i)
-    CALL writefield("ulon",f_buf1_i)
-    CALL writefield("ulat",f_buf2_i)
-
-    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
-      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
-      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
-      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
-    ENDIF
-    
-    ! geopotential
-    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
-    CALL writefield("p",f_buf_p)
-    CALL writefield("phi",f_buf_i)
-    CALL writefield("theta",f_buf1_i) ! potential temperature
-    CALL writefield("pk",f_buf2_i) ! Exner pressure
-  
-  
-  END SUBROUTINE write_output_fields
-  
-  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
-    USE field_mod
-    USE pression_mod
-    USE exner_mod
-    USE geopotential_mod
-    USE theta2theta_rhodz_mod
-    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
-         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
-    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
-         phi(:,:), phis(:), ps(:), pks(:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       ps = f_ps(ind)
-       p  = f_p(ind)
-       CALL compute_pression(ps,p,0)
-       pk = f_pk(ind)
-       pks = f_pks(ind)
-       CALL compute_exner(ps,p,pks,pk,0)
-       theta_rhodz = f_theta_rhodz(ind)
-       theta = f_theta(ind)
-       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
-       phis = f_phis(ind)
-       phi = f_phi(ind)
-       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
-    END DO
-
-  END SUBROUTINE thetarhodz2geopot
-  
-  SUBROUTINE un2ulonlat(f_u, f_u3d, f_ulon, f_ulat)
-    USE field_mod
-    USE wind_mod
-    TYPE(t_field), POINTER :: f_u(:), & ! IN  : normal velocity components on edges
-         f_u3d(:), f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
-    REAL(rstd),POINTER :: u(:,:), u3d(:,:,:), ulon(:,:), ulat(:,:)
-    INTEGER :: ind
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       u=f_u(ind)
-       u3d=f_u3d(ind)
-       CALL compute_wind_centered(u,u3d)
-       ulon=f_ulon(ind)
-       ulat=f_ulat(ind)
-       CALL compute_wind_centered_lonlat_compound(u3d, ulon, ulat)
-    END DO
-  END SUBROUTINE un2ulonlat
-
-  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_TV(:)
-    TYPE(t_field), POINTER :: f_q(:)
-    TYPE(t_field), POINTER :: f_T(:)
-    
-    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
-    INTEGER :: ind
-    
-    DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-       Tv=f_Tv(ind)
-       q=f_q(ind)
-       T=f_T(ind)
-       T=Tv/(1+0.608*q(:,:,1))
-    END DO
-    
-  END SUBROUTINE Tv2T
-  
-END MODULE caldyn_gcm_opt_mod
Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 150)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 151)
@@ -1,24 +1,24 @@
 MODULE check_conserve_mod
-	USE icosa 
-	IMPLICIT NONE 
-
-	PRIVATE 
-     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
-     TYPE(t_field),POINTER,SAVE :: f_vort(:)
-     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
+  USE icosa 
+  IMPLICIT NONE 
+
+  PRIVATE 
+    TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
+    TYPE(t_field),POINTER,SAVE :: f_vort(:)
+    TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
    
-     PUBLIC compute_conserve 
-     REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
-     REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
-     REAL(rstd),SAVE:: ztot,mtot
+  PUBLIC init_check_conserve, check_conserve 
+    REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
+    REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
+    REAL(rstd),SAVE:: ztot,mtot
 	
      
-Contains 
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE allocate_check_conserve
-    USE icosa
-    IMPLICIT NONE
+CONTAINS 
+
+!---------------------------------------------------------------------
+
+  SUBROUTINE init_check_conserve
+  USE icosa
+  IMPLICIT NONE
     CALL allocate_field(f_pk,field_t,type_real,llm)
     CALL allocate_field(f_p,field_t,type_real,llm+1)
@@ -26,9 +26,73 @@
     CALL allocate_field(f_rhodz,field_t,type_real,llm)
     CALL allocate_field(f_vort,field_z,type_real,llm)
-   END SUBROUTINE allocate_check_conserve
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_mass_conserve(f_ps,f_dps)
+  END SUBROUTINE init_check_conserve
+
+  SUBROUTINE check_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
+  USE icosa 
+  USE pression_mod
+  USE vorticity_mod
+  USE caldyn_gcm_mod
+  USE exner_mod 
+  USE mpipara
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_dps(:)
+    TYPE(t_field),POINTER :: f_ue(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_phis(:)
+    INTEGER::it 
+    REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
+    INTEGER::ind 
+
+    etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
+    mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
+
+    CALL pression(f_ps,f_p)
+
+    DO ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+        p=f_p(ind) 
+       rhodz=f_rhodz(ind)
+      CALL compute_rhodz(p,rhodz) 
+    END DO 
+
+    CALL vorticity(f_ue,f_vort)
+    CALL check_mass_conserve(f_ps,f_dps)
+    CALL check_PV 
+    CALL exner(f_ps,f_p,f_pks,f_pk)
+    CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
+       
+     IF ( it == 0  ) Then 
+       ztot0 = ztot
+       mtot0 = mtot
+       etot0 = etot 
+       ang0  = ang 
+       stot0 = stot 
+     END IF 
+
+     rmsv=SQRT(rmsv/mtot) 
+     ztot=ztot/ztot0 ; mtot=mtot/mtot0 
+     etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
+     rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
+
+
+     IF (is_mpi_root) THEN 
+       OPEN(134,file="checkconsicosa.txt",position='append')
+       WRITE(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
+       WRITE(134,*)"=================================================="
+       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+
+4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
+     ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
+      ,f10.6,e13.6,5f10.3/)     
+       close(134) 
+     END IF 
+  END SUBROUTINE check_conserve
+  
+!---------------------------------------------------------------------
+
+  SUBROUTINE check_mass_conserve(f_ps,f_dps)
   USE icosa
   IMPLICIT NONE
@@ -38,5 +102,5 @@
     INTEGER :: ind,i,j,ij  
 
-   DO ind=1,ndomain
+    DO ind=1,ndomain
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -49,5 +113,5 @@
             mtot=mtot+ps(ij)*Ai(ij)
 	       rmsdpdt=rmsdpdt+dps(ij)*dps(ij)
-	     ENDIF
+	  ENDIF
         ENDDO
       ENDDO
@@ -57,9 +121,9 @@
 !---------------------------------------------------------------------
 
-   SUBROUTINE check_EN(f_ue,f_theta_rhodz,f_phis) 
-    USE icosa
-    USE pression_mod
-    USE vorticity_mod
-    IMPLICIT NONE
+  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis) 
+  USE icosa
+  USE pression_mod
+  USE vorticity_mod
+  IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ue(:)
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
@@ -74,22 +138,23 @@
 
 
-      Do ind=1,ndomain 
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-	 ue=f_ue(ind) 
-         p=f_p(ind)
-         rhodz=f_rhodz(ind) 
-         theta_rhodz=f_theta_rhodz(ind) 
-	    pk=f_pk(ind) 
-	    phis=f_phis(ind) 
-        CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
-	 END DO 
-  END SUBROUTINE CHECK_EN 
-
-   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
-	 USE icosa
-      USE disvert_mod
-      USE wind_mod
-     IMPLICIT NONE
+    DO ind=1,ndomain 
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ue=f_ue(ind) 
+      p=f_p(ind)
+      rhodz=f_rhodz(ind) 
+      theta_rhodz=f_theta_rhodz(ind) 
+      pk=f_pk(ind) 
+      phis=f_phis(ind) 
+      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
+    END DO 
+ 
+  END SUBROUTINE check_en 
+
+  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
+  USE icosa
+  USE disvert_mod
+  USE wind_mod
+  IMPLICIT NONE
     INTEGER,INTENT(IN)::ind 
     REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
@@ -110,89 +175,88 @@
 
   
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-	         masse(ij,l) = rhodz(ij,l)*Ai(ij) 
-              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-             IF (domain(ind)%own(i,j)) THEN
-	         stot = stot + Ai(ij)*theta_rhodz(ij,l)
-             END IF 
-          END DO 
-         END DO
-	   END DO 
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-	        KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
-                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
-                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
-                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
-                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
-                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
-            END IF 
-	     END DO
-         END DO 
-        END DO  
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-	      etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
-           END IF 
-          END DO 
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+	  ij=(j-1)*iim+i
+	  masse(ij,l) = rhodz(ij,l)*Ai(ij) 
+          theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+          IF (domain(ind)%own(i,j)) THEN
+	    stot = stot + Ai(ij)*theta_rhodz(ij,l)
+          END IF 
         END DO 
-       END DO 
+      END DO
+    END DO 
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+	    KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
+                                      le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
+                                      le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
+                                      le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
+                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
+                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+            rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
+          END IF 
+        END DO
+      END DO 
+    END DO  
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+	  ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+	    etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
+          END IF 
+        END DO 
+      END DO 
+    END DO 
 
 !------------------------------ ANGULAR VELOCITY 
-	CALL compute_wind_centered(u,ucenter)
-     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
-	 rad = radius 
-      radomeg = rad*omega
-
-      DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-             CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
-                ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
-           END IF 
-         END DO
-       END DO
-      END DO    
+    CALL compute_wind_centered(u,ucenter)
+    CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
+    rad = radius 
+    radomeg = rad*omega
+
+    DO l = 1, llm
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+	  ij=(j-1)*iim+i
+          IF (domain(ind)%own(i,j)) THEN
+            CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
+            ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
+          END IF 
+        END DO
+      END DO
+    END DO    
           
-	END SUBROUTINE compute_en
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_PV 
-    USE icosa
-    IMPLICIT NONE
-
-    REAL(rstd), POINTER :: vort(:,:)
-    REAL(rstd), POINTER :: rhodz(:,:) 
-    INTEGER :: ind
-
-    Do ind=1,ndomain 
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      vort=f_vort(ind)
-      rhodz=f_rhodz(ind) 
-      CALL compute_PV(vort,rhodz) 
-    ENDDO
+  END SUBROUTINE compute_en
+
+!---------------------------------------------------------------------
+
+  SUBROUTINE check_PV 
+  USE icosa
+  IMPLICIT NONE
+     REAL(rstd), POINTER :: vort(:,:)
+     REAL(rstd), POINTER :: rhodz(:,:) 
+     INTEGER :: ind
+
+     DO ind=1,ndomain 
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       vort=f_vort(ind)
+       rhodz=f_rhodz(ind) 
+       CALL compute_PV(vort,rhodz) 
+     ENDDO
  
-   END SUBROUTINE check_PV 
-
-	SUBROUTINE compute_PV(vort,rhodz) 
-	  USE icosa
-      USE disvert_mod
-     IMPLICIT NONE 
+  END SUBROUTINE check_PV 
+
+  SUBROUTINE compute_PV(vort,rhodz) 
+  USE icosa
+  USE disvert_mod
+  IMPLICIT NONE 
     REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
@@ -201,32 +265,34 @@
     INTEGER :: i,j,ij,l,ij2
 
-	hv1 = 0.0 ; hv2 = 0.0 
+    hv1 = 0.0 ; hv2 = 0.0 
    
-     DO l = 1,llm
+    DO l = 1,llm
       DO j=jj_begin+1,jj_end-1
-       DO i=ii_begin+1,ii_end-1
-           ij=(j-1)*iim+i
-            hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
-              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-      
-            qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
-            
-            hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
-              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-                       
-            qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
-
-            ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
-	  ENDDO
-         ENDDO
+        DO i=ii_begin+1,ii_end-1
+          ij=(j-1)*iim+i
+          hv1 =  Riv2(ij,vup)        * rhodz(ij,l)           &
+            + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
+            + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
+    
+          qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
+          
+          hv2  =  Riv2(ij,vdown)     * rhodz(ij,l)          &
+            + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
+            + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
+                     
+          qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
+
+          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
+
         ENDDO
-    END SUBROUTINE compute_PV
-!---------------------------------------------------------------------
-
-    SUBROUTINE compute_rhodz(p,rhodz) 
-	USE icosa 
-	IMPLICIT NONE 
+      ENDDO
+    ENDDO
+
+  END SUBROUTINE compute_PV
+!---------------------------------------------------------------------
+
+   SUBROUTINE compute_rhodz(p,rhodz) 
+   USE icosa 
+   IMPLICIT NONE 
      REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
      REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
@@ -244,71 +310,7 @@
      ENDDO
 
-	END SUBROUTINE compute_rhodz
+   END SUBROUTINE compute_rhodz
 !====================================================================
 
-	SUBROUTINE compute_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
-	 USE icosa 
-      USE pression_mod
-      USE vorticity_mod
-      USE caldyn_gcm_mod
-      USE exner_mod 
-      USE mpipara
-     TYPE(t_field),POINTER :: f_ps(:)
-     TYPE(t_field),POINTER :: f_dps(:)
-     TYPE(t_field),POINTER :: f_ue(:)
-     TYPE(t_field),POINTER :: f_theta_rhodz(:)
-     TYPE(t_field),POINTER :: f_phis(:)
-     INTEGER::it 
-     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
-     INTEGER::ind 
-
-     etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
-     mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
-
-     CALL allocate_check_conserve
-
-     CALL pression(f_ps,f_p)
-
-     Do ind=1,ndomain 
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-         p=f_p(ind) 
-        rhodz=f_rhodz(ind)
-       CALL compute_rhodz(p,rhodz) 
-      END DO 
-
-     CALL vorticity(f_ue,f_vort)
-     CALL check_mass_conserve(f_ps,f_dps)
-     CALL check_PV 
-     CALL exner(f_ps,f_p,f_pks,f_pk)
-     CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
-    	
-	   IF ( it == 0  ) Then 
-	   ztot0 = ztot
-           mtot0 = mtot
-           etot0 = etot 
-           ang0  = ang 
-           stot0 = stot 
-        END IF 
-
-      rmsv=SQRT(rmsv/mtot) 
-      ztot=ztot/ztot0 ; mtot=mtot/mtot0 
-      etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
-      rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
-
-
-	  IF (is_mpi_root) THEN 
-       OPEN(134,file="checkconsicosa.txt",position='append')
-       write(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-       write(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
-       write(134,*)"=================================================="
-       write(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-
-4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
-      ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
-       ,f10.6,e13.6,5f10.3/)     
-       close(134) 
-       END IF 
-    END SUBROUTINE compute_conserve
 
  END MODULE check_conserve_mod
Index: codes/icosagcm/trunk/src/check_conserve.f90~
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90~	(revision 150)
+++ 	(revision )
@@ -1,315 +1,0 @@
-MODULE check_conserve_mod
-	USE icosa 
-	IMPLICIT NONE 
-
-	PRIVATE 
-     TYPE(t_field),POINTER,SAVE :: f_pk(:), f_pks(:), f_p(:)
-     TYPE(t_field),POINTER,SAVE :: f_vort(:)
-     TYPE(t_field),POINTER,SAVE :: f_rhodz(:) 
-   
-     PUBLIC compute_conserve 
-     REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
-     REAL(rstd),SAVE::rmsdpdt,etot,ang,stot,rmsv
-     REAL(rstd),SAVE:: ztot,mtot
-	
-     
-Contains 
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE allocate_check_conserve
-    USE icosa
-    IMPLICIT NONE
-    CALL allocate_field(f_pk,field_t,type_real,llm)
-    CALL allocate_field(f_p,field_t,type_real,llm+1)
-    CALL allocate_field(f_pks,field_t,type_real)
-    CALL allocate_field(f_rhodz,field_t,type_real,llm)
-    CALL allocate_field(f_vort,field_z,type_real,llm)
-   END SUBROUTINE allocate_check_conserve
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_mass_conserve(f_ps,f_dps)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_dps(:)
-    REAL(rstd),POINTER :: ps(:),dps(:) 
-    INTEGER :: ind,i,j,ij  
-
-   DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      dps=f_dps(ind) 
-      ps=f_ps(ind)
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          IF (domain(ind)%own(i,j)) THEN
-            mtot=mtot+ps(ij)*Ai(ij)
-	       rmsdpdt=rmsdpdt+dps(ij)*dps(ij)
-	     ENDIF
-        ENDDO
-      ENDDO
-    ENDDO
-  END SUBROUTINE check_mass_conserve
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_EN(f_ue,f_theta_rhodz,f_phis) 
-    USE icosa
-    USE pression_mod
-    USE vorticity_mod
-    IMPLICIT NONE
-    TYPE(t_field), POINTER :: f_ue(:)
-    TYPE(t_field), POINTER :: f_theta_rhodz(:)
-    TYPE(t_field), POINTER :: f_phis(:)
-    REAL(rstd), POINTER :: ue(:,:)
-    REAL(rstd), POINTER :: p(:,:) 
-    REAL(rstd), POINTER :: theta_rhodz(:,:) 
-    REAL(rstd), POINTER :: pk(:,:) 
-    REAL(rstd), POINTER :: phis(:) 
-    REAL(rstd), POINTER :: rhodz(:,:) 
-    INTEGER :: ind
-
-
-      Do ind=1,ndomain 
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-	    ue=f_ue(ind) 
-         p=f_p(ind)
-         rhodz=f_rhodz(ind) 
-         theta_rhodz=f_theta_rhodz(ind) 
-	    pk=f_pk(ind) 
-	    phis=f_phis(ind) 
-        CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
-	 END DO 
-  END SUBROUTINE CHECK_EN 
-
-   SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
-	 USE icosa
-      USE disvert_mod
-      USE wind_mod
-     IMPLICIT NONE
-    INTEGER,INTENT(IN)::ind 
-    REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
-    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: phis(iim*jjm)
-    REAL(rstd):: theta(iim*jjm,llm)
-    REAL(rstd)::KE(iim*jjm,llm) 
-    REAL(rstd) :: ucenter(iim*jjm,3,llm)
-    REAL(rstd) :: ulon(iim*jjm,llm)
-    REAL(rstd) :: ulat(iim*jjm,llm)
-
-    REAL(rstd)::masse(iim*jjm,llm)  
-    REAL(rstd)::rad,radomeg,lat,lon
-    INTEGER :: i,j,ij,l,ij2
-
-  
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-	         masse(ij,l) = rhodz(ij,l)*Ai(ij) 
-              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-             IF (domain(ind)%own(i,j)) THEN
-	         stot = stot + Ai(ij)*theta_rhodz(ij,l)
-             END IF 
-          END DO 
-         END DO
-	   END DO 
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-	        KE(ij,l)=   1/(4*Ai(ij))*(le(ij+u_right)*de(ij+u_right)*u(ij+u_right,l)**2 +   &
-                                     le(ij+u_rup)*de(ij+u_rup)*u(ij+u_rup,l)**2 +           &
-                                     le(ij+u_lup)*de(ij+u_lup)*u(ij+u_lup,l)**2 +           &
-                                     le(ij+u_left)*de(ij+u_left)*u(ij+u_left,l)**2 +        &
-                                     le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +     &
-                                     le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-             rmsv = rmsv + 2*masse(ij,l)*KE(ij,l) 
-            END IF 
-	     END DO
-         END DO 
-        END DO  
-
-       DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-	      etot = etot + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
-           END IF 
-          END DO 
-        END DO 
-       END DO 
-
-!------------------------------ ANGULAR VELOCITY 
-	CALL compute_wind_centered(u,ucenter)
-     CALL compute_wind_centered_lonlat_compound(ucenter,ulon,ulat)
-	 rad = radius 
-      radomeg = rad*omega
-
-      DO l = 1, llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-	        ij=(j-1)*iim+i
-           IF (domain(ind)%own(i,j)) THEN
-             CALL xyz2lonlat(xyz_i(ij,:),lon,lat) 
-                ang=ang+rad*cos(lat)*masse(ij,l)*(ulon(ij,l)+ radomeg*cos(lat))
-           END IF 
-         END DO
-       END DO
-      END DO    
-          
-	END SUBROUTINE compute_en
-
-!---------------------------------------------------------------------
-
-   SUBROUTINE check_PV 
-    USE icosa
-    IMPLICIT NONE
-
-    REAL(rstd), POINTER :: vort(:,:)
-    REAL(rstd), POINTER :: rhodz(:,:) 
-    INTEGER :: ind
-
-    Do ind=1,ndomain 
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      vort=f_vort(ind)
-      rhodz=f_rhodz(ind) 
-      CALL compute_PV(vort,rhodz) 
-    ENDDO
- 
-   END SUBROUTINE check_PV 
-
-	SUBROUTINE compute_PV(vort,rhodz) 
-	  USE icosa
-      USE disvert_mod
-     IMPLICIT NONE 
-    REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
-    REAL(rstd)::qv1,qv2 
-    REAL(rstd)::hv1,hv2  
-    INTEGER :: i,j,ij,l,ij2
-
-	hv1 = 0.0 ; hv2 = 0.0 
-     
-     DO l = 1,llm
-      DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-           ij=(j-1)*iim+i
-            hv1 =  Riv2(ij,vup)       * rhodz(ij,l)           &
-              + Riv2(ij+t_rup,vldown) * rhodz(ij+t_rup,l)     &
-              + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
-      
-            qv1 = ( vort(ij+z_up,l)+fv(ij+z_up) )/hv1
-            
-            hv2 =  Riv2(ij,vdown)     * rhodz(ij,l)          &
-              + Riv2(ij+t_ldown,vrup) * rhodz(ij+t_ldown,l)  &
-              + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
-                       
-            qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
-
-            ztot= ztot + (qv1*qv1*hv1 + qv2*qv2*hv2)
-	    ENDDO
-       ENDDO
-      ENDDO    
-    END SUBROUTINE compute_PV
-!---------------------------------------------------------------------
-
-    SUBROUTINE compute_rhodz(p,rhodz) 
-	USE icosa 
-	IMPLICIT NONE 
-     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
-     REAL(rstd),INTENT(OUT):: rhodz(iim*jjm,llm)
-     REAL(rstd) ::mass(iim*jjm,llm)
-     INTEGER :: i,j,ij,l 
-	
-     DO l = 1, llm
-      DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-         ij=(j-1)*iim+i
-         mass(ij,l) = ( p(ij,l) - p(ij,l+1) ) * Ai(ij)/g
-         rhodz(ij,l) = mass(ij,l) / Ai(ij)
-       ENDDO
-      ENDDO
-     ENDDO
-
-	END SUBROUTINE compute_rhodz
-
-!=============================================================================
-
-	SUBROUTINE compute_conserve(f_ps,f_dps,f_ue,f_theta_rhodz,f_phis,it) 
-	 USE icosa 
-      USE pression_mod
-      USE vorticity_mod
-      USE caldyn_gcm_mod
-      USE exner_mod 
-      USE mpipara
-     TYPE(t_field),POINTER :: f_ps(:)
-     TYPE(t_field),POINTER :: f_dps(:)
-     TYPE(t_field),POINTER :: f_ue(:)
-     TYPE(t_field),POINTER :: f_theta_rhodz(:)
-     TYPE(t_field),POINTER :: f_phis(:)
-     INTEGER::it 
-     REAL(rstd),POINTER :: p(:,:),rhodz(:,:) 
-     INTEGER::ind 
-
-     etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
-     mtot=0.0 ; rmsdpdt=0.0 ; ztot = 0.0 
-
-     CALL allocate_check_conserve
-
-     CALL pression(f_ps,f_p)
-
-     Do ind=1,ndomain 
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-         p=f_p(ind) 
-        rhodz=f_rhodz(ind)
-       CALL compute_rhodz(p,rhodz) 
-      END DO 
-
-     CALL vorticity(f_ue,f_vort)
-     CALL check_mass_conserve(f_ps,f_dps)
-     CALL check_PV 
-     CALL exner(f_ps,f_p,f_pks,f_pk)
-     CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
-    	
-	   IF ( it == 0  ) Then 
-		 ztot0 = ztot
-           mtot0 = mtot
-           etot0 = etot 
-           ang0  = ang 
-           stot0 = stot 
-        END IF 
-
-      rmsv=SQRT(rmsv/mtot) 
-	 ztot=ztot/ztot0 ; mtot=mtot/mtot0 
-      etot=etot/etot0 ; ang=ang/ang0 ; stot=stot/stot0 
-	 rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
-
-
-	  IF (is_mpi_root) THEN 
-       OPEN(134,file="checkconsicosa.txt",position='append')
-       write(134,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-       write(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang 
-       write(134,*)"=================================================="
-       write(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
-
-4000   FORMAT(10x,'masse',4x,'rmsdpdt',7x,'energie',2x,'enstrophie'  &
-      ,2x,'entropie',3x,'rmsv',4x,'mt.ang',/,'GLOB  '                &
-       ,f10.6,e13.6,5f10.3/)     
-       close(134) 
-       END IF 
-    END SUBROUTINE compute_conserve
-
- END MODULE check_conserve_mod
Index: codes/icosagcm/trunk/src/diffff
===================================================================
--- codes/icosagcm/trunk/src/diffff	(revision 150)
+++ 	(revision )
@@ -1,493 +1,0 @@
-Index: time.f90
-===================================================================
---- time.f90	(revision 142)
-+++ time.f90	(working copy)
-@@ -10,11 +10,23 @@
-   REAL(rstd),SAVE :: write_period
-   INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax
-   
-+  INTEGER,SAVE :: day_step,ndays
-+  REAL(rstd),SAVE :: jD_ref,jH_ref
-+  INTEGER,SAVE :: day_ini,day_end,annee_ref,day_ref
-+  REAL(rstd),SAVE::start_time
-+  CHARACTER(LEN=255) :: time_style
-+  INTEGER,SAVE:: an, mois, jour
-+  REAL(rstd),SAVE:: heure
-+  CHARACTER (LEN=10):: calend
-+
-   PUBLIC create_time_counter_header, update_time_counter, close_time_counter, init_time,  &
--         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax
-+         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax, & 
-+day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, &
-+            calend,time_style
- 
- 
- 
-+
- CONTAINS
-   
-   SUBROUTINE init_time
-@@ -23,22 +35,18 @@
-   USE mpipara
-   IMPLICIT NONE
-   REAL(rstd) :: run_length
--  
-+
-+
-+     time_style='dcmip'
-+   CALL getin('time_style',time_style)
-+
-+   IF (TRIM(time_style)=='dcmip')  Then
-     dt=90.
-     CALL getin('dt',dt)
- 
-     itaumax=100
-     CALL getin('itaumax',itaumax)
- 
--    itau_adv=1
--    CALL getin('itau_adv',itau_adv)
--
--    itau_dissip=1
--    CALL getin('itau_dissip',itau_dissip)
--
--    itau_physics=1
--    CALL getin('itau_physics',itau_physics)
--
-     run_length=dt*itaumax
-     CALL getin('run_length',run_length)
-     itaumax=run_length/dt
-@@ -53,6 +61,17 @@
-     write_period=write_period/scale_factor
-     itau_out=FLOOR(.5+write_period/dt)
-     IF (is_mpi_root) PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
-+    ENDIF 
-+
-+    itau_adv=1
-+    CALL getin('itau_adv',itau_adv)
-+
-+    itau_dissip=1
-+    CALL getin('itau_dissip',itau_dissip)
-+
-+    itau_physics=1
-+    CALL getin('itau_physics',itau_physics)
-+
-     
-     CALL create_time_counter_header
-     
-Index: caldyn_gcm.f90
-===================================================================
---- caldyn_gcm.f90	(revision 142)
-+++ caldyn_gcm.f90	(working copy)
-@@ -10,7 +10,7 @@
-   TYPE(t_field),POINTER :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
-   TYPE(t_field),POINTER :: f_buf_v(:), f_buf_s(:), f_buf_p(:)
- 
--  PUBLIC init_caldyn, caldyn, write_output_fields
-+  PUBLIC init_caldyn, caldyn, write_output_fields,un2ulonlat
- 
-   INTEGER :: caldyn_hydrostat, caldyn_conserv
- 
-@@ -708,17 +708,17 @@
-     str_pression=ADJUSTL(str_pression)
-     
-     CALL writefield("ps",f_ps)
--    CALL writefield("dps",f_dps)
--    CALL writefield("phis",f_phis)
--    CALL vorticity(f_u,f_buf_v)
--    CALL writefield("vort",f_buf_v)
-+!    CALL writefield("dps",f_dps)
-+!    CALL writefield("phis",f_phis)
-+!    CALL vorticity(f_u,f_buf_v)
-+!    CALL writefield("vort",f_buf_v)
- 
--    CALL w_omega(f_ps, f_u, f_buf_i)
--    CALL writefield('omega', f_buf_i)
--    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
--      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
--      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
--    ENDIF
-+!    CALL w_omega(f_ps, f_u, f_buf_i)
-+!    CALL writefield('omega', f_buf_i)
-+!    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
-+!      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
-+!      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
-+!    ENDIF
-     
-     ! Temperature
-     CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ;
-@@ -753,10 +753,10 @@
-     
-     ! geopotential
-     CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
--    CALL writefield("p",f_buf_p)
--    CALL writefield("phi",f_buf_i)
--    CALL writefield("theta",f_buf1_i) ! potential temperature
--    CALL writefield("pk",f_buf2_i) ! Exner pressure
-+!    CALL writefield("p",f_buf_p)
-+!    CALL writefield("phi",f_buf_i)
-+     CALL writefield("theta",f_buf1_i) ! potential temperature
-+!    CALL writefield("pk",f_buf2_i) ! Exner pressure
-   
-   
-   END SUBROUTINE write_output_fields
-Index: advect_tracer.f90
-===================================================================
---- advect_tracer.f90	(revision 142)
-+++ advect_tracer.f90	(working copy)
-@@ -57,7 +57,7 @@
-     CALL transfert_request(f_q,req_i1)
-     CALL transfert_request(f_rhodz,req_i1)
- 	
--    IF (is_mpi_root) PRINT *, 'Advection scheme'
-+!    IF (is_mpi_root) PRINT *, 'Advection scheme'
- 
- !    DO ind=1,ndomain
- !       CALL swap_dimensions(ind)
-Index: dissip_gcm.f90
-===================================================================
---- dissip_gcm.f90	(revision 142)
-+++ dissip_gcm.f90	(working copy)
-@@ -93,7 +93,8 @@
-       rayleigh_friction_type=2
-       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
-    CASE DEFAULT
--      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
-+      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
-+        ' in dissip_gcm.f90/init_dissip'
-       STOP
-    END SELECT
- 
-Index: physics.f90
-===================================================================
---- physics.f90	(revision 142)
-+++ physics.f90	(working copy)
-@@ -7,7 +7,8 @@
- 
-   SUBROUTINE init_physics
-   USE icosa
--  USE physics_dcmip_mod, init_physics_dcmip=>init_physics
-+  USE physics_dcmip_mod,init_physics_dcmip=>init_physics
-+  USE physics_dry_mod
-   IMPLICIT NONE
-     
-     CALL getin("physics",physics_type)
-@@ -17,36 +18,57 @@
-     
-       CASE ('dcmip')
-         CALL init_physics_dcmip
-+
-+      CASE ('lmd')
-+	CALL init_physics_dry
-       
-       CASE DEFAULT
--         PRINT*, 'Bad selector for variable physics <',physics_type, &
-+         PRINT*, 'Bad selector for variable physics init <',physics_type, &
-               '> options are <none>, <dcmip>,'
--         STOP
-+
-     END SELECT
-     
-   END SUBROUTINE init_physics
- 
--  SUBROUTINE physics(it,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-+  SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-   USE icosa
-+  USE physics_dry_mod
-   USE physics_dcmip_mod, physics_dcmip=>physics
-+  USE etat0_mod
-+  USE etat0_heldsz_mod
-   IMPLICIT NONE
-     INTEGER, INTENT(IN)   :: it
-+    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
-     TYPE(t_field),POINTER :: f_phis(:)
-     TYPE(t_field),POINTER :: f_ps(:)
-     TYPE(t_field),POINTER :: f_theta_rhodz(:)
-     TYPE(t_field),POINTER :: f_ue(:)
-     TYPE(t_field),POINTER :: f_q(:)
-+    LOGICAL:: firstcall,lastcall
-     
-     SELECT CASE(TRIM(physics_type))
-       CASE ('none')
-+
-+	SELECT CASE(TRIM(etat0_type))
-+	CASE('heldsz') 
-+	CALL transfert_request(f_ps,req_i1)
-+	CALL transfert_request(f_theta_rhodz,req_i1)
-+	CALL transfert_request(f_ue,req_e1)
-+	CALL held_saurez(f_ps,f_theta_rhodz,f_ue) 
-+	CASE DEFAULT
-+	PRINT*,"NO PHYSICAL PACAKAGE USED" 
-+	END SELECT
-     
-       CASE ('dcmip')
-         CALL physics_dcmip(it, f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-+
-+      CASE ('dry')
-+        CALL physics_dry(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
-       
-       CASE DEFAULT
-          PRINT*, 'Bad selector for variable physics <',physics_type, &
-               '> options are <none>, <dcmip>,'
--         STOP
-+  STOP
-     END SELECT
-     
-   END SUBROUTINE physics
-Index: etat0.f90
-===================================================================
---- etat0.f90	(revision 142)
-+++ etat0.f90	(working copy)
-@@ -1,4 +1,5 @@
- MODULE etat0_mod
-+    CHARACTER(len=255),SAVE :: etat0_type
- 
- CONTAINS
-   
-@@ -11,6 +12,10 @@
-     USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
-     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
-     USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0  
-+    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
-+    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 
-+    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0 
-+
-     IMPLICIT NONE
-     TYPE(t_field),POINTER :: f_ps(:)
-     TYPE(t_field),POINTER :: f_phis(:)
-@@ -18,7 +23,6 @@
-     TYPE(t_field),POINTER :: f_u(:)
-     TYPE(t_field),POINTER :: f_q(:)
-     
--    CHARACTER(len=255) :: etat0_type
-     etat0_type='jablonowsky06'
-     CALL getin("etat0",etat0_type)
-     
-@@ -27,6 +31,9 @@
-        CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('academic')
-        CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-+    CASE ('heldsz')
-+	print*,"heldsz test case"
-+       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip1')
-        CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
-@@ -37,6 +44,12 @@
-        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-      CASE ('dcmip5')
-        CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-+     CASE ('readnf_start') 
-+	  print*,"readnf_start used"	
-+       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
-+	CASE ('readnf_hz') 
-+	  print*,"readnf_hz used"
-+       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
-    CASE DEFAULT
-        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
-             '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
-Index: timeloop_gcm.f90
-===================================================================
---- timeloop_gcm.f90	(revision 142)
-+++ timeloop_gcm.f90	(working copy)
-@@ -17,7 +17,11 @@
-   USE advect_tracer_mod
-   USE physics_mod
-   USE mpipara
--  
-+!-------------------------------
-+  USE IOIPSL
-+  USE maxicosa
-+  USE check_conserve_mod
-+!------------------------------
-   IMPLICIT NONE
-   TYPE(t_field),POINTER :: f_phis(:)
- !  TYPE(t_field),POINTER :: f_theta(:)
-@@ -46,6 +50,10 @@
-   INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
-   CHARACTER(len=255) :: scheme_name
-   LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
-+  CHARACTER(len=7) :: first 
-+  REAL(rstd),SAVE :: jD_cur, jH_cur
-+  REAL(rstd),SAVE :: start_time
-+
- !  INTEGER :: itaumax
- !  REAL(rstd) ::write_period
- !  INTEGER    :: itau_out
-@@ -84,6 +92,58 @@
-   CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
-   CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
- 
-+!----------------------------------------------------
-+  IF (TRIM(time_style)=='lmd')  Then
-+
-+   day_step=180
-+   CALL getin('day_step',day_step)
-+
-+   ndays=1
-+   CALL getin('ndays',ndays)
-+
-+   dt = daysec/REAL(day_step)
-+   itaumax = ndays*day_step
-+
-+   calend = 'earth_360d'
-+   CALL getin('calend', calend)
-+
-+   day_ini = 0
-+   CALL getin('day_ini',day_ini)
-+
-+   day_end = 0
-+   CALL getin('day_end',day_end)
-+
-+   annee_ref = 1998
-+   CALL getin('annee_ref',annee_ref)
-+
-+   start_time = 0
-+   CALL getin('start_time',start_time) 
-+
-+   write_period=0
-+   CALL getin('write_period',write_period)
-+      
-+   write_period=write_period/scale_factor
-+   itau_out=FLOOR(write_period/dt)
-+   
-+   PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out 
-+
-+  mois = 1 ; heure = 0.
-+  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
-+  jH_ref = jD_ref - int(jD_ref) 
-+  jD_ref = int(jD_ref) 
-+  
-+  CALL ioconf_startdate(INT(jD_ref),jH_ref) 
-+  write(*,*)'annee_ref, mois, day_ref, heure, jD_ref'
-+  write(*,*)annee_ref, mois, day_ref, heure, jD_ref
-+  write(*,*)"ndays,day_step,itaumax,dt======>"
-+  write(*,*)ndays,day_step,itaumax,dt
-+  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
-+  write(*,*)'jD_ref+jH_ref,an, mois, jour, heure'
-+  write(*,*)jD_ref+jH_ref,an, mois, jour, heure
-+  day_end = day_ini + ndays 
-+ END IF 
-+!----------------------------------------------------
-+
-   scheme_name='runge_kutta'
-   CALL getin('scheme',scheme_name)
- 
-@@ -122,7 +182,6 @@
- !  CALL allocate_field(f_dum2,field_u,type_real,llm)
- !  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
- !  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
--
- !  CALL allocate_field(f_theta,field_t,type_real,llm)
- !  CALL allocate_field(f_dtheta,field_t,type_real,llm)
- 
-@@ -132,7 +191,11 @@
-   CALL init_advect_tracer
-   CALL init_physics
-   
--  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-+!============================ SET 0 to all field 
-+  CALL initial0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-+! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-+  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-+!================================================
-   CALL writefield("phis",f_phis,once=.TRUE.)
-   CALL transfert_request(f_q,req_i1) 
- 
-@@ -153,17 +216,18 @@
-      rhodz=f_rhodz(ind); ps=f_ps(ind)
-      CALL compute_rhodz(.FALSE., ps, rhodz)   
-   END DO
-+
-   
-   DO it=0,itaumax
- 
--    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
-     IF (mod(it,itau_out)==0 ) THEN
--      CALL writefield("q",f_q)
-+!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
-       CALL update_time_counter(dt*it)
-+      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
-     ENDIF
--    
--    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
- 
-+      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
-+
-     DO stage=1,nb_stage
-        CALL caldyn((stage==1) .AND. (MOD(it,itau_out)==0), &
-             f_phis,f_ps,f_theta_rhodz,f_u, f_q, &
-@@ -177,11 +241,11 @@
-           CALL  leapfrog_matsuno_scheme(stage)
-           
-           !      CASE ('leapfrog')
--          !        CALL leapfrog_scheme
-+          !      CALL leapfrog_scheme
-           !
-           !      CASE ('adam_bashforth')
--          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
--          !        CALL adam_bashforth_scheme
-+          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-+          !      CALL adam_bashforth_scheme
-        CASE DEFAULT
-           STOP
-        END SELECT
-@@ -192,7 +256,7 @@
- 
-     IF(MOD(it+1,itau_adv)==0) THEN
-        CALL transfert_request(f_wfluxt,req_i1) ! FIXME
--!       CALL transfert_request(f_hfluxt,req_e1) ! FIXME
-+!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
- 
-        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
-        fluxt_zero=.TRUE.
-@@ -207,11 +271,17 @@
-           CALL swap_geometry(ind)
-           rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
-           wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
--          CALL compute_rhodz(.FALSE., ps, rhodz)   
-+!!!!!          CALL compute_rhodz(.FALSE., ps, rhodz)   
-        END DO
--
-     END IF
--
-+!----------------------------------------------------
-+  jD_cur = jD_ref + day_ini - day_ref + it/day_step
-+  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
-+  jD_cur = jD_cur + int(jH_cur)
-+  jH_cur = jH_cur - int(jH_cur)
-+!	print*,"Just b4 phys"
-+    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
-+!----------------------------------------------------
- !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
-     ENDDO
-  
-@@ -399,7 +469,7 @@
-        DO j=jj_begin-dd,jj_end+dd
-           DO i=ii_begin-dd,ii_end+dd
-              ij=(j-1)*iim+i
--             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
-+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
-              IF(comp) THEN
-                 rhodz(ij,l) = m
-              ELSE
-@@ -414,7 +484,7 @@
-           PRINT *, 'Discrepancy between ps and rhodz detected', err
-           STOP
-        ELSE
--!          PRINT *, 'No discrepancy between ps and rhodz detected'
-+          PRINT *, 'No discrepancy between ps and rhodz detected'
-        END IF
-     END IF
- 
-Index: icosa_gcm.f90
-===================================================================
---- icosa_gcm.f90	(revision 142)
-+++ icosa_gcm.f90	(working copy)
-@@ -24,6 +24,7 @@
-   CALL compute_domain
-   CALL init_transfert
-   CALL init_writefield
-+
- !  CALL allocate_field(sum_ne,field_T,type_real)
- !  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
- !  
Index: codes/icosagcm/trunk/src/dimensions.f90
===================================================================
--- codes/icosagcm/trunk/src/dimensions.f90	(revision 150)
+++ codes/icosagcm/trunk/src/dimensions.f90	(revision 151)
@@ -44,4 +44,6 @@
     TYPE(t_domain),POINTER :: d
     
+    
+!$OMP MASTER
     d=>domain(ind)
      
@@ -79,5 +81,7 @@
      u_pos(:)=d%u_pos(:)
      z_pos(:)=d%z_pos(:)
-     
+
+!$OMP END MASTER     
+!$OMP BARRIER
    END SUBROUTINE swap_dimensions
    
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 150)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 151)
@@ -10,5 +10,5 @@
   TYPE(t_field),POINTER,SAVE :: f_dtheta_diss(:)
   TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz_diss(:)
-  
+  TYPE(t_message) :: req_due, req_dtheta 
  
   INTEGER,SAVE :: nitergdiv=1
@@ -29,5 +29,5 @@
 !  INTEGER,SAVE :: itau_dissip
   REAL,SAVE    :: dtdissip
-  
+ 
   PUBLIC init_dissip, dissip
 CONTAINS
@@ -57,4 +57,5 @@
   USE mpi_mod
   USE mpipara
+  USE transfert_mod
   USE time_mod
   
@@ -96,6 +97,5 @@
       IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
    CASE DEFAULT
-      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), & 
-        ' in dissip_gcm.f90/init_dissip'
+      IF (is_mpi_root) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
       STOP
    END SELECT
@@ -116,4 +116,7 @@
     CALL allocate_field(f_theta,field_t,type_real)
     CALL allocate_field(f_dtheta,field_t,type_real)
+    
+    CALL init_message(f_due_diss1,req_e1_vect,req_due)
+    CALL init_message(f_dtheta_diss,req_i1,req_dtheta)
 
     tau_graddiv(:)=5000
@@ -135,44 +138,4 @@
 
     CALL getin("niterdivgrad",niterdivgrad)
-  
-!   CALL create_request(field_u,req_dissip)
-
-!   DO ind=1,ndomain
-!     DO i=ii_begin,ii_end
-!       CALL request_add_point(ind,i,jj_begin-1,req_dissip,rup)
-!       CALL request_add_point(ind,i+1,jj_begin-1,req_dissip,lup)
-!     ENDDO
-
-!     DO j=jj_begin,jj_end
-!       CALL request_add_point(ind,ii_end+1,j,req_dissip,left)
-!       CALL request_add_point(ind,ii_end+1,j-1,req_dissip,lup)
-!     ENDDO    
-!   
-!     DO i=ii_begin,ii_end
-!       CALL request_add_point(ind,i,jj_end+1,req_dissip,ldown)
-!       CALL request_add_point(ind,i-1,jj_end+1,req_dissip,rdown)
-!     ENDDO    
-
-!     DO j=jj_begin,jj_end
-!       CALL request_add_point(ind,ii_begin-1,j,req_dissip,right)
-!       CALL request_add_point(ind,ii_begin-1,j+1,req_dissip,rdown)
-!     ENDDO   
-
-!     DO i=ii_begin+1,ii_end-1
-!       CALL request_add_point(ind,i,jj_begin,req_dissip,right)
-!       CALL request_add_point(ind,i,jj_end,req_dissip,right)
-!     ENDDO
-!   
-!     DO j=jj_begin+1,jj_end-1
-!       CALL request_add_point(ind,ii_begin,j,req_dissip,rup)
-!       CALL request_add_point(ind,ii_end,j,req_dissip,rup)
-!     ENDDO   
-
-!     CALL request_add_point(ind,ii_begin+1,jj_begin,req_dissip,left)
-!     CALL request_add_point(ind,ii_begin,jj_begin+1,req_dissip,ldown)
-!     CALL request_add_point(ind,ii_begin+1,jj_end,req_dissip,left)
-!     CALL request_add_point(ind,ii_end,jj_begin+1,req_dissip,ldown)
-!    
-!   ENDDO
 
 
@@ -214,5 +177,5 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradiv(u,du,1)
+          CALL compute_gradiv(u,du,1,1)
           u=du
         ENDDO
@@ -238,5 +201,5 @@
 
       IF (using_mpi) THEN 
-			  CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
+        CALL MPI_ALLREDUCE(dumax,dumax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
         dumax=dumax1
       ENDIF
@@ -291,5 +254,5 @@
           u=f_u(ind)
           du=f_du(ind)
-          CALL compute_gradrot(u,du,1)
+          CALL compute_gradrot(u,du,1,1)
           u=du
         ENDDO
@@ -362,9 +325,8 @@
           theta=f_theta(ind)
           dtheta=f_dtheta(ind)
-          CALL compute_divgrad(theta,dtheta,1)
+          CALL compute_divgrad(theta,dtheta,1,1)
           theta=dtheta
         ENDDO
       ENDDO
-!      CALL writefield("divgrad",f_dtheta)
 
       CALL transfert_request(f_dtheta,req_i1)
@@ -398,14 +360,8 @@
     ENDDO 
 
-!    CALL writefield("divgrad",f_dtheta)
     IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
   
     cdivgrad=dthetamax**(-1./niterdivgrad) 
     IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
-    
-
-
-
-
 
      
@@ -442,4 +398,5 @@
   USE trace
   USE time_mod
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -458,17 +415,19 @@
     INTEGER :: ind, shear
     INTEGER :: l,i,j,n
+
+!$OMP BARRIER
     
     CALL trace_start("dissip")
     CALL gradiv(f_ue,f_due_diss1)
     CALL gradrot(f_ue,f_due_diss2)
-    CALL theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
-    CALL divgrad(f_theta,f_dtheta_diss)
-    CALL theta2theta_rhodz(f_ps,f_dtheta_diss,f_dtheta_rhodz_diss)
-    
-    IF(rayleigh_friction_type>0) THEN
-       CALL pression(f_ps, f_p)
-       CALL exner(f_ps, f_p, f_pks, f_pk)
-       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
-    END IF
+    
+    CALL divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz_diss)
+    
+! later for openmp    
+!    IF(rayleigh_friction_type>0) THEN
+!       CALL pression(f_ps, f_p)
+!       CALL exner(f_ps, f_p, f_pks, f_pk)
+!       CALL geopotential(f_phis,f_pks,f_pk,f_theta,f_phi)
+!    END IF
 
     DO ind=1,ndomain
@@ -480,6 +439,6 @@
       dtheta_rhodz=f_dtheta_rhodz(ind)
       dtheta_rhodz_diss=f_dtheta_rhodz_diss(ind)
-            
-      DO l=1,llm
+
+      DO l=ll_begin,ll_end
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
@@ -495,21 +454,27 @@
       ENDDO
 
-      IF(rayleigh_friction_type>0) THEN
-         phi=f_phi(ind)
-         ue=f_ue(ind)
-         DO l=1,llm
-            DO j=jj_begin,jj_end
-               DO i=ii_begin,ii_end
-                  n=(j-1)*iim+i
-                  CALL relax(t_right, u_right)
-                  CALL relax(t_lup,   u_lup)
-                  CALL relax(t_ldown, u_ldown)
-               ENDDO
-            ENDDO
-         END DO
-      END IF
+!      dtheta_rhodz=0
+!      due=0
+
+! later for openmp  
+!      IF(rayleigh_friction_type>0) THEN
+!         phi=f_phi(ind)
+!         ue=f_ue(ind)
+!         DO l=1,llm
+!            DO j=jj_begin,jj_end
+!               DO i=ii_begin,ii_end
+!                  n=(j-1)*iim+i
+!                  CALL relax(t_right, u_right)
+!                  CALL relax(t_lup,   u_lup)
+!                  CALL relax(t_ldown, u_ldown)
+!               ENDDO
+!            ENDDO
+!         END DO
+!      END IF
    END DO
 
    CALL trace_end("dissip")
+
+!$OMP BARRIER
 
     CONTAINS
@@ -524,6 +489,5 @@
         z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
         IF(z>zh) THEN  ! relax only in the sponge layer z>zh
-!           PRINT *, 'z>zh : z,zh,l',z,zh,l
-!           STOP
+
            nn = n+shift_u
            CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
@@ -548,4 +512,5 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -554,5 +519,5 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it
+    INTEGER :: it,i,j,l,ij
        
     CALL trace_start("gradiv")
@@ -563,10 +528,20 @@
       ue=f_ue(ind)
       due=f_due(ind) 
-      due=ue
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+             due(ij+u_right,l)=ue(ij+u_right,l)
+             due(ij+u_lup,l)=ue(ij+u_lup,l)
+             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
+          ENDDO
+        ENDDO
+      ENDDO
     ENDDO
 
     DO it=1,nitergdiv
         
-      CALL transfert_request(f_due,req_e1_vect)
+      CALL send_message(f_due,req_due)
+      CALL wait_message(req_due)
         
       DO ind=1,ndomain
@@ -574,5 +549,5 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradiv(due,due,llm)
+        CALL compute_gradiv(due,due,ll_begin,ll_end)
       ENDDO
     ENDDO
@@ -586,4 +561,5 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ue(:)
@@ -592,5 +568,5 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it
+    INTEGER :: it,i,j,l,ij
        
     CALL trace_start("gradrot")
@@ -601,10 +577,20 @@
       ue=f_ue(ind)
       due=f_due(ind) 
-      due=ue
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+             due(ij+u_right,l)=ue(ij+u_right,l)
+             due(ij+u_lup,l)=ue(ij+u_lup,l)
+             due(ij+u_ldown,l)=ue(ij+u_ldown,l)
+          ENDDO
+        ENDDO
+      ENDDO
     ENDDO
 
     DO it=1,nitergrot
         
-      CALL transfert_request(f_due,req_e1_vect)
+      CALL send_message(f_due,req_due)
+      CALL wait_message(req_due)
         
       DO ind=1,ndomain
@@ -612,5 +598,5 @@
         CALL swap_geometry(ind)
         due=f_due(ind) 
-        CALL compute_gradrot(due,due,llm)
+        CALL compute_gradrot(due,due,ll_begin,ll_end)
       ENDDO
 
@@ -624,4 +610,5 @@
   USE icosa
   USE trace
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_theta(:)
@@ -650,5 +637,5 @@
         CALL swap_geometry(ind)
         dtheta=f_dtheta(ind) 
-        CALL compute_divgrad(dtheta,dtheta,llm)
+        CALL compute_divgrad(dtheta,dtheta,ll_begin,ll_end)
       ENDDO
 
@@ -659,97 +646,84 @@
   END SUBROUTINE divgrad
    
-  
-
-!  SUBROUTINE compute_dissip(ue,due,ps,theta_rhodz,dtheta_rhodz)
-!  USE icosa
-!  USE theta2theta_rhodz_mod
-!  IMPLICIT NONE
-!    REAL(rstd)         :: ue(3*iim*jjm,llm)
-!    REAL(rstd)         :: due(3*iim*jjm,llm)
-!    REAL(rstd)         :: ps(iim*jjm)
-!    REAL(rstd)         :: theta_rhodz(iim*jjm,llm)
-!    REAL(rstd)         :: dtheta_rhodz(iim*jjm,llm)
-!
-!    REAL(rstd),SAVE,ALLOCATABLE :: theta(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: du_dissip(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_dissip(:,:)
-!    REAL(rstd),SAVE,ALLOCATABLE :: dtheta_rhodz_dissip(:,:)
-!
-!    INTEGER :: ind
-!    INTEGER :: l,i,j,n
-!
-!!$OMP BARRIER
-!!$OMP MASTER
-!      ALLOCATE(theta(iim*jjm,llm))
-!      ALLOCATE(du_dissip(3*iim*jjm,llm))
-!      ALLOCATE(dtheta_dissip(iim*jjm,llm))
-!      ALLOCATE(dtheta_rhodz_dissip(iim*jjm,llm))
-!!$OMP END MASTER
-!!$OMP BARRIER
-!    
-!      CALL gradiv(ue,du_dissip,llm)
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_graddiv(l)*0.5
-!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_graddiv(l)*0.5
-!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_graddiv(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO
-!      
-!      CALL gradrot(ue,du_dissip,llm)
-!
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            due(n+u_right,l)=due(n+u_right,l)+du_dissip(n+u_right,l)/tau_gradrot(l)*0.5
-!            due(n+u_lup,l)=due(n+u_lup,l)+du_dissip(n+u_lup,l)/tau_gradrot(l)*0.5
-!            due(n+u_ldown,l)=due(n+u_ldown,l)+du_dissip(n+u_ldown,l)/tau_gradrot(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO 
-!      
-!      CALL compute_theta_rhodz2theta(ps,theta_rhodz,theta,1)
-!      CALL divgrad(theta,dtheta_dissip,llm)
-!      CALL compute_theta2theta_rhodz(ps,dtheta_dissip,dtheta_rhodz_dissip,0)
-!
-!      DO l=1,llm
-!!$OMP DO
-!        DO j=jj_begin,jj_end
-!          DO i=ii_begin,ii_end
-!            n=(j-1)*iim+i
-!            dtheta_rhodz(n,l)=dtheta_rhodz(n,l)+dtheta_rhodz_dissip(n,l)/tau_divgrad(l)*0.5
-!          ENDDO
-!        ENDDO
-!      ENDDO 
-!
-!!$OMP BARRIER
-!!$OMP MASTER
-!      DEALLOCATE(theta)
-!      DEALLOCATE(du_dissip)
-!      DEALLOCATE(dtheta_dissip)
-!      DEALLOCATE(dtheta_rhodz_dissip)
-!!$OMP END MASTER
-!!$OMP BARRIER      
-!
-!  END SUBROUTINE compute_dissip
-  
-  
-  SUBROUTINE compute_gradiv(ue,gradivu_e,ll)
-  USE icosa
+  SUBROUTINE divgrad_theta_rhodz(f_ps,f_theta_rhodz,f_dtheta_rhodz)
+  USE icosa
+  USE trace
+  USE omp_para
+  USE disvert_mod
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,ll)
-    REAL(rstd) :: divu_i(iim*jjm,ll)
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_dtheta_rhodz(:)
+    
+    REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: dtheta_rhodz(:,:)
+
+    INTEGER :: ind
+    INTEGER :: it,i,j,l,ij
+
+    CALL trace_start("divgrad")
+       
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind) 
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+
+    DO it=1,niterdivgrad
+        
+      CALL send_message(f_dtheta_rhodz,req_dtheta)
+      CALL wait_message(req_dtheta)
+      DO ind=1,ndomain
+        CALL swap_dimensions(ind)
+        CALL swap_geometry(ind)
+        dtheta_rhodz=f_dtheta_rhodz(ind) 
+        CALL compute_divgrad(dtheta_rhodz,dtheta_rhodz,ll_begin,ll_end)
+      ENDDO
+
+    ENDDO
+
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      dtheta_rhodz=f_dtheta_rhodz(ind) 
+    
+      DO  l = ll_begin, ll_end
+        DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+            ij=(j-1)*iim+i
+            dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * ( ap(l)-ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/g
+          ENDDO
+        ENDDO
+      ENDDO
+    ENDDO
+
+
+    CALL trace_end("divgrad")
+
+  END SUBROUTINE divgrad_theta_rhodz  
+  
+  SUBROUTINE compute_gradiv(ue,gradivu_e,llb,lle)
+  USE icosa
+  IMPLICIT NONE
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: gradivu_e(iim*3*jjm,llm)
+    REAL(rstd) :: divu_i(iim*jjm,llb:lle)
     
     INTEGER :: i,j,n,l
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -765,5 +739,5 @@
     ENDDO
     
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -781,5 +755,5 @@
     ENDDO
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -795,16 +769,17 @@
   END SUBROUTINE compute_gradiv
   
-  SUBROUTINE compute_divgrad(theta,divgrad_i,ll)
+  SUBROUTINE compute_divgrad(theta,divgrad_i,llb,lle)
   USE icosa
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,ll)
-    REAL(rstd)  :: grad_e(3*iim*jjm,ll)
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: divgrad_i(iim*jjm,llm)
+    REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
 
     INTEGER :: i,j,n,l
 
        
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -823,5 +798,5 @@
     
     
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -837,5 +812,5 @@
     ENDDO
    
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -849,15 +824,16 @@
 
     
-  SUBROUTINE compute_gradrot(ue,gradrot_e,ll)
+  SUBROUTINE compute_gradrot(ue,gradrot_e,llb,lle)
   USE icosa
   IMPLICIT NONE
-    INTEGER,INTENT(IN)     :: ll
-    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,ll)
-    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,ll)
-    REAL(rstd) :: rot_v(2*iim*jjm,ll)
+    INTEGER,INTENT(IN)     :: llb
+    INTEGER,INTENT(IN)     :: lle
+    REAL(rstd),INTENT(IN)  :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: gradrot_e(iim*3*jjm,llm)
+    REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
 
     INTEGER :: i,j,n,l
      
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin-1,jj_end+1
         DO i=ii_begin-1,ii_end+1
@@ -876,5 +852,5 @@
     ENDDO                              
   
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -891,5 +867,5 @@
     ENDDO
 
-    DO l=1,ll
+    DO l=llb,lle
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
Index: codes/icosagcm/trunk/src/domain.f90
===================================================================
--- codes/icosagcm/trunk/src/domain.f90	(revision 150)
+++ codes/icosagcm/trunk/src/domain.f90	(revision 151)
@@ -71,4 +71,5 @@
   SUBROUTINE create_domain
   USE grid_param
+  USE mpipara
   IMPLICIT NONE
   INTEGER :: ind,nf,ni,nj,i,j
@@ -98,7 +99,7 @@
           d%ii_end_glo=d%ii_begin_glo+d%ii_nb-1
  
-          IF (ni/=nsplit_i) THEN 
+          IF (ni/=1) THEN 
             d%ii_nb=d%ii_nb+1
-            d%ii_end_glo=d%ii_end_glo+1
+            d%ii_begin_glo=d%ii_begin_glo-1
           ENDIF
          
@@ -115,7 +116,7 @@
           d%jj_end_glo=d%jj_begin_glo+d%jj_nb-1
 
-          IF (nj/=nsplit_j) THEN 
+          IF (nj/=1) THEN 
             d%jj_nb=d%jj_nb+1
-            d%jj_end_glo=d%jj_end_glo+1
+            d%jj_begin_glo=d%jj_begin_glo-1
           ENDIF
 
@@ -142,4 +143,7 @@
           ALLOCATE(d%own(d%iim,d%jjm))
           ALLOCATE(d%ne(0:5,d%iim,d%jjm))
+          
+          IF (is_mpi_root) PRINT *,"Domain ",ind," : size ",d%ii_nb,"x",d%jj_nb
+          
         END DO
       END DO
@@ -261,5 +265,4 @@
             d%neighbour(:,k,i,j)=cell_glo(ind2)%xyz(:)
 
-!            d%ne(k,i,j)=vertex_glo(ii,jj,nf)%ne(k)
             d%ne(k,i,j)=1-2*MOD(k,2)
 
@@ -298,5 +301,5 @@
       edge_glo(e)%assign_pos=k
       edge_glo(e)%assign_delta=delta
-!      PRINT*,"Assign_edge",ind_d,ind,i,j,k,e
+
      END  SUBROUTINE assign_edge
          
Index: codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90	(revision 150)
+++ codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90	(revision 151)
@@ -18,7 +18,6 @@
 CONTAINS
 
-	SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
+  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
   USE icosa
-  USE caldyn_mod
   USE write_field  
   USE maxicosa
@@ -98,6 +97,4 @@
 !----------------------------------------------------
 !------------- OUTPUT OF Variables 
-	CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
-	CALL writefield("buf1",f_buf1_i)
   END SUBROUTINE etat0
 
Index: codes/icosagcm/trunk/src/dynetat0_hz_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_hz_mod.f90	(revision 150)
+++ codes/icosagcm/trunk/src/dynetat0_hz_mod.f90	(revision 151)
@@ -2,5 +2,4 @@
   USE genmod
   USE icosa
-  USE caldyn_gcm_mod 
   IMPLICIT NONE 
   PRIVATE
@@ -23,4 +22,5 @@
   USE write_field  
   USE maxicosa
+  USE wind_mod
         IMPLICIT NONE
      TYPE(t_domain),POINTER :: d 
@@ -97,5 +97,5 @@
 	ENDDO
 !------------- OUTPUT OF Variables 
-	CALL un2ulonlat(f_u,f_buf_i3,f_buf1_i,f_buf2_i) 
+	CALL un2ulonlat(f_u,f_buf1_i,f_buf2_i) 
 	CALL writefield("buf1",f_buf1_i)
   END SUBROUTINE etat0
Index: codes/icosagcm/trunk/src/exner.f90
===================================================================
--- codes/icosagcm/trunk/src/exner.f90	(revision 150)
+++ codes/icosagcm/trunk/src/exner.f90	(revision 151)
@@ -50,5 +50,4 @@
        
        ! for llm layer
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -61,5 +60,4 @@
        ! for other layer   
        DO l = llm-1 , 2 , -1
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
@@ -75,5 +73,4 @@
        
        ! for first layer
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -88,5 +85,4 @@
        
        DO l = 2, llm
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
@@ -99,5 +95,4 @@
     ELSE ! Simple calculation of Exner pressure based on centered average
        ! surface : pks
-       !$OMP DO
        DO j=jj_begin-offset,jj_end+offset
           DO i=ii_begin-offset,ii_end+offset
@@ -109,5 +104,4 @@
        ! 3D : pk
        DO l = 1, llm
-          !$OMP DO
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
Index: codes/icosagcm/trunk/src/geometry.f90
===================================================================
--- codes/icosagcm/trunk/src/geometry.f90	(revision 150)
+++ codes/icosagcm/trunk/src/geometry.f90	(revision 151)
@@ -91,4 +91,5 @@
   IMPLICIT NONE
     INTEGER,INTENT(IN) :: ind
+!$OMP MASTER
     Ai=geom%Ai(ind)
     xyz_i=geom%xyz_i(ind)
@@ -111,5 +112,6 @@
     bi=geom%bi(ind)
     fv=geom%fv(ind)
-    
+!$OMP END MASTER
+!$OMP BARRIER    
   END SUBROUTINE swap_geometry
   
@@ -213,15 +215,8 @@
           ENDDO
         ENDDO
-!        i=ii_begin ; j=jj_begin ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j) 
-!        i=ii_begin ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
-!        i=ii_end   ; j=jj_begin ; n=(j-1)*iim+i ;   xyz_i(n,:)=domain(ind)%xyz(:,i,j)
-!        PRINT *,"Pb ?? : " 
-!        PRINT *, xyz_i(n,:), domain(ind)%xyz(:,i,j), norm(xyz_i(n,:)- domain(ind)%xyz(:,i,j))
-!        i=ii_end   ; j=jj_end   ; n=(j-1)*iim+i ; xyz_i(n,:)=domain(ind)%xyz(:,i,j)
         
       ENDDO
 
       IF (check) THEN
-!        sum=sum/(ndomain*ii_nb*jj_nb)
         PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
       ENDIF 
@@ -287,8 +282,4 @@
             ne(n,k+1)=d%ne(k,i,j)
           ENDDO
-                  
-!          xyz_i(n,:)=d%xyz(:,i,j) 
-!          xyz_v(n+z_up,:)=d%vertex(:,vup-1,i,j) 
-!          xyz_v(n+z_down,:)=d%vertex(:,vdown-1,i,j) 
 
           vect(:,1)=xyz_v(n+z_rup,:)
@@ -309,23 +300,11 @@
           bi(n)=0. 
         
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),de(n+u_right))
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),de(n+u_lup))
-!          CALL dist_cart(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),de(n+u_ldown))
-
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
           CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
           
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,right-1,i,j),0.5,xyz_e(n+u_right,:))
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,lup-1,i,j),0.5,xyz_e(n+u_lup,:))
-!          CALL div_arc_bis(d%xyz(:,i,j),d%neighbour(:,ldown-1,i,j),0.5,xyz_e(n+u_ldown,:))
-
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
           CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
-         
-!          CALL dist_cart(d%vertex(:,vrdown-1,i,j),d%vertex(:,vrup-1,i,j),le(n+u_right))
-!          CALL dist_cart(d%vertex(:,vup-1,i,j),d%vertex(:,vlup-1,i,j),le(n+u_lup))
-!          CALL dist_cart(d%vertex(:,vldown-1,i,j),d%vertex(:,vdown-1,i,j),le(n+u_ldown))
 
           CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
@@ -335,6 +314,4 @@
           Ai(n)=0
           DO k=0,5
-!            CALL surf_triangle(d%xyz(:,i,j),d%neighbour(:,k,i,j),d%neighbour(:,MOD((k+1+6),6),i,j),surf_v(k+1))
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,MOD((k-1+6),6),i,j),d%vertex(:,k,i,j),surf(k+1))
             CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
             CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
@@ -414,6 +391,4 @@
           
           DO k=0,5
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,k,i,j),S1)
-!            CALL surf_triangle(d%xyz(:,i,j),d%vertex(:,k,i,j),d%neighbour(:,MOD(k+1+6,6),i,j),S2)
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
             CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
@@ -516,6 +491,5 @@
     CALL transfert_request(geom%centroid,req_i1)
     CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
-!    PRINT *,"Surf triangle : ",S*20/(4*Pi)
-  
+ 
   END SUBROUTINE set_geometry
   
Index: codes/icosagcm/trunk/src/icosa_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 150)
+++ codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 151)
@@ -6,5 +6,7 @@
   USE wind_mod
   USE mpipara
+  USE omp_para
   USE vertical_interp_mod
+  USE trace
   IMPLICIT NONE
   
@@ -21,30 +23,11 @@
   CALL init_earth_const 
   CALL init_grid_param
+  CALL init_omp_para
   CALL compute_metric
   CALL compute_domain
   CALL init_transfert
   CALL init_writefield
+  CALL init_trace
 
-!  CALL allocate_field(sum_ne,field_T,type_real)
-!  CALL allocate_field_glo(sum_ne_glo,field_T,type_real)
-!  
-! DO ind=1,ndomain
-!   CALL swap_dimensions(ind)
-!   pt_sum_ne=sum_ne(ind)
-!   DO j=jj_begin,jj_end
-!     DO i=ii_begin,ii_end    
-!       n=(j-1)*iim+i
-!       pt_sum_ne(n)=domloc_glo_ind(ind)
-!     ENDDO
-!   ENDDO
-! ENDDO
-
-! CALL WriteField("domain",sum_ne)
-! CALL WriteField_mpi("domain",sum_ne)
-! CALL transfert_request(sum_ne,req_i1)
-! CALL WriteField_mpi("domain",sum_ne)
-! CALL close_files
-! CALL finalize_mpipara
-! STOP
   
   CALL compute_geometry
@@ -80,8 +63,13 @@
   
   CALL WriteField("Ai",geom%Ai)
-!  CALL WriteField("sum_ne",sum_ne)
+
   IF (is_mpi_root) CALL write_apbp
   CALL init_time
+
+  CALL init_timeloop
+
+!$OMP PARALLEL
   CALL timeloop
+!$OMP END PARALLEL
 
   CALL close_files
Index: codes/icosagcm/trunk/src/mpi_mod.F90
===================================================================
--- codes/icosagcm/trunk/src/mpi_mod.F90	(revision 150)
+++ codes/icosagcm/trunk/src/mpi_mod.F90	(revision 151)
@@ -11,4 +11,5 @@
   INTEGER :: MPI_INFO_NULL
   INTEGER :: MPI_STATUS_SIZE 
+  INTEGER,PARAMETER :: MPI_ADDRESS_KIND=KIND(INTEGER)
 #endif
 
@@ -40,8 +41,14 @@
  END
 
+ SUBROUTINE MPI_ISSEND
+ END
+
  SUBROUTINE MPI_IRECV
  END
 
  SUBROUTINE MPI_WAITALL
+ END
+
+ SUBROUTINE MPI_TESTALL
  END
 
@@ -51,4 +58,10 @@
  SUBROUTINE MPI_ALLGATHER
  END
+ 
+ SUBROUTINE MPI_TYPE_EXTENT
+ END
+ 
+ SUBROUTINE MPI_ALLOC_MEM
+ END
 
 #endif
Index: codes/icosagcm/trunk/src/mpipara.F90
===================================================================
--- codes/icosagcm/trunk/src/mpipara.F90	(revision 150)
+++ codes/icosagcm/trunk/src/mpipara.F90	(revision 151)
@@ -8,4 +8,8 @@
   LOGICAL,SAVE :: using_mpi
   LOGICAL,SAVE :: is_mpi_root
+  
+  INTERFACE allocate_mpi_buffer
+    MODULE PROCEDURE allocate_mpi_buffer_r2, allocate_mpi_buffer_r3,allocate_mpi_buffer_r4
+  END INTERFACE allocate_mpi_buffer
 
 CONTAINS
@@ -48,4 +52,67 @@
    END SUBROUTINE  finalize_mpipara
    
+
+  SUBROUTINE allocate_mpi_buffer_r2(buffer,length)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+  IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:)
+    INTEGER,INTENT(IN)  :: length
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length /))
+
+   END SUBROUTINE allocate_mpi_buffer_r2
+
+  SUBROUTINE allocate_mpi_buffer_r3(buffer,length,dim3)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+    IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:,:)
+    INTEGER,INTENT(IN)  :: length
+    INTEGER,INTENT(IN)  :: dim3
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size*dim3
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length,dim3 /))
+    
+   END SUBROUTINE allocate_mpi_buffer_r3
+
+  SUBROUTINE allocate_mpi_buffer_r4(buffer,length,dim3,dim4)
+  USE ISO_C_BINDING
+  USE mpi_mod
+  USE prec
+  IMPLICIT NONE
+    REAL(rstd), POINTER :: buffer(:,:,:)
+    INTEGER,INTENT(IN)  :: length
+    INTEGER,INTENT(IN)  :: dim3
+    INTEGER,INTENT(IN)  :: dim4
+
+    TYPE(C_PTR)         :: base_ptr
+    INTEGER(KIND=MPI_ADDRESS_KIND) :: size
+    INTEGER :: real_size,ierr
+    
+    CALL MPI_TYPE_EXTENT(MPI_REAL8, real_size, ierr)
+    size=length*real_size*dim3*dim4
+    
+    CALL MPI_ALLOC_MEM(size,MPI_INFO_NULL,base_ptr,ierr)
+    CALL C_F_POINTER(base_ptr, buffer, (/ length, dim3, dim4 /))
+    
+   END SUBROUTINE allocate_mpi_buffer_r4
    
 END MODULE mpipara
Index: codes/icosagcm/trunk/src/omp_para.F90
===================================================================
--- codes/icosagcm/trunk/src/omp_para.F90	(revision 151)
+++ codes/icosagcm/trunk/src/omp_para.F90	(revision 151)
@@ -0,0 +1,100 @@
+MODULE omp_para
+
+  INTEGER,SAVE :: omp_size
+  INTEGER,SAVE :: omp_rank
+!$OMP THREADPRIVATE(omp_rank)
+
+  LOGICAL,SAVE :: omp_first
+  LOGICAL,SAVE :: omp_last
+  LOGICAL,SAVE :: omp_master
+!$OMP THREADPRIVATE(omp_first, omp_last,omp_master)
+
+  INTEGER,SAVE :: ll_begin
+  INTEGER,SAVE :: ll_beginp1
+  INTEGER,SAVE :: ll_end
+  INTEGER,SAVE :: ll_endm1
+  INTEGER,SAVE :: ll_endp1
+!$OMP THREADPRIVATE(ll_begin,ll_beginp1,ll_end,ll_endm1,ll_endp1)
+  LOGICAL,SAVE :: using_openmp
+
+CONTAINS
+
+
+  SUBROUTINE init_omp_para
+  USE grid_param
+#ifdef CPP_USING_OPENMP
+  USE omp_lib
+#endif
+  IMPLICIT NONE
+  INTEGER :: ll_nb,i
+
+#ifdef CPP_USING_OPENMP
+  using_openmp=.TRUE.
+#else
+  using_openmp=.FALSE.
+#endif
+
+  IF (using_openmp) THEN    
+!$OMP PARALLEL PRIVATE(ll_nb,i)
+  
+!$OMP MASTER
+#ifdef CPP_USING_OPENMP
+    omp_size=OMP_GET_NUM_THREADS()
+#endif
+!$OMP END MASTER
+!$OMP BARRIER
+#ifdef CPP_USING_OPENMP
+    omp_rank=OMP_GET_THREAD_NUM()
+#endif
+    omp_first=.FALSE.
+    omp_last=.FALSE.
+    omp_master=.FALSE.
+    
+    IF (omp_rank==0) THEN
+      omp_first=.TRUE.
+      omp_master=.TRUE.
+    ENDIF
+    
+    IF (omp_rank==omp_size-1) omp_last=.TRUE.
+    
+    ll_end=0
+    DO i=0,omp_rank
+      ll_begin=ll_end+1
+      ll_nb=llm/omp_size
+      IF (MOD(llm,omp_size)>i) ll_nb=ll_nb+1
+      ll_end=ll_begin+ll_nb-1
+    ENDDO
+    
+    ll_beginp1=ll_begin
+    ll_endp1=ll_end
+    ll_endm1=ll_end
+
+    IF (omp_first) ll_beginp1=ll_begin+1
+    IF (omp_last) ll_endp1=ll_endp1+1
+    IF (omp_last) ll_endm1=ll_endm1-1
+
+!$OMP END PARALLEL
+
+   ELSE
+     omp_size=1
+     omp_rank=0
+     omp_first=.TRUE.
+     omp_last=.TRUE.
+     omp_master=.TRUE.
+     ll_begin=1
+     ll_beginp1=2
+     ll_end=llm
+     ll_endm1=llm-1
+     ll_endp1=llm+1
+   ENDIF
+   
+  END SUBROUTINE init_omp_para
+
+  
+END MODULE omp_para
+    
+    
+      
+             
+  
+    
Index: codes/icosagcm/trunk/src/pression.f90
===================================================================
--- codes/icosagcm/trunk/src/pression.f90	(revision 150)
+++ codes/icosagcm/trunk/src/pression.f90	(revision 151)
@@ -33,5 +33,4 @@
 
     DO    l    = 1, llm+1
-!$OMP DO
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
Index: codes/icosagcm/trunk/src/theta_rhodz.f90
===================================================================
--- codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 150)
+++ codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 151)
@@ -83,14 +83,8 @@
     REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)
 
-!$OMP BARRIER
-!$OMP MASTER
-      ALLOCATE( p(iim*jjm,llm+1))
-!$OMP END MASTER
-!$OMP BARRIER
-
+    ALLOCATE( p(iim*jjm,llm+1))
     CALL compute_pression(ps,p,offset)
     
     DO    l    = 1, llm
-!$OMP DO
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
@@ -101,9 +95,5 @@
     ENDDO
 
-!$OMP BARRIER
-!$OMP MASTER
-      DEALLOCATE( p)
-!$OMP END MASTER
-!$OMP BARRIER    
+    DEALLOCATE( p)
 
   END SUBROUTINE compute_theta2theta_rhodz
@@ -120,9 +110,6 @@
     REAL(rstd),SAVE,ALLOCATABLE :: p(:,:)
 
-!$OMP BARRIER
-!$OMP MASTER
-      ALLOCATE( p(iim*jjm,llm+1))
-!$OMP END MASTER
-!$OMP BARRIER
+    ALLOCATE( p(iim*jjm,llm+1))
+
     CALL compute_pression(ps,p,offset)
     
@@ -136,9 +123,5 @@
     ENDDO
 
-!$OMP BARRIER
-!$OMP MASTER
-      DEALLOCATE( p)
-!$OMP END MASTER
-!$OMP BARRIER    
+    DEALLOCATE( p)
     
   END SUBROUTINE compute_theta_rhodz2theta
@@ -193,5 +176,4 @@
         DO i=ii_begin-offset,ii_end+offset
           ij=(j-1)*iim+i
-!          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp
           theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
         ENDDO
Index: codes/icosagcm/trunk/src/time.f90
===================================================================
--- codes/icosagcm/trunk/src/time.f90	(revision 150)
+++ codes/icosagcm/trunk/src/time.f90	(revision 151)
@@ -119,4 +119,6 @@
   INTEGER :: status
   REAL(rstd) ::time_array(1)
+
+!$OMP MASTER
     time_array(1)=time
   
@@ -126,4 +128,5 @@
       status=NF90_SYNC(ncid)
     ENDIF
+!$OMP END MASTER
 
   END SUBROUTINE update_time_counter    
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 150)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 151)
@@ -1,14 +1,33 @@
 MODULE timeloop_gcm_mod
-
+  USE transfert_mod
+  USE icosa
   PRIVATE
   
-  PUBLIC :: timeloop
+  PUBLIC :: init_timeloop, timeloop
 
   INTEGER, PARAMETER :: euler=1, rk4=2, mlf=3
   INTEGER  :: itau_sync=10
 
+  TYPE(t_message) :: req_ps0, req_theta_rhodz0, req_u0, req_q0
+
+  TYPE(t_field),POINTER :: f_phis(:)
+  TYPE(t_field),POINTER :: f_q(:)
+  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
+  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
+  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
+  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
+  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
+  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
+  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)  
+
+  INTEGER :: nb_stage, matsuno_period, scheme    
+
+  REAL(rstd),SAVE :: jD_cur, jH_cur
+  REAL(rstd),SAVE :: start_time
+
 CONTAINS
   
-  SUBROUTINE timeloop
+  SUBROUTINE init_timeloop
   USE icosa
   USE dissip_gcm_mod
@@ -19,77 +38,34 @@
   USE physics_mod
   USE mpipara
-  USE IOIPSL
-  USE maxicosa
-  USE check_conserve_mod
+  USE omp_para
   USE trace
   USE transfert_mod
+  USE check_conserve_mod
+  USE ioipsl
   IMPLICIT NONE
-  TYPE(t_field),POINTER :: f_phis(:)
-!  TYPE(t_field),POINTER :: f_theta(:)
-  TYPE(t_field),POINTER :: f_q(:)
-  TYPE(t_field),POINTER :: f_dtheta(:), f_rhodz(:)
-  TYPE(t_field),POINTER :: f_ps(:),f_psm1(:), f_psm2(:)
-  TYPE(t_field),POINTER :: f_u(:),f_um1(:),f_um2(:)
-  TYPE(t_field),POINTER :: f_theta_rhodz(:),f_theta_rhodzm1(:),f_theta_rhodzm2(:)
-  TYPE(t_field),POINTER :: f_dps(:),f_dpsm1(:), f_dpsm2(:)
-  TYPE(t_field),POINTER :: f_du(:),f_dum1(:),f_dum2(:)
-  TYPE(t_field),POINTER :: f_dtheta_rhodz(:),f_dtheta_rhodzm1(:),f_dtheta_rhodzm2(:)
-  TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:), f_hfluxt(:), f_wfluxt(:)
-
-  REAL(rstd),POINTER :: phis(:)
-  REAL(rstd),POINTER :: q(:,:,:)
-  REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
-  REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
-  REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
-  REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
-  REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
-  REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
-  REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
-!  REAL(rstd) :: dt, run_length
-  INTEGER :: ind
-  INTEGER :: it,i,j,n, nb_stage, stage, matsuno_period, scheme
-  CHARACTER(len=255) :: scheme_name
-  LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
-  CHARACTER(len=7) :: first 
-  REAL(rstd),SAVE :: jD_cur, jH_cur
-  REAL(rstd),SAVE :: start_time
-  LOGICAL, PARAMETER :: check=.FALSE.
-!  INTEGER :: itaumax
-!  REAL(rstd) ::write_period
-!  INTEGER    :: itau_out
-  
-!  dt=90.
-!  CALL getin('dt',dt)
-  
-!  itaumax=100
-!  CALL getin('itaumax',itaumax)
-
-!  run_length=dt*itaumax
-!  CALL getin('run_length',run_length)
-!  itaumax=run_length/dt
-!  PRINT *,'itaumax=',itaumax
-!  dt=dt/scale_factor
-
-! Trends
-  CALL allocate_field(f_dps,field_t,type_real)
-  CALL allocate_field(f_du,field_u,type_real,llm)
-  CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
+
+    CHARACTER(len=255) :: scheme_name
+
+    CALL allocate_field(f_dps,field_t,type_real)
+    CALL allocate_field(f_du,field_u,type_real,llm)
+    CALL allocate_field(f_dtheta_rhodz,field_t,type_real,llm)
 ! Model state at current time step (RK/MLF/Euler)
-  CALL allocate_field(f_phis,field_t,type_real)
-  CALL allocate_field(f_ps,field_t,type_real)
-  CALL allocate_field(f_u,field_u,type_real,llm)
-  CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_phis,field_t,type_real)
+    CALL allocate_field(f_ps,field_t,type_real)
+    CALL allocate_field(f_u,field_u,type_real,llm)
+    CALL allocate_field(f_theta_rhodz,field_t,type_real,llm)
 ! Model state at previous time step (RK/MLF)
-  CALL allocate_field(f_psm1,field_t,type_real)
-  CALL allocate_field(f_um1,field_u,type_real,llm)
-  CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
+    CALL allocate_field(f_psm1,field_t,type_real)
+    CALL allocate_field(f_um1,field_u,type_real,llm)
+    CALL allocate_field(f_theta_rhodzm1,field_t,type_real,llm)
 ! Tracers
-  CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
-  CALL allocate_field(f_rhodz,field_t,type_real,llm)
+    CALL allocate_field(f_q,field_t,type_real,llm,nqtot)
+    CALL allocate_field(f_rhodz,field_t,type_real,llm)
 ! Mass fluxes
-  CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
-  CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
-  CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
-  CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+    CALL allocate_field(f_hflux,field_u,type_real,llm)    ! instantaneous mass fluxes
+    CALL allocate_field(f_wflux,field_t,type_real,llm+1)  ! computed by caldyn
+    CALL allocate_field(f_hfluxt,field_u,type_real,llm)   ! mass "fluxes" accumulated in time
+    CALL allocate_field(f_wfluxt,field_t,type_real,llm+1)
+
 
 !----------------------------------------------------
@@ -145,93 +121,123 @@
 !----------------------------------------------------
 
-  scheme_name='runge_kutta'
-  CALL getin('scheme',scheme_name)
-
-  SELECT CASE (TRIM(scheme_name))
-  CASE('euler')
-     scheme=euler
-     nb_stage=1
-  CASE ('runge_kutta')
-     scheme=rk4
-     nb_stage=4
-  CASE ('leapfrog_matsuno')
-     scheme=mlf
-     matsuno_period=5
-     CALL getin('matsuno_period',matsuno_period)
-     nb_stage=matsuno_period+1
+
+
+
+
+    scheme_name='runge_kutta'
+    CALL getin('scheme',scheme_name)
+
+    SELECT CASE (TRIM(scheme_name))
+      CASE('euler')
+        scheme=euler
+        nb_stage=1
+      CASE ('runge_kutta')
+        scheme=rk4
+        nb_stage=4
+      CASE ('leapfrog_matsuno')
+        scheme=mlf
+        matsuno_period=5
+        CALL getin('matsuno_period',matsuno_period)
+        nb_stage=matsuno_period+1
      ! Model state 2 time steps ago (MLF)
-     CALL allocate_field(f_psm2,field_t,type_real)
-     CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
-     CALL allocate_field(f_um2,field_u,type_real,llm)
-  CASE default
-     PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
-          ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
-     STOP
-  END SELECT
+        CALL allocate_field(f_psm2,field_t,type_real)
+        CALL allocate_field(f_theta_rhodzm2,field_t,type_real,llm)
+        CALL allocate_field(f_um2,field_u,type_real,llm)
+    CASE default
+       PRINT*,'Bad selector for variable scheme : <', TRIM(scheme_name),             &
+            ' > options are <euler>, <runge_kutta>, <leapfrog_matsuno>'
+       STOP
+    END SELECT
+
+
+    CALL init_dissip
+    CALL init_caldyn
+    CALL init_guided
+    CALL init_advect_tracer
+    CALL init_check_conserve
+    CALL init_physics
+    CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
+
+    CALL transfert_request(f_phis,req_i0) 
+    CALL transfert_request(f_phis,req_i1) 
+    CALL writefield("phis",f_phis,once=.TRUE.)
+
+    CALL init_message(f_ps,req_i0,req_ps0)
+    CALL init_message(f_theta_rhodz,req_i0,req_theta_rhodz0)
+    CALL init_message(f_u,req_e0_vect,req_u0)
+    CALL init_message(f_q,req_i0,req_q0)
+
+  END SUBROUTINE init_timeloop
   
-!  write_period=0
-!  CALL getin('write_period',write_period)
-!  write_period=write_period/scale_factor
-!  itau_out=FLOOR(.5+write_period/dt)
-!  PRINT *, 'Output frequency (scaled) set to ',write_period, ' : itau_out = ',itau_out
- 
-! Trends at previous time steps needed only by Adams-Bashforth
-!  CALL allocate_field(f_dpsm1,field_t,type_real)
-!  CALL allocate_field(f_dpsm2,field_t,type_real)
-!  CALL allocate_field(f_dum1,field_u,type_real,llm)
-!  CALL allocate_field(f_dum2,field_u,type_real,llm)
-!  CALL allocate_field(f_dtheta_rhodzm1,field_t,type_real,llm)
-!  CALL allocate_field(f_dtheta_rhodzm2,field_t,type_real,llm)
-!  CALL allocate_field(f_theta,field_t,type_real,llm)
-!  CALL allocate_field(f_dtheta,field_t,type_real,llm)
-
-  CALL init_dissip
-  CALL init_caldyn
-  CALL init_guided
-  CALL init_advect_tracer
-  CALL init_physics
-!========================================= INITIALIZATION 
-! CALL dynetat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-  CALL etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
-  CALL writefield("phis",f_phis,once=.TRUE.)
-  CALL transfert_request(f_q,req_i1) 
-
+  SUBROUTINE timeloop
+  USE icosa
+  USE dissip_gcm_mod
+  USE caldyn_mod
+  USE etat0_mod
+  USE guided_mod
+  USE advect_tracer_mod
+  USE physics_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  USE transfert_mod
+  USE check_conserve_mod
+  IMPLICIT NONE  
+    REAL(rstd),POINTER :: phis(:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: ps(:) ,psm1(:), psm2(:)
+    REAL(rstd),POINTER :: u(:,:) , um1(:,:), um2(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:), theta_rhodz(:,:) , theta_rhodzm1(:,:), theta_rhodzm2(:,:)
+    REAL(rstd),POINTER :: dps(:), dpsm1(:), dpsm2(:)
+    REAL(rstd),POINTER :: du(:,:), dum1(:,:), dum2(:,:)
+    REAL(rstd),POINTER :: dtheta_rhodz(:,:),dtheta_rhodzm1(:,:),dtheta_rhodzm2(:,:)
+    REAL(rstd),POINTER :: hflux(:,:),wflux(:,:),hfluxt(:,:),wfluxt(:,:)
+
+    INTEGER :: ind
+    INTEGER :: it,i,j,n,  stage
+    CHARACTER(len=255) :: scheme_name
+    LOGICAL :: fluxt_zero(ndomain) ! set to .TRUE. to start accumulating fluxes in time
+    LOGICAL, PARAMETER :: check=.FALSE.
+
+!$OMP BARRIER
   DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(.TRUE., ps,rhodz) ! save rhodz for transport scheme before dynamics update ps
+     CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
   END DO
   fluxt_zero=.TRUE.
 
-  ! check that rhodz is consistent with ps
-  CALL transfert_request(f_rhodz,req_i1)
-  CALL transfert_request(f_ps,req_i1)
-  DO ind=1,ndomain
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind)
-     rhodz=f_rhodz(ind); ps=f_ps(ind)
-     CALL compute_rhodz(.FALSE., ps, rhodz)   
-  END DO
-
-  CALL transfert_request(f_phis,req_i0) 
-  CALL transfert_request(f_phis,req_i1) 
-  CALL transfert_request(f_phis,req_i1) 
+! check that rhodz is consistent with ps
+!  CALL transfert_request(f_rhodz,req_i1)
+!  CALL transfert_request(f_ps,req_i1)
+!  DO ind=1,ndomain
+!     CALL swap_dimensions(ind)
+!     CALL swap_geometry(ind)
+!     rhodz=f_rhodz(ind); ps=f_ps(ind)
+!     CALL compute_rhodz(.FALSE., ps, rhodz)   
+!  END DO
+  
 
   DO it=0,itaumax
     IF (MOD(it,itau_sync)==0) THEN
-      CALL transfert_request(f_ps,req_i0)
-      CALL transfert_request(f_theta_rhodz,req_i0) 
-      CALL transfert_request(f_u,req_e0_vect)
-      CALL transfert_request(f_q,req_i0) 
+      CALL send_message(f_ps,req_ps0)
+      CALL send_message(f_theta_rhodz,req_theta_rhodz0) 
+      CALL send_message(f_u,req_u0)
+      CALL send_message(f_q,req_q0) 
+      CALL wait_message(req_ps0)
+      CALL wait_message(req_theta_rhodz0) 
+      CALL wait_message(req_u0)
+      CALL wait_message(req_q0) 
     ENDIF
     
+!    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
     IF (mod(it,itau_out)==0 ) THEN
-!      IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
+      CALL writefield("q",f_q)
       CALL update_time_counter(dt*it)
-      CALL compute_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
+      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     ENDIF
-
-      CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
+    
+    CALL guided(it*dt,f_ps,f_theta_rhodz,f_u,f_q)
 
     DO stage=1,nb_stage
@@ -246,10 +252,11 @@
        CASE (mlf)
           CALL  leapfrog_matsuno_scheme(stage)
+          
           !      CASE ('leapfrog')
-          !      CALL leapfrog_scheme
+          !        CALL leapfrog_scheme
           !
           !      CASE ('adam_bashforth')
-          !      CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
-          !      CALL adam_bashforth_scheme
+          !        CALL dissip(f_u,f_du,f_ps,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+          !        CALL adam_bashforth_scheme
        CASE DEFAULT
           STOP
@@ -263,6 +270,4 @@
 
     IF(MOD(it+1,itau_adv)==0) THEN
-!       CALL transfert_request(f_wfluxt,req_i1) ! FIXME
-!      CALL transfert_request(f_hfluxt,req_e1) ! FIXME
 
        CALL advect_tracer(f_hfluxt,f_wfluxt,f_u, f_q,f_rhodz)  ! update q and rhodz after RK step
@@ -271,27 +276,26 @@
        ! FIXME : check that rhodz is consistent with ps
        IF (check) THEN
-         CALL transfert_request(f_rhodz,req_i1)
-         CALL transfert_request(f_ps,req_i1)
-         CALL transfert_request(f_dps,req_i1) ! FIXME
-         CALL transfert_request(f_wflux,req_i1) ! FIXME
          DO ind=1,ndomain
             CALL swap_dimensions(ind)
             CALL swap_geometry(ind)
-            rhodz=f_rhodz(ind); ps=f_ps(ind); dps=f_dps(ind); 
-            wflux=f_wflux(ind); wfluxt=f_wfluxt(ind)
+            rhodz=f_rhodz(ind); ps=f_ps(ind);
             CALL compute_rhodz(.FALSE., ps, rhodz)   
          END DO
        ENDIF
+
     END IF
+
+
+
 !----------------------------------------------------
-  jD_cur = jD_ref + day_ini - day_ref + it/day_step
-  jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
-  jD_cur = jD_cur + int(jH_cur)
-  jH_cur = jH_cur - int(jH_cur)
-!	print*,"Just b4 phys"
-    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
-!----------------------------------------------------
+!    jD_cur = jD_ref + day_ini - day_ref + it/day_step
+!    jH_cur = jH_ref + start_time + mod(it,day_step)/float(day_step)
+!    jD_cur = jD_cur + int(jH_cur)
+!    jH_cur = jH_cur - int(jH_cur)
+!    CALL physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
+
 !    CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u, f_q)
     ENDDO
+
  
   CONTAINS
@@ -310,10 +314,13 @@
          ps=f_ps(ind) ; dps=f_dps(ind) ; 
 
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             ps(ij)=ps(ij)+dt*dps(ij)
-           ENDDO
-         ENDDO
+         IF (omp_first) THEN
+           DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               ps(ij)=ps(ij)+dt*dps(ij)
+             ENDDO
+           ENDDO
+         ENDIF
+         
          hflux=f_hflux(ind);     hfluxt=f_hfluxt(ind)
          wflux=f_wflux(ind);     wfluxt=f_wfluxt(ind)
@@ -324,5 +331,5 @@
        du=f_du(ind) ; dtheta_rhodz=f_dtheta_rhodz(ind)
 
-       DO l=1,llm
+       DO l=ll_begin,ll_end
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -342,4 +349,5 @@
 
     SUBROUTINE RK_scheme(stage)
+      USE caldyn_gcm_mod
       IMPLICIT NONE
       INTEGER :: ind, stage
@@ -351,4 +359,38 @@
 
       tau = dt*coef(stage)
+
+      DO ind=1,ndomain
+         CALL swap_dimensions(ind)
+         CALL swap_geometry(ind)
+         ps=f_ps(ind)   
+         psm1=f_psm1(ind) 
+         dps=f_dps(ind) 
+         
+         IF (stage==1) THEN ! first stage : save model state in XXm1
+           IF (omp_first) THEN
+             DO j=jj_begin,jj_end
+               DO i=ii_begin,ii_end
+                 ij=(j-1)*iim+i
+                 psm1(ij)=ps(ij)
+               ENDDO
+             ENDDO
+           ENDIF
+         END IF
+         
+         ! updates are of the form x1 := x0 + tau*f(x1)
+           IF (omp_first) THEN
+             DO j=jj_begin,jj_end
+               DO i=ii_begin,ii_end
+                 ij=(j-1)*iim+i
+                 ps(ij)=psm1(ij)+tau*dps(ij)
+               ENDDO
+             ENDDO
+           ENDIF
+       
+       ENDDO
+
+       CALL send_message(f_ps,req_ps)      
+
+
       
       DO ind=1,ndomain
@@ -360,13 +402,6 @@
          
          IF (stage==1) THEN ! first stage : save model state in XXm1
-
-           DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-               ij=(j-1)*iim+i
-               psm1(ij)=ps(ij)
-             ENDDO
-           ENDDO
               
-           DO l=1,llm
+         DO l=ll_begin,ll_end
              DO j=jj_begin,jj_end
                DO i=ii_begin,ii_end
@@ -382,12 +417,6 @@
          END IF
          ! updates are of the form x1 := x0 + tau*f(x1)
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             ps(ij)=psm1(ij)+tau*dps(ij)
-           ENDDO
-         ENDDO
          
-         DO l=1,llm
+         DO l=ll_begin,ll_end
            DO j=jj_begin,jj_end
              DO i=ii_begin,ii_end
@@ -468,5 +497,4 @@
 
         
-!        IF (MOD(it,matsuno_period+1)==0) THEN
         IF (stage==1) THEN
           psm1(:)=ps(:) ; um1(:,:)=u(:,:) ; theta_rhodzm1(:,:)=theta_rhodz(:,:)
@@ -475,5 +503,4 @@
           theta_rhodz(:,:)=theta_rhodzm1(:,:)+tau*dtheta_rhodz(:,:)
 
-!        ELSE IF (MOD(it,matsuno_period+1)==1) THEN
         ELSE IF (stage==2) THEN
 
@@ -538,4 +565,5 @@
     USE icosa
     USE disvert_mod
+    USE omp_para
     LOGICAL, INTENT(IN) :: comp ! .TRUE. to compute, .FALSE. to check
     REAL(rstd), INTENT(IN) :: ps(iim*jjm)
@@ -544,16 +572,12 @@
     INTEGER :: l,i,j,ij,dd
     err=0.
-    IF(comp) THEN 
-       dd=1
-    ELSE
-!       dd=-1
-       dd=0
-    END IF
-
-    DO l = 1, llm
+
+    dd=0
+
+    DO l = ll_begin, ll_end
        DO j=jj_begin-dd,jj_end+dd
           DO i=ii_begin-dd,ii_end+dd
              ij=(j-1)*iim+i
-             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g
+             m = ( ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij) )/g 
              IF(comp) THEN
                 rhodz(ij,l) = m
@@ -570,5 +594,5 @@
           STOP
        ELSE
-          PRINT *, 'No discrepancy between ps and rhodz detected'
+!          PRINT *, 'No discrepancy between ps and rhodz detected'
        END IF
     END IF
@@ -578,4 +602,5 @@
   SUBROUTINE accumulate_fluxes(hflux,wflux, hfluxt,wfluxt, tau,fluxt_zero)
     USE icosa
+    USE omp_para
     REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm,llm), wflux(iim*jjm,llm+1)
     REAL(rstd), INTENT(INOUT) :: hfluxt(3*iim*jjm,llm), wfluxt(iim*jjm,llm+1)
@@ -585,7 +610,8 @@
 
     IF(fluxt_zero) THEN
-!       PRINT *, 'Accumulating fluxes (first)'
+
        fluxt_zero=.FALSE.
-       DO l=1,llm
+
+       DO l=ll_begin,ll_end
          DO j=jj_begin-1,jj_end+1
            DO i=ii_begin-1,ii_end+1
@@ -598,5 +624,5 @@
        ENDDO
 
-       DO l=1,llm+1
+       DO l=ll_begin,ll_endp1
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -607,6 +633,6 @@
        ENDDO
     ELSE
-!       PRINT *, 'Accumulating fluxes (next)'
-       DO l=1,llm
+
+       DO l=ll_begin,ll_end
          DO j=jj_begin-1,jj_end+1
            DO i=ii_begin-1,ii_end+1
@@ -619,5 +645,5 @@
        ENDDO
 
-       DO l=1,llm+1
+       DO l=ll_begin,ll_endp1
          DO j=jj_begin,jj_end
            DO i=ii_begin,ii_end
@@ -629,4 +655,5 @@
 
     END IF
+
   END SUBROUTINE accumulate_fluxes
   
Index: codes/icosagcm/trunk/src/trace.F90
===================================================================
--- codes/icosagcm/trunk/src/trace.F90	(revision 150)
+++ codes/icosagcm/trunk/src/trace.F90	(revision 151)
@@ -1,7 +1,21 @@
 MODULE trace
 
-
+  INTEGER,SAVE :: markId
+  
 CONTAINS
 
+  SUBROUTINE init_trace
+  IMPLICIT NONE
+#ifdef VTRACE
+#include <vt_user.inc>
+#endif 
+
+#ifdef VTRACE
+     VT_MARKER_DEF("marker", VT_MARKER_TYPE_HINT, markId)
+#endif
+  
+  END SUBROUTINE init_trace
+  
+  
   SUBROUTINE trace_start(name)
   IMPLICIT NONE
@@ -11,7 +25,9 @@
 #endif 
 
+!$OMP MASTER
 #ifdef VTRACE
      VT_USER_START(name)
 #endif
+!$OMP END MASTER
 
   END SUBROUTINE trace_start    
@@ -25,9 +41,26 @@
     CHARACTER(LEN=*),INTENT(IN) :: name
 
+!$OMP MASTER
 #ifdef VTRACE
      VT_USER_END(name)
 #endif
+!$OMP END MASTER
 
   END SUBROUTINE trace_end    
 
+  SUBROUTINE Marker(name)
+  IMPLICIT NONE
+  CHARACTER(LEN=*),INTENT(IN) :: name
+#ifdef VTRACE
+#include <vt_user.inc>
+#endif 
+
+!$OMP MASTER
+#ifdef VTRACE
+     VT_MARKER(markId,name)
+#endif
+!$OMP END MASTER
+
+  END SUBROUTINE Marker
+
 END MODULE trace
Index: codes/icosagcm/trunk/src/transfert.F90
===================================================================
--- codes/icosagcm/trunk/src/transfert.F90	(revision 150)
+++ codes/icosagcm/trunk/src/transfert.F90	(revision 151)
@@ -3,8 +3,12 @@
 #ifdef CPP_USING_MPI
   USE transfert_mpi_mod, ONLY : init_transfert, transfert_request=>transfert_request_mpi, req_i1,req_e1_vect, &
-                                req_e1_scal, req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field
+                                req_e1_scal, req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field, &
+                                t_message,init_message=>init_message_mpi,transfert_message=>transfert_message_mpi,  &
+                                send_message=>send_message_mpi,test_message=>test_message_mpi,wait_message=>wait_message_mpi,barrier
 #else 
-  USE transfert_mpi_mod, ONLY : init_transfert, transfert_request, req_i1,req_e1_vect, &
-                                req_e1_scal,request_add_point, create_request, gather_field
+  USE transfert_mpi_mod, ONLY : init_transfert, transfert_request=>transfert_request_seq, req_i1,req_e1_vect, &
+                                req_e1_scal,req_i0, req_e0_vect, req_e0_scal, request_add_point, create_request, gather_field, &
+                                t_message,init_message=>init_message_seq,transfert_message=>transfert_message_seq,  &
+                                send_message=>send_message_seq,test_message=>test_message_seq,wait_message=>wait_message_seq,barrier
 #endif
   
Index: codes/icosagcm/trunk/src/transfert_mpi.f90
===================================================================
--- codes/icosagcm/trunk/src/transfert_mpi.f90	(revision 150)
+++ codes/icosagcm/trunk/src/transfert_mpi.f90	(revision 151)
@@ -1,4 +1,5 @@
 MODULE transfert_mpi_mod
 USE genmod
+USE field_mod
   
   TYPE array
@@ -13,4 +14,10 @@
     REAL,POINTER    :: buffer_r4(:,:,:)
   END TYPE array
+  
+  TYPE t_buffer
+    REAL,POINTER    :: r2(:)
+    REAL,POINTER    :: r3(:,:)
+    REAL,POINTER    :: r4(:,:,:)
+  END TYPE t_buffer    
     
   TYPE t_request
@@ -41,7 +48,17 @@
   TYPE(t_request),POINTER :: req_e0_vect(:)
   
+  TYPE t_message
+    TYPE(t_request), POINTER :: request(:)
+    INTEGER :: nreq
+    INTEGER, POINTER :: mpi_req(:)
+    INTEGER, POINTER :: status(:,:)
+    TYPE(t_buffer),POINTER :: buffers(:) 
+    TYPE(t_field),POINTER :: field(:)
+    LOGICAL :: completed
+    LOGICAL :: pending
+  END TYPE t_message
   
 CONTAINS
-
+  
   SUBROUTINE init_transfert
   USE domain_mod
@@ -70,20 +87,4 @@
       DO j=jj_begin,jj_end+1
         CALL request_add_point(ind,ii_begin-1,j,req_i1)
-      ENDDO    
-    
-      DO i=ii_begin,ii_end
-!        CALL request_add_point(ind,i,jj_begin,req_i1)
-      ENDDO
-
-      DO j=jj_begin,jj_end
-!        CALL request_add_point(ind,ii_end,j,req_i1)
-      ENDDO    
-    
-      DO i=ii_begin,ii_end
-!        CALL request_add_point(ind,i,jj_end,req_i1)
-      ENDDO    
-
-      DO j=jj_begin,jj_end
-!        CALL request_add_point(ind,ii_begin,j,req_i1)
       ENDDO    
     
@@ -142,19 +143,4 @@
       ENDDO   
 
-      DO i=ii_begin+1,ii_end-1
-!        CALL request_add_point(ind,i,jj_begin,req_e1_scal,right)
-!        CALL request_add_point(ind,i,jj_end,req_e1_scal,right)
-      ENDDO
-    
-      DO j=jj_begin+1,jj_end-1
-!        CALL request_add_point(ind,ii_begin,j,req_e1_scal,rup)
-!        CALL request_add_point(ind,ii_end,j,req_e1_scal,rup)
-      ENDDO   
-
-!      CALL request_add_point(ind,ii_begin+1,jj_begin,req_e1_scal,left)
-!      CALL request_add_point(ind,ii_begin,jj_begin+1,req_e1_scal,ldown)
-!      CALL request_add_point(ind,ii_begin+1,jj_end,req_e1_scal,left)
-!      CALL request_add_point(ind,ii_end,jj_begin+1,req_e1_scal,ldown)
-
     ENDDO
 
@@ -211,18 +197,4 @@
       ENDDO   
 
-      DO i=ii_begin+1,ii_end-1
-!        CALL request_add_point(ind,i,jj_begin,req_e1_vect,right)
-!        CALL request_add_point(ind,i,jj_end,req_e1_vect,right)
-      ENDDO
-    
-      DO j=jj_begin+1,jj_end-1
-!        CALL request_add_point(ind,ii_begin,j,req_e1_vect,rup)
-!        CALL request_add_point(ind,ii_end,j,req_e1_vect,rup)
-      ENDDO   
-
-!      CALL request_add_point(ind,ii_begin+1,jj_begin,req_e1_vect,left)
-!      CALL request_add_point(ind,ii_begin,jj_begin+1,req_e1_vect,ldown)
-!      CALL request_add_point(ind,ii_begin+1,jj_end,req_e1_vect,left)
-!      CALL request_add_point(ind,ii_end,jj_begin+1,req_e1_vect,ldown)
   
     ENDDO  
@@ -322,5 +294,5 @@
     target_j=>req%target_j
     target_sign=>req%target_sign
-!    req%max_size=req%max_size*2
+
     ALLOCATE(req%src_domain(req%max_size*2))
     ALLOCATE(req%src_ind(req%max_size*2))
@@ -610,4 +582,532 @@
 
 
+  SUBROUTINE init_message_seq(field, request, message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE mpi_mod
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_request),POINTER :: request(:)
+    TYPE(t_message) :: message
+
+!$OMP MASTER    
+    message%request=>request
+!$OMP END MASTER    
+!$OMP BARRIER    
+
+  END SUBROUTINE init_message_seq
+
+  SUBROUTINE send_message_seq(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+
+    CALL transfert_request_seq(field,message%request)
+    
+  END SUBROUTINE send_message_seq
+  
+  SUBROUTINE test_message_seq(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+  END SUBROUTINE  test_message_seq
+  
+   
+  SUBROUTINE wait_message_seq(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+    
+  END SUBROUTINE wait_message_seq    
+
+  SUBROUTINE transfert_message_seq(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+
+   CALL send_message_seq(field,message)
+    
+  END SUBROUTINE transfert_message_seq    
+    
+    
+  SUBROUTINE init_message_mpi(field,request, message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE mpi_mod
+  IMPLICIT NONE
+  
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_request),POINTER :: request(:)
+    TYPE(t_message) :: message
+
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind
+    INTEGER :: dim3,dim4
+
+!$OMP MASTER
+    message%request=>request
+    nreq=sum(request(:)%nsend)+sum(request(:)%nrecv)
+    message%nreq=nreq
+    ALLOCATE(message%mpi_req(nreq))
+    ALLOCATE(message%buffers(nreq))
+    ALLOCATE(message%status(MPI_STATUS_SIZE,nreq))
+    
+    message%pending=.FALSE.
+    message%completed=.FALSE.
+    
+    IF (field(1)%data_type==type_real) THEN
+
+      IF (field(1)%ndim==2) THEN
+      
+        ireq=0
+        DO ind=1,ndomain
+          req=>request(ind)
+      
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r2,send%size)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r2,recv%size)
+          ENDDO
+        
+        ENDDO
+      
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+    
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval3d,2)
+          req=>request(ind)
+ 
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r3,send%size,dim3)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r3,recv%size,dim3)
+
+          ENDDO
+        
+        ENDDO
+
+
+      ELSE  IF (field(1)%ndim==4) THEN
+   
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval4d,2)
+          dim4=size(field(ind)%rval4d,3)
+          req=>request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r4,send%size,dim3,dim4)
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            CALL allocate_mpi_buffer(message%buffers(ireq)%r4,recv%size,dim3,dim4)
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+    ENDIF
+!$OMP END MASTER
+!$OMP BARRIER    
+  END SUBROUTINE init_message_mpi
+  
+  SUBROUTINE barrier
+  USE mpi_mod
+  USE mpipara
+  IMPLICIT NONE
+    
+    CALL MPI_BARRIER(comm_icosa,ierr)
+    
+  END SUBROUTINE barrier  
+    
+  SUBROUTINE transfert_message_mpi(field,message)
+  USE field_mod
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+    
+    CALL send_message_mpi(field,message)
+    CALL wait_message_mpi(message)
+    
+  END SUBROUTINE transfert_message_mpi
+
+
+  SUBROUTINE send_message_mpi(field,message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_field),POINTER :: field(:)
+    TYPE(t_message) :: message
+    REAL(rstd),POINTER :: rval2d(:) 
+    REAL(rstd),POINTER :: rval3d(:,:) 
+    REAL(rstd),POINTER :: rval4d(:,:,:) 
+    REAL(rstd),POINTER :: buffer_r2(:) 
+    REAL(rstd),POINTER :: buffer_r3(:,:) 
+    REAL(rstd),POINTER :: buffer_r4(:,:,:) 
+    INTEGER,POINTER :: value(:) 
+    INTEGER,POINTER :: sgn(:) 
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER, ALLOCATABLE :: mpi_req(:)
+    INTEGER, ALLOCATABLE :: status(:,:)
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind,n,l,m
+    INTEGER :: dim3,dim4
+
+!$OMP BARRIER
+
+    CALL trace_start("transfert_mpi")
+
+    nreq=message%nreq
+    message%field=>field
+
+!$OMP MASTER
+    message%completed=.FALSE.
+    message%pending=.TRUE.
+!$OMP END MASTER
+    
+    IF (field(1)%data_type==type_real) THEN
+      IF (field(1)%ndim==2) THEN
+
+        ireq=0
+        DO ind=1,ndomain
+          rval2d=>field(ind)%rval2d
+        
+          req=>message%request(ind)
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r2=>message%buffers(ireq)%r2
+            value=>send%value
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,send%size
+              buffer_r2(n)=rval2d(value(n))
+            ENDDO
+            
+            CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r2,send%size,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r2=>message%buffers(ireq)%r2
+!$OMP MASTER
+            CALL MPI_IRECV(buffer_r2,recv%size,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+        ENDDO
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+      
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval3d,2)
+          rval3d=>field(ind)%rval3d
+          req=>message%request(ind)
+ 
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r3=>message%buffers(ireq)%r3
+            value=>send%value
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+              DO n=1,send%size
+                buffer_r3(n,:)=rval3d(value(n),:)
+              ENDDO
+
+             CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r3,send%size*dim3,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r3=>message%buffers(ireq)%r3
+!$OMP MASTER           
+            CALL MPI_IRECV(buffer_r3,recv%size*dim3,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+
+          ENDDO
+        
+        ENDDO
+
+      ELSE  IF (field(1)%ndim==4) THEN
+    
+        ireq=0
+        DO ind=1,ndomain
+          dim3=size(field(ind)%rval4d,2)
+          dim4=size(field(ind)%rval4d,3)
+          rval4d=>field(ind)%rval4d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+            send=>req%send(isend)
+            buffer_r4=>message%buffers(ireq)%r4
+            value=>send%value
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,send%size
+               buffer_r4(n,:,:)=rval4d(value(n),:,:)
+            ENDDO
+
+           CALL trace_out
+
+!$OMP MASTER
+            CALL MPI_ISSEND(buffer_r4,send%size*dim3*dim4,MPI_REAL8,send%rank,ind,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r4=>message%buffers(ireq)%r4
+!$OMP MASTER           
+            CALL MPI_IRECV(buffer_r4,recv%size*dim3*dim4,MPI_REAL8,recv%rank,recv%domain,comm_icosa, message%mpi_req(ireq),ierr)
+!$OMP END MASTER
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+      
+    ENDIF
+
+    CALL trace_end("transfert_mpi")
+!$OMP BARRIER
+    
+  END SUBROUTINE send_message_mpi
+  
+  SUBROUTINE test_message_mpi(message)
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+    
+    INTEGER :: ierr
+!$OMP MASTER
+     IF (.NOT. message%pending) RETURN
+!$OMP END MASTER
+
+!$OMP MASTER
+     IF (.NOT. message%completed) CALL MPI_TESTALL(message%nreq,message%mpi_req,message%completed,message%status,ierr)
+!$OMP END MASTER
+  END SUBROUTINE  test_message_mpi
+  
+   
+  SUBROUTINE wait_message_mpi(message)
+  USE field_mod
+  USE domain_mod
+  USE mpi_mod
+  USE mpipara
+  USE omp_para
+  USE trace
+  IMPLICIT NONE
+    TYPE(t_message) :: message
+
+    TYPE(t_field),POINTER :: field(:)
+    REAL(rstd),POINTER :: rval2d(:) 
+    REAL(rstd),POINTER :: rval3d(:,:) 
+    REAL(rstd),POINTER :: rval4d(:,:,:) 
+    REAL(rstd),POINTER :: buffer_r2(:) 
+    REAL(rstd),POINTER :: buffer_r3(:,:) 
+    REAL(rstd),POINTER :: buffer_r4(:,:,:) 
+    INTEGER,POINTER :: value(:) 
+    INTEGER,POINTER :: sgn(:) 
+    TYPE(ARRAY),POINTER :: recv,send 
+    TYPE(t_request),POINTER :: req
+    INTEGER, ALLOCATABLE :: mpi_req(:)
+    INTEGER, ALLOCATABLE :: status(:,:)
+    INTEGER :: irecv,isend
+    INTEGER :: ireq,nreq
+    INTEGER :: ind,n,l,m
+    INTEGER :: dim3,dim4
+
+!$OMP BARRIER
+
+    CALL trace_start("transfert_mpi")
+
+    IF (.NOT. message%pending) RETURN
+    
+    field=>message%field
+    nreq=message%nreq
+    
+    IF (field(1)%data_type==type_real) THEN
+      IF (field(1)%ndim==2) THEN
+
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval2d=>field(ind)%rval2d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+     
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r2=>message%buffers(ireq)%r2
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval2d(value(n))=buffer_r2(n)*sgn(n)  
+            ENDDO        
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+      
+      
+      ELSE  IF (field(1)%ndim==3) THEN
+
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval3d=>field(ind)%rval3d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r3=>message%buffers(ireq)%r3
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+            
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval3d(value(n),:)=buffer_r3(n,:)*sgn(n)  
+            ENDDO  
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+
+      ELSE  IF (field(1)%ndim==4) THEN
+!$OMP MASTER
+        IF (.NOT. message%completed) CALL MPI_WAITALL(nreq,message%mpi_req,message%status,ierr)
+!$OMP END MASTER
+!$OMP BARRIER
+
+        ireq=0        
+        DO ind=1,ndomain
+          rval4d=>field(ind)%rval4d
+          req=>message%request(ind)
+
+          DO isend=1,req%nsend
+            ireq=ireq+1
+          ENDDO
+        
+          DO irecv=1,req%nrecv
+            ireq=ireq+1
+            recv=>req%recv(irecv)
+            buffer_r4=>message%buffers(ireq)%r4
+            value=>recv%value
+            sgn=>recv%sign
+
+            CALL trace_in
+
+!$OMP DO SCHEDULE(STATIC)
+            DO n=1,recv%size
+              rval4d(value(n),:,:)=buffer_r4(n,:,:)*sgn(n) 
+            ENDDO
+
+            CALL trace_out
+
+          ENDDO
+        
+        ENDDO
+      
+      ENDIF      
+      
+    ENDIF
+
+!$OMP MASTER
+    message%pending=.FALSE.
+!$OMP END MASTER
+
+    CALL trace_end("transfert_mpi")
+!$OMP BARRIER
+    
+  END SUBROUTINE wait_message_mpi
+
+
   SUBROUTINE transfert_request_mpi(field,request)
   USE field_mod
@@ -637,5 +1137,5 @@
 
     CALL trace_start("transfert_mpi")
-    
+
     IF (field(1)%data_type==type_real) THEN
       IF (field(1)%ndim==2) THEN
@@ -675,7 +1175,4 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -739,7 +1236,4 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -803,7 +1297,4 @@
         
         CALL MPI_WAITALL(nreq,mpi_req,status,ierr)
-!        DO ind=1,ndomain
-!          field(ind)%rval2d(:)=0.
-!        ENDDO
         
         DO ind=1,ndomain
@@ -837,5 +1328,5 @@
   END SUBROUTINE transfert_request_mpi
    
-  SUBROUTINE transfert_request(field,request)
+  SUBROUTINE transfert_request_seq(field,request)
   USE field_mod
   USE domain_mod
@@ -875,5 +1366,5 @@
     ENDDO
     
-  END SUBROUTINE transfert_request
+  END SUBROUTINE transfert_request_seq
   
   
@@ -949,5 +1440,19 @@
         
   END SUBROUTINE gather_field
-  
+
+   
+  SUBROUTINE trace_in
+  USE trace
+  IMPLICIT NONE
+  
+    CALL trace_start("transfert_buffer")
+  END SUBROUTINE trace_in              
+
+  SUBROUTINE trace_out
+  USE trace
+  IMPLICIT NONE
+  
+    CALL trace_end("transfert_buffer")
+  END SUBROUTINE trace_out              
 
 END MODULE transfert_mpi_mod
Index: codes/icosagcm/trunk/src/wind.f90
===================================================================
--- codes/icosagcm/trunk/src/wind.f90	(revision 150)
+++ codes/icosagcm/trunk/src/wind.f90	(revision 151)
@@ -3,4 +3,23 @@
 CONTAINS
 
+  SUBROUTINE un2ulonlat(f_u, f_ulon, f_ulat)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_u(:) ! IN  : normal velocity components on edges
+    TYPE(t_field), POINTER :: f_ulon(:), f_ulat(:) ! OUT : velocity reconstructed at hexagons
+    
+    REAL(rstd),POINTER :: u(:,:),  ulon(:,:), ulat(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       u=f_u(ind)
+       ulon=f_ulon(ind)
+       ulat=f_ulat(ind)
+       CALL compute_un2ulonlat(u,ulon, ulat)
+    END DO
+
+  END SUBROUTINE un2ulonlat
 
  
@@ -214,4 +233,18 @@
  END SUBROUTINE compute_wind_centered_lonlat_compound
 
+ SUBROUTINE compute_un2ulonlat(un, ulon, ulat)
+  USE icosa  
+    
+  IMPLICIT NONE
+  REAL(rstd),INTENT(IN)  :: un(3*iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: ulon(iim*jjm,llm)
+  REAL(rstd),INTENT(OUT) :: ulat(iim*jjm,llm)
+
+  REAL(rstd)             :: uc(iim*jjm,3,llm)
+    
+  CALL compute_wind_centered(un,uc)  
+  CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
+ 
+ END SUBROUTINE compute_un2ulonlat
 
 END MODULE wind_mod
Index: codes/icosagcm/trunk/src/write_field.f90
===================================================================
--- codes/icosagcm/trunk/src/write_field.f90	(revision 150)
+++ codes/icosagcm/trunk/src/write_field.f90	(revision 151)
@@ -75,5 +75,6 @@
       
       TYPE(t_field),POINTER :: field_glo(:)
-      
+
+!$OMP MASTER      
       IF(PRESENT(once)) THEN
          once_=once
@@ -101,4 +102,5 @@
       
       CALL deallocate_field_glo(field_glo)
+!$OMP END MASTER
       
    END SUBROUTINE writefield
@@ -333,6 +335,4 @@
     USE domain_mod
     USE field_mod
-!    USE dimensions
-!    USE geometry
     IMPLICIT NONE
       CHARACTER(LEN=*),INTENT(IN) :: name_in
@@ -482,8 +482,4 @@
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
-!        ncell=ncell+n
-
-!     ENDDO
      
      ELSE IF (Field(ind_b)%field_type==field_Z) THEN
@@ -604,9 +600,4 @@
           ENDDO
         ENDIF
-        
-!        PRINT *,NF90_STRERROR(status)
-!        ncell=ncell+n
-!
-!     ENDDO
      
      ENDIF
@@ -712,9 +703,8 @@
             ENDDO
           ENDDO
-!           DO l=1,size(field(ind)%rval3d,2)
-            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
+          status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,  &
                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-!           ENDDO
           DEALLOCATE(field_val3d)
+
         ELSE IF (field(1)%ndim==4) THEN
 
@@ -732,13 +722,11 @@
               ENDDO
              ENDDO
-!            DO l=1,size(field(ind)%rval4d,2)
+
             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d(:,l), &
                                 start=(/ displ+ncell,l,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO         
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
         ncell=ncell+n
 
@@ -810,9 +798,9 @@
             ENDDO
           ENDDO
-!           DO l=1,size(field(ind)%rval3d,2)
+
            status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(1),Field_val3d,    &
                                start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval3d,2),1 /))
-!           ENDDO
           DEALLOCATE(field_val3d)
+
         ELSE IF (field(1)%ndim==4) THEN
 
@@ -836,14 +824,10 @@
             ENDDO
 
-!            DO l=1,size(field(ind)%rval4d,2)
-
             status=NF90_PUT_VAR(FieldId(Index),FieldVarId(index)%nc_id(q),Field_val3d,  &
                                 start=(/ displ+ncell,1,FieldIndex(Index) /), count=(/n,size(field(ind)%rval4d,2),1 /))
-!            ENDDO
             DEALLOCATE(field_val3d)
           ENDDO
         ENDIF
         
-!        PRINT *,NF90_STRERROR(status)
         ncell=ncell+n
 
@@ -1159,6 +1143,4 @@
     USE metric
     USE spherical_geom_mod
-!    USE dimensions
-!    USE geometry
     IMPLICIT NONE
       CHARACTER(LEN=*),INTENT(IN) :: name_in
@@ -1292,16 +1274,5 @@
       
         status = NF90_ENDDEF(ncid)      
-
-!        ncell=1
-!        DO ind=ind_b,ind_e
-!          d=>domain_type(ind)
-  
-!          n=0
-!          DO j=d%jj_begin-halo_size,d%jj_end+halo_size
-!            DO i=d%ii_begin-halo_size,d%ii_end+halo_size
-!              IF (single .OR. d%assign_domain(i,j)==ind) n=n+1
-!            ENDDO
-!          ENDDO
-          
+         
          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
           
@@ -1330,7 +1301,5 @@
          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
   
-!         ncell=ncell+n
          DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!      END DO 
 
     ELSE IF (Field(ind_b)%field_type==field_Z) THEN
@@ -1416,22 +1385,4 @@
         status = NF90_ENDDEF(ncid)      
 
-!        ncell=1
-!        DO ind=ind_b,ind_e
-!          d=>domain_type(ind)
-!  
-!          n=0
-!          DO j=d%jj_begin+1,d%jj_end
-!            DO i=d%ii_begin,d%ii_end-1
-!              n=n+1
-!            ENDDO
-!          ENDDO
-!
-!          DO j=d%jj_begin,d%jj_end-1
-!            DO i=d%ii_begin+1,d%ii_end
-!              n=n+1
-!            ENDDO
-!          ENDDO
-
-!         ALLOCATE(lon(n),lat(n),bounds_lon(0:nvert-1,n),bounds_lat(0:nvert-1,n))
          ALLOCATE(lon(ncell),lat(ncell),bounds_lon(0:nvert-1,ncell),bounds_lat(0:nvert-1,ncell))
           
@@ -1446,9 +1397,5 @@
                CALL xyz2lonlat(d%vertex(:,vdown,i,j),lon(n),lat(n))
                lon(n)=lon(n)*180/Pi
- !              IF (lon(n)<0) lon(n)=lon(n)+360
                lat(n)=lat(n)*180/Pi
-!               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!               CALL xyz2lonlat(xyz_i(nij+t_ldown,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!               CALL xyz2lonlat(xyz_i(nij+t_rdown,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
 
                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
@@ -1459,5 +1406,4 @@
                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
                ENDDO
              ENDDO
@@ -1468,12 +1414,7 @@
                nij=(j-1)*d%iim+i
                n=n+1
-!              CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
                CALL xyz2lonlat(d%vertex(:,vup,i,j),lon(n),lat(n))
                lon(n)=lon(n)*180/Pi
-!              IF (lon(n)<0) lon(n)=lon(n)+360
                lat(n)=lat(n)*180/Pi
-!              CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
-!              CALL xyz2lonlat(xyz_i(nij+t_rup,:)/radius,bounds_lon(1,n), bounds_lat(1,n))
-!              CALL xyz2lonlat(xyz_i(nij+t_lup,:)/radius,bounds_lon(2,n), bounds_lat(2,n))
                CALL xyz2lonlat(d%xyz(:,i,j),bounds_lon(0,n), bounds_lat(0,n))
                CALL xyz2lonlat(d%xyz(:,i,j+1),bounds_lon(1,n), bounds_lat(1,n))
@@ -1483,5 +1424,4 @@
                  bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                  bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                 IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
                ENDDO
              ENDDO
@@ -1494,12 +1434,7 @@
          status=NF90_PUT_VAR(ncid,bounds_latId,REAL(bounds_lat,r8),start=(/ 1,1 /),count=(/ nvert,ncell /))
   
-!          ncell=ncell+n
           DEALLOCATE(lon,lat,bounds_lon,bounds_lat)
-!      END DO          
-    ENDIF
-
-
-      
-!      status = NF90_CLOSE(ncid)
+      ENDIF
+
 
     END SUBROUTINE Create_Header_gen
@@ -1783,5 +1718,4 @@
               CALL xyz2lonlat(xyz_v(nij+z_down,:)/radius,lon(n),lat(n))
               lon(n)=lon(n)*180/Pi
- !             IF (lon(n)<0) lon(n)=lon(n)+360
               lat(n)=lat(n)*180/Pi
               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
@@ -1792,5 +1726,4 @@
                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
               ENDDO
             ENDDO
@@ -1803,5 +1736,4 @@
               CALL xyz2lonlat(xyz_v(nij+z_up,:)/radius,lon(n),lat(n))
               lon(n)=lon(n)*180/Pi
-!              IF (lon(n)<0) lon(n)=lon(n)+360
               lat(n)=lat(n)*180/Pi
               CALL xyz2lonlat(xyz_i(nij,:)/radius,bounds_lon(0,n), bounds_lat(0,n))
@@ -1812,5 +1744,4 @@
                 bounds_lat(k,n)=bounds_lat(k,n)*180/Pi
                 bounds_lon(k,n)=bounds_lon(k,n)*180/Pi
-!                IF (bounds_lon(k,n)<0) bounds_lon(k,n)=bounds_lon(k,n)+360
               ENDDO
             ENDDO
@@ -1829,8 +1760,5 @@
 
 
-      
-!      status = NF90_CLOSE(ncid)
-
-    end subroutine Create_Header_mpi  
+   END SUBROUTINE Create_Header_mpi  
    
    SUBROUTINE Close_files
