Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 21)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 22)
@@ -1,93 +1,58 @@
-      MODULE advect_mod
-   USE icosa
-	IMPLICIT NONE 
-
-
-  TYPE(t_field),POINTER :: f_normal(:)
-  TYPE(t_field),POINTER :: f_tangent(:)
-  TYPE(t_field),POINTER ::  f_gradq3d(:)
-  REAL(rstd),POINTER :: gradq3d(:,:,:) 
-  REAL(rstd),POINTER :: normal(:,:)
-  REAL(rstd),POINTER :: tangent(:,:)
-			
+MODULE advect_mod
+  USE icosa
+  IMPLICIT NONE 
+
 CONTAINS
-	 
-	
-  SUBROUTINE allocate_advect
-  USE field_mod
-  USE domain_mod
-  USE dimensions
-  USE geometry
-  USE metric
-  IMPLICIT NONE
-
-    CALL allocate_field(f_normal,field_u,type_real,3)
-    CALL allocate_field(f_tangent,field_u,type_real,3)
-    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)
-
-  END SUBROUTINE allocate_advect
-!==========================================================================
-  SUBROUTINE swap_advect(ind)
-  USE domain_mod
-  USE dimensions
-  USE geometry
-  USE metric
-  IMPLICIT NONE
-  INTEGER,INTENT(IN) :: ind
-  
-      normal=f_normal(ind)
-      tangent=f_tangent(ind)
-      gradq3d = f_gradq3d(ind) 
-
-  END SUBROUTINE swap_advect
-!==========================================================================
-    
-  SUBROUTINE init_advect
-  USE domain_mod
-  USE dimensions
-  USE geometry
-  USE metric
-  USE vector
-  IMPLICIT NONE
-   INTEGER :: ind,i,j,n
-   
-   CALL allocate_advect
-   
+
+  !==========================================================================
+
+  SUBROUTINE init_advect(normal,tangent)
+    USE domain_mod
+    USE dimensions
+    USE geometry
+    USE metric
+    USE vector
+    IMPLICIT NONE
+    REAL(rstd),INTENT(OUT) :: normal(3*iim*jjm,3)
+    REAL(rstd),INTENT(OUT) :: tangent(3*iim*jjm,3)
+
+    INTEGER :: i,j,n
+
     DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        n=(j-1)*iim+i
-            
-        CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:))
-        normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right)
-         
-        CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:))
-         normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup)
-        
-        CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))	 
-        normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown)
-                      
-        tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:) 
-        tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50)
-             
-        tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:) 
-        tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50)
-     
-        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
-     ENDDO      
+       DO i=ii_begin-1,ii_end+1
+          n=(j-1)*iim+i
+
+          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),normal(n+u_right,:))
+          normal(n+u_right,:)=normal(n+u_right,:)/sqrt(sum(normal(n+u_right,:)**2)+1e-50)*ne(n,right)
+
+          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),normal(n+u_lup,:))
+          normal(n+u_lup,:)=normal(n+u_lup,:)/sqrt(sum(normal(n+u_lup,:)**2)+1e-50)*ne(n,lup)
+
+          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),normal(n+u_ldown,:))	 
+          normal(n+u_ldown,:)=normal(n+u_ldown,:)/sqrt(sum(normal(n+u_ldown,:)**2)+1e-50)*ne(n,ldown)
+
+          tangent(n+u_right,:)=xyz_v(n+z_rup,:)-xyz_v(n+z_rdown,:) 
+          tangent(n+u_right,:)=tangent(n+u_right,:)/sqrt(sum(tangent(n+u_right,:)**2)+1e-50)
+
+          tangent(n+u_lup,:)=xyz_v(n+z_lup,:)-xyz_v(n+z_up,:) 
+          tangent(n+u_lup,:)=tangent(n+u_lup,:)/sqrt(sum(tangent(n+u_lup,:)**2)+1e-50)
+
+          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
+    ENDDO
+
   END SUBROUTINE init_advect
