Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 843)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 844)
@@ -7,4 +7,5 @@
   USE compute_theta_mod, ONLY : compute_theta
   USE compute_caldyn_kv_mod, ONLY : compute_caldyn_kv
+  USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis
   IMPLICIT NONE
   PRIVATE
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 843)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 844)
@@ -13,6 +13,5 @@
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_caldyn_Coriolis, &
-       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
+  PUBLIC :: compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
        compute_caldyn_solver, compute_caldyn_fast
 
@@ -434,209 +433,4 @@
 
   END SUBROUTINE compute_caldyn_fast
-
-  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
-    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) :: 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) :: Ftheta(3*iim*jjm,llm)  ! 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")
-
-    IF(dysl_caldyn_coriolis) THEN
-
-#include "../kernels_hex/coriolis.k90"
-
-    ELSE
-#define FTHETA(ij) Ftheta(ij,1)
-
-    DO l=ll_begin, ll_end
-       ! compute theta flux
-       DO iq=1,nqdyn
-       !DIR$ SIMD
-          DO ij=ij_begin_ext,ij_end_ext      
-             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
-                                  * hflux(ij+u_right,l)
-             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
-                  * hflux(ij+u_lup,l)
-             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*FTHETA(ij+u_right)    +  &
-                  ne_rup*FTHETA(ij+u_rup)        +  &  
-                  ne_lup*FTHETA(ij+u_lup)        +  &  
-                  ne_left*FTHETA(ij+u_left)      +  &  
-                  ne_ldown*FTHETA(ij+u_ldown)    +  &
-                  ne_rdown*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,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
-                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
-                  wee(ij+u_right,5,2)*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,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
-                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
-                  wee(ij+u_lup,5,2)*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,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
-                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
-                  wee(ij+u_ldown,5,2)*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
-          ENDDO
-       ENDDO
-
-    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,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
-                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
-               .5*wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
-                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
-             uu_lup = &
-                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
-               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
-               .5*wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
-                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
-             uu_ldown = &
-                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
-               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
-               .5*wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
-                  wee(ij+u_ldown,5,2)*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
-          ENDDO
-       ENDDO
-
-    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,1)*hflux(ij+u_rup,l)+                           &
-                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
-                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
-                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
-                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
-                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
-                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
-                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
-                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
-                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
-             uu_lup = &
-                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
-                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
-                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
-                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
-                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
-                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
-                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
-                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
-                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
-                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
-             uu_ldown = &
-                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
-                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
-                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
-                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
-                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
-                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
-                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
-                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
-                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
-                  wee(ij+u_ldown,5,2)*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
-#undef FTHETA
-
-    END IF ! dysl
-
-    CALL trace_end("compute_caldyn_Coriolis")
-
-  END SUBROUTINE compute_caldyn_Coriolis
 
   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hv, hflux,Kv,du, zero)
Index: codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 844)
+++ codes/icosagcm/devel/src/dynamics/compute_caldyn_Coriolis.F90	(revision 844)
@@ -0,0 +1,219 @@
+MODULE compute_caldyn_Coriolis_mod
+  USE grid_param, ONLY : llm
+  IMPLICIT NONE
+  PRIVATE
+
+  PUBLIC :: compute_caldyn_Coriolis
+
+CONTAINS
+
+  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
+    USE icosa
+    USE trace
+    USE caldyn_vars_mod
+    USE omp_para, ONLY : ll_begin, ll_end
+    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) :: 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) :: Ftheta(3*iim*jjm,llm)  ! 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")
+
+    IF(dysl_caldyn_coriolis) THEN
+
+#include "../kernels_hex/coriolis.k90"
+
+    ELSE
+#define FTHETA(ij) Ftheta(ij,1)
+
+    DO l=ll_begin, ll_end
+       ! compute theta flux
+       DO iq=1,nqdyn
+       !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext      
+             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
+                                  * hflux(ij+u_right,l)
+             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
+                  * hflux(ij+u_lup,l)
+             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*FTHETA(ij+u_right)    +  &
+                  ne_rup*FTHETA(ij+u_rup)        +  &  
+                  ne_lup*FTHETA(ij+u_lup)        +  &  
+                  ne_left*FTHETA(ij+u_left)      +  &  
+                  ne_ldown*FTHETA(ij+u_ldown)    +  &
+                  ne_rdown*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,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
+                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
+                  wee(ij+u_right,5,2)*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,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
+                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
+                  wee(ij+u_lup,5,2)*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,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
+                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
+                  wee(ij+u_ldown,5,2)*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
+          ENDDO
+       ENDDO
+
+    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,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
+                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
+               .5*wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
+                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
+               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
+               .5*wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
+                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
+               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
+               .5*wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
+                  wee(ij+u_ldown,5,2)*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
+          ENDDO
+       ENDDO
+
+    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,1)*hflux(ij+u_rup,l)+                           &
+                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
+                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
+                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
+                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
+                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
+                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
+                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
+                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
+                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
+             uu_lup = &
+                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
+                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
+                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
+                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
+                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
+                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
+                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
+                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
+                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
+                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
+             uu_ldown = &
+                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
+                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
+                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
+                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
+                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
+                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
+                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
+                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
+                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
+                  wee(ij+u_ldown,5,2)*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
+#undef FTHETA
+
+    END IF ! dysl
+
+    CALL trace_end("compute_caldyn_Coriolis")
+
+  END SUBROUTINE compute_caldyn_Coriolis
+  
+END MODULE compute_caldyn_Coriolis_mod
