File: /nfs/home/0/users/jenkins/mfix.git/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 Number: 1                                                  C
11     !  Purpose: Mods for MFIX 2.0 (old name CALC_PHYSPROP)                 C
12     !  Author: M. Syamlal                                 Date: 23-APR-96  C
13     !  Reviewer:                                          Date: dd-mmm-yy  C
14     !                                                                      C
15     !  Revision Number: 2                                                  C
16     !  Purpose: allow SI                                                   C
17     !  Author: S. Dartevelle                              Date: 01-Jul-02  C
18     !  Reviewer:                                          Date: dd-mmm-yy  C
19     !                                                                      C
20     !  Literature/Document References:                                     C
21     !    Perry, R.H., and C.H. Chilton, Chemical Engineer's Handbook, 5th  C
22     !      edition, McGraw-Hill Kogakusha, Tokyo, 1973.                    C
23     !                                                                      C
24     !  Variables referenced: None                                          C
25     !  Variables modified: None                                            C
26     !  Local variables: None                                               C
27     !                                                                      C
28     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
29           SUBROUTINE PHYSICAL_PROP(IER, LEVEL)
30     
31           use compar
32           use funits
33           use geometry
34           use indices
35           use mpi_utility
36           use param1
37           use physprop
38     
39           use coeff, only: DENSITY  ! Density
40           use coeff, only: SP_HEAT  ! Specific heat
41           use coeff, only: PSIZE    ! Particle diameter
42     
43           implicit none
44     
45     ! Dummy arguments
46     !-----------------------------------------------------------------------
47     ! Global error Flag.
48           INTEGER, intent(inout) :: IER
49           INTEGER, intent(in) :: LEVEL
50     
51     ! Local variables
52     !-----------------------------------------------------------------------
53     ! Arrays for storing errors:
54     ! 100 - Negative gas phase density
55     ! 101 - Negative solids phase density
56     ! 10x - Unclassified
57     ! 900 - Invalid temperature in calc_CpoR
58           INTEGER :: Err_l(0:numPEs-1)  ! local
59           INTEGER :: Err_g(0:numPEs-1)  ! global
60     
61     ! Initialize error flags.
62           Err_l = 0
63     
64     ! Calculate density only. This is invoked several times within iterate,
65     ! making it the most frequently called.
66           if(LEVEL == 0) then
67              if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
68              if(any(DENSITY(1:DIMENSION_M))) CALL PHYSICAL_PROP_ROs
69     
70     ! Calculate everything except density. This is called at the start of
71     ! each iteration.
72           elseif(LEVEL == 1) then
73              if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
74     !         if(any(DENSITY(1:DIMENSION_M))) CALL PHYSICAL_PROP_ROs
75              if(any(SP_HEAT(1:DIMENSION_M))) CALL PHYSICAL_PROP_CPs
76              if(any(PSIZE(1:DIMENSION_M))) CALL PHYSICAL_PROP_Dp
77     
78     
79     ! Calculate everything. This is invoked via calc_coeff_all as part of
80     ! the initialization (before starting the time march) and at the start
81     ! of each step step thereafter.
82           elseif(LEVEL == 2) then
83              if(DENSITY(0)) CALL PHYSICAL_PROP_ROg
84              if(SP_HEAT(0)) CALL PHYSICAL_PROP_CPg
85              if(any(DENSITY(1:DIMENSION_M))) CALL PHYSICAL_PROP_ROs
86              if(any(SP_HEAT(1:DIMENSION_M))) CALL PHYSICAL_PROP_CPs
87              if(any(PSIZE(1:DIMENSION_M))) CALL PHYSICAL_PROP_Dp
88           endif
89     
90     
91     ! In case of negative density force exit from the physical property
92     ! calculation routine and reduce the time step
93           CALL global_all_sum(Err_l, Err_g)
94           IER = maxval(Err_g)
95           if(IER == 0) return
96     
97     
98     ! Error handeling. - Local.
99     !-----------------------------------------------------------------------
100     ! An invalid temperature was found by calc_CpoR. This is a fatal run-
101     ! time error and forces a call to MFIX_EXIT.
102           IF(IER == 901 .OR. IER == 902) then
103              if(myPE == PE_IO) then
104                 write(*,2000) IER
105                 write(UNIT_LOG,2000) IER
106              endif
107              CALL MFIX_EXIT(myPE)
108           ENDIF
109     
110           return
111     
112      2000 FORMAT(/1X,70('*')/' From: PHYSICAL_PROP',/' Fatal Error 2000:', &
113              ' calc_CpoR reporetd an invalid temperature: 0x0', I3/,       &
114              'See Cp.log for details. Calling MFIX_EXIT.',/1X,70('*')/)
115     
116           contains
117     
118     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
119     !                                                                      C
120     !  Subroutine: PHYSICAL_PROP_ROg                                       C
121     !  Purpose: Calculate the gas phase density.                           C
122     !                                                                      C
123     !  Author: J. Musser                                  Date: 28-JUN-13  C
124     !  Reviewer:                                          Date:            C
125     !                                                                      C
126     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
127           SUBROUTINE PHYSICAL_PROP_ROg
128     
129     
130     ! Global variables:
131     !-----------------------------------------------------------------------
132     ! Gas phase species mass fractions.
133           use fldvar, only: X_g
134     ! Gas phase temperature.
135           use fldvar, only: T_g
136     ! Gas phase density (compressible).
137           use fldvar, only: RO_g
138     ! Gas phase pressure.
139           use fldvar, only: P_g
140     ! Gas phase volume fraction.
141           use fldvar, only: EP_g
142     ! Gas phase material density.
143           use fldvar, only: ROP_g
144     ! Maximum value for molecular weight (divided by one)
145           use toleranc, only: OMW_MAX
146     ! Run time flag for generating negative gas density log files
147           use run, only: REPORT_NEG_DENSITY
148           use functions
149     ! Equation of State - GAS
150           use eos, only: EOSG
151     
152           implicit none
153     
154     ! Local Variables:
155     !-----------------------------------------------------------------------
156     ! Average molecular weight
157           DOUBLE PRECISION :: MW
158     
159     ! Loop indicies
160           INTEGER :: IJK   ! Computational cell
161     
162     ! Flag to write log header
163           LOGICAL wHeader
164     
165     ! Initialize:
166           wHeader = .TRUE.
167     
168     ! Average molecular weight: Xg1/Mw1 + Xg2/Mw2 + Xg3/Mw3 + ....
169           IF(.NOT.database_read) call read_database0(IER)
170     
171           IJK_LP: DO IJK = IJKSTART3, IJKEND3
172              IF(WALL_AT(IJK)) cycle IJK_LP
173              IF (MW_AVG == UNDEFINED) THEN
174     ! Calculating the average molecular weight of the fluid.
175                 MW = SUM(X_G(IJK,:NMAX(0))/MW_G(:NMAX(0)))
176                 MW = ONE/MAX(MW,OMW_MAX)
177                 MW_MIX_G(IJK) = MW
178     ! Calculate the fluid density and bulk density
179                 RO_G(IJK) = EOSG(MW,P_G(IJK),T_G(IJK))
180                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
181              ELSE
182                 RO_G(IJK) = EOSG(MW_AVG,P_G(IJK),T_G(IJK))
183                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
184              ENDIF
185     
186              IF(RO_G(IJK) < ZERO) THEN
187                 Err_l(myPE) = 100
188                 IF(REPORT_NEG_DENSITY)CALL ROgErr_LOG(IJK, wHeader)
189              ENDIF
190           ENDDO IJK_LP
191     
192     
193           RETURN
194           END SUBROUTINE PHYSICAL_PROP_ROg
195     
196     
197     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
198     !                                                                      C
199     !  Subroutine: PHYSICAL_PROP_ROs                                       C
200     !  Purpose: Calculate solids phase (variable) density.                 C
201     !                                                                      C
202     !  Author: J. Musser                                  Date: 28-JUN-13  C
203     !  Reviewer:                                          Date:            C
204     !                                                                      C
205     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
206           SUBROUTINE PHYSICAL_PROP_ROs
207     
208     ! Global variables:
209     !-----------------------------------------------------------------------
210     ! Solid phase species mass fractions.
211           use fldvar, only: X_s
212     ! Solid phase density (variable).
213           use fldvar, only: ROP_s, RO_s
214     ! Baseline/Unreaced solids density
215           use physprop, only: BASE_ROs
216     ! Initial mass fraction of inert species
217           use physprop, only: X_S0
218     ! Index of inert solids phase species.
219           use physprop, only: INERT_SPECIES
220     ! Flag for variable solids density.
221           use run, only: SOLVE_ROs
222     ! Run time flag for generating negative density log files.
223           use run, only: REPORT_NEG_DENSITY
224     ! Minimum solids volume fraction
225           use toleranc, only: DIL_EP_s
226           use functions
227     ! Equation of State - Solid
228           use eos, only: EOSS
229     
230           implicit none
231     
232     ! Local Variables:
233     !-----------------------------------------------------------------------
234     ! Loop indicies
235           INTEGER :: IJK   ! Computational cell
236           INTEGER :: M     ! Solids phase
237     ! Index of inert species
238           INTEGER :: IIS
239     ! Flag to write log header
240           LOGICAL wHeader
241     
242           DOUBLE PRECISION :: minROPs
243     
244           M_LP: DO M=1, MMAX
245              IF(.NOT.SOLVE_ROs(M)) cycle M_LP
246     ! Initialize header flag.
247              wHeader = .TRUE.
248     ! Set the index of the inert species
249              IIS = INERT_SPECIES(M)
250     ! Calculate the minimum solids denisty.
251              minROPs = BASE_ROs(M)*(10.0d0*DIL_EP_s)
252     
253     ! Calculate the solids denisty over all cells.
254              IJK_LP: DO IJK = IJKSTART3, IJKEND3
255                 IF(WALL_AT(IJK)) cycle IJK_LP
256                 IF(ROP_s(IJK,M) > minROPs) THEN
257                    RO_S(IJK,M) = EOSS(BASE_ROs(M), X_s0(M,IIS),            &
258                       X_s(IJK,M,IIS))
259                 ELSE
260                    RO_s(IJK,M) = BASE_ROs(M)
261                 ENDIF
262     
263     ! Report errors.
264                 IF(RO_S(IJK,M) <= ZERO) THEN
265                    Err_l(myPE) = 101
266                    IF(REPORT_NEG_DENSITY) CALL ROsErr_LOG(IJK, M, wHeader)
267                 ENDIF
268              ENDDO IJK_LP
269           ENDDO M_LP
270     
271     
272           RETURN
273           END SUBROUTINE PHYSICAL_PROP_ROs
274     
275     
276     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
277     !                                                                      C
278     !  Subroutine: PHYSICAL_PROP_CPg                                       C
279     !  Purpose: Calculate the gas phase constant pressure specific heat.   C
280     !                                                                      C
281     !  Author: J. Musser                                  Date: 28-JUN-13  C
282     !  Reviewer:                                          Date:            C
283     !                                                                      C
284     ! Notes:                                                               C
285     !  > Unit conversion: 1 cal = 4.183925 J                               C
286     !                                                                      C
287     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
288           SUBROUTINE PHYSICAL_PROP_CPg
289     
290     ! Global Variables:
291     !-----------------------------------------------------------------------
292     ! Universal gas constant in cal/mol.K
293           use constant, only: RGAS => GAS_CONST_cal
294     ! Maximum value for molecular weight (divided by one)
295           use toleranc, only: OMW_MAX
296     ! Gas phase species mass fractions.
297           use fldvar, only: X_g
298     ! Gas phase temperature.
299           use fldvar, only: T_g
300     ! Units: CGS/SI
301           use run, only: UNITS
302           use functions
303     
304           use read_thermochemical, only: calc_CpoR
305     
306           implicit none
307     
308     ! Local Variables:
309     !-----------------------------------------------------------------------
310     ! Species specific heat.
311           DOUBLE PRECISION :: lCp
312     ! Error flag returned from calc_CpoR
313           INTEGER :: lCP_Err
314           INTEGER :: gCP_Err
315     ! Loop indicies
316           INTEGER :: IJK   ! Computational cell
317           INTEGER :: N     ! Species index
318     
319     !-----------------------------------------------------------------------
320     
321     ! Ensure that the database was read. This *should* have been caught by
322     ! check_data_05 but this call remains to prevent an accident.
323           IF(.NOT.database_read) CALL read_database0(IER)
324     
325           gCP_Err = 0
326           lCP_Err = 0
327     
328           IJK_LP: DO IJK = IJKSTART3, IJKEND3
329              IF(WALL_AT(IJK)) CYCLE IJK_LP
330     ! Calculating an average specific heat for the fluid.
331              C_PG(IJK) = ZERO
332              DO N = 1, NMAX(0)
333                 lCp = calc_CpoR(T_G(IJK), 0, N, lCP_Err)
334                 C_PG(IJK) = C_PG(IJK) + X_g(IJK,N) * lCp * RGAS / MW_g(N)
335                 gCP_Err = max(gCP_Err, lCP_Err)
336              ENDDO
337           ENDDO IJK_LP
338     
339     ! The database calculation always returns cal/g.K thus the following
340     ! conversion is needed if using SI units. ** Vector operation
341           IF(UNITS == 'SI') C_PG = 4.183925d3 * C_PG  !in J/kg K
342     
343     ! Increment the error to 900+ to invoke fatal exit.
344           IF(gCP_Err /= 0) Err_l(myPE) = 800 + gCP_Err
345     
346           RETURN
347           END SUBROUTINE PHYSICAL_PROP_CPg
348     
349     
350     
351     
352     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
353     !                                                                      C
354     !  Subroutine: PHYSICAL_PROP_CPs                                       C
355     !  Purpose: Calculate solids phase constant pressure specific heat.    C
356     !                                                                      C
357     !  Author: J. Musser                                  Date: 28-JUN-13  C
358     !  Reviewer:                                          Date:            C
359     !                                                                      C
360     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
361           SUBROUTINE PHYSICAL_PROP_CPs
362     
363     ! Universal gas constant in cal/mol.K
364           use constant, only: RGAS => GAS_CONST_cal
365     ! Units: CGS/SI
366           use run, only: UNITS
367           use toleranc, only: OMW_MAX
368           use fldvar, only: T_s
369           use fldvar, only: X_s
370           use functions
371           use read_thermochemical, only: calc_CpoR
372     
373           implicit none
374     
375           DOUBLE PRECISION :: lCp
376     
377     ! Loop indicies
378           INTEGER :: IJK   ! Computational cell
379           INTEGER :: M     ! Solids phase
380           INTEGER :: N     ! Species index
381     
382     ! Local error flag indicating that the Cp is out of range.
383           INTEGER :: lCP_Err
384     
385     ! Ensure that the database was read. This *should* have been caught by
386     ! check_data_05 but this call remains to prevent an accident.
387           IF(.NOT.database_read) CALL read_database0(IER)
388     
389           lCP_Err = 0
390     
391           M_LP: DO M=1, MMAX
392              IJK_LP: DO IJK = IJKSTART3, IJKEND3
393                 IF(WALL_AT(IJK)) CYCLE IJK_LP
394     ! Calculating an average specific heat for the fluid.
395                 C_PS(IJK, M) = ZERO
396     
397                 DO N = 1, NMAX(M)
398                    lCp = calc_CpoR(T_S(IJK,M), M, N, lCP_Err)
399                    C_PS(IJK,M) = C_PS(IJK,M) + X_s(IJK,M,N) * &
400                       (lCp * RGAS / MW_s(M,N))
401                 ENDDO
402     
403              ENDDO IJK_LP
404           ENDDO M_LP
405     
406     ! The database calculation always returns cal/g.K thus the following
407     ! conversion is needed if using SI units. ** Vector operation
408           IF(UNITS == 'SI') C_PS = 4.183925d3 * C_PS  !in J/kg K
409     
410           END SUBROUTINE PHYSICAL_PROP_CPs
411     
412     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
413     !                                                                      C
414     !  Subroutine: PHYSICAL_PROP_CPs                                       C
415     !  Purpose: Calculate solids phase constant pressure specific heat.    C
416     !                                                                      C
417     !  Author: J. Musser                                  Date: 28-JUN-13  C
418     !  Reviewer:                                          Date:            C
419     !                                                                      C
420     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
421           SUBROUTINE PHYSICAL_PROP_Dp
422     
423           use run, only: CALL_DQMOM
424     
425           use fldvar, only: scalar
426           use scalars, only: phase4scalar
427     
428           use fldvar, only: ROP_s
429           use fldvar, only: RO_s
430           use fldvar, only: D_p, EP_S
431           use functions
432     
433           implicit none
434     
435     ! Loop indicies
436           INTEGER :: IJK   ! Computational cell
437           INTEGER :: M     ! Solids phase
438     
439     ! Map from true index to map.
440           INTEGER :: lM
441     
442           IF(.NOT.CALL_DQMOM) return
443     
444           M_LP: DO M=1, MMAX
445     
446              lM = phase4scalar(M) ! Map from scalar eq to solids phase
447     
448              IF(.NOT.PSIZE(M)) CYCLE M_LP
449              IJK_LP: DO IJK = IJKSTART3, IJKEND3
450                 IF(WALL_AT(IJK)) CYCLE IJK_LP
451     
452                 IF(EP_s(IJK,lM) > small_number) D_p(IJK,M)= Scalar(IJK,lM)
453     
454              ENDDO IJK_LP
455           ENDDO M_LP
456     
457           return
458           END SUBROUTINE PHYSICAL_PROP_Dp
459     
460     
461     
462     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
463     !                                                                      C
464     !  Subroutine: NegROg_LOG                                              C
465     !  Purpose: Record information about the location and conditions that  C
466     !           resulted in a negative gas phase density.                  C
467     !                                                                      C
468     !  Author: J. Musser                                  Date: 28-JUN-13  C
469     !  Reviewer:                                          Date:            C
470     !                                                                      C
471     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
472           SUBROUTINE ROgErr_LOG(IJK, tHeader)
473     
474     ! Simulation time
475           use run, only: TIME
476     ! Gas phase temperature.
477           use fldvar, only: T_g
478     ! Gas phase density (compressible).
479           use fldvar, only: RO_g
480     ! Gas phase pressure.
481           use fldvar, only: P_g
482           use cutcell
483     
484           INTEGER, intent(in) :: IJK
485           LOGICAL, intent(inout) :: tHeader
486     
487           LOGICAL :: lExists
488           CHARACTER(LEN=32) :: lFile
489           INTEGER, parameter :: lUnit = 4868
490           LOGICAL, save :: fHeader = .TRUE.
491     
492     
493           lFile = '';
494           if(numPEs > 1) then
495              write(lFile,"('ROgErr_',I4.4,'.log')") myPE
496           else
497              write(lFile,"('ROgErr.log')")
498           endif
499           inquire(file=trim(lFile),exist=lExists)
500           if(lExists) then
501              open(lUnit,file=trim(adjustl(lFile)),                         &
502                 status='old', position='append')
503           else
504              open(lUnit,file=trim(adjustl(lFile)), status='new')
505           endif
506     
507           if(fHeader) then
508              write(lUnit,1000)
509              fHeader = .FALSE.
510           endif
511     
512           if(tHeader) then
513              write(lUnit,"(/2x,'Simulation time: ',g12.5)") TIME
514              tHeader = .FALSE.
515           endif
516     
517           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
518           write(lUnit,"(6x,A,1X,g12.5)",ADVANCE='NO') 'RO_g:', RO_g(IJK)
519           write(lUnit,"(2x,A,1X,g12.5)",ADVANCE='NO') 'P_g:', P_g(IJK)
520           write(lUnit,"(2x,A,1X,g12.5)") 'T_g:', T_g(IJK)
521           if(CARTESIAN_GRID) then
522              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
523              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
524              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
525                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
526           endif
527     
528           close(lUnit)
529     
530           RETURN
531     
532      1000 FORMAT(2X,'One or more cells have reported a negative gas',      &
533              ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
534              ' lower UR_FAC(1) in mfix.dat.')
535     
536      1001 FORMAT(/4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
537     
538           END SUBROUTINE ROgErr_LOG
539     
540     
541     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
542     !                                                                      !
543     !  Subroutine: NegROs_LOG                                              !
544     !  Purpose: Record information about the location and conditions that  !
545     !           resulted in a negative solids phase density.               !
546     !                                                                      !
547     !  Author: J. Musser                                  Date: 09-Oct-13  !
548     !  Reviewer:                                          Date:            !
549     !                                                                      !
550     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
551           SUBROUTINE ROsErr_LOG(IJK, M, tHeader)
552     
553     ! Global variables:
554     !-----------------------------------------------------------------------
555     ! Simulation time
556           use run, only: TIME
557     ! Solid phase species mass fractions.
558           use fldvar, only: X_s
559     ! Solid phase density (variable).
560           use fldvar, only: RO_s
561     ! Baseline/Unreaced solids density
562           use physprop, only: BASE_ROs
563     ! Initial mass fraction of inert species
564           use physprop, only: X_S0
565     ! Index of inert solids phase species.
566           use physprop, only: INERT_SPECIES
567     
568     ! Full Access to cutcell data.
569           use cutcell
570     
571     
572     ! Local Variables:
573     !-----------------------------------------------------------------------
574     ! Fluid cell and solids phase indices
575           INTEGER, intent(in) :: IJK, M
576     ! Flag to output header
577           LOGICAL, intent(inout) :: tHeader
578     ! Local aliase for inert species index
579           INTEGER :: N
580     ! Local file values.
581           LOGICAL :: lExists
582           CHARACTER(LEN=32) :: lFile
583           INTEGER, parameter :: lUnit = 4868
584           LOGICAL, save :: fHeader = .TRUE.
585     
586     
587           lFile = '';
588           if(numPEs > 1) then
589              write(lFile,"('ROsErr_',I4.4,'.log')") myPE
590           else
591              write(lFile,"('ROsErr.log')")
592           endif
593           inquire(file=trim(lFile),exist=lExists)
594           if(lExists) then
595              open(lUnit,file=trim(adjustl(lFile)),                         &
596                 status='old', position='append')
597           else
598              open(lUnit,file=trim(adjustl(lFile)), status='new')
599           endif
600     
601           if(fHeader) then
602              write(lUnit,1000)
603              fHeader = .FALSE.
604           endif
605     
606           if(tHeader) then
607              write(lUnit,"(/2x,'Simulation time: ',g12.5,'  Phase: ',I2)")&
608                  TIME, M
609              tHeader = .FALSE.
610           endif
611     
612           N = INERT_SPECIES(M)
613           write(lUnit,1001) IJK, I_OF(IJK), J_OF(IJK), K_OF(IJK)
614           write(lUnit,"(6x,A,1X,g12.5)",advance='no') 'RO_s:', RO_s(IJK,M)
615           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'Base:', BASE_ROs(M)
616           write(lUnit,"(2x,A,1X,g12.5)",advance='no') 'X_s0:', X_s0(M,N)
617           write(lUnit,"(2x,A,1X,g12.5)") 'X_s:', X_s(IJK,M,N)
618     
619           if(CARTESIAN_GRID) then
620              write(lUnit,"(6x,A,1X,L1)",ADVANCE='NO') 'Cut Cell:', CUT_CELL_AT(IJK)
621              write(lUnit,"(6x,A,1X,L1)") 'Small Cell:', SMALL_CELL_AT(IJK)
622              write(lUnit,"(6x,'Coordinates (E/N/T): ',1X,3(2x, g17.8))") &
623                 xg_e(I_OF(IJK)), yg_n(J_of(ijk)), zg_t(k_of(ijk))
624           endif
625     
626           close(lUnit)
627     
628           RETURN
629     
630      1000 FORMAT(2X,'One or more cells have reported a negative gas',      &
631              ' density (RO_g(IJK)). If',/2x,'this is a persistent issue,', &
632              ' lower UR_FAC(1) in mfix.dat.')
633     
634      1001 FORMAT( 4X,'IJK: ',I8,7X,'I: ',I4,'  J: ',I4,'  K: ',I4)
635     
636           END SUBROUTINE ROsErr_LOG
637     
638           END SUBROUTINE PHYSICAL_PROP
639