Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 141)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 148)
@@ -11,8 +11,9 @@
   !==========================================================================
 
-  SUBROUTINE init_advect(normal,tangent)
+  SUBROUTINE init_advect(normal,tangent,one_over_sqrt_leng)
     IMPLICIT NONE
     REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
     REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
 
     INTEGER :: i,j,n
@@ -39,5 +40,9 @@
           tangent(n+u_ldown,:)=xyz_v(n+z_down,:)-xyz_v(n+z_ldown,:)
           tangent(n+u_ldown,:)=tangent(n+u_ldown,:)/sqrt(sum(tangent(n+u_ldown,:)**2)+1e-50)
-       END DO
+
+          one_over_sqrt_leng(n) = 1./sqrt(max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
+                                          sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
+                                          sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2)) )
+       ENDDO
     ENDDO
 
@@ -46,15 +51,19 @@
   !=======================================================================================
 
-  SUBROUTINE compute_gradq3d(qi,gradq3d)
+  SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d)
+    USE trace
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng(iim*jjm)
     REAL(rstd),INTENT(OUT) :: gradq3d(iim*jjm,llm,3) 
     REAL(rstd) :: maxq,minq,minq_c,maxq_c 
-    REAL(rstd) :: alphamx,alphami,alpha	,maggrd,leng
+    REAL(rstd) :: alphamx,alphami,alpha	,maggrd
     REAL(rstd) :: leng1,leng2
     REAL(rstd) :: arr(2*iim*jjm)
+    REAL(rstd) :: ar(iim*jjm)
     REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
-    REAL(rstd) :: ar
     INTEGER :: i,j,n,k,ind,l
+
+    CALL trace_start("compute_gradq3d")
 
     ! TODO : precompute ar, drop arr as output argument of gradq ?
@@ -67,22 +76,29 @@
           DO i=ii_begin-1,ii_end+1
              n=(j-1)*iim+i   
-             CALL gradq(n,n+t_rup,n+t_lup,n+z_up,qi(:,l),gradtri(n+z_up,l,:),arr(n+z_up))
-             CALL gradq(n,n+t_ldown,n+t_rdown,n+z_down,qi(:,l),gradtri(n+z_down,l,:),arr(n+z_down))
+             CALL gradq(n,l,n+t_rup,n+t_lup,n+z_up,qi,gradtri(n+z_up,l,:),arr(n+z_up))
+             CALL gradq(n,l,n+t_ldown,n+t_rdown,n+z_down,qi,gradtri(n+z_down,l,:),arr(n+z_down))
           END DO
        END DO
     END DO
 
-    !      Do l =1,llm
-    DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          gradq3d(n,:,:) = gradtri(n+z_up,:,:) + gradtri(n+z_down,:,:)   +   &
-               gradtri(n+z_rup,:,:) +  gradtri(n+z_ldown,:,:)   +   & 
-               gradtri(n+z_lup,:,:)+    gradtri(n+z_rdown,:,:)  
-          ar = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup)
-          gradq3d(n,:,:) = gradq3d(n,:,:)/(ar+1.e-50) 
-       END DO
-    END DO
-    !	 END DO
+    DO j=jj_begin-1,jj_end+1
+      DO i=ii_begin-1,ii_end+1
+         n=(j-1)*iim+i   
+         ar(n) = arr(n+z_up)+arr(n+z_lup)+arr(n+z_ldown)+arr(n+z_down)+arr(n+z_rdown)+arr(n+z_rup)+1.e-50
+      ENDDO
+    ENDDO
+      
+    DO k=1,3
+      DO l =1,llm
+        DO j=jj_begin,jj_end
+           DO i=ii_begin,ii_end
+              n=(j-1)*iim+i
+              gradq3d(n,l,k) = ( gradtri(n+z_up,l,k) + gradtri(n+z_down,l,k)   +   &
+                                 gradtri(n+z_rup,l,k) +  gradtri(n+z_ldown,l,k)   +   & 
+                                 gradtri(n+z_lup,l,k)+    gradtri(n+z_rdown,l,k) ) / ar(n)  
+           END DO
+        END DO
+      END DO
+    ENDDO
 
     !============================================================================================= LIMITING 