-!=======================================================================================
-
-
-
-!=======================================================================================
-  SUBROUTINE advect1(qi)
-  USE domain_mod
-  USE dimensions
-  USE geometry
-  USE metric
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: qi(iim*jjm,llm)
+
+  !=======================================================================================
+
+  SUBROUTINE compute_gradq3d(qi,gradq3d)
+    USE domain_mod
+    USE dimensions
+    USE geometry
+    USE metric
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN)  :: qi(iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: gradq3d(iim*jm,llm,3) 
     REAL(rstd) :: maxq,minq,minq_c,maxq_c 
     REAL(rstd) :: alphamx,alphami,alpha	,maggrd,leng
@@ -97,72 +62,77 @@
     REAL(rstd) :: ar
     INTEGER :: i,j,n,k,ind,l
-!========================================================================================== GRADIENT 
-       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   
-          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))
-         END DO
-        END DO
-       END DO
-
-!      Do l =1,llm
-        DO j=jj_begin,jj_end
+    !========================================================================================== GRADIENT 
+    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   
+             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))
+          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
+
+    !============================================================================================= LIMITING 
+    ! GO TO 120 
+    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
- 
-!============================================================================================= LIMITING 
-! GO TO 120 
-  Do l =1,llm
-   DO j=jj_begin,jj_end
-     DO i=ii_begin,ii_end
-       n=(j-1)*iim+i
-       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)
-       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))
-       minq = min(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))
-       alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) )
-       alphamx = max(alphamx,0.0)
-       alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l))
-       alphami = max(alphami,0.0) 
-       alpha   = min(alphamx,alphami,1.0)
-       gradq3d(n,l,:) = alpha*gradq3d(n,l,:)
-      END DO
-    END DO
-   END DO
-	 END SUBROUTINE ADVECT1		 
-!===================================================================================================	 
-	 SUBROUTINE advect2(qi,him,ue,he,bigt)
-  USE domain_mod
-  USE dimensions
-  USE geometry
-  USE metric
-  IMPLICIT NONE
+             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)
+             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))
+             minq = min(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))
+             alphamx = (maxq - qi(n,l)) ; alphamx = alphamx/(maxq_c - qi(n,l) )
+             alphamx = max(alphamx,0.0)
+             alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l))
+             alphami = max(alphami,0.0) 
+             alpha   = min(alphamx,alphami,1.0)
+             gradq3d(n,l,:) = alpha*gradq3d(n,l,:)
+          END DO
+       END DO
+    END DO
+  END SUBROUTINE compute_gradq3d
+
+  !===================================================================================================	 
+  SUBROUTINE compute_advect_horiz(normal,tangent,qi,him,ue,he,bigt)
+    USE domain_mod
+    USE dimensions
+    USE geometry
+    USE metric
+    IMPLICIT NONE
+    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*jm,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)
-    REAL(rstd),INTENT(IN)::bigt
+    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) 
     REAL(rstd) :: v_e(3*iim*jjm,llm,3)
     REAL(rstd) :: up_e
-    
+
     REAL(rstd) :: qe(3*iim*jjm,llm)
     REAL(rstd) :: ed(3)    
@@ -172,132 +142,132 @@
 
 
-!go to  123
-   DO l = 1,llm
+    !go to  123
+    DO l = 1,llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             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,:)
+
+             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(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
+
+             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(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
+
+             cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i 
+             cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt 
+             cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 
+          ENDDO
+       ENDDO
+    END DO
+    !123	 continue           
+    !==========================================================================================================
+    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
+             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) 
+             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) 
+             endif
+             if (ne(n,lup)*ue(n+u_lup,l)>0) then 
+                ed = cc(n+u_lup,l,:) - xyz_i(n,:)
+                qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)
+             else
+                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
+                qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
+             endif
+             if (ne(n,ldown)*ue(n+u_ldown,l)>0) then 
+                ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
+                qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed) 
+             else
+                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
+                qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
+             endif
+          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 
+          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
+
     DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
-        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,:)
-  
-                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(n+u_lup,l,:)= ue(n+u_lup,l)*normal(n+u_lup,:) + up_e*tangent(n+u_lup,:)
-    
-                  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(n+u_ldown,l,:)= ue(n+u_ldown,l)*normal(n+u_ldown,:) + up_e*tangent(n+u_ldown,:)
-
-          cc(n+u_right,l,:) = xyz_e(n+u_right,:) - v_e(n+u_right,l,:)*0.5*bigt    !! redge is mid point of edge i 
-          cc(n+u_lup,l,:)   = xyz_e(n+u_lup,:)   - v_e(n+u_lup,l,:)*0.5*bigt 
-          cc(n+u_ldown,l,:) = xyz_e(n+u_ldown,:) - v_e(n+u_ldown,l,:)*0.5*bigt 
-       ENDDO
-     ENDDO  
-    END DO
-!123	 continue           
-!==========================================================================================================
-    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
-	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) 
-	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) 
-	endif
-	if (ne(n,lup)*ue(n+u_lup,l)>0) then 
-	   ed = cc(n+u_lup,l,:) - xyz_i(n,:)
-  	 qe(n+u_lup,l)=qi(n,l)+sum2(gradq3d(n,l,:),ed)
-	else
-	   ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
-	 qe(n+u_lup,l)=qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
-	endif
-	if (ne(n,ldown)*ue(n+u_ldown,l)>0) then 
-	  ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
-	 qe(n+u_ldown,l)=qi(n,l)+ sum2(gradq3d(n,l,:),ed) 
-	else
-	 ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
-	qe(n+u_ldown,l)=qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
-	endif
-	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 
-	    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
-		 		
-     DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-        n=(j-1)*iim+i  
-
-	 dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) &
-                       + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) &
-                       + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) ) 
-
-         dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       & 
-                         +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    &
-                         - flxx(n+u_left,:) - flxx(n+u_rdown,:) )   
-         him_old(:) = him(n,:) 
-	 him(n,:) = him(n,:) + dhi(n,:) 
-         qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)  
-                       
-       END DO
-     END DO
-
-	CONTAINS
-!====================================================================================
-	  REAL FUNCTION sum2(a1,a2)
-        IMPLICIT NONE
-        REAL,INTENT(IN):: a1(3), a2(3)
-		sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
-	   END FUNCTION sum2
- 
-		END SUBROUTINE advect2 
-!==========================================================================
+          n=(j-1)*iim+i  
+
+          dhi(n,:)= -(1/Ai(n))*(he(n+u_right,:)*ne(n,right) + he(n+u_lup,:)*ne(n,lup) &
+               + he(n+u_ldown,:)*ne(n,ldown) + he(n+u_rup,:)*ne(n,rup) &
+               + he(n+u_left,:)*ne(n,left) +  he(n+u_rdown,:)*ne(n,rdown) ) 
+
+          dqi(n,:)= -(1/Ai(n))*(flxx(n+u_right,:)+flxx(n+u_lup,:)       & 
+               +flxx(n+u_ldown,:) - flxx(n+u_rup,:)    &
+               - flxx(n+u_left,:) - flxx(n+u_rdown,:) )   
+          him_old(:) = him(n,:) 
+          him(n,:) = him(n,:) + dhi(n,:) 
+          qi(n,:) = (qi(n,:)*him_old(:) + dqi(n,:))/him(n,:)  
+
+       END DO
+    END DO
+
+  CONTAINS
+    !====================================================================================
+    REAL FUNCTION sum2(a1,a2)
+      IMPLICIT NONE
+      REAL,INTENT(IN):: a1(3), a2(3)
+      sum2 = a1(1)*a2(1)+a1(2)*a2(2)+a1(3)*a2(3)
+    END FUNCTION sum2
+
+  END SUBROUTINE compute_advect_horiz
+  !==========================================================================
   SUBROUTINE gradq(n0,n1,n2,n3,q,dq,det)
