read_database.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: read_database(Ier)                                     C
!  Purpose: read thermochemical database                               C
!                                                                      C
!  Author: M. Syamlal                                 Date: 21-OCT-05  C
!                                                                      C
!  Modification 1: J. Musser                          Date: 02-May-11  C
!  Purpose: Provided support for DEM access to database.               C
!                                                                      C
!  Modification 2: J. Musser                          Date: 02-Oct-12  C
!  Purpose: Calls to READ_DATABASE were moved to check_gas_phase and   C
!  check_solids_common_all during input data integrity checks for the  C
!  gas and solids phases.                                              C
!  Rather than looping through all species for each phase, the model   C
!  (TFM/DEM), phase index, species index, species name, and molecular  C
!  weight are passed as dummy arguments so that only information for   C
!  referenced species (lName) is obtained.                             C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE read_database(lM, lN, lName, lMW)

      USE param
      USE param1
      USE physprop
      USE constant
      USE compar
      USE rxns
      USE funits
      USE discretelement
      USE des_rxns
      USE read_thermochemical, only: read_therm, calc_icpor, therm
      use run, only: reinitializing
      use error_manager

      IMPLICIT NONE

! Phase and species indices
      INTEGER, INTENT(IN) :: lM, lN
! Species name from data file.
      CHARACTER(len=*), INTENT(IN) :: lName
! Species molecular weight from the data file (if any)
      DOUBLE PRECISION, INTENT(INOUT) :: lMW

! Molecular weight read from database.
      DOUBLE PRECISION dbMW

! Error message returned from Read_Therm and sent to calling routine
      INTEGER IER
! File unit of databases searched
      INTEGER FUNIT
! Loop counter for checking all three locations for data
      INTEGER FILE
! Input/Output error status ( 0 is no error)
      INTEGER IOS

! Identifies if an error was detected.
      LOGICAL ErrorFlag

! Tcom +/- SMALL_NUMBER: This is done so that the specific heat
! polynomial can be evaluated at Tcom with the high and low
! coefficients.
      DOUBLE PRECISION :: xTc

! Various integrations of the specific heat polynomials:
      DOUBLE PRECISION :: ICpoR_TrL  ! 0.0 --> Tref using Alow
      DOUBLE PRECISION :: ICpoR_TcL  ! 0.0 --> Tcom using Alow
      DOUBLE PRECISION :: ICpoR_TcH  ! 0.0 --> Tcom using Ahigh

      LOGICAL :: testCp = .false.
! Database being searched.
      CHARACTER(len=256) :: DB

#ifdef BURCAT_THR
      therm = burcat_thr
#endif

      CALL init_err_msg('READ_DATABASE')

! Initialize the file unit to be used.
      funit = unit_dat  ! .dat file unit
! Read data from mfix.dat or from BURCAT.THR in run directory or
! mfix directory.
      file = 0
      db_lp: DO
         file = file + 1
! Check for thermochemical data in the mfix.dat file.
         IF(file == 1) THEN
           OPEN(convert='BIG_ENDIAN',unit=funit, file='mfix.dat', status='OLD', iostat= ios)
           IF(ios /= 0) cycle db_lp
           db=''; WRITE(db,1000) 'mfix.dat'
! Read thermochemical data from the BURCAT.THR database in the local
! run directory.
         ELSEIF(file == 2) THEN
            OPEN(convert='BIG_ENDIAN',unit=funit,file=trim(therm), status='OLD', iostat= ios)
            IF(ios /= 0) cycle db_lp
            db=''; WRITE(db,1000) trim(therm)
          ELSE
            EXIT db_lp
          ENDIF

         rewind(unit=funit)

! Initialize the error flag
         ier = 0

         CALL read_therm(funit, lname, thigh(lm,ln), tlow(lm,ln),      &
            tcom(lm,ln), dbmw, ahigh(:,lm,ln), alow(:,lm,ln),          &
            hfrefor(lm,ln), ier)

         IF(ier == 0) THEN
! If the user did not supply a value for the gas phase molecular weight
! in the mfix.dat file, use the value from the database.
            IF(lmw == undefined) lmw = dbmw
! There are a number of species with Tlow as 300, for which the
! following calculation will produce an error because T_ref = 298.  So
! slightly extend validity of the correlation.
            IF(abs(tlow(lm,ln)-t_ref)<=2.0d0 .AND. &
               tlow(lm,ln) > t_ref) tlow(lm,ln) = t_ref

! Initialize the reference integrals.
            icpor_l(lm,ln) = zero
            icpor_h(lm,ln) = zero

