File: N:\mfix\model\thermochemical\read_thermochemical_mod.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: read_thermochemical_mod.f                              C
4     !  Purpose: Read thermochemical data                                   C
5     !                                                                      C
6     !  Revision Number:                                                    C
7     !  Purpose:                                                            C
8     !  Author:                                            Date: dd-mmm-yy  C
9     !  Reviewer:                                          Date: dd-mmm-yy  C
10     !                                                                      C
11     !  Literature/Document References: None                                C
12     !                                                                      C
13     !  Variables referenced: None                                          C
14     !  Variables modified: None                                            C
15     !                                                                      C
16     !  Local variables: None                                               C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19     
20     #include "version.inc"
21     
22     MODULE read_thermochemical
23     
24       IMPLICIT NONE
25     
26     ! Full path to Burcat and Ruscic database
27       CHARACTER(len=142) :: THERM = 'BURCAT.THR'
28     
29     CONTAINS
30     
31     !      Program Test; CALL Read_Therm_tester; END Program Test
32     
33           SUBROUTINE READ_Therm_tester
34           Implicit none
35           DOUBLE PRECISION Ahigh(7), Alow(7)
36           DOUBLE PRECISION Thigh, Tlow, Tcom, MW
37           DOUBLE PRECISION Cp1, h1, h2, Hf298oR
38           CHARACTER(LEN=132) :: PATH
39           CHARACTER(LEN=18) :: SPECIES
40           integer funit, IER
41           CHARACTER(len=255) FILENAME
42           LOGICAL LocalCopy
43     
44           SPECIES = 'CH4'
45           PATH = '.'
46           funit = 5
47     
48           INQUIRE(FILE=TRIM(THERM),EXIST=LocalCopy)
49           IF(LocalCopy)Then
50             OPEN(CONVERT='BIG_ENDIAN',UNIT=funit,FILE=TRIM(THERM))
51           ELSE
52             FILENAME = './BURCAT.THR'
53             OPEN(CONVERT='BIG_ENDIAN',UNIT=funit,FILE=TRIM(FILENAME), ERR=500)
54           ENDIF
55     
56      !      Call Read_Therm(PATH, 'N2', Thigh, Tlow, Tcom, Ahigh, Alow, Hf298oR)
57           Call Read_Therm(funit, SPECIES, Thigh, Tlow, Tcom, MW, Ahigh, &
58              Alow, Hf298oR, IER)
59           IF(IER /= 0) GOTO 200
60     
61           print *, SPECIES
62           print *, Thigh, Tlow, Tcom, MW, Hf298oR*1.987207
63     
64     !      print *, Hf298oR
65     !      T = 300
66     !      DO i = 1, 12
67     !        Cp1 = calc_CpoR(T, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
68     !        T = T + 100
69     !        print *, T, Cp1
70     !      ENDDO
71     
72     !      Cp1 = calc_CpoR(8D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
73     !      h1 = calc_H0oR(4D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
74     !      h2 = calc_H0oR(12D2, Thigh, Tlow, Tcom, Ahigh, Alow)*1.987207
75           print *, Cp1, h1, h2
76           CLOSE(UNIT=funit)
77           ERROR_STOP
78     200   PRINT *, 'READ_Therm_tester: Species ', &
79              TRIM(SPECIES), ' not found in Database!'
80           ERROR_STOP
81     500   PRINT *, 'READ_Therm_tester: Cannot Open file ', TRIM(THERM), '!'
82           PRINT *, 'Check path or copy mfix/model/thermochemical/', &
83              TRIM(THERM), ' into run directory'
84           ERROR_STOP
85           END Subroutine READ_Therm_tester
86     
87     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv  C
88     !                                                                        C
89     !     Module name: READ_Therm()                                          C
90     !     Purpose: Read Thermo coefficients from Burcat and Ruscic (2005)    C
91     !     Author: M. Syamlal                                 Date: 30-SEP-05 C
92     !                                                                        C
93     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  C
94     !
95           SUBROUTINE READ_Therm(funit, Sp, Thigh, Tlow, Tcom, MW, Ahigh, &
96              Alow, Hf298oR, IER)
97     !
98     !-----------------------------------------------
99     !     M o d u l e s
100     !-----------------------------------------------
101     
102           IMPLICIT NONE
103     
104           CHARACTER(*) SP
105     
106     !     holds one line in the input file
107           CHARACTER(len=80) :: LINE_STRING
108     
109           CHARACTER(len=18) :: SPECIES, ss
110           INTEGER i, funit, IER
111           DOUBLE PRECISION Ahigh(7), Alow(7), Hf298oR
112           DOUBLE PRECISION Thigh, Tlow, Tcom, MW
113     
114     
115     !     Associate a simple species name with that in BURCAT.THR file
116           INTEGER, PARAMETER :: Max = 3
117           CHARACTER(LEN=18), DIMENSION(2,Max) :: SPECIES_ALIAS
118           SPECIES_ALIAS = RESHAPE ( (/ &
119     !        common name             BURCAT.THR name
120     !        123456789012345678     123456789012345678
121             'O2                ',  'O2 REF ELEMENT    ', &
122             'N2                ',  'N2  REF ELEMENT   ', &
123             'CH4               ',  'CH4   ANHARMONIC  ' &
124            /), (/2,Max/))
125     
126     
127           IER = 0
128           SPECIES = SP
129           DO i = 1, MAX
130             IF(TRIM(SP) == TRIM(SPECIES_ALIAS(1,I)))THEN
131               SPECIES = SPECIES_ALIAS(2,I)
132             ENDIF
133           ENDDO
134     
135     
136           LINE_STRING = '                '
137           DO WHILE(LINE_STRING(1:11) /= 'THERMO DATA')
138             READ(UNIT=funit, FMT='(A)',ERR=100,END=100)LINE_STRING
139           END DO
140     
141           ss = '                 '
142           call trimTab(SPECIES)
143           DO WHILE(TRIM(ss) /= TRIM(SPECIES))
144             READ(UNIT=funit, FMT='(A)',ERR=100,END=100)LINE_STRING
145             ss = LINE_STRING(1:18)
146             call trimTab(ss)
147     
148           END DO
149     !      print *, LINE_STRING
150     
151     
152           call get_values(LINE_STRING, Tlow, Thigh, MW)
153     
154     ! Tcom is almost always 1000K, however there are a few species where
155     ! this value is too high and causes a problem (e.g., liquid water).
156     ! Therefore, set Tcom = Thigh when Thigh < 1000K.
157           Tcom = min(1.0d3, Thigh)
158           READ(UNIT=funit, FMT='(5E15.0)',ERR=300,END=300)Ahigh(1:5)
159           READ(UNIT=funit, FMT='(5E15.0)',ERR=300,END=300)Ahigh(6:7), Alow(1:3)
160           READ(UNIT=funit, FMT='(5E15.0)',ERR=300,END=300)Alow(4:7), Hf298oR
161     
162     
163     
164           RETURN
165           ! species not found or THERMO DATA not found!
166     100   IER = 1
167           RETURN
168     
169     300   PRINT *, 'READ_Therm: Error reading coefficients for Species ', &
170              TRIM(LINE_STRING(1:18))
171           ERROR_STOP
172     
173           END SUBROUTINE READ_Therm
174     
175     !**********************************************************************!
176     ! Function: calc_CpoR                                                  !
177     ! Purpose: Evaluate the polynomial form of the specific heat.          !
178     !                                                                      !
179     !**********************************************************************!
180           PURE DOUBLE PRECISION FUNCTION  calc_CpoR(T, M, N)
181     
182     ! Polynomial coefficients
183           use physprop, only: Ahigh  ! for T in [Tcom, Thigh]
184           use physprop, only: Alow   ! for T in [Tlow, Tcom)
185           use physprop, only: Thigh  ! Upper bound of use
186           use physprop, only: Tlow   ! Lower bound of use
187           use physprop, only: Tcom   ! Switch from low to high coeffs
188     
189           implicit none
190     
191     ! Dummy Arguments:
192     !-----------------------------------------------------------------------
193     ! Evaluation temperaure (K)
194           DOUBLE PRECISION, intent(in) :: T
195     ! Phase index.
196           INTEGER, intent(in) :: M
197     ! Species index.
198           INTEGER, intent(in) :: N
199     
200     ! Local Variables:
201     !-----------------------------------------------------------------------
202     ! Bounded temperature.
203           DOUBLE PRECISION :: xT
204     
205     ! Bound the temperature to the valid region.
206           xT = max(min(T,Thigh(M,N)),Tlow(M,N))
207     
208     ! Evaluate the polynomial form.
209           IF(T < Tcom(M,N))THEN
210             calc_CpoR = calc_CpoR0(xT, Alow(1:5,M,N))
211           ELSE
212             calc_CpoR = calc_CpoR0(xT, Ahigh(1:5,M,N))
213           ENDIF
214     
215           RETURN
216           END Function calc_CpoR
217     
218     
219     !**********************************************************************!
220     ! Function: calc_CpoR0                                                 !
221     ! Purpose: Evaluate the polynomial form of the specific heat.          !
222     !                                                                      !
223     !**********************************************************************!
224           PURE DOUBLE PRECISION FUNCTION  calc_CpoR0(T, A)
225     
226           implicit none
227     
228     ! Dummy Arguments:
229     !-----------------------------------------------------------------------
230     ! Evaluation temperaure (K)
231           DOUBLE PRECISION, intent(in) :: T
232     ! Polynomial coefficients.
233           DOUBLE PRECISION, intent(in) :: A(1:5)
234     
235     ! Evaluate the polynomial.
236           calc_CpoR0 = (((A(5)*T +A(4))*T + A(3))*T + A(2))*T + A(1)
237     
238           RETURN
239           END Function calc_CpoR0
240     
241     
242     !**********************************************************************!
243     ! Function: calc_ICpoR                                                 !
244     ! Purpose: Integrate the polynomial form of the specific heat.         !
245     !                                                                      !
246     !**********************************************************************!
247           DOUBLE PRECISION FUNCTION calc_ICpoR(T, M, N, IER)
248     
249           use physprop, only: Ahigh
250           use physprop, only: Thigh
251           use physprop, only: ICpoR_h
252           use physprop, only: Alow
253           use physprop, only: Tlow
254           use physprop, only: ICpoR_l
255           use physprop, only: Tcom
256     
257           implicit none
258     
259     ! Dummy Arguments:
260     !-----------------------------------------------------------------------
261     ! Evaluation temperaure (K)
262           DOUBLE PRECISION, intent(in) :: T
263     ! Phase index.
264           INTEGER, intent(in) :: M
265     ! Species index.
266           INTEGER, intent(in) :: N
267     ! Error Flag.
268           INTEGER, intent(inout) :: IER
269     
270     ! Local Variables.
271     !-----------------------------------------------------------------------
272           DOUBLE PRECISION :: xT
273     !-----------------------------------------------------------------------
274     
275     ! Initialize the bounded temperature and error flag.
276           xT = T
277           IER = 0
278     
279     ! Verify that the temperature is in a valid range.
280           if(T > Thigh(M,N)) THEN
281             xT = Thigh(M,N)
282             IER = 101
283           elseif(T < Tlow(M,N)) THEN
284             xT = Tlow(M,N)
285             IER = 102
286           endif
287     
288     ! Integrate the polynomial from 0.0 to T.
289           if (xT < Tcom(M,N)) then
290             calc_ICpoR = calc_ICpoR0(xT, Alow(1:5,M,N),  ICpoR_l(M,N))
291           else
292             calc_ICpoR = calc_ICpoR0(xT, Ahigh(1:5,M,N), ICpoR_h(M,N))
293           endif
294     
295           RETURN
296           END FUNCTION calc_ICpoR
297     
298     
299     !**********************************************************************!
300     ! Function: calc_ICpoR                                                 !
301     ! Purpose: Integrate the polynomial form of the specific heat.         !
302     !                                                                      !
303     !**********************************************************************!
304           PURE DOUBLE PRECISION FUNCTION calc_ICpoR0(T, A, REF_ICpoR)
305     
306           implicit none
307     
308     ! Dummy Arguments:
309     !-----------------------------------------------------------------------
310     ! Evaluation temperaure (K)
311           DOUBLE PRECISION, intent(in) :: T
312     ! Polynomial coefficients.
313           DOUBLE PRECISION, intent(in) :: A(1:5)
314     ! Referece Integral
315           DOUBLE PRECISION, intent(in) :: REF_ICpoR
316     
317     ! Local Variables.
318     !-----------------------------------------------------------------------
319     ! Integral of specific heat polynomial (from 0 to T) over T
320           DOUBLE PRECISION ICpoRT
321     
322     !-----------------------------------------------------------------------
323     
324           ICpoRT = (((A(5)*T/5.0d0 + A(4)/4.0d0)*T + A(3)/3.0d0)*T +       &
325              A(2)/2.0d0)*T + A(1)
326     
327           calc_ICpoR0 = T*ICpoRT - REF_ICpoR
328     
329           RETURN
330           END FUNCTION calc_ICpoR0
331     
332     
333     
334     !**********************************************************************!
335     ! Function: calc_H0oR                                                  !
336     ! Purpose: Calculate the heat of formation from the first six poly-    !
337     !          nomial coefficients.                                        !
338     !                                                                      !
339     ! >>> This function is currently unused.                               !
340     !                                                                      !
341     !**********************************************************************!
342           PURE DOUBLE PRECISION FUNCTION calc_H0oR(T, Th, Tl, Tc, Ah, Al)
343     
344           implicit none
345     
346     ! Dummy Arguments:
347     !-----------------------------------------------------------------------
348     ! Polynomial coefficients. (High/Low)
349           DOUBLE PRECISION, intent(in) :: Ah(7), Al(7)
350     ! Temperature ranges of polynomials.
351           DOUBLE PRECISION, intent(in) :: Th   ! Max temp (for Ahigh)
352           DOUBLE PRECISION, intent(in) :: Tl   ! Min temp (for Alow)
353           DOUBLE PRECISION, intent(in) :: Tc   ! switch from low to high
354     ! Evaluation temperaure (K)
355           DOUBLE PRECISION, intent(in) :: T
356     
357     ! Local Variables.
358     !-----------------------------------------------------------------------
359     
360           !ICp = calc_ICpoR(T, Th, Tl, Tc, Ah, Al)
361           !If (T < Tc) then
362           !  calc_H0oR = ICp + Al(6)
363           !else
364           !  calc_H0oR = ICp + Ah(6)
365           !endif
366     
367           calc_H0oR = 0.0
368     
369           return
370           END FUNCTION calc_H0oR
371     
372     
373     !**********************************************************************!
374     ! SUBROUTINE: replaceTab                                               !
375     ! Purpose: Replace all instances of a tab with a single space.         !
376     !**********************************************************************!
377           SUBROUTINE replaceTab(C)
378     
379           implicit none
380     
381     ! Dummy Arguments:
382     !-----------------------------------------------------------------------
383     ! Incoming string that will have tabs removed.
384           CHARACTER(len=*) :: C
385     
386     ! Local Variables:
387     !-----------------------------------------------------------------------
388     ! Loop counter
389           INTEGER :: I
390     
391           DO I = 1, len(C)
392             IF(C(I:I) == CHAR(9)) C(I:I)=' '
393           ENDDO
394     
395           RETURN
396           END SUBROUTINE replaceTab
397     
398     
399     !**********************************************************************!
400     ! SUBROUTINE: trimTab                                                  !
401     ! Purpose: Search a string for the first instance of a tab. The        !
402     !          location of the tab and all remaing string entries are      !
403     !          replaced with blank spaces.                                 !
404     !**********************************************************************!
405           SUBROUTINE trimTab(C)
406     
407           implicit none
408     
409     ! Dummy Arguments:
410     !-----------------------------------------------------------------------
411     ! Incoming string that will have tabs removed.
412           CHARACTER(len=*) :: C
413     
414     ! Local Variables:
415     !-----------------------------------------------------------------------
416     ! Loop counter
417           INTEGER :: I
418     ! Logical indicating that a tab was located.
419           LOGICAL :: tabFound
420     
421     ! Initialize flag
422           tabFound = .FALSE.
423     
424     ! Look at each entry of the string. Once a tab is located, the rest of
425     ! the string is replaced by blank spaces.
426           DO I = 1, len(C)
427             IF(C(I:I) == CHAR(9) ) tabFound = .TRUE.
428             if(tabFound) C(I:I)=' '
429           ENDDO
430     
431           RETURN
432           END SUBROUTINE trimTab
433     
434         END MODULE read_thermochemical
435