!--------------------------------------------------------------------------
   !---------------------------- compute_rhodz ----------------------------------
   IF(comp) THEN
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            m = ( AP(l,ij)-AP(l+1,ij) + (BP(l,ij)-BP(l+1,ij))*ps(ij) )/g
            rhodz(l,ij)=m
         END DO
      END DO
      !$OMP END DO
   ELSE
      err=0.
      !$OMP DO SCHEDULE(STATIC)
      DO ij = 1, primal_num
         !DIR$ SIMD
         DO l = 1, llm
            m = ( AP(l,ij)-AP(l+1,ij) + (BP(l,ij)-BP(l+1,ij))*ps(ij) )/g
            err = MAX(err, ABS(m-rhodz(l,ij)))
         END DO
      END DO
      !$OMP END DO
      IF(err>1e-10) THEN
         PRINT *, 'Discrepancy between ps and rhodz detected', err
         STOP
      END IF
   END IF
   !---------------------------- compute_rhodz ----------------------------------
   !--------------------------------------------------------------------------