-  USE geometry
-  USE metric
-  USE dimensions
-  IMPLICIT NONE
+    USE geometry
+    USE metric
+    USE dimensions
+    IMPLICIT NONE
     INTEGER, INTENT(IN) :: n0,n1,n2,n3
     REAL,INTENT(IN)     :: q(iim*jjm)
@@ -313,27 +283,28 @@
     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)-q(n0)
     dq(2) = q(n2)-q(n0)
     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 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
-!==========================================================================
-	  SUBROUTINE determinant(a1,a2,a3,det)
-        IMPLICIT NONE
-        REAL, DIMENSION(3) :: a1, a2,a3
-        REAL ::  det,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    
-	END MODULE  
+  !==========================================================================
+  SUBROUTINE determinant(a1,a2,a3,det)
+    IMPLICIT NONE
+    REAL, DIMENSION(3) :: a1, a2,a3
+    REAL ::  det,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
Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 21)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 22)
@@ -4,400 +4,422 @@
   INTEGER,PARAMETER::iapp_tracvl= 3
   REAL(rstd),SAVE :: dt
-  
+
+  TYPE(t_field),POINTER :: f_normal(:)
+  TYPE(t_field),POINTER :: f_tangent(:)
+  TYPE(t_field),POINTER :: f_gradq3d(:)
+
   PUBLIC init_advect_tracer, advect_tracer
 
 CONTAINS
-  
+
   SUBROUTINE init_advect_tracer(dt_in)
-  USE advect_mod
-  IMPLICIT NONE
+    USE advect_mod
+    IMPLICIT NONE
     REAL(rstd),INTENT(IN) :: dt_in
-    
+    REAL(rstd),POINTER :: tangent(:,:)
+    REAL(rstd),POINTER :: normal(:,:)
+
     dt=dt_in
-    
-    CALL init_advect
-    
+    CALL allocate_field(f_normal,field_u,type_real,3)
+    CALL allocate_field(f_tangent,field_u,type_real,3)
+    CALL allocate_field(f_gradq3d,field_t,type_real,llm,3)
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       normal=f_normal(ind)
+       tangent=f_tangent(ind)
+       CALL init_advect(normal,tangent)
+    END DO
+
   END SUBROUTINE init_advect_tracer
