Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 174)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 186)
@@ -49,11 +49,13 @@
   !=======================================================================================
 
-  SUBROUTINE compute_gradq3d(qi,one_over_sqrt_leng,gradq3d)
+  SUBROUTINE compute_gradq3d(qi_in,one_over_sqrt_leng_in,gradq3d_out,xyz_i,xyz_v)
     USE trace
     USE omp_para
     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),INTENT(IN)  :: qi_in(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)  :: one_over_sqrt_leng_in(iim*jjm)
+    REAL(rstd),INTENT(IN)  :: xyz_i(iim*jjm,3)
+    REAL(rstd),INTENT(IN)  :: xyz_v(2*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: gradq3d_out(iim*jjm,llm,3) 
     REAL(rstd) :: maxq,minq,minq_c,maxq_c 
     REAL(rstd) :: alphamx,alphami,alpha	,maggrd
@@ -63,6 +65,16 @@
     REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
     INTEGER :: ij,k,ind,l
-
-    CALL trace_start("compute_gradq3d")
+    REAL(rstd)  :: qi(iim*jjm,llm)
+    REAL(rstd)  :: one_over_sqrt_leng(iim*jjm)
+    REAL(rstd) :: gradq3d(iim*jjm,llm,3) 
+    REAL(rstd) :: detx,dety,detz,det
+    REAL(rstd) :: A(3,3), a11,a12,a13,a21,a22,a23,a31,a32,a33
+    REAL(rstd) :: x1,x2,x3
+    REAL(rstd) :: dq(3)
+
+    qi=qi_in
+    one_over_sqrt_leng=one_over_sqrt_leng_in
+    
+    CALL trace_start("compute_gradq3d1")
 
     ! TODO : precompute ar, drop arr as output argument of gradq ?
@@ -71,16 +83,152 @@
     ! Compute gradient at triangles solving a linear system
     ! arr = area of triangle joining centroids of hexagons
+!     DO l = ll_begin,ll_end 
+!!$SIMD
+!      DO ij=ij_begin_ext,ij_end_ext
+!!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
+!!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
+!             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
+!             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
+!       END DO
+!    END DO
+
      DO l = ll_begin,ll_end 
 !$SIMD
       DO ij=ij_begin_ext,ij_end_ext
-             CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,:),arr(ij+z_up))
-             CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,:),arr(ij+z_down))
-       END DO
+!       CALL gradq(ij,l,ij+t_rup,ij+t_lup,ij+z_up,qi,gradtri(ij+z_up,l,1),gradtri(ij+z_up,l,2),gradtri(ij+z_up,l,3),arr(ij+z_up))
+
+        
+        A(1,1)=xyz_i(ij+t_rup,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_rup,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_rup,3)-xyz_i(ij,3) 
+        A(2,1)=xyz_i(ij+t_lup,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_lup,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_lup,3)-xyz_i(ij,3) 
+        A(3,1)=xyz_v(ij+z_up,1);               A(3,2)= xyz_v(ij+z_up,2);             A(3,3)=xyz_v(ij+z_up,3)
+    
+        dq(1) = qi(ij+t_rup,l)-qi(ij,l)
+        dq(2) = qi(ij+t_lup,l)-qi(ij,l)
+        dq(3) = 0.0
+
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         det =  x1 - x2 + x3
+                 
+!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
+
+         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detx =  x1 - x2 + x3
+        
+!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         dety =  x1 - x2 + x3
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detz =  x1 - x2 + x3
+
+        gradtri(ij+z_up,l,1) = detx
+        gradtri(ij+z_up,l,2) = dety
+        gradtri(ij+z_up,l,3) = detz
+        arr(ij+z_up) = det
+        
+      ENDDO
+      
+      DO ij=ij_begin_ext,ij_end_ext
+
+
+!        CALL gradq(ij,l,ij+t_ldown,ij+t_rdown,ij+z_down,qi,gradtri(ij+z_down,l,1),gradtri(ij+z_down,l,2),gradtri(ij+z_down,l,3),arr(ij+z_down))
+
+        A(1,1)=xyz_i(ij+t_ldown,1)-xyz_i(ij,1);  A(1,2)=xyz_i(ij+t_ldown,2)-xyz_i(ij,2); A(1,3)=xyz_i(ij+t_ldown,3)-xyz_i(ij,3) 
+        A(2,1)=xyz_i(ij+t_rdown,1)-xyz_i(ij,1);  A(2,2)=xyz_i(ij+t_rdown,2)-xyz_i(ij,2); A(2,3)=xyz_i(ij+t_rdown,3)-xyz_i(ij,3) 
+        A(3,1)=xyz_v(ij+z_down,1);               A(3,2)= xyz_v(ij+z_down,2);             A(3,3)=xyz_v(ij+z_down,3)
+ 
+        dq(1) = qi(ij+t_ldown,l)-qi(ij,l)
+        dq(2) = qi(ij+t_rdown,l)-qi(ij,l)
+        dq(3) = 0.0
+
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         det =  x1 - x2 + x3
+                 
+!        CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),detx)
+
+         a11=dq(1)  ; a12=dq(2)  ; a13=dq(3)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detx =  x1 - x2 + x3
+        
+!        CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dety)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=dq(1)  ; a22=dq(2)  ; a23=dq(3)
+         a31=A(1,3) ; a32=A(2,3) ; a33=A(3,3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         dety =  x1 - x2 + x3
+
+!        CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),detz)
+
+         a11=A(1,1) ; a12=A(2,1) ; a13=A(3,1)
+         a21=A(1,2) ; a22=A(2,2) ; a23=A(3,2)
+         a31=dq(1)  ; a32=dq(2)  ; a33=dq(3)
+
+         x1 =  a11 * (a22 * a33 - a23 * a32)
+         x2 =  a12 * (a21 * a33 - a23 * a31)
+         x3 =  a13 * (a21 * a32 - a22 * a31)
+         detz =  x1 - x2 + x3
+
+         gradtri(ij+z_down,l,1) = detx
+         gradtri(ij+z_down,l,2) = dety
+         gradtri(ij+z_down,l,3) = detz
+         arr(ij+z_down) = det
+
+      END DO
     END DO
 
 !$SIMD
