physical_prop.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP                                           C
!  Purpose: Calculate the indicated physical properties that vary      C
!           with time if directed to do so by the corresponding flag   C
!                                                                      C
!  Author: M. Syamlal                                 Date: 17-JUL-92  C
!  Reviewer: P. Nicoletti                             Date: 11-DEC-92  C
!                                                                      C
!  Revision: Mods for MFIX 2.0 (old name CALC_PHYSPROP)                C
!  Author: M. Syamlal                                 Date: 23-APR-96  C
!                                                                      C
!                                                                      C
!  Literature/Document References:                                     C
!    Perry, R.H., and C.H. Chilton, Chemical Engineer's Handbook, 5th  C
!      edition, McGraw-Hill Kogakusha, Tokyo, 1973.                    C
!                                                                      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop(IER, LEVEL)

! Modules
!---------------------------------------------------------------------//
      use compar, only: numpes, mype, pe_io
      use funits, only: unit_log
      use exit, only: mfix_exit
      use mpi_utility, only: global_all_sum
      use physprop, only: smax
      use coeff, only: density
      use coeff, only: sp_heat
      use coeff, only: psize
      implicit none

! Dummy arguments
!---------------------------------------------------------------------//
! Global error Flag.
      INTEGER, intent(inout) :: IER
      INTEGER, intent(in) :: LEVEL

! Local variables
!---------------------------------------------------------------------//
! Arrays for storing errors:
! 100 - Negative gas phase density
! 101 - Negative solids phase density
! 10x - Unclassified
! 900 - Invalid temperature in calc_CpoR
      INTEGER :: Err_l(0:numpes-1)  ! local
      INTEGER :: Err_g(0:numpes-1)  ! global
! loop index
      INTEGER :: M
!......................................................................!

! Initialize error flags.
      err_l = 0

! Calculate density only. This is invoked several times within iterate,
! making it the most frequently called.
      if(level == 0) then
         if(density(0)) CALL physical_prop_rog
           DO m=1,smax
              if(density(m)) CALL physical_prop_ros(m)
           ENDDO

! Calculate everything except density. This is called at the start of
! each iteration.
      elseif(level == 1) then
         if(sp_heat(0)) CALL physical_prop_cpg
         DO m=1,smax
            if(sp_heat(m)) CALL physical_prop_cps(m)
            if(psize(m)) CALL physical_prop_dp(m)
         ENDDO


! Calculate everything. This is invoked via calc_coeff_all as part of
! the initialization (before starting the time march) and at the start
! of each step step thereafter.
      elseif(level == 2) then
         if(density(0)) CALL physical_prop_rog
         if(sp_heat(0)) CALL physical_prop_cpg
         DO m=1,smax
            if(density(m)) CALL physical_prop_ros(m)
            if(sp_heat(m)) CALL physical_prop_cps(m)
            if(psize(m)) CALL physical_prop_dp(m)
         ENDDO
      endif


! In case of negative density force exit from the physical property
! calculation routine and reduce the time step
      CALL global_all_sum(err_l, err_g)
      ier = maxval(err_g)
      if(ier == 0) return


! Error handeling. - Local.
! An invalid temperature was found by calc_CpoR. This is a fatal run-
! time error and forces a call to MFIX_EXIT.
      IF(ier == 901 .OR. ier == 902) then
! Temperature is now bounded so this error will not occur.
         if(mype == pe_io) then
            write(*,2000) ier
            write(unit_log,2000) ier
         endif
         CALL mfix_exit(mype)
      ENDIF

      return

 2000 FORMAT(/1x,70('*')/' From: PHYSICAL_PROP',/' Fatal Error 2000:', &
         ' calc_CpoR reported an invalid temperature: 0x0', i3/,       &
         'See Cp.log for details. Calling MFIX_EXIT.',/1x,70('*')/)

      CONTAINS

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP_ROg                                       C
!  Purpose: Calculate the gas phase density.                           C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop_rog

