Index: codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_dmass.k90	(revision 658)
@@ -3,4 +3,5 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
+      !DIR$ SIMD
       DO l = 1, llm
          convm(l,ij) = mass_dbk(l,ij) * dmass_col(ij)
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_fast.k90	(revision 658)
@@ -6,4 +6,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
             berni(l,ij) = pk(l,ij)
@@ -16,4 +17,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
             berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
@@ -24,4 +26,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
             berni(l,ij) = .5*(geopot(l,ij)+geopot(l+1,ij))
@@ -39,4 +42,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          due = .5*(theta(l,ij_left,1)+theta(l,ij_right,1))*(pk(l,ij_right)-pk(l,ij_left)) + berni(l,ij_right)-berni(l,ij_left)
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_mil.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_mil.k90	(revision 658)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_mil.k90	(revision 658)
@@ -0,0 +1,16 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- caldyn_mil ----------------------------------
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      l=1
+      m_il(l,ij) = .5*rhodz(l,ij)
+      !DIR$ SIMD
+      DO l = 2, llm
+         m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij))
+      END DO
+      l=llm+1
+      m_il(l,ij) = .5*rhodz(l-1,ij)
+   END DO
+   !$OMP END DO
+   !---------------------------- caldyn_mil ----------------------------------
+   !--------------------------------------------------------------------------
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90	(revision 658)
@@ -7,4 +7,5 @@
       l=1
       w_il(l,ij) = 2.*W(l,ij)/rhodz(kup,ij)
+      !DIR$ SIMD
       DO l = 2, llm
          kdown = l-1
@@ -34,4 +35,5 @@
       v_el(l,edge) = .5*le_de(edge)*(u(kup,edge)+u(kdown,edge)) ! checked
       G_el(l,edge) = v_el(l,edge)*W_el - DePhil(l,edge)*W2_el
+      !DIR$ SIMD
       DO l = 2, llm
          kdown = l-1
