Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 173)
+++ codes/icosagcm/trunk/src/advect.f90	(revision 174)
@@ -17,32 +17,30 @@
     REAL(rstd),INTENT(OUT) :: one_over_sqrt_leng(iim*jjm)
 
-    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)
-
-          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
+    INTEGER :: ij
+
+!$SIMD
+      DO ij=ij_begin,ij_end
+
+          CALL cross_product2(xyz_v(ij+z_rdown,:),xyz_v(ij+z_rup,:),normal(ij+u_right,:))
+          normal(ij+u_right,:)=normal(ij+u_right,:)/sqrt(sum(normal(ij+u_right,:)**2)+1e-50)*ne(ij,right)
+
+          CALL cross_product2(xyz_v(ij+z_up,:),xyz_v(ij+z_lup,:),normal(ij+u_lup,:))
+          normal(ij+u_lup,:)=normal(ij+u_lup,:)/sqrt(sum(normal(ij+u_lup,:)**2)+1e-50)*ne(ij,lup)
+
+          CALL cross_product2(xyz_v(ij+z_ldown,:),xyz_v(ij+z_down,:),normal(ij+u_ldown,:))	 
+          normal(ij+u_ldown,:)=normal(ij+u_ldown,:)/sqrt(sum(normal(ij+u_ldown,:)**2)+1e-50)*ne(ij,ldown)
+
+          tangent(ij+u_right,:)=xyz_v(ij+z_rup,:)-xyz_v(ij+z_rdown,:) 
+          tangent(ij+u_right,:)=tangent(ij+u_right,:)/sqrt(sum(tangent(ij+u_right,:)**2)+1e-50)
+
+          tangent(ij+u_lup,:)=xyz_v(ij+z_lup,:)-xyz_v(ij+z_up,:) 
+          tangent(ij+u_lup,:)=tangent(ij+u_lup,:)/sqrt(sum(tangent(ij+u_lup,:)**2)+1e-50)
+
+          tangent(ij+u_ldown,:)=xyz_v(ij+z_down,:)-xyz_v(ij+z_ldown,:)
+          tangent(ij+u_ldown,:)=tangent(ij+u_ldown,:)/sqrt(sum(tangent(ij+u_ldown,:)**2)+1e-50)
+
+          one_over_sqrt_leng(ij) = 1./sqrt(max(sum((xyz_v(ij+z_up,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_down,:) - xyz_i(ij,:))**2), &
+                                          sum((xyz_v(ij+z_rup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_rdown,:) - xyz_i(ij,:))**2),  &
+                                          sum((xyz_v(ij+z_lup,:) - xyz_i(ij,:))**2),sum((xyz_v(ij+z_ldown,:) - xyz_i(ij,:))**2)) )
     ENDDO
 
@@ -64,5 +62,5 @@
     REAL(rstd) :: ar(iim*jjm)
     REAL(rstd) :: gradtri(2*iim*jjm,llm,3) 
-    INTEGER :: i,j,n,k,ind,l
+    INTEGER :: ij,k,ind,l
 
     CALL trace_start("compute_gradq3d")
@@ -74,29 +72,23 @@
     ! arr = area of triangle joining centroids of hexagons
      DO l = ll_begin,ll_end 
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             n=(j-1)*iim+i   
-             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
+!$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
     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
+!$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
     ENDDO
       
     DO k=1,3
       DO l =ll_begin,ll_end
-        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
+!$SIMD
+      DO ij=ij_begin,ij_end
+              gradq3d(ij,l,k) = ( gradtri(ij+z_up,l,k) + gradtri(ij+z_down,l,k)   +   &
+                                 gradtri(ij+z_rup,l,k) +  gradtri(ij+z_ldown,l,k)   +   & 
+                                 gradtri(ij+z_lup,l,k)+    gradtri(ij+z_rdown,l,k) ) / ar(ij)  
         END DO
       END DO
@@ -105,22 +97,20 @@
     !============================================================================================= LIMITING 
     DO l =ll_begin,ll_end
-       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,:))
+!$SIMD
+      DO ij=ij_begin,ij_end
+             maggrd =  dot_product(gradq3d(ij,l,:),gradq3d(ij,l,:))
              maggrd = sqrt(maggrd) 
