Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 294)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 295)
@@ -26,4 +26,5 @@
   SUBROUTINE init_advect_tracer
     USE advect_mod
+    USE omp_para
     REAL(rstd),POINTER :: tangent(:,:)
     REAL(rstd),POINTER :: normal(:,:)
@@ -48,5 +49,5 @@
        tangent=f_tangent(ind)
        sqrt_leng=f_sqrt_leng(ind)
-       CALL init_advect(normal,tangent,sqrt_leng)
+       IF (is_omp_level_master) CALL init_advect(normal,tangent,sqrt_leng)
     END DO
 
@@ -238,5 +239,5 @@
 
 !--> flush dzqw, adzqw
-!!$OMP BARRIER
+!$OMP BARRIER
 
     ! minmod-limited slope of q
@@ -258,5 +259,5 @@
 
     ! 0 slope in top and bottom layers
-    IF (omp_first) THEN
+    IF (is_omp_first_level) THEN
       DO ij=ijb,ije
            dzq(ij,1)=0.
@@ -264,5 +265,5 @@
     ENDIF
       
-    IF (omp_last) THEN
+    IF (is_omp_last_level) THEN
       DO ij=ijb,ije
           dzq(ij,llm)=0.
@@ -271,5 +272,5 @@
 
 !---> flush dzq
-!!$OMP BARRIER  
+!$OMP BARRIER  
 
     ! sigw = fraction of mass that leaves level l/l+1
@@ -290,5 +291,5 @@
     END DO
     ! wq = 0 at top and bottom
-    IF (omp_first) THEN
+    IF (is_omp_first_level) THEN
        DO ij=ijb,ije
             wq(ij,1)=0.
@@ -296,5 +297,5 @@
     ENDIF
     
-    IF (omp_last) THEN
+    IF (is_omp_last_level) THEN
       DO ij=ijb,ije
             wq(ij,llm+1)=0.
@@ -303,5 +304,5 @@
 
 ! --> flush wq
-!!$OMP BARRIER
+!$OMP BARRIER
 
 
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 294)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 295)
@@ -33,4 +33,5 @@
     USE exner_mod
     USE mpipara
+    USE omp_para
     IMPLICIT NONE
     CHARACTER(len=255) :: def
@@ -50,5 +51,5 @@
        STOP
     END SELECT
-    IF (is_mpi_root) PRINT *, 'caldyn_conserv=',def
+    IF (is_master) PRINT *, 'caldyn_conserv=',def
 
     CALL allocate_caldyn
@@ -72,5 +73,5 @@
     CALL allocate_field(f_qv,field_z,type_real,llm) 
   
-    CALL allocate_field(f_buf_i,   field_t,type_real,llm)
+    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
     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
@@ -102,5 +103,5 @@
     INTEGER :: ind,i,j,ij,l
 
-    IF (omp_first) THEN
+    IF (is_omp_first_level) THEN
        DO ind=1,ndomain
           IF (.NOT. assigned_domain(ind)) CYCLE
@@ -128,5 +129,5 @@
     ENDIF
 
-!    !$OMP BARRIER
+    !$OMP BARRIER
   END SUBROUTINE caldyn_BC
    
@@ -143,4 +144,5 @@
     USE omp_para
     USE output_field_mod
+    USE checksum_mod
     IMPLICIT NONE
     LOGICAL,INTENT(IN)    :: write_out
@@ -291,8 +293,8 @@
     END SELECT
 
-!!$OMP BARRIER
+!$OMP BARRIER
     IF (write_out) THEN
 
-       IF (is_mpi_root) PRINT *,'CALL write_output_fields'
+       IF (is_master) PRINT *,'CALL write_output_fields'
 
 ! ---> for openMP test to fix later
@@ -450,17 +452,18 @@
 
 !!! Compute exner function and geopotential
-       DO l = 1,llm
+       IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
+         DO l = 1,llm
 !         !$OMP DO SCHEDULE(STATIC) 
           !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext         
-             p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-             !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
-             exner_ik = cpp * (p_ik/preff) ** kappa
-             pk(ij,l) = exner_ik
+            DO ij=ij_begin_ext,ij_end_ext         
+               p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
+               !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
+               exner_ik = cpp * (p_ik/preff) ** kappa
+               pk(ij,l) = exner_ik
              ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-             geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
+               geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
           ENDDO
-       ENDDO
-
+         ENDDO
+       ENDIF
     ELSE 
        ! We are using a Lagrangian vertical coordinate
@@ -471,5 +474,6 @@
        IF(boussinesq) THEN ! compute only geopotential : pressure pk will be computed in compute_caldyn_horiz
           ! specific volume 1 = dphi/g/rhodz
-          DO l = 1,llm
+         IF (is_omp_level_master) THEN ! no openMP on vertical due to dependency
+           DO l = 1,llm
 !             !$OMP DO SCHEDULE(STATIC) 
              !DIR$ SIMD
@@ -477,40 +481,47 @@
                 geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
              ENDDO
-          ENDDO
+           ENDDO
+         ENDIF
        ELSE ! non-Boussinesq, compute geopotential and Exner pressure
           ! uppermost layer
-          !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext         
-             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-          END DO
-          ! other layers
-          DO l = llm-1, 1, -1
-
-!          !$OMP DO SCHEDULE(STATIC) 
-             !DIR$ SIMD
-             DO ij=ij_begin_ext,ij_end_ext         
-                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
-             END DO
-          END DO
+         IF (is_omp_level_master) THEN  ! no openMP on vertical due to dependency
+
+           !DIR$ SIMD
+           DO ij=ij_begin_ext,ij_end_ext         
+              pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
+           END DO
+           ! other layers
+           DO l = llm-1, 1, -1
+
+!           !$OMP DO SCHEDULE(STATIC) 
+              !DIR$ SIMD
+              DO ij=ij_begin_ext,ij_end_ext         
+                 pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
+              END DO
+           END DO
           ! surface pressure (for diagnostics)
-          DO ij=ij_begin_ext,ij_end_ext         
-             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-          END DO
+           DO ij=ij_begin_ext,ij_end_ext         
+              ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
+           END DO
 
           ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-          DO l = 1,llm
+           DO l = 1,llm
 
 !             !$OMP DO SCHEDULE(STATIC) 
              !DIR$ SIMD
-             DO ij=ij_begin_ext,ij_end_ext         
-                p_ik = pk(ij,l)
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
-                pk(ij,l) = exner_ik
-             ENDDO
-          ENDDO
+              DO ij=ij_begin_ext,ij_end_ext         
+                 p_ik = pk(ij,l)
+                 exner_ik = cpp * (p_ik/preff) ** kappa
+                 geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+                 pk(ij,l) = exner_ik
+              ENDDO
+            ENDDO
+          ENDIF
        END IF
 
     END IF
+
+!ym flush geopot
+!$OMP BARRIER
 
   CALL trace_end("compute_geopot")
@@ -799,21 +810,24 @@
   CALL trace_start("compute_caldyn_vert")
 
-!!$OMP BARRIER   
+!$OMP BARRIER   
 !!! cumulate mass flux convergence from top to bottom
-  DO  l = llm-1, 1, -1
+  IF (is_omp_level_master) THEN
+    DO  l = llm-1, 1, -1
 !    IF (caldyn_conserv==energy) CALL test_message(req_qu) 
 
 !!$OMP DO SCHEDULE(STATIC) 
 !DIR$ SIMD
-    DO ij=ij_begin,ij_end         
-        convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+      DO ij=ij_begin,ij_end         
+          convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+      ENDDO
     ENDDO
-  ENDDO
-
-! IMPLICIT FLUSH on convm
+  ENDIF
+
+!$OMP BARRIER
+! FLUSH on convm
 !!!!!!!!!!!!!!!!!!!!!!!!!
 
 ! compute dps
-  IF (omp_first) THEN
+  IF (is_omp_first_level) THEN
 !DIR$ SIMD
      DO ij=ij_begin,ij_end         
@@ -834,4 +848,7 @@
   ENDDO
 
+ !--> flush wflux  
+ !$OMP BARRIER
+
   DO l=ll_begin,ll_endm1
 !DIR$ SIMD
@@ -847,4 +864,5 @@
     ENDDO
   ENDDO
+
   
 ! Compute vertical transport
@@ -859,5 +877,5 @@
 
  !--> flush wwuu  
- ! !$OMP BARRIER
+ !$OMP BARRIER
 
 ! Add vertical transport to du
@@ -1018,8 +1036,11 @@
        ps = f_ps(ind)
        p  = f_p(ind)
+!$OMP BARRIER
        CALL compute_pression(ps,p,0)
        pk = f_pk(ind)
        pks = f_pks(ind)
+!$OMP BARRIER
        CALL compute_exner(ps,p,pks,pk,0)
+!$OMP BARRIER
        theta_rhodz = f_theta_rhodz(ind)
        theta = f_theta(ind)
Index: codes/icosagcm/trunk/src/check_conserve.f90
===================================================================
--- codes/icosagcm/trunk/src/check_conserve.f90	(revision 294)
+++ codes/icosagcm/trunk/src/check_conserve.f90	(revision 295)
@@ -9,10 +9,6 @@
    
   PUBLIC init_check_conserve, check_conserve 
-    REAL(rstd),SAVE:: mtot0,ztot0,etot0,ang0,stot0,rmsv0
-!$OMP THREADPRIVATE(mtot0,ztot0,etot0,ang0,stot0,rmsv0)
-    REAL(rstd),SAVE:: etot,ang,stot,rmsv
-!$OMP THREADPRIVATE(etot,ang,stot,rmsv)
-    REAL(rstd),SAVE:: ztot
-!$OMP THREADPRIVATE(ztot)	
+    REAL(rstd),SAVE:: mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0
+!$OMP THREADPRIVATE(mtot0,ztot0,etot0,angtot0,stot0,rmsvtot0)
      
 CONTAINS 
@@ -36,5 +32,6 @@
   USE caldyn_gcm_mod
   USE exner_mod 
-  USE mpipara, ONLY : is_mpi_root, comm_icosa
+  USE mpipara, ONLY : is_mpi_master, comm_icosa
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -48,8 +45,7 @@
     INTEGER::ind,ierr
     REAL(rstd) :: mtot, rmsdpdt
-
-    etot=0.0; ang=0.0;stot=0.0;rmsv=0.0 
-    ztot = 0.0 
-
+    REAL(rstd) :: etot, stot, angtot, rmsvtot, ztot
+
+    CALL transfert_request(f_ue,req_e1_vect)
     CALL pression(f_ps,f_p)
 
@@ -65,34 +61,33 @@
     CALL vorticity(f_ue,f_vort)
     CALL check_mass_conserve(f_ps,f_dps,mtot,rmsdpdt)
-    CALL check_PV 
+    CALL check_PV(ztot) 
     CALL exner(f_ps,f_p,f_pks,f_pk)
-    CALL check_EN(f_ue,f_theta_rhodz,f_phis) 
-
-    IF (is_mpi_root) THEN 
-!$OMP MASTER       
-       IF ( it == itau0  ) Then 
+    CALL check_EN(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot) 
+
+    IF (is_master) THEN 
+       IF ( it == itau0  ) THEN 
           ztot0 = ztot
           mtot0 = mtot
           etot0 = etot 
-          ang0  = ang 
+          angtot0  = angtot 
           stot0 = stot 
        END IF
 
-       rmsv=SQRT(rmsv/mtot) 
+       rmsvtot=SQRT(rmsvtot/mtot) 
        ztot=ztot/ztot0-1. ; mtot=mtot/mtot0-1. 
-       etot=etot/etot0-1. ; ang=ang/ang0-1. ; stot=stot/stot0-1. 
+       etot=etot/etot0-1. ; angtot=angtot/angtot0-1. ; stot=stot/stot0-1. 
        rmsdpdt= daysec*1.e-2*sqrt(rmsdpdt/ncell_glo)  
        
        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,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
+       WRITE(134,*)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot 
        WRITE(134,*)"=================================================="
-       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsv,ang
+       WRITE(*,4000)mtot,rmsdpdt,etot,ztot,stot,rmsvtot,angtot
        
 4000   FORMAT(10x,'masse',5x,'rmsdpdt',5x,'energie',5x,'enstrophie'  &
             ,5x,'entropie',5x,'rmsv',5x,'mt.ang',/,'GLOB  '                &
-            ,e10.3,e13.6,5e10.3/)     
+            ,e10.3,e13.6,5e13.3/)     
        close(134)
-!$OMP END MASTER 
+      
     END IF
   END SUBROUTINE check_conserve
@@ -105,4 +100,5 @@
   USE transfert_omp_mod
   USE icosa
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -117,5 +113,5 @@
     mloc=0.0; rmsloc=0.0
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -141,5 +137,6 @@
 !$OMP END MASTER
     ELSE
-      
+      mtot=mloc
+      rmsdpdt=rmsloc
     ENDIF
     
@@ -148,12 +145,16 @@
 !---------------------------------------------------------------------
 
-  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis) 
+  SUBROUTINE check_en(f_ue,f_theta_rhodz,f_phis, etot, stot, angtot, rmsvtot)  
   USE icosa
   USE pression_mod
   USE vorticity_mod
+  USE mpi_mod
+  USE mpipara
+  USE transfert_omp_mod
   IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ue(:)
     TYPE(t_field), POINTER :: f_theta_rhodz(:)
     TYPE(t_field), POINTER :: f_phis(:)
+    REAL(rstd),INTENT(OUT) :: etot, stot, angtot, rmsvtot
     REAL(rstd), POINTER :: ue(:,:)
     REAL(rstd), POINTER :: p(:,:) 
@@ -162,6 +163,9 @@
     REAL(rstd), POINTER :: phis(:) 
     REAL(rstd), POINTER :: rhodz(:,:) 
+    REAL(rstd) :: e, s, ang, rmsv
+    REAL(rstd) :: e_mpi, s_mpi, ang_mpi, rmsv_mpi
     INTEGER :: ind
 
+    e = 0 ; s = 0 ; ang = 0 ; rmsv = 0
 
     DO ind=1,ndomain 
@@ -175,13 +179,32 @@
       pk=f_pk(ind) 
       phis=f_phis(ind) 
