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