read_thermochemical_mod.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Module name: read_thermochemical_mod.f                              C
!  Purpose: Read thermochemical data                                   C
!                                                                      C
!  Revision Number:                                                    C
!  Purpose:                                                            C
!  Author:                                            Date: dd-mmm-yy  C
!  Reviewer:                                          Date: dd-mmm-yy  C
!                                                                      C
!  Literature/Document References: None                                C
!                                                                      C
!  Variables referenced: None                                          C
!  Variables modified: None                                            C
!                                                                      C
!  Local variables: None                                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

#include "version.inc"

MODULE read_thermochemical

  IMPLICIT NONE

! Full path to Burcat and Ruscic database
  CHARACTER(len=142) :: therm = 'BURCAT.THR'

CONTAINS

!      Program Test; CALL Read_Therm_tester; END Program Test

      SUBROUTINE read_therm_tester
      Implicit none
      DOUBLE PRECISION Ahigh(7), Alow(7)
      DOUBLE PRECISION Thigh, Tlow, Tcom, MW
      DOUBLE PRECISION Cp1, h1, h2, Hf298oR
      CHARACTER(LEN=132) :: PATH
      CHARACTER(LEN=18) :: SPECIES
      integer funit, IER
      CHARACTER(len=255) FILENAME
      LOGICAL LocalCopy

      species = 'CH4'
      path = '.'
      funit = 5

      INQUIRE(file=trim(therm),exist=localcopy)
      IF(localcopy)Then
        OPEN(convert='BIG_ENDIAN',unit=funit,file=trim(therm))
      ELSE
        filename = './BURCAT.THR'
        OPEN(convert='BIG_ENDIAN',unit=funit,file=trim(filename), err=500)
      ENDIF

 !      Call Read_Therm(PATH, 'N2', Thigh, Tlow, Tcom, Ahigh, Alow, Hf298oR)
      Call read_therm(funit, species, thigh, tlow, tcom, mw, ahigh, &
         alow, hf298or, ier)
      IF(ier /= 0) GOTO 200

      print *, species
      print *, thigh, tlow, tcom, mw, hf298or*1.987207

