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