Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 733)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 734)
@@ -51,5 +51,5 @@
     REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
     REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
-    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:)
+    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:)
 
 ! temporary shared variable
@@ -74,4 +74,5 @@
       CALL init_message(f_u,req_e1_vect,req_u)
       CALL init_message(f_qu,req_e1_scal,req_qu)
+      IF(caldyn_kinetic==kinetic_consistent) CALL init_message(f_Kv,req_z1_scal,req_Kv)
       IF(.NOT.hydrostatic) THEN
          CALL init_message(f_geopot,req_i1,req_geopot)
@@ -140,9 +141,18 @@
        qv=f_qv(ind)
        CALL compute_pvort_only(u,mass,qu,qv)
+       IF(caldyn_kinetic==kinetic_consistent) THEN
+          Kv=f_Kv(ind)
+          CALL compute_caldyn_Kv(u,Kv)
+       END IF
     ENDDO
     
     CALL send_message(f_qu,req_qu) ! COM03
     CALL wait_message(req_qu) ! COM03
-       
+
+    IF(caldyn_kinetic==kinetic_consistent) THEN
+       CALL send_message(f_Kv,req_Kv)
+       CALL wait_message(req_Kv)
+    END IF
+
     DO ind=1,ndomain
        IF (.NOT. assigned_domain(ind)) CYCLE
@@ -159,5 +169,6 @@
 
        IF(hydrostatic) THEN
-          CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.)
+          Kv=f_Kv(ind)
+          CALL compute_caldyn_slow_hydro(u,mass,hflux,Kv,du, .TRUE.)
        ELSE
           W = f_W(ind)
Index: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 733)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 734)
@@ -13,5 +13,5 @@
   LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
 
-  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, &
+  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Kv, compute_caldyn_Coriolis, &
        compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
        compute_caldyn_solver, compute_caldyn_fast
@@ -99,4 +99,48 @@
 
   END SUBROUTINE compute_pvort_only
+
+  SUBROUTINE compute_caldyn_kv(ue, Kv)
+    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: Kv(2*iim*jjm,llm)
+    REAL(rstd) :: ue2(3*iim*jjm), dem2(3*iim*jjm), r2_Av(2*iim*jjm), rad2
+    INTEGER :: ij,l, u_up, u_down
+    
+    u_up    = t_lup + u_right
+    u_down  = t_rdown + u_left
+    
+    rad2=radius**2
+
+    !DIR$ SIMD
+    DO ij=ij_begin_ext,ij_end_ext
+       dem2(ij+u_right) = de(ij+u_right)**(-2)
+       dem2(ij+u_lup)   = de(ij+u_lup)**(-2)
+       dem2(ij+u_ldown) = de(ij+u_ldown)**(-2)
+       r2_Av(ij+z_up)   = rad2*(1./Av(ij+z_up))
+       r2_Av(ij+z_down) = rad2*(1./Av(ij+z_down))
+    END DO
+
+    DO l=ll_begin,ll_end
+       ! compute squared normal component from 1-form
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          ue2(ij+u_right) = dem2(ij+u_right)* (ue(ij+u_right,l)**2)
+          ue2(ij+u_lup)   = dem2(ij+u_lup)  * (ue(ij+u_lup,l)**2)
+          ue2(ij+u_ldown) = dem2(ij+u_ldown)* (ue(ij+u_ldown,l)**2)
+       END DO
+       ! average squared normal component to vertices
+       !DIR$ SIMD
+       DO ij=ij_begin_ext,ij_end_ext
+          Kv(ij+z_up,l) = r2_Av(ij+z_up)*(                       &
+                               S1(ij,vup)*ue2(ij+u_rup) +        &
+                               S2(ij,vup)*ue2(ij+u_lup) +        &
+                               S2(ij+t_lup,vrdown)*ue2(ij+u_up))
+
+          Kv(ij+z_down,l) = r2_Av(ij+z_down)*(                   &
+                               S1(ij,vdown)*ue2(ij+u_ldown) +    &
+                               S2(ij,vdown)*ue2(ij+u_rdown) +    &
+                               S2(ij+t_rdown,vlup)*ue2(ij+u_down) )
+       ENDDO
+    ENDDO
+  END SUBROUTINE compute_caldyn_kv
 
   SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
@@ -663,7 +707,7 @@
                   wee(ij+u_ldown,5,2)*hflux(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
+             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
@@ -679,7 +723,8 @@
   END SUBROUTINE compute_caldyn_Coriolis
 
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero)
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,Kv,du, zero)
     LOGICAL, INTENT(IN) :: zero
     REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
+    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
@@ -719,14 +764,26 @@
        ENDDO
        ! Compute Bernoulli=kinetic energy
-       !DIR$ SIMD
-       DO ij=ij_begin,ij_end         
-          BERNI(ij) = &
-               1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
-               le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
-               le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
-               le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
-               le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
-               le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
-       ENDDO
+       IF(caldyn_kinetic==kinetic_trisk) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end         
+             BERNI(ij) = &
+                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
+                                le_de(ij+u_rup)*u(ij+u_rup,l)**2     +    &
+                                le_de(ij+u_lup)*u(ij+u_lup,l)**2     +    &
+                                le_de(ij+u_left)*u(ij+u_left,l)**2   +    &
+                                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
+                                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )  
+          ENDDO
+       ELSE
+          !DIR$ SIMD
+          DO ij=ij_begin,ij_end
+             BERNI(ij) = Riv(ij,vup)   *Kv(ij+z_up,l)    + &
+                         Riv(ij,vlup)  *Kv(ij+z_lup,l)   + &
+                         Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
+                         Riv(ij,vdown) *Kv(ij+z_down,l)  + &
+                         Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
+                         Riv(ij,vrup)  *Kv(ij+z_rup,l)
+          END DO
+       END IF
        ! Compute du=-grad(Bernoulli)
        IF(zero) THEN
Index: codes/icosagcm/devel/src/dynamics/caldyn_vars.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 733)
+++ codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 734)
@@ -26,5 +26,5 @@
                            f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:)
 
-  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w
+  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w, req_Kv
 
 END MODULE caldyn_vars_mod