-  
-  
-  SUBROUTINE advect_tracer(f_ps,f_u, f_q)
-  USE icosa
-  USE advect_mod
-  USE disvert_mod
-  IMPLICIT NONE
-  TYPE(t_field),POINTER :: f_ps(:)
-  TYPE(t_field),POINTER :: f_u(:)
-  TYPE(t_field),POINTER :: f_q(:)
-  REAL(rstd),POINTER :: q(:,:,:)
-  REAL(rstd),POINTER :: u(:,:) 
-  REAL(rstd),POINTER :: ps(:) 
-  REAL(rstd),POINTER :: massflx(:,:)
-  REAL(rstd),POINTER :: rhodz(:,:)
-  TYPE(t_field),POINTER,SAVE :: f_massflxc(:)
-  TYPE(t_field),POINTER,SAVE :: f_massflx(:)
-  TYPE(t_field),POINTER,SAVE :: f_uc(:)
-  TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)
-  TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
-  REAL(rstd),POINTER,SAVE :: massflxc(:,:)
-  REAL(rstd),POINTER,SAVE :: uc(:,:)
-  REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)
-  REAL(rstd):: bigt 
-  INTEGER :: ind,it,iapptrac,i,j,l,ij
-  INTEGER,SAVE :: iadvtr=0
-  LOGICAL,SAVE:: first=.TRUE. 
-!------------------------------------------------------sarvesh
+
+  SUBROUTINE advect_tracer(f_ps,f_u,f_q)
+    USE icosa
+    USE advect_mod
+    USE disvert_mod
+    IMPLICIT NONE
+    TYPE(t_field),POINTER :: f_ps(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_q(:)
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: u(:,:) 
+    REAL(rstd),POINTER :: ps(:) 
+    REAL(rstd),POINTER :: massflx(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:)
+    TYPE(t_field),POINTER,SAVE :: f_massflxc(:)
+    TYPE(t_field),POINTER,SAVE :: f_massflx(:)
+    TYPE(t_field),POINTER,SAVE :: f_uc(:)
+    TYPE(t_field),POINTER,SAVE :: f_rhodzm1(:)
+    TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
+    REAL(rstd),POINTER,SAVE :: massflxc(:,:)
+    REAL(rstd),POINTER,SAVE :: uc(:,:)
+    REAL(rstd),POINTER,SAVE :: rhodzm1(:,:)
+    REAL(rstd):: bigt
+    INTEGER :: ind,it,i,j,l,ij
+    INTEGER,SAVE :: iadvtr=0
+    LOGICAL,SAVE:: first=.TRUE. 
+    !------------------------------------------------------sarvesh
     IF ( first ) THEN 
-      CALL allocate_field(f_rhodz,field_t,type_real,llm) 
-      CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
-      CALL allocate_field(f_massflxc,field_u,type_real,llm) 
-      CALL allocate_field(f_massflx,field_u,type_real,llm) 
-      CALL allocate_field(f_uc,field_u,type_real,llm) 
-      first = .FALSE.
-    END IF 
-
-    DO ind=1,ndomain
-      CALL swap_dimensions(ind)
-      CALL swap_geometry(ind)
-      rhodz=f_rhodz(ind)
-      massflx=f_massflx(ind)
-      ps=f_ps(ind)
-      u=f_u(ind)
-      
-      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) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/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
-            massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
-            massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
-            massflx(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
-    
-    ENDDO
-    
-
-    
-	IF ( iadvtr == 0 ) THEN 
-	 DO ind=1,ndomain         
-         CALL swap_dimensions(ind)
-         CALL swap_geometry(ind)
-	   rhodz=f_rhodz(ind)  
-	   rhodzm1 = f_rhodzm1(ind) 
-	   massflxc = f_massflxc(ind) 
-	   rhodzm1 = rhodz
-	   massflxc = 0.0 
-	   uc = f_uc(ind)
-	   uc = 0.0 
-	 END DO
-         CALL transfert_request(f_rhodzm1,req_i1)   !---->
- 	 CALL transfert_request(f_massflxc,req_e1) !---->
-         CALL transfert_request(f_massflxc,req_e1) !------>
-	 CALL transfert_request(f_uc,req_e1) !---->
-         CALL transfert_request(f_uc,req_e1) 
-        END IF
-
-	iadvtr = iadvtr + 1 
-	iapptrac = iadvtr  
-
-      DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-	    massflx=f_massflx(ind)
-	    rhodzm1 = f_rhodzm1(ind) 
-	    massflxc = f_massflxc(ind) 
-	    massflxc = massflxc + massflx 
-	     uc = f_uc(ind) 
-	     u  = f_u(ind) 
-	     uc = uc + u 
+       CALL allocate_field(f_rhodz,field_t,type_real,llm) 
+       CALL allocate_field(f_rhodzm1,field_t,type_real,llm) 
+       CALL allocate_field(f_massflxc,field_u,type_real,llm) 
+       CALL allocate_field(f_massflx,field_u,type_real,llm) 
+       CALL allocate_field(f_uc,field_u,type_real,llm) 
+       first = .FALSE.
+    END IF
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       rhodz=f_rhodz(ind)
+       massflx=f_massflx(ind)
+       ps=f_ps(ind)
+       u=f_u(ind)
+
+       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) = (ap(l) - ap(l+1) + (bp(l)-bp(l+1))*ps(ij))/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
+                massflx(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
+                massflx(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
+                massflx(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
+    ENDDO
+
+    IF ( iadvtr == 0 ) THEN 
+       DO ind=1,ndomain         
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          rhodz=f_rhodz(ind)  
+          rhodzm1 = f_rhodzm1(ind) 
+          massflxc = f_massflxc(ind)  ! accumulated mass fluxes
+          uc = f_uc(ind)              ! time-integrated normal velocity to compute back-trajectory (Miura)
+          rhodzm1 = rhodz
+          massflxc = 0.0 
+          uc = 0.0 
        END DO
-
-     IF ( iadvtr == iapp_tracvl ) THEN 
-	 bigt = dt*iapp_tracvl
-	 DO ind=1,ndomain
-           CALL swap_dimensions(ind)
-           CALL swap_geometry(ind)
-	      uc = f_uc(ind) 
-	      uc = uc/real(iapp_tracvl) 
-         END DO
+       CALL transfert_request(f_rhodzm1,req_i1)   !---->
+       CALL transfert_request(f_massflxc,req_e1) !---->
+       CALL transfert_request(f_massflxc,req_e1) !------>
+       CALL transfert_request(f_uc,req_e1) !---->
+       CALL transfert_request(f_uc,req_e1) 
+    END IF
+
+    iadvtr = iadvtr + 1 
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       massflx  = f_massflx(ind)
+       massflxc = f_massflxc(ind) 
+       uc = f_uc(ind) 
+       u  = f_u(ind) 
+       massflxc = massflxc + massflx 
+       uc = uc + u 
+    END DO
+
+    IF ( iadvtr == iapp_tracvl ) THEN 
+       bigt = dt*iapp_tracvl
+       DO ind=1,ndomain
+          CALL swap_dimensions(ind)
+          CALL swap_geometry(ind)
+          uc = f_uc(ind) 
+          uc = uc/real(iapp_tracvl) 
+          massflxc = f_massflx(ind)
+          massflxc = massflxc*dt
+          ! now massflx is time-integrated
+       END DO
 
        CALL vlsplt(f_q,f_rhodzm1,f_massflxc,2.0,f_uc,bigt) 
-            iadvtr = 0
-     END IF
+       iadvtr = 0
+    END IF
   END SUBROUTINE advect_tracer
-  
-!==============================================================================  
-  SUBROUTINE advtrac(massflx,wgg)
-  USE domain_mod
-  USE dimensions
-  USE grid_param
-  USE geometry
-  USE metric
-  USE disvert_mod
-  IMPLICIT NONE
-    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm)
-    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm)
-     
-    INTEGER :: i,j,ij,l
-    REAL(rstd) :: convm(iim*jjm,llm) 
-   
-     DO l = 1, llm
-      DO j=jj_begin,jj_end
-       DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  &
-                                ne(ij,rup)*massflx(ij+u_rup,l)       +  &
-                                ne(ij,lup)*massflx(ij+u_lup,l)       +  &
-                                ne(ij,left)*massflx(ij+u_left,l)     +  &
-                                ne(ij,ldown)*massflx(ij+u_ldown,l)   +  &
-                                ne(ij,rdown)*massflx(ij+u_rdown,l))
-        ENDDO
-      ENDDO
-     ENDDO
-
-    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
-
-    !!! Compute vertical velocity
-  DO l = 1,llm-1
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
-        wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )) 
-      ENDDO
-    ENDDO
-  ENDDO
-
-    DO j=jj_begin,jj_end
-     DO i=ii_begin,ii_end
-      ij=(j-1)*iim+i
-      wgg(ij,1)  = 0.
-     ENDDO
-   ENDDO
- END SUBROUTINE advtrac	
-
-	SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt) 
-  USE field_mod
-  USE domain_mod
-  USE dimensions
-  USE grid_param
-  USE geometry
-  USE metric
-  USE advect_mod
-  IMPLICIT NONE
-
-  TYPE(t_field),POINTER :: f_q(:)
-  TYPE(t_field),POINTER :: f_u(:)
-  TYPE(t_field),POINTER :: f_rhodz(:) 
-  TYPE(t_field),POINTER :: f_massflx(:)
-  TYPE(t_field),POINTER,SAVE :: f_wg(:)
-  TYPE(t_field),POINTER,SAVE :: f_zm(:)
-  TYPE(t_field),POINTER,SAVE :: f_zq(:)
-
-  REAL(rstd)::bigt,dt
-  REAL(rstd),POINTER :: q(:,:,:)
-  REAL(rstd),POINTER :: u(:,:)
-  REAL(rstd),POINTER :: rhodz(:,:) 
-  REAL(rstd),POINTER :: massflx(:,:) 
-  REAL(rstd),POINTER,SAVE :: wg(:,:)
-  REAL(rstd),POINTER,SAVE::zq(:,:,:) 
-  REAL(rstd),POINTER,SAVE::zm(:,:) 
-
-  REAL(rstd)::pente_max
-  LOGICAL,SAVE::first = .true. 
-  INTEGER :: i,ij,l,j,ind,k
-  REAL(rstd) :: zzpbar, zzw 
-  REAL::qvmax,qvmin 
-  
-     IF ( first ) THEN 
+
+  SUBROUTINE vlsplt(f_q,f_rhodz,f_massflx,pente_max,f_u,bigt)
+    USE field_mod
+    USE domain_mod
+    USE dimensions
+    USE grid_param
+    USE geometry
+    USE metric
+    USE advect_mod
+    IMPLICIT NONE
+
+    TYPE(t_field),POINTER :: f_q(:)
+    TYPE(t_field),POINTER :: f_u(:)
+    TYPE(t_field),POINTER :: f_rhodz(:) 
+    TYPE(t_field),POINTER :: f_massflx(:)
+
+    TYPE(t_field),POINTER,SAVE :: f_wg(:)
+    TYPE(t_field),POINTER,SAVE :: f_zm(:)
+    TYPE(t_field),POINTER,SAVE :: f_zq(:)
+
+    REAL(rstd)::bigt,dt
+    REAL(rstd),POINTER :: q(:,:,:)
+    REAL(rstd),POINTER :: u(:,:)
+    REAL(rstd),POINTER :: rhodz(:,:) 
+    REAL(rstd),POINTER :: massflx(:,:) 
+    REAL(rstd),POINTER :: wg(:,:)
+    REAL(rstd),POINTER :: zq(:,:,:) 
+    REAL(rstd),POINTER :: zm(:,:) 
+    REAL(rstd),POINTER :: tangent(:,:)
+    REAL(rstd),POINTER :: normal(:,:)
+    REAL(rstd),POINTER :: gradq3d(:,:,:)
+
+    REAL(rstd)::pente_max
+    LOGICAL,SAVE::first = .true. 
+    INTEGER :: i,ij,l,j,ind,k
+    REAL(rstd) :: zzpbar, zzw 
+    REAL::qvmax,qvmin 
+
+    IF ( first ) THEN 
        CALL allocate_field(f_wg,field_t,type_real,llm)
        CALL allocate_field(f_zm,field_t,type_real,llm)
        CALL allocate_field(f_zq,field_t,type_real,llm,nqtot) 
        first = .FALSE.
