File: RELATIVE:/../../../mfix.git/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 vtp
95           use run, only: ENERGY_EQ
96           use discretelement, only: S_TIME
97           use discretelement, only: DES_POS_NEW, DES_VEL_NEW, DES_USR_VAR
98           use discretelement, only: DES_USR_VAR, DES_USR_VAR_SIZE
99           use des_thermo, only: DES_T_s_NEW
100           use discretelement, only: DES_RADIUS
101           use discretelement, only: USE_COHESION, PostCohesive
102           use des_rxns, only: DES_X_s
103           use run, only: ANY_SPECIES_EQ
104           use param, only: DIMENSION_N_S
105           USE mfix_pic, only: des_stat_wt, mppic
106     
107           use error_manager
108     
109           IMPLICIT NONE
110     
111           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: DES_DIAMETER !(PARTICLES)
112           CHARACTER(len=10) :: lNoP
113           CHARACTER(len=24) :: sTIMEc
114           INTEGER :: N
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           CALL VTP_WRITE_ELEMENT('<PointData Scalars="Diameter" &
143              &Vectors="Velocity">')
144     
145           ALLOCATE(DES_DIAMETER(SIZE(DES_RADIUS)))
146           DES_DIAMETER(:) = 2.0d0*DES_RADIUS(:)
147           CALL VTP_WRITE_DATA('Diameter', DES_DIAMETER)
148           DEALLOCATE(DES_DIAMETER)
149     
150           CALL VTP_WRITE_DATA('Velocity', DES_VEL_NEW)
151     
152           IF(DES_USR_VAR_SIZE > 0) &
153              CALL VTP_WRITE_DATA('User Defined Var', DES_USR_VAR)
154     
155           IF(PARTICLE_ORIENTATION) &
156              CALL VTP_WRITE_DATA('Orientation', ORIENTATION)
157     
158     !      IF(MPPIC) CALL VTP_WRITE_DATA('Statwt', DES_STAT_WT)
159     
160           IF(ENERGY_EQ) &
161              CALL VTP_WRITE_DATA('Temperature', DES_T_s_NEW)
162     
163           IF(ANY_SPECIES_EQ) THEN
164              DO N=1, DIMENSION_N_S
165                 CALL VTP_WRITE_DATA(trim(iVar('X_s',N)), DES_X_s(:,N))
166              ENDDO
167           ENDIF
168     
169           IF(USE_COHESION) &
170              CALL VTP_WRITE_DATA('CohesiveForce', PostCohesive)
171     
172           CALL VTP_WRITE_ELEMENT('</PointData>')
173     
174     ! Open/Close the unused VTP tags.
175     !----------------------------------------------------------------------/
176           CALL VTP_WRITE_ELEMENT('<CellData></CellData>')
177           CALL VTP_WRITE_ELEMENT('<Verts></Verts>')
178           CALL VTP_WRITE_ELEMENT('<Lines></Lines>')
179           CALL VTP_WRITE_ELEMENT('<Strips></Strips>')
180           CALL VTP_WRITE_ELEMENT('<Polys></Polys>')
181     
182     ! Close all the opened tags:
183     !----------------------------------------------------------------------/
184           CALL VTP_WRITE_ELEMENT('</Piece>')
185           CALL VTP_WRITE_ELEMENT('</PolyData>')
186           CALL VTP_WRITE_ELEMENT('</VTKFile>')
187     
188           CALL VTP_CLOSE_FILE
189     
190     ! Add the new VTP file to the PVD file for time association.
191           CALL ADD_VTP_TO_PVD
192     
193           RETURN
194           END SUBROUTINE WRITE_DES_VTP
195     
196     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
197     !
198     !  Module name: WRITE_DES_TECPLOT
199     !  Purpose: Writing DES output in TECPLOT format
200     !
201     !  Revision: For parallel runs added distributed and single IO
202     !  Comment: In earlier version the time instances are keep appended to
203     !           the tecplot file. This will make tecplot file so large for
204     !           large simulations. Hence seperate files are written for
205     !           each instances
206     !  Author : Pradeep G.
207     !
208     !
209     !
210     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
211           SUBROUTINE WRITE_DES_TECPLOT
212     
213           USE param
214           USE param1
215           USE parallel
216           USE fldvar
217           USE discretelement
218           USE run
219           USE geometry
220           USE physprop
221           USE sendrecv
222           USE des_bc
223           USE compar
224           USE cdist
225           USE desmpi
226           USE mpi_comm_des
227           USE mpi_utility
228           USE functions
229           IMPLICIT NONE
230     
231     !-----------------------------------------------
232     ! Local Variables
233     !-----------------------------------------------
234     ! file units for
235     ! Pradeep remove parameter from following variables
236           INTEGER:: DES_DATA,DES_EX,DES_EPS
237     
238     ! output file for basic DES variables including: position, velocity,
239     ! radius, density, mark (flag)
240           CHARACTER(LEN=50)     :: FNAME_DATA
241     
242     ! output file for extra DES variables including:
243     ! solids time (S_TIME), maximum neighbor count, maximum overlap
244     ! granular energy and granular temperature
245           CHARACTER(LEN=50)     :: FNAME_EXTRA
246     
247     ! output file for axial solids volume fraction and granular temp
248           CHARACTER(LEN=50)     :: FNAME_EPS
249     
250     ! dummy indices
251           INTEGER L, K
252     
253     ! index to track accounted for particles
254           INTEGER PC
255     
256     ! Variables related to gathering info at PE_IO
257           integer llocalcnt,lglocnt,lgathercnts(0:numpes-1),lproc,ltotvar,lcount
258           double precision, dimension(:,:), allocatable :: ltemp_array
259     
260           INTEGER :: wDIMN
261     
262     ! Set output dimnensions
263           wDIMN = merge(2,3,NO_K)
264     ! set the total variable based on dimension
265           ltotvar = merge(8,9,NO_K)
266     
267     ! set the file name and unit number and open file
268           des_data = 2000
269           des_ex = 2100
270           des_eps = 2200
271           if (bdist_io) then
272              write(fname_data,'(A,"_DES_DATA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
273              write(fname_extra,'(A,"_DES_EXTRA",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
274              write(fname_eps,'(A,"_DES_EPS",I4.4,"_",I4.4,".dat")') trim(run_name),tecplot_findex,mype
275              open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
276              open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
277              open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
278           else
279              if(mype.eq.pe_io) then
280                 write(fname_data,'(A,"_DES_DATA_",I4.4,".dat")') trim(run_name),tecplot_findex
281                 write(fname_extra,'(A,"_DES_EXTRA_",I4.4,".dat")') trim(run_name),tecplot_findex
282                 write(fname_eps,'(A,"_DES_EPS_",I4.4,".dat")') trim(run_name),tecplot_findex
283                 open(convert='big_endian',unit=des_data,file=fname_data,status='new',err=999)
284                 open(convert='big_endian',unit=des_ex,file=fname_extra,status='new',err=999)
285                 open(convert='big_endian',unit=des_eps,file=fname_eps,status='new',err=999)
286              end if
287           end if
288           tecplot_findex = tecplot_findex + 1
289     
290     ! write header
291           if (bdist_io .or. mype .eq. pe_io) then
292              if(DO_K) then
293                 write (des_data,'(A)') &
294                 'variables = "x" "y" "z" "vx" "vy" "vz" "rad" "den" "mark"'
295              else
296                 write (des_data, '(A)') &
297                 'variables = "x" "y" "vx" "vy" "omega" "rad" "den" "mark"'
298              endif
299              write (des_data, "(A,F15.7,A)")'zone T ="', s_time, '"'
300           end if
301     
302     ! write data in ditributed mode or as a single file
303           if(bdist_io) then
304              pc = 1
305              do l = 1,max_pip
306                 if(pc.gt.pip) exit
307                 if(is_nonexistent(l)) cycle
308                 pc = pc+1
309                 if(is_ghost(l) .or. is_entering_ghost(l) .or. is_exiting_ghost(l)) cycle
310                 if(DO_K) then
311                    write (des_data, '(8(2x,es12.5))')&
312                       (des_pos_new(k,l),k=1,wDIMN),(des_vel_new(k,l),k=1,wDIMN), &
313                        des_radius(l),ro_sol(l)
314                 else
315                    write (des_data, '(7(2x,es12.5))')&
316                       (des_pos_new(k,l),k=1,wDIMN), (des_vel_new(k,l),k=1,wDIMN), &
317                       omega_new(1,l), des_radius(l), ro_sol(l)
318                 endif
319             end do
320     
321           else ! if bdist_io
322     ! set parameters required for gathering info at PEIO and write as single file
323              lglocnt = 10
324              llocalcnt = pip - ighost_cnt
325              call global_sum(llocalcnt,lglocnt)
326              allocate (dprocbuf(llocalcnt),drootbuf(lglocnt),iprocbuf(llocalcnt),irootbuf(lglocnt))
327              allocate (ltemp_array(lglocnt,ltotvar))
328              igath_sendcnt = llocalcnt
329              lgathercnts = 0
330              lgathercnts(mype) = llocalcnt
331              call global_sum(lgathercnts,igathercnts)
332              idispls(0) = 0
333              do lproc = 1,numpes-1
334                 idispls(lproc) = idispls(lproc-1) + igathercnts(lproc-1)
335              end do
336     
337     ! gather information from all processor
338              lcount = 1
339              do k = 1,wDIMN
340                 call des_gather(des_pos_new(k,:))
341                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
342              end do
343              do k = 1,wDIMN
344                 call des_gather(des_vel_new(k,:))
345                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
346              end do
347              if(NO_K) then
348                 call des_gather(omega_new(1,:))
349                 ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
350              end if
351              call des_gather(des_radius)
352              ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
353              call des_gather(ro_sol)
354              ltemp_array(:,lcount) = drootbuf(:); lcount=lcount+1
355     
356     ! write the data into file
357              if (mype.eq.pe_io) then
358                 if(DO_K) then
359                    do l =1,lglocnt
360                       write (des_data,'(8(2x,es12.5),I5)') (ltemp_array(l,k),k=1,8),int(ltemp_array(l,9))
361                    end do
362                 else
363                    do l =1,lglocnt
364                       write (des_data,'(7(2x,es12.5),I5)') (ltemp_array(l,k),k=1,7),int(ltemp_array(l,8))
365                    end do
366                 end if
367              end if
368              deallocate (dprocbuf,drootbuf,iprocbuf,irootbuf,ltemp_array)
369           end if
370     ! close the files
371           if (mype.eq.pe_io .or. bdist_io) then
372              close(des_data)
373              close(des_ex)
374              close(des_eps)
375           end if
376           return
377     
378       999 write(*,"(/1x,70('*'),//,a,/,a,/1x,70('*'))")&
379              ' From: write_des_tecplot ',&
380              ' message: error opening des tecplot file. terminating run.'
381     
382     
383           END SUBROUTINE WRITE_DES_TECPLOT
384     !-----------------------------------------------
385     
386     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
387     !
388     !  Module name: WRITE_DES_BEDHEIGHT
389     !  Purpose: Writing DES output on bed height.
390     
391     !  WARNING: This code is out-of-date and should be modified for consistency
392     !  with current DEM version.  Also this routine will be fairly specific
393     !  to a user needs and should probably be tailored as such
394     
395     !
396     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
397           SUBROUTINE WRITE_DES_BEDHEIGHT
398     
399           USE param
400           USE param1
401           USE parallel
402           USE fldvar
403           USE discretelement
404           USE run
405           USE geometry
406           USE physprop
407           USE sendrecv
408           USE des_bc
409           IMPLICIT NONE
410     
411     !-----------------------------------------------
412     ! Local Variables
413     !-----------------------------------------------
414     ! logical that identifies that the data file has been created
415     ! and is already opened (initial checks/writes)
416           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
417     ! logical used for testing is the data files already exists
418           LOGICAL :: F_EXISTS
419     ! output file for the bed height data
420           CHARACTER(LEN=50)     :: FNAME_BH
421     ! file unit for the bed height data
422           INTEGER, PARAMETER :: BH_UNIT = 2010
423     ! dummy index values
424           INTEGER I, M
425     ! variables for bed height calculation
426           INTEGER, SAVE :: tcount = 1
427           DOUBLE PRECISION :: height_avg, height_rms
428           DOUBLE PRECISION, PARAMETER :: tmin = 5.d0
429           DOUBLE PRECISION, DIMENSION(5000), SAVE :: bed_height_time
430     
431     !-----------------------------------------------
432     ! Functions
433     !-----------------------------------------------
434     !-----------------------------------------------
435     
436     ! after tmin start storing bed height. after enough measurements
437     ! have been taken (i.e. tcount > 20) start to calculate a running
438     ! average bed height and running rms bed height for solids phase 1 only
439           height_avg = zero
440           height_rms = zero
441     
442           if(time.gt.tmin) then
443              if(tcount.le.5000)  then
444                 bed_height_time(tcount) = bed_height(1)
445                 !dt_time(tcount) = DT
446                 tcount = tcount + 1
447     
448                 if(tcount.gt.20)  then
449                    do i = 1, tcount-1,1
450                       height_avg = height_avg + bed_height_time(i)!*dt_time(i)
451                    enddo
452                    height_avg = height_avg/(tcount-1)
453                    do i = 1, tcount-1,1
454                       height_rms = height_rms + ((bed_height_time(i)&
455                            &-height_avg)**2)!*dt_time(i)
456                    enddo
457     
458                    height_rms = sqrt(height_rms/(tcount-1))
459                 endif
460              endif
461           endif
462     
463           FNAME_BH = TRIM(RUN_NAME)//'_DES_BEDHEIGHT.dat'
464           IF(FIRST_PASS) THEN
465              F_EXISTS = .FALSE.
466              INQUIRE(FILE=FNAME_BH,EXIST=F_EXISTS)
467     ! If the file does not exist, then create it with the necessary
468     ! header information.
469              IF (.NOT.F_EXISTS) THEN
470                 OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,&
471                    FORM="formatted",STATUS="new")
472              ELSE
473     ! To prevent overwriting existing files accidently, exit if the file
474     ! exists and this is a NEW run.
475                 IF(RUN_TYPE .EQ. 'NEW') THEN
476                    WRITE(*,1000)
477                    WRITE(UNIT_LOG, 1000)
478                    CALL MFIX_EXIT(myPE)
479                 ELSE
480     ! Open the file for appending of new data (RESTART_1 Case)
481                    OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,POSITION="append")
482                 ENDIF
483              ENDIF
484              FIRST_PASS = .FALSE.
485           ELSE
486     ! Open the file and mark for appending
487              OPEN(CONVERT='BIG_ENDIAN',UNIT=BH_UNIT,FILE=FNAME_BH,POSITION="append")
488           ENDIF
489     
490           WRITE(BH_UNIT, '(10(2X,E20.12))') s_time, &
491              (bed_height(M), M=1,DES_MMAX), height_avg, height_rms
492     ! Close the file and keep
493           CLOSE(BH_UNIT, STATUS="KEEP")
494     
495           RETURN
496     
497      1000 FORMAT(/1X,70('*')//, ' From: WRITE_DES_BEDHEIGHT',/,&
498              ' Message: bed_height.dat already exists in the run',&
499              ' directory.',/10X, 'Run terminated to prevent',&
500              ' accidental overwriting of files.',/1X,70('*')/)
501     
502           END SUBROUTINE WRITE_DES_BEDHEIGHT
503     !-----------------------------------------------
504     
505     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
506     !
507     !  Module name: WRITE_DES_THETA
508     !  Purpose: The following code writes out des_theta to a file for each
509     !  ijk cell in the system each time des_granular_temperature is called.
510     !
511     !
512     !
513     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
514           SUBROUTINE WRITE_DES_THETA
515     
516           USE param
517           USE param1
518           USE parallel
519           USE fldvar
520           USE discretelement
521           USE run
522           USE geometry
523           USE physprop
524           USE sendrecv
525           USE des_bc
526           USE functions
527           IMPLICIT NONE
528     
529     !-----------------------------------------------
530     ! Local Variables
531     !-----------------------------------------------
532     ! indices
533           INTEGER I, J, K, IJK
534     !
535           INTEGER M, NP
536     ! logical that identifies that the data file has been created
537     ! and is already opened (initial checks/writes)
538           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
539     ! logical used for testing is the data file already exists
540           LOGICAL :: F_EXISTS
541     ! file unit for the granular temperature data
542           INTEGER, PARAMETER :: GT_UNIT = 2020
543     ! output file for the granular temperature data
544           CHARACTER(LEN=50)  :: FNAME_GT
545     !-----------------------------------------------
546     
547           FNAME_GT = TRIM(RUN_NAME)//'_DES_THETA.dat'
548           IF (FIRST_PASS) THEN
549              F_EXISTS = .FALSE.
550              INQUIRE(FILE=FNAME_GT,EXIST=F_EXISTS)
551     
552              IF (.NOT.F_EXISTS) THEN
553     ! If the file does not exist, then create it with the necessary
554     ! header information.
555                 OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,STATUS='NEW')
556              ELSE
557                 IF(RUN_TYPE .EQ. 'NEW') THEN
558     ! If the run is new and the GT file already exists replace it with a
559     ! new file.
560     !               OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,STATUS='REPLACE')
561     ! Prevent overwriting an existing file by exiting if the file exists
562     ! and this is a NEW run.
563                    WRITE(*,1001) FNAME_GT
564                    WRITE(UNIT_LOG,1001) FNAME_GT
565                    CALL MFIX_EXIT(myPE)
566                 ELSE
567     ! Open the file for appending of new data (RESTART_1 Case)
568                    OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT, FILE=FNAME_GT, POSITION='APPEND')
569                 ENDIF
570              ENDIF
571              FIRST_PASS =  .FALSE.
572           ELSE
573     ! Open file and mark for appending
574              OPEN(CONVERT='BIG_ENDIAN',UNIT=GT_UNIT,FILE=FNAME_GT,POSITION='APPEND')
575           ENDIF   ! endif (first_pass)
576     
577           WRITE(GT_UNIT,*) ''
578           WRITE(GT_UNIT,'(A6,ES24.16)') 'Time=', S_TIME
579           WRITE(GT_UNIT,'(A6,2X,3(A6,2X),A8)',ADVANCE="NO") 'IJK', &
580              'I', 'J', 'K', 'NP'
581           DO M = 1,DES_MMAX
582              WRITE(GT_UNIT,'(7X,A6,I1)',ADVANCE="NO") 'THETA_',M
583           ENDDO
584           WRITE(GT_UNIT,*) ''
585           DO IJK = IJKSTART3, IJKEND3
586              IF(FLUID_AT(IJK)) THEN
587                 I = I_OF(IJK)
588                 J = J_OF(IJK)
589                 K = K_OF(IJK)
590                 NP = PINC(IJK)
591                 WRITE(GT_UNIT,'(I6,2X,3(I6,2X),I8,(2X,ES15.5))') &
592                    IJK, I, J, K, NP, (DES_THETA(IJK,M), M = 1,DES_MMAX)
593              ENDIF
594           ENDDO
595     
596     ! Close the file and keep
597           CLOSE(GT_UNIT, STATUS='KEEP')
598     
599           RETURN
600     
601      1001 FORMAT(/1X,70('*')//, ' From: WRITE_DES_THETA',/,&
602              ' Message: ', A, ' already exists in the run',/10X,&
603              'directory. Run terminated to prevent accidental overwriting',&
604              /10X,'of files.',/1X,70('*')/)
605     
606           END SUBROUTINE WRITE_DES_THETA
607     !-----------------------------------------------
608