-      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis)
+      CALL compute_en(ind,ue,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
     END DO 
- 
+
+    IF (using_mpi) THEN
+      CALL reduce_sum_omp(e, e_mpi)
+      CALL reduce_sum_omp(s, s_mpi)
+      CALL reduce_sum_omp(ang, ang_mpi)
+      CALL reduce_sum_omp(rmsv, rmsv_mpi)
+!$OMP MASTER    
+      CALL MPI_REDUCE(e_mpi, etot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+      CALL MPI_REDUCE(s_mpi, stot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+      CALL MPI_REDUCE(ang_mpi, angtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+      CALL MPI_REDUCE(rmsv_mpi, rmsvtot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+!$OMP END MASTER
+    ELSE
+      etot=e
+      stot=s
+      angtot=ang
+      rmsvtot=rmsv
+    ENDIF 
+
   END SUBROUTINE check_en 
 
-  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis)
+  SUBROUTINE compute_en(ind,u,p,rhodz,theta_rhodz,pk,phis, e, s, ang, rmsv)
   USE icosa
   USE disvert_mod
   USE wind_mod
+  USE omp_para
   IMPLICIT NONE
     INTEGER,INTENT(IN)::ind 
@@ -192,4 +215,9 @@
     REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
+    REAL(rstd),INTENT(INOUT) :: e
+    REAL(rstd),INTENT(INOUT) :: s
+    REAL(rstd),INTENT(INOUT) :: ang
+    REAL(rstd),INTENT(INOUT) :: rmsv
+    
     REAL(rstd):: theta(iim*jjm,llm)
     REAL(rstd)::KE(iim*jjm,llm) 
@@ -204,5 +232,5 @@
 
   
-    DO l = 1, llm
+    DO l = ll_begin, ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -211,5 +239,5 @@
           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)
+	    s = s + Ai(ij)*theta_rhodz(ij,l)
           END IF 
         END DO 
@@ -217,5 +245,5 @@
     END DO 
 
-    DO l = 1, llm
+    DO l = ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -234,10 +262,10 @@
     END DO  
 
-    DO l = 1, llm
+    DO l = ll_begin,ll_end
       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))
+	    e = e + masse(ij,l)*(phis(ij)+theta(ij,l)*pk(ij,l)+KE(ij,l))
           END IF 
         END DO 
@@ -251,5 +279,5 @@
     radomeg = rad*omega
 
-    DO l = 1, llm
+    DO l = ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -267,11 +295,17 @@
 !---------------------------------------------------------------------
 
-  SUBROUTINE check_PV 
-  USE icosa
-  IMPLICIT NONE
+  SUBROUTINE check_PV(ztot)
+  USE icosa
+  USE mpi_mod
+  USE mpipara
+  USE transfert_omp_mod
+  IMPLICIT NONE
+     REAL(rstd),INTENT(OUT) :: ztot
      REAL(rstd), POINTER :: vort(:,:)
      REAL(rstd), POINTER :: rhodz(:,:) 
      INTEGER :: ind
-
+     REAL(rstd) :: z,z_mpi
+
+     z=0
      DO ind=1,ndomain 
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -280,15 +314,26 @@
        vort=f_vort(ind)
        rhodz=f_rhodz(ind) 
-       CALL compute_PV(vort,rhodz) 
+       CALL compute_PV(vort,rhodz,z) 
      ENDDO
+
+    IF (using_mpi) THEN
+      CALL reduce_sum_omp(z, z_mpi)
+!$OMP MASTER    
+      CALL MPI_REDUCE(z_mpi, ztot, 1, MPI_REAL8, MPI_SUM, 0, comm_icosa, ierr)
+!$OMP END MASTER
+    ELSE
+      ztot=z
+    ENDIF 
  
   END SUBROUTINE check_PV 
 
-  SUBROUTINE compute_PV(vort,rhodz) 
+  SUBROUTINE compute_PV(vort,rhodz,z) 
   USE icosa
   USE disvert_mod
+  USE omp_para
   IMPLICIT NONE 
     REAL(rstd),INTENT(IN) :: vort(2*iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 
+    REAL(rstd),INTENT(INOUT) :: z
     REAL(rstd)::qv1,qv2 
     REAL(rstd)::hv1,hv2  
@@ -297,5 +342,5 @@
     hv1 = 0.0 ; hv2 = 0.0 
    
-    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
@@ -313,8 +358,9 @@
           qv2 =( vort(ij+z_down,l)+fv(ij+z_down) )/hv2
 
-          ztot = ztot + (qv1*qv1*hv1 + qv2*qv2*hv2) 
+          z = z + (qv1*qv1*hv1 + qv2*qv2*hv2) 
 
         ENDDO
       ENDDO
+      
     ENDDO
 
Index: codes/icosagcm/trunk/src/checksum.f90
===================================================================
--- codes/icosagcm/trunk/src/checksum.f90	(revision 295)
+++ codes/icosagcm/trunk/src/checksum.f90	(revision 295)
@@ -0,0 +1,66 @@
+MODULE checksum_mod
+
+
+! simple function to perform checksum on T field
+! works only on 1 process
+! usefull to check openMP synchronisation problem
+! need to be enhanced...
+
+CONTAINS
+
+  SUBROUTINE checksum(field)
+  USE icosa
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: field(:)
+    INTEGER :: intval(2)
+    INTEGER :: ind,i,j,ij,l
+    INTEGER :: tot_sum
+        
+
+!$OMP BARRIER
+!$OMP MASTER
+    tot_sum=0
+    DO ind=1,ndomain
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+       
+      IF (field(ind)%ndim==2) THEN
+        DO j=jj_begin,jj_end
+           DO i=ii_begin,ii_end
+              ij=(j-1)*iim+i
+              IF (domain(ind)%own(i,j)) THEN
+                intval=transfer(field(ind)%rval2d(ij),intval,2)  
+                tot_sum=tot_sum+intval(1)+intval(2)
+              ENDIF
+           ENDDO
+        ENDDO
+   
+      ELSE IF (field(ind)%ndim==3) THEN
+        
+        DO l=1,size(field(ind)%rval3d,2)
+          DO j=jj_begin,jj_end
+             DO i=ii_begin,ii_end
+                ij=(j-1)*iim+i
+                IF (domain(ind)%own(i,j)) THEN
+                  intval=transfer(field(ind)%rval3d(ij,l),intval,2)  
+                  tot_sum=tot_sum+intval(1)+intval(2)
+                ENDIF
+             ENDDO
+          ENDDO
+        ENDDO
+     
+      ENDIF
+    ENDDO
+!$OMP END MASTER
+!$OMP BARRIER    
+    
+!$OMP MASTER
+    PRINT*,"CheckSum Field",field(1)%name,tot_sum
+!$OMP END MASTER
+  
+  
+  
+  END SUBROUTINE checksum
+  
+END MODULE checksum_mod
+  
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 294)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 295)
@@ -71,4 +71,5 @@
   USE time_mod
   USE transfert_omp_mod
+  USE omp_para
   IMPLICIT NONE
   
@@ -100,15 +101,15 @@
    CASE('none')
       rayleigh_friction_type=0
-      IF (is_mpi_root) PRINT *, 'No Rayleigh friction'
+      IF (is_master) PRINT *, 'No Rayleigh friction'
    CASE('dcmip2_schaer_noshear')
       rayleigh_friction_type=1
       rayleigh_shear=0
-      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
+      IF (is_master) PRINT *, 'Rayleigh friction : Schaer-like mountain without shear DCMIP2.1'
    CASE('dcmip2_schaer_shear')
       rayleigh_shear=1
       rayleigh_friction_type=2
-      IF (is_mpi_root) PRINT *, 'Rayleigh friction : Schaer-like mountain with shear DCMIP2.2'
+      IF (is_master) 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_master) PRINT *, 'Bad selector : rayleigh_friction_type =', TRIM(rayleigh_friction_key), ' in dissip_gcm.f90/init_dissip'
       STOP
    END SELECT
@@ -119,5 +120,5 @@
       rayleigh_tau = rayleigh_tau / scale_factor
       IF(rayleigh_tau<=0) THEN
-         IF (is_mpi_root) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
+         IF (is_master) PRINT *, 'Forbidden : negative value for rayleigh_friction_tau =',rayleigh_tau
          STOP
       END IF
@@ -156,7 +157,8 @@
     cdivgrad=-1
     cgradrot=-1
-    
+
+!$OMP BARRIER
+!$OMP MASTER
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -175,6 +177,6 @@
       ENDDO        
     ENDDO
-  
-
+!$OMP END MASTER  
+!$OMP BARRIER
 
     DO it=1,20
@@ -184,5 +186,5 @@
         CALL transfert_request(f_u,req_e1_vect)
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
@@ -197,5 +199,5 @@
       
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -214,5 +216,5 @@
 
       IF (using_mpi) THEN 
-        CALL reduce_sum_omp(dumax,dumax1)
+        CALL reduce_max_omp(dumax,dumax1)
 !$OMP MASTER        
         CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
@@ -220,10 +222,10 @@
         CALL bcast_omp(dumax)  
       ELSE
-        CALL allreduce_sum_omp(dumax,dumax1)
+        CALL allreduce_max_omp(dumax,dumax1)
         dumax=dumax1
       ENDIF  
 			
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -232,18 +234,19 @@
         u=du/dumax
       ENDDO
-      IF (is_mpi_root) PRINT *,"gradiv : it :",it ,": dumax",dumax
+      IF (is_master) PRINT *,"gradiv : it :",it ,": dumax",dumax
 
     ENDDO 
-    IF (is_mpi_root) PRINT *,"gradiv : dumax",dumax
-    IF (is_mpi_root) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
+    IF (is_master) PRINT *,"gradiv : dumax",dumax
+    IF (is_master) PRINT *, 'mean T-cell edge size (km)', 1.45*radius/iim_glo/1000., &
                               'effective T-cell half-edge size (km)', dumax**(-.5/nitergdiv)/1000
-    IF (is_mpi_root)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', &
+    IF (is_master)  PRINT *, 'Max. time step assuming c=340 m/s and Courant number=2.8 :', &
                                2.8/340.*dumax**(-.5/nitergdiv)
 
     cgraddiv=dumax**(-1./nitergdiv)
-    IF (is_mpi_root) PRINT *,"cgraddiv : ",cgraddiv
-
+    IF (is_master) PRINT *,"cgraddiv : ",cgraddiv
+
+!$OMP BARRIER
+!$OMP MASTER
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -262,4 +265,6 @@
       ENDDO        
     ENDDO
+!$OMP END MASTER  
+!$OMP BARRIER
 
 
@@ -270,5 +275,5 @@
         CALL transfert_request(f_u,req_e1_vect)
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
@@ -283,5 +288,5 @@
       
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -300,5 +305,5 @@
 
       IF (using_mpi) THEN 