-      END IF 
-
-     DO ind=1,ndomain
-       CALL swap_dimensions(ind)
-       CALL swap_geometry(ind)
-         q=f_q(ind) 
-         rhodz=f_rhodz(ind) 
-         zq=f_zq(ind)
-	 zm=f_zm(ind) 
-         zm = rhodz ; zq = q  
-         wg = f_wg(ind)
-	 wg = 0.0
-         massflx=f_massflx(ind) 
-       CALL advtrac(massflx,wg) 
-     END DO
-
-!   CALL transfert_request(f_wg,req_i1) 
-
-      DO ind=1,ndomain
-        CALL swap_dimensions(ind)
-        CALL swap_geometry(ind)
-         zq=f_zq(ind)
-	 zm=f_zm(ind) 
-         wg=f_wg(ind)
-	 wg=wg*0.5
+    END IF
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       q=f_q(ind) 
+       rhodz=f_rhodz(ind) 
+       zq=f_zq(ind)
+       zm=f_zm(ind) 
+       zm = rhodz ; zq = q  
+       wg = f_wg(ind)
+       wg = 0.0
+       massflx=f_massflx(ind) 
+       CALL advtrac(massflx,wg) ! compute vertical mass fluxes
+    END DO
+
+    !   CALL transfert_request(f_wg,req_i1) 
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       zq=f_zq(ind)
+       zm=f_zm(ind) 
+       wg=f_wg(ind)
+       wg=wg*0.5
        DO k = 1, nqtot
