File: N:\mfix\model\physical_prop.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: PHYSICAL_PROP                                           C
4     !  Purpose: Calculate the indicated physical properties that vary      C
5     !           with time if directed to do so by the corresponding flag   C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 17-JUL-92  C
8     !  Reviewer: P. Nicoletti                             Date: 11-DEC-92  C
9     !                                                                      C
10     !  Revision: Mods for MFIX 2.0 (old name CALC_PHYSPROP)                C
11     !  Author: M. Syamlal                                 Date: 23-APR-96  C
12     !                                                                      C
13     !                                                                      C
14     !  Literature/Document References:                                     C
15     !    Perry, R.H., and C.H. Chilton, Chemical Engineer's Handbook, 5th  C
16     !      edition, McGraw-Hill Kogakusha, Tokyo, 1973.                    C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE PHYSICAL_PROP(IER, LEVEL)
21     
22     ! Modules
23     !---------------------------------------------------------------------//
24           use compar, only: numPEs, myPe, pe_io
25           use funits, only: unit_log
26           use exit, only: mfix_exit
27           use mpi_utility, only: global_all_sum
28           use physprop, only: smax
29           use coeff, only: DENSITY
30           use coeff, only: SP_HEAT
31           use coeff, only: PSIZE
32           implicit none
33     
34     ! Dummy arguments
35     !---------------------------------------------------------------------//
36     ! Global error Flag.
37           INTEGER, intent(inout) :: IER
38           INTEGER, intent(in) :: LEVEL
39     
40     ! Local variables
41     !---------------------------------------------------------------------//
42     ! Arrays for storing errors:
43     ! 100 - Negative gas phase density
44     ! 101 - Negative solids phase density
45     ! 10x - Unclassified
46     ! 900 - Invalid temperature in calc_CpoR
47           INTEGER :: Err_l(0:numPEs-1)  ! local
48           INTEGER :: Err_g(0:numPEs-1)  ! global
49     ! loop index
50           INTEGER :: M
51     !......................................................................!
52     
53     ! Initialize error flags.
54           Err_l = 0
55     
56     ! Calculate density only. This is invoked several times within iterate,
57     ! making it the most frequently called.
58           if(LEVEL == 0) then
59              if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
60                DO M=1,SMAX
61                   if(DENSITY(M)) CALL PHYSICAL_PROP_ROs(M)
62                ENDDO
63     
64     ! Calculate everything except density. This is called at the start of
65     ! each iteration.
66           elseif(LEVEL == 1) then
67              if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
68              DO M=1,SMAX
69                 if(SP_HEAT(M)) CALL PHYSICAL_PROP_CPs(M)
70                 if(PSIZE(M)) CALL PHYSICAL_PROP_Dp(M)
71              ENDDO
72     
73     
74     ! Calculate everything. This is invoked via calc_coeff_all as part of
75     ! the initialization (before starting the time march) and at the start
76     ! of each step step thereafter.
77           elseif(LEVEL == 2) then
78              if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
79              if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
80              DO M=1,SMAX
81                 if(DENSITY(M)) CALL PHYSICAL_PROP_ROs(M)
82                 if(SP_HEAT(M)) CALL PHYSICAL_PROP_CPs(M)
83                 if(PSIZE(M)) CALL PHYSICAL_PROP_Dp(M)
84              ENDDO
85           endif
86     
87     
88     ! In case of negative density force exit from the physical property
89     ! calculation routine and reduce the time step
90           CALL global_all_sum(Err_l, Err_g)
91           IER = maxval(Err_g)
92           if(IER == 0) return
93     
94     
95     ! Error handeling. - Local.
96     ! An invalid temperature was found by calc_CpoR. This is a fatal run-
97     ! time error and forces a call to MFIX_EXIT.
98           IF(IER == 901 .OR. IER == 902) then
99     ! Temperature is now bounded so this error will not occur.
100              if(myPE == PE_IO) then
101                 write(*,2000) IER
102                 write(UNIT_LOG,2000) IER
103              endif
104              CALL MFIX_EXIT(myPE)
105           ENDIF
106     
107           return
108     
109      2000 FORMAT(/1X,70('*')/' From: PHYSICAL_PROP',/' Fatal Error 2000:', &
110              ' calc_CpoR reported an invalid temperature: 0x0', I3/,       &
111              'See Cp.log for details. Calling MFIX_EXIT.',/1X,70('*')/)
112     
113           CONTAINS
114     
115     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
116     !                                                                      C
117     !  Subroutine: PHYSICAL_PROP_ROg                                       C
118     !  Purpose: Calculate the gas phase density.                           C
119     !                                                                      C
120     !  Author: J. Musser                                  Date: 28-JUN-13  C
121     !                                                                      C
122     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
123           SUBROUTINE PHYSICAL_PROP_ROg
124     
125     ! Global variables:
126     !---------------------------------------------------------------------//
127           use compar, only: ijkstart3, ijkend3
128     ! Equation of State - GAS
129           use eos, only: EOSG
130     ! Gas phase species mass fractions.
131           use fldvar, only: X_g
132     ! Gas phase temperature.
133           use fldvar, only: T_g
134     ! Gas phase density (compressible).
135           use fldvar, only: RO_g
136     ! Gas phase pressure.
137           use fldvar, only: P_g
138     ! Gas phase volume fraction.
139           use fldvar, only: EP_g
140     ! Gas phase material density.
141           use fldvar, only: ROP_g
142     
143           use functions, only: wall_at
144     ! Run time flag for generating negative gas density log files
145           use output, only: REPORT_NEG_DENSITY
146           use param1, only: undefined, zero, one
147           use physprop, only: database_read, mw_g, mw_mix_g, mw_avg, nmax
148     ! Maximum value for molecular weight (divided by one)
149           use toleranc, only: OMW_MAX
150     ! Detect NaN in density for compressible flows.
151           use utilities, only: mfix_isnan
152     ! invoke user defined quantity
153           USE usr_prop, only: usr_rog, calc_usr_prop
154           USE usr_prop, only: gas_density
155           IMPLICIT NONE
156     
157     ! Local Variables:
158     !---------------------------------------------------------------------//
159     ! Average molecular weight
160           DOUBLE PRECISION :: MW
161     ! Loop indicies
162           INTEGER :: IJK   ! Computational cell
163     ! Flag to write log header
164           LOGICAL :: wHeader
165     !......................................................................!
166     
167     ! Ensure that the database was read. This *should* have been caught by
168     ! check_gas_phase but this call remains to prevent an accident.
169           IF(.NOT.database_read) call read_database0(IER)
170     
171     ! User-defined function
172           IF(USR_ROg) THEN
173              CALL CALC_USR_PROP(Gas_Density,lm=0,lerr=err_l)
174              RETURN
175           ENDIF
176     
177     ! Initialize:
178           wHeader = .TRUE.
179     
180           IJK_LP: DO IJK = IJKSTART3, IJKEND3
181              IF(WALL_AT(IJK)) cycle IJK_LP
182              IF (MW_AVG == UNDEFINED) THEN
183     ! Calculating the average molecular weight of the fluid.
184                 MW = SUM(X_G(IJK,:NMAX(0))/MW_G(:NMAX(0)))
185                 MW = ONE/MAX(MW,OMW_MAX)
186                 MW_MIX_G(IJK) = MW
187     ! Calculate the fluid density and bulk density
188                 RO_G(IJK) = EOSG(MW,P_G(IJK),T_G(IJK))
189                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
190              ELSE
191                 RO_G(IJK) = EOSG(MW_AVG,P_G(IJK),T_G(IJK))
192                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
193              ENDIF
194     
195              IF(RO_G(IJK) < ZERO .or. mfix_isnan(RO_G(IJK))) THEN
196                 Err_l(myPE) = 100
197                 IF(REPORT_NEG_DENSITY) CALL ROgErr_LOG(IJK, wHeader)
198              ENDIF
199           ENDDO IJK_LP
200     
201           RETURN
202           END SUBROUTINE PHYSICAL_PROP_ROg
203     
204     
205     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
206     !                                                                      C
207     !  Subroutine: PHYSICAL_PROP_ROs                                       C
208     !  Purpose: Calculate solids phase (variable) density.                 C
209     !                                                                      C
210     !  Author: J. Musser                                  Date: 28-JUN-13  C
211     !                                                                      C
212     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
213           SUBROUTINE PHYSICAL_PROP_ROs(M)
214     
215     ! Global variables:
216     !---------------------------------------------------------------------//
217           use compar, only: ijkstart3, ijkend3
218     ! Equation of State - Solid
219           use eos, only: EOSS
220     ! Solid phase species mass fractions.
221           use fldvar, only: X_s, ROP_s, RO_s
222     ! Solid phase density (variable).
223           use fldvar, only: ROP_s, RO_s
224     
225           use functions, only: wall_at
226     ! Run time flag for generating negative density log files.
227           use output, only: REPORT_NEG_DENSITY
228           use param1, only: undefined, zero
229     ! Baseline/Unreaced solids density
230           use physprop, only: RO_s0
231     ! Initial mass fraction of inert species
232           use physprop, only: X_S0
233     ! Index of inert solids phase species.
234           use physprop, only: INERT_SPECIES
235     ! Inert solids phase species mass fraction in dilute region.
236           use physprop, only: DIL_INERT_X_VSD
237     ! Factor to define dilute region where DIL_INERT_X_VSD is used
238           use physprop, only: DIL_FACTOR_VSD
239     ! Flag for variable solids density.
240           use run, only: SOLVE_ROs
241     ! Minimum solids volume fraction
242           use toleranc, only: DIL_EP_s
243     ! invoke user defined quantity
244           USE usr_prop, only: usr_ros, calc_usr_prop
245           USE usr_prop, only: solids_density
246           IMPLICIT NONE
247     
248     ! Dummy arguments
249     !---------------------------------------------------------------------//
250     ! solids phase index
251           INTEGER, INTENT(IN) :: M
252     
253     ! Local Variables:
254     !---------------------------------------------------------------------//
255     ! Loop indicies
256           INTEGER :: IJK
257     ! Index of inert species
258           INTEGER :: IIS
259     ! Flag to write log header
260           LOGICAL :: wHeader
261     ! Minimum bulk density
262           DOUBLE PRECISION :: minROPs
263     !......................................................................!
264     
265     ! User defined function
266           IF(USR_ROs(m)) THEN
267              CALL CALC_USR_PROP(Solids_Density,lm=M,lerr=err_l)
268              RETURN
269           ENDIF
270     
271           IF(SOLVE_ROs(M)) THEN
272     ! Initialize header flag.
273              wHeader = .TRUE.
274     ! Set the index of the inert species
275              IIS = INERT_SPECIES(M)
276     ! Calculate the minimum solids denisty.
277              minROPs = RO_s0(M)*(DIL_FACTOR_VSD*DIL_EP_s)
278     
279     ! Calculate the solids denisty over all cells.
280              IJK_LP: DO IJK = IJKSTART3, IJKEND3
281                 IF(WALL_AT(IJK)) cycle IJK_LP
282                 IF(ROP_s(IJK,M) > minROPs) THEN
283                    RO_S(IJK,M) = EOSS(RO_s0(M), X_s0(M,IIS), &
284                       X_s(IJK,M,IIS))
285                 ELSE
286     !               RO_s(IJK,M) = RO_s0(M)
287                    RO_S(IJK,M) = EOSS(RO_s0(M), X_s0(M,IIS), &
288                       DIL_INERT_X_VSD(M))
289                 ENDIF
290     
291     ! Report errors.
292                 IF(RO_S(IJK,M) <= ZERO) THEN
293                    Err_l(myPE) = 101
294                    IF(REPORT_NEG_DENSITY) CALL ROsErr_LOG(IJK, M, wHeader)
295                 ENDIF
296              ENDDO IJK_LP
297           ENDIF
298     
299           RETURN
300           END SUBROUTINE PHYSICAL_PROP_ROs
301     
302     
303     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
304     !                                                                      C
305     !  Subroutine: PHYSICAL_PROP_CPg                                       C
306     !  Purpose: Calculate the gas phase constant pressure specific heat.   C
307     !                                                                      C
308     !  Author: J. Musser                                  Date: 28-JUN-13  C
309     !                                                                      C
310     ! Notes:                                                               C
311     !  > Unit conversion: 1 cal = 4.183925 J                               C
312     !                                                                      C
313     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
314           SUBROUTINE PHYSICAL_PROP_CPg
315     
316     ! Global Variables:
317     !---------------------------------------------------------------------//
318           use compar, only: ijkstart3, ijkend3
319     ! Universal gas constant in cal/mol.K
320           use constant, only: RGAS => GAS_CONST_cal
321     ! Gas phase species mass fractions.
322           use fldvar, only: X_g
323     ! Gas phase temperature.
324           use fldvar, only: T_g
325     
326           use functions, only: wall_at
327           use output, only: REPORT_NEG_SPECIFICHEAT
328           use param1, only: undefined, zero
329           use physprop, only: database_read
330           use physprop, only: MW_g, c_pg, nmax
331     ! Function to calculate Cp over gas constant R
332           use read_thermochemical, only: calc_CpoR
333     ! Units: CGS/SI
334           use run, only: UNITS
335     ! invoke user defined quantity
336           USE usr_prop, only: usr_cpg, calc_usr_prop
337           USE usr_prop, only: gas_specificheat
338           IMPLICIT NONE
339     
340     ! Local Variables:
341     !---------------------------------------------------------------------//
342     ! Species specific heat.
343           DOUBLE PRECISION :: lCp
344     ! Loop indicies
345           INTEGER :: IJK, NN
346     ! Flag to write log header
347           LOGICAL :: wHeader
348     ! Error flag returned from calc_CpoR
349           INTEGER :: lCP_Err
350           INTEGER :: gCP_Err
351     !......................................................................!
352     
353     ! Ensure that the database was read. This *should* have been caught by
354     ! check_gas_phase but this call remains to prevent an accident.
355           IF(.NOT.database_read) CALL read_database0(IER)
356     
357     ! User defined function
358           IF(USR_CPg) THEN
359              CALL CALC_USR_PROP(Gas_SpecificHeat,lm=0,lerr=err_l)
360              RETURN
361           ENDIF
362     
363     ! Initialize
364           wHeader = .TRUE.
365           gCP_Err = 0
366           lCP_Err = 0
367           IJK_LP: DO IJK = IJKSTART3, IJKEND3
368              IF(WALL_AT(IJK)) CYCLE IJK_LP
369     ! Calculating an average specific heat for the fluid.
370              C_PG(IJK) = ZERO
371              DO NN = 1, NMAX(0)
372                 lCp = calc_CpoR(T_G(IJK), 0, NN)
373     ! calc_cpor used to return error flag (101/102) on out of bound
374     ! temperature. however, temperature is now bounded in routine.
375                 gCP_Err = max(gCP_Err, lCP_Err)
376                 C_PG(IJK) = C_PG(IJK) + X_g(IJK,NN) * lCp * RGAS / MW_g(NN)
377              ENDDO
378     
379     ! report errors
380              IF(C_PG(IJK) <= ZERO) THEN
381     !            gCP_err = 103
382                 IF(REPORT_NEG_SPECIFICHEAT) CALL CPgErr_LOG(IJK, wHeader)
383              ENDIF
384     
385           ENDDO IJK_LP
386     
387     ! The database calculation always returns cal/g.K thus the following
388     ! conversion is needed if using SI units.
389     ! Only convert useful part of the array. When the array is re-indexed,
390     ! elements past IJKEND3 are UNDEFINED and would keep growing if
391     ! multiplied, leading to overflow.
392           IF(UNITS == 'SI') THEN
393              DO IJK = IJKSTART3, IJKEND3
394                C_PG(IJK) = 4.183925d3 * C_PG(IJK)
395              ENDDO
396           ENDIF
397     
398     ! Increment the error to 900+ to invoke fatal exit.
399           IF(gCP_Err /= 0) Err_l(myPE) = 800 + gCP_Err
400     
401           RETURN
402           END SUBROUTINE PHYSICAL_PROP_CPg
403     
404     
405     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
406     !                                                                      C
407     !  Subroutine: PHYSICAL_PROP_CPs                                       C
408     !  Purpose: Calculate solids phase constant pressure specific heat.    C
409     !                                                                      C
410     !  Author: J. Musser                                  Date: 28-JUN-13  C
411     !                                                                      C
412     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
413           SUBROUTINE PHYSICAL_PROP_CPs(M)
414     
415     ! Global Variables:
416     !---------------------------------------------------------------------//
417           use compar, only: ijkstart3, ijkend3
418     ! Universal gas constant in cal/mol.K
419           use constant, only: RGAS => GAS_CONST_cal
420     ! Solids temperature
421           use fldvar, only: T_s
422     ! Solids species mass fractions
423           use fldvar, only: X_s
424     
425           use functions, only: wall_at
426           use output, only: REPORT_NEG_SPECIFICHEAT
427           use param1, only: undefined, zero
428           use physprop, only: database_read
429           use physprop, only: MW_s, c_ps, nmax
430     ! Function to calculate Cp over gas constant R
431           use read_thermochemical, only: calc_CpoR
432     ! Units: CGS/SI
433           use run, only: UNITS
434     ! invoke user defined quantity
435           USE usr_prop, only: usr_cps, calc_usr_prop
436           USE usr_prop, only: solids_specificheat
437           IMPLICIT NONE
438     
439     ! Dummy arguments
440     !---------------------------------------------------------------------//
441     ! solids phase index
442           INTEGER, INTENT(IN) :: M
443     
444     ! Local variables:
445     !---------------------------------------------------------------------//
446     ! Local value for species specific heat (Cp)
447           DOUBLE PRECISION :: lCp
448     ! Loop indicies
449           INTEGER :: IJK, NN
450     ! Flag to write log header
451           LOGICAL :: wHeader
452     ! Error flag returned from calc_CpoR
453           INTEGER :: lCP_Err
454           INTEGER :: gCP_Err
455     !......................................................................!
456     
457     ! Ensure that the database was read. This *should* have been caught by
458     ! check_solids_common_all but this call remains to prevent an accident.
459           IF(.NOT.database_read) CALL read_database0(IER)
460     
461     ! User defined function
462           IF(USR_CPs(M)) THEN
463              CALL CALC_USR_PROP(Solids_SpecificHeat,lm=M,lerr=err_l)
464              RETURN
465           ENDIF
466     
467     ! Initialize
468           wHeader = .TRUE.
469           gCP_Err = 0
470           lCP_Err = 0
471           IJK_LP: DO IJK = IJKSTART3, IJKEND3
472              IF(WALL_AT(IJK)) CYCLE IJK_LP
473     ! Calculating an average specific heat for the fluid.
474              C_PS(IJK, M) = ZERO
475              DO NN = 1, NMAX(M)
476                 lCp = calc_CpoR(T_S(IJK,M), M, NN)
477     ! calc_cpor used to return error flag (101/102) based on out of bound
478     ! temperature. however, temperature is now bounded in routine
479                 gCP_Err = max(gCP_Err, lCP_Err)
480                 C_PS(IJK,M) = C_PS(IJK,M) + X_s(IJK,M,NN) * &
481                    (lCp * RGAS / MW_s(M,NN))
482              ENDDO
483     
484     ! report errors
485              IF(C_PS(IJK,M) <= ZERO) THEN
486     !            gCP_err = 104
487                 IF(REPORT_NEG_SPECIFICHEAT) CALL CPsErr_LOG(IJK, M, wHeader)
488              ENDIF
489     
490           ENDDO IJK_LP
491     
492     ! The database calculation always returns cal/g.K thus the following
493     ! conversion is needed if using SI units.
494     ! Only convert useful part of the array. When the array is re-indexed,
495     ! elements past IJKEND3 are UNDEFINED and would keep growing if
496     ! multiplied, leading to overflow.
497           IF(UNITS == 'SI') THEN
498              DO IJK = IJKSTART3, IJKEND3
499                C_PS(IJK,M) = 4.183925d3 * C_PS(IJK,M)
500              ENDDO
501           ENDIF
502     
503     ! Increment the error to 900+ to invoke fatal exit.
504           IF(gCP_Err /= 0) Err_l(myPE) = 800 + gCP_Err
505     
506           RETURN
507           END SUBROUTINE PHYSICAL_PROP_CPs
508     
509     
510     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
511     !                                                                      C
512     !  Subroutine: PHYSICAL_PROP_DP                                        C
513     !  Purpose: Calculate solids phase diameter.                           C
514     !                                                                      C
515     !  Author: J. Musser                                  Date: 28-JUN-13  C
516     !                                                                      C
517     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
518           SUBROUTINE PHYSICAL_PROP_Dp(M)
519     
520     ! Modules
521     !---------------------------------------------------------------------//
522           use compar, only: ijkstart3, ijkend3
523           use scalars, only: phase4scalar
524           use fldvar, only: scalar
525           use fldvar, only: D_p, EP_S
526           use functions, only: wall_at
527           use param1, only: small_number
528           use run, only: CALL_DQMOM
529           implicit none
530     
531     ! Dummy arguments
532     !---------------------------------------------------------------------//
533     ! solids phase index
534           INTEGER, INTENT(IN) :: M
535     
536     ! Local variables
537     !---------------------------------------------------------------------//
538     ! Loop indicies
539           INTEGER :: IJK   ! Computational cell
540     ! Map from true index to map.
541           INTEGER :: lM
542     !......................................................................!
543     
544           IF(.NOT.CALL_DQMOM) return
545     
546           lM = phase4scalar(M) ! Map from scalar eq to solids phase
547     
548           IJK_LP: DO IJK = IJKSTART3, IJKEND3
549              IF(WALL_AT(IJK)) CYCLE IJK_LP
550              IF(EP_s(IJK,lM) > small_number) D_p(IJK,M)= Scalar(IJK,lM)
551           ENDDO IJK_LP
552     
553           RETURN
554           END SUBROUTINE PHYSICAL_PROP_Dp
555     
556     
557           END SUBROUTINE PHYSICAL_PROP
558     
559     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
560     !                                                                      C
561     !  Subroutine: ROgErr_LOG                                              C
562     !  Purpose: Record information about the location and conditions that  C
563     !           resulted in a negative gas phase density.                  C
564     !                                                                      C
565     !  Author: J. Musser                                  Date: 28-JUN-13  C
566     !                                                                      C
567     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
568           SUBROUTINE ROgErr_LOG(IJK, tHeader)
569     
570     ! Global Variables:
571     !---------------------------------------------------------------------//
572           use compar, only: mype, numPEs
573     ! Cutcell data.
574           use cutcell, only: cartesian_grid
575           use cutcell, only: cut_cell_at, small_cell_at
576           use cutcell, only: xg_e, yg_n, zg_t
577     ! Simulation time
578           use run, only: TIME
579     ! Gas phase temperature.
580           use fldvar, only: T_g
581     ! Gas phase density (compressible).
582           use fldvar, only: RO_g
583     !
584           use indices, only: i_of, j_of, k_of
585     ! Gas phase pressure.
586           use fldvar, only: P_g
587           IMPLICIT NONE
588     
589     ! Dummy arguments
590     !---------------------------------------------------------------------//
591           INTEGER, intent(in) :: IJK
592           LOGICAL, intent(inout) :: tHeader
593     
594     ! Local variables:
595     !---------------------------------------------------------------------//
596           LOGICAL :: lExists
597           CHARACTER(LEN=255) :: lFile
598           INTEGER, parameter :: lUnit = 4868
599           LOGICAL, save :: fHeader = .TRUE.
600     !......................................................................!
601     
602           lFile = '';
603           if(numPEs > 1) then
604              write(lFile,"('ROgErr_',I4.4,'.log')") myPE
605           else
606              write(lFile,"('ROgErr.log')")
607           endif
608           inquire(file=trim(lFile),exist=lExists)
609           if(lExists) then
610              open(lUnit,file=trim(adjustl(lFile)),  &
611                 status='old', position='append', convert='big_endian')
612           else
613              open(lUnit,file=trim(adjustl(lFile)), status='new',&
614                   convert='big_endian')
615           endif
616     
617           if(fHeader) then
618              write(lUnit,1000)
619              fHeader = .FALSE.
620           endif
621     
622           if(tHeader) then
623              write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
624              tHeader = .FALSE.
625           endif
626     
627           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
628           write(lUnit,"(6x,A,1X,g12.5)",ADVANCE='NO') 'RO_g:', RO_g(IJK)
629           write(lUnit,"(2x,A,1X,g12.5)",ADVANCE='NO') 'P_g:', P_g(IJK)
630           write(lUnit,"(2x,A,1X,g12.5)") 'T_g:', T_g(IJK)
631           if(CARTESIAN_GRID) then
632              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
633              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
634              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
635                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
636           endif
637     
638           close(lUnit)
639     
640           RETURN
641     
642      1000 FORMAT(2X,'One or more cells have reported a negative/NaN gas',  &
643              ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
644              ' lower UR_FAC(1) in mfix.dat.')
645     
646      1001 FORMAT(/4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
647     
648           END SUBROUTINE ROgErr_LOG
649     
650     
651     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
652     !                                                                      !
653     !  Subroutine: ROsErr_LOG                                              !
654     !  Purpose: Record information about the location and conditions that  !
655     !           resulted in a negative solids phase density.               !
656     !                                                                      !
657     !  Author: J. Musser                                  Date: 09-Oct-13  !
658     !                                                                      !
659     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
660           SUBROUTINE ROsErr_LOG(IJK, M, tHeader)
661     
662     ! Global variables:
663     !---------------------------------------------------------------------//
664           use compar, only: mype, numPEs
665     ! Cutcell data.
666           use cutcell, only: cartesian_grid
667           use cutcell, only: cut_cell_at, small_cell_at
668           use cutcell, only: xg_e, yg_n, zg_t
669     ! Solid phase species mass fractions.
670           use fldvar, only: X_s
671     ! Solid phase density (variable).
672           use fldvar, only: RO_s
673     !
674           use indices, only: i_of, j_of, k_of
675     ! Baseline/Unreaced solids density
676           use physprop, only: RO_s0
677     ! Initial mass fraction of inert species
678           use physprop, only: X_S0
679     ! Index of inert solids phase species.
680           use physprop, only: INERT_SPECIES
681     ! Simulation time
682           use run, only: TIME
683           IMPLICIT NONE
684     
685     ! Dummy arguments
686     !---------------------------------------------------------------------//
687     ! Fluid cell and solids phase indices
688           INTEGER, intent(in) :: IJK, M
689     ! Flag to output header
690           LOGICAL, intent(inout) :: tHeader
691     
692     ! Local Variables:
693     !---------------------------------------------------------------------//
694     ! Local aliase for inert species index
695           INTEGER :: NN
696     ! Local file values.
697           LOGICAL :: lExists
698           CHARACTER(LEN=255) :: lFile
699           INTEGER, parameter :: lUnit = 4868
700           LOGICAL, save :: fHeader = .TRUE.
701     !......................................................................!
702     
703           lFile = '';
704           if(numPEs > 1) then
705              write(lFile,"('ROsErr_',I4.4,'.log')") myPE
706           else
707              write(lFile,"('ROsErr.log')")
708           endif
709           inquire(file=trim(lFile),exist=lExists)
710           if(lExists) then
711              open(lUnit,file=trim(adjustl(lFile)), &
712                 status='old', position='append',convert='big_endian')
713           else
714              open(lUnit,file=trim(adjustl(lFile)), status='new',&
715                   convert='big_endian')
716           endif
717     
718           if(fHeader) then
719              write(lUnit,1000)
720              fHeader = .FALSE.
721           endif
722     
723           if(tHeader) then
724              write(lUnit,"(/2x,'Simulation time: ',g12.5,'  Phase: ',I2)")&
725                  TIME, M
726              tHeader = .FALSE.
727           endif
728     
729           NN = INERT_SPECIES(M)
730           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
731           write(lUnit,"(6x,A,1X,g12.5)",advance='no') 'RO_s:', RO_s(IJK,M)
732           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'Base:', RO_s0(M)
733           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'X_s0:', X_s0(M,NN)
734           write(lUnit,"(2x,A,1X,g12.5)") 'X_s:', X_s(IJK,M,NN)
735     
736           if(CARTESIAN_GRID) then
737              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
738              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
739              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
740                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
741           endif
742     
743           close(lUnit)
744     
745           RETURN
746     
747      1000 FORMAT(2X,'One or more cells have reported a negative solids',   &
748              ' density (RO_s(IJK)).')
749     
750      1001 FORMAT( 4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
751     
752           END SUBROUTINE ROsErr_LOG
753     
754     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
755     !                                                                      !
756     !  Subroutine: CPgErr_LOG                                              !
757     !  Purpose: Record information about the location and conditions that  !
758     !           resulted in a zero or negative gas specific heat.          !
759     !                                                                      !
760     !                                                                      !
761     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
762           SUBROUTINE CPgErr_LOG(IJK, tHeader)
763     
764     ! Global variables:
765     !---------------------------------------------------------------------//
766           use compar, only: mype, numPEs
767     ! Cutcell data.
768           use cutcell, only: cartesian_grid
769           use cutcell, only: cut_cell_at, small_cell_at
770           use cutcell, only: xg_e, yg_n, zg_t
771     ! Gas phase species mass fractions.
772           use fldvar, only: T_g, X_g, EP_g
773     ! Gas phase density (variable).
774           use fldvar, only: RO_g
775     !
776           use indices, only: i_of, j_of, k_of
777     
778           use physprop, only: MW_g, C_pg, nmax
779     ! Simulation time
780           use run, only: TIME
781           IMPLICIT NONE
782     
783     ! Dummy arguments
784     !---------------------------------------------------------------------//
785     ! Fluid cell and solids phase indices
786           INTEGER, intent(in) :: IJK
787     ! Flag to output header
788           LOGICAL, intent(inout) :: tHeader
789     
790     ! Local Variables:
791     !---------------------------------------------------------------------//
792     ! Local aliase for inert species index
793           INTEGER :: NN
794     ! Local file values.
795           LOGICAL :: lExists
796           CHARACTER(LEN=255) :: lFile
797           INTEGER, parameter :: lUnit = 4868
798           LOGICAL, save :: fHeader = .TRUE.
799           CHARACTER(LEN=7) :: X_gN
800     !......................................................................!
801     
802           lFile = '';
803           if(numPEs > 1) then
804              write(lFile,"('CPgErr_',I4.4,'.log')") myPE
805           else
806              write(lFile,"('CPgErr.log')")
807           endif
808           inquire(file=trim(lFile),exist=lExists)
809           if(lExists) then
810              open(lUnit,file=trim(adjustl(lFile)), &
811                 status='old', position='append',convert='big_endian')
812           else
813              open(lUnit,file=trim(adjustl(lFile)), status='new',&
814                   convert='big_endian')
815           endif
816     
817           if(fHeader) then
818              write(lUnit,1000)
819              fHeader = .FALSE.
820           endif
821     
822           if(tHeader) then
823              write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
824              tHeader = .FALSE.
825           endif
826     
827           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
828           write(lUnit,"(6x,A,1X,g12.5)",advance='no') 'C_PG:', C_PG(IJK)
829           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'EP_G:', EP_g(IJK)
830           write(lUnit,"(2x,A,1X,g12.5)") 'T_G:', T_g(IJK)
831           write(lUnit,"(6x,A,1X)",advance='no') 'X_gN:'
832           DO NN = 1,NMAX(0)
833              write(lUnit,"(g12.5,2X)",advance='no') X_g(IJK,NN)
834           ENDDO
835           write(lUnit,fmt='(/)')
836     
837           if(CARTESIAN_GRID) then
838              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
839              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
840              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
841                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
842           endif
843     
844           close(lUnit)
845     
846           RETURN
847     
848      1000 FORMAT(2X,'One or more cells have reported a negative or zero',&
849              ' gas specific',/2X,'heat (C_pg(IJK)).')
850     
851      1001 FORMAT( 4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
852     
853           END SUBROUTINE CPgErr_LOG
854     
855     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
856     !                                                                      !
857     !  Subroutine: CPsErr_LOG                                              !
858     !  Purpose: Record information about the location and conditions that  !
859     !           resulted in a zero or negative solids specific heat.       !
860     !                                                                      !
861     !                                                                      !
862     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
863           SUBROUTINE CPsErr_LOG(IJK, M, tHeader)
864     
865     ! Global variables:
866     !---------------------------------------------------------------------//
867           use compar, only: mype, numPEs
868     ! Cutcell data.
869           use cutcell, only: cartesian_grid
870           use cutcell, only: cut_cell_at, small_cell_at
871           use cutcell, only: xg_e, yg_n, zg_t
872     ! Solid phase species mass fractions.
873           use fldvar, only: T_s, X_s, EP_S
874     ! Solid phase density (variable).
875           use fldvar, only: RO_s
876     !
877           use indices, only: i_of, j_of, k_of
878     
879           use physprop, only: MW_s, C_ps, nmax
880     ! Simulation time
881           use run, only: TIME
882           IMPLICIT NONE
883     
884     ! Dummy arguments
885     !---------------------------------------------------------------------//
886     ! Fluid cell and solids phase indices
887           INTEGER, intent(in) :: IJK, M
888     ! Flag to output header
889           LOGICAL, intent(inout) :: tHeader
890     
891     ! Local Variables:
892     !---------------------------------------------------------------------//
893     ! Local aliase for inert species index
894           INTEGER :: NN
895     ! Local file values.
896           LOGICAL :: lExists
897           CHARACTER(LEN=255) :: lFile
898           INTEGER, parameter :: lUnit = 4868
899           LOGICAL, save :: fHeader = .TRUE.
900           CHARACTER(LEN=7) :: X_sN
901     !......................................................................!
902     
903           lFile = '';
904           if(numPEs > 1) then
905              write(lFile,"('CPsErr_',I4.4,'.log')") myPE
906           else
907              write(lFile,"('CPsErr.log')")
908           endif
909           inquire(file=trim(lFile),exist=lExists)
910           if(lExists) then
911              open(lUnit,file=trim(adjustl(lFile)), &
912                 status='old', position='append',convert='big_endian')
913           else
914              open(lUnit,file=trim(adjustl(lFile)), status='new',&
915                   convert='big_endian')
916           endif
917     
918           if(fHeader) then
919              write(lUnit,1000)
920              fHeader = .FALSE.
921           endif
922     
923           if(tHeader) then
924              write(lUnit,"(/2x,'Simulation time: ',g12.5,'  Phase: ',I2)")&
925                  TIME, M
926              tHeader = .FALSE.
927           endif
928     
929           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
930           write(lUnit,"(6x,A,1X,g12.5)",advance='no') 'C_PS:', C_PS(IJK,M)
931           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'EP_S:', EP_s(IJK,M)
932           write(lUnit,"(2x,A,1X,g12.5)") 'T_S:', T_s(IJK,M)
933           write(lUnit,"(6x,A,1X)",advance='no') 'X_sN:'
934           DO NN = 1,NMAX(M)
935              write(lUnit,"(g12.5,2X)",advance='no') X_s(IJK,M,NN)
936           ENDDO
937           write(lUnit,fmt='(/)')
938     
939           if(CARTESIAN_GRID) then
940              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
941              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
942              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
943                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
944           endif
945     
946           close(lUnit)
947     
948           RETURN
949     
950      1000 FORMAT(2X,'One or more cells have reported a negative or zero',&
951              ' solids specific',/2X,'heat (C_ps(IJK)).')
952     
953      1001 FORMAT( 4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
954     
955           END SUBROUTINE CPsErr_LOG
956     
957