-        CALL reduce_sum_omp(dumax,dumax1)
+        CALL reduce_max_omp(dumax,dumax1)
 !$OMP MASTER        
         CALL MPI_ALLREDUCE(dumax1,dumax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
@@ -306,5 +311,5 @@
         CALL bcast_omp(dumax)  
       ELSE
-        CALL allreduce_sum_omp(dumax,dumax1)
+        CALL allreduce_max_omp(dumax,dumax1)
         dumax=dumax1
       ENDIF  
@@ -312,5 +317,5 @@
       
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -320,16 +325,17 @@
       ENDDO
       
-      IF (is_mpi_root) PRINT *,"gradrot : it :",it ,": dumax",dumax
+      IF (is_master) PRINT *,"gradrot : it :",it ,": dumax",dumax
 
     ENDDO 
-    IF (is_mpi_root) PRINT *,"gradrot : dumax",dumax
+    IF (is_master) PRINT *,"gradrot : dumax",dumax
   
     cgradrot=dumax**(-1./nitergrot) 
-    IF (is_mpi_root) PRINT *,"cgradrot : ",cgradrot
+    IF (is_master) PRINT *,"cgradrot : ",cgradrot
    
 
 
+!$OMP BARRIER
+!$OMP MASTER
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -344,4 +350,6 @@
       ENDDO  
     ENDDO
+!$OMP END MASTER
+!$OMP BARRIER
 
     DO it=1,20
@@ -351,5 +359,5 @@
         CALL transfert_request(f_theta,req_i1)
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           CALL swap_dimensions(ind)
           CALL swap_geometry(ind)
@@ -364,5 +372,5 @@
       
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -379,5 +387,5 @@
 
       IF (using_mpi) THEN 
-        CALL reduce_sum_omp(dthetamax ,dthetamax1)
+        CALL reduce_max_omp(dthetamax ,dthetamax1)
 !$OMP MASTER        
         CALL MPI_ALLREDUCE(dthetamax1,dthetamax,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
@@ -385,12 +393,12 @@
         CALL bcast_omp(dthetamax)  
       ELSE
-        CALL allreduce_sum_omp(dthetamax,dthetamax1)
+        CALL allreduce_max_omp(dthetamax,dthetamax1)
         dumax=dumax1
       ENDIF  
       
-      IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
+      IF (is_master) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
 
       DO ind=1,ndomain
-        IF (.NOT. assigned_domain(ind)) CYCLE
+        IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -401,8 +409,8 @@
     ENDDO 
 
-    IF (is_mpi_root) PRINT *,"divgrad : divgrad",dthetamax
+    IF (is_master) PRINT *,"divgrad : divgrad",dthetamax
   
     cdivgrad=dthetamax**(-1./niterdivgrad) 
-    IF (is_mpi_root) PRINT *,"cdivgrad : ",cdivgrad
+    IF (is_master) PRINT *,"cdivgrad : ",cdivgrad
 
      
@@ -431,9 +439,9 @@
        dtdissip=itau_dissip*dt
     ELSE
-       IF (is_mpi_root) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
+       IF (is_master) PRINT *,"No dissipation time set, setting itau_dissip to 1000000000"
        itau_dissip=100000000
     END IF
     itau_dissip=MAX(1,itau_dissip)
-    IF (is_mpi_root) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
+    IF (is_master) PRINT *,"mintau ",mintau,"itau_dissip",itau_dissip," dtdissip ",dtdissip
 
   END SUBROUTINE init_dissip 
@@ -713,5 +721,5 @@
 
     CALL trace_start("divgrad")
-       
+
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
Index: codes/icosagcm/trunk/src/domain.f90
===================================================================
--- codes/icosagcm/trunk/src/domain.f90	(revision 294)
+++ codes/icosagcm/trunk/src/domain.f90	(revision 295)
@@ -449,11 +449,11 @@
     
     ind=0
-    DO rank=0,omp_size-1
-      nb_domain=ndomain/omp_size
-      IF ( rank < MOD(ndomain,omp_size) ) nb_domain=nb_domain+1
+    DO rank=0,omp_domain_size-1
+      nb_domain=ndomain/omp_domain_size
+      IF ( rank < MOD(ndomain,omp_domain_size) ) nb_domain=nb_domain+1
     
       DO i=1,nb_domain
        ind=ind+1
-       IF (rank==omp_rank) THEN 
+       IF (rank==omp_domain_rank) THEN 
          assigned_domain(ind)=.TRUE.
          PRINT *,"Rank ",mpi_rank," task ",rank," local domain : ",ind," global domain ",domloc_glo_ind(ind)
Index: codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_gcm_mod.f90	(revision 294)
+++ 	(revision )
@@ -1,364 +1,0 @@
-MODULE dynetat0_gcm_mod 
-  USE genmod
-  USE icosa
-  USE caldyn_gcm_mod 
- 	IMPLICIT NONE 
-	  PRIVATE
-
-       PUBLIC  etat0
-	 INTEGER,SAVE::ncell
-	 TYPE(t_field),POINTER:: f_iu(:)
-      TYPE(t_field),POINTER:: f_iv(:) 
-	 TYPE(t_field),POINTER:: f_iue(:)
-      TYPE(t_field),POINTER:: f_ive(:) 
-	 REAL(rstd),POINTER :: iu(:,:),iv(:,:)
-      REAL(rstd),POINTER :: iue(:,:),ive(:,:) 
-
-	  
-CONTAINS
-
-  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
-  USE icosa
-  USE write_field  
-  USE maxicosa
-	IMPLICIT NONE
-	TYPE(t_domain),POINTER :: d 
-	TYPE(t_field),POINTER:: f_ps(:)
-	TYPE(t_field),POINTER:: f_phis(:)
-	TYPE(t_field),POINTER:: f_u(:)
-	TYPE(t_field),POINTER:: f_q(:)
-	TYPE(t_field),POINTER:: f_theta_rhodz(:) 
-     TYPE(t_field),POINTER::  f_buf_i3(:), f_buf1_i(:), f_buf2_i(:)
-     REAL(rstd),POINTER :: ps(:)
-     REAL(rstd),POINTER :: phis(:)
-     REAL(rstd),POINTER :: theta_rhodz(:,:)
-     REAL(rstd),POINTER :: u(:,:) 
-     REAL(rstd),POINTER :: q(:,:,:)
-	REAL(rstd):: maxff,minff,maxuu,minuu
-	INTEGER :: ind
- 
-	CALL allocate_field(f_iu,field_t,type_real,llm) 
-     CALL allocate_field(f_iv,field_t,type_real,llm)
-	CALL allocate_field(f_iue,field_u,type_real,llm) 
-     CALL allocate_field(f_ive,field_u,type_real,llm)
-	CALL allocate_field(f_u,field_u,type_real,llm) 
-	CALL allocate_field(f_buf1_i,field_t,type_real,llm) 
-	CALL allocate_field(f_buf2_i,field_t,type_real,llm) 
-	CALL allocate_field(f_buf_i3,field_u,type_real,3,llm) 
-  
-	PRINT*,"IN NETCDF READ"
-!------------------------------------zero 
-	DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-	  iu = f_iu(ind) 
-	  iv = f_iv(ind) 
-      iue = f_iue(ind) 
-	 ive = f_ive(ind) 	
-	  iu = 0.0 
-	  iv = 0.0 
-        u = f_u(ind) 
-	   u = 0.0 
-	 iue = 0.0	
-	 ive = 0.0 
-     END DO 
-!--------------------------------------------
-	 ncell = 0
-     DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-	 d => domain_glo(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      q=f_q(ind)
-      iu=f_iu(ind) 
-	 iv=f_iv(ind) 
-      CALL compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q)
-	ENDDO
-
-     CALL transfert_request(f_ps,req_i1)
-	CALL transfert_request(f_phis,req_i1)
-	CALL transfert_request(f_theta_rhodz,req_i1)
-	CALL transfert_request(f_q,req_i1)
-	CALL transfert_request(f_iu,req_i1)
-	CALL transfert_request(f_iv,req_i1)
-!------------------------------------------
-	DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-	 u=f_u(ind)
-	 iu=f_iu(ind) 
-	 iv=f_iv(ind) 
-	 iue=f_iue(ind) 
-	 ive=f_ive(ind) 
-  	 CALL compute_dynetatu(iu,iv,iue,ive,u)
-	ENDDO
-!----------------------------------------------------
-!------------- OUTPUT OF Variables 
-  END SUBROUTINE etat0
-
-!==================================================================
-  SUBROUTINE compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q) 
-   use icosa
-   use netcdf
-   use wind_mod 
-   USE disvert_mod
- 	IMPLICIT NONE 
-   TYPE(t_domain),POINTER :: d 
-   CHARACTER*20::dimname 
-   REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
-   REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
-   REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-   REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
-   REAL(rstd),ALLOCATABLE :: mass(:,:)   ! mass    
-   REAL(rstd),ALLOCATABLE :: rhodz(:,:)   ! mass density  
-   REAL(rstd),ALLOCATABLE :: theta(:,:) 
-   REAL(rstd),ALLOCATABLE :: p(:,:)  ! pression
-   REAL(rstd),POINTER :: iu(:,:),iv(:,:)
-   REAL(rstd),POINTER :: icops(:)
-   REAL(rstd),POINTER :: icophis(:)
-   REAL(rstd),POINTER :: icou(:,:),icov(:,:)
-   REAL(rstd),POINTER :: icotheta(:,:)
-   REAL(rstd),POINTER :: icoq(:,:,:)
-
-   INTEGER length,iq,ind,l
-   PARAMETER (length = 100)
-   REAL tab_cntrl(length) ! tableau des parametres du run
-   INTEGER::ierr,nid,ncid,nvarid,dimid,nind
-   INTEGER::ncells  
-   INTEGER::halo_size,i,j,k,ij
-   LOGICAL::single 
-   INTEGER::nDims,nVars,nGlobalAtts,unlimdimid
-   INTEGER:: len
-   CHARACTER(LEN=200):: iqq 
-   
-!	OPEN NETCDF FILE 
-	 ierr = NF90_OPEN ("start_icosa25.nc",NF90_NOWRITE,nid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)'dynetat0: with file start_icosa.nc'
-        write(*,*)' ierr = ', ierr
-        STOP 
-      ENDIF
-
-     ierr= nf90_inquire(nid,nDims,nVars,nGlobalAtts,unlimdimid)
-	IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)'Problem in inquire'
-        write(*,*)' ierr = ', ierr
-        STOP 
-      ENDIF
-
-!	PRINT*,"nDims,nVars,nGlobalAtts,unlimdimid"
-!	PRINT*,nDims,nVars,nGlobalAtts,unlimdimid
-
-
-	 ierr = NF90_INQ_DIMID(nid,"ncells",dimid)
-	 IF (ierr .NE. NF90_NOERR ) THEN 
-	   write(*,*)'ncells is not present in start_icosa.nc'
-        write(*,*)' ierr = ', ierr
-	   STOP 
-	 ENDIF 
-	   
-	 ierr = nf90_inquire_dimension(nid,dimid,dimname,ncells)
-	 IF (ierr .NE. NF90_NOERR ) THEN 
-	   write(*,*)'ncells  in start_icosa.nc'
-        write(*,*)' ierr = ', ierr
-	   STOP 
-	 ENDIF
-
-	  ALLOCATE(icops(ncells))
-	  ALLOCATE(icophis(ncells))
-	  ALLOCATE(icou(ncells,llm))
-	  ALLOCATE(icov(ncells,llm))
-	  ALLOCATE(icotheta(ncells,llm))
-	  ALLOCATE(icoq(ncells,llm,nqtot))
-	  ALLOCATE(p(iim*jjm,llm+1))
-	  ALLOCATE(theta(iim*jjm,llm)) 
-       ALLOCATE(mass(iim*jjm,llm))   ! mass    
-       ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
-!============================================================
-      ierr = NF90_INQ_VARID(nid, "phisinit", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: phisinit is absent"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid, icophis)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN PHIS"
-         STOP
-      ENDIF
-!==============================================================
-	  ierr = NF90_INQ_VARID(nid, "ps", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ps is absent"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid, icops)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN PS"
-         STOP
-      ENDIF
-!================================================================
-	  ierr = NF90_INQ_VARID(nid, "teta", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: teta is not available in start.nc"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icotheta)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN Teta"
-         STOP
-      ENDIF
-!================================================================
-	DO iq = 1,nqtot   
-		write(iqq,*)INT(iq)
-		iqq=ADJUSTL(iqq)  
-     	ierr = NF90_INQ_VARID(nid,"q"//iqq, nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ","q"//iqq,"not here"
-	   write(*,*)' ierr = ', ierr
-!         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,iq))
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN Q"
-!         STOP
-      ENDIF
-	END DO
-!================================================================
-	GO TO 121
-	  ierr = NF90_INQ_VARID(nid, "q01", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: q1 is not available in start.nc"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,1))
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN Q01"
-         STOP
-      ENDIF
-121	CONTINUE 
-!================================================================
-	  ierr = NF90_INQ_VARID(nid, "ucov", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ucov is not available in start.nc"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icou)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN ucov"
-         STOP
-      ENDIF
-	PRINT*,"UCOV is read using start_icosa.nc" 
-!================================================================
-	  ierr = NF90_INQ_VARID(nid, "vcov", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: PROBLEM in VCOV"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icov)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN vcov"
-         STOP
-      ENDIF
-!================================================================
-			iu = 0.0 ; iv = 0.0 
-             DO j=d%jj_begin,d%jj_end
-              DO i=d%ii_begin,d%ii_end
-                  k=d%iim*(j-1)+i
-                 IF (d%assign_domain(i,j)==ind ) THEN
-                     ncell=ncell+1
-		           phis(k)=  0.0  !icophis(ncell) 
-	                ps(k)= icops(ncell) 
-	                theta(k,:) = icotheta(ncell,:) 
-	                q(k,:,1)= icoq(ncell,:,1) 
-	                iu(k,:) = icou(ncell,:)
-				 iv(k,:) = icov(ncell,:) 
-	         ENDIF
-                ENDDO
-              ENDDO
-
-    DO    l    = 1, llm+1
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          p(ij,l) = ap(l) + bp(l) * ps(ij)
-        ENDDO
-      ENDDO
-    ENDDO
-
-   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
-
-    DO    l    = 1, llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          theta_rhodz(ij,l) = theta(ij,l)*rhodz(ij,l)
-        ENDDO
-      ENDDO
-    ENDDO
-
-	  DEALLOCATE(icops)
-	  DEALLOCATE(icophis)
-	  DEALLOCATE(icou)
-	  DEALLOCATE(icov)
-	  DEALLOCATE(icotheta)
-	  DEALLOCATE(p)
-	  DEALLOCATE(theta) 
-       DEALLOCATE(mass)   ! mass    
-       DEALLOCATE(rhodz)   !
-	END SUBROUTINE compute_dynetat0
-
-!==================================================================
-	  SUBROUTINE compute_dynetatu(iu,iv,iue,ive,u) 
-   use icosa
-   use wind_mod 
-    	IMPLICIT NONE 
-   CHARACTER*20::dimname 
-   REAL(rstd),INTENT(OUT):: u(3*iim*jjm,llm)
-   REAL(rstd) :: iu(iim*jjm,llm),iv(iim*jjm,llm)
-   REAL(rstd) :: iue(3*iim*jjm,llm),ive(3*iim*jjm,llm)
-   INTEGER::halo_size,i,j,k,ij,l
-
-
-  Do l = 1, llm
-   DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-         k=iim*(j-1)+i
-	   iue(k+u_right,l)=0.5*(iu(k,l)+iu(k+t_right,l))
-        iue(k+u_lup,l)  =0.5*(iu(k,l)+iu(k+t_lup,l))
-        iue(k+u_ldown,l)=0.5*(iu(k,l)+iu(k+t_ldown,l)) 
-!------------------------------------------------------ 
-	   ive(k+u_right,l)=0.5*(iv(k,l)+iv(k+t_right,l))
-        ive(k+u_lup,l)  =0.5*(iv(k,l)+iv(k+t_lup,l))
-        ive(k+u_ldown,l)=0.5*(iv(k,l)+iv(k+t_ldown,l)) 
-	 END DO
-    END DO 
-  END DO 
-	CALL compute_wind_perp_from_lonlat_compound(iue,ive,u) 
-
-	END SUBROUTINE compute_dynetatu
-
-
- END MODULE dynetat0_gcm_mod 
Index: codes/icosagcm/trunk/src/dynetat0_hz_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/dynetat0_hz_mod.f90	(revision 294)
+++ 	(revision )
@@ -1,348 +1,0 @@
-MODULE dynetat0_hz_mod 
-  USE genmod
-  USE icosa
-  IMPLICIT NONE 
-  PRIVATE
-
-       PUBLIC  etat0
-       INTEGER,SAVE::ncell
-       TYPE(t_field),POINTER:: f_iu(:)
-       TYPE(t_field),POINTER:: f_iv(:) 
-       TYPE(t_field),POINTER:: f_iue(:)
-       TYPE(t_field),POINTER:: f_ive(:) 
-       REAL(rstd),POINTER :: iu(:,:),iv(:,:)
-       REAL(rstd),POINTER :: iue(:,:),ive(:,:) 
-
- 
-CONTAINS
-
-    SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
-  USE icosa
-  USE caldyn_mod
-  USE write_field  
-  USE maxicosa
-  USE wind_mod
-        IMPLICIT NONE
-     TYPE(t_domain),POINTER :: d 
-     TYPE(t_field),POINTER:: f_ps(:)
-     TYPE(t_field),POINTER:: f_phis(:)
-     TYPE(t_field),POINTER:: f_u(:)
-     TYPE(t_field),POINTER:: f_q(:)
-     TYPE(t_field),POINTER:: f_theta_rhodz(:) 
-     TYPE(t_field),POINTER::  f_buf_i3(:), f_buf1_i(:), f_buf2_i(:)
-     REAL(rstd),POINTER :: ps(:)
-     REAL(rstd),POINTER :: phis(:)
-     REAL(rstd),POINTER :: theta_rhodz(:,:)
-     REAL(rstd),POINTER :: u(:,:) 
-     REAL(rstd),POINTER :: q(:,:,:)
-     REAL(rstd):: maxff,minff,maxuu,minuu
-     INTEGER :: ind
- 
-     CALL allocate_field(f_iu,field_t,type_real,llm) 
-     CALL allocate_field(f_iv,field_t,type_real,llm)
-     CALL allocate_field(f_iue,field_u,type_real,llm) 
-     CALL allocate_field(f_ive,field_u,type_real,llm)
-     CALL allocate_field(f_u,field_u,type_real,llm) 
-     CALL allocate_field(f_buf1_i,field_t,type_real,llm) 
-     CALL allocate_field(f_buf2_i,field_t,type_real,llm) 
-     CALL allocate_field(f_buf_i3,field_u,type_real,3,llm) 
-  
-     PRINT*,"IN NETCDF READ"
-!------------------------------------zero 
-	DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-	  iu = f_iu(ind) 
-	  iv = f_iv(ind) 
-      iue = f_iue(ind) 
-	 ive = f_ive(ind) 	
-	   u = f_u(ind) 
-	  iu = 0.0 
-	  iv = 0.0 
-	   u = 0.0 
-	 iue = 0.0	
-	 ive = 0.0 
-     END DO 
-!--------------------------------------------
-     ncell = 0
-     DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      d => domain_glo(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      q=f_q(ind)
-      iu=f_iu(ind) 
-      iv=f_iv(ind) 
-      CALL compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q)
-      ENDDO
-
-     CALL transfert_request(f_ps,req_i1)
-	CALL transfert_request(f_phis,req_i1)
-	CALL transfert_request(f_theta_rhodz,req_i1)
-	CALL transfert_request(f_q,req_i1)
-	CALL transfert_request(f_iu,req_i1)
-	CALL transfert_request(f_iv,req_i1)
-!------------------------------------------
-	DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-	 u=f_u(ind)
-	 iu=f_iu(ind) 
-	 iv=f_iv(ind) 
-	 iue=f_iue(ind) 
-	 ive=f_ive(ind) 
-  	 CALL compute_dynetatu(iu,iv,iue,ive,u)
-	ENDDO
-!------------- OUTPUT OF Variables 
-	CALL un2ulonlat(f_u,f_buf1_i,f_buf2_i) 
-	CALL writefield("buf1",f_buf1_i)
-  END SUBROUTINE etat0
-
-!==================================================================
-  SUBROUTINE compute_dynetat0(ind,d,ps,phis,theta_rhodz,iu,iv,q) 
-   use icosa
-   use netcdf
-   use wind_mod 
-   USE disvert_mod
- 	IMPLICIT NONE 
-   TYPE(t_domain),POINTER :: d 
-   CHARACTER*20::dimname 
-   REAL(rstd), INTENT(OUT) :: ps(iim*jjm)
-   REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
-   REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-   REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
-   REAL(rstd),ALLOCATABLE :: mass(:,:)   ! mass    
-   REAL(rstd),ALLOCATABLE :: rhodz(:,:)   ! mass density  
-   REAL(rstd),ALLOCATABLE :: theta(:,:) 
-   REAL(rstd),ALLOCATABLE :: p(:,:)  ! pression
-   REAL(rstd),POINTER :: iu(:,:),iv(:,:)
-   REAL(rstd),POINTER :: icops(:)
-   REAL(rstd),POINTER :: icophis(:)
-   REAL(rstd),POINTER :: icou(:,:),icov(:,:)
-   REAL(rstd),POINTER :: icotheta(:,:)
-   REAL(rstd),POINTER :: icoq(:,:,:)
-   
-   INTEGER length,iq,ind,l
-   PARAMETER (length = 1)
-   CHARACTER(LEN=200):: iqq ! tableau des parametres du run
-   INTEGER::ierr,nid,ncid,nvarid,dimid,nind
-   INTEGER::ncells  
-   INTEGER::halo_size,i,j,k,ij
-   LOGICAL::single 
-   INTEGER::nDims,nVars,nGlobalAtts,unlimdimid
-   INTEGER:: len
-   
-!	OPEN NETCDF FILE 
-	 ierr = NF90_OPEN ("icosa_hz30.nc",NF90_NOWRITE,nid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)'dynetat0: with file icosa_hz30.nc'
-        write(*,*)' ierr = ', ierr
-        STOP 
-      ENDIF
-
-     ierr= nf90_inquire(nid,nDims,nVars,nGlobalAtts,unlimdimid)
-	IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)'Problem in inquire'
-        write(*,*)' ierr = ', ierr
-        STOP 
-      ENDIF
-
-	PRINT*,"nDims,nVars,nGlobalAtts,unlimdimid"
-	PRINT*,nDims,nVars,nGlobalAtts,unlimdimid
-
-
-	 ierr = NF90_INQ_DIMID(nid,"ncells",dimid)
-	 IF (ierr .NE. NF90_NOERR ) THEN 
-	   write(*,*)'ncells is not present in hzdy_30.nc'
-        write(*,*)' ierr = ', ierr
-	   STOP 
-	 ENDIF 
-	   
-	 ierr = nf90_inquire_dimension(nid,dimid,dimname,ncells)
-	 IF (ierr .NE. NF90_NOERR ) THEN 
-	   write(*,*)'ncells  in hzdy_30.nc'
-        write(*,*)' ierr = ', ierr
-	   STOP 
-	 ENDIF
-
-	  ALLOCATE(icops(ncells))
-	  ALLOCATE(icophis(ncells))
-	  ALLOCATE(icou(ncells,llm))
-	  ALLOCATE(icov(ncells,llm))
-	  ALLOCATE(icotheta(ncells,llm))
-	  ALLOCATE(icoq(ncells,llm,nqtot))
-       ALLOCATE(p(iim*jjm,llm+1))
-	  ALLOCATE(theta(iim*jjm,llm)) 
-       ALLOCATE(mass(iim*jjm,llm))   ! mass    
-       ALLOCATE(rhodz(iim*jjm,llm))   ! mass density   
-!============================================================
-      ierr = NF90_INQ_VARID(nid, "phis", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: phis is absent"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid, icophis)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN PHIS"
-         STOP
-      ENDIF
-!==============================================================
-	  ierr = NF90_INQ_VARID(nid, "ps", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ps is absent"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid, icops)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN PS"
-         STOP
-      ENDIF
-!================================================================
-	  ierr = NF90_INQ_VARID(nid, "theta", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: teta is not available in start.nc"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icotheta)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN Teta"
-         STOP
-      ENDIF
-!================================================================
-	DO iq = 1,nqtot   
-		write(iqq,*)INT(iq)
-		iqq=ADJUSTL(iqq)  
-     	ierr = NF90_INQ_VARID(nid,"q"//iqq, nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ","q"//iqq,"not here"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icoq(:,:,iq))
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN Q"
-         STOP
-      ENDIF
-	END DO
-!===============================================================
-	 ierr = NF90_INQ_VARID(nid, "ulon", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ulon is not in file"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-	
-      ierr = NF90_GET_VAR(nid, nvarid,icou)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN ucov"
-         STOP
-      ENDIF
-!================================================================
-	  ierr = NF90_INQ_VARID(nid, "ulat", nvarid)
-      IF (ierr .NE. NF90_NOERR) THEN
-        write(*,*)"dynetat0: ulat is not available in start.nc"
-	   write(*,*)' ierr = ', ierr
-         STOP       
-      ENDIF
-
-      ierr = NF90_GET_VAR(nid, nvarid,icov)
-      IF (ierr .NE. NF90_NOERR) THEN
-         write(*,*)"dynetat0: PROBLEM IN vlat"
-         STOP
-      ENDIF
-!================================================================
-			iu = 0.0 ; iv = 0.0 
-             DO j=d%jj_begin,d%jj_end
-              DO i=d%ii_begin,d%ii_end
-                  k=d%iim*(j-1)+i
-                 IF (d%assign_domain(i,j)==ind ) THEN
-                     ncell=ncell+1
-				 phis(k)= icophis(ncell) 
-	                ps(k)= icops(ncell) 
-	                theta(k,:) = icotheta(ncell,:) 
-	                q(k,:,:)= icoq(ncell,:,:) 
-	                iu(k,:) = icou(ncell,:)
-				 iv(k,:) = icov(ncell,:)
-			   ENDIF
-                ENDDO
-              ENDDO
-
-    DO    l    = 1, llm+1
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          p(ij,l) = ap(l) + bp(l) * ps(ij)
-        ENDDO
-      ENDDO
-    ENDDO
-
-   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
-
-    DO    l    = 1, llm
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-          theta_rhodz(ij,l) = theta(ij,l)*rhodz(ij,l)
-        ENDDO
-      ENDDO
-    ENDDO
-
-	  DEALLOCATE(icops)
-	  DEALLOCATE(icophis)
-	  DEALLOCATE(icou)
-	  DEALLOCATE(icov)
-	  DEALLOCATE(icotheta)
-	  DEALLOCATE(p)
-	  DEALLOCATE(theta) 
-       DEALLOCATE(mass)   ! mass    
-       DEALLOCATE(rhodz)   !
-	END SUBROUTINE compute_dynetat0
-!==================================================================
-	  SUBROUTINE compute_dynetatu(iu,iv,iue,ive,u) 
-   use icosa
-   use wind_mod 
-    	IMPLICIT NONE 
-   CHARACTER*20::dimname 
-   REAL(rstd),INTENT(OUT):: u(3*iim*jjm,llm)
-   REAL(rstd) :: iu(iim*jjm,llm),iv(iim*jjm,llm)
-   REAL(rstd) :: iue(3*iim*jjm,llm),ive(3*iim*jjm,llm)
-   INTEGER::halo_size,i,j,k,ij,l
-
-
-  Do l = 1, llm
-   DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-         k=iim*(j-1)+i
-	   iue(k+u_right,l)=0.5*(iu(k,l)+iu(k+t_right,l))
-        iue(k+u_lup,l)  =0.5*(iu(k,l)+iu(k+t_lup,l))
-        iue(k+u_ldown,l)=0.5*(iu(k,l)+iu(k+t_ldown,l)) 
-!------------------------------------------------------ 
-	   ive(k+u_right,l)=0.5*(iv(k,l)+iv(k+t_right,l))
-        ive(k+u_lup,l)  =0.5*(iv(k,l)+iv(k+t_lup,l))
-        ive(k+u_ldown,l)=0.5*(iv(k,l)+iv(k+t_ldown,l)) 
-	 END DO
-    END DO 
-  END DO 
-	CALL compute_wind_perp_from_lonlat_compound(iue,ive,u) 
-
-	END SUBROUTINE compute_dynetatu
-
-
- END MODULE dynetat0_hz_mod 
Index: codes/icosagcm/trunk/src/etat0.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0.f90	(revision 294)
+++ codes/icosagcm/trunk/src/etat0.f90	(revision 295)
@@ -25,6 +25,4 @@
     USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>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 
     USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0  
 
@@ -80,5 +78,5 @@
     CASE ('dcmip3')
        CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-     CASE ('dcmip4')
+    CASE ('dcmip4')
         IF(nqtot<2) THEN
            IF (is_mpi_root)  THEN
@@ -87,11 +85,5 @@
            STOP
         END IF
-       CALL etat0_dcmip4(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) 
+        CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE DEFAULT
        PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
@@ -114,4 +106,5 @@
 
   SUBROUTINE etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_q)
+    USE theta2theta_rhodz_mod
     IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -119,4 +112,5 @@
     TYPE(t_field),POINTER :: f_phis(:)
     TYPE(t_field),POINTER :: f_theta_rhodz(:)
+    TYPE(t_field),POINTER :: f_temp(:)
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
@@ -126,4 +120,5 @@
     REAL(rstd),POINTER :: phis(:)
     REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd),POINTER :: u(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
@@ -138,12 +133,21 @@
       phis=f_phis(ind)
       theta_rhodz=f_theta_rhodz(ind)
+      temp=f_temp(ind)
       u=f_u(ind)
       q=f_q(ind)
-      CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
+
+      IF( TRIM(etat0_type)=='williamson91.6' ) THEN 
+       CALL compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
+      ELSE
+       CALL compute_etat0_collocated(ps,mass, phis, temp, u, q)
+      ENDIF
     ENDDO
+    
+    IF( TRIM(etat0_type)/='williamson91.6' ) CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+    
+    
   END SUBROUTINE etat0_collocated
 
-  SUBROUTINE compute_etat0_collocated(ps,mass, phis, theta_rhodz, u, q)
-    USE theta2theta_rhodz_mod
+  SUBROUTINE compute_etat0_collocated(ps,mass, phis, temp_i, u, q)
     USE wind_mod
     USE etat0_jablonowsky06_mod, ONLY : compute_jablonowsky06 => compute_etat0
@@ -154,12 +158,15 @@
     REAL(rstd),INTENT(INOUT) :: mass(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
-    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: temp_i(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm,nqtot)
 
-    REAL(rstd) :: temp_i(iim*jjm,llm)
+    REAL(rstd) :: lon_i(iim*jjm)
+    REAL(rstd) :: lat_i(iim*jjm)
     REAL(rstd) :: ulon_i(iim*jjm,llm)
     REAL(rstd) :: ulat_i(iim*jjm,llm)
 
+    REAL(rstd) :: lon_e(3*iim*jjm)
+    REAL(rstd) :: lat_e(3*iim*jjm)
     REAL(rstd) :: ps_e(3*iim*jjm)
     REAL(rstd) :: mass_e(3*iim*jjm,llm)
@@ -183,15 +190,8 @@
        CALL compute_dcmip5(3*iim*jjm,lon_e,lat_e, phis_e, ps_e, temp_e, ulon_e, ulat_e, q_e)
     CASE('williamson91.6')
-       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), theta_rhodz(:,1), ulon_i(:,1), ulat_i(:,1))
+       CALL compute_w91_6(iim*jjm,lon_i,lat_i, phis, mass(:,1), temp_i(:,1), ulon_i(:,1), ulat_i(:,1))
        CALL compute_w91_6(3*iim*jjm,lon_e,lat_e, phis_e, mass_e(:,1), temp_e(:,1), ulon_e(:,1), ulat_e(:,1))
     END SELECT
 
-    SELECT CASE (TRIM(etat0_type))
-    CASE('williamson91.6')
-       ! Do nothing
-    CASE DEFAULT
-       CALL compute_temperature2theta_rhodz(ps,temp_i,theta_rhodz,0)    
-    END SELECT
- 
     CALL compute_wind_perp_from_lonlat_compound(ulon_e, ulat_e, u)
 
Index: codes/icosagcm/trunk/src/etat0_academic.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_academic.f90	(revision 294)
+++ codes/icosagcm/trunk/src/etat0_academic.f90	(revision 295)
@@ -134,6 +134,9 @@
      
 
+!$OMP BARRIER
     CALL compute_pression(ps,p,1)     
+!$OMP BARRIER
     CALL compute_exner(ps,p,pks,pk,1)  
+!$OMP BARRIER
     CALL compute_geopotential(phis,pks,pk,theta,phi,1)
 
Index: codes/icosagcm/trunk/src/etat0_dcmip2.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 294)
+++ codes/icosagcm/trunk/src/etat0_dcmip2.f90	(revision 295)
@@ -14,4 +14,5 @@
   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
   USE icosa
+  USE theta2theta_rhodz_mod
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -20,9 +21,10 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
-    
+        
+    TYPE(t_field),POINTER,SAVE :: f_temp(:)
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: phis(:)
     REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: Temp(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
 
@@ -45,4 +47,6 @@
     PRINT *, 'Orographic gravity-wave test case :', TRIM(etat0_type)
 
+    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 
+
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -53,12 +57,16 @@
       phis=f_phis(ind)
       u=f_u(ind)
-      theta_rhodz=f_theta_rhodz(ind)
       q=f_q(ind)
-      CALL compute_etat0_DCMIP2(icase,ps,phis,u,theta_rhodz,q)
+      temp=f_temp(ind)
+      CALL compute_etat0_DCMIP2(icase,ps,phis,u,Temp,q)
     ENDDO
+    
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
             
+    CALL deallocate_field(f_temp)
+   
   END SUBROUTINE etat0
   
-  SUBROUTINE compute_etat0_DCMIP2(icase, ps, phis, u, theta_rhodz,q)
+  SUBROUTINE compute_etat0_DCMIP2(icase, ps, phis, u, temp,q)
     USE icosa
     USE disvert_mod, ONLY : ap,bp
@@ -71,7 +79,7 @@
     REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
     REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm)
-    REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd), INTENT(OUT) :: temp(iim*jjm,llm)
     REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm)
-    REAL(rstd) :: ulon(3*iim*jjm,llm), ulat(3*iim*jjm,llm), temp(iim*jjm,llm)
+    REAL(rstd) :: ulon(3*iim*jjm,llm), ulat(3*iim*jjm,llm)
     
     INTEGER :: i,j,l,ij
@@ -90,5 +98,4 @@
        END DO
     END DO
-    CALL compute_temperature2theta_rhodz(ps,temp,theta_rhodz,1)
     
     ! Edges : ulon,ulat
Index: codes/icosagcm/trunk/src/etat0_dcmip3.f90
===================================================================
--- codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 294)
+++ codes/icosagcm/trunk/src/etat0_dcmip3.f90	(revision 295)
@@ -18,4 +18,5 @@
   SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
   USE icosa
+  USE theta2theta_rhodz_mod
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -24,12 +25,15 @@
     TYPE(t_field),POINTER :: f_u(:)
     TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER,SAVE :: f_temp(:)
     
     REAL(rstd),POINTER :: ps(:)
     REAL(rstd),POINTER :: phis(:)
     REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
+    REAL(rstd),POINTER :: Temp(:,:)
     REAL(rstd),POINTER :: q(:,:,:)
 
     INTEGER :: ind
+
+    CALL allocate_field(f_temp,field_t,type_real,llm,name='temp') 
 
     DO ind=1,ndomain
@@ -41,13 +45,16 @@
       phis=f_phis(ind)
       u=f_u(ind)
-      theta_rhodz=f_theta_rhodz(ind)
       q=f_q(ind)
-      CALL compute_etat0_DCMIP3(ps,phis,u,theta_rhodz,q)
+      temp=f_temp(ind)
+      CALL compute_etat0_DCMIP3(ps,phis,u,Temp,q)
     ENDDO
+
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+    CALL deallocate_field(f_temp)
             
   END SUBROUTINE etat0
   
 
-  SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, theta_rhodz,q)
+  SUBROUTINE compute_etat0_DCMIP3(ps, phis, u, temp,q)
   USE icosa
   USE pression_mod
@@ -68,10 +75,9 @@
   REAL(rstd), INTENT(OUT) :: phis(iim*jjm)
   REAL(rstd), INTENT(OUT) :: u(3*iim*jjm,llm)
-  REAL(rstd), INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
+  REAL(rstd), INTENT(OUT) :: Temp(iim*jjm,llm)
   REAL(rstd), INTENT(OUT) :: q(iim*jjm,llm,nqtot)
   
   REAL(rstd) :: Ts(iim*jjm)
   REAL(rstd) :: s(iim*jjm)
-  REAL(rstd) :: T(iim*jjm,llm)
   REAL(rstd) :: p(iim*jjm,llm+1)
   REAL(rstd) :: theta(iim*jjm,llm)
@@ -125,5 +131,7 @@
   END DO
   
+!$OMP BARRIER
   CALL compute_pression(ps,p,0)
+!$OMP BARRIER
   
   DO l=1,llm
@@ -134,5 +142,5 @@
            IF(use_dcmip_routine) THEN
               CALL test3_gravity_wave(lon_i(ij),lat_i(ij),pp,dummy,0, &
-                   dummy,dummy,dummy,T(ij,l),dummy,dummy,dummy,dummy)
+                   dummy,dummy,dummy,Temp(ij,l),dummy,dummy,dummy,dummy)
            ELSE
               pspsk=(pp/ps(ij))**kappa
@@ -142,5 +150,5 @@
               thetap = dtheta *sin(2*Pi*zz/Lz) * s(ij)                ! perturbation pot. temp.
               theta(ij,l) = thetab + thetap
-              T(ij,l) = theta(ij,l)* ((pp/peq)**kappa)
+              Temp(ij,l) = theta(ij,l)* ((pp/peq)**kappa)
               ! T(ij,l) = Ts(ij)*pspsk / ( Ts(ij) / GG * ( pspsk-1) +1)  ! background temp.
            END IF
@@ -149,9 +157,9 @@
   ENDDO
   
-  IF(use_dcmip_routine) THEN
-     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
-  ELSE
-     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
-  END IF
+!  IF(use_dcmip_routine) THEN
+!     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
+!  ELSE
+!     CALL compute_temperature2theta_rhodz(ps,T,theta_rhodz,0)
+!  END IF
   
   pp=peq
Index: codes/icosagcm/trunk/src/exner.f90
===================================================================
--- codes/icosagcm/trunk/src/exner.f90	(revision 294)
+++ codes/icosagcm/trunk/src/exner.f90	(revision 295)
@@ -22,4 +22,5 @@
     INTEGER :: ind
 
+!$OMP BARRIER
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -32,4 +33,5 @@
       CALL compute_exner(ps, p, pks, pk, 0)
     ENDDO
+!$OMP BARRIER
   
   END SUBROUTINE exner
@@ -39,4 +41,5 @@
   USE disvert_mod
   USE pression_mod
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: ps(iim*jjm)
@@ -52,66 +55,74 @@
        !! Compute Alpha and Beta
        
+       IF (is_omp_level_master) THEN
        ! for llm layer
-       DO j=jj_begin-offset,jj_end+offset
-          DO i=ii_begin-offset,ii_end+offset
-             ij=(j-1)*iim+i
-             alpha(ij,llm) = 0.
-             beta (ij,llm) = 1./ (1+ 2*kappa)
-          ENDDO
-       ENDDO
+         DO j=jj_begin-offset,jj_end+offset
+            DO i=ii_begin-offset,ii_end+offset
+               ij=(j-1)*iim+i
+               alpha(ij,llm) = 0.
+               beta (ij,llm) = 1./ (1+ 2*kappa)
+            ENDDO
+         ENDDO
        
        ! for other layer   
-       DO l = llm-1 , 2 , -1
+         DO l = llm-1 , 2 , -1
+            DO j=jj_begin-offset,jj_end+offset
+               DO i=ii_begin-offset,ii_end+offset
+                  ij=(j-1)*iim+i
+                  delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
+                  alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
+                  beta (ij,l)  =   p(ij,l  ) / delta   
+               ENDDO
+            ENDDO
+         ENDDO
+       
+         !! Compute pk
+       
+         ! for first layer
+         DO j=jj_begin-offset,jj_end+offset
+            DO i=ii_begin-offset,ii_end+offset
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+               pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    &
+                    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
+            ENDDO
+         ENDDO
+       
+       ! for other layers
+       
+         DO l = 2, llm
+            DO j=jj_begin-offset,jj_end+offset
+               DO i=ii_begin-offset,ii_end+offset
+                  ij=(j-1)*iim+i
+                  pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
+               ENDDO
+            ENDDO
+         ENDDO
+
+       ENDIF
+
+    ELSE ! Simple calculation of Exner pressure based on centered average
+       ! surface : pks
+       IF (is_omp_level_master) THEN
+
+         DO j=jj_begin-offset,jj_end+offset
+            DO i=ii_begin-offset,ii_end+offset
+               ij=(j-1)*iim+i
+               pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
+            ENDDO
+         ENDDO
+
+       ENDIF
+
+         ! 3D : pk
+       DO l = 1, llm
           DO j=jj_begin-offset,jj_end+offset
              DO i=ii_begin-offset,ii_end+offset
-                ij=(j-1)*iim+i
-                delta = p(ij,l)* (1+2*kappa) + p(ij,l+1)* ( beta(ij,l+1)- (1+2*kappa) )
-                alpha(ij,l)  = - p(ij,l+1) / delta * alpha(ij,l+1)
-                beta (ij,l)  =   p(ij,l  ) / delta   
+               ij=(j-1)*iim+i
+               pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
              ENDDO
           ENDDO
        ENDDO
-       
-       !! Compute pk
-       
-       ! for first layer
-       DO j=jj_begin-offset,jj_end+offset
-          DO i=ii_begin-offset,ii_end+offset
-             ij=(j-1)*iim+i
-             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-             pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )     /    &
-                  (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-(1.+2*kappa) )* p(ij,2) )
-          ENDDO
-       ENDDO
-       
-       ! for other layers
-       
-       DO l = 2, llm
-          DO j=jj_begin-offset,jj_end+offset
-             DO i=ii_begin-offset,ii_end+offset
-                ij=(j-1)*iim+i
-                pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
-             ENDDO
-          ENDDO
-       ENDDO
-       
-    ELSE ! Simple calculation of Exner pressure based on centered average
-       ! surface : pks
-       DO j=jj_begin-offset,jj_end+offset
-          DO i=ii_begin-offset,ii_end+offset
-             ij=(j-1)*iim+i
-             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
-          ENDDO
-       ENDDO
 
