Index: /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 664)
+++ /codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 665)
@@ -46,5 +46,5 @@
 SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
   FIELD_MASS   :: m_ik, theta   ! IN*2
-  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL*5 
+  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL
   DBL :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff
   INTEGER :: iter
@@ -83,9 +83,9 @@
 END SUBROUTINE compute_NH_geopot
 
-SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
+SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, w_il,berni,gradPhi2,DePhil,v_el,G_el,F_el, hflux,du,dPhi,dW)
   FIELD_U      :: u, hflux, du   ! IN, OUT, OUT
-  FIELD_MASS   :: rhodz, berni   ! IN, LOCAL
-  FIELD_GEOPOT :: Phi,W,dPhi,dW, w_il, gradPhi2  ! IN,IN, OUT,OUT, LOCAL 
-  FIELD_UL     :: DePhil, v_el, G_el, F_el ! LOCAL
+  FIELD_MASS   :: rhodz, berni   ! IN, BUF
+  FIELD_GEOPOT :: Phi,W,dPhi,dW, w_il, gradPhi2  ! IN,IN, OUT,OUT, BUF*2
+  FIELD_UL     :: DePhil, v_el, G_el, F_el ! BUF*4
   DECLARE_INDICES
   DECLARE_EDGES
@@ -96,9 +96,9 @@
 END SUBROUTINE compute_caldyn_slow_NH
 
-SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi,dW,du)
+SUBROUTINE compute_caldyn_solver(tau,rhodz,theta, berni,pres,m_il, pk,geopot,W,dPhi,dW,du)
   DBL, INTENT(IN) :: tau
-  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, LOCAL
+  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, BUF*2
   FIELD_THETA  :: theta                  ! IN
-  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, LOCAL 
+  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, BUF
   FIELD_U      :: du                     ! OUT
   DECLARE_INDICES
@@ -113,8 +113,8 @@
 END SUBROUTINE compute_caldyn_solver
 
-SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
-  FIELD_MASS   :: mass, eta_dot, wcov, W_etadot
-  FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW
-  FIELD_U      :: du
+SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW)
+  FIELD_MASS   :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3
+  FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2
+  FIELD_U      :: du ! INOUT
   DECLARE_INDICES
   DBL :: w_ij, wflux_ij
Index: /codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 664)
+++ /codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 665)
@@ -58,8 +58,10 @@
   !----------------------------      PROFILING      --------------------------
   
-  SUBROUTINE init_trace()
+  SUBROUTINE init_trace() BINDC(init_trace)
+    !$OMP MASTER
     time_spent(:)=0.
     bytes(:)=0
     nb_calls(:)=0
+    !$OMP END MASTER
   END SUBROUTINE init_trace
 
@@ -67,28 +69,36 @@
     INTEGER :: id
     DBL :: total_spent
+    !$OMP MASTER
     total_spent=SUM(time_spent)
-    PRINT *, '========================= Performance metrics ========================='
-    PRINT *, 'Total time spent in instrumented code (seconds) :', total_spent
-    PRINT *, 'Name, #calls, %time, microsec/call, MB/sec'    
-    DO id=1,nb_routines
-       IF(nb_calls(id)>0) PRINT *, id_name(id), nb_calls(id), INT(100.*time_spent(id)/total_spent), &
-            INT(1e6*time_spent(id)/nb_calls(id)), INT(1e-6*bytes(id)/time_spent(id))
-    END DO
-    CALL init_trace()
+    IF(total_spent>.01) THEN
+       PRINT *, '========================= Performance metrics ========================='
+       PRINT *, 'Total time spent in instrumented code (seconds) :', total_spent
+       PRINT *, 'Name, #calls, %time, microsec/call, MB/sec'    
+       DO id=1,nb_routines
+          IF(nb_calls(id)>0) PRINT *, id_name(id), nb_calls(id), INT(100.*time_spent(id)/total_spent), &
+               INT(1e6*time_spent(id)/nb_calls(id)), INT(1e-6*bytes(id)/time_spent(id))
+       END DO
+       CALL init_trace()
+    END IF
+    !$OMP END MASTER
   END SUBROUTINE print_trace
 
   SUBROUTINE enter_trace(id, nbytes)
     INTEGER :: id, nbytes
+    !$OMP MASTER
     current_id = id
     bytes(id) = bytes(id) + nbytes
     nb_calls(id)=nb_calls(id)+1
     start_time = OMP_GET_WTIME()
+    !$OMP END MASTER
   END SUBROUTINE enter_trace
 
   SUBROUTINE exit_trace()
     DBL :: elapsed
