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