-       ! 3D : pk
-       DO l = 1, llm
-          DO j=jj_begin-offset,jj_end+offset
-             DO i=ii_begin-offset,ii_end+offset
-                ij=(j-1)*iim+i
-                pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
-             ENDDO
-          ENDDO
-       ENDDO
     END IF
     
Index: codes/icosagcm/trunk/src/field.f90
===================================================================
--- codes/icosagcm/trunk/src/field.f90	(revision 294)
+++ codes/icosagcm/trunk/src/field.f90	(revision 295)
@@ -48,4 +48,5 @@
   SUBROUTINE allocate_field(field,field_type,data_type,dim1,dim2,name)
   USE domain_mod
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: field(:)
@@ -64,5 +65,5 @@
 
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
 
       IF(PRESENT(name)) THEN
@@ -181,4 +182,5 @@
   SUBROUTINE deallocate_field(field)
   USE domain_mod
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: field(:)
@@ -188,5 +190,5 @@
 !$OMP BARRIER
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
 
       data_type=field(ind)%data_type
Index: codes/icosagcm/trunk/src/geometry.f90
===================================================================
--- codes/icosagcm/trunk/src/geometry.f90	(revision 294)
+++ codes/icosagcm/trunk/src/geometry.f90	(revision 295)
@@ -166,4 +166,5 @@
     USE vector
     USE transfert_mod