+    !$OMP MASTER
     elapsed = OMP_GET_WTIME()-start_time
     IF(elapsed<0.) elapsed=0.
     time_spent(current_id) = time_spent(current_id) + elapsed
+    !$OMP END MASTER
   END SUBROUTINE exit_trace
 
Index: /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 664)
+++ /codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 665)
@@ -45,12 +45,15 @@
     FIELD_GEOPOT :: wflux, w, geopot, &               ! DIAG, INOUT
          dPhi_fast, dPhi_slow, dW_fast, dW_slow       ! OUT
-    FIELD_U      :: u,du_fast,du_slow,hflux,Ftheta,qu ! INOUT,OUT,OUT,DIAG
+    FIELD_U      :: u,du_fast,du_slow,hflux,qu        ! INOUT,OUT,OUT,DIAG
     FIELD_Z      :: qv                                ! DIAG
     FIELD_PS     :: ps,dmass_col,mass_col             ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
-    FIELD_UL     :: wwuu
+    ! buffers for fields that need to be shared among OpenMP threads
+    ! vertical size is allocated as llm+1 even if only llm is needed
+    FIELD_GEOPOT :: bufm1, bufm2, bufm3
+    FIELD_UL     :: bufu1, bufu2, bufu3,bufu4
+
     DBL          :: time1,time2
     INTEGER :: ij
     
-    CALL init_trace()
     !  CALL CPU_TIME(time1)
     time1=OMP_GET_WTIME()
@@ -72,7 +75,7 @@
        CALL compute_pvort_only(rhodz,u,qv,qu)
        CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
-       CALL compute_coriolis(hflux,theta,qu,Ftheta, drhodz,dtheta_rhodz,du_slow)
+       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
        IF(caldyn_eta == eta_mass) THEN
-          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu)
+          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
        END IF
        !$OMP END PARALLEL
@@ -88,12 +91,12 @@
        END DO
        CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
-       CALL compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi_fast,dW_fast,du_fast)
+       CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast)
        CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
        CALL compute_pvort_only(rhodz,u,qv,qu)
-       CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux,du_slow,dPhi_slow,dW_slow) 
-       CALL compute_coriolis(hflux,theta,qu,Ftheta, drhodz,dtheta_rhodz,du_slow)
+       CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow) 
+       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
        IF(caldyn_eta == eta_mass) THEN
-          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu)
-          CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow,dPhi_slow,dW_slow)
+          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
+          CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow)
        END IF
     END IF
@@ -117,6 +120,6 @@
     FIELD_MASS   :: rhodz, pk, berni                 ! IN, DIAG
     FIELD_THETA  :: theta_rhodz, theta               ! IN, DIAG
-    FIELD_U      :: u,hflux,qu,Ftheta                ! INOUT,DIAG*3
-    FIELD_GEOPOT :: geopot, w, wflux, w_il           ! IN, INOUT, DIAG, DIAG
+    FIELD_U      :: u,hflux,qu                       ! INOUT,DIAG*2
+    FIELD_GEOPOT :: geopot, w, wflux                 ! IN, INOUT, DIAG
     FIELD_Z      :: qv                               ! DIAG
     DOUBLE2(       primal_num, max_nb_stage) :: dmass_col    ! OUT
@@ -126,9 +129,11 @@
     DOUBLE3(llm+1, primal_num, max_nb_stage) :: &
          dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT
-    FIELD_UL     :: DePhil, v_el, G_el, wwuu         ! DIAG*4
+    ! buffers for fields that need to be shared among OpenMP threads
+    ! vertical size is allocated as llm+1 even if only llm is needed
+    FIELD_GEOPOT :: bufm1, bufm2, bufm3
+    FIELD_UL     :: bufu1, bufu2, bufu3,bufu4
     DBL       :: time1,time2
     INTEGER :: step, stage, ij
     
-    CALL init_trace()
     !CALL CPU_TIME(time1)
     time1=OMP_GET_WTIME()
@@ -154,8 +159,8 @@
              CALL compute_pvort_only(rhodz,u,qv,qu)
              CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
-             CALL compute_coriolis(hflux,theta,qu,Ftheta, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
+             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
              IF(caldyn_eta == eta_mass) THEN
                 CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
-                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu)
+                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1)
              END IF
              
@@ -168,14 +173,16 @@
              !$OMP END DO
              CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
-             CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, &
+             CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, &
                   dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
              CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
              CALL compute_pvort_only(rhodz,u,qv,qu)
-             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
-             CALL compute_coriolis(hflux,theta,qu,Ftheta, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
+             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux, &
+                  du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
+             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
              IF(caldyn_eta == eta_mass) THEN
                 CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
-                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu)
-                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
+                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1)
+                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, &
+                     du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
              END IF
           END IF ! NH
