SUBMODULE(compute_caldyn_mod) compute_caldyn_Coriolis_mod
   USE prec, ONLY: rstd
   USE caldyn_vars_mod, ONLY: caldyn_conserv, conserv_energy, conserv_enstrophy, conserv_gassmann
   USE grid_param
   USE earth_const
   USE disvert_mod
   USE omp_para
   USE trace
   IMPLICIT NONE

CONTAINS

   MODULE SUBROUTINE compute_caldyn_coriolis(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du)
      USE grid_param, ONLY: grid_type, grid_ico
      USE geometry, ONLY: Ai, wee
      REAL(rstd), INTENT(IN)    :: hflux(:, :), theta(:, :, :), qu(:, :)
      REAL(rstd), INTENT(OUT)   :: Ftheta(:, :), convm(:, :), dtheta_rhodz(:, :, :)
      REAL(rstd), INTENT(INOUT) :: du(:, :)
      IF (grid_type == grid_ico) THEN
         CALL compute_caldyn_coriolis_hex_(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du, Ai, wee)
      ELSE
         CALL compute_caldyn_coriolis_unst_(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du)
      END IF
   END SUBROUTINE compute_caldyn_coriolis

!--------------------------------------------------------------

   SUBROUTINE compute_caldyn_coriolis_unst_(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du)
      USE ISO_C_BINDING, only: C_DOUBLE, C_FLOAT
      USE data_unstructured_mod, ONLY: enter_trace, exit_trace, &
                                       id_coriolis, left, right, primal_deg, primal_edge, primal_ne, &
                                       trisk_deg, trisk, wee

      REAL(C_DOUBLE), DIMENSION(llm, edge_num)     :: hflux, Ftheta, qu, du ! IN, BUF, IN, INOUT
      REAL(C_DOUBLE), DIMENSION(llm, primal_num)  :: convm                 ! BUF
      REAL(C_DOUBLE), DIMENSION(llm, primal_num, nqdyn) :: theta, dtheta_rhodz   ! 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
      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 compute_caldyn_coriolis_unst_

   SUBROUTINE compute_caldyn_Coriolis_hex_(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du, Ai, wee)
      USE icosa, hide_Ai => Ai, hide_wee => wee
      REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm, llm)  ! hflux in kg/s
      REAL(rstd), INTENT(IN)    :: theta(iim*jjm, llm, nqdyn) ! active scalars
      REAL(rstd), INTENT(IN)    :: qu(3*iim*jjm, llm)
      REAL(rstd), INTENT(OUT)   :: Ftheta(3*iim*jjm, llm)  ! potential temperature flux
      REAL(rstd), INTENT(OUT)   :: convm(iim*jjm, llm)  ! mass flux convergence
      REAL(rstd), INTENT(OUT)   :: dtheta_rhodz(iim*jjm, llm, nqdyn)
      REAL(rstd), INTENT(INOUT) :: du(3*iim*jjm, llm)
      REAL(rstd), INTENT(IN)    :: Ai(iim*jjm)
      REAL(rstd), INTENT(IN)    :: wee(3*iim*jjm, 10)

      REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
      INTEGER :: ij, iq, l, kdown
      CALL trace_start("compute_caldyn_Coriolis")

      !$acc data async &
      !$acc present(hflux(:,:), theta(:,:,:), qu(:,:), Ftheta(:,:), convm(:,:), dtheta_rhodz(:,:,:), du(:,:), Ai(:), wee(:,:))

      DO iq = 1, nqdyn
         !$acc parallel loop collapse(2) present(Ftheta, theta) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               Ftheta(ij + u_right, l) = .5*(theta(ij, l, iq) + theta(ij + t_right, l, iq))*hflux(ij + u_right, l)
               Ftheta(ij + u_lup, l) = .5*(theta(ij, l, iq) + theta(ij + t_lup, l, iq))*hflux(ij + u_lup, l)
               Ftheta(ij + u_ldown, l) = .5*(theta(ij, l, iq) + theta(ij + t_ldown, l, iq))*hflux(ij + u_ldown, l)
            END DO
         END DO
         !$acc parallel loop collapse(2) present(Ftheta) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               divF = 0.
               divF = divF + Ftheta(ij + u_rup, l)*ne_rup
               divF = divF + Ftheta(ij + u_lup, l)*ne_lup
               divF = divF + Ftheta(ij + u_left, l)*ne_left
               divF = divF + Ftheta(ij + u_ldown, l)*ne_ldown
               divF = divF + Ftheta(ij + u_rdown, l)*ne_rdown
               divF = divF + Ftheta(ij + u_right, l)*ne_right
               dtheta_rhodz(ij, l, iq) = -divF/Ai(ij)
            END DO
         END DO
      END DO ! iq

      !$acc parallel loop collapse(2) present(Ai) async
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            divF = 0.
            divF = divF + hflux(ij + u_rup, l)*ne_rup
            divF = divF + hflux(ij + u_lup, l)*ne_lup
            divF = divF + hflux(ij + u_left, l)*ne_left
            divF = divF + hflux(ij + u_ldown, l)*ne_ldown
            divF = divF + hflux(ij + u_rdown, l)*ne_rdown
            divF = divF + hflux(ij + u_right, l)*ne_right
            convm(ij, l) = -divF/Ai(ij)
         END DO
      END DO
      !
      SELECT CASE (caldyn_conserv)
      CASE (conserv_energy) ! energy-conserving TRiSK
         !$acc parallel loop collapse(2) present(qu, wee, du) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_right, 1)*hflux(ij + u_rup, l)* &
                          (qu(ij + u_right, l) + qu(ij + u_rup, l))
               du_trisk = du_trisk + wee(ij + u_right, 2)*hflux(ij + u_lup, l)* &
                          (qu(ij + u_right, l) + qu(ij + u_lup, l))
               du_trisk = du_trisk + wee(ij + u_right, 3)*hflux(ij + u_left, l)* &
                          (qu(ij + u_right, l) + qu(ij + u_left, l))
               du_trisk = du_trisk + wee(ij + u_right, 4)*hflux(ij + u_ldown, l)* &
                          (qu(ij + u_right, l) + qu(ij + u_ldown, l))
               du_trisk = du_trisk + wee(ij + u_right, 5)*hflux(ij + u_rdown, l)* &
                          (qu(ij + u_right, l) + qu(ij + u_rdown, l))
               du_trisk = du_trisk + wee(ij + u_right, 6)*hflux(ij + t_right + u_ldown, l)* &
                          (qu(ij + u_right, l) + qu(ij + t_right + u_ldown, l))
               du_trisk = du_trisk + wee(ij + u_right, 7)*hflux(ij + t_right + u_rdown, l)* &
                          (qu(ij + u_right, l) + qu(ij + t_right + u_rdown, l))
               du_trisk = du_trisk + wee(ij + u_right, 8)*hflux(ij + t_right + u_right, l)* &
                          (qu(ij + u_right, l) + qu(ij + t_right + u_right, l))
               du_trisk = du_trisk + wee(ij + u_right, 9)*hflux(ij + t_right + u_rup, l)* &
                          (qu(ij + u_right, l) + qu(ij + t_right + u_rup, l))
               du_trisk = du_trisk + wee(ij + u_right, 10)*hflux(ij + t_right + u_lup, l)* &
                          (qu(ij + u_right, l) + qu(ij + t_right + u_lup, l))
               du(ij + u_right, l) = du(ij + u_right, l) + .5*du_trisk
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_lup, 1)*hflux(ij + u_left, l)* &
                          (qu(ij + u_lup, l) + qu(ij + u_left, l))
               du_trisk = du_trisk + wee(ij + u_lup, 2)*hflux(ij + u_ldown, l)* &
                          (qu(ij + u_lup, l) + qu(ij + u_ldown, l))
               du_trisk = du_trisk + wee(ij + u_lup, 3)*hflux(ij + u_rdown, l)* &
                          (qu(ij + u_lup, l) + qu(ij + u_rdown, l))
               du_trisk = du_trisk + wee(ij + u_lup, 4)*hflux(ij + u_right, l)* &
                          (qu(ij + u_lup, l) + qu(ij + u_right, l))
               du_trisk = du_trisk + wee(ij + u_lup, 5)*hflux(ij + u_rup, l)* &
                          (qu(ij + u_lup, l) + qu(ij + u_rup, l))
               du_trisk = du_trisk + wee(ij + u_lup, 6)*hflux(ij + t_lup + u_right, l)* &
                          (qu(ij + u_lup, l) + qu(ij + t_lup + u_right, l))
               du_trisk = du_trisk + wee(ij + u_lup, 7)*hflux(ij + t_lup + u_rup, l)* &
                          (qu(ij + u_lup, l) + qu(ij + t_lup + u_rup, l))
               du_trisk = du_trisk + wee(ij + u_lup, 8)*hflux(ij + t_lup + u_lup, l)* &
                          (qu(ij + u_lup, l) + qu(ij + t_lup + u_lup, l))
               du_trisk = du_trisk + wee(ij + u_lup, 9)*hflux(ij + t_lup + u_left, l)* &
                          (qu(ij + u_lup, l) + qu(ij + t_lup + u_left, l))
               du_trisk = du_trisk + wee(ij + u_lup, 10)*hflux(ij + t_lup + u_ldown, l)* &
                          (qu(ij + u_lup, l) + qu(ij + t_lup + u_ldown, l))
               du(ij + u_lup, l) = du(ij + u_lup, l) + .5*du_trisk
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_ldown, 1)*hflux(ij + u_rdown, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + u_rdown, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 2)*hflux(ij + u_right, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + u_right, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 3)*hflux(ij + u_rup, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + u_rup, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 4)*hflux(ij + u_lup, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + u_lup, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 5)*hflux(ij + u_left, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + u_left, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 6)*hflux(ij + t_ldown + u_lup, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_lup, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 7)*hflux(ij + t_ldown + u_left, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_left, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 8)*hflux(ij + t_ldown + u_ldown, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_ldown, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 9)*hflux(ij + t_ldown + u_rdown, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_rdown, l))
               du_trisk = du_trisk + wee(ij + u_ldown, 10)*hflux(ij + t_ldown + u_right, l)* &
                          (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_right, l))
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + .5*du_trisk
            END DO
         END DO
      CASE (conserv_enstrophy) ! enstrophy-conserving TRiSK
         !$acc parallel loop collapse(2) present(wee, du) async
         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_right, 1)*hflux(ij + u_rup, l)
               du_trisk = du_trisk + wee(ij + u_right, 2)*hflux(ij + u_lup, l)
               du_trisk = du_trisk + wee(ij + u_right, 3)*hflux(ij + u_left, l)
               du_trisk = du_trisk + wee(ij + u_right, 4)*hflux(ij + u_ldown, l)
               du_trisk = du_trisk + wee(ij + u_right, 5)*hflux(ij + u_rdown, l)
               du_trisk = du_trisk + wee(ij + u_right, 6)*hflux(ij + t_right + u_ldown, l)
               du_trisk = du_trisk + wee(ij + u_right, 7)*hflux(ij + t_right + u_rdown, l)
               du_trisk = du_trisk + wee(ij + u_right, 8)*hflux(ij + t_right + u_right, l)
               du_trisk = du_trisk + wee(ij + u_right, 9)*hflux(ij + t_right + u_rup, l)
               du_trisk = du_trisk + wee(ij + u_right, 10)*hflux(ij + t_right + u_lup, l)
               du(ij + u_right, l) = du(ij + u_right, l) + du_trisk*qu(ij + u_right, l)
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_lup, 1)*hflux(ij + u_left, l)
               du_trisk = du_trisk + wee(ij + u_lup, 2)*hflux(ij + u_ldown, l)
               du_trisk = du_trisk + wee(ij + u_lup, 3)*hflux(ij + u_rdown, l)
               du_trisk = du_trisk + wee(ij + u_lup, 4)*hflux(ij + u_right, l)
               du_trisk = du_trisk + wee(ij + u_lup, 5)*hflux(ij + u_rup, l)
               du_trisk = du_trisk + wee(ij + u_lup, 6)*hflux(ij + t_lup + u_right, l)
               du_trisk = du_trisk + wee(ij + u_lup, 7)*hflux(ij + t_lup + u_rup, l)
               du_trisk = du_trisk + wee(ij + u_lup, 8)*hflux(ij + t_lup + u_lup, l)
               du_trisk = du_trisk + wee(ij + u_lup, 9)*hflux(ij + t_lup + u_left, l)
               du_trisk = du_trisk + wee(ij + u_lup, 10)*hflux(ij + t_lup + u_ldown, l)
               du(ij + u_lup, l) = du(ij + u_lup, l) + du_trisk*qu(ij + u_lup, l)
               du_trisk = 0.
               du_trisk = du_trisk + wee(ij + u_ldown, 1)*hflux(ij + u_rdown, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 2)*hflux(ij + u_right, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 3)*hflux(ij + u_rup, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 4)*hflux(ij + u_lup, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 5)*hflux(ij + u_left, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 6)*hflux(ij + t_ldown + u_lup, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 7)*hflux(ij + t_ldown + u_left, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 8)*hflux(ij + t_ldown + u_ldown, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 9)*hflux(ij + t_ldown + u_rdown, l)
               du_trisk = du_trisk + wee(ij + u_ldown, 10)*hflux(ij + t_ldown + u_right, l)
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + du_trisk*qu(ij + u_ldown, l)
            END DO
         END DO
      END SELECT

      !$acc end data
      CALL trace_end("compute_caldyn_Coriolis")
   END SUBROUTINE compute_caldyn_Coriolis_hex_

   SUBROUTINE compute_caldyn_Coriolis_manual(hflux, theta, qu, Ftheta, convm, dtheta_rhodz, du)
      USE icosa
      USE caldyn_vars_mod
      REAL(rstd), INTENT(IN)    :: hflux(3*iim*jjm, llm)  ! hflux in kg/s
      REAL(rstd), INTENT(IN)    :: theta(iim*jjm, llm, nqdyn) ! active scalars
      REAL(rstd), INTENT(IN)    :: qu(3*iim*jjm, llm)
      REAL(rstd), INTENT(IN)    :: Ftheta(3*iim*jjm, llm)  ! ignored in favor of local buffer
      REAL(rstd), INTENT(OUT)   :: convm(iim*jjm, llm)  ! mass flux convergence
      REAL(rstd), INTENT(OUT)   :: dtheta_rhodz(iim*jjm, llm, nqdyn)
      REAL(rstd), INTENT(INOUT) :: du(3*iim*jjm, llm)
      REAL(rstd)  :: buf_Ftheta(3*iim*jjm)  ! local buffer for potential temperature flux
      REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
      INTEGER :: ij, iq, l, kdown

      CALL trace_start("compute_caldyn_Coriolis")

      DO l = ll_begin, ll_end
         ! compute theta flux
         DO iq = 1, nqdyn
            !DIR$ SIMD
            DO ij = ij_begin_ext, ij_end_ext
               buf_Ftheta(ij + u_right) = 0.5*(theta(ij, l, iq) + theta(ij + t_right, l, iq)) &
                                          *hflux(ij + u_right, l)
               buf_Ftheta(ij + u_lup) = 0.5*(theta(ij, l, iq) + theta(ij + t_lup, l, iq)) &
                                        *hflux(ij + u_lup, l)
               buf_Ftheta(ij + u_ldown) = 0.5*(theta(ij, l, iq) + theta(ij + t_ldown, l, iq)) &
                                          *hflux(ij + u_ldown, l)
            END DO
            ! horizontal divergence of fluxes
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               ! dtheta_rhodz =  -div(flux.theta)
               dtheta_rhodz(ij, l, iq) = &
                  -1./Ai(ij)*(ne_right*buf_Ftheta(ij + u_right) + &
                              ne_rup*buf_Ftheta(ij + u_rup) + &
                              ne_lup*buf_Ftheta(ij + u_lup) + &
                              ne_left*buf_Ftheta(ij + u_left) + &
                              ne_ldown*buf_Ftheta(ij + u_ldown) + &
                              ne_rdown*buf_Ftheta(ij + u_rdown))
            END DO
         END DO

         !DIR$ SIMD
         DO ij = ij_begin, ij_end
            ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
            convm(ij, l) = -1./Ai(ij)*(ne_right*hflux(ij + u_right, l) + &
                                       ne_rup*hflux(ij + u_rup, l) + &
                                       ne_lup*hflux(ij + u_lup, l) + &
                                       ne_left*hflux(ij + u_left, l) + &
                                       ne_ldown*hflux(ij + u_ldown, l) + &
                                       ne_rdown*hflux(ij + u_rdown, l))
         END DO ! ij
      END DO ! llm