+    USE omp_para
 
     IMPLICIT NONE
@@ -186,5 +187,5 @@
     
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
@@ -212,4 +213,5 @@
   USE vector
   USE getin_mod
+  USE omp_para
   IMPLICIT NONE
     INTEGER :: nb_it=0
@@ -226,5 +228,5 @@
     
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
       d=>domain(ind)
       CALL swap_dimensions(ind)
@@ -241,5 +243,5 @@
 
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
       d=>domain(ind)
       CALL swap_dimensions(ind)
@@ -271,5 +273,5 @@
       sum=0
       DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
         CALL swap_dimensions(ind)
         CALL swap_geometry(ind)
@@ -311,4 +313,5 @@
   USE transfert_mod
   USE getin_mod
+  USE omp_para
   IMPLICIT NONE
 
@@ -330,5 +333,5 @@
  
     DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
+      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
       d=>domain(ind)
       CALL swap_dimensions(ind)
@@ -559,5 +562,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)
+!    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
  
   END SUBROUTINE set_geometry
Index: codes/icosagcm/trunk/src/geopotential_mod.f90
===================================================================
--- codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 294)
+++ codes/icosagcm/trunk/src/geopotential_mod.f90	(revision 295)
@@ -35,4 +35,5 @@
   SUBROUTINE compute_geopotential(phis,pks,pk,theta,phi,offset)
   USE icosa
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: phis(iim*jjm)
@@ -51,21 +52,28 @@
    
   ! for first layer
-   DO j=jj_begin-offset,jj_end+offset
-     DO i=ii_begin-offset,ii_end+offset
-       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
+
+! flush pks, pk thetha, phis
+!$OMP BARRIER
+   IF(is_omp_level_master) THEN
      DO j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          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,1 ) = phis( ij ) + theta(ij,1) * ( pks(ij) - pk(ij,1) )
        ENDDO
      ENDDO
-   ENDDO
+          
+    ! for other layers
+     DO l = 2, llm
+       DO j=jj_begin-offset,jj_end+offset
+         DO i=ii_begin-offset,ii_end+offset
+           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
+   ENDIF
+! flush phi
+!$OMP BARRIER
    
   END SUBROUTINE compute_geopotential
Index: codes/icosagcm/trunk/src/icosa_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 294)
+++ codes/icosagcm/trunk/src/icosa_gcm.f90	(revision 295)
@@ -19,6 +19,6 @@
   CALL xios_init
   CALL init_earth_const 
-  CALL init_grid_param(is_mpi_root)
-  CALL init_omp_para
+  CALL init_grid_param(is_mpi_master)
+  CALL init_omp_para(is_mpi_master)
   CALL compute_metric
   CALL compute_domain
@@ -28,4 +28,5 @@
 
 !$OMP PARALLEL  
+  CALL switch_omp_no_distrib_level
   CALL compute_geometry
   CALL check_total_area
@@ -43,5 +44,5 @@
   
   CALL timeloop
-
+  CALL switch_omp_no_distrib_level
 !$OMP END PARALLEL
 
Index: codes/icosagcm/trunk/src/kinetic.f90
===================================================================
--- codes/icosagcm/trunk/src/kinetic.f90	(revision 294)
+++ codes/icosagcm/trunk/src/kinetic.f90	(revision 295)
@@ -30,4 +30,5 @@
   SUBROUTINE compute_kinetic(ue, Ki)
   USE icosa
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
@@ -35,5 +36,5 @@
     INTEGER :: i,j,ij,l
     