! Calculate the integral of specific heat from zero to Tref using the
! Alow coefficients.
               icpor_trl = calc_icpor(t_ref, lm, ln, ier)
! Calculate the integral of specific heat from zero to Tcom using the
! Alow coefficients.
            xtc = tcom(lm,ln)-small_number
            icpor_tcl = calc_icpor(xtc, lm, ln, ier)
! Calculate the integral of specific heat from zero to Tcom using the
! Ahigh coefficients.
            xtc = tcom(lm,ln)+small_number
            icpor_tch = calc_icpor(xtc, lm, ln, ier)
! Store the integrals in global variables.
            icpor_l(lm,ln) = icpor_trl
            icpor_h(lm,ln) = icpor_tch - (icpor_tcl - icpor_trl)
         ENDIF

         errorflag = .true.
         IF(ier == 0) THEN
            IF(.NOT.reinitializing)THEN
               WRITE(err_msg,1001) trim(adjustl(db)), 'Found!'
               CALL flush_err_msg(header=.false., footer=.false.)
            ENDIF
            if(testcp) CALL writecp(lm, ln, lname, lmw)
            errorflag = .false.
            EXIT db_lp
        ELSEIF(.NOT.reinitializing) THEN
            WRITE(err_msg,1001) trim(adjustl(db)), 'Not found.'
            CALL flush_err_msg(header=.false., footer=.false.)
        ENDIF

      ENDDO db_lp

! Error message control.
!-----------------------------------------------------------------------
! Write remaining error message if needed.
      IF(errorflag) THEN
         CALL flush_err_msg(header=.false.)
         WRITE(err_msg,1010) trim(lname), trim(therm)
         CALL flush_err_msg(abort=.true.)
      ENDIF

      CALL finl_err_msg

      RETURN


! Messages
!-----------------------------------------------------------------------
 1000 FORMAT('Checking ',a)
 1001 FORMAT(8x,a,1x,' :: ',a)

! Error Flags
!-----------------------------------------------------------------------
 1010 FORMAT('Message 1010: Species "',a,'" was not matched to any ',  &
         'entry in the',/'thermochemical databases.',2/,'SUGGESTION: ',&
         'Search the database for the exact species name. The ',/      &
         'species names are case sensitive and should match the names',&
         ' in',/'BURCAT.THR exactly excluding trailing blanks and ',   &
         'tabs. Also verify',/'that the data section in the mfix.dat ',&
         'file (if any) is below a line',/'that starts with THERMO ',  &
         'DATA.',2/'Database location:', /a)

      END SUBROUTINE read_database

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: READ_DATABASE0(IER)                                    C
!  Purpose: Provides legacy support for rrates files.                  C
!                                                                      C
!  Author: J. Musser                                  Date: 02-Oct-12  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE read_database0()

      USE compar
      USE constant
      USE des_rxns
      USE discretelement
      USE exit, only: mfix_exit
      USE funits
      USE param
      USE param1
      USE physprop
      USE rxns

      IMPLICIT NONE

! Loop indices for mass phase and species
      INTEGER M, NN
! Loop counter for continuum and discrete species
      INTEGER Nsp, DES_Nsp

! Return to the calling routine if the database has already been called.
      IF(database_read)RETURN

! Set the flag identifying that the database has been read.
      database_read = .true.

! Initialize counters
      nsp = 0
      des_nsp = 0

! Read species data for the gas phase.
!-----------------------------------------------------------------------
      DO nn = 1, nmax(0)
         nsp = nsp + 1
! If a species name was not specified in mfix.dat, flag error and exit.
          IF(species_name(nsp) == undefined_c) THEN
            WRITE(*,1010) nn         ! screen
            IF(dmp_log) WRITE(unit_log,1010) nn  ! log file
             CALL mfix_exit(mype)
          ENDIF
! Read the database.
         CALL read_database(0, nn, species_name(nsp), mw_g(nn))
       ENDDO

! Read species data for the continuum solids phases.
!-----------------------------------------------------------------------
! Skip reading the database for the continuum solids phase if the
! simulation is only employing discrete solids.
      IF(.NOT.discrete_element .OR. des_continuum_hybrid)THEN
          DO m = 1, mmax
            DO nn = 1, nmax(m)
                nsp = nsp + 1
! If a species name was not specified in mfix.dat, flag error and exit.
                IF(species_name(nsp) == undefined_c)THEN
                  WRITE(*,1011)'continuum', m, nn ! screen
                  IF(dmp_log) WRITE(unit_log,1011)'continuum', m, nn
                   CALL mfix_exit(mype)
                ENDIF
               CALL read_database(m, nn, species_name(nsp), mw_s(m,nn))
             ENDDO   ! N=1, NMAX(M)
          ENDDO   ! M=1, MMAX
      ENDIF

      RETURN

