Index: codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 734)
+++ codes/icosagcm/devel/src/dynamics/caldyn_gcm.F90	(revision 735)
@@ -61,9 +61,11 @@
     CASE('energy')
        caldyn_conserv=conserv_energy
+    CASE('energy_gassmann')
+       caldyn_conserv=conserv_gassmann
     CASE('enstrophy')
        caldyn_conserv=conserv_enstrophy
     CASE DEFAULT
        IF (is_mpi_root) PRINT *,'Bad selector for variable caldyn_conserv : <', &
-            TRIM(def),'> options are <energy>, <enstrophy>'
+            TRIM(def),'> options are <energy>, <energy_gassmann>, <enstrophy>'
        STOP
     END SELECT
@@ -156,9 +158,10 @@
 
   SUBROUTINE allocate_caldyn
-    CALL allocate_field(f_qu,field_u,type_real,llm) 
-    CALL allocate_field(f_qv,field_z,type_real,llm) 
-    CALL allocate_field(f_Kv,field_z,type_real,llm) 
+    CALL allocate_field(f_qu,field_u,type_real,llm, name='qu')
+    CALL allocate_field(f_qv,field_z,type_real,llm, name='qv')
+    CALL allocate_field(f_Kv,field_z,type_real,llm, name='Kv')
+    CALL allocate_field(f_hv,field_z,type_real,llm, name='hv')
     CALL allocate_field(f_pk,    field_t,type_real,llm,  name='pk')
-    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')    
+    CALL allocate_field(f_wwuu,  field_u,type_real,llm+1,name='wwuu')
     CALL allocate_field(f_planetvel,  field_u,type_real, name='planetvel') ! planetary velocity at r=a
     IF(.NOT.hydrostatic) THEN
Index: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 734)
+++ codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90	(revision 735)
@@ -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(:,:), Kv(:,:)
+    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:), hv(:,:)
 
 ! temporary shared variable
@@ -140,5 +140,6 @@
        qu=f_qu(ind)
        qv=f_qv(ind)
-       CALL compute_pvort_only(u,mass,qu,qv)
+       hv=f_hv(ind)
+       CALL compute_pvort_only(u,mass,qu,qv,hv)
        IF(caldyn_kinetic==kinetic_consistent) THEN
           Kv=f_Kv(ind)
@@ -151,4 +152,5 @@
 
     IF(caldyn_kinetic==kinetic_consistent) THEN
+       CALL transfert_request(f_hv,req_z1_scal)
        CALL send_message(f_Kv,req_Kv)
        CALL wait_message(req_Kv)
@@ -169,6 +171,7 @@
 
        IF(hydrostatic) THEN
+          hv=f_hv(ind)
           Kv=f_Kv(ind)
-          CALL compute_caldyn_slow_hydro(u,mass,hflux,Kv,du, .TRUE.)
+          CALL compute_caldyn_slow_hydro(u,mass,hv, 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 734)
+++ codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90	(revision 735)
@@ -50,9 +50,10 @@
   END SUBROUTINE compute_theta
 
-  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv)
+  SUBROUTINE compute_pvort_only(u,rhodz,qu,qv,hv_)
     REAL(rstd),INTENT(IN)  :: u(iim*3*jjm,llm)
     REAL(rstd),INTENT(INOUT) :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: qu(iim*3*jjm,llm)
     REAL(rstd),INTENT(OUT) :: qv(iim*2*jjm,llm)
+    REAL(rstd),INTENT(OUT) :: hv_(iim*2*jjm,llm)
 
     INTEGER :: ij,l
@@ -76,4 +77,5 @@
                + Riv2(ij+t_lup,vrdown) * rhodz(ij+t_lup,l)
           qv(ij+z_up,l) = ( etav+fv(ij+z_up) )/hv
+          hv_(ij+z_up,l) = hv
           
           etav = 1./Av(ij+z_down)*(  ne_ldown * u(ij+u_ldown,l)   &
@@ -84,4 +86,5 @@
                + Riv2(ij+t_rdown,vlup) * rhodz(ij+t_rdown,l)
           qv(ij+z_down,l) =( etav+fv(ij+z_down) )/hv
+          hv_(ij+z_down,l) = hv
        ENDDO
 
@@ -668,4 +671,48 @@
        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
 
@@ -723,8 +770,9 @@
   END SUBROUTINE compute_caldyn_Coriolis
 
-  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,Kv,du, zero)
+  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hv, 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)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
     REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
     REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
@@ -751,16 +799,33 @@
        !  Compute mass fluxes
        IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu) 
-       !DIR$ SIMD
-       DO ij=ij_begin_ext,ij_end_ext
-          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
-          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
-          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
-          uu_right= uu_right*le_de(ij+u_right)
-          uu_lup  = uu_lup  *le_de(ij+u_lup)
-          uu_ldown= uu_ldown*le_de(ij+u_ldown)
-          hflux(ij+u_right,l)=uu_right
-          hflux(ij+u_lup,l)  =uu_lup
-          hflux(ij+u_ldown,l)=uu_ldown
-       ENDDO
+
+       IF(caldyn_kinetic==kinetic_trisk) THEN
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
+             uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
+             uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
+             uu_right= uu_right*le_de(ij+u_right)
+             uu_lup  = uu_lup  *le_de(ij+u_lup)
+             uu_ldown= uu_ldown*le_de(ij+u_ldown)
+             hflux(ij+u_right,l)=uu_right
+             hflux(ij+u_lup,l)  =uu_lup
+             hflux(ij+u_ldown,l)=uu_ldown
+          ENDDO
+       ELSE ! mass flux deriving from consistent kinetic energy
+          !DIR$ SIMD
+          DO ij=ij_begin_ext,ij_end_ext
+             uu_right=0.5*(hv(ij+z_rup,l)+hv(ij+z_rdown,l))*u(ij+u_right,l)
+             uu_lup=0.5*(hv(ij+z_up,l)+hv(ij+z_lup,l))*u(ij+u_lup,l)
+             uu_ldown=0.5*(hv(ij+z_ldown,l)+hv(ij+z_down,l))*u(ij+u_ldown,l)
+             uu_right= uu_right*le_de(ij+u_right)
+             uu_lup  = uu_lup  *le_de(ij+u_lup)
+             uu_ldown= uu_ldown*le_de(ij+u_ldown)
+             hflux(ij+u_right,l)=uu_right
+             hflux(ij+u_lup,l)  =uu_lup
+             hflux(ij+u_ldown,l)=uu_ldown
+          ENDDO
+       END IF
+
        ! Compute Bernoulli=kinetic energy
        IF(caldyn_kinetic==kinetic_trisk) THEN
Index: codes/icosagcm/devel/src/dynamics/caldyn_vars.f90
===================================================================
--- codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 734)
+++ codes/icosagcm/devel/src/dynamics/caldyn_vars.f90	(revision 735)
@@ -6,5 +6,5 @@
   SAVE
 
-  INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, &
+  INTEGER, PARAMETER :: conserv_energy=1, conserv_enstrophy=2, conserv_gassmann=3, &
        kinetic_trisk=1, kinetic_consistent=2, &
        caldyn_vert_noncons=1, caldyn_vert_cons=2, &
@@ -22,5 +22,5 @@
 
   ! temporary shared variables for caldyn
-  TYPE(t_field),POINTER :: f_qu(:), f_qv(:), f_Kv(:), &
+  TYPE(t_field),POINTER :: f_qu(:), f_qv(:), f_Kv(:), f_hv(:), &
                            f_pk(:),f_wwuu(:),f_planetvel(:), &
                            f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:)