-             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))
-             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) )
+             maxq_c = qi(ij,l) + maggrd*one_over_sqrt_leng(ij)
+             minq_c = qi(ij,l) - maggrd*one_over_sqrt_leng(ij)
+             maxq = max(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
+                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
+             minq = min(qi(ij,l),qi(ij+t_right,l),qi(ij+t_lup,l),qi(ij+t_rup,l),qi(ij+t_left,l), &
+                  qi(ij+t_rdown,l),qi(ij+t_ldown,l))
+             alphamx = (maxq - qi(ij,l)) ; alphamx = alphamx/(maxq_c - qi(ij,l) )
              alphamx = max(alphamx,0.0)
-             alphami = (minq - qi(n,l)); alphami = alphami/(minq_c - qi(n,l))
+             alphami = (minq - qi(ij,l)); alphami = alphami/(minq_c - qi(ij,l))
              alphami = max(alphami,0.0) 
              alpha   = min(alphamx,alphami,1.0)
-             gradq3d(n,l,:) = alpha*gradq3d(n,l,:)
-          END DO
+             gradq3d(ij,l,:) = alpha*gradq3d(ij,l,:)
        END DO
     END DO
@@ -141,5 +131,5 @@
 
     REAL(rstd) :: v_e(3), up_e, qe, ed(3)    
-    INTEGER :: i,j,n,l
+    INTEGER :: ij,l
 
     CALL trace_start("compute_backward_traj")
@@ -149,53 +139,51 @@
     ! reconstruct tangential wind then 3D wind at edge then cc = edge midpoint - u*tau
     DO l = ll_begin,ll_end
-       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)         &
+!$SIMD
+      DO ij=ij_begin,ij_end
+             up_e =1/de(ij+u_right)*(                                                   &
+                  wee(ij+u_right,1,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
+                  wee(ij+u_right,2,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
+                  wee(ij+u_right,3,1)*le(ij+u_left)*ue(ij+u_left,l)+                      &
+                  wee(ij+u_right,4,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                    &
+                  wee(ij+u_right,5,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    & 
+                  wee(ij+u_right,1,2)*le(ij+t_right+u_ldown)*ue(ij+t_right+u_ldown,l)+    &
+                  wee(ij+u_right,2,2)*le(ij+t_right+u_rdown)*ue(ij+t_right+u_rdown,l)+    &
+                  wee(ij+u_right,3,2)*le(ij+t_right+u_right)*ue(ij+t_right+u_right,l)+    &
+                  wee(ij+u_right,4,2)*le(ij+t_right+u_rup)*ue(ij+t_right+u_rup,l)+        &
+                  wee(ij+u_right,5,2)*le(ij+t_right+u_lup)*ue(ij+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(ij+u_right,l)*normal(ij+u_right,:) + up_e*tangent(ij+u_right,:)
+             cc(ij+u_right,l,:) = xyz_e(ij+u_right,:) - v_e*tau
+
+             up_e=1/de(ij+u_lup)*(                                                      &
+                  wee(ij+u_lup,1,1)*le(ij+u_left)*ue(ij+u_left,l)+                        &
+                  wee(ij+u_lup,2,1)*le(ij+u_ldown)*ue(ij+u_ldown,l)+                      &
+                  wee(ij+u_lup,3,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                      &
+                  wee(ij+u_lup,4,1)*le(ij+u_right)*ue(ij+u_right,l)+                      &
+                  wee(ij+u_lup,5,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                          & 
+                  wee(ij+u_lup,1,2)*le(ij+t_lup+u_right)*ue(ij+t_lup+u_right,l)+          &
+                  wee(ij+u_lup,2,2)*le(ij+t_lup+u_rup)*ue(ij+t_lup+u_rup,l)+              &
+                  wee(ij+u_lup,3,2)*le(ij+t_lup+u_lup)*ue(ij+t_lup+u_lup,l)+              &
+                  wee(ij+u_lup,4,2)*le(ij+t_lup+u_left)*ue(ij+t_lup+u_left,l)+            &
+                  wee(ij+u_lup,5,2)*le(ij+t_lup+u_ldown)*ue(ij+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(ij+u_lup,l)*normal(ij+u_lup,:) + up_e*tangent(ij+u_lup,:)
+             cc(ij+u_lup,l,:) = xyz_e(ij+u_lup,:) - v_e*tau
+            
+
+             up_e=1/de(ij+u_ldown)*(                                                    &
+                  wee(ij+u_ldown,1,1)*le(ij+u_rdown)*ue(ij+u_rdown,l)+                    &
+                  wee(ij+u_ldown,2,1)*le(ij+u_right)*ue(ij+u_right,l)+                    &
+                  wee(ij+u_ldown,3,1)*le(ij+u_rup)*ue(ij+u_rup,l)+                        &
+                  wee(ij+u_ldown,4,1)*le(ij+u_lup)*ue(ij+u_lup,l)+                        &
+                  wee(ij+u_ldown,5,1)*le(ij+u_left)*ue(ij+u_left,l)+                      & 
+                  wee(ij+u_ldown,1,2)*le(ij+t_ldown+u_lup)*ue(ij+t_ldown+u_lup,l)+        &
+                  wee(ij+u_ldown,2,2)*le(ij+t_ldown+u_left)*ue(ij+t_ldown+u_left,l)+      &
+                  wee(ij+u_ldown,3,2)*le(ij+t_ldown+u_ldown)*ue(ij+t_ldown+u_ldown,l)+    &
+                  wee(ij+u_ldown,4,2)*le(ij+t_ldown+u_rdown)*ue(ij+t_ldown+u_rdown,l)+    &
+                  wee(ij+u_ldown,5,2)*le(ij+t_ldown+u_right)*ue(ij+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
+             v_e = ue(ij+u_ldown,l)*normal(ij+u_ldown,:) + up_e*tangent(ij+u_ldown,:)
+             cc(ij+u_ldown,l,:) = xyz_e(ij+u_ldown,:) - v_e*tau
        ENDDO
     END DO
@@ -220,5 +208,5 @@
     REAL(rstd) :: dq,dmass,qe,ed(3), newmass
     REAL(rstd) :: qflux(3*iim*jjm,llm)
-    INTEGER :: i,j,n,k,l
+    INTEGER :: ij,k,l
 
     CALL trace_start("compute_advect_horiz")
@@ -228,61 +216,61 @@
     ! ne*hfluxt>0 iff outgoing
     DO l = ll_begin,ll_end
-       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)*hfluxt(n+u_right,l)>0) then 
-                ed = cc(n+u_right,l,:) - xyz_i(n,:)
-                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
-             else
-                ed = cc(n+u_right,l,:) - xyz_i(n+t_right,:)
-                qe = qi(n+t_right,l)+sum2(gradq3d(n+t_right,l,:),ed) 
-             endif
-             qflux(n+u_right,l) = hfluxt(n+u_right,l)*qe
+
+!$SIMD
+      DO ij=ij_begin_ext,ij_end_ext
+
+             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) 
+             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) 
+             ENDIF
+             qflux(ij+u_right,l) = hfluxt(ij+u_right,l)*qe
              
-             if (ne(n,lup)*hfluxt(n+u_lup,l)>0) then 
-                ed = cc(n+u_lup,l,:) - xyz_i(n,:)
-                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed)
-             else
-                ed = cc(n+u_lup,l,:) - xyz_i(n+t_lup,:)
-                qe = qi(n+t_lup,l)+sum2(gradq3d(n+t_lup,l,:),ed) 
-             endif
-             qflux(n+u_lup,l) = hfluxt(n+u_lup,l)*qe 
-             
-             if (ne(n,ldown)*hfluxt(n+u_ldown,l)>0) then 
-                ed = cc(n+u_ldown,l,:) - xyz_i(n,:)
-                qe = qi(n,l)+sum2(gradq3d(n,l,:),ed) 
-             else
-                ed = cc(n+u_ldown,l,:) - xyz_i(n+t_ldown,:)
-                qe = qi(n+t_ldown,l)+sum2(gradq3d(n+t_ldown,l,:),ed) 
-             endif
-             qflux(n+u_ldown,l) = hfluxt(n+u_ldown,l)*qe
-          end do
-       end do
+             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)
+             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) 
+             ENDIF
+             qflux(ij+u_lup,l) = hfluxt(ij+u_lup,l)*qe 
+
+             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) 
+             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) 
+             ENDIF
+             qflux(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)*qe
+       END DO
     END DO
 
     ! update q and, if update_mass, update mass
     DO l = ll_begin,ll_end
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             n=(j-1)*iim+i  
+!$SIMD
+      DO ij=ij_begin,ij_end
              ! sign convention as Ringler et al. (2010) eq. 21
-             dmass =   hfluxt(n+u_right,l)*ne(n,right)   &
-                     + hfluxt(n+u_lup,l)  *ne(n,lup)     &
-                     + hfluxt(n+u_ldown,l)*ne(n,ldown)   &
-                     + hfluxt(n+u_rup,l)  *ne(n,rup)     &
-                     + hfluxt(n+u_left,l) *ne(n,left)    &
-                     + hfluxt(n+u_rdown,l)*ne(n,rdown)
-
-             dq    =   qflux(n+u_right,l) *ne(n,right)   &
-                     + qflux(n+u_lup,l)   *ne(n,lup)     &
-                     + qflux(n+u_ldown,l) *ne(n,ldown)   &
-                     + qflux(n+u_rup,l)   *ne(n,rup)     &
-                     + qflux(n+u_left,l)  *ne(n,left)    &
-                     + qflux(n+u_rdown,l) *ne(n,rdown)
-
-             newmass = mass(n,l) - dmass/Ai(n)
-             qi(n,l) = (qi(n,l)*mass(n,l) - dq/Ai(n) ) / newmass
-             IF(update_mass) mass(n,l) = newmass
-          END DO
+             dmass =   hfluxt(ij+u_right,l)*ne(ij,right)   &
+                     + hfluxt(ij+u_lup,l)  *ne(ij,lup)     &
+                     + hfluxt(ij+u_ldown,l)*ne(ij,ldown)   &
+                     + hfluxt(ij+u_rup,l)  *ne(ij,rup)     &
+                     + hfluxt(ij+u_left,l) *ne(ij,left)    &
+                     + hfluxt(ij+u_rdown,l)*ne(ij,rdown)
+
+             dq    =   qflux(ij+u_right,l) *ne(ij,right)   &
+                     + qflux(ij+u_lup,l)   *ne(ij,lup)     &
+                     + qflux(ij+u_ldown,l) *ne(ij,ldown)   &
+                     + qflux(ij+u_rup,l)   *ne(ij,rup)     &
+                     + qflux(ij+u_left,l)  *ne(ij,left)    &
+                     + qflux(ij+u_rdown,l) *ne(ij,rdown)
+
+             
+             newmass = mass(ij,l) - dmass/Ai(ij)
+             qi(ij,l) = (qi(ij,l)*mass(ij,l) - dq/Ai(ij) ) / newmass
+             IF(update_mass) mass(ij,l) = newmass
+
        END DO
     END DO
Index: codes/icosagcm/trunk/src/advect_tracer.f90
===================================================================
--- codes/icosagcm/trunk/src/advect_tracer.f90	(revision 173)
+++ codes/icosagcm/trunk/src/advect_tracer.f90	(revision 174)
@@ -10,5 +10,5 @@
   TYPE(t_field),POINTER :: f_one_over_sqrt_leng(:)
 
-  TYPE(t_message) :: req_u, req_wfluxt, req_q, req_rhodz, req_gradq3d
+  TYPE(t_message) :: req_u, req_cc, req_wfluxt, req_q, req_rhodz, req_gradq3d
 
   REAL(rstd), PARAMETER :: pente_max=2.0 ! for vlz
@@ -81,4 +81,5 @@
       first=.FALSE.
       CALL init_message(f_u,req_e1_vect,req_u)
+      CALL init_message(f_cc,req_e1_scal,req_cc)
       CALL init_message(f_wfluxt,req_i1,req_wfluxt)
       CALL init_message(f_q,req_i1,req_q)
@@ -124,4 +125,6 @@
     END DO
 
+    CALL send_message(f_cc,req_cc)
+
 
     ! horizontal transport - split in two to place transfer of gradq3d
@@ -138,4 +141,5 @@
 
        CALL send_message(f_gradq3d,req_gradq3d)
+       CALL wait_message(req_cc)
        CALL wait_message(req_gradq3d)
 
@@ -205,17 +209,18 @@
 
     REAL(rstd) :: dzqmax, newmass, sigw, qq, w
-    INTEGER :: i,ij,l,j
+    INTEGER :: i,ij,l,j,ijb,ije
 
     CALL trace_start("vlz")
+     
+     ijb=((jj_begin-halo)-1)*iim+ii_begin-halo
+     ije = ((jj_end+halo)-1)*iim+ii_end+halo
 
     ! finite difference of q
 
      DO l=ll_beginp1,ll_end
-       DO j=jj_begin-halo,jj_end+halo
-          DO i=ii_begin-halo,ii_end+halo
-             ij=(j-1)*iim+i
-             dzqw(ij,l)=q(ij,l)-q(ij,l-1)
-             adzqw(ij,l)=abs(dzqw(ij,l))
-          ENDDO
+!$SIMD
+       DO ij=ijb,ije
+         dzqw(ij,l)=q(ij,l)-q(ij,l-1)
+         adzqw(ij,l)=abs(dzqw(ij,l))
        ENDDO
     ENDDO
@@ -228,15 +233,13 @@
 
      DO l=ll_beginp1,ll_endm1
-       DO j=jj_begin-halo,jj_end+halo
-          DO i=ii_begin-halo,ii_end+halo
-             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) )
-                dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) )
-                dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b)
-             ELSE
-                dzq(ij,l)=0.
-             ENDIF
-          ENDDO
+!$SIMD
+       DO ij=ijb,ije 
+         IF(dzqw(ij,l)*dzqw(ij,l+1).gt.0.) THEN
+             dzq(ij,l) = 0.5*( dzqw(ij,l)+dzqw(ij,l+1) )
+             dzqmax    = pente_max * min( adzqw(ij,l),adzqw(ij,l+1) )
+             dzq(ij,l) = sign( min(abs(dzq(ij,l)),dzqmax) , dzq(ij,l) )  ! NB : sign(a,b)=a*sign(b)
+          ELSE
+             dzq(ij,l)=0.
+          ENDIF
        ENDDO
     ENDDO