-      DO ij=ij_begin,ij_end
-         ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
+    DO ij=ij_begin,ij_end
+       ar(ij) = arr(ij+z_up)+arr(ij+z_lup)+arr(ij+z_ldown)+arr(ij+z_down)+arr(ij+z_rdown)+arr(ij+z_rup)+1.e-50
     ENDDO
+    CALL trace_end("compute_gradq3d1")
+    CALL trace_start2("compute_gradq3d2")
       
     DO k=1,3
@@ -94,4 +242,6 @@
       END DO
     ENDDO
+    CALL trace_end2("compute_gradq3d2")
+    CALL trace_start("compute_gradq3d3")
 
     !============================================================================================= LIMITING 
@@ -99,5 +249,6 @@
 !$SIMD
       DO ij=ij_begin,ij_end
-             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
+!             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
+             maggrd = gradq3d(ij,l,1)*gradq3d(ij,l,1) + gradq3d(ij,l,2)*gradq3d(ij,l,2) + gradq3d(ij,l,3)*gradq3d(ij,l,3) 
              maggrd = sqrt(maggrd) 
              maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
@@ -112,9 +263,81 @@
              alphami = max(alphami,0.0) 
              alpha   = min(alphamx,alphami,1.0)
-             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
+!             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
+             gradq3d(ij,l,1) = alpha*gradq3d(ij,l,1)
+             gradq3d(ij,l,2) = alpha*gradq3d(ij,l,2)
+             gradq3d(ij,l,3) = alpha*gradq3d(ij,l,3)
        END DO
     END DO
 
-  CALL trace_end("compute_gradq3d")
+  CALL trace_end("compute_gradq3d3")
+  
+  gradq3d_out=gradq3d
+  
+  CONTAINS
+
+    SUBROUTINE gradq(n0,l,n1,n2,n3,q,dq1,dq2,dq3,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
+    REAL(rstd), INTENT(OUT)    :: dq1,dq2,dq3,det
+    REAL(rstd)    :: dq(3)
+    
+    REAL(rstd)  ::detx,dety,detz
+    INTEGER    :: info
+    INTEGER    :: IPIV(3)
+
+    REAL(rstd) :: A(3,3)
+    REAL(rstd) :: B(3)
+    
+    ! TODO : replace A by A1,A2,A3
+    A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3) 
+    A(2,1)=xyz_i(n2,1)-xyz_i(n0,1);  A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3) 
+    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,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)
+!    CALL determinant(A(:,1),A(:,2),A(:,3),det)
+!    CALL determinant(dq,A(:,2),A(:,3),detx)
+!    CALL determinant(A(:,1),dq,A(:,3),dety)
+!    CALL determinant(A(:,1),A(:,2),dq,detz)
+!    dq(1) = detx
+!    dq(2) = dety
+!    dq(3) = detz 
+
+    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),det)
+    CALL determinant(dq(1),dq(2),dq(3),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),dq1)
+    CALL determinant(A(1,1),A(2,1),A(3,1),dq(1),dq(2),dq(3),A(1,3),A(2,3),A(3,3),dq2)
+    CALL determinant(A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),dq(1),dq(2),dq(3),dq3)
+
+  END SUBROUTINE gradq
+
+  !==========================================================================
+!  PURE SUBROUTINE determinant(a1,a2,a3,det)
+!    IMPLICIT NONE
+!    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
+!    REAL(rstd), INTENT(OUT) :: det
+!    REAL(rstd) ::  x1,x2,x3
+!    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
+!    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
+!    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
+!    det =  x1 - x2 + x3
+!  END SUBROUTINE determinant
+
+   SUBROUTINE determinant(a11,a12,a13,a21,a22,a23,a31,a32,a33,det)
+    IMPLICIT NONE
+    REAL(rstd), INTENT(IN) :: a11,a12,a13,a21,a22,a23,a31,a32,a33
+    REAL(rstd), INTENT(OUT) :: det
+    REAL(rstd) ::  x1,x2,x3
+    x1 =  a11 * (a22 * a33 - a23 * a32)
+    x2 =  a12 * (a21 * a33 - a23 * a31)
+    x3 =  a13 * (a21 * a32 - a22 * a31)
+    det =  x1 - x2 + x3
+  END SUBROUTINE determinant
+
+
   END SUBROUTINE compute_gradq3d
 