! Global variables:
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
! Equation of State - GAS
      use eos, only: eosg
! Gas phase species mass fractions.
      use fldvar, only: x_g
! Gas phase temperature.
      use fldvar, only: t_g
! Gas phase density (compressible).
      use fldvar, only: ro_g
! Gas phase pressure.
      use fldvar, only: p_g
! Gas phase volume fraction.
      use fldvar, only: ep_g
! Gas phase material density.
      use fldvar, only: rop_g

      use functions, only: wall_at
! Run time flag for generating negative gas density log files
      use output, only: report_neg_density
      use param1, only: undefined, zero, one
      use physprop, only: database_read, mw_g, mw_mix_g, mw_avg, nmax
! Maximum value for molecular weight (divided by one)
      use toleranc, only: omw_max
! Detect NaN in density for compressible flows.
      use utilities, only: mfix_isnan
! invoke user defined quantity
      USE usr_prop, only: usr_rog, calc_usr_prop
      USE usr_prop, only: gas_density
      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Average molecular weight
      DOUBLE PRECISION :: MW
! Loop indicies
      INTEGER :: IJK   ! Computational cell
! Flag to write log header
      LOGICAL :: wHeader
!......................................................................!

! Ensure that the database was read. This *should* have been caught by
! check_gas_phase but this call remains to prevent an accident.
      IF(.NOT.database_read) call read_database0(ier)

! User-defined function
      IF(usr_rog) THEN
         CALL calc_usr_prop(gas_density,lm=0,lerr=err_l)
         RETURN
      ENDIF

! Initialize:
      wheader = .true.

      ijk_lp: DO ijk = ijkstart3, ijkend3
         IF(wall_at(ijk)) cycle ijk_lp
         IF (mw_avg == undefined) THEN
! Calculating the average molecular weight of the fluid.
            mw = sum(x_g(ijk,:nmax(0))/mw_g(:nmax(0)))
            mw = one/max(mw,omw_max)
            mw_mix_g(ijk) = mw
! Calculate the fluid density and bulk density
            ro_g(ijk) = eosg(mw,p_g(ijk),t_g(ijk))
            rop_g(ijk) = ro_g(ijk)*ep_g(ijk)
         ELSE
            ro_g(ijk) = eosg(mw_avg,p_g(ijk),t_g(ijk))
            rop_g(ijk) = ro_g(ijk)*ep_g(ijk)
         ENDIF

         IF(ro_g(ijk) < zero .or. mfix_isnan(ro_g(ijk))) THEN
            err_l(mype) = 100
            IF(report_neg_density) CALL rogerr_log(ijk, wheader)
         ENDIF
      ENDDO ijk_lp

      RETURN
      END SUBROUTINE physical_prop_rog


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP_ROs                                       C
!  Purpose: Calculate solids phase (variable) density.                 C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop_ros(M)

! Global variables:
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
! Equation of State - Solid
      use eos, only: eoss
! Solid phase species mass fractions.
      use fldvar, only: x_s, rop_s, ro_s
! Solid phase density (variable).
      use fldvar, only: rop_s, ro_s

      use functions, only: wall_at
! Run time flag for generating negative density log files.
      use output, only: report_neg_density
      use param1, only: undefined, zero
! Baseline/Unreaced solids density
      use physprop, only: ro_s0
! Initial mass fraction of inert species
      use physprop, only: x_s0
! Index of inert solids phase species.
      use physprop, only: inert_species
! Inert solids phase species mass fraction in dilute region.
      use physprop, only: dil_inert_x_vsd
! Factor to define dilute region where DIL_INERT_X_VSD is used
      use physprop, only: dil_factor_vsd
! Flag for variable solids density.
      use run, only: solve_ros
! Minimum solids volume fraction
      use toleranc, only: dil_ep_s