@@ -245,17 +248,12 @@
     ! 0 slope in top and bottom layers
     IF (omp_first) THEN
-      DO j=jj_begin-halo,jj_end+halo
-         DO i=ii_begin-halo,ii_end+halo
-           ij=(j-1)*iim+i
+      DO ij=ijb,ije
            dzq(ij,1)=0.
-         ENDDO
       ENDDO
     ENDIF
       
     IF (omp_last) THEN
-      DO j=jj_begin-halo,jj_end+halo
-         DO i=ii_begin-halo,ii_end+halo
+      DO ij=ijb,ije
           dzq(ij,llm)=0.
-         ENDDO
       ENDDO
     ENDIF
@@ -267,7 +265,6 @@
     ! then amount of q leaving level l/l+1 = wq = w * qq
      DO l=ll_beginp1,ll_end
-       DO j=jj_begin-halo,jj_end+halo
-          DO i=ii_begin-halo,ii_end+halo
-             ij=(j-1)*iim+i
+!$SIMD
+       DO ij=ijb,ije
              w = fac*wfluxt(ij,l)
              IF(w>0.) THEN  ! upward transport, upwind side is at level l 
@@ -279,23 +276,16 @@
              ENDIF
              wq(ij,l) = w*qq
-          ENDDO
        ENDDO
     END DO
     ! wq = 0 at top and bottom
     IF (omp_first) THEN
-      DO j=jj_begin-halo,jj_end+halo
-         DO i=ii_begin-halo,ii_end+halo
-            ij=(j-1)*iim+i
+       DO ij=ijb,ije
             wq(ij,1)=0.
-         ENDDO
       END DO
     ENDIF
     
     IF (omp_last) THEN
-      DO j=jj_begin-halo,jj_end+halo
-         DO i=ii_begin-halo,ii_end+halo
-            ij=(j-1)*iim+i
+      DO ij=ijb,ije
             wq(ij,llm+1)=0.
-         ENDDO
       END DO
     ENDIF
@@ -307,11 +297,9 @@
     ! update q, mass is updated only after all q's have been updated
     DO l=ll_begin,ll_end
-       DO j=jj_begin-halo,jj_end+halo
-          DO i=ii_begin-halo,ii_end+halo
-             ij=(j-1)*iim+i
+!$SIMD
+       DO ij=ijb,ije
              newmass = mass(ij,l) + fac*(wfluxt(ij,l)-wfluxt(ij,l+1))
              q(ij,l) = ( q(ij,l)*mass(ij,l) + wq(ij,l)-wq(ij,l+1) ) / newmass
              IF(update_mass) mass(ij,l)=newmass
-          ENDDO
        ENDDO
     END DO
Index: codes/icosagcm/trunk/src/caldyn_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 173)
+++ codes/icosagcm/trunk/src/caldyn_gcm.f90	(revision 174)
@@ -97,19 +97,16 @@
           wwuu=f_wwuu(ind)
           
-          DO j=jj_begin-1,jj_end+1
-             DO i=ii_begin-1,ii_end+1
-                ij=(j-1)*iim+i
-                ! lower BCs : geopot=phis, wflux=0, wwuu=0
-                geopot(ij,1) = phis(ij)
-                wflux(ij,1) = 0.
-                wwuu(ij+u_right,1)=0   
-                wwuu(ij+u_lup,1)=0   
-                wwuu(ij+u_ldown,1)=0
-                ! top BCs : wflux=0, wwuu=0
-                wflux(ij,llm+1)  = 0.
-                wwuu(ij+u_right,llm+1)=0   
-                wwuu(ij+u_lup,llm+1)=0   
-                wwuu(ij+u_ldown,llm+1)=0
-             ENDDO
+          DO ij=ij_begin_ext,ij_end_ext
+              ! lower BCs : geopot=phis, wflux=0, wwuu=0
+              geopot(ij,1) = phis(ij)
+              wflux(ij,1) = 0.
+              wwuu(ij+u_right,1)=0   
+              wwuu(ij+u_lup,1)=0   
+              wwuu(ij+u_ldown,1)=0
+              ! top BCs : wflux=0, wwuu=0
+              wflux(ij,llm+1)  = 0.
+              wwuu(ij+u_right,llm+1)=0   
+              wwuu(ij+u_lup,llm+1)=0   
+              wwuu(ij+u_ldown,llm+1)=0
           ENDDO
        END DO
@@ -303,5 +300,4 @@
 !  REAL(rstd) :: qv(2*iim*jjm,llm)     ! potential velocity  
   
-
   CALL trace_start("compute_pvort")  
 
@@ -316,11 +312,9 @@
      DO l = ll_begin,ll_end
         CALL test_message(req_u) 
-        DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-              ij=(j-1)*iim+i
-              m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
-              rhodz(ij,l) = m
-              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-           ENDDO
+!DIR$ SIMD
+        DO ij=ij_begin_ext,ij_end_ext
+           m = ( mass_dak(l)+ps(ij)*mass_dbk(l) )/g
+           rhodz(ij,l) = m
+           theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
         ENDDO
      ENDDO
@@ -328,10 +322,8 @@
      DO l = ll_begin,ll_end
         CALL test_message(req_u) 
-        DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-              ij=(j-1)*iim+i
-              theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
-           ENDDO
-        ENDDO
+!DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+         theta(ij,l) = theta_rhodz(ij,l)/rhodz(ij,l)
+       ENDDO
      ENDDO
   END IF
@@ -341,9 +333,6 @@
 !!! Compute shallow-water potential vorticity
   DO l = ll_begin,ll_end
