MODULE caldyn_unstructured_mod
   USE ISO_C_BINDING
   USE OMP_LIB
   USE earth_const
   USE data_unstructured_mod
   IMPLICIT NONE
   SAVE
   PRIVATE

   PUBLIC :: compute_NH_geopot, compute_caldyn_slow_NH, compute_caldyn_solver, &
             compute_caldyn_vert_NH, compute_geopot, compute_caldyn_slow_hydro, &
             caldyn_vert, compute_coriolis, compute_theta, compute_pvort_only, &
             compute_caldyn_fast, gradient, div, &
             compute_scalar_laplacian, compute_curl_laplacian

CONTAINS

!----------------------------- Non-Hydrostatic -----------------------------

   SUBROUTINE compute_NH_geopot(tau, m_ik, m_il, theta, W_il, Phi_il)
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: m_ik, theta   ! IN*2
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL
      REAL(C_DOUBLE) :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff
      INTEGER :: iter
      LOGICAL :: debug_hevi_solver
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: p_ik, A_ik, C_ik
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: R_il, B_il, D_il, x_il

      debug_hevi_solver = .FALSE.
!$OMP MASTER
      debug_hevi_solver = debug_hevi_solver_
!$OMP END MASTER

      CALL enter_trace(id_NH_geopot, 8*llm*((7)*primal_num + (0)*dual_num + (0)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- compute_NH_geopot ----------------------------------
      tau2_g = tau*tau/g
      g2 = g*g
      gm2 = 1./g2
      vreff = Treff*cpp/preff*kappa
      gamma = 1./(1.-kappa)
      !$OMP BARRIER
      ! compute Phi_star
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = 1, llm + 1
            Phi_star_il(l, ij) = Phi_il(l, ij) + tau*g2*(W_il(l, ij)/m_il(l, ij) - tau)
         END DO
      END DO
      !$OMP END DO
      ! Newton-Raphson iteration : Phi_il contains current guess value
      DO iter = 1, 2 ! 2 iterations should be enough
         ! Compute pressure, A_ik
         SELECT CASE (caldyn_thermo)
         CASE (thermo_theta)
            !$OMP DO SCHEDULE(STATIC)
            DO ij = 1, primal_num
               DO l = 1, llm
                  rho_ij = (g*m_ik(l, ij))/(Phi_il(l + 1, ij) - Phi_il(l, ij))
                  X_ij = (cpp/preff)*kappa*theta(l, ij)*rho_ij
                  p_ik(l, ij) = preff*(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
               END DO
            END DO
            !$OMP END DO
         CASE (thermo_entropy)
            !$OMP DO SCHEDULE(STATIC)
            DO ij = 1, primal_num
               DO l = 1, llm
                  rho_ij = (g*m_ik(l, ij))/(Phi_il(l + 1, ij) - Phi_il(l, ij))
                  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
               END DO
            END DO
            !$OMP END DO
         CASE DEFAULT
            PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
            STOP
         END SELECT
         ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
         ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            ! Compute residual R_il and B_il
            ! bottom interface l=1
            ml_g2 = gm2*m_il(1, ij)
            B_il(1, ij) = A_ik(1, ij) + ml_g2 + tau2_g*rho_bot
            R_il(1, ij) = ml_g2*(Phi_il(1, ij) - Phi_star_il(1, ij)) &
                          + tau2_g*(p_ik(1, ij) - pbot + rho_bot*(Phi_il(1, ij) - Phi_bot))
            DO l = 2, llm
               ! inner interfaces
               ml_g2 = gm2*m_il(l, ij)
               B_il(l, ij) = A_ik(l, ij) + A_ik(l - 1, ij) + ml_g2
               R_il(l, ij) = ml_g2*(Phi_il(l, ij) - Phi_star_il(l, ij)) &
                             + tau2_g*(p_ik(l, ij) - p_ik(l - 1, ij))
               ! consistency check : if Wil=0 and initial state is in hydrostatic balance
               ! then Phi_star_il(l,ij) = Phi_il(l,ij) - tau^2*g^2
               ! and residual = tau^2*(ml+(1/g)dl_pi)=0
            END DO
            ! top interface l=llm+1
            ml_g2 = gm2*m_il(llm + 1, ij)
            B_il(llm + 1, ij) = A_ik(llm + 1 - 1, ij) + ml_g2
            R_il(llm + 1, ij) = ml_g2*(Phi_il(llm + 1, ij) - Phi_star_il(llm + 1, ij)) &
                                + tau2_g*(ptop - p_ik(llm + 1 - 1, ij))
            !
            ! Forward sweep :
            ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
            ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
            X_ij = 1./B_il(1, ij)
            C_ik(1, ij) = -A_ik(1, ij)*X_ij
            D_il(1, ij) = R_il(1, ij)*X_ij
            DO l = 2, llm
               X_ij = 1./(B_il(l, ij) + A_ik(l - 1, ij)*C_ik(l - 1, ij))
               C_ik(l, ij) = -A_ik(l, ij)*X_ij
               D_il(l, ij) = (R_il(l, ij) + A_ik(l - 1, ij)*D_il(l - 1, ij))*X_ij
            END DO
            X_ij = 1./(B_il(llm + 1, ij) + A_ik(llm + 1 - 1, ij)*C_ik(llm + 1 - 1, ij))
            D_il(llm + 1, ij) = (R_il(llm + 1, ij) + A_ik(llm + 1 - 1, ij)*D_il(llm + 1 - 1, ij))*X_ij
            ! Back substitution :
            ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
            ! + Newton-Raphson update
            ! top interface l=llm+1
            x_il(llm + 1, ij) = D_il(llm + 1, ij)
            Phi_il(llm + 1, ij) = Phi_il(llm + 1, ij) - x_il(llm + 1, ij)
            DO l = llm, 1, -1
               ! Back substitution at lower interfaces
               x_il(l, ij) = D_il(l, ij) - C_ik(l, ij)*x_il(l + 1, ij)
               Phi_il(l, ij) = Phi_il(l, ij) - x_il(l, ij)
            END DO
         END DO
         !$OMP END DO
         IF (debug_hevi_solver) THEN
            PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)), MAXVAL(ABS(B_il))
            PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)), MAXVAL(ABS(D_il))
            DO l = 1, llm + 1
               WRITE (*, '(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter, l, MAXVAL(ABS(x_il(l, :)))
            END DO
            DO l = 1, llm + 1
               WRITE (*, '(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter, l, MAXVAL(ABS(R_il(l, :)))
            END DO
         END IF
      END DO ! Newton-Raphson
      !$OMP BARRIER
      debug_hevi_solver = .FALSE.
      !---------------------------- compute_NH_geopot ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()

!$OMP MASTER
      debug_hevi_solver_ = debug_hevi_solver
!$OMP END MASTER

   END SUBROUTINE compute_NH_geopot

   SUBROUTINE compute_caldyn_slow_NH(u, rhodz, Phi, W, w_il, berni, gradPhi2, DePhil, v_el, G_el, F_el, hflux, du, dPhi, dW)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: u, hflux, du   ! IN, OUT, OUT
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: rhodz, berni   ! IN, BUF
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: Phi, W, dPhi, dW, w_il, gradPhi2  ! IN,IN, OUT,OUT, BUF*2
      REAL(C_DOUBLE), DIMENSION(llm + 1, edge_num)     :: DePhil, v_el, G_el, F_el ! BUF*4
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      REAL(C_DOUBLE) :: W_el, W2_el, gPhi2, dP, divG, u2, uu
      CALL enter_trace(id_slow_NH, 8*llm*((5)*primal_num + (0)*dual_num + (3)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_slow_NH ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         kdown = 1
         kup = 1
         l = 1
         w_il(l, ij) = 2.*W(l, ij)/rhodz(kup, ij)
         !DIR$ SIMD
         DO l = 2, llm
            kdown = l - 1
            kup = l
            w_il(l, ij) = 2.*W(l, ij)/(rhodz(kdown, ij) + rhodz(kup, ij))
         END DO
         kdown = llm
         kup = llm
         l = llm + 1
         w_il(l, ij) = 2.*W(l, ij)/rhodz(kdown, ij)
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         kdown = 1
         kup = 1
         l = 1
         ! compute DePhi, v_el, G_el, F_el
         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
         ! while DePhil, W_el and F_el do not
         W_el = .5*(W(l, ij_right) + W(l, ij_left))
         DePhil(l, edge) = 1.*(Phi(l, ij_right) - Phi(l, ij_left))
         F_el(l, edge) = DePhil(l, edge)*W_el
         W2_el = .5*le_de(edge)*(W(l, ij_left)*w_il(l, ij_left) + W(l, ij_right)*w_il(l, ij_right))
         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
            kup = l
            ! compute DePhi, v_el, G_el, F_el
            ! v_el, W2_el and therefore G_el incorporate metric factor le_de
            ! while DePhil, W_el and F_el do not
            W_el = .5*(W(l, ij_right) + W(l, ij_left))
            DePhil(l, edge) = 1.*(Phi(l, ij_right) - Phi(l, ij_left))
            F_el(l, edge) = DePhil(l, edge)*W_el
            W2_el = .5*le_de(edge)*(W(l, ij_left)*w_il(l, ij_left) + W(l, ij_right)*w_il(l, ij_right))
            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
         END DO
         kdown = llm
         kup = llm
         l = llm + 1
         ! compute DePhi, v_el, G_el, F_el
         ! v_el, W2_el and therefore G_el incorporate metric factor le_de
         ! while DePhil, W_el and F_el do not
         W_el = .5*(W(l, ij_right) + W(l, ij_left))
         DePhil(l, edge) = 1.*(Phi(l, ij_right) - Phi(l, ij_left))
         F_el(l, edge) = DePhil(l, edge)*W_el
         W2_el = .5*le_de(edge)*(W(l, ij_left)*w_il(l, ij_left) + W(l, ij_right)*w_il(l, ij_right))
         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
      END DO
      !$OMP END DO
      ! compute GradPhi2, dPhi, dW
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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
      ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
      !$OMP BARRIER
      ! Compute berni at scalar points
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_left = left(edge)
         ij_right = right(edge)
         !DIR$ SIMD
         DO l = 1, llm
            ! Compute mass flux and grad(berni)
            uu = .5*(rhodz(l, ij_left) + rhodz(l, ij_right))*u(l, edge) - .5*(F_el(l, edge) + F_el(l + 1, edge))
            hflux(l, edge) = le_de(edge)*uu
            du(l, edge) = 1.*(berni(l, ij_left) - berni(l, ij_right))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- caldyn_slow_NH ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_slow_NH

   SUBROUTINE compute_caldyn_solver(tau, rhodz, theta, berni, pres, m_il, pk, geopot, W, dPhi, dW, du)
      USE grid_param
      REAL(C_DOUBLE), INTENT(IN) :: tau
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: rhodz, pk, berni, pres    ! IN, OUT, BUF*2
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn)  :: theta                  ! IN
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: geopot, W, dPhi, dW, m_il ! INOUT,INOUT, OUT,OUT, BUF
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: du                     ! OUT
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff

      !--------------------------------------------------------------------------
      !---------------------------- 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 ----------------------------------
      !--------------------------------------------------------------------------
      IF (tau > 0) THEN ! solve implicit problem for geopotential
         CALL compute_NH_geopot(tau, rhodz, m_il, theta, W, geopot)
      END IF
      CALL enter_trace(id_solver, 8*llm*((7)*primal_num + (0)*dual_num + (1)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_solver ----------------------------------
      !
      ! Compute pressure (pres) and Exner function (pk)
      ! kappa = R/Cp
      ! 1-kappa = Cv/Cp
      ! Cp/Cv = 1/(1-kappa)
      gamma = 1./(1.-kappa)
      vreff = Rd*Treff/preff ! reference specific volume
      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 = 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)
               ! = (p/p0)^(1-kappa)
               pres(l, ij) = preff*(X_ij**gamma) ! pressure
               ! Compute Exner function (needed by compute_caldyn_fast)
               ! other formulae possible if exponentiation is slow
               pk(l, ij) = cpp*((pres(l, ij)/preff)**kappa) ! Exner
            END DO
         END DO
         !$OMP END DO
      CASE (thermo_entropy)
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               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
            END DO
         END DO
         !$OMP END DO
      CASE DEFAULT
         STOP
      END SELECT
      ! We need a barrier here because we compute pres above and do a vertical difference below
      !$OMP BARRIER
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         l = 1
         ! Lower BC
         dW(l, ij) = (1./g)*(pbot - rho_bot*(geopot(l, ij) - Phi_bot) - pres(l, ij)) - m_il(l, ij)
         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)
            W(l, ij) = W(l, ij) + tau*dW(l, ij) ! update W
            dPhi(l, ij) = g*g*W(l, ij)/m_il(l, ij)
         END DO
         l = llm + 1
         ! Top BC
         dW(l, ij) = (1./g)*(pres(l - 1, ij) - ptop) - m_il(l, ij)
         W(l, ij) = W(l, ij) + tau*dW(l, ij) ! update W
         dPhi(l, ij) = g*g*W(l, ij)/m_il(l, ij)
      END DO
      !$OMP END DO
      ! We need a barrier here because we update W above and do a vertical average below
      !$OMP BARRIER
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            ! compute du = -0.5*g^2.grad(w^2)
            berni(l, ij) = (-.25*g*g)*((W(l, ij)/m_il(l, ij))**2 + (W(l + 1, ij)/m_il(l + 1, ij))**2)
         END DO
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- caldyn_solver ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_solver

   SUBROUTINE compute_caldyn_vert_NH(mass, geopot, W, wflux, eta_dot, wcov, W_etadot, du, dPhi, dW)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: geopot, W, wflux, dPhi, dW ! IN*3, INOUT*2
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: du ! INOUT
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: w_ij, wflux_ij
      CALL enter_trace(id_vert_NH, 8*llm*((6)*primal_num + (0)*dual_num + (1)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_vert_NH ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         l = 1
         w_ij = (W(l, ij) + .5*W(l + 1, ij))/mass(l, ij)
         wflux_ij = .5*(wflux(l, ij) + wflux(l + 1, ij))
         W_etadot(l, ij) = wflux_ij*w_ij
         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)
            wflux_ij = .5*(wflux(l, ij) + wflux(l + 1, ij))
            W_etadot(l, ij) = wflux_ij*w_ij
            eta_dot(l, ij) = wflux_ij/mass(l, ij)
            wcov(l, ij) = w_ij*(geopot(l + 1, ij) - geopot(l, ij))
         END DO
         l = llm
         w_ij = (.5*W(l, ij) + W(l + 1, ij))/mass(l, ij)
         wflux_ij = .5*(wflux(l, ij) + wflux(l + 1, ij))
         W_etadot(l, ij) = wflux_ij*w_ij
         eta_dot(l, ij) = wflux_ij/mass(l, ij)
         wcov(l, ij) = w_ij*(geopot(l + 1, ij) - geopot(l, ij))
      END DO
      !$OMP END DO
      ! add NH term to du
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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))
         END DO
      END DO
      !$OMP END DO
      ! add NH terms to dW, dPhi
      ! FIXME : TODO top and bottom
      !$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))
         END DO
      END DO
      !$OMP END DO
      ! We need a barrier here because we compute W_etadot above and do a vertical difference below
      !$OMP BARRIER
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         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)
         END DO
         l = llm + 1
         dW(l, ij) = dW(l, ij) + W_etadot(l - 1, ij)
      END DO
      !$OMP END DO
      !---------------------------- caldyn_vert_NH ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_vert_NH