! invoke user defined quantity
      USE usr_prop, only: usr_ros, calc_usr_prop
      USE usr_prop, only: solids_density
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local Variables:
!---------------------------------------------------------------------//
! Loop indicies
      INTEGER :: IJK
! Index of inert species
      INTEGER :: IIS
! Flag to write log header
      LOGICAL :: wHeader
! Minimum bulk density
      DOUBLE PRECISION :: minROPs
!......................................................................!

! User defined function
      IF(usr_ros(m)) THEN
         CALL calc_usr_prop(solids_density,lm=m,lerr=err_l)
         RETURN
      ENDIF

      IF(solve_ros(m)) THEN
! Initialize header flag.
         wheader = .true.
! Set the index of the inert species
         iis = inert_species(m)
! Calculate the minimum solids denisty.
         minrops = ro_s0(m)*(dil_factor_vsd*dil_ep_s)

! Calculate the solids denisty over all cells.
         ijk_lp: DO ijk = ijkstart3, ijkend3
            IF(wall_at(ijk)) cycle ijk_lp
            IF(rop_s(ijk,m) > minrops) THEN
               ro_s(ijk,m) = eoss(ro_s0(m), x_s0(m,iis), &
                  x_s(ijk,m,iis))
            ELSE
!               RO_s(IJK,M) = RO_s0(M)
               ro_s(ijk,m) = eoss(ro_s0(m), x_s0(m,iis), &
                  dil_inert_x_vsd(m))
            ENDIF

! Report errors.
            IF(ro_s(ijk,m) <= zero) THEN
               err_l(mype) = 101
               IF(report_neg_density) CALL roserr_log(ijk, m, wheader)
            ENDIF
         ENDDO ijk_lp
      ENDIF

      RETURN
      END SUBROUTINE physical_prop_ros


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP_CPg                                       C
!  Purpose: Calculate the gas phase constant pressure specific heat.   C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
! Notes:                                                               C
!  > Unit conversion: 1 cal = 4.183925 J                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop_cpg

! Global Variables:
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
! Universal gas constant in cal/mol.K
      use constant, only: rgas => gas_const_cal
! Gas phase species mass fractions.
      use fldvar, only: x_g
! Gas phase temperature.
      use fldvar, only: t_g

      use functions, only: wall_at
      use output, only: report_neg_specificheat
      use param1, only: undefined, zero
      use physprop, only: database_read
      use physprop, only: mw_g, c_pg, nmax
! Function to calculate Cp over gas constant R
      use read_thermochemical, only: calc_cpor
! Units: CGS/SI
      use run, only: units
! invoke user defined quantity
      USE usr_prop, only: usr_cpg, calc_usr_prop
      USE usr_prop, only: gas_specificheat
      IMPLICIT NONE

! Local Variables:
!---------------------------------------------------------------------//
! Species specific heat.
      DOUBLE PRECISION :: lCp
! Loop indicies
      INTEGER :: IJK, NN
! Flag to write log header
      LOGICAL :: wHeader
! Error flag returned from calc_CpoR
      INTEGER :: lCP_Err
      INTEGER :: gCP_Err
!......................................................................!

! Ensure that the database was read. This *should* have been caught by
! check_gas_phase but this call remains to prevent an accident.
      IF(.NOT.database_read) CALL read_database0(ier)

! User defined function
      IF(usr_cpg) THEN
         CALL calc_usr_prop(gas_specificheat,lm=0,lerr=err_l)
         RETURN
      ENDIF

! Initialize
      wheader = .true.
      gcp_err = 0
      lcp_err = 0
      ijk_lp: DO ijk = ijkstart3, ijkend3
         IF(wall_at(ijk)) cycle ijk_lp
! Calculating an average specific heat for the fluid.
         c_pg(ijk) = zero
         DO nn = 1, nmax(0)
            lcp = calc_cpor(t_g(ijk), 0, nn)