-
-    DO j=jj_begin-1,jj_end+1
-       DO i=ii_begin-1,ii_end+1
-          ij=(j-1)*iim+i
-          
+!DIR$ SIMD
+     DO ij=ij_begin_ext,ij_end_ext         
           etav= 1./Av(ij+z_up)*(  ne_rup        * u(ij+u_rup,l)        * de(ij+u_rup)         &
                                 + ne_left * u(ij+t_rup+u_left,l) * de(ij+t_rup+u_left)  &
@@ -366,14 +355,11 @@
           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
 
-        ENDDO
       ENDDO
 
-      DO j=jj_begin,jj_end
-         DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
-            qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
-            qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
-            qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
-         END DO
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end
+         qu(ij+u_right,l) = 0.5*(qv(ij+z_rdown,l)+qv(ij+z_rup,l))  
+         qu(ij+u_lup,l) = 0.5*(qv(ij+z_up,l)+qv(ij+z_lup,l))  
+         qu(ij+u_ldown,l) = 0.5*(qv(ij+z_ldown,l)+qv(ij+z_down,l))  
       END DO
       
@@ -407,14 +393,12 @@
        DO l = 1,llm
           !$OMP DO SCHEDULE(STATIC) 
-          DO j=jj_begin-1,jj_end+1
-             DO i=ii_begin-1,ii_end+1
-                ij=(j-1)*iim+i
-                p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
-                !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
-                exner_ik = cpp * (p_ik/preff) ** kappa
-                pk(ij,l) = exner_ik
-                ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
-                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
-             ENDDO
+!DIR$ SIMD
+         DO ij=ij_begin_ext,ij_end_ext         
+            p_ik = ptop + mass_ak(l) + mass_bk(l)*ps(ij) ! FIXME : leave ps for the moment ; change ps to Ms later
+            !         p_ik = ptop + g*(mass_ak(l)+ mass_bk(l)*ps(i,j))
+            exner_ik = cpp * (p_ik/preff) ** kappa
+            pk(ij,l) = exner_ik
+            ! specific volume v = kappa*theta*pi/p = dphi/g/rhodz
+            geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik
           ENDDO
        ENDDO
@@ -427,26 +411,19 @@
 
        ! uppermost layer
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
-             pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
-          END DO
+!DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext         
+         pk(ij,llm) = ptop + (.5*g)*rhodz(ij,llm)
        END DO
        ! other layers
        DO l = llm-1, 1, -1
           !$OMP DO SCHEDULE(STATIC) 
-          DO j=jj_begin-1,jj_end+1
-             DO i=ii_begin-1,ii_end+1
-                ij=(j-1)*iim+i
-                pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
-             END DO
+!DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext         
+             pk(ij,l) = pk(ij,l+1) + (.5*g)*(rhodz(ij,l)+rhodz(ij,l+1))
           END DO
        END DO
        ! surface pressure (for diagnostics)
-       DO j=jj_begin-1,jj_end+1
-          DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
-             ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
-          END DO
+       DO ij=ij_begin_ext,ij_end_ext         
+         ps(ij) = pk(ij,1) + (.5*g)*rhodz(ij,1)
        END DO
 
@@ -454,9 +431,7 @@
           DO l = 1,llm
              !$OMP DO SCHEDULE(STATIC) 
-             DO j=jj_begin-1,jj_end+1
-                DO i=ii_begin-1,ii_end+1
-                   ij=(j-1)*iim+i
-                   geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
-                ENDDO
+!DIR$ SIMD
+             DO ij=ij_begin_ext,ij_end_ext         
+               geopot(ij,l+1) = geopot(ij,l) + g*rhodz(ij,l)
              ENDDO
           ENDDO
@@ -464,12 +439,10 @@
           DO l = 1,llm
              !$OMP DO SCHEDULE(STATIC) 
-             DO j=jj_begin-1,jj_end+1
-                DO i=ii_begin-1,ii_end+1
-                   ij=(j-1)*iim+i
-                   p_ik = pk(ij,l)
-                   exner_ik = cpp * (p_ik/preff) ** kappa
-                   geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
-                   pk(ij,l) = exner_ik
-                ENDDO
+!DIR$ SIMD
+              DO ij=ij_begin_ext,ij_end_ext         
+                p_ik = pk(ij,l)
+                exner_ik = cpp * (p_ik/preff) ** kappa
+                geopot(ij,l+1) = geopot(ij,l) + (g*kappa)*rhodz(ij,l)*theta(ij,l)*exner_ik/p_ik 
+                pk(ij,l) = exner_ik
              ENDDO
           ENDDO
@@ -514,7 +487,6 @@
 !!!  Compute mass and theta fluxes
     IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-    DO j=jj_begin-1,jj_end+1
-      DO i=ii_begin-1,ii_end+1
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext         
         hflux(ij+u_right,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)*le(ij+u_right)
         hflux(ij+u_lup,l)=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)*le(ij+u_lup)
@@ -524,11 +496,9 @@
         Ftheta(ij+u_lup,l)=0.5*(theta(ij,l)+theta(ij+t_lup,l))*hflux(ij+u_lup,l)
         Ftheta(ij+u_ldown,l)=0.5*(theta(ij,l)+theta(ij+t_ldown,l))*hflux(ij+u_ldown,l)
-      ENDDO
     ENDDO
 
 !!! compute horizontal divergence of fluxes
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i  
+!DIR$ SIMD
+    DO ij=ij_begin,ij_end         
         ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
         convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
@@ -547,5 +517,4 @@
                                  ne_ldown*Ftheta(ij+u_ldown,l)    +  &  
                                  ne_rdown*Ftheta(ij+u_rdown,l))
-      ENDDO
     ENDDO
 
@@ -560,8 +529,7 @@
 
         DO l=ll_begin,ll_end
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
-
+!DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+    
                 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
                      wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
@@ -601,5 +569,4 @@
                 du(ij+u_ldown,l) = .5*uu/de(ij+u_ldown)
                 
-             ENDDO
           ENDDO
        ENDDO
@@ -608,7 +575,6 @@
   
         DO l=ll_begin,ll_end
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
+!DIR$ SIMD
+          DO ij=ij_begin,ij_end         
 
                 uu = wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
@@ -649,5 +615,4 @@
                 du(ij+u_ldown,l) = qu(ij+u_ldown,l)*uu/de(ij+u_ldown)
      
-             ENDDO
           ENDDO
        ENDDO
@@ -661,7 +626,6 @@
 
        DO l=ll_begin,ll_end
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
+!DIR$ SIMD
+          DO ij=ij_begin,ij_end         
                 
                 berni(ij,l) = pk(ij,l) + &
@@ -673,5 +637,4 @@
                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
                 pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
-             ENDDO
           ENDDO
        ENDDO
@@ -680,7 +643,6 @@
 
        DO l=ll_begin,ll_end
-          DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-                ij=(j-1)*iim+i
+!DIR$ SIMD
+            DO ij=ij_begin,ij_end         
                 
                 berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
@@ -691,6 +653,4 @@
                                      le(ij+u_ldown)*de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
                                      le(ij+u_rdown)*de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-       
-             ENDDO
           ENDDO
        ENDDO
@@ -700,7 +660,6 @@
 !!! Add gradients of Bernoulli and Exner functions to du
     DO l=ll_begin,ll_end
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
+!DIR$ SIMD
+       DO ij=ij_begin,ij_end         
         
              du(ij+u_right,l) = du(ij+u_right,l) + 1/de(ij+u_right) * (             &
@@ -720,5 +679,4 @@
                                    + ne_ldown*berni(ij,l)+ne_rup*berni(ij+t_ldown,l) )
 
-          ENDDO
        ENDDO
     ENDDO
@@ -750,4 +708,5 @@
     REAL(rstd) :: p_ik, exner_ik
 
+
   CALL trace_start("compute_caldyn_vert")
 
@@ -757,9 +716,7 @@
     IF (caldyn_conserv==energy) CALL test_message(req_qu) 
 !$OMP DO SCHEDULE(STATIC) 
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+    DO ij=ij_begin,ij_end         
         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
-      ENDDO
     ENDDO
   ENDDO
@@ -770,10 +727,8 @@
 ! compute dps
   IF (omp_first) THEN
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+     DO ij=ij_begin,ij_end         
         ! dps/dt = -int(div flux)dz
         dps(ij) = convm(ij,1) * g 
-      ENDDO
     ENDDO
   ENDIF
@@ -782,29 +737,23 @@
   DO l=ll_beginp1,ll_end
     IF (caldyn_conserv==energy) CALL test_message(req_qu) 
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end         
         ! w = int(z,ztop,div(flux)dz) + B(eta)dps/dt
         ! => w>0 for upward transport
         wflux( ij, l ) = bp(l) * convm( ij, 1 ) - convm( ij, l ) 
-      ENDDO
     ENDDO
   ENDDO
 
   DO l=ll_begin,ll_endm1
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+     DO ij=ij_begin,ij_end         
         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  -   0.5 * (  wflux(ij,l+1) * (theta(ij,l) + theta(ij,l+1)))
-      ENDDO
     ENDDO
   ENDDO
   
   DO l=ll_beginp1,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end         
         dtheta_rhodz(ij, l ) = dtheta_rhodz(ij, l )  +   0.5 * ( wflux(ij,l  ) * (theta(ij,l-1) + theta(ij,l) ) )
-      ENDDO
     ENDDO
   ENDDO
@@ -812,11 +761,9 @@
 ! Compute vertical transport
   DO l=ll_beginp1,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+      DO ij=ij_begin,ij_end         
         wwuu(ij+u_right,l) = 0.5*( wflux(ij,l) + wflux(ij+t_right,l)) * (u(ij+u_right,l) - u(ij+u_right,l-1))
         wwuu(ij+u_lup,l) = 0.5* ( wflux(ij,l) + wflux(ij+t_lup,l)) * (u(ij+u_lup,l) - u(ij+u_lup,l-1))
         wwuu(ij+u_ldown,l) = 0.5*( wflux(ij,l) + wflux(ij+t_ldown,l)) * (u(ij+u_ldown,l) - u(ij+u_ldown,l-1))
-       ENDDO
      ENDDO
    ENDDO
@@ -827,11 +774,9 @@
 ! Add vertical transport to du
   DO l=ll_begin,ll_end
