, ele, dtheta_t)             do i=0,ele_ngi(velocity,ele)                 du_dz(0,i)=dot_product(grad_theta_gi(:,i),grav_at_quads(:,i)) ! Divide this by rho_0 for non-Boussinesq?             end do             