! calc_cpor used to return error flag (101/102) on out of bound
! temperature. however, temperature is now bounded in routine.
            gcp_err = max(gcp_err, lcp_err)
            c_pg(ijk) = c_pg(ijk) + x_g(ijk,nn) * lcp * rgas / mw_g(nn)
         ENDDO

! report errors
         IF(c_pg(ijk) <= zero) THEN
!            gCP_err = 103
            IF(report_neg_specificheat) CALL cpgerr_log(ijk, wheader)
         ENDIF

      ENDDO ijk_lp

! The database calculation always returns cal/g.K thus the following
! conversion is needed if using SI units.
! Only convert useful part of the array. When the array is re-indexed,
! elements past IJKEND3 are UNDEFINED and would keep growing if
! multiplied, leading to overflow.
      IF(units == 'SI') THEN
         DO ijk = ijkstart3, ijkend3
           c_pg(ijk) = 4.183925d3 * c_pg(ijk)
         ENDDO
      ENDIF

! Increment the error to 900+ to invoke fatal exit.
      IF(gcp_err /= 0) err_l(mype) = 800 + gcp_err

      RETURN
      END SUBROUTINE physical_prop_cpg


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP_CPs                                       C
!  Purpose: Calculate solids phase constant pressure specific heat.    C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop_cps(M)

! Global Variables:
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
! Universal gas constant in cal/mol.K
      use constant, only: rgas => gas_const_cal
! Solids temperature
      use fldvar, only: t_s
! Solids species mass fractions
      use fldvar, only: x_s

      use functions, only: wall_at
      use output, only: report_neg_specificheat
      use param1, only: undefined, zero
      use physprop, only: database_read
      use physprop, only: mw_s, c_ps, nmax
! Function to calculate Cp over gas constant R
      use read_thermochemical, only: calc_cpor
! Units: CGS/SI
      use run, only: units
! invoke user defined quantity
      USE usr_prop, only: usr_cps, calc_usr_prop
      USE usr_prop, only: solids_specificheat
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables:
!---------------------------------------------------------------------//
! Local value for species specific heat (Cp)
      DOUBLE PRECISION :: lCp
! Loop indicies
      INTEGER :: IJK, NN
! Flag to write log header
      LOGICAL :: wHeader
! Error flag returned from calc_CpoR
      INTEGER :: lCP_Err
      INTEGER :: gCP_Err
!......................................................................!

! Ensure that the database was read. This *should* have been caught by
! check_solids_common_all but this call remains to prevent an accident.
      IF(.NOT.database_read) CALL read_database0(ier)

! User defined function
      IF(usr_cps(m)) THEN
         CALL calc_usr_prop(solids_specificheat,lm=m,lerr=err_l)
         RETURN
      ENDIF

! Initialize
      wheader = .true.
      gcp_err = 0
      lcp_err = 0
      ijk_lp: DO ijk = ijkstart3, ijkend3
         IF(wall_at(ijk)) cycle ijk_lp
! Calculating an average specific heat for the fluid.
         c_ps(ijk, m) = zero
         DO nn = 1, nmax(m)
            lcp = calc_cpor(t_s(ijk,m), m, nn)
! calc_cpor used to return error flag (101/102) based on out of bound
! temperature. however, temperature is now bounded in routine
            gcp_err = max(gcp_err, lcp_err)
            c_ps(ijk,m) = c_ps(ijk,m) + x_s(ijk,m,nn) * &
               (lcp * rgas / mw_s(m,nn))
         ENDDO

! report errors
         IF(c_ps(ijk,m) <= zero) THEN
!            gCP_err = 104
            IF(report_neg_specificheat) CALL cpserr_log(ijk, m, wheader)
         ENDIF

      ENDDO ijk_lp