-    DO j=jj_begin,jj_end
-      DO i=ii_begin,ii_end
-        ij=(j-1)*iim+i
+!DIR$ SIMD
+     DO ij=ij_begin,ij_end         
         du(ij+u_right, l )   = du(ij+u_right,l)  - (wwuu(ij+u_right,l+1)+ wwuu(ij+u_right,l)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
         du(ij+u_lup, l )     = du(ij+u_lup,l)    - (wwuu(ij+u_lup,l+1)  + wwuu(ij+u_lup,l))   / (rhodz(ij,l)+rhodz(ij+t_lup,l))
         du(ij+u_ldown, l )   = du(ij+u_ldown,l)  - (wwuu(ij+u_ldown,l+1)+ wwuu(ij+u_ldown,l)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
-      ENDDO
     ENDDO
   ENDDO      
Index: codes/icosagcm/trunk/src/dimensions.f90
===================================================================
--- codes/icosagcm/trunk/src/dimensions.f90	(revision 173)
+++ codes/icosagcm/trunk/src/dimensions.f90	(revision 174)
@@ -9,5 +9,9 @@
     INTEGER :: ii_nb
     INTEGER :: jj_nb
-    
+    INTEGER :: ij_begin
+    INTEGER :: ij_end
+    INTEGER :: ij_begin_ext
+    INTEGER :: ij_end_ext
+        
     INTEGER :: t_right
     INTEGER :: t_rup
@@ -56,5 +60,10 @@
      ii_nb=d%ii_nb
      jj_nb=d%jj_nb
-     
+    
+     ij_begin=(jj_begin-1)*iim+ii_begin
+     ij_end=(jj_end-1)*iim+ii_end 
+     ij_begin_ext=((jj_begin-1)-1)*iim+ii_begin-1
+     ij_end_ext = ((jj_end+1)-1)*iim+ii_end+1
+        
      t_right=d%t_right
      t_rup=d%t_rup
Index: codes/icosagcm/trunk/src/dissip_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 173)
+++ codes/icosagcm/trunk/src/dissip_gcm.f90	(revision 174)
@@ -77,7 +77,9 @@
   CHARACTER(len=255)    :: rayleigh_friction_key
   REAL(rstd)            :: mintau
+  INTEGER               :: seed_size
+  INTEGER,ALLOCATABLE   :: seed(:)
   
             
-  INTEGER :: i,j,n,ind,it,iter
+  INTEGER :: i,j,ij,ind,it,iter
 
    rayleigh_friction_key='none'
@@ -142,6 +144,4 @@
     cdivgrad=-1
     cgradrot=-1
-  
-    CALL RANDOM_SEED()
     
     DO ind=1,ndomain
@@ -152,14 +152,13 @@
       DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
+          ij=(j-1)*iim+i   
           CALL RANDOM_NUMBER(r)
-          u(n+u_right)=r-0.5
+          u(ij+u_right)=r-0.5
           CALL RANDOM_NUMBER(r)
-          u(n+u_lup)=r-0.5
+          u(ij+u_lup)=r-0.5
           CALL RANDOM_NUMBER(r)
-          u(n+u_ldown)=r-0.5
-       ENDDO
-      ENDDO
-        
+          u(ij+u_ldown)=r-0.5
+        ENDDO
+      ENDDO        
     ENDDO
   
@@ -191,8 +190,8 @@
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
-            n=(j-1)*iim+i
-            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
-            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
-            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
+            ij=(j-1)*iim+i   
+            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
+            if (le(ij+u_lup)>1e-100)   dumax=MAX(dumax,ABS(du(ij+u_lup)))
+            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
           ENDDO
         ENDDO
@@ -228,16 +227,15 @@
       u=f_u(ind)
 
-     DO j=jj_begin,jj_end
+      DO j=jj_begin,jj_end
         DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
+          ij=(j-1)*iim+i   
           CALL RANDOM_NUMBER(r)
-          u(n+u_right)=r-0.5
+          u(ij+u_right)=r-0.5
           CALL RANDOM_NUMBER(r)
-          u(n+u_lup)=r-0.5
+          u(ij+u_lup)=r-0.5
           CALL RANDOM_NUMBER(r)
-          u(n+u_ldown)=r-0.5
-       ENDDO
-      ENDDO
-        
+          u(ij+u_ldown)=r-0.5
+        ENDDO
+      ENDDO        
     ENDDO
 
@@ -268,8 +266,8 @@
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
-            n=(j-1)*iim+i
-            if (le(n+u_right)>1e-100) dumax=MAX(dumax,ABS(du(n+u_right)))
-            if (le(n+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(n+u_lup)))
-            if (le(n+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(n+u_ldown)))
+            ij=(j-1)*iim+i   
+            if (le(ij+u_right)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_right)))
+            if (le(ij+u_lup)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_lup)))
+            if (le(ij+u_ldown)>1e-100) dumax=MAX(dumax,ABS(du(ij+u_ldown)))
           ENDDO
         ENDDO
@@ -304,12 +302,11 @@
       theta=f_theta(ind)
  
-       DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            n=(j-1)*iim+i
-            CALL RANDOM_NUMBER(r)
-            theta(n)=r-0.5
-         ENDDO
-        ENDDO
-        
+      DO j=jj_begin,jj_end
+        DO i=ii_begin,ii_end
+          ij=(j-1)*iim+i   
+          CALL RANDOM_NUMBER(r)
+          theta(ij)=r-0.5
+        ENDDO
+      ENDDO  
     ENDDO
 
@@ -339,13 +336,15 @@
         DO j=jj_begin,jj_end
           DO i=ii_begin,ii_end
-            n=(j-1)*iim+i
-            dthetamax=MAX(dthetamax,ABS(dtheta(n)))
+            ij=(j-1)*iim+i   
+            dthetamax=MAX(dthetamax,ABS(dtheta(ij)))
           ENDDO
         ENDDO
       ENDDO
+      
       IF (using_mpi) THEN
         CALL MPI_ALLREDUCE(dthetamax,dthetamax1,1,MPI_REAL8,MPI_MAX,comm_icosa,ierr)
 	dthetamax=dthetamax1
       ENDIF
+      
       IF (is_mpi_root) PRINT *,"divgrad : it :",it ,": dthetamax",dthetamax
 
@@ -417,5 +416,5 @@
 
     INTEGER :: ind, shear
-    INTEGER :: l,i,j,n
+    INTEGER :: l,ij
 
 !$OMP BARRIER
@@ -444,14 +443,12 @@
 
       DO l=ll_begin,ll_end
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            n=(j-1)*iim+i
-
-            due(n+u_right,l) = -0.5*( due_diss1(n+u_right,l)/tau_graddiv(l) + due_diss2(n+u_right,l)/tau_gradrot(l))*itau_dissip 
-            due(n+u_lup,l)   = -0.5*( due_diss1(n+u_lup,l)  /tau_graddiv(l) + due_diss2(n+u_lup,l)  /tau_gradrot(l))*itau_dissip
-            due(n+u_ldown,l) = -0.5*( due_diss1(n+u_ldown,l)/tau_graddiv(l) + due_diss2(n+u_ldown,l)/tau_gradrot(l))*itau_dissip
-
-            dtheta_rhodz(n,l) = -0.5*dtheta_rhodz_diss(n,l)/tau_divgrad(l)*itau_dissip
-          ENDDO
+!$SIMD
+        DO ij=ij_begin,ij_end
+
+            due(ij+u_right,l) = -0.5*( due_diss1(ij+u_right,l)/tau_graddiv(l) + due_diss2(ij+u_right,l)/tau_gradrot(l))*itau_dissip 
+            due(ij+u_lup,l)   = -0.5*( due_diss1(ij+u_lup,l)  /tau_graddiv(l) + due_diss2(ij+u_lup,l)  /tau_gradrot(l))*itau_dissip
+            due(ij+u_ldown,l) = -0.5*( due_diss1(ij+u_ldown,l)/tau_graddiv(l) + due_diss2(ij+u_ldown,l)/tau_gradrot(l))*itau_dissip
+
+            dtheta_rhodz(ij,l) = -0.5*dtheta_rhodz_diss(ij,l)/tau_divgrad(l)*itau_dissip
         ENDDO
       ENDDO
@@ -490,8 +487,8 @@
         LOGICAL :: hybrid_eta
         INTEGER :: shift_u, shift_t, zcoords, nn
-        z = (phi(n,l)+phi(n+shift_t,l))/(2.*g)
+        z = (phi(ij,l)+phi(ij+shift_t,l))/(2.*g)
         IF(z>zh) THEN  ! relax only in the sponge layer z>zh
 
