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