! The database calculation always returns cal/g.K thus the following
! conversion is needed if using SI units.
! Only convert useful part of the array. When the array is re-indexed,
! elements past IJKEND3 are UNDEFINED and would keep growing if
! multiplied, leading to overflow.
      IF(units == 'SI') THEN
         DO ijk = ijkstart3, ijkend3
           c_ps(ijk,m) = 4.183925d3 * c_ps(ijk,m)
         ENDDO
      ENDIF

! Increment the error to 900+ to invoke fatal exit.
      IF(gcp_err /= 0) err_l(mype) = 800 + gcp_err

      RETURN
      END SUBROUTINE physical_prop_cps


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: PHYSICAL_PROP_DP                                        C
!  Purpose: Calculate solids phase diameter.                           C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE physical_prop_dp(M)

! Modules
!---------------------------------------------------------------------//
      use compar, only: ijkstart3, ijkend3
      use scalars, only: phase4scalar
      use fldvar, only: scalar
      use fldvar, only: d_p, ep_s
      use functions, only: wall_at
      use param1, only: small_number
      use run, only: call_dqmom
      implicit none

! Dummy arguments
!---------------------------------------------------------------------//
! solids phase index
      INTEGER, INTENT(IN) :: M

! Local variables
!---------------------------------------------------------------------//
! Loop indicies
      INTEGER :: IJK   ! Computational cell
! Map from true index to map.
      INTEGER :: lM
!......................................................................!

      IF(.NOT.call_dqmom) return

      lm = phase4scalar(m) ! Map from scalar eq to solids phase

      ijk_lp: DO ijk = ijkstart3, ijkend3
         IF(wall_at(ijk)) cycle ijk_lp
         IF(ep_s(ijk,lm) > small_number) d_p(ijk,m)= scalar(ijk,lm)
      ENDDO ijk_lp

      RETURN
      END SUBROUTINE physical_prop_dp


      END SUBROUTINE physical_prop

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: ROgErr_LOG                                              C
!  Purpose: Record information about the location and conditions that  C
!           resulted in a negative gas phase density.                  C
!                                                                      C
!  Author: J. Musser                                  Date: 28-JUN-13  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE rogerr_log(IJK, tHeader)

! Global Variables:
!---------------------------------------------------------------------//
      use compar, only: mype, numpes
! Cutcell data.
      use cutcell, only: cartesian_grid
      use cutcell, only: cut_cell_at, small_cell_at
      use cutcell, only: xg_e, yg_n, zg_t
! Simulation time
      use run, only: time
! Gas phase temperature.
      use fldvar, only: t_g
! Gas phase density (compressible).
      use fldvar, only: ro_g
!
      use indices, only: i_of, j_of, k_of
! Gas phase pressure.
      use fldvar, only: p_g
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
      INTEGER, intent(in) :: IJK
      LOGICAL, intent(inout) :: tHeader

! Local variables:
!---------------------------------------------------------------------//
      LOGICAL :: lExists
      CHARACTER(LEN=255) :: lFile
      INTEGER, parameter :: lUnit = 4868
      LOGICAL, save :: fHeader = .true.