-           nn = n+shift_u
+           nn = ij+shift_u
            CALL xyz2lonlat(xyz_e(nn,:),lon,lat)
            zcoords = 1 ; hybrid_eta = .FALSE. ! use z instead of p or hyam/hybm
@@ -522,5 +519,5 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it,i,j,l,ij
+    INTEGER :: it,l,ij
        
     CALL trace_start("gradiv")
@@ -532,11 +529,9 @@
       due=f_due(ind) 
       DO  l = ll_begin, ll_end
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
+!$SIMD
+        DO ij=ij_begin,ij_end
              due(ij+u_right,l)=ue(ij+u_right,l)
              due(ij+u_lup,l)=ue(ij+u_lup,l)
              due(ij+u_ldown,l)=ue(ij+u_ldown,l)
-          ENDDO
         ENDDO
       ENDDO
@@ -571,5 +566,5 @@
     REAL(rstd),POINTER    :: due(:,:)
     INTEGER :: ind
-    INTEGER :: it,i,j,l,ij
+    INTEGER :: it,l,ij
        
     CALL trace_start("gradrot")
@@ -581,11 +576,9 @@
       due=f_due(ind) 
       DO  l = ll_begin, ll_end
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
+!$SIMD
+        DO ij=ij_begin,ij_end
              due(ij+u_right,l)=ue(ij+u_right,l)
              due(ij+u_lup,l)=ue(ij+u_lup,l)
              due(ij+u_ldown,l)=ue(ij+u_ldown,l)
-          ENDDO
         ENDDO
       ENDDO
@@ -663,5 +656,5 @@
 
     INTEGER :: ind
-    INTEGER :: it,i,j,l,ij
+    INTEGER :: it,l,ij
 
     CALL trace_start("divgrad")
@@ -674,9 +667,7 @@
       dtheta_rhodz=f_dtheta_rhodz(ind) 
       DO  l = ll_begin, ll_end
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
+!$SIMD
+        DO ij=ij_begin,ij_end
             dtheta_rhodz(ij,l) = theta_rhodz(ij,l) / mass(ij,l)
-          ENDDO
         ENDDO
       ENDDO
@@ -702,9 +693,7 @@
     
       DO  l = ll_begin, ll_end
-        DO j=jj_begin,jj_end
-          DO i=ii_begin,ii_end
-            ij=(j-1)*iim+i
+!$SIMD
+        DO ij=ij_begin,ij_end
             dtheta_rhodz(ij,l) = dtheta_rhodz(ij,l) * mass(ij,l)
-          ENDDO
         ENDDO
       ENDDO
@@ -725,44 +714,37 @@
     REAL(rstd) :: divu_i(iim*jjm,llb:lle)
     
-    INTEGER :: i,j,n,l
+    INTEGER :: ij,l
 
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          divu_i(n,l)=1./Ai(n)*(ne(n,right)*ue(n+u_right,l)*le(n+u_right)  +  &
-                             ne(n,rup)*ue(n+u_rup,l)*le(n+u_rup)        +  &  
-                             ne(n,lup)*ue(n+u_lup,l)*le(n+u_lup)        +  &  
-                             ne(n,left)*ue(n+u_left,l)*le(n+u_left)     +  &  
-                             ne(n,ldown)*ue(n+u_ldown,l)*le(n+u_ldown)  +  &  
-                             ne(n,rdown)*ue(n+u_rdown,l)*le(n+u_rdown))
-        ENDDO
+!$SIMD
+      DO ij=ij_begin,ij_end
+          divu_i(ij,l)=1./Ai(ij)*(ne(ij,right)*ue(ij+u_right,l)*le(ij+u_right)  +  &
+                             ne(ij,rup)*ue(ij+u_rup,l)*le(ij+u_rup)        +  &  
+                             ne(ij,lup)*ue(ij+u_lup,l)*le(ij+u_lup)        +  &  
+                             ne(ij,left)*ue(ij+u_left,l)*le(ij+u_left)     +  &  
+                             ne(ij,ldown)*ue(ij+u_ldown,l)*le(ij+u_ldown)  +  &  
+                             ne(ij,rdown)*ue(ij+u_rdown,l)*le(ij+u_rdown))
       ENDDO
     ENDDO
     
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
+!$SIMD
+      DO ij=ij_begin,ij_end
  
-          n=(j-1)*iim+i
- 
-          gradivu_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*divu_i(n,l)+ ne(n+t_right,left)*divu_i(n+t_right,l) )       
-
-          gradivu_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*divu_i(n,l)+ ne(n+t_lup,rdown)*divu_i(n+t_lup,l))       
-    
-          gradivu_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*divu_i(n,l)+ne(n+t_ldown,rup)*divu_i(n+t_ldown,l) )
-
-        ENDDO
+          gradivu_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*divu_i(ij,l)+ ne(ij+t_right,left)*divu_i(ij+t_right,l) )       
+
+          gradivu_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*divu_i(ij,l)+ ne(ij+t_lup,rdown)*divu_i(ij+t_lup,l))       
+    
+          gradivu_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*divu_i(ij,l)+ne(ij+t_ldown,rup)*divu_i(ij+t_ldown,l) )
+
       ENDDO
     ENDDO
 
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          gradivu_e(n+u_right,l)=-gradivu_e(n+u_right,l)*cgraddiv
-          gradivu_e(n+u_lup,l)=-gradivu_e(n+u_lup,l)*cgraddiv
-          gradivu_e(n+u_ldown,l)=-gradivu_e(n+u_ldown,l)*cgraddiv
-        ENDDO
+!$SIMD
+      DO ij=ij_begin,ij_end
+          gradivu_e(ij+u_right,l)=-gradivu_e(ij+u_right,l)*cgraddiv
+          gradivu_e(ij+u_lup,l)=-gradivu_e(ij+u_lup,l)*cgraddiv
+          gradivu_e(ij+u_ldown,l)=-gradivu_e(ij+u_ldown,l)*cgraddiv
       ENDDO
     ENDDO
@@ -780,20 +762,17 @@
     REAL(rstd)  :: grad_e(3*iim*jjm,llb:lle)
 
-    INTEGER :: i,j,n,l
+    INTEGER :: ij,l
 
        
     DO l=llb,lle
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-   
-          n=(j-1)*iim+i
+!$SIMD
+      DO ij=ij_begin_ext,ij_end_ext
  
-          grad_e(n+u_right,l)=-1/de(n+u_right)*(ne(n,right)*theta(n,l)+ ne(n+t_right,left)*theta(n+t_right,l) )        
-
-          grad_e(n+u_lup,l)=-1/de(n+u_lup)*(ne(n,lup)*theta(n,l)+ ne(n+t_lup,rdown)*theta(n+t_lup,l ))       
-    
-          grad_e(n+u_ldown,l)=-1/de(n+u_ldown)*(ne(n,ldown)*theta(n,l)+ne(n+t_ldown,rup)*theta(n+t_ldown,l) )
-
-        ENDDO
+          grad_e(ij+u_right,l)=-1/de(ij+u_right)*(ne(ij,right)*theta(ij,l)+ ne(ij+t_right,left)*theta(ij+t_right,l) )        
+
+          grad_e(ij+u_lup,l)=-1/de(ij+u_lup)*(ne(ij,lup)*theta(ij,l)+ ne(ij+t_lup,rdown)*theta(ij+t_lup,l ))       
+    
+          grad_e(ij+u_ldown,l)=-1/de(ij+u_ldown)*(ne(ij,ldown)*theta(ij,l)+ne(ij+t_ldown,rup)*theta(ij+t_ldown,l) )
+
       ENDDO
     ENDDO
@@ -801,23 +780,19 @@
     
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          divgrad_i(n,l)=1./Ai(n)*(ne(n,right)*grad_e(n+u_right,l)*le(n+u_right)  +  &
-                             ne(n,rup)*grad_e(n+u_rup,l)*le(n+u_rup)              +  &  
-                             ne(n,lup)*grad_e(n+u_lup,l)*le(n+u_lup)              +  &  
-                             ne(n,left)*grad_e(n+u_left,l)*le(n+u_left)           +  &  
-                             ne(n,ldown)*grad_e(n+u_ldown,l)*le(n+u_ldown)        +  &  
-                             ne(n,rdown)*grad_e(n+u_rdown,l)*le(n+u_rdown))
-        ENDDO
+!$SIMD
+      DO ij=jj_begin,ij_end
+
+          divgrad_i(ij,l)=1./Ai(ij)*(ne(ij,right)*grad_e(ij+u_right,l)*le(ij+u_right)  +  &
+                             ne(ij,rup)*grad_e(ij+u_rup,l)*le(ij+u_rup)              +  &  
+                             ne(ij,lup)*grad_e(ij+u_lup,l)*le(ij+u_lup)              +  &  
+                             ne(ij,left)*grad_e(ij+u_left,l)*le(ij+u_left)           +  &  
+                             ne(ij,ldown)*grad_e(ij+u_ldown,l)*le(ij+u_ldown)        +  &  
+                             ne(ij,rdown)*grad_e(ij+u_rdown,l)*le(ij+u_rdown))
       ENDDO
     ENDDO
    
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-          divgrad_i(n,l)=-divgrad_i(n,l)*cdivgrad
-        ENDDO
+      DO ij=ij_begin,ij_end
+          divgrad_i(ij,l)=-divgrad_i(ij,l)*cdivgrad
       ENDDO
     ENDDO