@@ -64,18 +66,143 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
-      DO l = 1, llm+1
-         gPhi2=0.
-         dP=0.
-         divG=0
-         DO iedge = 1, primal_deg(ij)
-            edge = primal_edge(iedge,ij)
-            gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2
-            dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge)
-            divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de
-         END DO
-         gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
-         dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
-         dW(l,ij) = (-1./Ai(ij))*divG
-      END DO
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm+1
+            gPhi2=0.
+            dP=0.
+            divG=0
+            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
+            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
+            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
+            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
+            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
+            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
+            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
+            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
+            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
+            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
+            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
+            dW(l,ij) = (-1./Ai(ij))*divG
+         END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         !DIR$ SIMD
+         DO l = 1, llm+1
+            gPhi2=0.
+            dP=0.
+            divG=0
+            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
+            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
+            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
+            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
+            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
+            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
+            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
+            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
+            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2
+            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5)
+            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de
+            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
+            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
+            dW(l,ij) = (-1./Ai(ij))*divG
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         le_de6 = le_de(edge6)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         sign6 = primal_ne(6,ij)
+         !DIR$ SIMD
+         DO l = 1, llm+1
+            gPhi2=0.
+            dP=0.
+            divG=0
+            gPhi2 = gPhi2 + le_de1*DePhil(l,edge1)**2
+            dP = dP + le_de1*DePhil(l,edge1)*v_el(l,edge1)
+            divG = divG + sign1*G_el(l,edge1) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de2*DePhil(l,edge2)**2
+            dP = dP + le_de2*DePhil(l,edge2)*v_el(l,edge2)
+            divG = divG + sign2*G_el(l,edge2) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de3*DePhil(l,edge3)**2
+            dP = dP + le_de3*DePhil(l,edge3)*v_el(l,edge3)
+            divG = divG + sign3*G_el(l,edge3) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de4*DePhil(l,edge4)**2
+            dP = dP + le_de4*DePhil(l,edge4)*v_el(l,edge4)
+            divG = divG + sign4*G_el(l,edge4) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de5*DePhil(l,edge5)**2
+            dP = dP + le_de5*DePhil(l,edge5)*v_el(l,edge5)
+            divG = divG + sign5*G_el(l,edge5) ! -div(G_el), G_el already has le_de
+            gPhi2 = gPhi2 + le_de6*DePhil(l,edge6)**2
+            dP = dP + le_de6*DePhil(l,edge6)*v_el(l,edge6)
+            divG = divG + sign6*G_el(l,edge6) ! -div(G_el), G_el already has le_de
+            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
+            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
+            dW(l,ij) = (-1./Ai(ij))*divG
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm+1
+            gPhi2=0.
+            dP=0.
+            divG=0
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               gPhi2 = gPhi2 + le_de(edge)*DePhil(l,edge)**2
+               dP = dP + le_de(edge)*DePhil(l,edge)*v_el(l,edge)
+               divG = divG + primal_ne(iedge,ij)*G_el(l,edge) ! -div(G_el), G_el already has le_de
+            END DO
+            gradPhi2(l,ij) = 1./(2.*Ai(ij)) * gPhi2
+            dPhi(l,ij) = gradPhi2(l,ij)*w_il(l,ij) - 1./(2.*Ai(ij))*dP
+            dW(l,ij) = (-1./Ai(ij))*divG
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
@@ -85,12 +212,80 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
-      DO l = 1, llm
-         u2=0.
-         DO iedge = 1, primal_deg(ij)
-            edge = primal_edge(iedge,ij)
-            u2 = u2 + le_de(edge)*u(l,edge)**2
-         END DO
-         berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
-      END DO
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         !DIR$ SIMD
+         DO l = 1, llm
+            u2=0.
+            u2 = u2 + le_de1*u(l,edge1)**2
+            u2 = u2 + le_de2*u(l,edge2)**2
+            u2 = u2 + le_de3*u(l,edge3)**2
+            u2 = u2 + le_de4*u(l,edge4)**2
+            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
+         END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         !DIR$ SIMD
+         DO l = 1, llm
+            u2=0.
+            u2 = u2 + le_de1*u(l,edge1)**2
+            u2 = u2 + le_de2*u(l,edge2)**2
+            u2 = u2 + le_de3*u(l,edge3)**2
+            u2 = u2 + le_de4*u(l,edge4)**2
+            u2 = u2 + le_de5*u(l,edge5)**2
+            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         le_de6 = le_de(edge6)
+         !DIR$ SIMD
+         DO l = 1, llm
+            u2=0.
+            u2 = u2 + le_de1*u(l,edge1)**2
+            u2 = u2 + le_de2*u(l,edge2)**2
+            u2 = u2 + le_de3*u(l,edge3)**2
+            u2 = u2 + le_de4*u(l,edge4)**2
+            u2 = u2 + le_de5*u(l,edge5)**2
+            u2 = u2 + le_de6*u(l,edge6)**2
+            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            u2=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               u2 = u2 + le_de(edge)*u(l,edge)**2
+            END DO
+            berni(l,ij) = 1./(4.*Ai(ij)) * u2 - .25*( gradPhi2(l,ij)*w_il(l,ij)**2 + gradPhi2(l+1,ij)*w_il(l+1,ij)**2 )
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
@@ -99,4 +294,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          ! Compute mass flux and grad(berni)
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_slow_hydro.k90	(revision 658)
@@ -5,4 +5,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          uu = .5*(rhodz(l,ij_left)+rhodz(l,ij_right))*u(l,edge)
@@ -13,12 +14,80 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
-      DO l = 1, llm
-         ke=0.d0
-         DO iedge = 1, primal_deg(ij)
-            edge = primal_edge(iedge,ij)
-            ke = ke + le_de(edge)*u(l,edge)**2
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         !DIR$ SIMD
+         DO l = 1, llm
+            ke=0.d0
+            ke = ke + le_de1*u(l,edge1)**2
+            ke = ke + le_de2*u(l,edge2)**2
+            ke = ke + le_de3*u(l,edge3)**2
+            ke = ke + le_de4*u(l,edge4)**2
+            BERNI(l,ij)=ke*(.25/Ai(ij))
          END DO
-         BERNI(l,ij)=ke*(.25/Ai(ij))
-      END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         !DIR$ SIMD
+         DO l = 1, llm
+            ke=0.d0
+            ke = ke + le_de1*u(l,edge1)**2
+            ke = ke + le_de2*u(l,edge2)**2
+            ke = ke + le_de3*u(l,edge3)**2
+            ke = ke + le_de4*u(l,edge4)**2
+            ke = ke + le_de5*u(l,edge5)**2
+            BERNI(l,ij)=ke*(.25/Ai(ij))
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         le_de6 = le_de(edge6)
+         !DIR$ SIMD
+         DO l = 1, llm
+            ke=0.d0
+            ke = ke + le_de1*u(l,edge1)**2
+            ke = ke + le_de2*u(l,edge2)**2
+            ke = ke + le_de3*u(l,edge3)**2
+            ke = ke + le_de4*u(l,edge4)**2
+            ke = ke + le_de5*u(l,edge5)**2
+            ke = ke + le_de6*u(l,edge6)**2
+            BERNI(l,ij)=ke*(.25/Ai(ij))
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            ke=0.d0
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               ke = ke + le_de(edge)*u(l,edge)**2
+            END DO
+            BERNI(l,ij)=ke*(.25/Ai(ij))
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
@@ -28,4 +97,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 1, llm
             du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
