function glm_psi_flux(psi, b_cc) result(gpf)
real, dimension(:,:), intent(in) :: psi (is it rank 2 ?)
real, dimension(size(psi,1),size(psi,2))              :: gpf

    if(has_ion) then

       gpf(1,:) = chspeed**2*(b_cc(xdim,:))

    else
       gpf(1,:) = zero

    end if
    apart from being intent(in) psi does not seem have any role