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