@@ -38,4 +108,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 1, llm
             du(l,edge) = du(l,edge) + 1.*(berni(l,ij_left)-berni(l,ij_right)) ! minus gradient
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_solver.k90	(revision 658)
@@ -1,19 +1,4 @@
    !--------------------------------------------------------------------------
    !---------------------------- caldyn_solver ----------------------------------
-   !$OMP DO SCHEDULE(STATIC)
-   DO ij = 1, primal_num
-      l=1
-      m_il(l,ij) = .5*rhodz(l,ij)
-      DO l = 2, llm
-         m_il(l,ij) = .5*(rhodz(l,ij)+rhodz(l-1,ij))
-      END DO
-      l=llm+1
-      m_il(l,ij) = .5*rhodz(l-1,ij)
-   END DO
-   !$OMP END DO
-   !
-   IF(tau>0) THEN ! solve implicit problem for geopotential
-      CALL compute_NH_geopot(tau,PHI_BOT_VAR, rhodz, m_il, theta, W, geopot)
-   END IF
    !
    ! Compute pressure (pres) and Exner function (pk)
@@ -23,12 +8,15 @@
    gamma = 1./(1.-kappa)
    vreff = Rd*Treff/preff ! reference specific volume
-   Cvd = cpp-Rd
+   Cvd = 1./(cpp-Rd)
+   Rd_preff = kappa*cpp/preff
    SELECT CASE(caldyn_thermo)
    CASE(thermo_theta)
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
-            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
-            X_ij = (cpp/preff)*kappa*theta(l,ij,1)*rho_ij
+            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij))
+            rho_ij = rho_ij*g*rhodz(l,ij)
+            X_ij = Rd_preff*theta(l,ij,1)*rho_ij
             ! kappa.theta.rho = p/exner
             ! => X = (p/p0)/(exner/Cp)
@@ -44,7 +32,9 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
-            rho_ij = g*rhodz(l,ij)/(geopot(l+1,ij)-geopot(l,ij))
-            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))/Cvd )
+            rho_ij = 1./(geopot(l+1,ij)-geopot(l,ij))
+            rho_ij = rho_ij*g*rhodz(l,ij)
+            T_ij = Treff*exp( (theta(l,ij,1)+Rd*log(vreff*rho_ij))*Cvd )
             pres(l,ij) = rho_ij*Rd*T_ij
             pk(l,ij) = T_ij
@@ -64,4 +54,5 @@
       W(l,ij) = W(l,ij)+tau*dW(l,ij) ! update W
       dPhi(l,ij) = g*g*W(l,ij)/m_il(l,ij)
+      !DIR$ SIMD
       DO l = 2, llm
          dW(l,ij) = (1./g)*(pres(l-1,ij)-pres(l,ij)) - m_il(l,ij)
@@ -80,4 +71,5 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
+      !DIR$ SIMD
       DO l = 1, llm
          ! compute du = -0.5*g^2.grad(w^2)
@@ -90,4 +82,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          du(l,edge) = 1.*(berni(l,ij_left)-berni(l,ij_right))
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_vert.k90	(revision 658)
@@ -4,4 +4,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 2, llm
             dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) + 0.5*(theta(l,ij,iq)+theta(l-1,ij,iq))*wflux(l,ij)
@@ -11,4 +12,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm-1
             dtheta_rhodz(l,ij,iq) = dtheta_rhodz(l,ij,iq) - 0.5*(theta(l,ij,iq)+theta(l+1,ij,iq))*wflux(l+1,ij)
@@ -23,4 +25,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 2, llm
             wwuu(l,edge) = .25*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)+u(l-1,edge)) ! Fu
@@ -34,4 +37,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 1, llm
             dFu_deta = wwuu(l+1,edge)-wwuu(l,edge) ! d/deta (F*u)
@@ -46,4 +50,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 2, llm
             wwuu(l,edge) = .5*(wflux(l,ij_left)+wflux(l,ij_right))*(u(l,edge)-u(l-1,edge))
@@ -57,4 +62,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 1, llm
             du(l,edge) = du(l,edge) - (wwuu(l,edge)+wwuu(l+1,edge)) / (rhodz(l,ij_left)+rhodz(l,ij_right))
Index: codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/caldyn_vert_NH.k90	(revision 658)
@@ -9,4 +9,5 @@
       eta_dot(l,ij) = wflux_ij / mass(l,ij)
       wcov(l,ij) = w_ij*(geopot(l+1,ij)-geopot(l,ij))
+      !DIR$ SIMD
       DO l = 2, llm-1
          w_ij = .5*( W(l,ij)+W(l+1,ij) )/mass(l,ij)