@@ -835,48 +810,42 @@
     REAL(rstd) :: rot_v(2*iim*jjm,llb:lle)
 
-    INTEGER :: i,j,n,l
+    INTEGER :: ij,l
      
     DO l=llb,lle
-      DO j=jj_begin-1,jj_end+1
-        DO i=ii_begin-1,ii_end+1
-          n=(j-1)*iim+i
-        
-          rot_v(n+z_up,l)= 1./Av(n+z_up)*(  ne(n,rup)*ue(n+u_rup,l)*de(n+u_rup)                     &
-                                + ne(n+t_rup,left)*ue(n+t_rup+u_left,l)*de(n+t_rup+u_left)          &
-                                - ne(n,lup)*ue(n+u_lup,l)*de(n+u_lup) ) 
+!$SIMD
+      DO ij=ij_begin_ext,ij_end_ext
+        
+          rot_v(ij+z_up,l)= 1./Av(ij+z_up)*(  ne(ij,rup)*ue(ij+u_rup,l)*de(ij+u_rup)                     &
+                                + ne(ij+t_rup,left)*ue(ij+t_rup+u_left,l)*de(ij+t_rup+u_left)          &
+                                - ne(ij,lup)*ue(ij+u_lup,l)*de(ij+u_lup) ) 
                               
-          rot_v(n+z_down,l) = 1./Av(n+z_down)*( ne(n,ldown)*ue(n+u_ldown,l)*de(n+u_ldown)                 &
-                                     + ne(n+t_ldown,right)*ue(n+t_ldown+u_right,l)*de(n+t_ldown+u_right)  &
-                                     - ne(n,rdown)*ue(n+u_rdown,l)*de(n+u_rdown) )
-  
-        ENDDO
+          rot_v(ij+z_down,l) = 1./Av(ij+z_down)*( ne(ij,ldown)*ue(ij+u_ldown,l)*de(ij+u_ldown)                 &
+                                     + ne(ij+t_ldown,right)*ue(ij+t_ldown+u_right,l)*de(ij+t_ldown+u_right)  &
+                                     - ne(ij,rdown)*ue(ij+u_rdown,l)*de(ij+u_rdown) )
+  
       ENDDO
     ENDDO                              
   
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-  
-          gradrot_e(n+u_right,l)=1/le(n+u_right)*ne(n,right)*(rot_v(n+z_rdown,l)-rot_v(n+z_rup,l)) 
+!$SIMD
+      DO ij=ij_begin,ij_end
+  
+          gradrot_e(ij+u_right,l)=1/le(ij+u_right)*ne(ij,right)*(rot_v(ij+z_rdown,l)-rot_v(ij+z_rup,l)) 
          
-          gradrot_e(n+u_lup,l)=1/le(n+u_lup)*ne(n,lup)*(rot_v(n+z_up,l)-rot_v(n+z_lup,l)) 
-        
-          gradrot_e(n+u_ldown,l)=1/le(n+u_ldown)*ne(n,ldown)*(rot_v(n+z_ldown,l)-rot_v(n+z_down,l))
-        
-        ENDDO
+          gradrot_e(ij+u_lup,l)=1/le(ij+u_lup)*ne(ij,lup)*(rot_v(ij+z_up,l)-rot_v(ij+z_lup,l)) 
+        
+          gradrot_e(ij+u_ldown,l)=1/le(ij+u_ldown)*ne(ij,ldown)*(rot_v(ij+z_ldown,l)-rot_v(ij+z_down,l))
+        
       ENDDO
     ENDDO
 
     DO l=llb,lle
-      DO j=jj_begin,jj_end
-        DO i=ii_begin,ii_end
-          n=(j-1)*iim+i
-    
-          gradrot_e(n+u_right,l)=-gradrot_e(n+u_right,l)*cgradrot       
-          gradrot_e(n+u_lup,l)=-gradrot_e(n+u_lup,l)*cgradrot
-          gradrot_e(n+u_ldown,l)=-gradrot_e(n+u_ldown,l)*cgradrot
-        
-        ENDDO
+!$SIMD
+      DO ij=ij_begin,ij_end
+    
+          gradrot_e(ij+u_right,l)=-gradrot_e(ij+u_right,l)*cgradrot       
+          gradrot_e(ij+u_lup,l)=-gradrot_e(ij+u_lup,l)*cgradrot
+          gradrot_e(ij+u_ldown,l)=-gradrot_e(ij+u_ldown,l)*cgradrot
+        
       ENDDO
     ENDDO  
Index: codes/icosagcm/trunk/src/domain.f90
===================================================================
--- codes/icosagcm/trunk/src/domain.f90	(revision 173)
+++ codes/icosagcm/trunk/src/domain.f90	(revision 174)
@@ -72,7 +72,10 @@
   USE grid_param
   USE mpipara
+  USE IOIPSL
   IMPLICIT NONE
   INTEGER :: ind,nf,ni,nj,i,j
   INTEGER :: quotient, rest
+  INTEGER :: halo_i,halo_j
+  
   TYPE(t_domain),POINTER :: d
  
@@ -81,4 +84,9 @@
     ALLOCATE(domglo_rank(ndomain_glo))
     ALLOCATE(domglo_loc_ind(ndomain_glo))
+   
+    halo_i=0
+    CALL getin("halo_i",halo_i)
+    halo_j=1
+    CALL getin("halo_j",halo_j)
 
     ind=0
@@ -122,8 +130,8 @@
 
 
-          d%iim=d%ii_nb+2*halo
-          d%jjm=d%jj_nb+2*halo
-          d%ii_begin=halo+1
-          d%jj_begin=halo+1
+          d%iim=d%ii_nb+2*halo  + halo_i*2
+          d%jjm=d%jj_nb+2*halo  + halo_j*2
+          d%ii_begin=halo+1  + halo_i
+          d%jj_begin=halo+1  + halo_j
           d%ii_end=d%ii_begin+d%ii_nb-1
           d%jj_end=d%jj_begin+d%jj_nb-1
Index: codes/icosagcm/trunk/src/domain_param.f90
===================================================================
--- codes/icosagcm/trunk/src/domain_param.f90	(revision 173)
+++ codes/icosagcm/trunk/src/domain_param.f90	(revision 174)
@@ -3,5 +3,5 @@
 INTEGER :: nsplit_i
 INTEGER :: nsplit_j
-INTEGER :: halo=3
+INTEGER :: halo=1
 
 INTEGER, PARAMETER :: default_nsplit_i=1
Index: codes/icosagcm/trunk/src/timeloop_gcm.f90
===================================================================
--- codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 173)
+++ codes/icosagcm/trunk/src/timeloop_gcm.f90	(revision 174)
@@ -336,9 +336,7 @@
              ps=f_ps(ind) ; dps=f_dps(ind) ;              
              IF (omp_first) THEN
-                DO j=jj_begin,jj_end
-                   DO i=ii_begin,ii_end
-                      ij=(j-1)*iim+i
-                      ps(ij)=ps(ij)+dt*dps(ij)
-                   ENDDO
+!$SIMD
+                DO ij=ij_begin,ij_end
+                    ps(ij)=ps(ij)+dt*dps(ij)
                 ENDDO
              ENDIF 
@@ -346,9 +344,7 @@
              mass=f_mass(ind) ; dmass=f_dmass(ind) ;              
              DO l=1,llm
-                DO j=jj_begin,jj_end
-                   DO i=ii_begin,ii_end
-                      ij=(j-1)*iim+i
-                      mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
-                   ENDDO
+!$SIMD
+                DO ij=ij_begin,ij_end
+                    mass(ij,l)=mass(ij,l)+dt*dmass(ij,l)
                 ENDDO
              END DO
@@ -364,12 +360,10 @@
 
        DO l=ll_begin,ll_end
