!--------------------------------------------------------------------------
   !---------------------------- caldyn_vert ----------------------------------
   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
   IF(caldyn_vert_variant == caldyn_vert_cons) THEN
      ! conservative vertical transport of momentum : (F/m)du/deta = 1/m (d/deta(Fu)-u.dF/deta)
      DO l = ll_beginp1, ll_end
         !DIR$ SIMD
         DO ij=ij_begin, ij_end
            wwuu(ij+u_right,l) = .25*(wflux(ij,l)+wflux(ij+t_right,l))*(u(ij+u_right,l)+u(ij+u_right,l-1)) ! Fu
            wwuu(ij+u_lup,l) = .25*(wflux(ij,l)+wflux(ij+t_lup,l))*(u(ij+u_lup,l)+u(ij+u_lup,l-1)) ! Fu
            wwuu(ij+u_ldown,l) = .25*(wflux(ij,l)+wflux(ij+t_ldown,l))*(u(ij+u_ldown,l)+u(ij+u_ldown,l-1)) ! Fu
         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
            dFu_deta = wwuu(ij+u_right,l+1)-wwuu(ij+u_right,l) ! d/deta (F*u)
            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_right,l+1)-(wflux(ij,l)+wflux(ij+t_right,l))) ! d/deta(F)
            du(ij+u_right,l) = du(ij+u_right,l) - (dFu_deta-u(ij+u_right,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_right,l))) ! (F/m)du/deta
            dFu_deta = wwuu(ij+u_lup,l+1)-wwuu(ij+u_lup,l) ! d/deta (F*u)
            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_lup,l+1)-(wflux(ij,l)+wflux(ij+t_lup,l))) ! d/deta(F)
            du(ij+u_lup,l) = du(ij+u_lup,l) - (dFu_deta-u(ij+u_lup,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))) ! (F/m)du/deta
            dFu_deta = wwuu(ij+u_ldown,l+1)-wwuu(ij+u_ldown,l) ! d/deta (F*u)
            dF_deta = .5*(wflux(ij,l+1)+wflux(ij+t_ldown,l+1)-(wflux(ij,l)+wflux(ij+t_ldown,l))) ! d/deta(F)
            du(ij+u_ldown,l) = du(ij+u_ldown,l) - (dFu_deta-u(ij+u_ldown,l)*dF_deta) / (.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))) ! (F/m)du/deta
         END DO
      END DO
   ELSE
      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
            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
   END IF
   !---------------------------- caldyn_vert ----------------------------------
   !--------------------------------------------------------------------------