Index: codes/icosagcm/trunk/src/advect.f90
===================================================================
--- codes/icosagcm/trunk/src/advect.f90	(revision 151)
+++ 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