-         DO j=jj_begin,jj_end
-           DO i=ii_begin,ii_end
-             ij=(j-1)*iim+i
+!$SIMD
+         DO ij=ij_begin,ij_end
              u(ij+u_right,l)=u(ij+u_right,l)+dt*du(ij+u_right,l)
              u(ij+u_lup,l)=u(ij+u_lup,l)+dt*du(ij+u_lup,l)
              u(ij+u_ldown,l)=u(ij+u_ldown,l)+dt*du(ij+u_ldown,l)
              theta_rhodz(ij,l)=theta_rhodz(ij,l)+dt*dtheta_rhodz(ij,l)
-           ENDDO
          ENDDO
        ENDDO
@@ -401,18 +395,14 @@
                
                IF (stage==1) THEN ! first stage : save model state in XXm1
-                  DO j=jj_begin,jj_end
-                     DO i=ii_begin,ii_end
-                        ij=(j-1)*iim+i
-                        psm1(ij)=ps(ij)
-                     ENDDO
-                  ENDDO
+!$SIMD
+                 DO ij=ij_begin,ij_end
+                   psm1(ij)=ps(ij)
+                 ENDDO
                ENDIF
             
                ! updates are of the form x1 := x0 + tau*f(x1)
-               DO j=jj_begin,jj_end
-                  DO i=ii_begin,ii_end
-                     ij=(j-1)*iim+i
-                     ps(ij)=psm1(ij)+tau*dps(ij)
-                  ENDDO
+!$SIMD
+               DO ij=ij_begin,ij_end
+                   ps(ij)=psm1(ij)+tau*dps(ij)
                ENDDO
             ENDDO
@@ -428,10 +418,8 @@
             IF (stage==1) THEN ! first stage : save model state in XXm1
                DO l=ll_begin,ll_end
-                  DO j=jj_begin,jj_end
-                     DO i=ii_begin,ii_end
-                        ij=(j-1)*iim+i
-                        massm1(ij,l)=mass(ij,l)
-                     ENDDO
-                  ENDDO
+!$SIMD
+                 DO ij=ij_begin,ij_end
+                    massm1(ij,l)=mass(ij,l)
+                 ENDDO
                ENDDO
             END IF
@@ -439,9 +427,7 @@
             ! updates are of the form x1 := x0 + tau*f(x1)
             DO l=ll_begin,ll_end
-               DO j=jj_begin,jj_end
-                  DO i=ii_begin,ii_end
-                     ij=(j-1)*iim+i
-                     mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
-                  ENDDO
+!$SIMD
+               DO ij=ij_begin,ij_end
+                   mass(ij,l)=massm1(ij,l)+tau*dmass(ij,l)
                ENDDO
             ENDDO
@@ -462,12 +448,10 @@
          IF (stage==1) THEN ! first stage : save model state in XXm1
            DO l=ll_begin,ll_end
-             DO j=jj_begin,jj_end
-               DO i=ii_begin,ii_end
-                 ij=(j-1)*iim+i
+!$SIMD
+                DO ij=ij_begin,ij_end
                  um1(ij+u_right,l)=u(ij+u_right,l)
                  um1(ij+u_lup,l)=u(ij+u_lup,l)
                  um1(ij+u_ldown,l)=u(ij+u_ldown,l)
                  theta_rhodzm1(ij,l)=theta_rhodz(ij,l)
-               ENDDO
              ENDDO
            ENDDO
@@ -475,12 +459,10 @@
 
          DO l=ll_begin,ll_end
-           DO j=jj_begin,jj_end
-             DO i=ii_begin,ii_end
-               ij=(j-1)*iim+i
+!$SIMD
+           DO ij=ij_begin,ij_end
                u(ij+u_right,l)=um1(ij+u_right,l)+tau*du(ij+u_right,l)
                u(ij+u_lup,l)=um1(ij+u_lup,l)+tau*du(ij+u_lup,l)
                u(ij+u_ldown,l)=um1(ij+u_ldown,l)+tau*du(ij+u_ldown,l)
                theta_rhodz(ij,l)=theta_rhodzm1(ij,l)+tau*dtheta_rhodz(ij,l)
-             ENDDO
            ENDDO
          ENDDO
@@ -564,11 +546,9 @@
 
        DO l=ll_begin,ll_end
-         DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
+!$SIMD
+         DO ij=ij_begin_ext,ij_end_ext
              hfluxt(ij+u_right,l) = tau*hflux(ij+u_right,l)
              hfluxt(ij+u_lup,l) = tau*hflux(ij+u_lup,l)
              hfluxt(ij+u_ldown,l) = tau*hflux(ij+u_ldown,l)
-           ENDDO
          ENDDO
        ENDDO
@@ -576,9 +556,7 @@
        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
           DO l=ll_begin,ll_endp1
-             DO j=jj_begin,jj_end
-                DO i=ii_begin,ii_end
-                   ij=(j-1)*iim+i
-                   wfluxt(ij,l) = tau*wflux(ij,l)
-                ENDDO
+!$SIMD
+             DO ij=ij_begin,ij_end
+                wfluxt(ij,l) = tau*wflux(ij,l)
              ENDDO
           ENDDO
@@ -588,11 +566,9 @@
 
        DO l=ll_begin,ll_end
-         DO j=jj_begin-1,jj_end+1
-           DO i=ii_begin-1,ii_end+1
-             ij=(j-1)*iim+i
+!$SIMD
+         DO ij=ij_begin_ext,ij_end_ext
              hfluxt(ij+u_right,l) = hfluxt(ij+u_right,l)+tau*hflux(ij+u_right,l)
              hfluxt(ij+u_lup,l) = hfluxt(ij+u_lup,l)+tau*hflux(ij+u_lup,l)
              hfluxt(ij+u_ldown,l) = hfluxt(ij+u_ldown,l)+tau*hflux(ij+u_ldown,l)
-           ENDDO
          ENDDO
        ENDDO
@@ -600,9 +576,7 @@
        IF(caldyn_eta==eta_mass) THEN ! no need for vertical fluxes if eta_lag
           DO l=ll_begin,ll_endp1
-             DO j=jj_begin,jj_end
-                DO i=ii_begin,ii_end
-                   ij=(j-1)*iim+i
+!$SIMD
+             DO ij=ij_begin,ij_end
                    wfluxt(ij,l) = wfluxt(ij,l)+tau*wflux(ij,l)
-                ENDDO
              ENDDO
           ENDDO
@@ -613,34 +587,34 @@
   END SUBROUTINE accumulate_fluxes
   
-  FUNCTION maxval_i(p)
-    USE icosa
-    IMPLICIT NONE
-    REAL(rstd), DIMENSION(iim*jjm) :: p
-    REAL(rstd) :: maxval_i
-    INTEGER :: j, ij
-    
-    maxval_i=p((jj_begin-1)*iim+ii_begin)
-    
-    DO j=jj_begin-1,jj_end+1
-       ij=(j-1)*iim
-       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
-    END DO
-  END FUNCTION maxval_i
-
-  FUNCTION maxval_ik(p)
-    USE icosa
-    IMPLICIT NONE
-    REAL(rstd) :: p(iim*jjm, llm)
-    REAL(rstd) :: maxval_ik(llm)
-    INTEGER :: l,j, ij
-    
-    DO l=1,llm
-       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
-       DO j=jj_begin-1,jj_end+1
-          ij=(j-1)*iim
-          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
-       END DO
-    END DO
-  END FUNCTION maxval_ik
+!  FUNCTION maxval_i(p)
+!    USE icosa
+!    IMPLICIT NONE
+!    REAL(rstd), DIMENSION(iim*jjm) :: p
+!    REAL(rstd) :: maxval_i
+!    INTEGER :: j, ij
+!    
+!    maxval_i=p((jj_begin-1)*iim+ii_begin)
+!    
+!    DO j=jj_begin-1,jj_end+1
+!       ij=(j-1)*iim
+!       maxval_i = MAX(maxval_i, MAXVAL(p(ij+ii_begin:ij+ii_end)))
+!    END DO
+!  END FUNCTION maxval_i
+
+!  FUNCTION maxval_ik(p)
+!    USE icosa
+!    IMPLICIT NONE
+!    REAL(rstd) :: p(iim*jjm, llm)
+!    REAL(rstd) :: maxval_ik(llm)
+!    INTEGER :: l,j, ij
+!    
+!    DO l=1,llm
+!       maxval_ik(l)=p((jj_begin-1)*iim+ii_begin,l)
+!       DO j=jj_begin-1,jj_end+1
+!          ij=(j-1)*iim
+!          maxval_ik(l) = MAX(maxval_ik(l), MAXVAL(p(ij+ii_begin:ij+ii_end,l)))
+!       END DO
+!    END DO
+!  END FUNCTION maxval_ik
 
 END MODULE timeloop_gcm_mod
