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