!      print *, Hf298oR
!      T = 300
!      DO i = 1, 12
!        Cp1 = calc_CpoR(T, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
!        T = T + 100
!        print *, T, Cp1
!      ENDDO

!      Cp1 = calc_CpoR(8D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
!      h1 = calc_H0oR(4D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
!      h2 = calc_H0oR(12D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
      print *, cp1, h1, h2
      CLOSE(unit=funit)
      error_stop
200   print *, 'READ_Therm_tester: Species ', &
         trim(species), ' not found in Database!'
      error_stop
500   print *, 'READ_Therm_tester: Cannot Open file ', trim(therm), '!'
      print *, 'Check path or copy mfix/model/thermochemical/', &
         trim(therm), ' into run directory'
      error_stop
      END Subroutine read_therm_tester

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv  C
!                                                                        C
!     Module name: READ_Therm()                                          C
!     Purpose: Read Thermo coefficients from Burcat and Ruscic (2005)    C
!     Author: M. Syamlal                                 Date: 30-SEP-05 C
!                                                                        C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  C
!
      SUBROUTINE read_therm(funit, Sp, Thigh, Tlow, Tcom, MW, Ahigh, &
         alow, hf298or, ier)
!
!-----------------------------------------------
!     M o d u l e s
!-----------------------------------------------

      IMPLICIT NONE

      CHARACTER(*) SP

!     holds one line in the input file
      CHARACTER(len=80) :: LINE_STRING

      CHARACTER(len=18) :: SPECIES, ss
      INTEGER i, funit, IER
      DOUBLE PRECISION Ahigh(7), Alow(7), Hf298oR
      DOUBLE PRECISION Thigh, Tlow, Tcom, MW


!     Associate a simple species name with that in BURCAT.THR file
      INTEGER, PARAMETER :: Max = 3
      CHARACTER(LEN=18), DIMENSION(2,Max) :: SPECIES_ALIAS
      species_alias = reshape( (/ &
!        common name             BURCAT.THR name
!        123456789012345678     123456789012345678
        'O2                ',  'O2 REF ELEMENT    ', &
        'N2                ',  'N2  REF ELEMENT   ', &
        'CH4               ',  'CH4   ANHARMONIC  ' &
       /), (/2,max/))


      ier = 0
      species = sp
      DO i = 1, max
        IF(trim(sp) == trim(species_alias(1,i)))THEN
          species = species_alias(2,i)
        ENDIF
      ENDDO


      line_string = '                '
      DO WHILE(line_string(1:11) /= 'THERMO DATA')
        READ(unit=funit, fmt='(A)',err=100,end=100)line_string
      END DO

      ss = '                 '
      call trimtab(species)
      DO WHILE(trim(ss) /= trim(species))
        READ(unit=funit, fmt='(A)',err=100,end=100)line_string
        ss = line_string(1:18)
        call trimtab(ss)

      END DO
!      print *, LINE_STRING


      call get_values(line_string, tlow, thigh, mw)

! Tcom is almost always 1000K, however there are a few species where
! this value is too high and causes a problem (e.g., liquid water).
! Therefore, set Tcom = Thigh when Thigh < 1000K.
      tcom = min(1.0d3, thigh)
      READ(unit=funit, fmt='(5E15.0)',err=300,end=300)ahigh(1:5)
      READ(unit=funit, fmt='(5E15.0)',err=300,end=300)ahigh(6:7), alow(1:3)
      READ(unit=funit, fmt='(5E15.0)',err=300,end=300)alow(4:7), hf298or



      RETURN
      ! species not found or THERMO DATA not found!
100   ier = 1
      RETURN

300   print *, 'READ_Therm: Error reading coefficients for Species ', &
         trim(line_string(1:18))
      error_stop

      END SUBROUTINE read_therm

!**********************************************************************!
! Function: calc_CpoR                                                  !
! Purpose: Evaluate the polynomial form of the specific heat.          !
!                                                                      !
!**********************************************************************!
      PURE DOUBLE PRECISION FUNCTION  calc_cpor(T, M, N)

! Polynomial coefficients
      use physprop, only: ahigh  ! for T in [Tcom, Thigh]
      use physprop, only: alow   ! for T in [Tlow, Tcom)
      use physprop, only: thigh  ! Upper bound of use
      use physprop, only: tlow   ! Lower bound of use
      use physprop, only: tcom   ! Switch from low to high coeffs

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Evaluation temperaure (K)
      DOUBLE PRECISION, intent(in) :: T
! Phase index.
      INTEGER, intent(in) :: M
! Species index.
      INTEGER, intent(in) :: N

! Local Variables:
!-----------------------------------------------------------------------
! Bounded temperature.
      DOUBLE PRECISION :: xT

! Bound the temperature to the valid region.
      xt = max(min(t,thigh(m,n)),tlow(m,n))

! Evaluate the polynomial form.
      IF(t < tcom(m,n))THEN
        calc_cpor = calc_cpor0(xt, alow(1:5,m,n))
      ELSE
        calc_cpor = calc_cpor0(xt, ahigh(1:5,m,n))
      ENDIF

      RETURN
      END Function calc_cpor


!**********************************************************************!
! Function: calc_CpoR0                                                 !
! Purpose: Evaluate the polynomial form of the specific heat.          !
!                                                                      !
!**********************************************************************!
      PURE DOUBLE PRECISION FUNCTION  calc_cpor0(T, A)

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Evaluation temperaure (K)
      DOUBLE PRECISION, intent(in) :: T
! Polynomial coefficients.
      DOUBLE PRECISION, intent(in) :: A(1:5)

! Evaluate the polynomial.
      calc_cpor0 = (((a(5)*t +a(4))*t + a(3))*t + a(2))*t + a(1)

      RETURN
      END Function calc_cpor0


!**********************************************************************!
! Function: calc_ICpoR                                                 !
! Purpose: Integrate the polynomial form of the specific heat.         !
!                                                                      !
!**********************************************************************!
      DOUBLE PRECISION FUNCTION calc_icpor(T, M, N, IER)

      use physprop, only: ahigh
      use physprop, only: thigh
      use physprop, only: icpor_h
      use physprop, only: alow
      use physprop, only: tlow
      use physprop, only: icpor_l
      use physprop, only: tcom

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Evaluation temperaure (K)
      DOUBLE PRECISION, intent(in) :: T
! Phase index.
      INTEGER, intent(in) :: M
! Species index.
      INTEGER, intent(in) :: N
! Error Flag.
      INTEGER, intent(inout) :: IER

! Local Variables.
!-----------------------------------------------------------------------
      DOUBLE PRECISION :: xT
!-----------------------------------------------------------------------

! Initialize the bounded temperature and error flag.
      xt = t
      ier = 0

! Verify that the temperature is in a valid range.
      if(t > thigh(m,n)) THEN
        xt = thigh(m,n)
        ier = 101
      elseif(t < tlow(m,n)) THEN
        xt = tlow(m,n)
        ier = 102
      endif

! Integrate the polynomial from 0.0 to T.
      if (xt < tcom(m,n)) then
        calc_icpor = calc_icpor0(xt, alow(1:5,m,n),  icpor_l(m,n))
      else
        calc_icpor = calc_icpor0(xt, ahigh(1:5,m,n), icpor_h(m,n))
      endif

      RETURN
      END FUNCTION calc_icpor


!**********************************************************************!
! Function: calc_ICpoR                                                 !
! Purpose: Integrate the polynomial form of the specific heat.         !
!                                                                      !
!**********************************************************************!
      PURE DOUBLE PRECISION FUNCTION calc_icpor0(T, A, REF_ICpoR)

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Evaluation temperaure (K)
      DOUBLE PRECISION, intent(in) :: T
! Polynomial coefficients.
      DOUBLE PRECISION, intent(in) :: A(1:5)
! Referece Integral
      DOUBLE PRECISION, intent(in) :: REF_ICpoR

! Local Variables.
!-----------------------------------------------------------------------
! Integral of specific heat polynomial (from 0 to T) over T
      DOUBLE PRECISION ICpoRT

!-----------------------------------------------------------------------

      icport = (((a(5)*t/5.0d0 + a(4)/4.0d0)*t + a(3)/3.0d0)*t +       &
         a(2)/2.0d0)*t + a(1)

      calc_icpor0 = t*icport - ref_icpor

      RETURN
      END FUNCTION calc_icpor0



!**********************************************************************!
! Function: calc_H0oR                                                  !
! Purpose: Calculate the heat of formation from the first six poly-    !
!          nomial coefficients.                                        !
!                                                                      !
! >>> This function is currently unused.                               !
!                                                                      !
!**********************************************************************!
      PURE DOUBLE PRECISION FUNCTION calc_h0or(T, Th, Tl, Tc, Ah, Al)

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Polynomial coefficients. (High/Low)
      DOUBLE PRECISION, intent(in) :: Ah(7), Al(7)
! Temperature ranges of polynomials.
      DOUBLE PRECISION, intent(in) :: Th   ! Max temp (for Ahigh)
      DOUBLE PRECISION, intent(in) :: Tl   ! Min temp (for Alow)
      DOUBLE PRECISION, intent(in) :: Tc   ! switch from low to high
! Evaluation temperaure (K)
      DOUBLE PRECISION, intent(in) :: T

! Local Variables.
!-----------------------------------------------------------------------

      !ICp = calc_ICpoR(T, Th, Tl, Tc, Ah, Al)
      !If (T < Tc) then
      !  calc_H0oR = ICp + Al(6)
      !else
      !  calc_H0oR = ICp + Ah(6)
      !endif

      calc_h0or = 0.0

      return
      END FUNCTION calc_h0or


!**********************************************************************!
! SUBROUTINE: replaceTab                                               !
! Purpose: Replace all instances of a tab with a single space.         !
!**********************************************************************!
      SUBROUTINE replacetab(C)

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Incoming string that will have tabs removed.
      CHARACTER(len=*) :: C

! Local Variables:
!-----------------------------------------------------------------------
! Loop counter
      INTEGER :: I

      DO i = 1, len(c)
        IF(c(i:i) == char(9)) c(i:i)=' '
      ENDDO

      RETURN
      END SUBROUTINE replacetab


!**********************************************************************!
! SUBROUTINE: trimTab                                                  !
! Purpose: Search a string for the first instance of a tab. The        !
!          location of the tab and all remaing string entries are      !
!          replaced with blank spaces.                                 !
!**********************************************************************!
      SUBROUTINE trimtab(C)

      implicit none

! Dummy Arguments:
!-----------------------------------------------------------------------
! Incoming string that will have tabs removed.
      CHARACTER(len=*) :: C

! Local Variables:
!-----------------------------------------------------------------------
! Loop counter
      INTEGER :: I
! Logical indicating that a tab was located.
      LOGICAL :: tabFound

! Initialize flag
      tabfound = .false.

! Look at each entry of the string. Once a tab is located, the rest of
! the string is replaced by blank spaces.
      DO i = 1, len(c)
        IF(c(i:i) == char(9) ) tabfound = .true.
        if(tabfound) c(i:i)=' '
      ENDDO

      RETURN
      END SUBROUTINE trimtab

    END MODULE read_thermochemical