@@ -94,9 +110,11 @@
              maggrd =  dot_product(gradq3d(n,l,:),gradq3d(n,l,:))
              maggrd = sqrt(maggrd) 
-             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
-                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
-                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2))
-             maxq_c = qi(n,l) + maggrd*sqrt(leng)
-             minq_c = qi(n,l) - maggrd*sqrt(leng)
+!             leng = max(sum((xyz_v(n+z_up,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_down,:) - xyz_i(n,:))**2), &
+!                  sum((xyz_v(n+z_rup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_rdown,:) - xyz_i(n,:))**2),  &
+!                  sum((xyz_v(n+z_lup,:) - xyz_i(n,:))**2),sum((xyz_v(n+z_ldown,:) - xyz_i(n,:))**2))
+!             maxq_c = qi(n,l) + maggrd*sqrt(leng(n))
+!             minq_c = qi(n,l) - maggrd*sqrt(leng(n))
+             maxq_c = qi(n,l) + maggrd*one_over_sqrt_leng(n)
+             minq_c = qi(n,l) - maggrd*one_over_sqrt_leng(n)
              maxq = max(qi(n,l),qi(n+t_right,l),qi(n+t_lup,l),qi(n+t_rup,l),qi(n+t_left,l), &
                   qi(n+t_rdown,l),qi(n+t_ldown,l))
@@ -112,8 +130,11 @@
        END DO
     END DO
+
+  CALL trace_end("compute_gradq3d")
   END SUBROUTINE compute_gradq3d
 
   ! Backward trajectories, for use with Miura approach
   SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
+    USE trace
     IMPLICIT NONE
     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
@@ -125,4 +146,6 @@
     REAL(rstd) :: v_e(3), up_e, qe, ed(3)    
     INTEGER :: i,j,n,l
+
+    CALL trace_start("compute_backward_traj")
 
     ! TODO : compute normal displacement ue*tau as hfluxt / mass(upwind) then reconstruct tangential displacement
@@ -181,4 +204,7 @@
        ENDDO
     END DO
+
+    CALL trace_end("compute_backward_traj")
+
   END SUBROUTINE compute_backward_traj
 
@@ -186,4 +212,5 @@
   ! Slope-limited van Leer approach with hexagons
   SUBROUTINE compute_advect_horiz(update_mass,hfluxt,cc,gradq3d, mass,qi) 
+    USE trace
     IMPLICIT NONE
     LOGICAL, INTENT(IN)       :: update_mass
@@ -197,4 +224,6 @@
     REAL(rstd) :: qflux(3*iim*jjm,llm)
     INTEGER :: i,j,n,k,l
+
+    CALL trace_start("compute_advect_horiz")
 
     ! evaluate tracer field at point cc using piecewise linear reconstruction
@@ -262,4 +291,5 @@
     END DO
 
+    CALL trace_end("compute_advect_horiz")
   CONTAINS
 
@@ -275,8 +305,8 @@
 
   !==========================================================================
-  PURE SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
-    IMPLICIT NONE
-    INTEGER, INTENT(IN) :: n0,n1,n2,n3
-    REAL(rstd), INTENT(IN)     :: q(iim*jjm)
+  PURE SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq,det)
+    IMPLICIT NONE
+    INTEGER, INTENT(IN) :: n0,l,n1,n2,n3
+    REAL(rstd), INTENT(IN)     :: q(iim*jjm,llm)
     REAL(rstd), INTENT(OUT)    :: dq(3), det
     
@@ -293,6 +323,6 @@
     A(3,1)=xyz_v(n3,1);              A(3,2)= xyz_v(n3,2);            A(3,3)=xyz_v(n3,3)
 
-    dq(1) = q(n1)-q(n0)
-    dq(2) = q(n2)-q(n0)
+    dq(1) = q(n1,l)-q(n0,l)
+    dq(2) = q(n2,l)-q(n0,l)
     dq(3) = 0.0
     !    CALL DGESV(3,1,A,3,IPIV,dq(:),3,info)
