File: N:\mfix\model\des\write_des_data.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
2     !                                                                      !
3     !  Subroutine: WRITE_DES_DATA                                          !
4     !  Purpose: Writing DES output in Paraview format                      !
5     !                                                                      !
6     !                                                                      !
7     !  Author: Jay Boyalakuntla                           Date: 26-Jul-06  !
8     !  Reviewer: Sreekanth Pannala                        Date: 31-Oct-06  !
9     !                                                                      !
10     !  Reviewer: J. Musser                                Date: 20-Apr-10  !
11     !  Comments: Split original subroutine into one for ParaView *.vtp     !
12     !  files, and a second for TECPLOT files *.dat.                        !
13     !                                                                      !
14     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
15     
16           SUBROUTINE WRITE_DES_DATA
17     
18     !-----------------------------------------------
19     ! Modules
20     !-----------------------------------------------
21           USE param
22           USE param1
23           USE parallel
24           USE fldvar
25           USE discretelement
26           USE run
27           USE geometry
28           USE physprop
29           USE sendrecv
30           USE des_bc
31     
32           use error_manager
33     
34           IMPLICIT NONE
35     !-----------------------------------------------
36     ! Local Variables
37     !-----------------------------------------------
38     
39     !-----------------------------------------------
40     ! Functions
41     !-----------------------------------------------
42     
43     !-----------------------------------------------
44     
45           IF (TRIM(DES_OUTPUT_TYPE) .EQ. 'TECPLOT') THEN
46              CALL WRITE_DES_TECPLOT
47           ELSE
48              CALL WRITE_DES_VTP
49           ENDIF
50     
51     ! Invoke at own risk
52     
53     ! Granular temperature subroutine should be called/calculated when
54     ! writing DES data
55     
56           IF (DES_CALC_BEDHEIGHT) THEN
57              CALL CALC_DES_BEDHEIGHT
58              CALL WRITE_DES_BEDHEIGHT
59           ENDIF
60     
61           IF (.FALSE.) CALL DES_GRANULAR_TEMPERATURE
62           IF (.FALSE.) CALL WRITE_DES_THETA
63     
64           RETURN
65           END SUBROUTINE WRITE_DES_DATA
66     !-----------------------------------------------
67     
68     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
69     !                                                                      !
70     !  Subroutine: WRITE_DES_VTP                                           !
71     !  Purpose: Writing DES output in Paraview format                      !
72     !                                                                      !
73     !                                                                      !
74     !  Reviewer: Rahul Garg                               Dare: 01-Aug-07  !
75     !  Comments: Added one more output file containing averaged bed height !
76     !                                                                      !
77     !  Revision : For parallel runs added cumulative and parallel IO       !
78     !  Author   : Pradeep G.                                               !
79     !                                                                      !
80     !  NOTE: If the system starts with zero particles, ParaView may have   !
81     !  trouble showing the results. To view the results in the current     !
82     !  version of ParaView, Version 3.6.1:                                 !
83     !    i - load the *.vtp files                                          !
84     !   ii - filter with glyph (Filters > Common > Glyph)                  !
85     !        a - change glyph to sphere                                    !
86     !        b - change scale mode to scalar                               !
87     !        c - check the "Edit" Set Scale Factor Box                     !
88     !        d - change the value to 1.0                                   !
89     !                                                                      !
90     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
91     
92           SUBROUTINE WRITE_DES_VTP
93     
94           use des_rxns, only: DES_X_s
95           use des_thermo, only: DES_T_s
96           use discretelement, only: DES_POS_NEW, DES_VEL_NEW, DES_USR_VAR
97           use discretelement, only: DES_RADIUS
98           use discretelement, only: DES_USR_VAR, DES_USR_VAR_SIZE
99           use discretelement, only: S_TIME
100           use discretelement, only: USE_COHESION, PostCohesive
101           use error_manager
102           use mfix_pic, only: des_stat_wt, mppic
103           use param
104           use param1, only: zero
105           use run, only: ANY_SPECIES_EQ
106           use run, only: ENERGY_EQ
107           use vtp
108     
109           IMPLICIT NONE
110     
111           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP_PAR
112           CHARACTER(len=10) :: lNoP
113           CHARACTER(len=24) :: sTIMEc
114           INTEGER :: NN
115     
116           sTIMEc=''; WRITE(sTIMEc,"(ES24.16)") S_TIME
117     
118     ! This routine opens the VTP file and calculates send/recv information.
119     ! It returns back the number of points as a string.
120           CALL VTP_OPEN_FILE(lNoP)
121     
122     ! Standard VTP header information:
123     !----------------------------------------------------------------------/
124           CALL VTP_WRITE_ELEMENT('<?xml version="1.0"?>')
125           CALL VTP_WRITE_ELEMENT('<!-- Time ='//sTIMEc//'s -->')
126           CALL VTP_WRITE_ELEMENT('<VTKFile type="PolyData" version="0.1" &
127              &byte_order="LittleEndian">')
128           CALL VTP_WRITE_ELEMENT('<PolyData>')
129     
130           CALL VTP_WRITE_ELEMENT('<Piece NumberOfPoints="'//lNoP//'" &
131              &NumberOfVerts="0" NumberOfLines="0" NumberOfStrips="0" &
132              &NumberOfPolys="0">')
133     
134     ! Points are the particle identified by position:
135     !----------------------------------------------------------------------/
136           CALL VTP_WRITE_ELEMENT('<Points>')
137           CALL VTP_WRITE_DATA('Position', DES_POS_NEW)
138           CALL VTP_WRITE_ELEMENT('</Points>')
139     
140     ! PointData are individual particle properties:
141     !----------------------------------------------------------------------/
142           ALLOCATE(TMP_PAR(SIZE(DES_RADIUS)))
143     
144           CALL VTP_WRITE_ELEMENT('<PointData Scalars="Diameter" &
145              &Vectors="Velocity">')
146     
147           CALL GET_DIAMETER(TMP_PAR)
148           CALL VTP_WRITE_DATA('Diameter', TMP_PAR)
149     
150           CALL VTP_WRITE_DATA('Velocity', DES_VEL_NEW)
151     
152           CALL VTP_WRITE_DATA('ID', iGLOBAL_ID)
153     
154           DO NN=1, DES_USR_VAR_SIZE
155              CALL VTP_WRITE_DATA(trim(iVar('USR_VAR',NN)),DES_USR_VAR(NN,:))
156           ENDDO
157     
158           IF(PARTICLE_ORIENTATION) &
159              CALL VTP_WRITE_DATA('Orientation', ORIENTATION)
160     
161           IF(MPPIC) THEN
162              TMP_PAR = 1.0d0 - EPG_P
163              CALL VTP_WRITE_DATA('EPs', TMP_PAR)
164           ENDIF
165     
166           IF(ENERGY_EQ) &
167              CALL VTP_WRITE_DATA('Temperature', DES_T_s)
168     
169           IF(ANY_SPECIES_EQ) THEN
170              DO NN=1, DIMENSION_N_S
171                 CALL VTP_WRITE_DATA(trim(iVar('X_s',NN)), DES_X_s(:,NN))
172              ENDDO
173           ENDIF
174     
175           IF(USE_COHESION) &
176              CALL VTP_WRITE_DATA('CohesiveForce', PostCohesive)
177     
178           CALL VTP_WRITE_ELEMENT('</PointData>')
179           DEALLOCATE(TMP_PAR)
180     
181     ! Open/Close the unused VTP tags.
182     !----------------------------------------------------------------------/
183           CALL VTP_WRITE_ELEMENT('<CellData></CellData>')
184           CALL VTP_WRITE_ELEMENT('<Verts></Verts>')
185           CALL VTP_WRITE_ELEMENT('<Lines></Lines>')
186           CALL VTP_WRITE_ELEMENT('<Strips></Strips>')
187           CALL VTP_WRITE_ELEMENT('<Polys></Polys>')
188     
189     ! Close all the opened tags:
190     !----------------------------------------------------------------------/
191           CALL VTP_WRITE_ELEMENT('</Piece>')
192           CALL VTP_WRITE_ELEMENT('</PolyData>')
193           CALL VTP_WRITE_ELEMENT('</VTKFile>')
194     
195           CALL VTP_CLOSE_FILE
196     
197     ! Add the new VTP file to the PVD file for time association.
198           CALL ADD_VTP_TO_PVD
199     
200           RETURN
201     
202           contains
203     
204     !......................................................................!
205     ! Purpose: Return the diamerter of DES particles.                      !
206     !......................................................................!
207           SUBROUTINE GET_DIAMETER(OUT_VAR)
208     
209           IMPLICIT NONE
210     
211           DOUBLE PRECISION, INTENT(OUT) :: OUT_VAR(:)
212           INTEGER :: LC
213     
214     ! Return the effective diamter of the parcle.
215           IF(MPPIC) THEN
216              DO LC=1,size(OUT_VAR)
217                 IF(DES_STAT_WT(LC) > ZERO) THEN
218                    OUT_VAR(LC) = 2.0d0*DES_RADIUS(LC)*&
219                       (DES_STAT_WT(LC)**(1./3.))
220                 ELSE
221                    OUT_VAR(LC) = 2.0d0*DES_RADIUS(LC)
222                 ENDIF
223              ENDDO
224           ELSE
225              OUT_VAR = 2.0d0*DES_RADIUS
226           ENDIF
227     
228           RETURN
229           END SUBROUTINE GET_DIAMETER
230     
231           END SUBROUTINE WRITE_DES_VTP
232     
233     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
234     !
235     !  Module name: WRITE_DES_TECPLOT
236     !  Purpose: Writing DES output in TECPLOT format
237     !
238     !  Revision: For parallel runs added distributed and single IO
239     !  Comment: In earlier version the time instances are keep appended to
240     !           the tecplot file. This will make tecplot file so large for
241     !           large simulations. Hence seperate files are written for
242     !           each instances
243     !  Author : Pradeep G.
244     !
245     !
246     !
247     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
248           SUBROUTINE WRITE_DES_TECPLOT
249     
250           USE param
251           USE param1
252           USE parallel
253           USE fldvar
254           USE discretelement
255           USE run
256           USE geometry
257           USE physprop
258           USE sendrecv
259           USE des_bc
260           USE compar
261           USE cdist
262           USE desmpi
263           USE mpi_comm_des
264           USE mpi_utility
265           USE functions
266           IMPLICIT NONE
267     
268     !-----------------------------------------------
269     ! Local Variables
270     !-----------------------------------------------
271     ! file units for
272     ! Pradeep remove parameter from following variables
273           INTEGER:: DES_DATA,DES_EX,DES_EPS
274     
275     ! output file for basic DES variables including: position, velocity,
276     ! radius, density, mark (flag)
277           CHARACTER(LEN=50)     :: FNAME_DATA
278     
279     ! output file for extra DES variables including:
280     ! solids time (S_TIME), maximum neighbor count, maximum overlap
281     ! granular energy and granular temperature
282           CHARACTER(LEN=50)     :: FNAME_EXTRA
283     
284     ! output file for axial solids volume fraction and granular temp
285           CHARACTER(LEN=50)     :: FNAME_EPS
286     
287     ! dummy indices
288           INTEGER L, K
289     
290     ! index to track accounted for particles
291           INTEGER PC
292     
293     ! Variables related to gathering info at PE_IO
294           integer llocalcnt,lglocnt,lgathercnts(0:numpes-1),lproc,ltotvar,lcount
295           double precision, dimension(:,:), allocatable :: ltemp_array
296     
297           INTEGER :: wDIMN
298     
299     ! Set output dimnensions
300           wDIMN = merge(2,3,NO_K)
301     ! set the total variable based on dimension
302           ltotvar = merge(8,9,NO_K)
303     
304     ! set the file name and unit number and open file
305           des_data = 2000
306           des_ex = 2100
307           des_eps = 2200
308           if (bdist_io) then
309              write(fname_data,'(A,"_DES_DATA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
310              write(fname_extra,'(A,"_DES_EXTRA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
311              write(fname_eps,'(A,"_DES_EPS",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
312              open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
313              open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
314              open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
315           else
316              if(mype.eq.pe_io) then
317                 write(fname_data,'(A,"_DES_DATA_",I4.4,".dat")') trim(run_name),tecplot_findex
318                 write(fname_extra,'(A,"_DES_EXTRA_",I4.4,".dat")') trim(run_name),tecplot_findex
319                 write(fname_eps,'(A,"_DES_EPS_",I4.4,".dat")') trim(run_name),tecplot_findex
320                 open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
321                 open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
322                 open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
323              end if
324           end if
325           tecplot_findex = tecplot_findex + 1
326     
327     ! write header
328           if (bdist_io .or. mype .eq. pe_io) then
329              if(DO_K) then
330                 write (des_data,'(A)') &
331                 'variables = "x" "y" "z" "vx" "vy" "vz" "rad" "den" "mark"'
332              else
333                 write (des_data, '(A)') &
334                 'variables = "x" "y" "vx" "vy" "omega" "rad" "den" "mark"'
335              endif
336              write (des_data, "(A,F15.7,A)")'zone T ="', s_time, '"'
337           end if
338     
339     ! write data in ditributed mode or as a single file
340           if(bdist_io) then
341              pc = 1
342              do l = 1,max_pip
343                 if(pc.gt.pip) exit
344                 if(is_nonexistent(l)) cycle
345                 pc = pc+1
346                 if(is_ghost(l) .or. is_entering_ghost(l) .or. is_exiting_ghost(l)) cycle
347                 if(DO_K) then
348                    write (des_data, '(8(2x,es12.5))')&
349                       (des_pos_new(l,k),k=1,wDIMN),(des_vel_new(l,k),k=1,wDIMN), &
350                        des_radius(l),ro_sol(l)
351                 else
352                    write (des_data, '(7(2x,es12.5))')&
353                       (des_pos_new(l,k),k=1,wDIMN), (des_vel_new(l,k),k=1,wDIMN), &
354                       omega_new(l,1), des_radius(l), ro_sol(l)
355                 endif
356             end do
357     
358           else ! if bdist_io
359     ! set parameters required for gathering info at PEIO and write as single file
360              lglocnt = 10
361              llocalcnt = pip - ighost_cnt
362              call global_sum(llocalcnt,lglocnt)
363              allocate (dprocbuf(llocalcnt),drootbuf(lglocnt),iprocbuf(llocalcnt),irootbuf(lglocnt))
364              allocate (ltemp_array(lglocnt,ltotvar))
365              igath_sendcnt = llocalcnt
366              lgathercnts = 0
367              lgathercnts(mype) = llocalcnt
368              call global_sum(lgathercnts,igathercnts)
369              idispls(0) = 0
370              do lproc = 1,numpes-1
371                 idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
372              end do
373     
374     ! gather information from all processor
375              lcount = 1
376              do k = 1,wDIMN
377                 call des_gather(des_pos_new(:,k))
378                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
379              end do
380              do k = 1,wDIMN
381                 call des_gather(des_vel_new(:,k))
382                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
383              end do
384              if(NO_K) then
385                 call des_gather(omega_new(:,1))
386                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
387              end if
388              call des_gather(des_radius)
389              ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
390              call des_gather(ro_sol)
391              ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
392     
393     ! write the data into file
394              if (mype.eq.pe_io) then
395                 if(DO_K) then
396                    do l =1,lglocnt
397                       write (des_data,'(8(2x,es12.5),I5)') (ltemp_array(l,k),k=1,8),int(ltemp_array(l,9))
398                    end do
399                 else
400                    do l =1,lglocnt
401                       write (des_data,'(7(2x,es12.5),I5)') (ltemp_array(l,k),k=1,7),int(ltemp_array(l,8))
402                    end do
403                 end if
404              end if
405              deallocate (dprocbuf,drootbuf,iprocbuf,irootbuf,ltemp_array)
406           end if
407     ! close the files
408           if (mype.eq.pe_io .or. bdist_io) then
409              close(des_data)
410              close(des_ex)
411              close(des_eps)
412           end if
413           return
414     
415       999 write(*,"(/1x,70('*'),//,a,/,a,/1x,70('*'))")&
416              ' From: write_des_tecplot ',&
417              ' message: error opening des tecplot file. terminating run.'
418     
419     
420           END SUBROUTINE WRITE_DES_TECPLOT
421     !-----------------------------------------------
422     
423     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
424     !
425     !  Module name: WRITE_DES_BEDHEIGHT
426     !  Purpose: Writing DES output on bed height.
427     
428     !  WARNING: This code is out-of-date and should be modified for consistency
429     !  with current DEM version.  Also this routine will be fairly specific
430     !  to a user needs and should probably be tailored as such
431     
432     !
433     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
434           SUBROUTINE WRITE_DES_BEDHEIGHT
435     
436           USE param
437           USE param1
438           USE parallel
439           USE fldvar
440           USE discretelement
441           USE run
442           USE geometry
443           USE physprop
444           USE sendrecv
445           USE des_bc
446           IMPLICIT NONE
447     
448     !-----------------------------------------------
449     ! Local Variables
450     !-----------------------------------------------
451     ! logical that identifies that the data file has been created
452     ! and is already opened (initial checks/writes)
453           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
454     ! logical used for testing is the data files already exists
455           LOGICAL :: F_EXISTS
456     ! output file for the bed height data
457           CHARACTER(LEN=50)     :: FNAME_BH
458     ! file unit for the bed height data
459           INTEGER, PARAMETER :: BH_UNIT = 2010
460     ! dummy index values
461           INTEGER I, M
462     ! variables for bed height calculation
463           INTEGER, SAVE :: tcount = 1
464           DOUBLE PRECISION :: height_avg, height_rms
465           DOUBLE PRECISION, PARAMETER :: tmin = 5.d0
466           DOUBLE PRECISION, DIMENSION(5000), SAVE :: bed_height_time
467     
468     !-----------------------------------------------
469     ! Functions
470     !-----------------------------------------------
471     !-----------------------------------------------
472     
473     ! after tmin start storing bed height. after enough measurements
474     ! have been taken (i.e. tcount > 20) start to calculate a running
475     ! average bed height and running rms bed height for solids phase 1 only
476           height_avg = zero
477           height_rms = zero
478     
479           if(time.gt.tmin) then
480              if(tcount.le.5000)  then
481                 bed_height_time(tcount) = bed_height(1)
482                 !dt_time(tcount) = DT
483                 tcount = tcount + 1
484     
485                 if(tcount.gt.20)  then
486                    do i = 1, tcount-1,1
487                       height_avg = height_avg + bed_height_time(i)!*dt_time(i)
488                    enddo
489                    height_avg = height_avg/(tcount-1)
490                    do i = 1, tcount-1,1
491                       height_rms = height_rms + ((bed_height_time(i)&
492                            &-height_avg)**2)!*dt_time(i)
493                    enddo
494     
495                    height_rms = sqrt(height_rms/(tcount-1))
496                 endif
497              endif
498           endif
499     
500           FNAME_BH = TRIM(RUN_NAME)//'_DES_BEDHEIGHT.dat'
501           IF(FIRST_PASS) THEN
502              F_EXISTS = .FALSE.
503              INQUIRE(FILE=FNAME_BH,EXIST=F_EXISTS)
504     ! If the file does not exist, then create it with the necessary
505     ! header information.
506              IF (.NOT.F_EXISTS) THEN
507                 OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,&
508                    FORM="formatted",STATUS="new")
509              ELSE
510     ! To prevent overwriting existing files accidently, exit if the file
511     ! exists and this is a NEW run.
512                 IF(RUN_TYPE .EQ. 'NEW') THEN
513                    WRITE(*,1000)
514                    WRITE(UNIT_LOG, 1000)
515                    CALL MFIX_EXIT(myPE)
516                 ELSE
517     ! Open the file for appending of new data (RESTART_1 Case)
518                    OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,&
519                         POSITION="append")
520                 ENDIF
521              ENDIF
522              FIRST_PASS = .FALSE.
523           ELSE
524     ! Open the file and mark for appending
525              OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,&
526                   POSITION="append")
527           ENDIF
528     
529           WRITE(BH_UNIT, '(10(2X,E20.12))') s_time, &
530              (bed_height(M), M=MMAX+1,MMAX+DES_MMAX), height_avg, height_rms
531     ! Close the file and keep
532           CLOSE(BH_UNIT, STATUS="KEEP")
533     
534           RETURN
535     
536      1000 FORMAT(/1X,70('*')//, ' From: WRITE_DES_BEDHEIGHT',/,&
537              ' Message: bed_height.dat already exists in the run',&
538              ' directory.',/10X, 'Run terminated to prevent',&
539              ' accidental overwriting of files.',/1X,70('*')/)
540     
541           END SUBROUTINE WRITE_DES_BEDHEIGHT
542     
543     
544     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
545     !
546     !  Module name: WRITE_DES_THETA
547     !  Purpose: The following code writes out theta_m for discrete
548     !  particles to a file for each ijk cell in the system each time
549     !  des_granular_temperature is called.
550     !
551     !
552     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
553           SUBROUTINE WRITE_DES_THETA
554     
555           USE param
556           USE param1
557           USE parallel
558           USE fldvar
559           USE discretelement
560           USE run
561           USE geometry
562           USE physprop
563           USE sendrecv
564           USE des_bc
565           USE functions
566           IMPLICIT NONE
567     
568     !-----------------------------------------------
569     ! Local Variables
570     !-----------------------------------------------
571     ! indices
572           INTEGER I, J, K, IJK
573     !
574           INTEGER M, NP
575     ! logical that identifies that the data file has been created
576     ! and is already opened (initial checks/writes)
577           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
578     ! logical used for testing is the data file already exists
579           LOGICAL :: F_EXISTS
580     ! file unit for the granular temperature data
581           INTEGER, PARAMETER :: GT_UNIT = 2020
582     ! output file for the granular temperature data
583           CHARACTER(LEN=50)  :: FNAME_GT
584     !-----------------------------------------------
585     
586           FNAME_GT = TRIM(RUN_NAME)//'_DES_THETA.dat'
587           IF (FIRST_PASS) THEN
588              F_EXISTS = .FALSE.
589              INQUIRE(FILE=FNAME_GT,EXIST=F_EXISTS)
590     
591              IF (.NOT.F_EXISTS) THEN
592     ! If the file does not exist, then create it with the necessary
593     ! header information.
594                 OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,&
595                      STATUS='NEW')
596              ELSE
597                 IF(RUN_TYPE .EQ. 'NEW') THEN
598     ! If the run is new and the GT file already exists replace it with a
599     ! new file.
600     !               OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,STATUS='REPLACE')
601     ! Prevent overwriting an existing file by exiting if the file exists
602     ! and this is a NEW run.
603                    WRITE(*,1001) FNAME_GT
604                    WRITE(UNIT_LOG,1001) FNAME_GT
605                    CALL MFIX_EXIT(myPE)
606                 ELSE
607     ! Open the file for appending of new data (RESTART_1 Case)
608                    OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT, FILE=FNAME_GT,&
609                         POSITION='APPEND')
610                 ENDIF
611              ENDIF
612              FIRST_PASS =  .FALSE.
613           ELSE
614     ! Open file and mark for appending
615              OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,&
616                   POSITION='APPEND')
617           ENDIF   ! endif (first_pass)
618     
619           WRITE(GT_UNIT,*) ''
620           WRITE(GT_UNIT,'(A6,ES24.16)') 'Time=', S_TIME
621           WRITE(GT_UNIT,'(A6,2X,3(A6,2X),A8)',ADVANCE="NO") 'IJK', &
622              'I', 'J', 'K', 'NP'
623           DO M = MMAX+1,DES_MMAX+MMAX
624              WRITE(GT_UNIT,'(7X,A6,I1)',ADVANCE="NO") 'THETA_',M
625           ENDDO
626           WRITE(GT_UNIT,*) ''
627           DO IJK = IJKSTART3, IJKEND3
628              IF(FLUID_AT(IJK)) THEN
629                 I = I_OF(IJK)
630                 J = J_OF(IJK)
631                 K = K_OF(IJK)
632                 NP = PINC(IJK)
633                 WRITE(GT_UNIT,'(I6,2X,3(I6,2X),I8,(2X,ES15.5))') &
634                    IJK, I, J, K, NP, (THETA_M(IJK,M), M = MMAX+1,DES_MMAX+MMAX)
635              ENDIF
636           ENDDO
637     
638     ! Close the file and keep
639           CLOSE(GT_UNIT, STATUS='KEEP')
640     
641           RETURN
642     
643      1001 FORMAT(/1X,70('*')//, ' From: WRITE_DES_THETA',/,&
644              ' Message: ', A, ' already exists in the run',/10X,&
645              'directory. Run terminated to prevent accidental overwriting',&
646              /10X,'of files.',/1X,70('*')/)
647     
648           END SUBROUTINE WRITE_DES_THETA
649     !-----------------------------------------------
650