Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 136)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 137)
@@ -114,6 +114,6 @@
   END SUBROUTINE compute_gradq3d
 
-  !===================================================================================================	 
-  SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,he,bigt)
+  ! Backward trajectories, used in Miura approach
+  SUBROUTINE compute_backward_traj(normal,tangent,ue,tau, cc)
     USE domain_mod
     USE dimensions
@@ -121,28 +121,16 @@
     USE metric
     IMPLICIT NONE
-    LOGICAL, INTENT(IN)      :: update_mass
     REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
     REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
-    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
-    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm)
     REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
-    REAL(rstd),INTENT(IN)    :: he(3*iim*jjm,llm) ! mass flux
-    REAL(rstd),INTENT(IN)    :: bigt
-
-    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 
-    REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
-    REAL(rstd) :: v_e(3*iim*jjm,llm,3)
-    REAL(rstd) :: up_e
-
-    REAL(rstd) :: qe(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT)   :: cc(3*iim*jjm,llm,3) ! start of backward trajectory
+    REAL(rstd),INTENT(IN)    :: tau
+
+    REAL(rstd) :: v_e(3), up_e
+    REAL(rstd) :: qe
     REAL(rstd) :: ed(3)    
-    REAL(rstd) :: flxx(3*iim*jjm,llm)
-    INTEGER :: i,j,n,k,l
-    REAL(rstd):: him_old(llm) 
-
-
-    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2
-    ! TODO : this does not depend on q and should be done only once, not for each tracer
+    INTEGER :: i,j,n,l
+
+    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
     DO l = 1,llm
        DO j=jj_begin-1,jj_end+1
@@ -162,4 +150,91 @@
                   wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         &
                   )