-    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/maxico.f90
===================================================================
--- codes/icosagcm/trunk/src/maxico.f90	(revision 294)
+++ 	(revision )
@@ -1,218 +1,0 @@
- MODULE maxicosa
-contains
-    SUBROUTINE maxis(f_fl,maxff,minff) 
-     USE icosa
-     IMPLICIT NONE
- INTEGER::i,j,ij,l,ind
- TYPE(t_field),POINTER::f_fl(:) 
- REAL(rstd),POINTER::fl(:)  
- REAL(rstd)::maxf,minf 
- REAL(rstd),INTENT(OUT)::maxff,minff 
- 
-    maxff=-1e50
-    minff=1e50 
-    DO ind=1,ndomain 
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind) 
-     fl=f_fl(ind) 
-     CALL compute_maxis(fl,maxf,minf)  
-     maxff=max(maxff,maxf) 
-     minff=min(minff,minf) 
-    ENDDO 
-   END SUBROUTINE maxis
-!------------------------------------------------------------
-	    SUBROUTINE maxiv(f_fu,maxuu,minuu) 
-     USE icosa
-     IMPLICIT NONE
- INTEGER::i,j,ij,l,ind
- TYPE(t_field),POINTER::f_fu(:) 
- REAL(rstd),POINTER::fu(:,:)  
- REAL(rstd)::maxu,minu 
- REAL(rstd),INTENT(OUT)::maxuu,minuu 
- 
-    maxuu=-1e50
-    minuu=1e50 
-    DO ind=1,ndomain 
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind) 
-     fu=f_fu(ind) 
-     CALL compute_maxiv(fu,maxu,minu)  
-	maxuu=max(maxuu,maxu) 
-	minuu=min(minuu,minu) 
-    ENDDO 
-	END SUBROUTINE maxiv
-
-    SUBROUTINE errornorms(f_ref,f_sim,norml1,norml2,normlf) 
-     USE icosa
-     IMPLICIT NONE
- INTEGER::i,j,ij,l,ind
- TYPE(t_field),POINTER::f_ref(:),f_sim(:) 
- REAL(rstd),POINTER::ref(:),sim(:)  
- REAL(rstd),INTENT(OUT)::norml1,norml2,normlf 
- REAL(rstd)::maxf,minf
- REAL(rstd)::maxff,minff
- 
-    norml1 = 0.0 ; norml2 = 0.0 ; normlf = 0.0
-    DO ind=1,ndomain 
-     CALL swap_dimensions(ind)
-     CALL swap_geometry(ind) 
-      ref=f_ref(ind) 
-      sim=f_sim(ind) 
-     CALL compute_maxis(sim,maxf,minf)  
-     maxff=max(maxff,maxf) 
-     minff=min(minff,minf) 
-    ENDDO 
-	END SUBROUTINE errornorms
-
-	SUBROUTINE compute_maxis(fl,maxf,minf) 
-	USE icosa
-	IMPLICIT NONE 
-	REAL(rstd),INTENT(IN)::fl(iim*jjm)
-	REAL(rstd),INTENT(OUT)::maxf,minf 
-	INTEGER::i,j,ij 
-
-	maxf=-1e50
-	minf= 1e50 
-     DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-           ij=(j-1)*iim+i
-          IF ( maxf .LT. fl(ij) ) Then 
-	      maxf = fl(ij) 
-	     ENDIF 
-	    IF ( minf .GT. fl(ij) ) Then 
-	     minf = fl(ij) 
-	    ENDIF 
-	   ENDDO
-      ENDDO
-	END SUBROUTINE compute_maxis
-!------------------------------------------------------
-		SUBROUTINE compute_maxiv(fu,maxu,minu) 
-	USE icosa
-	IMPLICIT NONE 
-	REAL(rstd),INTENT(IN)::fu(3*iim*jjm,llm) 
-	REAL(rstd),INTENT(OUT)::maxu,minu  
-	INTEGER::i,j,ij 
-
-     maxu=-1e50
-     minu=1e50 
-	DO j=jj_begin-1,jj_end+1
-         DO i=ii_begin-1,ii_end+1
-           ij=(j-1)*iim+i
-          IF ( maxu .LT. fu(ij+u_right,llm) ) Then 
-	      maxu = fu(ij+u_right,llm) 
-	     ENDIF 
-		IF ( maxu .LT. fu(ij+u_lup,llm) ) Then 
-	      maxu = fu(ij+u_lup,llm) 
-	     ENDIF
-	     IF ( maxu .LT. fu(ij+u_ldown,llm) ) Then 
-	      maxu = fu(ij+u_ldown,llm) 
-	     ENDIF
-!------------------------------
-		IF ( minu .GT. fu(ij+u_right,llm) ) Then 
-	      minu = fu(ij+u_right,llm) 
-	     ENDIF 
-		IF ( minu .GT. fu(ij+u_lup,llm) ) Then 
-	      minu = fu(ij+u_lup,llm) 
-	     ENDIF
-	     IF ( maxu .GT. fu(ij+u_ldown,llm) ) Then 
-	      minu = fu(ij+u_ldown,llm) 
-	     ENDIF 
-	   ENDDO
-      ENDDO
-	END SUBROUTINE compute_maxiv
-!--------------------------------------------------
-	Subroutine compute_l1(reference,simulated,nrml1)
-	use icosa
-	IMPLICIT NONE 
-	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
-	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
-	REAL(rstd),INTENT(OUT)::nrml1
-	REAL(rstd) :: temp1,temp2
-	INTEGER::i,j,ij 
-
-	temp1 = 0.0 ; temp2 = 0.0 
-	nrml1 = 0.0 
-	DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-	  temp1 = temp1 + abs(simulated(ij) - reference(ij))*Ai(ij) 
-	  temp2 = temp2 + abs(reference(ij))*Ai(ij) 
-	 END DO
-	END DO
-	  nrml1 = temp1/temp2 
-	END SUBROUTINE compute_l1 
-
-	Subroutine compute_l2(reference,simulated,nrml2)
-	use icosa
-	IMPLICIT NONE 
-	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
-	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
-	REAL(rstd),INTENT(OUT)::nrml2
-	REAL(rstd) :: temp1,temp2
-	INTEGER::i,j,ij 
-
-	temp1 = 0.0 ; temp2 = 0.0 
-	nrml2 = 0.0 
-	DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-          ij=(j-1)*iim+i
-	  temp1 = temp1 + (simulated(ij) - reference(ij))*(simulated(ij) - reference(ij))*Ai(ij) 
-	  temp2 = temp2 + reference(ij)*reference(ij)*Ai(ij) 
-	 END DO
-	END DO
-	  temp1 = sqrt(temp1) 
-	  temp2 = sqrt(temp2)  
-	  nrml2 = temp1/temp2 
-	END SUBROUTINE compute_l2 
-
-	Subroutine compute_lf(reference,simulated,nrmlf)
-	use icosa
-	IMPLICIT NONE 
-	REAL(rstd),INTENT(IN)::reference(iim*jjm) 
-	REAL(rstd),INTENT(IN)::simulated(iim*jjm) 
-	REAL(rstd),INTENT(OUT)::nrmlf
-	REAL(rstd):: difference(iim*jjm) 
-	REAL(rstd) :: temp1,temp2
-	INTEGER::i,j,ij 
-
-	temp1 = 0.0 ; temp2 = 0.0 
-	nrmlf = 0.0 
-	
-	temp1 = maxval(ABS(difference))
-	temp2 = maxval(ABS(reference)) 
-	  nrmlf = temp1/temp2 
-	END SUBROUTINE compute_lf 
-!================================================== initial0 
-	 SUBROUTINE initial0(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
-  USE icosa
-  IMPLICIT NONE
-    TYPE(t_field),POINTER :: f_ps(:)
-    TYPE(t_field),POINTER :: f_phis(:)
-    TYPE(t_field),POINTER :: f_theta_rhodz(:)
-    TYPE(t_field),POINTER :: f_u(:)
-    TYPE(t_field),POINTER :: f_q(:)
-  
-    REAL(rstd),POINTER :: ps(:)
-    REAL(rstd),POINTER :: phis(:)
-    REAL(rstd),POINTER :: theta_rhodz(:,:)
-    REAL(rstd),POINTER :: u(:,:)
-    REAL(rstd),POINTER :: q(:,:,:)
-    INTEGER :: ind
-       
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
-      phis=f_phis(ind)
-      theta_rhodz=f_theta_rhodz(ind)
-      u=f_u(ind)
-      q=f_q(ind)
-      q=0.0
-      ps=0.0
-	 phis=0.0
-      theta_rhodz = 0.0
-	 u = 0.0
-    ENDDO
-  END SUBROUTINE initial0
-
-	END MODULE maxicosa
Index: codes/icosagcm/trunk/src/omega.f90
===================================================================
--- codes/icosagcm/trunk/src/omega.f90	(revision 294)
+++ codes/icosagcm/trunk/src/omega.f90	(revision 295)
@@ -23,6 +23,10 @@
   END SUBROUTINE W_omega
 
+
+
   SUBROUTINE compute_omega(ps,u, w)
     USE disvert_mod, ONLY : ap,bp
+    USE omp_para
+    IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: u(iim*3*jjm,llm), ps(iim*jjm)
     REAL(rstd),INTENT(OUT):: w(iim*jjm,llm)
@@ -30,63 +34,68 @@
     REAL(rstd):: p(iim*jjm,llm+1), rhodz(iim*jjm,llm), Fe(iim*3*jjm,llm)
     REAL(rstd):: ugradps
-    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
+    
+    INTEGER :: i,j,l,ij
+
+!$OMP BARRIER    
+    IF (is_omp_level_master) THEN
+      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
+      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 mass flux
-    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
-             Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-             Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-             Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
-          ENDDO
-       ENDDO
-    ENDDO
+      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
+               Fe(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
+               Fe(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
+               Fe(ij+u_ldown,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)*le(ij+u_ldown)
+            ENDDO
+         ENDDO
+      ENDDO
     
 !!! mass flux convergence computation  
     
     ! horizontal convergence
-    DO l = 1, llm
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
-             convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
-                  ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
-                  ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
-                  ne(ij,left)*Fe(ij+u_left,l)     +  &  
-                  ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
-                  ne(ij,rdown)*Fe(ij+u_rdown,l))
-          ENDDO
-       ENDDO
-    ENDDO
+      DO l = 1, llm
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               ! convm = +div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
+               convm(ij,l)= 1./Ai(ij)*(ne(ij,right)*Fe(ij+u_right,l)   +  &
+                    ne(ij,rup)*Fe(ij+u_rup,l)       +  &  
+                    ne(ij,lup)*Fe(ij+u_lup,l)       +  &  
+                    ne(ij,left)*Fe(ij+u_left,l)     +  &  
+                    ne(ij,ldown)*Fe(ij+u_ldown,l)   +  &  
+                    ne(ij,rdown)*Fe(ij+u_rdown,l))
+            ENDDO
+         ENDDO
+      ENDDO
     
-    ! vertical integration from up to down
-    DO  l = llm-1, 1, -1
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-          ENDDO
-       ENDDO
-    ENDDO
-    convm(:,llm+1)=0.
+      ! vertical integration from up to down
+      DO  l = llm-1, 1, -1
+         DO j=jj_begin,jj_end
+            DO i=ii_begin,ii_end
+               ij=(j-1)*iim+i
+               convm(ij,l) = convm(ij,l) + convm(ij,l+1)
+            ENDDO
+         ENDDO
+      ENDDO
+      convm(:,llm+1)=0.
 
 !!! Compute dps
@@ -125,8 +134,7 @@
     ! -grad ps : ( ne(ij,ldown)*ps(ij,l) + ne(ij+t_ldown,rup)*ps(ij+t_ldown,l) ) ) / de(ij+u_ldown)
 
-    DO l = 1,llm
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             toto = 1
+       DO l = 1,llm
+         DO j=jj_begin,jj_end 
+           DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
              ugradps = &
@@ -140,6 +148,8 @@
              w( ij, l) =  ugradps - .5*(convm( ij,l+1)+convm(ij,l)) 
           ENDDO
-       ENDDO
-    ENDDO
+        ENDDO
+      ENDDO
+    ENDIF
+!$OMP BARRIER
     
   END SUBROUTINE compute_omega
Index: codes/icosagcm/trunk/src/omp_para.F90
===================================================================
--- codes/icosagcm/trunk/src/omp_para.F90	(revision 294)
+++ codes/icosagcm/trunk/src/omp_para.F90	(revision 295)
@@ -5,8 +5,8 @@
 !$OMP THREADPRIVATE(omp_rank)
 
-  LOGICAL,SAVE :: omp_first
-  LOGICAL,SAVE :: omp_last
-  LOGICAL,SAVE :: omp_master
-!$OMP THREADPRIVATE(omp_first, omp_last,omp_master)
+  LOGICAL,SAVE :: is_omp_first_level
+  LOGICAL,SAVE :: is_omp_last_level
+  LOGICAL,SAVE :: is_omp_master
+!$OMP THREADPRIVATE(is_omp_first_level, is_omp_last_level,is_omp_master)
 
   INTEGER,SAVE :: ll_begin
@@ -17,16 +17,57 @@
 !$OMP THREADPRIVATE(ll_begin,ll_beginp1,ll_end,ll_endm1,ll_endp1)
   LOGICAL,SAVE :: using_openmp
-
+  
+  INTEGER,SAVE :: omp_domain_size
+  INTEGER,SAVE :: omp_domain_rank
+  INTEGER,SAVE :: omp_level_size
+  INTEGER,SAVE :: omp_level_rank
+!$OMP THREADPRIVATE( omp_domain_size, omp_level_size,omp_domain_rank,omp_level_rank)
+  LOGICAL,SAVE :: is_omp_domain_master
+  LOGICAL,SAVE :: is_omp_level_master
+!$OMP THREADPRIVATE(is_omp_domain_master,is_omp_level_master )
+  
   LOGICAL,PARAMETER :: omp_by_domain=.TRUE. 
+  LOGICAL,SAVE :: is_master
+!$OMP THREADPRIVATE(is_master)
+
+
+   LOGICAL,SAVE :: is_omp_first_level_full
+   LOGICAL,SAVE :: is_omp_last_level_full
+   INTEGER,SAVE :: ll_begin_full
+   INTEGER,SAVE :: ll_beginp1_full
+   INTEGER,SAVE :: ll_end_full
+   INTEGER,SAVE :: ll_endm1_full
+   INTEGER,SAVE :: ll_endp1_full
+!$OMP THREADPRIVATE(is_omp_first_level_full,is_omp_last_level_full)   
+!$OMP THREADPRIVATE( ll_begin_full, ll_beginp1_full, ll_end_full, ll_endm1_full, ll_endp1_full)
+  PRIVATE :: is_omp_first_level_full,is_omp_last_level_full
+  PRIVATE :: ll_begin_full, ll_beginp1_full, ll_end_full, ll_endm1_full, ll_endp1_full 
+  
+
+   LOGICAL,SAVE :: is_omp_first_level_distrib
+   LOGICAL,SAVE :: is_omp_last_level_distrib
+   INTEGER,SAVE :: ll_begin_distrib
+   INTEGER,SAVE :: ll_beginp1_distrib
+   INTEGER,SAVE :: ll_end_distrib
+   INTEGER,SAVE :: ll_endm1_distrib
+   INTEGER,SAVE :: ll_endp1_distrib
+!$OMP THREADPRIVATE(is_omp_first_level_distrib,is_omp_last_level_distrib)   
+!$OMP THREADPRIVATE( ll_begin_distrib, ll_beginp1_distrib, ll_end_distrib, ll_endm1_distrib, ll_endp1_distrib)
+
+  PRIVATE :: is_omp_first_level_distrib,is_omp_last_level_distrib
+  PRIVATE :: ll_begin_distrib, ll_beginp1_distrib, ll_end_distrib, ll_endm1_distrib, ll_endp1_distrib 
+
 
 CONTAINS
 
 
-  SUBROUTINE init_omp_para
+  SUBROUTINE init_omp_para(is_mpi_master)
   USE grid_param
+  USE getin_mod
 #ifdef CPP_USING_OMP
   USE omp_lib
 #endif
   IMPLICIT NONE
+  LOGICAL, INTENT(IN) :: is_mpi_master
   INTEGER :: ll_nb,i
 
@@ -49,58 +90,85 @@
     omp_rank=OMP_GET_THREAD_NUM()
 #endif
-
-    IF (omp_by_domain) THEN
-      omp_first=.TRUE.
-      omp_last=.TRUE.
-      IF (omp_rank==0) THEN
-        omp_master=.TRUE.
-      ELSE
-        omp_master=.FALSE.
-      ENDIF
-     
-      ll_begin=1
-      ll_beginp1=2
-      ll_end=llm
-      ll_endm1=llm-1
-      ll_endp1=llm+1
-     
-    ELSE    
-    
-      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
-    
+    
+    is_omp_master=.FALSE.
+    is_master=.FALSE.
+      
+    IF (omp_rank==0) THEN
+      is_omp_master=.TRUE.
+      IF (is_mpi_master) is_master=.TRUE.
     ENDIF
+
+    omp_level_size=1 
+    CALL getin("omp_level_size",omp_level_size)
+    IF (MOD(omp_size,omp_level_size)/=0) THEN
+      IF (is_mpi_master) PRINT*,"omp_size /= omp_level_size x omp_domain_size  => disable omp threads on vertical layers"
+      omp_level_size=1
+    ENDIF
+    omp_domain_size=omp_size/omp_level_size
+    omp_domain_rank = omp_rank / omp_level_size    
+    omp_level_rank = MOD(omp_rank, omp_level_size)    
+    
+    IF (is_mpi_master) PRINT*,"omp_domain_size",omp_domain_size,"omp_domain_rank", omp_domain_rank
+    IF (is_mpi_master) PRINT*,"omp_level_size",omp_level_size,"omp_level_rank", omp_level_rank
+    
+    is_omp_first_level=.FALSE.
+    is_omp_last_level= .FALSE.
+    is_omp_domain_master=.FALSE.
+    is_omp_level_master=.FALSE.
+
+    IF (omp_domain_rank==0) is_omp_domain_master = .TRUE.
+    IF (omp_level_rank==0) is_omp_level_master = .TRUE.
+    IF (omp_level_rank==0) is_omp_first_level=.TRUE.
+    
+    IF (omp_level_rank==omp_level_size-1) is_omp_last_level=.TRUE.
+    
+    ll_end=0
+    DO i=0,omp_level_rank
+      ll_begin=ll_end+1
+      ll_nb=llm/omp_level_size
+      IF (MOD(llm,omp_level_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 (is_omp_first_level) ll_beginp1=ll_begin+1
+    IF (is_omp_last_level) ll_endp1=ll_endp1+1
+    IF (is_omp_last_level) ll_endm1=ll_endm1-1
+    
+    
+    
+    is_omp_first_level_distrib = is_omp_first_level
+    is_omp_last_level_distrib  = is_omp_last_level
+    ll_begin_distrib           = ll_begin
+    ll_beginp1_distrib         = ll_beginp1
+    ll_end_distrib             = ll_end
+    ll_endm1_distrib           = ll_endm1
+    ll_endp1_distrib           = ll_endp1
+    
+    is_omp_first_level_full = .TRUE.
+    is_omp_last_level_full  = .TRUE.
+    ll_begin_full           = 1
+    ll_beginp1_full         = 2
+    ll_end_full             = llm
+    ll_endm1_full           = llm-1
+    ll_endp1_full           = llm+1    
+
 !$OMP END PARALLEL
 
    ELSE
      omp_size=1
+     omp_level_size=1
+     omp_domain_size=1
      omp_rank=0
-     omp_first=.TRUE.
-     omp_last=.TRUE.
-     omp_master=.TRUE.
+     omp_level_rank=0
+     omp_domain_rank=0
+     is_omp_first_level=.TRUE.
+     is_omp_last_level=.TRUE.
+     is_omp_master=.TRUE.
+     is_omp_domain_master=.TRUE.
+     is_omp_level_master=.TRUE.
      ll_begin=1
      ll_beginp1=2
@@ -108,8 +176,53 @@
      ll_endm1=llm-1
      ll_endp1=llm+1
+     
+     is_omp_first_level_distrib = is_omp_first_level
+     is_omp_last_level_distrib  = is_omp_last_level
+     ll_begin_distrib           = ll_begin
+     ll_beginp1_distrib         = ll_beginp1
+     ll_end_distrib             = ll_end
+     ll_endm1_distrib           = ll_endm1
+     ll_endp1_distrib           = ll_endp1
+    
+     is_omp_first_level_full = .TRUE.
+     is_omp_last_level_full  = .TRUE.
+     ll_begin_full           = 1
+     ll_beginp1_full         = 2
+     ll_end_full             = llm
+     ll_endm1_full           = llm-1
+     ll_endp1_full           = llm+1    
+
    ENDIF
    
   END SUBROUTINE init_omp_para
 
+
+  SUBROUTINE switch_omp_distrib_level
+  IMPLICIT NONE
+     is_omp_first_level = is_omp_first_level_distrib 
+     is_omp_last_level  = is_omp_last_level_distrib 
+     ll_begin           = ll_begin_distrib 
+     ll_beginp1         = ll_beginp1_distrib 
+     ll_end             = ll_end_distrib 
+     ll_endm1           = ll_endm1_distrib 
+     ll_endp1           = ll_endp1_distrib 
+  
+  END SUBROUTINE switch_omp_distrib_level
+  
+
+  SUBROUTINE switch_omp_no_distrib_level
+  IMPLICIT NONE
+
+     is_omp_first_level = is_omp_first_level_full 
+     is_omp_last_level  = is_omp_last_level_full 
+     ll_begin           = ll_begin_full 
+     ll_beginp1         = ll_beginp1_full 
+     ll_end             = ll_end_full 
+     ll_endm1           = ll_endm1_full 
+     ll_endp1           = ll_endp1_full 
+  
+  END SUBROUTINE switch_omp_no_distrib_level
+  
+
   FUNCTION omp_in_parallel()
 #ifdef CPP_USING_OMP
Index: codes/icosagcm/trunk/src/physics.f90
===================================================================
--- codes/icosagcm/trunk/src/physics.f90	(revision 294)
+++ codes/icosagcm/trunk/src/physics.f90	(revision 295)
@@ -10,4 +10,5 @@
   TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)
   TYPE(t_field),POINTER :: f_dulon(:), f_dulat(:)
+  TYPE(t_field),POINTER :: f_temp(:)
 
   CHARACTER(LEN=255) :: physics_type
@@ -38,4 +39,5 @@
        CALL allocate_field(f_dulon,field_t,type_real,llm, name='dulon')
        CALL allocate_field(f_dulat,field_t,type_real,llm, name='dulat')
+       CALL allocate_field(f_temp,field_t,type_real,llm, name='temp')
        CALL init_pack_before ! Compute physics_inout%ngrid and offsets used by pack/unpack
        CALL init_physics_dcmip
@@ -102,4 +104,5 @@
     USE physics_interface_mod
     USE physics_dcmip_mod, ONLY : full_physics_dcmip => full_physics
+    USE theta2theta_rhodz_mod
     USE mpipara
     IMPLICIT NONE
@@ -111,4 +114,5 @@
     REAL(rstd),POINTER :: phis(:)
     REAL(rstd),POINTER :: ps(:)
+    REAL(rstd),POINTER :: temp(:,:)
     REAL(rstd),POINTER :: theta_rhodz(:,:)
     REAL(rstd),POINTER :: ue(:,:)
@@ -118,4 +122,6 @@
     INTEGER :: it, ind
 
+    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
+
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -124,8 +130,8 @@
        phis=f_phis(ind)
        ps=f_ps(ind)
-       theta_rhodz=f_theta_rhodz(ind)
+       temp=f_temp(ind)
        ue=f_ue(ind)
        q=f_q(ind)
-       CALL pack_physics(pack_info(ind), phis, ps, theta_rhodz, ue, q)
+       CALL pack_physics(pack_info(ind), phis, ps, temp, ue, q)
     END DO
 
@@ -134,5 +140,5 @@
        CALL full_physics_dcmip
     CASE DEFAULT
-       IF(is_mpi_root) PRINT *,'Internal error : illegal value of phys_type', phys_type
+       IF(is_mpi_master) PRINT *,'Internal error : illegal value of phys_type', phys_type
        STOP
     END SELECT
@@ -143,10 +149,11 @@
        CALL swap_geometry(ind)
        ps=f_ps(ind)
-       theta_rhodz=f_theta_rhodz(ind)
+       temp=f_temp(ind)
        q=f_q(ind)
        dulon=f_dulon(ind)
        dulat=f_dulat(ind)
-       CALL unpack_physics(pack_info(ind), ps, theta_rhodz, q, dulon, dulat)
+       CALL unpack_physics(pack_info(ind), ps, temp, q, dulon, dulat)
     END DO
+    CALL temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
 
     ! Transfer dulon, dulat
@@ -166,5 +173,5 @@
   END SUBROUTINE physics_column
 
-  SUBROUTINE pack_physics(info, phis, ps, theta_rhodz, ue, q)
+  SUBROUTINE pack_physics(info, phis, ps, temp, ue, q)
     USE icosa
     USE wind_mod
@@ -176,16 +183,16 @@
     REAL(rstd) :: phis(iim*jjm)
     REAL(rstd) :: ps(iim*jjm)
-    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
+    REAL(rstd) :: temp(iim*jjm,llm)
     REAL(rstd) :: ue(3*iim*jjm,llm)
     REAL(rstd) :: q(iim*jjm,llm,nqtot)
 
     REAL(rstd) :: p(iim*jjm,llm+1)
-    REAL(rstd) :: Temp(iim*jjm,llm)
     REAL(rstd) :: uc(iim*jjm,3,llm)
     REAL(rstd) :: ulon(iim*jjm,llm)
     REAL(rstd) :: ulat(iim*jjm,llm)
 
+!$OMP BARRIER
     CALL compute_pression(ps,p,0)
-    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
+!$OMP BARRIER
     CALL compute_wind_centered(ue,uc)
     CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
@@ -199,5 +206,5 @@
   END SUBROUTINE pack_physics
 
-  SUBROUTINE unpack_physics(info, ps,theta_rhodz, q, dulon, dulat)
+  SUBROUTINE unpack_physics(info, ps,temp, q, dulon, dulat)
     USE icosa
     USE physics_interface_mod
@@ -206,6 +213,5 @@
     TYPE(t_pack_info) :: info
     REAL(rstd) :: ps(iim*jjm)
-    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd) :: Temp(iim*jjm,llm)
+    REAL(rstd) :: temp(iim*jjm,llm)
     REAL(rstd) :: q(iim*jjm,llm,nqtot)
     REAL(rstd) :: dulon(iim*jjm,llm)
@@ -221,5 +227,5 @@
     q = q + physics_inout%dt_phys * dq
     Temp = Temp + physics_inout%dt_phys * dTemp
-    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
+!    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
   END SUBROUTINE unpack_physics
 
Index: codes/icosagcm/trunk/src/pression.f90
===================================================================
--- codes/icosagcm/trunk/src/pression.f90	(revision 294)
+++ codes/icosagcm/trunk/src/pression.f90	(revision 295)
@@ -13,4 +13,5 @@
     INTEGER :: ind
 
+!$OMP BARRIER
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -21,4 +22,5 @@
       CALL compute_pression(ps, p,0)
     ENDDO
+!$OMP BARRIER
   
   END SUBROUTINE pression
@@ -27,4 +29,5 @@
   USE icosa
   USE disvert_mod
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: ps(iim*jjm)
@@ -34,13 +37,15 @@
 
     IF(ap_bp_present) THEN
-    DO    l    = 1, llm+1
-      DO j=jj_begin-offset,jj_end+offset
-        DO i=ii_begin-offset,ii_end+offset
-          ij=(j-1)*iim+i
-          p(ij,l) = ap(l) + bp(l) * ps(ij)
+      DO    l    = ll_begin, ll_endp1
+!      DO    l    = 1, llm + 1
+        DO j=jj_begin-offset,jj_end+offset
+          DO i=ii_begin-offset,ii_end+offset
+            ij=(j-1)*iim+i
+            p(ij,l) = ap(l) + bp(l) * ps(ij)
+          ENDDO
         ENDDO
       ENDDO
-    ENDDO
     END IF
+
   END SUBROUTINE compute_pression
 
Index: codes/icosagcm/trunk/src/theta_rhodz.f90
===================================================================
--- codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 294)
+++ codes/icosagcm/trunk/src/theta_rhodz.f90	(revision 295)
@@ -1,6 +1,25 @@
 MODULE theta2theta_rhodz_mod