-	 CALL vlz(zq(:,:,k),2.,zm,wg)
+          CALL vlz(zq(:,:,k),2.,zm,wg)
        END DO
-      END DO
-
-	DO ind=1,ndomain
-          CALL swap_dimensions(ind)
-          CALL swap_geometry(ind)
-          CALL swap_advect(ind)
-	  zq=f_zq(ind)
-	  zq = f_zq(ind)
-          zm = f_zm(ind)
-	  massflx =f_massflx(ind)
-	  u = f_u(ind) 
-         DO k = 1,nqtot
-           CALL advect1(zq(:,:,k))
-       	   CALL advect2(zq(:,:,k),zm,u,massflx,bigt) 
-         END DO
-	END DO
-
-	DO ind=1,ndomain
-           CALL swap_dimensions(ind)
-           CALL swap_geometry(ind)
-	    q = f_q(ind)
-	   zq = f_zq(ind) 
-	   zm = f_zm(ind) 
-	   wg = f_wg(ind)
-          DO k = 1,nqtot
-	    CALL vlz(zq(:,:,k),2.,zm,wg)
-	  END DO 
-            q = zq 
-	END DO
-
-  END SUBROUTINE vlsplt 
-
-	
-     SUBROUTINE vlz(q,pente_max,masse,wgg)
-!c
-!c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
-!c
-!c    ********************************************************************
-!c     Shema  d'advection " pseudo amont " .
-!c    ********************************************************************
-      USE icosa
-      IMPLICIT NONE
-!c
-!c   Arguments:
-!c   ----------
-      REAL masse(iim*jjm,llm),pente_max
-      REAL q(iim*jjm,llm)
-      REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1) 
-      REAL dq(iim*jjm,llm)
-      INTEGER i,ij,l,j,ii
-!c
-      REAL wq(iim*jjm,llm+1),newmasse
-
-      REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax
-      REAL sigw
-
-      REAL      SSUM
-
-
-	 w(:,1:llm) = wgg(:,:)
-	 w(:,llm+1) = 0.0 
-
-!c    On oriente tout dans le sens de la pression c'est a dire dans le
-!c    sens de W
-
-       DO l=2,llm
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
-            dzqw(ij,l)=q(ij,l-1)-q(ij,l)
-            adzqw(ij,l)=abs(dzqw(ij,l))
-         ENDDO
-        ENDDO
-       ENDDO
-
-       DO l=2,llm-1
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
-           IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
+    END DO
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       zq=f_zq(ind)
+       zq = f_zq(ind)
+       zm = f_zm(ind)
+       massflx =f_massflx(ind)
+       u = f_u(ind)
+       tangent = f_tangent(ind)
+       normal = f_normal(ind)
+       gradq3d = f_gradq3d(ind)
+
+       DO k = 1,nqtot
+          CALL compute_gradq3d(zq(:,:,k),gradq3d)
+          CALL compute_advect_horiz(zq(:,:,k),zm,u,massflx,bigt) 
+       END DO
+    END DO
+
+    DO ind=1,ndomain
+       CALL swap_dimensions(ind)
+       CALL swap_geometry(ind)
+       q = f_q(ind)
+       zq = f_zq(ind) 
+       zm = f_zm(ind) 
+       wg = f_wg(ind)
+       DO k = 1,nqtot
+          CALL vlz(zq(:,:,k),2.,zm,wg)
+       END DO
+       q = zq 
+    END DO
+
+  END SUBROUTINE vlsplt
+
+  !==============================================================================  
+  SUBROUTINE advtrac(massflx,wgg)
+    USE domain_mod
+    USE dimensions
+    USE grid_param
+    USE geometry
+    USE metric
+    USE disvert_mod
+    IMPLICIT NONE
+    REAL(rstd),INTENT(IN) :: massflx(iim*3*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: wgg(iim*jjm,llm)
+
+    INTEGER :: i,j,ij,l
+    REAL(rstd) :: convm(iim*jjm,llm) 
+
+    DO l = 1, llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             ! divergence of horizontal flux
+             convm(ij,l)= 1/(Ai(ij))*(ne(ij,right)*massflx(ij+u_right,l)   +  &
+                  ne(ij,rup)*massflx(ij+u_rup,l)       +  &
+                  ne(ij,lup)*massflx(ij+u_lup,l)       +  &
+                  ne(ij,left)*massflx(ij+u_left,l)     +  &
+                  ne(ij,ldown)*massflx(ij+u_ldown,l)   +  &
+                  ne(ij,rdown)*massflx(ij+u_rdown,l))
+          ENDDO
+       ENDDO
+    ENDDO
+
+    ! accumulate divergence from top of model
+    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
+
+!!! Compute vertical velocity
+    DO l = 1,llm-1
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             wgg( ij, l+1 ) = (convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )) 
+          ENDDO
+       ENDDO
+    ENDDO
+
+    DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          wgg(ij,1)  = 0.
+       ENDDO
+    ENDDO
+  END SUBROUTINE advtrac
+
+  SUBROUTINE vlz(q,pente_max,masse,wgg)
+    !c
+    !c     Auteurs:   P.Le Van, F.Hourdin, F.Forget 
+    !c
+    !c    ********************************************************************
+    !c     Shema  d'advection " pseudo amont " .
+    !c    ********************************************************************
+    USE icosa
+    IMPLICIT NONE
+    !c
+    !c   Arguments:
+    !c   ----------
+    REAL masse(iim*jjm,llm),pente_max
+    REAL q(iim*jjm,llm)
+    REAL wgg(iim*jjm,llm),w(iim*jjm,llm+1) 
+    REAL dq(iim*jjm,llm)
+    INTEGER i,ij,l,j,ii
+    !c
+    REAL wq(iim*jjm,llm+1),newmasse
+
+    REAL dzq(iim*jjm,llm),dzqw(iim*jjm,llm),adzqw(iim*jjm,llm),dzqmax
+    REAL sigw
+
+    REAL      SSUM
+
+
+    w(:,1:llm) = -wgg(:,:)  ! w>0 for downward transport
+    w(:,llm+1) = 0.0 
+
+    !c    On oriente tout dans le sens de la pression c'est a dire dans le
+    !c    sens de W
+
+    DO l=2,llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             dzqw(ij,l)=q(ij,l-1)-q(ij,l)
+             adzqw(ij,l)=abs(dzqw(ij,l))
+          ENDDO
+       ENDDO
+    ENDDO
+
+    DO l=2,llm-1
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
                 dzq(ij,l)=0.5*(dzqw(ij,l)+dzqw(ij,l+1))