! Error Messages
!-----------------------------------------------------------------------
 1010 FORMAT(/1x,70('*')/, ' From: READ_DATABASE0',/, ' Message: ',    &
         'No SPECIES_NAME provided for gas phase species ',i3,'.',/' ',&
         'Check mfix.dat.',/1x,70('*')/)
 1011 FORMAT(/1x,70('*')/, ' From: READ_DATABASE0',/, ' Message: ',    &
         'No SPECIES_NAME provided for ',a,' solids phase ',i2,', ',/  &
         ' species ',i3,'.',/' Check mfix.dat.',/1x,70('*')/)

      END SUBROUTINE read_database0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: read_database(Ier)                                     C
!  Purpose: read thermochemical database                               C
!                                                                      C
!  Author: M. Syamlal                                 Date: 21-OCT-05  C
!                                                                      C
!  Modification 1: J. Musser                          Date: 02-May-11  C
!  Purpose: Provided support for DEM access to database.               C
!                                                                      C
!  Modification 2: J. Musser                          Date: 02-Oct-12  C
!  Purpose: Calls to READ_DATABASE were moved to CHECK_DATA_04/05      C
!  duing input data integrity checks for the gas and solids phases.    C
!  Rather than looping through all species for each phase, the model   C
!  (TFM/DEM), phase index, species index, species name, and molecular  C
!  weight are passed as dummy arguments so that only infomration for   C
!  referenced species (lName) is obtained.                             C                             C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE writecp(lM, lN, lName, lMW)

      use param1
      USE physprop
      USE compar
      USE read_thermochemical, only: calc_cpor, calc_icpor, calc_icpor0

! Universal gas constant in cal/mol.K
      use constant, only: rgas => gas_const_cal

      IMPLICIT NONE

! Phase and species indices
      INTEGER, INTENT(IN) :: lM, lN
! Species name from data file.
      CHARACTER(len=*), INTENT(IN) :: lName
! Species molecular weight from the data file (if any)
      DOUBLE PRECISION, INTENT(in) :: lMW

      INTEGER :: IER1, IER2, lc

      DOUBLE PRECISION :: T
      DOUBLE PRECISION :: lCP, lICP

      write(*,"(2/,3x,'Specific Heat report for ',A)")trim(lname)

      write(*,"(/12x,'Low',9x,'High')")
      write(*,"(6x,'T',3x,g12.5,2x,g12.5)") tlow(lm,ln), thigh(lm,ln)
      DO lc=1,5
         write(*,"(4x,'A(',I1,')',2(2x,g12.5))") lc, &
            alow(lc,lm,ln), ahigh(lc,lm,ln)
      ENDDO
      write(*,"('')")
      write(*,"(5x,'Tcom: ',g12.5)")tcom(lm,ln)
      write(*,"('')")

      write(*,"(5x,'Temperature',8x,'Cp',11x,'ICp')")

      t = tcom(lm,ln) - 100.0
      DO WHILE(t <= tcom(lm,ln) - small_number)

         ier1 = 0
         ier2 = 0

         write(*,"(7x,g12.5)",advance="NO") t
         lcp  = calc_cpor(t, lm, ln) * rgas / lmw
         licp = calc_icpor(t, lm, ln, ier2) * rgas / lmw
         write(*,"(2(3x,g12.5))",advance="NO")lcp, licp

         IF(ier1 /= 0) write(*,"(3x,'Cp Error!')",advance="NO")
         IF(ier2 /= 0) write(*,"(3x,'ICp Error!')",advance="NO")
         write(*,"('')")

         t = t + 5.0
      ENDDO


      t = tcom(lm,ln) + small_number
      DO WHILE(t <= tcom(lm,ln) + 100.0)

         ier1 = 0
         ier2 = 0

         write(*,"(7x,g12.5)",advance="NO") t
         lcp  = calc_cpor(t, lm, ln) * rgas / lmw
         licp = calc_icpor(t, lm, ln, ier2) * rgas / lmw
         write(*,"(2(3x,g12.5))",advance="NO")lcp, licp

         IF(ier1 /= 0) write(*,"(3x,'Cp Error!')",advance="NO")
         IF(ier2 /= 0) write(*,"(3x,'ICp Error!')",advance="NO")
         write(*,"('')")

         t = t + 5.0
      ENDDO

      write(*,"('')")
      END SUBROUTINE writecp