+  USE field_mod
+  
+  TYPE(t_field), POINTER, SAVE  :: f_p(:)
+  TYPE(t_field), POINTER, SAVE  :: f_pks(:)
+  TYPE(t_field), POINTER, SAVE  :: f_pk(:)
+
+  PRIVATE :: f_p,f_pk,f_pks 
 
 CONTAINS
   
+  SUBROUTINE init_theta2theta_rhodz
+  USE icosa
+  USE field_mod
+  IMPLICIT NONE
+    CALL allocate_field(f_p,field_t,type_real,llm+1,name='p')
+    CALL allocate_field(f_pk,field_t,type_real,llm,name='pk')
+    CALL allocate_field(f_pks,field_t,type_real,name='pks')
+    
+  END SUBROUTINE init_theta2theta_rhodz
+
+
+
   SUBROUTINE theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
   USE icosa
@@ -15,4 +34,5 @@
     INTEGER :: ind
 
+!$OMP BARRIER
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -24,4 +44,5 @@
       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
     ENDDO
+!$OMP BARRIER
   
   END SUBROUTINE theta_rhodz2theta
@@ -29,4 +50,6 @@
   SUBROUTINE theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
   USE icosa
+  USE pression_mod
+  USE exner_mod
   IMPLICIT NONE
     TYPE(t_field), POINTER :: f_ps(:)
@@ -37,17 +60,71 @@
     REAL(rstd), POINTER :: theta_rhodz(:,:)
     REAL(rstd), POINTER :: temp(:,:)
-    INTEGER :: ind
-
-    DO ind=1,ndomain
-      IF (.NOT. assigned_domain(ind)) CYCLE
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      ps=f_ps(ind)
+    REAL(rstd), POINTER :: p(:)
+    REAL(rstd), POINTER :: pk(:,:)
+    REAL(rstd), POINTER :: pks(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      p=f_p(ind)
+      pks=f_pks(ind)
+      pk=f_pk(ind)
       theta_rhodz=f_theta_rhodz(ind)
       temp=f_temp(ind)
-      CALL compute_theta_rhodz2temperature(ps, theta_rhodz,temp,0)
-    ENDDO
+
+!$OMP BARRIER
+      CALL compute_pression(ps,p,0)
+!$OMP BARRIER
+      CALL compute_exner(ps,p,pks,pk,0)
+!$OMP BARRIER
+      CALL compute_theta_rhodz2temperature(p, pk, theta_rhodz,temp,0)
+    ENDDO
+!$OMP BARRIER
   
   END SUBROUTINE theta_rhodz2temperature
+ 
+  SUBROUTINE temperature2theta_rhodz(f_ps,f_temp,f_theta_rhodz)
+  USE icosa
+  USE pression_mod
+  USE exner_mod
+  IMPLICIT NONE
+    TYPE(t_field), POINTER :: f_ps(:)
+    TYPE(t_field), POINTER :: f_theta_rhodz(:)
+    TYPE(t_field), POINTER :: f_temp(:)
+  
+    REAL(rstd), POINTER :: ps(:)
+    REAL(rstd), POINTER :: theta_rhodz(:,:)
+    REAL(rstd), POINTER :: temp(:,:)
+    REAL(rstd), POINTER :: p(:)
+    REAL(rstd), POINTER :: pk(:,:)
+    REAL(rstd), POINTER :: pks(:,:)
+    INTEGER :: ind
+
+    DO ind=1,ndomain
+      IF (.NOT. assigned_domain(ind)) CYCLE
+      CALL swap_dimensions(ind)
+      CALL swap_geometry(ind)
+      ps=f_ps(ind)
+      p=f_p(ind)
+      pks=f_pks(ind)
+      pk=f_pk(ind)
+      theta_rhodz=f_theta_rhodz(ind)
+      temp=f_temp(ind)
+
+!$OMP BARRIER
+      CALL compute_pression(ps,p,0)
+!$OMP BARRIER
+      CALL compute_exner(ps,p,pks,pk,0)
+!$OMP BARRIER
+      CALL compute_temperature2theta_rhodz(p, pk, temp, theta_rhodz, 0)
+    ENDDO
+!$OMP BARRIER
+  
+  END SUBROUTINE temperature2theta_rhodz
+ 
+ 
     
   SUBROUTINE theta2theta_rhodz(f_ps,f_theta,f_theta_rhodz)
@@ -63,4 +140,5 @@
     INTEGER :: ind
 
+!$OMP BARRIER
     DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
@@ -72,4 +150,5 @@
       CALL compute_theta2theta_rhodz(ps, theta, theta_rhodz,0)
     ENDDO
+!$OMP BARRIER
   
   END SUBROUTINE theta2theta_rhodz
@@ -77,5 +156,6 @@
   SUBROUTINE compute_theta2theta_rhodz(ps,theta, theta_rhodz,offset)
   USE icosa
-  USE pression_mod
+  USE disvert_mod
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: ps(iim*jjm)
@@ -84,20 +164,16 @@
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
-    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)
-!$OMP THREADPRIVATE(p)
-
-    ALLOCATE( p(iim*jjm,llm+1))
-    CALL compute_pression(ps,p,offset)
-    
-    DO    l    = 1, llm
-      DO j=jj_begin-offset,jj_end+offset
-        DO i=ii_begin-offset,ii_end+offset
-          ij=(j-1)*iim+i
-          theta_rhodz(ij,l) = theta(ij,l) * (p(ij,l)-p(ij,l+1))/g
-        ENDDO
-      ENDDO
-    ENDDO
-
-    DEALLOCATE( p)
+    
+!$OMP BARRIER
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          theta_rhodz(ij,l) = theta(ij,l) * ( (ap(l)-ap(l+1)) + ( bp(l)- bp(l+1))* ps(ij) )/g
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP BARRIER
+
 
   END SUBROUTINE compute_theta2theta_rhodz
@@ -105,5 +181,6 @@
   SUBROUTINE compute_theta_rhodz2theta(ps,theta_rhodz,theta,offset)
   USE icosa
-  USE pression_mod
+  USE disvert_mod
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: ps(iim*jjm)
@@ -112,42 +189,40 @@
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
-    REAL(rstd),SAVE,ALLOCATABLE :: p(:,:)
-!$OMP THREADPRIVATE(p)
-
-    ALLOCATE( p(iim*jjm,llm+1))
-
-    CALL compute_pression(ps,p,offset)
-    
-    DO    l    = 1, llm
-      DO j=jj_begin-offset,jj_end+offset
-        DO i=ii_begin-offset,ii_end+offset
-          ij=(j-1)*iim+i
-          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g)
-        ENDDO
-      ENDDO
-    ENDDO
-
-    DEALLOCATE( p)
+
+!$OMP BARRIER
+    DO    l    = ll_begin, ll_end
+      DO j=jj_begin-offset,jj_end+offset
+        DO i=ii_begin-offset,ii_end+offset
+          ij=(j-1)*iim+i
+          theta(ij,l) = theta_rhodz(ij,l) / ( (ap(l)-ap(l+1)) + ( bp(l)- bp(l+1))* ps(ij) )/g
+        ENDDO
+      ENDDO
+    ENDDO
+!$OMP BARRIER
+
     
   END SUBROUTINE compute_theta_rhodz2theta
 
-  SUBROUTINE compute_theta_rhodz2temperature(ps,theta_rhodz,temp,offset)
-  USE icosa
-  USE pression_mod
-  USE exner_mod
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
+
+
+
+
+
+  SUBROUTINE compute_theta_rhodz2temperature(p,pk,theta_rhodz,temp,offset)
+  USE icosa
+  USE pression_mod
+  USE exner_mod
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN) :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm)
     INTEGER,INTENT(IN) :: offset
     INTEGER :: i,j,ij,l
-    REAL(rstd) :: p(iim*jjm,llm+1)
-    REAL(rstd) :: pk(iim*jjm,llm)
-    REAL(rstd) :: pks(iim*jjm)
-
-    CALL compute_pression(ps,p,offset)
-    CALL compute_exner(ps,p,pks,pk,offset)
         
-    DO    l    = 1, llm
+! flush p
+!$OMP BARRIER
+    DO    l    = ll_begin, ll_end
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
@@ -157,25 +232,27 @@
       ENDDO
     ENDDO
+!$OMP BARRIER
+    
     
   END SUBROUTINE compute_theta_rhodz2temperature
 
-  SUBROUTINE compute_temperature2theta_rhodz(ps,temp,theta_rhodz,offset)
-  USE icosa
-  USE pression_mod
-  USE exner_mod
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
+  SUBROUTINE compute_temperature2theta_rhodz(p,pk,temp,theta_rhodz,offset)
+  USE icosa
+  USE pression_mod
+  USE exner_mod
+  USE omp_para
+  IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: p(iim*jjm,llm+1)
+    REAL(rstd),INTENT(IN)  :: pk(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
-    REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm)
-    INTEGER,INTENT(IN) :: offset
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: p(iim*jjm,llm+1)
-    REAL(rstd) :: pk(iim*jjm,llm)
-    REAL(rstd) :: pks(iim*jjm)
-
-    CALL compute_pression(ps,p,offset)
-    CALL compute_exner(ps,p,pks,pk,offset)
+    REAL(rstd),INTENT(IN)  :: temp(iim*jjm,llm)
+    INTEGER,INTENT(IN) :: offset
+    INTEGER :: i,j,ij,l
+
         
-    DO    l    = 1, llm
+! flush p
+!$OMP BARRIER
+
+    DO    l    = ll_begin, ll_end
       DO j=jj_begin-offset,jj_end+offset
         DO i=ii_begin-offset,ii_end+offset
@@ -185,4 +262,5 @@
       ENDDO
     ENDDO
+!$OMP BARRIER
     
   END SUBROUTINE compute_temperature2theta_rhodz
Index: codes/icosagcm/trunk/src/time.f90
===================================================================
--- codes/icosagcm/trunk/src/time.f90	(revision 294)
+++ codes/icosagcm/trunk/src/time.f90	(revision 295)
@@ -19,4 +19,7 @@
   INTEGER,SAVE    :: itau_out, itau_adv, itau_dissip, itau_physics, itaumax
 !$OMP THREADPRIVATE(itau_out, itau_adv, itau_dissip, itau_physics, itaumax)  
+
+  INTEGER,SAVE    :: itau_check_conserv
+!$OMP THREADPRIVATE(itau_check_conserv)  
   
   INTEGER,SAVE :: day_step,ndays
@@ -38,7 +41,8 @@
 
   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, & 
-day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref,an, mois, jour,heure, &
-            calend,time_style,itau0
+         dt, write_period, itau_out, itau_adv, itau_dissip, itau_physics, itaumax,  &
+         itau_check_conserv,  & 
+         day_step,ndays,jD_ref,jH_ref,day_ini,day_end,annee_ref,day_ref, &
+         an, mois, jour,heure, calend,time_style,itau0
 
 
@@ -82,11 +86,14 @@
     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)
-    
+
+    itau_check_conserv=HUGE(itau_check_conserv)
+    CALL getin('itau_check_conserv',itau_check_conserv)
+
     IF (is_mpi_root)  THEN
        PRINT *, 'itaumax=',itaumax
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 294)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 295)
@@ -39,4 +39,5 @@
   USE output_field_mod
   USE write_field
+  USE theta2theta_rhodz_mod
   IMPLICIT NONE
 
@@ -117,5 +118,5 @@
     END SELECT
 
-
+    CALL init_theta2theta_rhodz
     CALL init_dissip
     CALL init_caldyn
@@ -124,4 +125,5 @@
     CALL init_check_conserve
     CALL init_physics
+    
 
     CALL etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
@@ -157,4 +159,5 @@
   USE output_field_mod
   USE write_etat0_mod
+  USE checksum_mod
   IMPLICIT NONE  
     REAL(rstd),POINTER :: q(:,:,:)
@@ -172,4 +175,5 @@
     INTEGER :: stop_clock
     INTEGER :: rate_clock
+    INTEGER :: l
     
     
@@ -178,7 +182,8 @@
 !    CALL write_restart(f_ps,f_mass,f_phis,f_theta_rhodz,f_u,f_q) 
     
+    CALL switch_omp_distrib_level
     CALL caldyn_BC(f_phis, f_wflux) ! set constant values in first/last interfaces
-
-!!$OMP BARRIER
+    
+!$OMP BARRIER
   DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
@@ -189,7 +194,10 @@
         CALL compute_rhodz(.TRUE., ps, rhodz) ! save rhodz for transport scheme before dynamics update ps
      ELSE
-        rhodz(:,:)=mass(:,:)
+       DO l=ll_begin,ll_end
+         rhodz(:,l)=mass(:,l)
+       ENDDO
      END IF
   END DO
+!$OMP BARRIER
   fluxt_zero=.TRUE.
 
@@ -224,11 +232,9 @@
     ENDIF
 
-!$OMP MASTER    
-    IF (is_mpi_root) PRINT *,"It No :",It,"   t :",dt*It
-!$OMP END MASTER    
+    IF (is_master) PRINT *,"It No :",It,"   t :",dt*It
+
     IF (mod(it,itau_out)==0 ) THEN
       CALL update_time_counter(dt*it)
       CALL output_field("q",f_q)
-      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
     ENDIF
     
@@ -256,4 +262,6 @@
        
        IF(caldyn_eta==eta_mass) THEN
+!ym flush ps
+!$OMP BARRIER
           DO ind=1,ndomain
              IF (.NOT. assigned_domain(ind)) CYCLE
@@ -267,4 +275,5 @@
 !       CALL wait_message(req_mass)  
        CALL dissip(f_u,f_du,f_mass,f_phis, f_theta_rhodz,f_dtheta_rhodz)
+
 !       CALL send_message(f_mass,req_mass)
 !       CALL wait_message(req_mass)  
@@ -291,4 +300,8 @@
 
     CALL physics(it,f_phis, f_ps, f_theta_rhodz, f_u,  f_q)
+
+    IF (MOD(it,itau_check_conserv)==0) THEN
+      CALL check_conserve(f_ps,f_dps,f_u,f_theta_rhodz,f_phis,it)  
+    ENDIF
     
   ENDDO
@@ -324,5 +337,5 @@
           IF(caldyn_eta==eta_mass) THEN ! update ps
              ps=f_ps(ind) ; dps=f_dps(ind) ;              
-             IF (omp_first) THEN
+             IF (is_omp_first_level) THEN
 !$SIMD
                 DO ij=ij_begin,ij_end
@@ -376,5 +389,5 @@
       ! if mass coordinate, deal with ps first on one core
       IF(caldyn_eta==eta_mass) THEN
-         IF (omp_first) THEN
+         IF (is_omp_first_level) THEN
 
             DO ind=1,ndomain
Index: codes/icosagcm/trunk/src/transfert_mpi.f90
===================================================================
--- codes/icosagcm/trunk/src/transfert_mpi.f90	(revision 294)
+++ codes/icosagcm/trunk/src/transfert_mpi.f90	(revision 295)
@@ -1119,5 +1119,5 @@
 
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           
           rval2d=>field(ind)%rval2d
@@ -1152,5 +1152,5 @@
         
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           rval2d=>field(ind)%rval2d
           req=>message%request(ind)        
@@ -1192,5 +1192,5 @@
       
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
 
           dim3=size(field(ind)%rval3d,2)
@@ -1226,5 +1226,5 @@
          
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           dim3=size(field(ind)%rval3d,2)
           rval3d=>field(ind)%rval3d
@@ -1263,5 +1263,5 @@
     
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE
 
           dim3=size(field(ind)%rval4d,2)
@@ -1302,5 +1302,5 @@
         
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           
           dim3=size(field(ind)%rval4d,2)
@@ -1421,5 +1421,5 @@
         
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
           
           rval2d=>field(ind)%rval2d
@@ -1452,5 +1452,5 @@
         
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
 
           rval3d=>field(ind)%rval3d
@@ -1485,5 +1485,5 @@
                 
         DO ind=1,ndomain
-          IF (.NOT. assigned_domain(ind)) CYCLE
+          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
 
           rval4d=>field(ind)%rval4d
Index: codes/icosagcm/trunk/src/transfert_omp.f90
===================================================================
--- codes/icosagcm/trunk/src/transfert_omp.f90	(revision 294)
+++ codes/icosagcm/trunk/src/transfert_omp.f90	(revision 295)
@@ -34,5 +34,17 @@
   END INTERFACE 
 
-  PUBLIC bcast_omp, reduce_sum_omp, allreduce_sum_omp
+
+  INTERFACE reduce_max_omp
+    MODULE PROCEDURE reduce_max_omp_i,reduce_max_omp_i1,reduce_max_omp_i2,reduce_max_omp_i3,reduce_max_omp_i4, &
+                     reduce_max_omp_r,reduce_max_omp_r1,reduce_max_omp_r2,reduce_max_omp_r3,reduce_max_omp_r4
+  END INTERFACE 
+
+  INTERFACE allreduce_max_omp
+    MODULE PROCEDURE allreduce_max_omp_i,allreduce_max_omp_i1,allreduce_max_omp_i2,allreduce_max_omp_i3,allreduce_max_omp_i4, &
+                     allreduce_max_omp_r,allreduce_max_omp_r1,allreduce_max_omp_r2,allreduce_max_omp_r3,allreduce_max_omp_r4
+  END INTERFACE 
+
+
+  PUBLIC bcast_omp, reduce_sum_omp, allreduce_sum_omp, reduce_max_omp, allreduce_max_omp
   
 CONTAINS