@@ -29,4 +30,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          du(l,edge) = du(l,edge) - .5*(wcov(l,ij_left)+wcov(l,ij_right))*1.*(eta_dot(l,ij_right)-eta_dot(l,ij_left))
@@ -38,4 +40,5 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
+      !DIR$ SIMD
       DO l = 2, llm
          dPhi(l,ij)=dPhi(l,ij)-wflux(l,ij)*(geopot(l+1,ij)-geopot(l-1,ij))/(mass(l-1,ij)+mass(l,ij))
@@ -49,4 +52,5 @@
       l=1
       dW(l,ij) = dW(l,ij) - W_etadot(l,ij)
+      !DIR$ SIMD
       DO l = 2, llm
          dW(l,ij) = dW(l,ij) + W_etadot(l-1,ij) - W_etadot(l,ij)
Index: codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90	(revision 658)
@@ -4,4 +4,5 @@
    g2=g*g
    gm2 = 1./g2
+   vreff = Treff*cpp/preff*kappa
    gamma = 1./(1.-kappa)
    !$OMP BARRIER
@@ -35,7 +36,6 @@
             DO l = 1,llm
                rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij))
-               X_ij = Treff*exp(theta(l,ij)/cpp) ! theta = Tref.exp(s/Cp)
-               X_ij = (cpp/preff)*kappa*X_ij*rho_ij
-               p_ik(l,ij) = preff*(X_ij**gamma)
+               X_ij = log(vreff*rho_ij) + theta(l,ij)/cpp
+               p_ik(l,ij) = preff*exp(X_ij*gamma)
                c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho
                A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2
Index: codes/icosagcm/devel/src/kernels_unst/coriolis.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/coriolis.k90	(revision 658)
@@ -7,4 +7,5 @@
          ij_left = left(edge)
          ij_right = right(edge)
+         !DIR$ SIMD
          DO l = 1, llm
             Ftheta(l,edge) = .5*(theta(l,ij_left,iq)+theta(l,ij_right,iq))*hflux(l,edge)
@@ -14,12 +15,80 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
-         DO l = 1, llm
-            divF=0.
-            DO iedge = 1, primal_deg(ij)
-               edge = primal_edge(iedge,ij)
-               divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij)
-            END DO
-            dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
-         END DO
+         ! this VLOOP iterates over primal cell edges
+         SELECT CASE(primal_deg(ij))
+         CASE(4)
+            edge1 = primal_edge(1,ij)
+            edge2 = primal_edge(2,ij)
+            edge3 = primal_edge(3,ij)
+            edge4 = primal_edge(4,ij)
+            sign1 = primal_ne(1,ij)
+            sign2 = primal_ne(2,ij)
+            sign3 = primal_ne(3,ij)
+            sign4 = primal_ne(4,ij)
+            !DIR$ SIMD
+            DO l = 1, llm
+               divF=0.
+               divF = divF + Ftheta(l,edge1)*sign1
+               divF = divF + Ftheta(l,edge2)*sign2
+               divF = divF + Ftheta(l,edge3)*sign3
+               divF = divF + Ftheta(l,edge4)*sign4
+               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
+            END DO
+         CASE(5)
+            edge1 = primal_edge(1,ij)
+            edge2 = primal_edge(2,ij)
+            edge3 = primal_edge(3,ij)
+            edge4 = primal_edge(4,ij)
+            edge5 = primal_edge(5,ij)
+            sign1 = primal_ne(1,ij)
+            sign2 = primal_ne(2,ij)
+            sign3 = primal_ne(3,ij)
+            sign4 = primal_ne(4,ij)
+            sign5 = primal_ne(5,ij)
+            !DIR$ SIMD
+            DO l = 1, llm
+               divF=0.
+               divF = divF + Ftheta(l,edge1)*sign1
+               divF = divF + Ftheta(l,edge2)*sign2
+               divF = divF + Ftheta(l,edge3)*sign3
+               divF = divF + Ftheta(l,edge4)*sign4
+               divF = divF + Ftheta(l,edge5)*sign5
+               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
+            END DO
+         CASE(6)
+            edge1 = primal_edge(1,ij)
+            edge2 = primal_edge(2,ij)
+            edge3 = primal_edge(3,ij)
+            edge4 = primal_edge(4,ij)
+            edge5 = primal_edge(5,ij)
+            edge6 = primal_edge(6,ij)
+            sign1 = primal_ne(1,ij)
+            sign2 = primal_ne(2,ij)
+            sign3 = primal_ne(3,ij)
+            sign4 = primal_ne(4,ij)
+            sign5 = primal_ne(5,ij)
+            sign6 = primal_ne(6,ij)
+            !DIR$ SIMD
+            DO l = 1, llm
+               divF=0.
+               divF = divF + Ftheta(l,edge1)*sign1
+               divF = divF + Ftheta(l,edge2)*sign2
+               divF = divF + Ftheta(l,edge3)*sign3
+               divF = divF + Ftheta(l,edge4)*sign4
+               divF = divF + Ftheta(l,edge5)*sign5
+               divF = divF + Ftheta(l,edge6)*sign6
+               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
+            END DO
+         CASE DEFAULT
+            !DIR$ SIMD
+            DO l = 1, llm
+               divF=0.
+               DO iedge = 1, primal_deg(ij)
+                  edge = primal_edge(iedge,ij)
+                  divF = divF + Ftheta(l,edge)*primal_ne(iedge,ij)
+               END DO
+               dtheta_rhodz(l,ij,iq) = -divF / Ai(ij)
+            END DO
+         END SELECT
       END DO
       !$OMP END DO