@@ -222,8 +445,10 @@
              IF (ne(ij,right)*hfluxt(ij+u_right,l)>0) THEN 
                 ed = cc(ij+u_right,l,:) - xyz_i(ij,:)
-                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
+!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
+                 qe = qi(ij,l)+gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
              ELSE
                 ed = cc(ij+u_right,l,:) - xyz_i(ij+t_right,:)
-                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 
+!                qe = qi(ij+t_right,l)+sum2(gradq3d(ij+t_right,l,:),ed) 
+                qe = qi(ij+t_right,l) + gradq3d(ij+t_right,l,1)*ed(1)+gradq3d(ij+t_right,l,2)*ed(2)+gradq3d(ij+t_right,l,3)*ed(3) 
              ENDIF
              qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
@@ -231,8 +456,10 @@
              IF (ne(ij,lup)*hfluxt(ij+u_lup,l)>0) THEN 
                 ed = cc(ij+u_lup,l,:) - xyz_i(ij,:)
-                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
+!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed)
+                qe = qi(ij,l)  + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
              ELSE
                 ed = cc(ij+u_lup,l,:) - xyz_i(ij+t_lup,:)
-                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 
+!                qe = qi(ij+t_lup,l)+sum2(gradq3d(ij+t_lup,l,:),ed) 
+                qe = qi(ij+t_lup,l) + gradq3d(ij+t_lup,l,1)*ed(1)+gradq3d(ij+t_lup,l,2)*ed(2)+gradq3d(ij+t_lup,l,3)*ed(3) 
              ENDIF
              qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
@@ -240,8 +467,10 @@
              IF (ne(ij,ldown)*hfluxt(ij+u_ldown,l)>0) THEN 
                 ed = cc(ij+u_ldown,l,:) - xyz_i(ij,:)
-                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
+!                qe = qi(ij,l)+sum2(gradq3d(ij,l,:),ed) 
+                qe = qi(ij,l) + gradq3d(ij,l,1)*ed(1)+gradq3d(ij,l,2)*ed(2)+gradq3d(ij,l,3)*ed(3) 
              ELSE
                 ed = cc(ij+u_ldown,l,:) - xyz_i(ij+t_ldown,:)
-                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 
+!                qe = qi(ij+t_ldown,l)+sum2(gradq3d(ij+t_ldown,l,:),ed) 
+                qe = qi(ij+t_ldown,l)+gradq3d(ij+t_ldown,l,1)*ed(1)+gradq3d(ij+t_ldown,l,2)*ed(2)+gradq3d(ij+t_ldown,l,3)*ed(3) 
              ENDIF
              qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
@@ -289,47 +518,4 @@
   END SUBROUTINE compute_advect_horiz
 
-  !==========================================================================
-  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
-    
-    REAL(rstd)  ::detx,dety,detz
-    INTEGER    :: info
-    INTEGER    :: IPIV(3)
-
-    REAL(rstd) :: A(3,3)
-    REAL(rstd) :: B(3)
-
-    ! TODO : replace A by A1,A2,A3
-    A(1,1)=xyz_i(n1,1)-xyz_i(n0,1);  A(1,2)=xyz_i(n1,2)-xyz_i(n0,2); A(1,3)=xyz_i(n1,3)-xyz_i(n0,3) 
-    A(2,1)=xyz_i(n2,1)-xyz_i(n0,1);  A(2,2)=xyz_i(n2,2)-xyz_i(n0,2); A(2,3)=xyz_i(n2,3)-xyz_i(n0,3) 
-    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,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)
-    CALL determinant(A(:,1),A(:,2),A(:,3),det)
-    CALL determinant(dq,A(:,2),A(:,3),detx)
-    CALL determinant(A(:,1),dq,A(:,3),dety)
-    CALL determinant(A(:,1),A(:,2),dq,detz)
-    dq(1) = detx
-    dq(2) = dety
-    dq(3) = detz 
-  END SUBROUTINE gradq
-
-  !==========================================================================
-  PURE SUBROUTINE determinant(a1,a2,a3,det)
-    IMPLICIT NONE
-    REAL(rstd), DIMENSION(3), INTENT(IN) :: a1,a2,a3
-    REAL(rstd), INTENT(OUT) :: det
-    REAL(rstd) ::  x1,x2,x3
-    x1 =  a1(1) * (a2(2) * a3(3) - a2(3) * a3(2))
-    x2 =  a1(2) * (a2(1) * a3(3) - a2(3) * a3(1))
-    x3 =  a1(3) * (a2(1) * a3(2) - a2(2) * a3(1))
-    det =  x1 - x2 + x3
-  END SUBROUTINE determinant
 
 END MODULE advect_mod
