!--------------------------------------------------------------------------
   !---------------------------- compute_rhodz ----------------------------------
   IF(comp) THEN
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij=ij_begin, ij_end
            m = ( AP(ij,l)-AP(ij,l+1) + (BP(ij,l)-BP(ij,l+1))*ps(ij) )/g
            rhodz(ij,l)=m
         END DO
      END DO
   ELSE
      err=0.
      DO l = ll_begin, ll_end
         !DIR$ SIMD
         DO ij=ij_begin_ext, ij_end_ext
            m = ( AP(ij,l)-AP(ij,l+1) + (BP(ij,l)-BP(ij,l+1))*ps(ij) )/g
            err = MAX(err, ABS(m-rhodz(ij,l)))
         END DO
      END DO
      IF(err>1e-10) THEN
         PRINT *, 'Discrepancy between ps and rhodz detected', err
         STOP
      END IF
   END IF
   !---------------------------- compute_rhodz ----------------------------------
   !--------------------------------------------------------------------------