@@ -28,12 +97,80 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
-      DO l = 1, llm
-         divF=0.
-         DO iedge = 1, primal_deg(ij)
-            edge = primal_edge(iedge,ij)
-            divF = divF + hflux(l,edge)*primal_ne(iedge,ij)
-         END DO
-         convm(l,ij) = -divF / Ai(ij)
-      END DO
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            divF=0.
+            divF = divF + hflux(l,edge1)*sign1
+            divF = divF + hflux(l,edge2)*sign2
+            divF = divF + hflux(l,edge3)*sign3
+            divF = divF + hflux(l,edge4)*sign4
+            convm(l,ij) = -divF / Ai(ij)
+         END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            divF=0.
+            divF = divF + hflux(l,edge1)*sign1
+            divF = divF + hflux(l,edge2)*sign2
+            divF = divF + hflux(l,edge3)*sign3
+            divF = divF + hflux(l,edge4)*sign4
+            divF = divF + hflux(l,edge5)*sign5
+            convm(l,ij) = -divF / Ai(ij)
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         sign6 = primal_ne(6,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            divF=0.
+            divF = divF + hflux(l,edge1)*sign1
+            divF = divF + hflux(l,edge2)*sign2
+            divF = divF + hflux(l,edge3)*sign3
+            divF = divF + hflux(l,edge4)*sign4
+            divF = divF + hflux(l,edge5)*sign5
+            divF = divF + hflux(l,edge6)*sign6
+            convm(l,ij) = -divF / Ai(ij)
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            divF=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               divF = divF + hflux(l,edge)*primal_ne(iedge,ij)
+            END DO
+            convm(l,ij) = -divF / Ai(ij)
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
@@ -41,12 +178,71 @@
    !$OMP DO SCHEDULE(STATIC)
    DO edge = 1, edge_num
-      DO l = 1, llm
-         du_trisk=0.
-         DO itrisk = 1, trisk_deg(edge)
-            edge_trisk = trisk(itrisk,edge)
-            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
-         END DO
-         du(l,edge) = du(l,edge) + .5*du_trisk
-      END DO
+      ! this VLOOP iterates over the TRISK stencil
+      SELECT CASE(trisk_deg(edge))
+      CASE(10)
+         !DIR$ SIMD
+         DO l = 1, llm
+            du_trisk=0.
+            itrisk = 1
+            edge_trisk = trisk(1,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 2
+            edge_trisk = trisk(2,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 3
+            edge_trisk = trisk(3,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 4
+            edge_trisk = trisk(4,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 5
+            edge_trisk = trisk(5,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 6
+            edge_trisk = trisk(6,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 7
+            edge_trisk = trisk(7,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 8
+            edge_trisk = trisk(8,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 9
+            edge_trisk = trisk(9,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 10
+            edge_trisk = trisk(10,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            du(l,edge) = du(l,edge) + .5*du_trisk
+         END DO
+      CASE(4)
+         !DIR$ SIMD
+         DO l = 1, llm
+            du_trisk=0.
+            itrisk = 1
+            edge_trisk = trisk(1,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 2
+            edge_trisk = trisk(2,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 3
+            edge_trisk = trisk(3,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            itrisk = 4
+            edge_trisk = trisk(4,edge)
+            du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            du(l,edge) = du(l,edge) + .5*du_trisk
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            du_trisk=0.
+            DO itrisk = 1, trisk_deg(edge)
+               edge_trisk = trisk(itrisk,edge)
+               du_trisk = du_trisk + wee(itrisk,edge)*hflux(l,edge_trisk)*(qu(l,edge)+qu(l,edge_trisk))
+            END DO
+            du(l,edge) = du(l,edge) + .5*du_trisk
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
Index: codes/icosagcm/devel/src/kernels_unst/div.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/div.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/div.k90	(revision 658)
@@ -3,12 +3,95 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, primal_num
-      DO l = 1, llm
-         div_ij=0.
-         DO iedge = 1, primal_deg(ij)
-            edge = primal_edge(iedge,ij)
-            div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge)
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*u(l,edge1)
+            div_ij = div_ij + sign2*le_de2*u(l,edge2)
+            div_ij = div_ij + sign3*le_de3*u(l,edge3)
+            div_ij = div_ij + sign4*le_de4*u(l,edge4)
+            divu(l,ij) = div_ij / Ai(ij)
          END DO
-         divu(l,ij) = div_ij / Ai(ij)
-      END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*u(l,edge1)
+            div_ij = div_ij + sign2*le_de2*u(l,edge2)
+            div_ij = div_ij + sign3*le_de3*u(l,edge3)
+            div_ij = div_ij + sign4*le_de4*u(l,edge4)
+            div_ij = div_ij + sign5*le_de5*u(l,edge5)
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         le_de1 = le_de(edge1)
+         le_de2 = le_de(edge2)
+         le_de3 = le_de(edge3)
+         le_de4 = le_de(edge4)
+         le_de5 = le_de(edge5)
+         le_de6 = le_de(edge6)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         sign6 = primal_ne(6,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            div_ij = div_ij + sign1*le_de1*u(l,edge1)
+            div_ij = div_ij + sign2*le_de2*u(l,edge2)
+            div_ij = div_ij + sign3*le_de3*u(l,edge3)
+            div_ij = div_ij + sign4*le_de4*u(l,edge4)
+            div_ij = div_ij + sign5*le_de5*u(l,edge5)
+            div_ij = div_ij + sign6*le_de6*u(l,edge6)
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            div_ij=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               div_ij = div_ij + primal_ne(iedge,ij)*le_de(edge)*u(l,edge)
+            END DO
+            divu(l,ij) = div_ij / Ai(ij)
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
Index: codes/icosagcm/devel/src/kernels_unst/gradient.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/gradient.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/gradient.k90	(revision 658)
@@ -5,4 +5,5 @@
       ij_left = left(edge)
       ij_right = right(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          grad(l,edge) = 1.*(b(l,ij_right)-b(l,ij_left))
Index: codes/icosagcm/devel/src/kernels_unst/pvort_only.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/pvort_only.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/pvort_only.k90	(revision 658)
@@ -3,17 +3,71 @@
    !$OMP DO SCHEDULE(STATIC)
    DO ij = 1, dual_num
-      DO l = 1, llm
-         etav = 0.d0
-         DO iedge = 1, dual_deg(ij)
-            edge = dual_edge(iedge,ij)
-            etav = etav + dual_ne(iedge,ij)*u(l,edge)
+      ! this VLOOP iterates over dual cell edges
+      SELECT CASE(dual_deg(ij))
+      CASE(3)
+         edge1 = dual_edge(1,ij)
+         edge2 = dual_edge(2,ij)
+         edge3 = dual_edge(3,ij)
+         sign1 = dual_ne(1,ij)
+         sign2 = dual_ne(2,ij)
+         sign3 = dual_ne(3,ij)
+         vertex1 = dual_vertex(1,ij)
+         vertex2 = dual_vertex(2,ij)
+         vertex3 = dual_vertex(3,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            etav = etav + sign1*u(l,edge1)
+            etav = etav + sign2*u(l,edge2)
+            etav = etav + sign3*u(l,edge3)
+            hv=0.
+            hv = hv + Riv2(1,ij)*rhodz(l,vertex1)
+            hv = hv + Riv2(2,ij)*rhodz(l,vertex2)
+            hv = hv + Riv2(3,ij)*rhodz(l,vertex3)
+            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
          END DO
-         hv=0.
-         DO ivertex = 1, dual_deg(ij)
-            vertex = dual_vertex(ivertex,ij)
-            hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex)
+      CASE(4)
+         edge1 = dual_edge(1,ij)
+         edge2 = dual_edge(2,ij)
+         edge3 = dual_edge(3,ij)
+         edge4 = dual_edge(4,ij)
+         sign1 = dual_ne(1,ij)
+         sign2 = dual_ne(2,ij)
+         sign3 = dual_ne(3,ij)
+         sign4 = dual_ne(4,ij)
+         vertex1 = dual_vertex(1,ij)
+         vertex2 = dual_vertex(2,ij)
+         vertex3 = dual_vertex(3,ij)
+         vertex4 = dual_vertex(4,ij)
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            etav = etav + sign1*u(l,edge1)
+            etav = etav + sign2*u(l,edge2)
+            etav = etav + sign3*u(l,edge3)
+            etav = etav + sign4*u(l,edge4)
+            hv=0.
+            hv = hv + Riv2(1,ij)*rhodz(l,vertex1)
+            hv = hv + Riv2(2,ij)*rhodz(l,vertex2)
+            hv = hv + Riv2(3,ij)*rhodz(l,vertex3)
+            hv = hv + Riv2(4,ij)*rhodz(l,vertex4)
+            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
          END DO
-         qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
-      END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            etav = 0.d0
+            DO iedge = 1, dual_deg(ij)
+               edge = dual_edge(iedge,ij)
+               etav = etav + dual_ne(iedge,ij)*u(l,edge)
+            END DO
+            hv=0.
+            DO ivertex = 1, dual_deg(ij)
+               vertex = dual_vertex(ivertex,ij)
+               hv = hv + Riv2(ivertex,ij)*rhodz(l,vertex)
+            END DO
+            qv(l,ij) = (etav + fv(ij)*Av(ij) )/(hv*Av(ij))
+         END DO
+      END SELECT
    END DO
    !$OMP END DO
@@ -22,4 +76,5 @@
       ij_up = up(edge)
       ij_down = down(edge)
+      !DIR$ SIMD
       DO l = 1, llm
          qu(l,edge)=0.5d0*(qv(l,ij_down)+qv(l,ij_up))
Index: codes/icosagcm/devel/src/kernels_unst/theta.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 657)
+++ codes/icosagcm/devel/src/kernels_unst/theta.k90	(revision 658)
@@ -16,4 +16,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
             ! FIXME : formula below (used in DYNAMICO) is for dak, dbk based on pressure rather than mass
@@ -28,4 +29,5 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij = 1, primal_num
+         !DIR$ SIMD
          DO l = 1, llm
             theta(l,ij,iq) = theta_rhodz(l,ij,iq)/rhodz(l,ij)
Index: codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 657)
+++ codes/icosagcm/devel/src/unstructured/caldyn_unstructured.F90	(revision 658)
@@ -21,5 +21,4 @@
 #define DECLARE_VERTICES INTEGER VERTICES
 #define PHI_BOT(ij) Phi_bot
-#define PHI_BOT_VAR 0.
 
 #define BINDC_(thename) BIND(C, name=#thename)
@@ -40,13 +39,46 @@
 #define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
 
+#define START_TRACE(id,nprimal,ndual,nedge) CALL enter_trace(id, 8*llm*((nprimal)*primal_num+(ndual)*dual_num+(nedge)*edge_num) )
+#define STOP_TRACE CALL exit_trace()
+
 !----------------------------- Non-Hydrostatic -----------------------------
 
-SUBROUTINE compute_NH_geopot(tau,dummy, m_ik, m_il, theta, W_il, Phi_il)
-  FIELD_MASS   :: m_ik, theta, p_ik, A_ik, C_ik   ! IN*2,LOCAL*3
-  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il, R_il, x_il, B_il, D_il  ! IN,INOUT*2, LOCAL*5 
-  DBL :: tau, dummy, gamma, rho_ij, X_ij, Y_ij, wil, tau2_g, g2, gm2, ml_g2, c2_mik
-  DECLARE_INDICES
+SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
+  FIELD_MASS   :: m_ik, theta   ! IN*2
+  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL*5 
+  DBL :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff
   INTEGER :: iter
+  DECLARE_INDICES
+  DBL :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2
+#define COLUMN 0
+#if COLUMN 
+  DOUBLE1(llm)  :: pk, Ak, Ck
+  DOUBLE1(llm+1):: Rl, Bl, Dl, xl
+#define p_ik(l,ij) pk(l)
+#define A_ik(l,ij) Ak(l)
+#define C_ik(l,ij) Ck(l)
+#define R_il(l,ij) Rl(l)
+#define B_il(l,ij) Bl(l)
+#define D_il(l,ij) Dl(l)
+#define x_il(l,ij) xl(l)
+#else
+  FIELD_MASS :: p_ik, A_ik, C_ik
+  FIELD_GEOPOT :: R_il, B_il, D_il, x_il
+#endif
+
+  START_TRACE(id_NH_geopot, 7,0,0)
 #include "../kernels_unst/compute_NH_geopot.k90"
+  STOP_TRACE  
+
+#if COLUMN
+#undef p_ik
+#undef A_ik
+#undef C_ik
+#undef R_il
+#undef B_il
+#undef D_il
+#undef x_il
+#endif
+#undef COLUMN
 END SUBROUTINE compute_NH_geopot
 
@@ -59,5 +91,7 @@
   DECLARE_EDGES
   DBL :: W_el, W2_el, gPhi2, dP, divG, u2, uu
+  START_TRACE(id_slow_NH, 5,0,3)
 #include "../kernels_unst/caldyn_slow_NH.k90"
+  STOP_TRACE
 END SUBROUTINE compute_caldyn_slow_NH
 
@@ -69,6 +103,12 @@
   FIELD_U      :: du                     ! OUT
   DECLARE_INDICES
-  DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff
+  DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff
+#include "../kernels_unst/caldyn_mil.k90"
+  IF(tau>0) THEN ! solve implicit problem for geopotential
+    CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
+  END IF
+  START_TRACE(id_solver, 7,0,1)
 #include "../kernels_unst/caldyn_solver.k90"
+  STOP_TRACE
 END SUBROUTINE compute_caldyn_solver
 
@@ -79,5 +119,7 @@
   DECLARE_INDICES
   DBL :: w_ij, wflux_ij
+  START_TRACE(id_vert_NH, 6,0,1)
 #include "../kernels_unst/caldyn_vert_NH.k90"
+  STOP_TRACE
 END SUBROUTINE compute_caldyn_vert_NH
 
@@ -91,5 +133,7 @@
   DECLARE_INDICES
   DBL :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix
+  START_TRACE(id_geopot, 3,0,3)
 #include "../kernels_unst/compute_geopot.k90"
+  STOP_TRACE
 END SUBROUTINE compute_geopot
 !
@@ -102,5 +146,7 @@
   LOGICAL, PARAMETER :: zero=.TRUE.
   DBL :: ke, uu
+  START_TRACE(id_slow_hydro, 3,0,3)
 #include "../kernels_unst/caldyn_slow_hydro.k90"
+  STOP_TRACE
 END SUBROUTINE compute_caldyn_slow_hydro
 
@@ -118,7 +164,9 @@
   wwuu=0.
   !$OMP BARRIER
+  START_TRACE(id_vert, 5,0,3)
 #include "../kernels_unst/caldyn_wflux.k90"
 #include "../kernels_unst/caldyn_dmass.k90"
 #include "../kernels_unst/caldyn_vert.k90"
+  STOP_TRACE
 END SUBROUTINE caldyn_vert
 
@@ -130,5 +178,7 @@
   DECLARE_EDGES
   DBL :: divF, du_trisk
+  START_TRACE(id_coriolis, 3,4,0) ! primal, dual, edge
 #include "../kernels_unst/coriolis.k90"
+  STOP_TRACE
 END SUBROUTINE
 
@@ -139,5 +189,7 @@
   DECLARE_INDICES
   DBL :: m
+  START_TRACE(id_theta, 3,0,0) ! primal, dual, edge
 #include "../kernels_unst/theta.k90"
+  STOP_TRACE
 END SUBROUTINE
 
@@ -150,5 +202,7 @@
   DECLARE_VERTICES
   DBL :: etav, hv
+  START_TRACE(id_pvort_only, 1,1,2) ! primal, dual, edge
 #include "../kernels_unst/pvort_only.k90"
+  STOP_TRACE
 END SUBROUTINE compute_pvort_only
 
@@ -162,7 +216,7 @@
   DECLARE_EDGES
   DBL          :: due
-
+  START_TRACE(id_fast, 4,0,2) ! primal, dual, edge
 #include "../kernels_unst/caldyn_fast.k90"
-
+  STOP_TRACE
 END SUBROUTINE compute_caldyn_fast
 
Index: codes/icosagcm/devel/src/unstructured/data_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 657)
+++ codes/icosagcm/devel/src/unstructured/data_unstructured.F90	(revision 658)
@@ -44,12 +44,13 @@
   INTEGER(C_INT), BIND(C) :: comm_icosa
 
-  INTEGER, PARAMETER :: id_pvort_only=1, id_slow_hydro=2, id_fast=3, id_coriolis=4, id_theta=5, &
-       id_vert_NH=6, id_solver=7, id_slow_NH=8, id_NH_geopot=9, id_vert=10, id_NH_Phi_star=11, nb_routines=11
+  INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, &
+       id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, &
+       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, nb_routines=14 
   DBL, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel
   INTEGER, PRIVATE :: current_id, nb_calls(nb_routines), bytes(nb_routines) ! bytes read or written be each kernel
   CHARACTER(len = 10) :: id_name(nb_routines) = &
-       (/'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', &
-       'vert_NH   ', 'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert      ', &
-       'Phi_star  '/)
+       (/'dev1      ', 'dev2      ', &
+       'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', &
+       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ' /)
 
 CONTAINS
Index: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90
===================================================================
--- codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 657)
+++ codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90	(revision 658)
@@ -219,4 +219,5 @@
     !
     INTEGER :: l, ij
+    CALL enter_trace(id_update, 8*sz*(j+1) )
     !  PRINT *, clj(1:j,j)
     SELECT CASE(j)
@@ -263,4 +264,5 @@
        END DO
     END SELECT
+    CALL exit_trace()
   END SUBROUTINE update
   