!!! Compute potential vorticity (Coriolis) contribution to du
      SELECT CASE (caldyn_conserv)

      CASE (conserv_energy) ! energy-conserving TRiSK

         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               uu_right = &
                  wee(ij + u_right, 1)*hflux(ij + u_rup, l)*(qu(ij + u_right, l) + qu(ij + u_rup, l)) + &
                  wee(ij + u_right, 2)*hflux(ij + u_lup, l)*(qu(ij + u_right, l) + qu(ij + u_lup, l)) + &
                  wee(ij + u_right, 3)*hflux(ij + u_left, l)*(qu(ij + u_right, l) + qu(ij + u_left, l)) + &
                  wee(ij + u_right, 4)*hflux(ij + u_ldown, l)*(qu(ij + u_right, l) + qu(ij + u_ldown, l)) + &
                  wee(ij + u_right, 5)*hflux(ij + u_rdown, l)*(qu(ij + u_right, l) + qu(ij + u_rdown, l)) + &
                  wee(ij + u_right, 6)*hflux(ij + t_right + u_ldown, l)*(qu(ij + u_right, l) + qu(ij + t_right + u_ldown, l)) + &
                  wee(ij + u_right, 7)*hflux(ij + t_right + u_rdown, l)*(qu(ij + u_right, l) + qu(ij + t_right + u_rdown, l)) + &
                  wee(ij + u_right, 8)*hflux(ij + t_right + u_right, l)*(qu(ij + u_right, l) + qu(ij + t_right + u_right, l)) + &
                  wee(ij + u_right, 9)*hflux(ij + t_right + u_rup, l)*(qu(ij + u_right, l) + qu(ij + t_right + u_rup, l)) + &
                  wee(ij + u_right, 10)*hflux(ij + t_right + u_lup, l)*(qu(ij + u_right, l) + qu(ij + t_right + u_lup, l))
               uu_lup = &
                  wee(ij + u_lup, 1)*hflux(ij + u_left, l)*(qu(ij + u_lup, l) + qu(ij + u_left, l)) + &
                  wee(ij + u_lup, 2)*hflux(ij + u_ldown, l)*(qu(ij + u_lup, l) + qu(ij + u_ldown, l)) + &
                  wee(ij + u_lup, 3)*hflux(ij + u_rdown, l)*(qu(ij + u_lup, l) + qu(ij + u_rdown, l)) + &
                  wee(ij + u_lup, 4)*hflux(ij + u_right, l)*(qu(ij + u_lup, l) + qu(ij + u_right, l)) + &
                  wee(ij + u_lup, 5)*hflux(ij + u_rup, l)*(qu(ij + u_lup, l) + qu(ij + u_rup, l)) + &
                  wee(ij + u_lup, 6)*hflux(ij + t_lup + u_right, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_right, l)) + &
                  wee(ij + u_lup, 7)*hflux(ij + t_lup + u_rup, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_rup, l)) + &
                  wee(ij + u_lup, 8)*hflux(ij + t_lup + u_lup, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_lup, l)) + &
                  wee(ij + u_lup, 9)*hflux(ij + t_lup + u_left, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_left, l)) + &
                  wee(ij + u_lup, 10)*hflux(ij + t_lup + u_ldown, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_ldown, l))
               uu_ldown = &
                  wee(ij + u_ldown, 1)*hflux(ij + u_rdown, l)*(qu(ij + u_ldown, l) + qu(ij + u_rdown, l)) + &
                  wee(ij + u_ldown, 2)*hflux(ij + u_right, l)*(qu(ij + u_ldown, l) + qu(ij + u_right, l)) + &
                  wee(ij + u_ldown, 3)*hflux(ij + u_rup, l)*(qu(ij + u_ldown, l) + qu(ij + u_rup, l)) + &
                  wee(ij + u_ldown, 4)*hflux(ij + u_lup, l)*(qu(ij + u_ldown, l) + qu(ij + u_lup, l)) + &
                  wee(ij + u_ldown, 5)*hflux(ij + u_left, l)*(qu(ij + u_ldown, l) + qu(ij + u_left, l)) + &
                  wee(ij + u_ldown, 6)*hflux(ij + t_ldown + u_lup, l)*(qu(ij + u_ldown, l) + qu(ij + t_ldown + u_lup, l)) + &
                  wee(ij + u_ldown, 7)*hflux(ij + t_ldown + u_left, l)*(qu(ij + u_ldown, l) + qu(ij + t_ldown + u_left, l)) + &
                  wee(ij + u_ldown, 8)*hflux(ij + t_ldown + u_ldown, l)*(qu(ij + u_ldown, l) + qu(ij + t_ldown + u_ldown, l)) + &
                  wee(ij + u_ldown, 9)*hflux(ij + t_ldown + u_rdown, l)*(qu(ij + u_ldown, l) + qu(ij + t_ldown + u_rdown, l)) + &
                  wee(ij + u_ldown, 10)*hflux(ij + t_ldown + u_right, l)*(qu(ij + u_ldown, l) + qu(ij + t_ldown + u_right, l))
               du(ij + u_right, l) = du(ij + u_right, l) + .5*uu_right
               du(ij + u_lup, l) = du(ij + u_lup, l) + .5*uu_lup
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + .5*uu_ldown
            END DO
         END DO

      CASE (conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)

         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               uu_right = &
                  wee(ij + u_right, 1)*hflux(ij + u_rup, l)*qu(ij + t_right + u_lup, l) + &
                  wee(ij + u_right, 2)*hflux(ij + u_lup, l)*qu(ij + u_rup, l) + &
                  .5*wee(ij + u_right, 3)*hflux(ij + u_left, l)*(qu(ij + u_right, l) + qu(ij + u_left, l)) + &
                  wee(ij + u_right, 4)*hflux(ij + u_ldown, l)*qu(ij + u_rdown, l) + &
                  wee(ij + u_right, 5)*hflux(ij + u_rdown, l)*qu(ij + t_right + u_ldown, l) + &
                  wee(ij + u_right, 6)*hflux(ij + t_right + u_ldown, l)*qu(ij + u_rdown, l) + &
                  wee(ij + u_right, 7)*hflux(ij + t_right + u_rdown, l)*qu(ij + t_right + u_ldown, l) + &
                  .5*wee(ij + u_right, 8)*hflux(ij + t_right + u_right, l)* &
                  (qu(ij + u_right, l) + qu(ij + t_right + u_right, l)) + &
                  wee(ij + u_right, 9)*hflux(ij + t_right + u_rup, l)*qu(ij + t_right + u_lup, l) + &
                  wee(ij + u_right, 10)*hflux(ij + t_right + u_lup, l)*qu(ij + u_rup, l)
               uu_lup = &
                  wee(ij + u_lup, 1)*hflux(ij + u_left, l)*qu(ij + t_lup + u_ldown, l) + &
                  wee(ij + u_lup, 2)*hflux(ij + u_ldown, l)*qu(ij + u_left, l) + &
                  .5*wee(ij + u_lup, 3)*hflux(ij + u_rdown, l)*(qu(ij + u_lup, l) + qu(ij + u_rdown, l)) + &
                  wee(ij + u_lup, 4)*hflux(ij + u_right, l)*qu(ij + u_rup, l) + &
                  wee(ij + u_lup, 5)*hflux(ij + u_rup, l)*qu(ij + t_lup + u_right, l) + &
                  wee(ij + u_lup, 6)*hflux(ij + t_lup + u_right, l)*qu(ij + u_rup, l) + &
                  wee(ij + u_lup, 7)*hflux(ij + t_lup + u_rup, l)*qu(ij + t_lup + u_right, l) + &
                  .5*wee(ij + u_lup, 8)*hflux(ij + t_lup + u_lup, l)*(qu(ij + u_lup, l) + qu(ij + t_lup + u_lup, l)) + &
                  wee(ij + u_lup, 9)*hflux(ij + t_lup + u_left, l)*qu(ij + t_lup + u_ldown, l) + &
                  wee(ij + u_lup, 10)*hflux(ij + t_lup + u_ldown, l)*qu(ij + u_left, l)
               uu_ldown = &
                  wee(ij + u_ldown, 1)*hflux(ij + u_rdown, l)*qu(ij + t_ldown, l + u_right) + &
                  wee(ij + u_ldown, 2)*hflux(ij + u_right, l)*qu(ij + u_rdown, l) + &
                  .5*wee(ij + u_ldown, 3)*hflux(ij + u_rup, l)*(qu(ij + u_ldown, l) + qu(ij + u_rup, l)) + &
                  wee(ij + u_ldown, 4)*hflux(ij + u_lup, l)*qu(ij + u_left, l) + &
                  wee(ij + u_ldown, 5)*hflux(ij + u_left, l)*qu(ij + t_ldown + u_lup, l) + &
                  wee(ij + u_ldown, 6)*hflux(ij + t_ldown + u_lup, l)*qu(ij + u_left, l) + &
                  wee(ij + u_ldown, 7)*hflux(ij + t_ldown + u_left, l)*qu(ij + t_ldown + u_lup, l) + &
                  .5*wee(ij + u_ldown, 8)*hflux(ij + t_ldown + u_ldown, l)* &
                  (qu(ij + u_ldown, l) + qu(ij + t_ldown + u_ldown, l)) + &
                  wee(ij + u_ldown, 9)*hflux(ij + t_ldown + u_rdown, l)*qu(ij + t_ldown + u_right, l) + &
                  wee(ij + u_ldown, 10)*hflux(ij + t_ldown + u_right, l)*qu(ij + u_rdown, l)
               du(ij + u_right, l) = du(ij + u_right, l) + uu_right
               du(ij + u_lup, l) = du(ij + u_lup, l) + uu_lup
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + uu_ldown
            END DO
         END DO

      CASE (conserv_enstrophy) ! enstrophy-conserving TRiSK

         DO l = ll_begin, ll_end
            !DIR$ SIMD
            DO ij = ij_begin, ij_end
               uu_right = &
                  wee(ij + u_right, 1)*hflux(ij + u_rup, l) + &
                  wee(ij + u_right, 2)*hflux(ij + u_lup, l) + &
                  wee(ij + u_right, 3)*hflux(ij + u_left, l) + &
                  wee(ij + u_right, 4)*hflux(ij + u_ldown, l) + &
                  wee(ij + u_right, 5)*hflux(ij + u_rdown, l) + &
                  wee(ij + u_right, 6)*hflux(ij + t_right + u_ldown, l) + &
                  wee(ij + u_right, 7)*hflux(ij + t_right + u_rdown, l) + &
                  wee(ij + u_right, 8)*hflux(ij + t_right + u_right, l) + &
                  wee(ij + u_right, 9)*hflux(ij + t_right + u_rup, l) + &
                  wee(ij + u_right, 10)*hflux(ij + t_right + u_lup, l)
               uu_lup = &
                  wee(ij + u_lup, 1)*hflux(ij + u_left, l) + &
                  wee(ij + u_lup, 2)*hflux(ij + u_ldown, l) + &
                  wee(ij + u_lup, 3)*hflux(ij + u_rdown, l) + &
                  wee(ij + u_lup, 4)*hflux(ij + u_right, l) + &
                  wee(ij + u_lup, 5)*hflux(ij + u_rup, l) + &
                  wee(ij + u_lup, 6)*hflux(ij + t_lup + u_right, l) + &
                  wee(ij + u_lup, 7)*hflux(ij + t_lup + u_rup, l) + &
                  wee(ij + u_lup, 8)*hflux(ij + t_lup + u_lup, l) + &
                  wee(ij + u_lup, 9)*hflux(ij + t_lup + u_left, l) + &
                  wee(ij + u_lup, 10)*hflux(ij + t_lup + u_ldown, l)
               uu_ldown = &
                  wee(ij + u_ldown, 1)*hflux(ij + u_rdown, l) + &
                  wee(ij + u_ldown, 2)*hflux(ij + u_right, l) + &
                  wee(ij + u_ldown, 3)*hflux(ij + u_rup, l) + &
                  wee(ij + u_ldown, 4)*hflux(ij + u_lup, l) + &
                  wee(ij + u_ldown, 5)*hflux(ij + u_left, l) + &
                  wee(ij + u_ldown, 6)*hflux(ij + t_ldown + u_lup, l) + &
                  wee(ij + u_ldown, 7)*hflux(ij + t_ldown + u_left, l) + &
                  wee(ij + u_ldown, 8)*hflux(ij + t_ldown + u_ldown, l) + &
                  wee(ij + u_ldown, 9)*hflux(ij + t_ldown + u_rdown, l) + &
                  wee(ij + u_ldown, 10)*hflux(ij + t_ldown + u_right, l)

               du(ij + u_right, l) = du(ij + u_right, l) + uu_right*qu(ij + u_right, l)
               du(ij + u_lup, l) = du(ij + u_lup, l) + uu_lup*qu(ij + u_lup, l)
               du(ij + u_ldown, l) = du(ij + u_ldown, l) + uu_ldown*qu(ij + u_ldown, l)
            END DO
         END DO
      CASE DEFAULT
         STOP
      END SELECT

      CALL trace_end("compute_caldyn_Coriolis")

   END SUBROUTINE compute_caldyn_Coriolis_manual

END SUBMODULE compute_caldyn_Coriolis_mod