@@ -521,4 +533,268 @@
   END SUBROUTINE allreduce_sum_omp_r4  
   
+
+
+
+
+
+
+
+
+
+
+
+  SUBROUTINE reduce_max_omp_i(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN)  :: VarIn
+    INTEGER,INTENT(OUT) :: VarOut
+    INTEGER             :: VarIn_tmp(1)
+    INTEGER             :: VarOut_tmp(1)
+    
+    VarIn_tmp(1)=VarIn
+    CALL Check_buffer_i(1)   
+    CALL reduce_max_omp_igen(VarIn_tmp,Varout_tmp,1,buffer_i)
+    VarOut=VarOut_tmp(1)
+    
+  END SUBROUTINE reduce_max_omp_i
+
+  SUBROUTINE reduce_max_omp_i1(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:) :: VarOut
+    
+    CALL Check_buffer_i(size(VarIn))   
+    CALL reduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+   
+  END SUBROUTINE reduce_max_omp_i1
+  
+  
+  SUBROUTINE reduce_max_omp_i2(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:) :: VarOut
+
+    CALL Check_buffer_i(size(VarIn))   
+    CALL reduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE reduce_max_omp_i2
+
+
+  SUBROUTINE reduce_max_omp_i3(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:,:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:,:) :: VarOut
+    
+    CALL Check_buffer_i(size(VarIn))   
+    CALL reduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE reduce_max_omp_i3
+
+
+  SUBROUTINE reduce_max_omp_i4(VarIn, VarOut)
+    IMPLICIT NONE
+
+    INTEGER,INTENT(IN),DIMENSION(:,:,:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:,:,:) :: VarOut
+  
+    CALL Check_buffer_i(size(VarIn))   
+    CALL reduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE reduce_max_omp_i4
+
+
+  SUBROUTINE reduce_max_omp_r(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN)  :: VarIn
+    REAL,INTENT(OUT) :: VarOut
+    REAL             :: VarIn_tmp(1)
+    REAL             :: VarOut_tmp(1)
+    
+    VarIn_tmp(1)=VarIn
+    CALL Check_buffer_r(1)   
+    CALL reduce_max_omp_rgen(VarIn_tmp,Varout_tmp,1,buffer_r)
+    VarOut=VarOut_tmp(1)
+  
+  END SUBROUTINE reduce_max_omp_r
+
+  SUBROUTINE reduce_max_omp_r1(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL reduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+   
+  END SUBROUTINE reduce_max_omp_r1
+  
+  
+  SUBROUTINE reduce_max_omp_r2(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL reduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE reduce_max_omp_r2
+
+
+  SUBROUTINE reduce_max_omp_r3(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:,:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:,:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL reduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE reduce_max_omp_r3
+
+
+  SUBROUTINE reduce_max_omp_r4(VarIn, VarOut)
+    IMPLICIT NONE
+
+    REAL,INTENT(IN),DIMENSION(:,:,:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:,:,:) :: VarOut
+  
+    CALL Check_buffer_r(size(VarIn))   
+    CALL reduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE reduce_max_omp_r4  
+  
+  
+  
+  
+    SUBROUTINE allreduce_max_omp_i(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN)  :: VarIn
+    INTEGER,INTENT(OUT) :: VarOut
+    INTEGER             :: VarIn_tmp(1)
+    INTEGER             :: VarOut_tmp(1)
+    
+    VarIn_tmp(1)=VarIn
+    CALL Check_buffer_i(1)   
+    CALL allreduce_max_omp_igen(VarIn_tmp,Varout_tmp,1,buffer_i)
+    VarOut=VarOut_tmp(1)
+    
+  END SUBROUTINE allreduce_max_omp_i
+
+  SUBROUTINE allreduce_max_omp_i1(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:) :: VarOut
+    
+    CALL Check_buffer_i(size(VarIn))   
+    CALL allreduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+   
+  END SUBROUTINE allreduce_max_omp_i1
+  
+  
+  SUBROUTINE allreduce_max_omp_i2(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:) :: VarOut
+
+    CALL Check_buffer_i(size(VarIn))   
+    CALL allreduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE allreduce_max_omp_i2
+
+
+  SUBROUTINE allreduce_max_omp_i3(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    INTEGER,INTENT(IN),DIMENSION(:,:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:,:) :: VarOut
+    
+    CALL Check_buffer_i(size(VarIn))   
+    CALL allreduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE allreduce_max_omp_i3
+
+
+  SUBROUTINE allreduce_max_omp_i4(VarIn, VarOut)
+    IMPLICIT NONE
+
+    INTEGER,INTENT(IN),DIMENSION(:,:,:,:)  :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(:,:,:,:) :: VarOut
+  
+    CALL Check_buffer_i(size(VarIn))   
+    CALL allreduce_max_omp_igen(VarIn,Varout,Size(VarIn),buffer_i)
+  
+  END SUBROUTINE allreduce_max_omp_i4
+
+
+  SUBROUTINE allreduce_max_omp_r(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN)  :: VarIn
+    REAL,INTENT(OUT) :: VarOut
+    REAL             :: VarIn_tmp(1)
+    REAL             :: VarOut_tmp(1)
+    
+    VarIn_tmp(1)=VarIn
+    CALL Check_buffer_r(1)   
+    CALL allreduce_max_omp_rgen(VarIn_tmp,Varout_tmp,1,buffer_r)
+    VarOut=VarOut_tmp(1)
+  
+  END SUBROUTINE allreduce_max_omp_r
+
+  SUBROUTINE allreduce_max_omp_r1(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL allreduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+   
+  END SUBROUTINE allreduce_max_omp_r1
+  
+  
+  SUBROUTINE allreduce_max_omp_r2(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL allreduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE allreduce_max_omp_r2
+
+
+  SUBROUTINE allreduce_max_omp_r3(VarIn, VarOut)
+    IMPLICIT NONE
+  
+    REAL,INTENT(IN),DIMENSION(:,:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:,:) :: VarOut
+    
+    CALL Check_buffer_r(size(VarIn))   
+    CALL allreduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE allreduce_max_omp_r3
+
+
+  SUBROUTINE allreduce_max_omp_r4(VarIn, VarOut)
+    IMPLICIT NONE
+
+    REAL,INTENT(IN),DIMENSION(:,:,:,:)  :: VarIn
+    REAL,INTENT(OUT),DIMENSION(:,:,:,:) :: VarOut
+  
+    CALL Check_buffer_r(size(VarIn))   
+    CALL allreduce_max_omp_rgen(VarIn,Varout,Size(VarIn),buffer_r)
+  
+  END SUBROUTINE allreduce_max_omp_r4  
+
   
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -740,4 +1016,141 @@
   
   END SUBROUTINE allreduce_sum_omp_rgen
+
+
+
+
+
+
+  SUBROUTINE reduce_max_omp_igen(VarIn,VarOut,dimsize,Buff)
+  IMPLICIT NONE
+
+    INTEGER,INTENT(IN) :: dimsize
+    INTEGER,INTENT(IN),DIMENSION(dimsize) :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(dimsize) :: VarOut
+    INTEGER,INTENT(INOUT),DIMENSION(dimsize) :: Buff
+
+    INTEGER :: i
+
+  !$OMP MASTER
+    Buff(:)=VarIn(1)
+  !$OMP END MASTER
+  !$OMP BARRIER
+  
+  !$OMP CRITICAL     
+    DO i=1,dimsize
+      Buff(i)=MAX(Buff(i),VarIn(i))
+    ENDDO
+  !$OMP END CRITICAL
+  !$OMP BARRIER  
+  
+  !$OMP MASTER
+    DO i=1,dimsize
+      VarOut(i)=Buff(i)
+    ENDDO
+  !$OMP END MASTER
+  !$OMP BARRIER
+  
+  END SUBROUTINE reduce_max_omp_igen
+
+  SUBROUTINE reduce_max_omp_rgen(VarIn,VarOut,dimsize,Buff)
+  IMPLICIT NONE
+
+    INTEGER,INTENT(IN) :: dimsize
+    REAL,INTENT(IN),DIMENSION(dimsize) :: VarIn
+    REAL,INTENT(OUT),DIMENSION(dimsize) :: VarOut
+    REAL,INTENT(INOUT),DIMENSION(dimsize) :: Buff
+
+    INTEGER :: i
+
+  !$OMP MASTER
+    Buff(:)=VarIn(1)
+  !$OMP END MASTER
+  !$OMP BARRIER
+  
+  !$OMP CRITICAL     
+    DO i=1,dimsize
+      Buff(i)=MAX(Buff(i),VarIn(i))
+    ENDDO
+  !$OMP END CRITICAL
+  !$OMP BARRIER  
+  
+    DO i=1,dimsize
+      VarOut(i)=Buff(i)
+    ENDDO
+  !$OMP BARRIER
+  
+  END SUBROUTINE reduce_max_omp_rgen
+
+
+
+  SUBROUTINE allreduce_max_omp_igen(VarIn,VarOut,dimsize,Buff)
+  IMPLICIT NONE
+
+    INTEGER,INTENT(IN) :: dimsize
+    INTEGER,INTENT(IN),DIMENSION(dimsize) :: VarIn
+    INTEGER,INTENT(OUT),DIMENSION(dimsize) :: VarOut
+    INTEGER,INTENT(INOUT),DIMENSION(dimsize) :: Buff
+
+    INTEGER :: i
+
+  !$OMP MASTER
+    Buff(:)=VarIn(1)
+  !$OMP END MASTER
+  !$OMP BARRIER
+  
+  !$OMP CRITICAL     
+    DO i=1,dimsize
+      Buff(i)=MAX(Buff(i),VarIn(i))
+    ENDDO
+  !$OMP END CRITICAL
+  !$OMP BARRIER  
+  
+    DO i=1,dimsize
+      VarOut(i)=Buff(i)
+    ENDDO
+  !$OMP BARRIER
+  
+  END SUBROUTINE allreduce_max_omp_igen
+
+  SUBROUTINE allreduce_max_omp_rgen(VarIn,VarOut,dimsize,Buff)
+  IMPLICIT NONE
+
+    INTEGER,INTENT(IN) :: dimsize
+    REAL,INTENT(IN),DIMENSION(dimsize) :: VarIn
+    REAL,INTENT(OUT),DIMENSION(dimsize) :: VarOut
+    REAL,INTENT(INOUT),DIMENSION(dimsize) :: Buff
+
+    INTEGER :: i
+
+  !$OMP MASTER
+    Buff(:)=VarIn(1)
+  !$OMP END MASTER
+  !$OMP BARRIER
+  
+  !$OMP CRITICAL     
+    DO i=1,dimsize
+      Buff(i)=MAX(Buff(i),VarIn(i))
+    ENDDO
+  !$OMP END CRITICAL
+  !$OMP BARRIER  
+  
+    DO i=1,dimsize
+      VarOut(i)=Buff(i)
+    ENDDO
+
+  !$OMP BARRIER
+  
+  END SUBROUTINE allreduce_max_omp_rgen
+    
+
+
+
+
+
+
+
+
+
+
     
 END MODULE transfert_omp_mod
Index: codes/icosagcm/trunk/src/vertical_interp.f90
===================================================================
--- codes/icosagcm/trunk/src/vertical_interp.f90	(revision 294)
+++ codes/icosagcm/trunk/src/vertical_interp.f90	(revision 295)
@@ -21,4 +21,5 @@
   USE icosa
   USE pression_mod
+  USE omp_para
   IMPLICIT NONE
     TYPE(t_field),POINTER :: f_ps(:)
@@ -48,4 +49,5 @@
 
   SUBROUTINE compute_vertical_interp(p,in,out,pval)
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
@@ -56,19 +58,25 @@
     INTEGER :: i,j,ij,l
         
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
-        l=llm-1
-        DO WHILE(0.5*(p(ij,l)+p(ij,l+1))<pval .AND. l>1)
-          l=l-1
+!$OMP BARRIER    
+    IF (is_omp_level_master) THEN
+    
+      DO j=jj_begin-1,jj_end+1
+        DO i=ii_begin-1,ii_end+1
+          ij=(j-1)*iim+i
+          l=llm-1
+          DO WHILE(0.5*(p(ij,l)+p(ij,l+1))<pval .AND. l>1)
+            l=l-1
+          ENDDO
+          pmid=0.5*(p(ij,l)+p(ij,l+1))
+          pmidp1=0.5*(p(ij,l+1)+p(ij,l+2))
+
+          coeff=(pval-pmid)/(pmid-pmidp1)
+        
+          out(ij)=in(ij,l)+coeff*(in(ij,l)-in(ij,l+1))
         ENDDO
-        pmid=0.5*(p(ij,l)+p(ij,l+1))
-        pmidp1=0.5*(p(ij,l+1)+p(ij,l+2))
-
-        coeff=(pval-pmid)/(pmid-pmidp1)
-        
-        out(ij)=in(ij,l)+coeff*(in(ij,l)-in(ij,l+1))
       ENDDO
-    ENDDO
+   
+    ENDIF
+!$OMP BARRIER    
         
   END SUBROUTINE compute_vertical_interp
Index: codes/icosagcm/trunk/src/vorticity.f90
===================================================================
--- codes/icosagcm/trunk/src/vorticity.f90	(revision 294)
+++ codes/icosagcm/trunk/src/vorticity.f90	(revision 295)
@@ -29,4 +29,5 @@
   USE icosa
   USE disvert_mod
+  USE omp_para
   IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: ue(3*iim*jjm,llm)
@@ -34,5 +35,5 @@
     INTEGER :: i,j,ij,l
 
-    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
Index: codes/icosagcm/trunk/src/wind.f90
===================================================================
--- codes/icosagcm/trunk/src/wind.f90	(revision 294)
+++ codes/icosagcm/trunk/src/wind.f90	(revision 295)
@@ -47,5 +47,5 @@
   SUBROUTINE compute_wind_centered(ue,ucenter)
   USE icosa
-  
+  USE omp_para
   IMPLICIT NONE
   REAL(rstd) :: ue(3*iim*jjm,llm)
@@ -53,5 +53,5 @@
   INTEGER :: i,j,ij,l    
  
-    DO l=1,llm
+    DO l=ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -73,4 +73,5 @@
  SUBROUTINE compute_wind_on_edge(ue,uedge)
   USE icosa
+  USE omp_para
     
   IMPLICIT NONE
@@ -83,5 +84,5 @@
     CALL compute_tangential_compound(ue,ut)
   
-    DO l=1,llm
+    DO l=ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -100,4 +101,5 @@
  SUBROUTINE compute_tangential_compound(ue,ut)
   USE icosa  
+  USE omp_para
   IMPLICIT NONE
   REAL(rstd) :: ue(3*iim*jjm,llm)
@@ -105,5 +107,5 @@
   INTEGER :: i,j,l,ij
     
-  DO l=1,llm
+  DO l=ll_begin,ll_end
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
@@ -155,4 +157,5 @@
  SUBROUTINE compute_wind_lonlat_compound(u, ulon, ulat)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -164,5 +167,5 @@
     
   
-    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
@@ -184,4 +187,5 @@
   SUBROUTINE compute_wind_from_lonlat_compound(ulon, ulat, u)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -192,5 +196,5 @@
   INTEGER :: i,j,ij,l     
   
-    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
@@ -207,4 +211,5 @@
   SUBROUTINE compute_wind_centered_from_lonlat_compound(ulon, ulat, u)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -214,5 +219,5 @@
 
   INTEGER :: i,j,ij,l     
-  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
@@ -248,4 +253,5 @@
   SUBROUTINE compute_wind_perp_from_lonlat_compound(ulon, ulat, up)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -259,5 +265,5 @@
    CALL compute_wind_from_lonlat_compound(ulon, ulat, u)
 
-    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
@@ -297,4 +303,5 @@
  SUBROUTINE compute_wind_centered_lonlat_compound(uc, ulon, ulat)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -306,5 +313,5 @@
     
   
-    DO l=1,llm
+    DO l=ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -320,4 +327,5 @@
  SUBROUTINE compute_wind_centered_from_wind_lonlat_centered(ulon, ulat,uc)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -329,5 +337,5 @@
     
   
-    DO l=1,llm
+    DO l=ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
@@ -344,4 +352,5 @@
  SUBROUTINE compute_wind_perp_from_wind_centered(uc,un)
   USE icosa  
+  USE omp_para
     
   IMPLICIT NONE
@@ -352,5 +361,5 @@
     
   
-    DO l=1,llm
+    DO l=ll_begin,ll_end
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
