File: /nfs/home/0/users/jenkins/mfix.git/model/read_database.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: read_database(Ier)                                     C
4     !  Purpose: read thermochemical database                               C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 21-OCT-05  C
7     !                                                                      C
8     !  Modification 1: J. Musser                          Date: 02-May-11  C
9     !  Purpose: Provided support for DEM access to database.               C
10     !                                                                      C
11     !  Modification 2: J. Musser                          Date: 02-Oct-12  C
12     !  Purpose: Calls to READ_DATABASE were moved to CHECK_DATA_04/05      C
13     !  duing input data integrity checks for the gas and solids phases.    C
14     !  Rather than looping through all species for each phase, the model   C
15     !  (TFM/DEM), phase index, species index, species name, and molecular  C
16     !  weight are passed as dummy arguments so that only infomration for   C
17     !  referenced species (lName) is obtained.                             C                             C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE READ_DATABASE(lM, lN, lName, lMW)
21     
22           USE param
23           USE param1
24           USE physprop
25           USE constant
26           USE compar
27           USE rxns
28           USE funits
29           USE discretelement
30           USE des_rxns
31           USE read_thermochemical, only: read_therm, calc_ICpoR
32     
33           use error_manager
34     
35           IMPLICIT NONE
36     
37     ! Phase and species indices
38           INTEGER, INTENT(IN) :: lM, lN
39     ! Species name from data file.
40           CHARACTER(len=*), INTENT(IN) :: lName
41     ! Species molecular weight from the data file (if any)
42           DOUBLE PRECISION, INTENT(INOUT) :: lMW
43     
44     ! Molecular weight read from database.
45           DOUBLE PRECISION dbMW
46     
47     ! Error message returned from Read_Therm and sent to calling routine
48           INTEGER IER
49     ! File unit of databases searched
50           INTEGER FUNIT
51     ! Loop counter for checking all three locations for data
52           INTEGER FILE
53     ! Input/Output error status ( 0 is no error)
54           INTEGER IOS
55     
56     ! Identifies if an error was detected.
57           LOGICAL ErrorFlag
58     
59     ! File name of Burcat and Ruscic database
60           CHARACTER(len=10) :: THERM = 'BURCAT.THR'
61     ! Full path to model directory
62           INCLUDE 'mfix_directory_path.inc'
63     ! Full path to Burcat and Ruscic database
64           CHARACTER(len=147) FILENAME
65     
66     ! Tcom +/- SMALL_NUMBER: This is done so that the specific heat
67     ! polynomial can be evaluated at Tcom with the high and low
68     ! coefficients.
69           DOUBLE PRECISION :: xTc
70     
71     ! Various integrations of the specific heat polynomials:
72           DOUBLE PRECISION :: ICpoR_TrL  ! 0.0 --> Tref using Alow
73           DOUBLE PRECISION :: ICpoR_TcL  ! 0.0 --> Tcom using Alow
74           DOUBLE PRECISION :: ICpoR_TcH  ! 0.0 --> Tcom using Ahigh
75     
76           LOGICAL :: testCp = .FALSE.
77     ! Database being searched.
78           CHARACTER(len=256) :: DB
79     
80           CALL INIT_ERR_MSG('READ_DATABASE')
81     
82     ! Initialize the file unit to be used.
83           FUNIT = UNIT_DAT  ! .dat file unit
84     ! Read data from mfix.dat or from BURCAT.THR in run directory or
85     ! mfix directory.
86           FILE = 0
87           DB_LP: DO
88              FILE = FILE + 1
89     ! Check for thermochemical data in the mfix.dat file.
90              IF(FILE == 1) THEN
91                OPEN(UNIT=FUNIT, FILE='mfix.dat', STATUS='OLD', IOSTAT= IOS)
92                IF(IOS /= 0) CYCLE DB_LP
93                DB=''; WRITE(DB,1000) 'mfix.dat'
94     ! Read thermochemical data from the BURCAT.THR database in the local
95     ! run directory.
96              ELSEIF(FILE == 2) THEN
97                 OPEN(UNIT=FUNIT,FILE=TRIM(THERM), STATUS='OLD', IOSTAT= IOS)
98                 IF(IOS /= 0) CYCLE DB_LP
99                 DB=''; WRITE(DB,1000) TRIM(THERM)
100     ! Read thermochemical data from the BURCAT.THR database in the model
101     ! directory (model/thermochemical/BURCAT.THR).
102               ELSEIF(file == 3) THEN
103                 FILENAME = trim(MFIX_PATH)//'/thermochemical/'//TRIM(THERM)
104                 OPEN(UNIT=FUNIT,FILE=TRIM(FILENAME), STATUS='OLD',IOSTAT= IOS)
105                 IF(IOS /= 0) CYCLE DB_LP
106                 DB=''; WRITE(DB,1000) ('/thermochemical/'//TRIM(THERM))
107               ELSE
108                 EXIT DB_LP
109               ENDIF
110     
111              REWIND(UNIT=funit)
112     
113     ! Initialize the error flag
114              IER = 0
115     
116              CALL READ_THERM(FUNIT, lName, Thigh(lM,lN), Tlow(lM,lN),      &
117                 Tcom(lM,lN), dbMW, Ahigh(:,lM,lN), Alow(:,lM,lN),          &
118                 HfrefoR(lM,lN), IER)
119     
120              IF(IER == 0) THEN
121     ! If the user did not supply a value for the gas phase molecular weight
122     ! in the mfix.dat file, use the value from the database.
123                 IF(lMW == UNDEFINED) lMW = dbMW
124     ! There are a number of species with Tlow as 300, for which the
125     ! following calculation will produce an error because T_ref = 298.  So
126     ! slightly extend validity of the correlation.
127                 IF(ABS(Tlow(lM,lN)-T_ref)<=2.0D0 .AND. &
128                    Tlow(lM,lN) > T_ref) Tlow(lM,lN) = T_ref
129     
130     ! Initialize the reference integrals.
131                 ICpoR_l(lM,lN) = ZERO
132                 ICpoR_h(lM,lN) = ZERO
133     
134     ! Calculate the integral of specific heat from zero to Tref using the
135     ! Alow coefficients.
136                    ICpoR_TrL = calc_ICpoR(T_ref, lM, lN, IER)
137     ! Calculate the integral of specific heat from zero to Tcom using the
138     ! Alow coefficients.
139                 xTc = Tcom(lM,lN)-SMALL_NUMBER
140                 ICpoR_TcL = calc_ICpoR(xTc, lM, lN, IER)
141     ! Calculate the integral of specific heat from zero to Tcom using the
142     ! Ahigh coefficients.
143                 xTc = Tcom(lM,lN)+SMALL_NUMBER
144                 ICpoR_TcH = calc_ICpoR(xTc, lM, lN, IER)
145     ! Store the integrals in global variables.
146                 ICpoR_l(lM,lN) = ICpoR_TrL
147                 ICpoR_h(lM,lN) = ICpoR_TcH - (ICpoR_TcL - ICpoR_TrL)
148              ENDIF
149     
150              ErrorFlag = .TRUE.
151              IF(IER == 0) THEN
152                 WRITE(ERR_MSG,1001) trim(adjustl(DB)), 'Found!'
153                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
154                 if(testCP) CALL writeCp(lM, lN, lName, lMW)
155                 ErrorFlag = .FALSE.
156                 EXIT DB_LP
157             ELSE
158                 WRITE(ERR_MSG,1001) trim(adjustl(DB)), 'Not found.'
159                 CALL FLUSH_ERR_MSG(HEADER=.FALSE., FOOTER=.FALSE.)
160             ENDIF
161     
162           ENDDO DB_LP
163     
164     ! Error message control.
165     !-----------------------------------------------------------------------
166     ! Write remaining error message if needed.
167           IF(ErrorFlag) THEN
168              CALL FLUSH_ERR_MSG(HEADER=.FALSE.)
169              WRITE(ERR_MSG,1010) trim(lName), trim(FILENAME)
170              CALL FLUSH_ERR_MSG(ABORT=.TRUE.)
171           ENDIF
172     
173           CALL FINL_ERR_MSG
174     
175           RETURN
176     
177     
178     ! Messages
179     !-----------------------------------------------------------------------
180      1000 FORMAT('Checking ',A)
181      1001 FORMAT(8X,A,1X,' :: ',A)
182     
183     ! Error Flags
184     !-----------------------------------------------------------------------
185      1010 FORMAT('Message 1010: Species "',A,'" was not matched to any ',  &
186              'entry in the',/'thermochemical databases.',2/,'SUGGESTION: ',&
187              'Search the database for the exact species name. The ',/      &
188              'species names are case sensitive and should match the names',&
189              ' in',/'BURCAT.THR exactly excluding trailing blanks and ',   &
190              'tabs. Also verify',/'that the data section in the mfix.dat ',&
191              'file (if any) is below a line',/'that starts with THERMO ',  &
192              'DATA.',2/'Database location:', /A)
193     
194      1020 FORMAT('Message 1020: This routine was entered with an invalid ',&
195              'model argument.',/' Acceptable values are TFM and DEM.',/    &
196              ' Passed Argument: MODEL = ',A)
197     
198           END SUBROUTINE READ_DATABASE
199     
200     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
201     !                                                                      C
202     !  Module name: READ_DATABASE0(IER)                                    C
203     !  Purpose: Provides legacy support for rrates files.                  C
204     !                                                                      C
205     !  Author: J. Musser                                  Date: 02-Oct-12  C
206     !                                                                      C
207     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
208           SUBROUTINE READ_DATABASE0(Ier)
209     
210           USE param
211           USE param1
212           USE physprop
213           USE constant
214           USE compar
215           USE rxns
216           USE funits
217           USE discretelement
218           USE des_rxns
219     
220           IMPLICIT NONE
221     
222     ! Error message returned from Read_Therm and sent to calling routine
223           INTEGER IER
224     ! Loop indices for mass phase and species
225           INTEGER M, N
226     ! Loop counter for continuum and discrete species
227           INTEGER Nsp, DES_Nsp
228     
229     ! Return to the calling routine if the database has already been called.
230           IF(database_read)RETURN
231     
232     ! Set the flag identifying that the database has been read.
233           database_read = .TRUE.
234     
235     ! Initialize counters
236           Nsp = 0
237           DES_Nsp = 0
238     
239     ! Read species data for the gas phase.
240     !-----------------------------------------------------------------------
241           DO N = 1, NMAX(0)
242              Nsp = Nsp + 1
243     ! If a species name was not specified in mfix.dat, flag error and exit.
244               IF(SPECIES_NAME(Nsp) == UNDEFINED_C) THEN
245                 WRITE(*,1010) N         ! screen
246                 IF(DMP_LOG) WRITE(UNIT_LOG,1010) N  ! log file
247                  CALL MFIX_EXIT(mypE)
248               ENDIF
249     ! Read the database.
250              CALL READ_DATABASE(0, N, SPECIES_NAME(Nsp), MW_g(N))
251            ENDDO
252     
253     ! Read species data for the continuum solids phases.
254     !-----------------------------------------------------------------------
255     ! Skip reading the database for the continuum solids phase if the
256     ! simulation is only employing discrete solids.
257           IF(.NOT.DISCRETE_ELEMENT .OR. DES_CONTINUUM_HYBRID)THEN
258               DO M = 1, MMAX
259                 DO N = 1, NMAX(M)
260                     Nsp = Nsp + 1
261     ! If a species name was not specified in mfix.dat, flag error and exit.
262                     IF(SPECIES_NAME(Nsp) == UNDEFINED_C)THEN
263                       WRITE(*,1011)'continuum', M, N ! screen
264                       IF(DMP_LOG) WRITE(UNIT_LOG,1011)'continuum', M, N
265                        CALL MFIX_EXIT(mypE)
266                     ENDIF
267                    CALL READ_DATABASE(M, N, SPECIES_NAME(Nsp), MW_s(M,N))
268                  ENDDO   ! N=1, NMAX(M)
269               ENDDO   ! M=1, MMAX
270           ENDIF
271     
272           RETURN
273     
274     ! Error Messages
275     !-----------------------------------------------------------------------
276      1010 FORMAT(/1X,70('*')/, ' From: READ_DATABASE0',/, ' Message: ',    &
277              'No SPECIES_NAME provided for gas phase species ',I3,'.',/' ',&
278              'Check mfix.dat.',/1X,70('*')/)
279      1011 FORMAT(/1X,70('*')/, ' From: READ_DATABASE0',/, ' Message: ',    &
280              'No SPECIES_NAME provided for ',A,' solids phase ',I2,', ',/  &
281              ' species ',I3,'.',/' Check mfix.dat.',/1X,70('*')/)
282     
283           END SUBROUTINE READ_DATABASE0
284     
285     
286     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
287     !                                                                      C
288     !  Module name: read_database(Ier)                                     C
289     !  Purpose: read thermochemical database                               C
290     !                                                                      C
291     !  Author: M. Syamlal                                 Date: 21-OCT-05  C
292     !                                                                      C
293     !  Modification 1: J. Musser                          Date: 02-May-11  C
294     !  Purpose: Provided support for DEM access to database.               C
295     !                                                                      C
296     !  Modification 2: J. Musser                          Date: 02-Oct-12  C
297     !  Purpose: Calls to READ_DATABASE were moved to CHECK_DATA_04/05      C
298     !  duing input data integrity checks for the gas and solids phases.    C
299     !  Rather than looping through all species for each phase, the model   C
300     !  (TFM/DEM), phase index, species index, species name, and molecular  C
301     !  weight are passed as dummy arguments so that only infomration for   C
302     !  referenced species (lName) is obtained.                             C                             C
303     !                                                                      C
304     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
305           SUBROUTINE writeCp(lM, lN, lName, lMW)
306     
307           use param1
308           USE physprop
309           USE compar
310           USE read_thermochemical, only: calc_CpoR, calc_ICpoR, calc_ICpoR0
311     
312     ! Universal gas constant in cal/mol.K
313           use constant, only: RGAS => GAS_CONST_cal
314     
315           IMPLICIT NONE
316     
317     ! Phase and species indices
318           INTEGER, INTENT(IN) :: lM, lN
319     ! Species name from data file.
320           CHARACTER(len=*), INTENT(IN) :: lName
321     ! Species molecular weight from the data file (if any)
322           DOUBLE PRECISION, INTENT(in) :: lMW
323     
324           INTEGER :: IER1, IER2, lc
325     
326           DOUBLE PRECISION :: T
327           DOUBLE PRECISION :: lCP, lICP
328     
329           write(*,"(2/,3x,'Specific Heat report for ',A)")trim(lName)
330     
331           write(*,"(/12x,'Low',9x,'High')")
332           write(*,"(6x,'T',3x,g12.5,2x,g12.5)") Tlow(lM,lN), Thigh(lM,lN)
333           DO lc=1,5
334              write(*,"(4x,'A(',I1,')',2(2x,g12.5))") lc, &
335                 Alow(lc,lM,lN), Ahigh(lc,lM,lN)
336           ENDDO
337           write(*,"('')")
338           write(*,"(5x,'Tcom: ',g12.5)")Tcom(lM,lN)
339           write(*,"('')")
340     
341           write(*,"(5x,'Temperature',8x,'Cp',11x,'ICp')")
342     
343           T = Tcom(lM,lN) - 100.0
344           DO WHILE(T <= Tcom(lM,lN) - SMALL_NUMBER)
345     
346              IER1 = 0
347              IER2 = 0
348     
349              write(*,"(7x,g12.5)",ADVANCE="NO") T
350              lCP  = calc_CpoR(T, lM, lN, IER1) * RGAS / lMW
351              lICP = calc_ICpoR(T, lM, lN, IER2) * RGAS / lMW
352              write(*,"(2(3x,g12.5))",ADVANCE="NO")lCP, lICP
353     
354              IF(IER1 /= 0) write(*,"(3x,'Cp Error!')",ADVANCE="NO")
355              IF(IER2 /= 0) write(*,"(3x,'ICp Error!')",ADVANCE="NO")
356              write(*,"('')")
357     
358              T = T + 5.0
359           ENDDO
360     
361     
362           T = Tcom(lM,lN) + SMALL_NUMBER
363           DO WHILE(T <= Tcom(lM,lN) + 100.0)
364     
365              IER1 = 0
366              IER2 = 0
367     
368              write(*,"(7x,g12.5)",ADVANCE="NO") T
369              lCP  = calc_CpoR(T, lM, lN, IER1) * RGAS / lMW
370              lICP = calc_ICpoR(T, lM, lN, IER2) * RGAS / lMW
371              write(*,"(2(3x,g12.5))",ADVANCE="NO")lCP, lICP
372     
373              IF(IER1 /= 0) write(*,"(3x,'Cp Error!')",ADVANCE="NO")
374              IF(IER2 /= 0) write(*,"(3x,'ICp Error!')",ADVANCE="NO")
375              write(*,"('')")
376     
377              T = T + 5.0
378           ENDDO
379     
380           write(*,"('')")
381           END SUBROUTINE writeCp
382