!----------------------------- Hydrostatic -------------------------------

   SUBROUTINE compute_geopot(rhodz, theta, ps, pk, geopot)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: rhodz, pk   ! IN, OUT
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: theta      ! IN
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: geopot    ! IN(l=1)/OUT(l>1)
      REAL(C_DOUBLE), DIMENSION(primal_num)     :: ps        ! OUT
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix, Cp_ik
      CALL enter_trace(id_geopot, 8*llm*((3)*primal_num + (0)*dual_num + (3)*edge_num))
! #include "../kernels_unst/compute_geopot.k90"
      CALL exit_trace()
   END SUBROUTINE compute_geopot
!
   SUBROUTINE compute_caldyn_slow_hydro(rhodz, theta, u, berni, hflux, du)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: rhodz, berni! IN
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: theta      ! IN
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)     :: u, hflux, du ! IN, OUT, OUT
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      LOGICAL, PARAMETER :: zero = .TRUE.
      REAL(C_DOUBLE) :: ke, uu
      CALL enter_trace(id_slow_hydro, 8*llm*((3)*primal_num + (0)*dual_num + (3)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_slow_hydro ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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)
            hflux(l, edge) = uu*le_de(edge)
         END DO
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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
         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
      IF (zero) THEN
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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
            END DO
         END DO
         !$OMP END DO
      ELSE
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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
            END DO
         END DO
         !$OMP END DO
      END IF
      !---------------------------- caldyn_slow_hydro ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_slow_hydro

!---------------------------------- Generic ------------------------------

   SUBROUTINE caldyn_vert(convm, rhodz, theta, u, dmass_col, wflux, dtheta_rhodz, du, wwuu)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(primal_num)    :: dmass_col
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: convm, rhodz
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)     :: u, du
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: dtheta_rhodz, theta
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: wflux
      REAL(C_DOUBLE), DIMENSION(llm + 1, edge_num)     :: wwuu
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: dF_deta, dFu_deta
      wwuu = 0.
      !$OMP BARRIER
      CALL enter_trace(id_vert, 8*llm*((5)*primal_num + (0)*dual_num + (3)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_wflux ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         DO l = llm - 1, 1, -1
            ! cumulate mass flux convergence from top to bottom
            convm(l, ij) = convm(l, ij) + convm(l + 1, ij)
         END DO
         dmass_col(ij) = convm(1, ij)
         DO l = 2, llm
            ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
            wflux(l, ij) = mass_bl(l, ij)*dmass_col(ij) - convm(l, ij)
         END DO
      END DO
      !$OMP END DO
      ! make sure wflux is up to date
      !$OMP BARRIER
      !---------------------------- caldyn_wflux ----------------------------------
      !--------------------------------------------------------------------------

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_dmass ----------------------------------
      !$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)
         END DO
      END DO
      !$OMP END DO
      !---------------------------- caldyn_dmass ----------------------------------
      !--------------------------------------------------------------------------

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_vert ----------------------------------
      DO iq = 1, nqdyn
         !$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)
            END DO
         END DO
         !$OMP END DO
         !$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)
            END DO
         END DO
         !$OMP END DO
      END DO
      IF (caldyn_vert_variant == caldyn_vert_cons) THEN
         ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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
            END DO
         END DO
         !$OMP END DO
         ! make sure wwuu is up to date
         !$OMP BARRIER
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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)
               dF_deta = .5*(wflux(l + 1, ij_left) + wflux(l + 1, ij_right) - (wflux(l, ij_left) + wflux(l, ij_right))) ! d/deta(F)
               du(l, edge) = du(l, edge) - (dFu_deta - u(l, edge)*dF_deta)/(.5*(rhodz(l, ij_left) + rhodz(l, ij_right))) ! (F/m)du/deta
            END DO
         END DO
         !$OMP END DO
      ELSE
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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))
            END DO
         END DO
         !$OMP END DO
         ! make sure wwuu is up to date
         !$OMP BARRIER
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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))
            END DO
         END DO
         !$OMP END DO
      END IF
      !---------------------------- caldyn_vert ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE caldyn_vert

   SUBROUTINE compute_coriolis(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)     :: hflux, Ftheta, qu, du
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: convm
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: theta, dtheta_rhodz
      INTEGER, PARAMETER :: conserv_energy = 1, conserv_enstrophy = 2, caldyn_conserv = conserv_enstrophy ! FIXME
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      REAL(C_DOUBLE) :: divF, du_trisk
      CALL enter_trace(id_coriolis, 8*llm*((3)*primal_num + (4)*dual_num + (0)*edge_num)) ! primal, dual, edge

      !--------------------------------------------------------------------------
      !---------------------------- coriolis ----------------------------------
      !
      DO iq = 1, nqdyn
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            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)
            END DO
         END DO
         !$OMP END DO
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            ! 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
      END DO ! iq
      !
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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
      !
      SELECT CASE (caldyn_conserv)
      CASE (conserv_energy) ! energy-conserving TRiSK
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            ! this VLOOP iterates over the TRISK stencil
            SELECT CASE (trisk_deg(edge))
            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 (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 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
      CASE (conserv_enstrophy) ! enstrophy-conserving TRiSK
         !$OMP DO SCHEDULE(STATIC)
         DO edge = 1, edge_num
            ! this VLOOP iterates over the TRISK stencil
            SELECT CASE (trisk_deg(edge))
            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)
                  itrisk = 2
                  edge_trisk = trisk(2, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 3
                  edge_trisk = trisk(3, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 4
                  edge_trisk = trisk(4, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  du(l, edge) = du(l, edge) + du_trisk*qu(l, edge)
               END DO
            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)
                  itrisk = 2
                  edge_trisk = trisk(2, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 3
                  edge_trisk = trisk(3, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 4
                  edge_trisk = trisk(4, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 5
                  edge_trisk = trisk(5, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 6
                  edge_trisk = trisk(6, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 7
                  edge_trisk = trisk(7, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 8
                  edge_trisk = trisk(8, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 9
                  edge_trisk = trisk(9, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  itrisk = 10
                  edge_trisk = trisk(10, edge)
                  du_trisk = du_trisk + wee(itrisk, edge)*hflux(l, edge_trisk)
                  du(l, edge) = du(l, edge) + du_trisk*qu(l, edge)
               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)
                  END DO
                  du(l, edge) = du(l, edge) + du_trisk*qu(l, edge)
               END DO
            END SELECT
         END DO
         !$OMP END DO
      END SELECT
      !---------------------------- coriolis ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE

   SUBROUTINE compute_theta(mass_col, rhodz, theta_rhodz, theta)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(primal_num) :: mass_col
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: rhodz
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: theta, theta_rhodz
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: m
      CALL enter_trace(id_theta, 8*llm*((3)*primal_num + (0)*dual_num + (0)*edge_num)) ! primal, dual, edge

      !--------------------------------------------------------------------------
      !---------------------------- theta ----------------------------------
      IF (caldyn_eta == eta_mass) THEN ! Compute mass
         ! compute mass_col from rhodz
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            mass_col(ij) = 0.
            DO l = 1, llm
               mass_col(ij) = mass_col(ij) + rhodz(l, ij)
            END DO
         END DO
         !$OMP END DO
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               rhodz(l, ij) = MASS_DAK(l, ij) + mass_col(ij)*MASS_DBK(l, ij)
            END DO
         END DO
         !$OMP END DO
      END IF
      DO iq = 1, nqdyn
         !$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)
            END DO
         END DO
         !$OMP END DO
      END DO
      !---------------------------- theta ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE

   SUBROUTINE compute_pvort_only(rhodz, u, qv, qu)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: rhodz
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: u, qu
      REAL(C_DOUBLE), DIMENSION(llm, dual_num)    :: qv
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      INTEGER :: vertex1, vertex2, vertex3, vertex4, vertex5, vertex6
      REAL(C_DOUBLE) :: etav, hv
      CALL enter_trace(id_pvort_only, 8*llm*((1)*primal_num + (1)*dual_num + (2)*edge_num)) ! primal, dual, edge

      !--------------------------------------------------------------------------
      !---------------------------- pvort_only ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, dual_num
         ! 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 + R_iv2(1, ij)*rhodz(l, vertex1)
               hv = hv + R_iv2(2, ij)*rhodz(l, vertex2)
               hv = hv + R_iv2(3, ij)*rhodz(l, vertex3)
               qv(l, ij) = (etav + fv(ij)*Av(ij))/(hv*Av(ij))
            END DO
         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 + R_iv2(1, ij)*rhodz(l, vertex1)
               hv = hv + R_iv2(2, ij)*rhodz(l, vertex2)
               hv = hv + R_iv2(3, ij)*rhodz(l, vertex3)
               hv = hv + R_iv2(4, ij)*rhodz(l, vertex4)
               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 + R_iv2(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
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- pvort_only ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_pvort_only

   SUBROUTINE compute_caldyn_fast(tau, pk, berni, theta, geopot, du, u)
      USE grid_param
      REAL(C_DOUBLE), INTENT(IN) :: tau
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: pk, berni  ! INOUT, OUT
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn)  :: theta     ! IN
      REAL(C_DOUBLE), DIMENSION(llm + 1, primal_num) :: geopot    ! IN
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: u, du      ! INOUT,INOUT
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
         de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      REAL(C_DOUBLE) :: due, cp_ik, Phi_ik
      CALL enter_trace(id_fast, 8*llm*((4)*primal_num + (0)*dual_num + (2)*edge_num)) ! primal, dual, edge

      !--------------------------------------------------------------------------
      !---------------------------- caldyn_fast ----------------------------------
      !
      SELECT CASE (caldyn_thermo)
      CASE (thermo_boussinesq)
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               berni(l, ij) = pk(l, ij)
               ! from now on pk contains the vertically-averaged geopotential
               pk(l, ij) = .5*(geopot(l, ij) + geopot(l + 1, ij))
            END DO
         END DO
         !$OMP END DO
      CASE (thermo_theta)
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               Phi_ik = .5*(geopot(l, ij) + geopot(l + 1, ij))
               berni(l, ij) = Phi_ik
            END DO
         END DO
         !$OMP END DO
      CASE (thermo_entropy)
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               Phi_ik = .5*(geopot(l, ij) + geopot(l + 1, ij))
               berni(l, ij) = Phi_ik + pk(l, ij)*(cpp - theta(l, ij, 1)) ! Gibbs = Cp.T-Ts = T(Cp-s)
            END DO
         END DO
         !$OMP END DO
      CASE (thermo_variable_Cp)
         !$OMP DO SCHEDULE(STATIC)
         DO ij = 1, primal_num
            !DIR$ SIMD
            DO l = 1, llm
               Phi_ik = .5*(geopot(l, ij) + geopot(l + 1, ij))
               ! thermodynamics with variable Cp
               ! Cp(T) = Cp0 * (T/T0)^nu
               ! => h = Cp(T).T/(nu+1)
               cp_ik = cpp*(pk(ij, l)/Treff)**nu
               berni(l, ij) = Phi_ik + pk(l, ij)*(cp_ik/(nu + 1.) - theta(l, ij, 1)) ! Gibbs = h-Ts = T(Cp/(nu+1)-s)
            END DO
         END DO
         !$OMP END DO
      CASE DEFAULT
         PRINT *, 'Unsupported value of caldyn_thermo : ', caldyn_thermo ! FIXME
         STOP
      END SELECT
      !
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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)
            du(l, edge) = du(l, edge) - 1.*due
            u(l, edge) = u(l, edge) + tau*du(l, edge)
         END DO
      END DO
      !$OMP END DO
      !
      !---------------------------- caldyn_fast ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_caldyn_fast

!----------------------------- Unused -----------------------------

   SUBROUTINE gradient(b, grad) BIND(C, name='dynamico_gradient')
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: b
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: grad
      INTEGER ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown

      !--------------------------------------------------------------------------
      !---------------------------- gradient ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- gradient ----------------------------------
      !--------------------------------------------------------------------------
   END SUBROUTINE
!
   SUBROUTINE div(u, divu) BIND(C, name='dynamico_div')
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num) :: divu
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)    :: u
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      REAL(C_DOUBLE) :: div_ij
      !$OMP PARALLEL NUM_THREADS(nb_threads)

      !--------------------------------------------------------------------------
      !---------------------------- div ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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
         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
      !---------------------------- div ----------------------------------
      !--------------------------------------------------------------------------
      !$OMP END PARALLEL
   END SUBROUTINE

   SUBROUTINE compute_scalar_laplacian(b, grad, divu)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)   :: b, divu ! IN, OUT
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)      :: grad ! BUF
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      REAL(C_DOUBLE) :: div_ij
      CALL enter_trace(id_scalar_laplacian, 8*llm*((0)*primal_num + (0)*dual_num + (1)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- scalar_laplacian ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         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))
         END DO
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         ! 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*grad(l, edge1)
               div_ij = div_ij + sign2*le_de2*grad(l, edge2)
               div_ij = div_ij + sign3*le_de3*grad(l, edge3)
               div_ij = div_ij + sign4*le_de4*grad(l, edge4)
               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*grad(l, edge1)
               div_ij = div_ij + sign2*le_de2*grad(l, edge2)
               div_ij = div_ij + sign3*le_de3*grad(l, edge3)
               div_ij = div_ij + sign4*le_de4*grad(l, edge4)
               div_ij = div_ij + sign5*le_de5*grad(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*grad(l, edge1)
               div_ij = div_ij + sign2*le_de2*grad(l, edge2)
               div_ij = div_ij + sign3*le_de3*grad(l, edge3)
               div_ij = div_ij + sign4*le_de4*grad(l, edge4)
               div_ij = div_ij + sign5*le_de5*grad(l, edge5)
               div_ij = div_ij + sign6*le_de6*grad(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)*grad(l, edge)
               END DO
               divu(l, ij) = div_ij/Ai(ij)
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !---------------------------- scalar_laplacian ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_scalar_laplacian

   SUBROUTINE compute_curl_laplacian(u, qv, curlcurl)
      USE grid_param
      REAL(C_DOUBLE), DIMENSION(llm, dual_num)   :: qv ! BUF
      REAL(C_DOUBLE), DIMENSION(llm, edge_num)   :: u, curlcurl ! IN,OUT
      INTEGER :: ij, l, iq, iedge, edge, ivertex, vertex, ij_left, ij_right, ij_up, ij_down, itrisk, edge_trisk, kup, kdown
      REAL(C_DOUBLE) :: sign1, sign2, sign3, sign4, sign5, sign6, le_de1, le_de2, le_de3, le_de4, le_de5, le_de6, &
                        de1, de2, de3, de4, de5, de6, le1, le2, le3, le4, le5, le6
      INTEGER :: edge1, edge2, edge3, edge4, edge5, edge6, ij_up1, ij_up2, ij_up3, ij_up4, ij_up5, ij_up6, &
                 ij_down1, ij_down2, ij_down3, ij_down4, ij_down5, ij_down6
      INTEGER :: vertex1, vertex2, vertex3, vertex4, vertex5, vertex6
      REAL(C_DOUBLE) :: etav
      CALL enter_trace(id_scalar_laplacian, 8*llm*((0)*primal_num + (0)*dual_num + (1)*edge_num))

      !--------------------------------------------------------------------------
      !---------------------------- curl_laplacian ----------------------------------
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, dual_num
         ! 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)
            !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)
               qv(l, ij) = etav/Av(ij)
            END DO
         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)
            !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)
               qv(l, ij) = etav/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
               qv(l, ij) = etav/Av(ij)
            END DO
         END SELECT
      END DO
      !$OMP END DO
      !$OMP DO SCHEDULE(STATIC)
      DO edge = 1, edge_num
         ij_up = up(edge)
         ij_down = down(edge)
         !DIR$ SIMD
         DO l = 1, llm
            curlcurl(l, edge) = 1.*(qv(l, ij_down) - qv(l, ij_up))*(1./le_de(edge))
         END DO
      END DO
      !$OMP END DO
      !---------------------------- curl_laplacian ----------------------------------
      !--------------------------------------------------------------------------
      CALL exit_trace()
   END SUBROUTINE compute_curl_laplacian

END MODULE caldyn_unstructured_mod
