Index: codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/caldyn_slow_NH.k90	(revision 614)
+++ 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)