+             v_e = ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)
+             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e*tau
+
+             up_e=1/de(n+u_lup)*(                                                      &
+                  wee(n+u_lup,1,1)*le(n+u_left)*ue(n+u_left,l)+                        &
+                  wee(n+u_lup,2,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                      &
+                  wee(n+u_lup,3,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                      &
+                  wee(n+u_lup,4,1)*le(n+u_right)*ue(n+u_right,l)+                      &
+                  wee(n+u_lup,5,1)*le(n+u_rup)*ue(n+u_rup,l)+                          & 
+                  wee(n+u_lup,1,2)*le(n+t_lup+u_right)*ue(n+t_lup+u_right,l)+          &
+                  wee(n+u_lup,2,2)*le(n+t_lup+u_rup)*ue(n+t_lup+u_rup,l)+              &
+                  wee(n+u_lup,3,2)*le(n+t_lup+u_lup)*ue(n+t_lup+u_lup,l)+              &
+                  wee(n+u_lup,4,2)*le(n+t_lup+u_left)*ue(n+t_lup+u_left,l)+            &
+                  wee(n+u_lup,5,2)*le(n+t_lup+u_ldown)*ue(n+t_lup+u_ldown,l)           &
+                  )
+             v_e = ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
+             cc(n+u_lup,l,:) = xyz_e(n+u_lup,:) - v_e*tau
+
+             up_e=1/de(n+u_ldown)*(                                                    &
+                  wee(n+u_ldown,1,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    &
+                  wee(n+u_ldown,2,1)*le(n+u_right)*ue(n+u_right,l)+                    &
+                  wee(n+u_ldown,3,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
+                  wee(n+u_ldown,4,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
+                  wee(n+u_ldown,5,1)*le(n+u_left)*ue(n+u_left,l)+                      & 
+                  wee(n+u_ldown,1,2)*le(n+t_ldown+u_lup)*ue(n+t_ldown+u_lup,l)+        &
+                  wee(n+u_ldown,2,2)*le(n+t_ldown+u_left)*ue(n+t_ldown+u_left,l)+      &
+                  wee(n+u_ldown,3,2)*le(n+t_ldown+u_ldown)*ue(n+t_ldown+u_ldown,l)+    &
+                  wee(n+u_ldown,4,2)*le(n+t_ldown+u_rdown)*ue(n+t_ldown+u_rdown,l)+    &
+                  wee(n+u_ldown,5,2)*le(n+t_ldown+u_right)*ue(n+t_ldown+u_right,l)     &
+                  )
+             v_e = ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
+             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e*tau
+          ENDDO
+       ENDDO
+    END DO
+  END SUBROUTINE compute_backward_traj
+
+  ! Horizontal transport (S. Dubey, T. Dubos)
+  ! Slope-limited van Leer approach with hexagons
+  SUBROUTINE compute_advect_horiz(update_mass, normal,tangent,qi,gradq3d,him,ue,hfluxt,bigt)
+    USE domain_mod
+    USE dimensions
+    USE geometry
+    USE metric
+    IMPLICIT NONE
+    LOGICAL, INTENT(IN)      :: update_mass
+    REAL(rstd),INTENT(IN)    :: normal(3*iim*jjm,3)
+    REAL(rstd),INTENT(IN)    :: tangent(3*iim*jjm,3)
+    REAL(rstd),INTENT(INOUT) :: qi(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: gradq3d(iim*jjm,llm,3) 
+    REAL(rstd),INTENT(INOUT) :: him(iim*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: ue(iim*3*jjm,llm)
+    REAL(rstd),INTENT(IN)    :: hfluxt(3*iim*jjm,llm) ! mass flux
+    REAL(rstd),INTENT(IN)    :: bigt
+
+    REAL(rstd) :: dqi(iim*jjm,llm),dhi(iim*jjm,llm) 
+    REAL(rstd) :: cc(3*iim*jjm,llm,3)   ! barycenter of quadrilateral, where q is evaluated (1-point quadrature)
+    REAL(rstd) :: v_e(3*iim*jjm,llm,3)
+    REAL(rstd) :: up_e
+
+!    REAL(rstd) :: qe(3*iim*jjm,llm)
+    REAL(rstd) :: qe
+    REAL(rstd) :: ed(3)    
+    REAL(rstd) :: flxx(3*iim*jjm,llm)
+    INTEGER :: i,j,n,k,l
+    REAL(rstd):: him_old(llm) 
+
+
+    ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*dt/2
+    ! TODO : this does not depend on q and should be done only once, not for each tracer
+    DO l = 1,llm
+       DO j=jj_begin-1,jj_end+1
+          DO i=ii_begin-1,ii_end+1
+             n=(j-1)*iim+i
+
+             up_e =1/de(n+u_right)*(                                                   &
+                  wee(n+u_right,1,1)*le(n+u_rup)*ue(n+u_rup,l)+                        &
+                  wee(n+u_right,2,1)*le(n+u_lup)*ue(n+u_lup,l)+                        &
+                  wee(n+u_right,3,1)*le(n+u_left)*ue(n+u_left,l)+                      &
+                  wee(n+u_right,4,1)*le(n+u_ldown)*ue(n+u_ldown,l)+                    &
+                  wee(n+u_right,5,1)*le(n+u_rdown)*ue(n+u_rdown,l)+                    & 
+                  wee(n+u_right,1,2)*le(n+t_right+u_ldown)*ue(n+t_right+u_ldown,l)+    &
+                  wee(n+u_right,2,2)*le(n+t_right+u_rdown)*ue(n+t_right+u_rdown,l)+    &
+                  wee(n+u_right,3,2)*le(n+t_right+u_right)*ue(n+t_right+u_right,l)+    &
+                  wee(n+u_right,4,2)*le(n+t_right+u_rup)*ue(n+t_right+u_rup,l)+        &
+                  wee(n+u_right,5,2)*le(n+t_right+u_lup)*ue(n+t_right+u_lup,l)         &
+                  )
              v_e(n+u_right,l,:)= ue(n+u_right,l)*normal(n+u_right,:) + up_e*tangent(n+u_right,:)
 
@@ -208,9 +283,11 @@
              if (ne(n,right)*ue(n+u_right,l)>0) then 
                 ed = cc(n+u_right,l,:) - xyz_i(n,:)
-                qe(n+u_right,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
+                qe=qi(n,l)+sum2(gradq3d(n,l,:),ed) 
              else
                 ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
-                qe(n+u_right,l)=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
+                qe=qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
              endif
+             flxx(n+u_right,l) = hfluxt(n+u_right,:)*qe*ne(n,right)
+             
              if (ne(n,lup)*ue(n+u_lup,l)>0) then 
                 ed = cc(n+u_lup,l,:) - xyz_i(n,:)
@@ -220,4 +297,6 @@
                 qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
              endif
+             flxx(n+u_lup,:) = hfluxt(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup) 
+             
              if (ne(n,ldown)*ue(n+u_ldown,l)>0) then 
                 ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
@@ -227,4 +306,5 @@
                 qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
              endif
+             flxx(n+u_ldown,:) = hfluxt(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)
           end do
        end do
@@ -237,7 +317,4 @@
        DO i=ii_begin-1,ii_end+1
           n=(j-1)*iim+i 
-          flxx(n+u_right,:) = he(n+u_right,:)*qe(n+u_right,:)*ne(n,right)
-          flxx(n+u_lup,:) = he(n+u_lup,:)*qe(n+u_lup,:)*ne(n,lup) 
-          flxx(n+u_ldown,:) = he(n+u_ldown,:)*qe(n+u_ldown,:)*ne(n,ldown)
        ENDDO
     ENDDO