!......................................................................!

      lfile = '';
      if(numpes > 1) then
         write(lfile,"('ROgErr_',I4.4,'.log')") mype
      else
         write(lfile,"('ROgErr.log')")
      endif
      inquire(file=trim(lfile),exist=lexists)
      if(lexists) then
         open(lunit,file=trim(adjustl(lfile)),  &
            status='old', position='append', convert='big_endian')
      else
         open(lunit,file=trim(adjustl(lfile)), status='new',&
              convert='big_endian')
      endif

      if(fheader) then
         write(lunit,1000)
         fheader = .false.
      endif

      if(theader) then
         write(lunit,"(/2x,'Simulation time: ',g12.5)") time
         theader = .false.
      endif

      write(lunit,1001) ijk, i_of(ijk), j_of(ijk), k_of(ijk)
      write(lunit,"(6x,A,1X,g12.5)",advance='NO') 'RO_g:', ro_g(ijk)
      write(lunit,"(2x,A,1X,g12.5)",advance='NO') 'P_g:', p_g(ijk)
      write(lunit,"(2x,A,1X,g12.5)") 'T_g:', t_g(ijk)
      if(cartesian_grid) then
         write(lunit,"(6x,A,1X,L1)",advance='NO') 'Cut Cell:', cut_cell_at(ijk)
         write(lunit,"(6x,A,1X,L1)") 'Small Cell:', small_cell_at(ijk)
         write(lunit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
            xg_e(i_of(ijk)), yg_n(j_of(ijk)), zg_t(k_of(ijk))
      endif

      close(lunit)

      RETURN

 1000 FORMAT(2x,'One or more cells have reported a negative/NaN gas',  &
         ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
         ' lower UR_FAC(1) in mfix.dat.')

 1001 FORMAT(/4x,'IJK: ',i8,7x,'I: ',i4,'  J: ',i4,'  K: ',i4)

      END SUBROUTINE rogerr_log


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: ROsErr_LOG                                              !
!  Purpose: Record information about the location and conditions that  !
!           resulted in a negative solids phase density.               !
!                                                                      !
!  Author: J. Musser                                  Date: 09-Oct-13  !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE roserr_log(IJK, M, tHeader)

! Global variables:
!---------------------------------------------------------------------//
      use compar, only: mype, numpes
! Cutcell data.
      use cutcell, only: cartesian_grid
      use cutcell, only: cut_cell_at, small_cell_at
      use cutcell, only: xg_e, yg_n, zg_t
! Solid phase species mass fractions.
      use fldvar, only: x_s
! Solid phase density (variable).
      use fldvar, only: ro_s
!
      use indices, only: i_of, j_of, k_of
! Baseline/Unreaced solids density
      use physprop, only: ro_s0
! Initial mass fraction of inert species
      use physprop, only: x_s0
! Index of inert solids phase species.
      use physprop, only: inert_species
! Simulation time
      use run, only: time
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Fluid cell and solids phase indices
      INTEGER, intent(in) :: IJK, M
! Flag to output header
      LOGICAL, intent(inout) :: tHeader

! Local Variables:
!---------------------------------------------------------------------//
! Local aliase for inert species index
      INTEGER :: NN
! Local file values.
      LOGICAL :: lExists
      CHARACTER(LEN=255) :: lFile
      INTEGER, parameter :: lUnit = 4868
      LOGICAL, save :: fHeader = .true.
!......................................................................!

      lfile = '';
      if(numpes > 1) then
         write(lfile,"('ROsErr_',I4.4,'.log')") mype
      else
         write(lfile,"('ROsErr.log')")
      endif
      inquire(file=trim(lfile),exist=lexists)
      if(lexists) then
         open(lunit,file=trim(adjustl(lfile)), &
            status='old', position='append',convert='big_endian')
      else
         open(lunit,file=trim(adjustl(lfile)), status='new',&
              convert='big_endian')
      endif

      if(fheader) then
         write(lunit,1000)
         fheader = .false.
      endif

      if(theader) then
         write(lunit,"(/2x,'Simulation time: ',g12.5,'  Phase: ',I2)")&
             time, m
         theader = .false.
      endif

      nn = inert_species(m)
      write(lunit,1001) ijk, i_of(ijk), j_of(ijk), k_of(ijk)
      write(lunit,"(6x,A,1X,g12.5)",advance='no') 'RO_s:', ro_s(ijk,m)
      write(lunit,"(2x,A,1X,g12.5)",advance='no') 'Base:', ro_s0(m)
      write(lunit,"(2x,A,1X,g12.5)",advance='no') 'X_s0:', x_s0(m,nn)
      write(lunit,"(2x,A,1X,g12.5)") 'X_s:', x_s(ijk,m,nn)

      if(cartesian_grid) then
         write(lunit,"(6x,A,1X,L1)",advance='NO') 'Cut Cell:', cut_cell_at(ijk)
         write(lunit,"(6x,A,1X,L1)") 'Small Cell:', small_cell_at(ijk)
         write(lunit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
            xg_e(i_of(ijk)), yg_n(j_of(ijk)), zg_t(k_of(ijk))
      endif

      close(lunit)

      RETURN

 1000 FORMAT(2x,'One or more cells have reported a negative solids',   &
         ' density (RO_s(IJK)).')

 1001 FORMAT( 4x,'IJK: ',i8,7x,'I: ',i4,'  J: ',i4,'  K: ',i4)

      END SUBROUTINE roserr_log

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CPgErr_LOG                                              !
!  Purpose: Record information about the location and conditions that  !
!           resulted in a zero or negative gas specific heat.          !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE cpgerr_log(IJK, tHeader)

! Global variables:
!---------------------------------------------------------------------//
      use compar, only: mype, numpes
! Cutcell data.
      use cutcell, only: cartesian_grid
      use cutcell, only: cut_cell_at, small_cell_at
      use cutcell, only: xg_e, yg_n, zg_t
! Gas phase species mass fractions.
      use fldvar, only: t_g, x_g, ep_g
! Gas phase density (variable).
      use fldvar, only: ro_g
!
      use indices, only: i_of, j_of, k_of

      use physprop, only: mw_g, c_pg, nmax
! Simulation time
      use run, only: time
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Fluid cell and solids phase indices
      INTEGER, intent(in) :: IJK
! Flag to output header
      LOGICAL, intent(inout) :: tHeader

! Local Variables:
!---------------------------------------------------------------------//
! Local aliase for inert species index
      INTEGER :: NN
! Local file values.
      LOGICAL :: lExists
      CHARACTER(LEN=255) :: lFile
      INTEGER, parameter :: lUnit = 4868
      LOGICAL, save :: fHeader = .true.
      CHARACTER(LEN=7) :: X_gN
!......................................................................!

      lfile = '';
      if(numpes > 1) then
         write(lfile,"('CPgErr_',I4.4,'.log')") mype
      else
         write(lfile,"('CPgErr.log')")
      endif
      inquire(file=trim(lfile),exist=lexists)
      if(lexists) then
         open(lunit,file=trim(adjustl(lfile)), &
            status='old', position='append',convert='big_endian')
      else
         open(lunit,file=trim(adjustl(lfile)), status='new',&
              convert='big_endian')
      endif

      if(fheader) then
         write(lunit,1000)
         fheader = .false.
      endif

      if(theader) then
         write(lunit,"(/2x,'Simulation time: ',g12.5)") time
         theader = .false.
      endif

      write(lunit,1001) ijk, i_of(ijk), j_of(ijk), k_of(ijk)
      write(lunit,"(6x,A,1X,g12.5)",advance='no') 'C_PG:', c_pg(ijk)
      write(lunit,"(2x,A,1X,g12.5)",advance='no') 'EP_G:', ep_g(ijk)
      write(lunit,"(2x,A,1X,g12.5)") 'T_G:', t_g(ijk)
      write(lunit,"(6x,A,1X)",advance='no') 'X_gN:'
      DO nn = 1,nmax(0)
         write(lunit,"(g12.5,2X)",advance='no') x_g(ijk,nn)
      ENDDO
      write(lunit,fmt='(/)')

      if(cartesian_grid) then
         write(lunit,"(6x,A,1X,L1)",advance='NO') 'Cut Cell:', cut_cell_at(ijk)
         write(lunit,"(6x,A,1X,L1)") 'Small Cell:', small_cell_at(ijk)
         write(lunit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
            xg_e(i_of(ijk)), yg_n(j_of(ijk)), zg_t(k_of(ijk))
      endif

      close(lunit)

      RETURN

 1000 FORMAT(2x,'One or more cells have reported a negative or zero',&
         ' gas specific',/2x,'heat (C_pg(IJK)).')

 1001 FORMAT( 4x,'IJK: ',i8,7x,'I: ',i4,'  J: ',i4,'  K: ',i4)

      END SUBROUTINE cpgerr_log

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: CPsErr_LOG                                              !
!  Purpose: Record information about the location and conditions that  !
!           resulted in a zero or negative solids specific heat.       !
!                                                                      !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE cpserr_log(IJK, M, tHeader)

! Global variables:
!---------------------------------------------------------------------//
      use compar, only: mype, numpes
! Cutcell data.
      use cutcell, only: cartesian_grid
      use cutcell, only: cut_cell_at, small_cell_at
      use cutcell, only: xg_e, yg_n, zg_t
! Solid phase species mass fractions.
      use fldvar, only: t_s, x_s, ep_s
! Solid phase density (variable).
      use fldvar, only: ro_s
!
      use indices, only: i_of, j_of, k_of

      use physprop, only: mw_s, c_ps, nmax
! Simulation time
      use run, only: time
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Fluid cell and solids phase indices
      INTEGER, intent(in) :: IJK, M
! Flag to output header
      LOGICAL, intent(inout) :: tHeader

! Local Variables:
!---------------------------------------------------------------------//
! Local aliase for inert species index
      INTEGER :: NN
! Local file values.
      LOGICAL :: lExists
      CHARACTER(LEN=255) :: lFile
      INTEGER, parameter :: lUnit = 4868
      LOGICAL, save :: fHeader = .true.
      CHARACTER(LEN=7) :: X_sN
!......................................................................!

      lfile = '';
      if(numpes > 1) then
         write(lfile,"('CPsErr_',I4.4,'.log')") mype
      else
         write(lfile,"('CPsErr.log')")
      endif
      inquire(file=trim(lfile),exist=lexists)
      if(lexists) then
         open(lunit,file=trim(adjustl(lfile)), &
            status='old', position='append',convert='big_endian')
      else
         open(lunit,file=trim(adjustl(lfile)), status='new',&
              convert='big_endian')
      endif

      if(fheader) then
         write(lunit,1000)
         fheader = .false.
      endif

      if(theader) then
         write(lunit,"(/2x,'Simulation time: ',g12.5,'  Phase: ',I2)")&
             time, m
         theader = .false.
      endif

      write(lunit,1001) ijk, i_of(ijk), j_of(ijk), k_of(ijk)
      write(lunit,"(6x,A,1X,g12.5)",advance='no') 'C_PS:', c_ps(ijk,m)
      write(lunit,"(2x,A,1X,g12.5)",advance='no') 'EP_S:', ep_s(ijk,m)
      write(lunit,"(2x,A,1X,g12.5)") 'T_S:', t_s(ijk,m)
      write(lunit,"(6x,A,1X)",advance='no') 'X_sN:'
      DO nn = 1,nmax(m)
         write(lunit,"(g12.5,2X)",advance='no') x_s(ijk,m,nn)
      ENDDO
      write(lunit,fmt='(/)')

      if(cartesian_grid) then
         write(lunit,"(6x,A,1X,L1)",advance='NO') 'Cut Cell:', cut_cell_at(ijk)
         write(lunit,"(6x,A,1X,L1)") 'Small Cell:', small_cell_at(ijk)
         write(lunit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
            xg_e(i_of(ijk)), yg_n(j_of(ijk)), zg_t(k_of(ijk))
      endif

      close(lunit)

      RETURN

 1000 FORMAT(2x,'One or more cells have reported a negative or zero',&
         ' solids specific',/2x,'heat (C_ps(IJK)).')

 1001 FORMAT( 4x,'IJK: ',i8,7x,'I: ',i4,'  J: ',i4,'  K: ',i4)

      END SUBROUTINE cpserr_log