-           ELSE
+             ELSE
                 dzq(ij,l)=0.
-           ENDIF
-            dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
-            dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
-         ENDDO
-        ENDDO
-	ENDDO
-
-       DO l=2,llm-1
-        DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
-            dzq(ij,1)=0.
-           dzq(ij,llm)=0.
-         ENDDO
-	 ENDDO
-      ENDDO
-!c ---------------------------------------------------------------
-!c   .... calcul des termes d'advection verticale  .......
-!c ---------------------------------------------------------------
-
-!c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
-
-        DO l = 1,llm-1
-   	   DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
+             ENDIF
+             dzqmax=pente_max*min(adzqw(ij,l),adzqw(ij,l+1))
+             dzq(ij,l)=sign(min(abs(dzq(ij,l)),dzqmax),dzq(ij,l))
+          ENDDO
+       ENDDO
+    ENDDO
+
+    DO l=2,llm-1
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             dzq(ij,1)=0.
+             dzq(ij,llm)=0.
+          ENDDO
+       ENDDO
+    ENDDO
+    !c ---------------------------------------------------------------
+    !c   .... calcul des termes d'advection verticale  .......
+    !c ---------------------------------------------------------------
+
+    !c calcul de  - d( q   * w )/ d(sigma)    qu'on ajoute a  dq pour calculer dq
+
+    DO l = 1,llm-1
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
              ij=(j-1)*iim+i
              IF(w(ij,l+1).gt.0.) THEN
