!--------------------------------------------------------------------------
   !---------------------------- caldyn_vert ----------------------------------
   DO l = llm-1,1,-1
      DO ij=ij_omp_begin,ij_omp_end
         ! cumulate mass flux convergence from top to bottom
         convm(ij,l) = convm(ij,l) + convm(ij,l+1)
      END DO
   END DO
   DO ij=ij_omp_begin,ij_omp_end
      dmass_col(ij) = convm(ij,1)
   END DO
   DO l = 2,llm
      DO ij=ij_omp_begin,ij_omp_end
         ! Compute vertical mass flux (l=1,llm+1 set to zero at init)
         wflux(ij,l) = mass_bl(l) * dmass_col(ij) - convm(ij,l)
      END DO
   END DO
   ! make sure wflux is up to date
   !$OMP BARRIER
   DO l = ll_beginp1, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         wwuu(ij+u_right,l) = .5*(wflux(ij,l)+wflux(ij+t_right,l))*(u(ij+u_right,l)-u(ij+u_right,l-1))
         wwuu(ij+u_lup,l) = .5*(wflux(ij,l)+wflux(ij+t_lup,l))*(u(ij+u_lup,l)-u(ij+u_lup,l-1))
         wwuu(ij+u_ldown,l) = .5*(wflux(ij,l)+wflux(ij+t_ldown,l))*(u(ij+u_ldown,l)-u(ij+u_ldown,l-1))
      END DO
   END DO
   ! make sure wwuu is up to date
   !$OMP BARRIER
   DO l = ll_begin, ll_end
      !DIR$ SIMD
      DO ij=ij_begin, ij_end
         convm(ij,l) = mass_dbk(l) * dmass_col(ij) ! FIXME : we shoud prognose mass_col
         du(ij+u_right,l) = du(ij+u_right,l) - (wwuu(ij+u_right,l)+wwuu(ij+u_right,l+1)) / (rhodz(ij,l)+rhodz(ij+t_right,l))
         du(ij+u_lup,l) = du(ij+u_lup,l) - (wwuu(ij+u_lup,l)+wwuu(ij+u_lup,l+1)) / (rhodz(ij,l)+rhodz(ij+t_lup,l))
         du(ij+u_ldown,l) = du(ij+u_ldown,l) - (wwuu(ij+u_ldown,l)+wwuu(ij+u_ldown,l+1)) / (rhodz(ij,l)+rhodz(ij+t_ldown,l))
      END DO
   END DO
   DO iq=1,nqdyn
      DO l = ll_beginp1, ll_end
         !DIR$ SIMD
         DO ij=ij_begin, ij_end
            dtheta_rhodz(ij,l,iq) = dtheta_rhodz(ij,l,iq) + 0.5*(theta(ij,l,iq)+theta(ij,l-1,iq))*wflux(ij,l)
         END DO
      END DO
      DO l = ll_begin, ll_endm1
         !DIR$ SIMD
         DO ij=ij_begin, ij_end
            dtheta_rhodz(ij,l,iq) = dtheta_rhodz(ij,l,iq) - 0.5*(theta(ij,l,iq)+theta(ij,l+1,iq))*wflux(ij,l+1)
         END DO
      END DO
   END DO
   !---------------------------- caldyn_vert ----------------------------------
   !--------------------------------------------------------------------------