-             sigw=w(ij,l+1)/masse(ij,l+1)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
+                ! upwind only if downward transport
+                sigw=w(ij,l+1)/masse(ij,l+1)
+                wq(ij,l+1)=w(ij,l+1)*(q(ij,l+1)+0.5*(1.-sigw)*dzq(ij,l+1))
              ELSE
-             sigw=w(ij,l+1)/masse(ij,l)
-             wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
-            ENDIF
-           ENDDO
-         ENDDO
-	 END DO
-
-	  DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-             wq(ij,llm+1)=0.
-             wq(ij,1)=0.
-           ENDDO
-	  END DO
-
-      DO l=1,llm
-       DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
-            newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l))
-	    dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>>
-            q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse    
-	                 
-            masse(ij,l)=newmasse
-	   ENDDO
-        ENDDO
-      END DO	
-      RETURN
-      END SUBROUTINE vlz
+                ! upwind only if upward transport
+                sigw=w(ij,l+1)/masse(ij,l)
+                wq(ij,l+1)=w(ij,l+1)*(q(ij,l)-0.5*(1.+sigw)*dzq(ij,l))
+             ENDIF
+          ENDDO
+       ENDDO
+    END DO
+
+    DO j=jj_begin,jj_end
+       DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i
+          wq(ij,llm+1)=0.
+          wq(ij,1)=0.
+       ENDDO
+    END DO
+
+    DO l=1,llm
+       DO j=jj_begin,jj_end
+          DO i=ii_begin,ii_end
+             ij=(j-1)*iim+i
+             ! masse -= dw/dz but w>0 <=> downward
+             newmasse=masse(ij,l)+(w(ij,l+1)-w(ij,l)) 
+!             dq(ij,l) = (wq(ij,l+1)-wq(ij,l)) !================>>>>
+             q(ij,l)=(q(ij,l)*masse(ij,l)+wq(ij,l+1)-wq(ij,l))/newmasse    
+!             masse(ij,l)=newmasse
+          ENDDO
+       ENDDO
+    END DO
+    RETURN
+  END SUBROUTINE vlz
 
 END MODULE advect_tracer_mod
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 21)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 22)
@@ -438,4 +438,5 @@
         ij=(j-1)*iim+i  
 ! signe ? attention d (rho theta dz)
+        ! dtheta_rhodz =  -div(flux.theta)
         dtheta_rhodz(ij,l)=-1./Ai(ij)*(ne(ij,right)*Ftheta(ij+u_right,l)    +  &
                                  ne(ij,rup)*Ftheta(ij+u_rup,l)        +  &  
@@ -458,5 +459,5 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
-!signe ?
+        ! 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)       +  &  
@@ -498,4 +499,5 @@
     DO i=ii_begin,ii_end
       ij=(j-1)*iim+i
+      ! dps/dt = -int(div flux)dz
       dps(ij)=-convm(ij,1) * g 
     ENDDO
@@ -510,4 +512,6 @@
       DO i=ii_begin,ii_end
         ij=(j-1)*iim+i
+        ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
+        ! => w>0 for upward transport
         w( ij, l+1 ) = convm( ij, l+1 ) - bp(l+1) * convm( ij, 1 )
       ENDDO
@@ -682,7 +686,8 @@
     DO j=jj_begin,jj_end
       DO i=ii_begin,ii_end
+         ! ww>0 <=> upward transport
         ij=(j-1)*iim+i
         ww            =   0.5 * w(ij,l+1) * (theta(ij,l) + theta(ij,l+1) )
-        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww 
+        dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -  ww
         dtheta_rhodz(ij,l+1) = dtheta_rhodz(ij,l+1)  +  ww
       ENDDO
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 21)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 22)
@@ -145,5 +145,5 @@
 
 
-    DO it=1,500
+    DO it=1,20
 
       CALL transfert_request(f_u,req_dissip)
@@ -212,5 +212,5 @@
 
 
-    DO it=1,500
+    DO it=1,20
  
       CALL transfert_request(f_u,req_dissip)
@@ -279,5 +279,5 @@
     ENDDO
 
-    DO it=1,500
+    DO it=1,20
 
       CALL transfert_request(f_theta,req_i1)
