MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
convert.t
Go to the documentation of this file.
1 !=============================================================================
2 subroutine generate_plotfile
8 !-----------------------------------------------------------------------------
9 
10 if(mype==0.and.level_io>0) write(unitterm,*)'reset tree to fixed level=',level_io
11 if(level_io>0 .or. level_io_min.ne.1 .or. level_io_max.ne.nlevelshi) then
13 else if(.not. phys_req_diagonal) then
14  call getbc(global_time,0.d0,ps,iwstart,nwgc)
15 end if
16 
17 select case(convert_type)
18  case('tecplot','tecplotCC','tecline')
19  call tecplot(unitconvert)
20  case('tecplotmpi','tecplotCCmpi','teclinempi')
22  case('vtu','vtuCC')
24  case('vtumpi','vtuCCmpi')
26  case('vtuB','vtuBCC','vtuBmpi','vtuBCCmpi')
28  case('pvtumpi','pvtuCCmpi')
30  case('pvtuBmpi','pvtuBCCmpi')
32  case('vtimpi','vtiCCmpi')
34  case('onegrid','onegridmpi')
35  call onegrid(unitconvert)
36  case('oneblock','oneblockB')
37  call oneblock(unitconvert)
38  case('user','usermpi')
39  if (.not. associated(usr_special_convert)) then
40  call mpistop("usr_special_convert not defined")
41  else
43  end if
44  case default
45  call mpistop("Error in generate_plotfile: Unknown convert_type")
46 end select
47 
48 
49 ! output synthetic euv emission
50 if (ndim==3) then
53 endif
54 
55 end subroutine generate_plotfile
56 !=============================================================================
57 subroutine oneblock(qunit)
58 ! this is for turning an AMR run into a single block
59 ! the data will be all on selected level level_io
60 
61 ! this version should work for any dimension
62 ! only writes w_write selected 1:nw variables, also nwauxio
63 ! may use saveprim to switch to primitives
64 ! this version can not work on multiple CPUs
65 ! does not renormalize variables
66 
67 ! header info differs from onegrid below
68 
69 ! ASCII or binary output
70 
74 use mod_physics
76 integer, intent(in) :: qunit
77 
78 integer :: Morton_no,igrid,ix^D,ig^D,level
79 integer, pointer :: ig_to_igrid(:^d&,:)
80 logical :: fileopen,writeblk(max_blocks)
81 character(len=80) :: filename
82 integer :: filenr,ncells^D,ncellx^D,jg^D,jig^D
83 
84 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
85 character(len=1024) :: outfilehead
86 
87 double precision :: wval1,xval1
88 double precision, dimension({^D&1:1},1:nw+nwauxio) :: wval
89 double precision, dimension({^D&1:1},1:ndim) :: xval
90 double precision:: normconv(0:nw+nwauxio)
91 
92 integer :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
93 logical :: patchw(ixg^t)
94 !-----------------------------------------------------------------------------
95 
96 if(level_io<1)then
97  call mpistop('please specify level_io>0 for usage with oneblock')
98 end if
99 
100 if(autoconvert)then
101  call mpistop('Set autoconvert=F and convert oneblock data manually')
102 end if
103 
104 if(npe>1)then
105  if(mype==0) print *,'ONEBLOCK as yet to be parallelized'
106  call mpistop('npe>1, oneblock')
107 end if
108 
109 ! only variables selected by w_write will be written out
110 normconv(0:nw+nwauxio)=one
111 normconv(0) = length_convert_factor
112 normconv(1:nw) = w_convert_factor
113 writenw=count(w_write(1:nw))+nwauxio
114 iiw=0
115 do iw =1,nw
116  if (.not.w_write(iw))cycle
117  iiw=iiw+1
118  iwrite(iiw)=iw
119 end do
120 if(nwauxio>0)then
121  do iw =nw+1,nw+nwauxio
122  iiw=iiw+1
123  iwrite(iiw)=iw
124  end do
125 endif
126 
127 allocate(ig_to_igrid(ng^d(level_io),0:npe-1))
128 ig_to_igrid=-1
129 writeblk=.false.
130 do morton_no=morton_start(mype),morton_stop(mype)
131  igrid=sfc_to_igrid(morton_no)
132  level=node(plevel_,igrid)
133  ig^d=igrid_to_node(igrid,mype)%node%ig^d;
134  ig_to_igrid(ig^d,mype)=igrid
135  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
136  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
137  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
138  writeblk(igrid)=.true.
139  end if
140 end do
141 
142 call getheadernames(wnamei,xandwnamei,outfilehead)
143 ncells^d=0;
144 ncellx^d=ixmhi^d-ixmlo^d+1\
145 {do ig^d=1,ng^d(level_io)\}
146  igrid=ig_to_igrid(ig^d,mype)
147  if(writeblk(igrid)) go to 20
148 {end do\}
149 20 continue
150 jg^d=ig^d;
151 {
152 jig^dd=jg^dd;
153 do ig^d=1,ng^d(level_io)
154  jig^d=ig^d
155  igrid=ig_to_igrid(jig^dd,mype)
156  if(writeblk(igrid)) ncells^d=ncells^d+ncellx^d
157 end do
158 \}
159 
160 do iigrid=1,igridstail; igrid=igrids(iigrid)
161  if(.not.writeblk(igrid)) cycle
162  block=>ps(igrid)
163  if (nwauxio > 0) then
164  if (.not. associated(usr_aux_output)) then
165  call mpistop("usr_aux_output not defined")
166  else
167  call usr_aux_output(ixg^ll,ixm^ll^ladd1, &
168  ps(igrid)%w,ps(igrid)%x,normconv)
169  end if
170  end if
171 end do
172 
173 if (saveprim) then
174  do iigrid=1,igridstail; igrid=igrids(iigrid)
175  if (.not.writeblk(igrid)) cycle
176  block=>ps(igrid)
177  call phys_to_primitive(ixg^ll,ixg^ll^lsub1,ps(igrid)%w,ps(igrid)%x)
178  if(b0field) then
179  ! add background magnetic field B0 to B
180  ps(igrid)%w(ixg^t,iw_mag(:))=ps(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
181  end if
182  end do
183 else
184  do iigrid=1,igridstail; igrid=igrids(iigrid)
185  if (.not.writeblk(igrid)) cycle
186  block=>ps(igrid)
187  if (b0field) then
188  ! add background magnetic field B0 to B
189  if(phys_energy) &
190  ps(igrid)%w(ixg^t,iw_e)=ps(igrid)%w(ixg^t,iw_e)+0.5d0*sum(ps(igrid)%B0(ixg^t,:,0)**2,dim=ndim+1) &
191  + sum(ps(igrid)%w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
192  ps(igrid)%w(ixg^t,iw_mag(:))=ps(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
193  end if
194  end do
195 end if
196 
197 master_cpu_open : if (mype == 0) then
198  inquire(qunit,opened=fileopen)
199  if (.not.fileopen) then
200  ! generate filename
201  filenr=snapshotini
202  if (autoconvert) filenr=snapshotnext
203  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
204  select case(convert_type)
205  case("oneblock")
206  open(qunit,file=filename,status='unknown')
207  write(qunit,*) trim(outfilehead)
208  write(qunit,*) ncells^d
209  write(qunit,*) real(global_time*time_convert_factor)
210  case("oneblockB")
211  open(qunit,file=filename,form='unformatted',status='unknown')
212  write(qunit) outfilehead
213  write(qunit) ncells^d
214  write(qunit) real(global_time*time_convert_factor)
215  end select
216  end if
217 end if master_cpu_open
218 
219 {^ifthreed
220 do ig3=1,ng3(level_io)
221  do ix3=ixmlo3,ixmhi3}
222 
223  {^nooned
224  do ig2=1,ng2(level_io)
225  do ix2=ixmlo2,ixmhi2}
226 
227  do ig1=1,ng1(level_io)
228  igrid=ig_to_igrid(ig^d,mype)
229  if(.not.writeblk(igrid)) cycle
230  do ix1=ixmlo1,ixmhi1
231  master_write : if(mype==0) then
232  select case(convert_type)
233  case("oneblock")
234  write(qunit,fmt="(100(e14.6))") &
235  ps(igrid)%x(ix^d,1:ndim)*normconv(0),&
236  (ps(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
237  case("oneblockB")
238  write(qunit) real(ps(igrid)%x(ix^D,1:ndim)*normconv(0)),&
239  (real(ps(igrid)%w(ix^D,iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
240  end select
241  end if master_write
242  end do
243  end do
244  {^nooned
245  end do
246  end do}
247  {^ifthreed
248  end do
249 end do}
250 
251 close(qunit)
252 
253 end subroutine oneblock
254 !=============================================================================
255 subroutine onegrid(qunit)
257 ! this is for turning an AMR run into a single grid
258 ! this version should work for any dimension, can be in parallel
259 ! in 1D, should behave much like oneblock, except for header info
260 
261 ! only writes all 1:nw variables, no nwauxio
262 ! may use saveprim to switch to primitives
263 ! this version can work on multiple CPUs
264 ! does not renormalize variables
265 ! ASCII output
266 
269 use mod_physics
271 integer, intent(in) :: qunit
272 
273 integer :: itag,Morton_no,igrid,ix^D,iw
274 logical :: fileopen
275 character(len=80) :: filename
276 integer :: filenr
277 
278 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
279 character(len=1024) :: outfilehead
280 
281 !.. MPI variables ..
282 integer :: igrid_recv,ipe
283 double precision :: w_recv(ixg^t,1:nw),x_recv(ixg^t,1:ndim)
284 integer, allocatable :: intstatus(:,:)
285 
286 !-----------------------------------------------------------------------------
287 
288 if(nwauxio>0)then
289  if(mype==0) print *,'ONEGRID to be used without nwauxio'
290  call mpistop('nwauxio>0, onegrid')
291 end if
292 
293 if(saveprim)then
294  if(mype==0.and.nwaux>0) print *,'warning: ONEGRID used with saveprim, check auxiliaries'
295 end if
296 
297 
298 
299 master_cpu_open : if (mype == 0) then
300  call getheadernames(wnamei,xandwnamei,outfilehead)
301  write(outfilehead,'(a)') "#"//" "//trim(outfilehead)
302  inquire(qunit,opened=fileopen)
303  if (.not.fileopen) then
304  ! generate filename
305  filenr=snapshotini
306  if (autoconvert) filenr=snapshotnext
307  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
308  open(qunit,file=filename,status='unknown')
309  end if
310  write(qunit,"(a)")outfilehead
311  write(qunit,"(i7)") ( {^d&(ixmhi^d-ixmlo^d+1)*} )*(morton_stop(npe-1)-morton_start(0)+1)
312 end if master_cpu_open
313 
314 do morton_no=morton_start(mype),morton_stop(mype)
315  igrid=sfc_to_igrid(morton_no)
316  if(saveprim) call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
317  if (mype/=0)then
318  itag=morton_no
319  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
320  call mpi_send(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
321  itag=igrid
322  call mpi_send(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
323  else
324  {do ix^db=ixmlo^db,ixmhi^db\}
325  do iw=1,nw
326  if( dabs(ps(igrid)%w(ix^d,iw)) < 1.0d-32 ) ps(igrid)%w(ix^d,iw) = zero
327  enddo
328  write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix^d,1:ndim)&
329  ,ps(igrid)%w(ix^d,1:nw)
330  {end do\}
331  end if
332 end do
333 
334 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
335 
336 manycpu : if (npe>1) then
337  if (mype==0) then
338  loop_cpu : do ipe =1, npe-1
339  loop_morton : do morton_no=morton_start(ipe),morton_stop(ipe)
340  itag=morton_no
341  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
342  call mpi_recv(x_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
343  itag=igrid_recv
344  call mpi_recv(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
345  {do ix^db=ixmlo^db,ixmhi^db\}
346  do iw=1,nw
347  if( dabs(ps(igrid)%w(ix^d,iw)) < smalldouble ) ps(igrid)%w(ix^d,iw) = zero
348  enddo
349  write(qunit,fmt="(100(e14.6))") x_recv(ix^d,1:ndim)&
350  ,w_recv(ix^d,1:nw)
351  {end do\}
352  end do loop_morton
353  end do loop_cpu
354  end if
355 end if manycpu
356 
357 if (npe>1) then
358  call mpi_barrier(icomm,ierrmpi)
359  if(mype==0)deallocate(intstatus)
360 endif
361 
362 if(mype==0) close(qunit)
363 end subroutine onegrid
364 !============================================================================
365 subroutine tecplot(qunit)
367 ! output for tecplot (ASCII format)
368 ! not parallel, uses calc_grid to compute nwauxio variables
369 ! allows renormalizing using convert factors
370 
373 
374 integer, intent(in) :: qunit
375 
376 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
377 integer:: NumGridsOnLevel(1:nlevelshi)
378 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
379 
380 integer :: nodes, elems
381 
382 
383 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
384 
385 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
386 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
387 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
388 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
389 
390 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
391 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
392 double precision, dimension(0:nw+nwauxio) :: normconv
393 logical :: fileopen,first
394 character(len=80) :: filename
395 integer :: filenr
396 
397 !!! possible length conflict
398 character(len=1024) :: tecplothead
399 
400 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
401 character(len=1024) :: outfilehead
402 !-----------------------------------------------------------------------------
403 
404 if(npe>1)then
405  if(mype==0) print *,'tecplot not parallel, use tecplotmpi'
406  call mpistop('npe>1, tecplot')
407 end if
408 
409 if(nw/=count(w_write(1:nw)))then
410  if(mype==0) print *,'tecplot does not use w_write=F'
411  call mpistop('w_write, tecplot')
412 end if
413 
414 if(nocartesian)then
415  if(mype==0) print *,'tecplot with nocartesian'
416 endif
417 
418 inquire(qunit,opened=fileopen)
419 if (.not.fileopen) then
420  ! generate filename
421  filenr=snapshotini
422  if (autoconvert) filenr=snapshotnext
423  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
424  open(qunit,file=filename,status='unknown')
425 end if
426 
427 call getheadernames(wnamei,xandwnamei,outfilehead)
428 
429 write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
430 
431 write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
432 
433 numgridsonlevel(1:nlevelshi)=0
434 do level=levmin,levmax
435  numgridsonlevel(level)=0
436  do iigrid=1,igridstail; igrid=igrids(iigrid);
437  if (node(plevel_,igrid)/=level) cycle
438  numgridsonlevel(level)=numgridsonlevel(level)+1
439  end do
440 end do
441 
442 nx^d=ixmhi^d-ixmlo^d+1;
443 nxc^d=nx^d+1;
444 
445 {^ifoned
446 if(convert_type=='tecline') then
447  nodes=0
448  elems=0
449  do level=levmin,levmax
450  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
451  elems=elems + numgridsonlevel(level)*{nx^d*}
452  enddo
453 
454  write(qunit,"(a,i7,a,1pe12.5,a)") &
455  'ZONE T="all levels", I=',elems, &
456  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
457 
458  igonlevel=0
459  do iigrid=1,igridstail; igrid=igrids(iigrid);
460  block=>ps(igrid)
461  call calc_x(igrid,xc,xcc)
462  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
463  {do ix^db=ixccmin^db,ixccmax^db\}
464  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
465  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
466  !write(qunit,fmt="(100(e14.6))") x_TEC, w_TEC
467  write(qunit,fmt="(100(e24.16))") x_tec, w_tec
468  {end do\}
469  enddo
470  close(qunit)
471 else
472 }
473 do level=levmin,levmax
474  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
475  elemsonlevel=numgridsonlevel(level)*{nx^d*}
476  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
477  ! with the AMR grid LEVEL. Other options would be
478  ! let each grid define a zone: inefficient for TECPLOT internal workings
479  ! hence not implemented
480  ! let entire octree define 1 zone: no difference in interpolation
481  ! properties across TECPLOT zones detected as yet, hence not done
482  select case(convert_type)
483  case('tecplot')
484  ! in this option, we store the corner coordinates, as well as the corner
485  ! values of all variables (obtained by averaging). This allows POINT packaging,
486  ! and thus we can save full grid info by using one call to calc_grid
487  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
488  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
489  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
490  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
491  do iigrid=1,igridstail; igrid=igrids(iigrid);
492  if (node(plevel_,igrid)/=level) cycle
493  block=>ps(igrid)
494  call calc_x(igrid,xc,xcc)
495  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
496  ixc^l,ixcc^l,.true.)
497  {do ix^db=ixcmin^db,ixcmax^db\}
498  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
499  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
500  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
501  {end do\}
502  enddo
503  case('tecplotCC')
504  ! in this option, we store the corner coordinates, and the cell center
505  ! values of all variables. Due to this mix of corner/cell center, we must
506  ! use BLOCK packaging, and thus we have enormous overhead by using
507  ! calc_grid repeatedly to merely fill values of cell corner coordinates
508  ! and cell center values per dimension, per variable
509  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
510  if(nw+nwauxio==1)then
511  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
512  ! and just set [ndim+1]
513  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
514  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
515  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
516  ndim+1,']=CELLCENTERED), ZONETYPE=', &
517  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
518  else
519  if(ndim+nw+nwauxio<10) then
520  ! difference only in length of integer format specification for ndim+nw+nwauxio
521  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
522  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
523  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
524  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
525  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
526  else
527  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
528  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
529  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
530  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
531  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
532  endif
533  endif
534  do idim=1,ndim
535  first=(idim==1)
536  do iigrid=1,igridstail; igrid=igrids(iigrid);
537  if (node(plevel_,igrid)/=level) cycle
538  block=>ps(igrid)
539  call calc_x(igrid,xc,xcc)
540  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
541  ixc^l,ixcc^l,first)
542  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
543  enddo
544  enddo
545  do iw=1,nw+nwauxio
546  do iigrid=1,igridstail; igrid=igrids(iigrid);
547  if (node(plevel_,igrid)/=level) cycle
548  block=>ps(igrid)
549  call calc_x(igrid,xc,xcc)
550  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
551  ixc^l,ixcc^l,.true.)
552  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
553  enddo
554  enddo
555  case default
556  call mpistop('no such tecplot type')
557  end select
558  igonlevel=0
559  do iigrid=1,igridstail; igrid=igrids(iigrid);
560  if (node(plevel_,igrid)/=level) cycle
561  block=>ps(igrid)
562  igonlevel=igonlevel+1
563  call save_conntec(qunit,igrid,igonlevel)
564  end do
565 end do
566 {^ifoned endif}
567 
568 close(qunit)
569 
570 end subroutine tecplot
571 !=============================================================================
572 subroutine save_conntec(qunit,igrid,igonlevel)
574 ! this saves the basic line, quad and brick connectivity,
575 ! as used by TECPLOT file outputs for unstructured grid
576 
578 
579 integer, intent(in) :: qunit, igrid, igonlevel
580 
581 integer :: nx^D, nxC^D, ix^D
582 {^ifoned integer, external:: nodenumbertec1D \}
583 {^iftwod integer, external:: nodenumbertec2D \}
584 {^ifthreed integer, external:: nodenumbertec3D \}
585 !-----------------------------------------------------------------------------
586 nx^d=ixmhi^d-ixmlo^d+1;
587 nxc^d=nx^d+1;
588 
589 ! connectivity list
590 {do ix^db=1,nx^db\}
591  {^ifthreed
592  ! basic brick connectivity
593  write(qunit,'(8(i7,1x))') &
594  nodenumbertec3d(ix1, ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
595  nodenumbertec3d(ix1+1,ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
596  nodenumbertec3d(ix1+1,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
597  nodenumbertec3d(ix1 ,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
598  nodenumbertec3d(ix1 ,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
599  nodenumbertec3d(ix1+1,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
600  nodenumbertec3d(ix1+1,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
601  nodenumbertec3d(ix1 ,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid)
602  }
603  {^iftwod
604  ! basic quadrilateral connectivity
605  write(qunit,'(4(i7,1x))') &
606  nodenumbertec2d(ix1, ix2-1,nxc1,nxc2,igonlevel,igrid),&
607  nodenumbertec2d(ix1+1,ix2-1,nxc1,nxc2,igonlevel,igrid),&
608  nodenumbertec2d(ix1+1,ix2 ,nxc1,nxc2,igonlevel,igrid),&
609  nodenumbertec2d(ix1 ,ix2 ,nxc1,nxc2,igonlevel,igrid)
610  }
611  {^ifoned
612  ! basic line connectivity
613  write(qunit,'(2(i7,1x))') nodenumbertec1d(ix1,nxc1,igonlevel,igrid),&
614  nodenumbertec1d(ix1+1,nxc1,igonlevel,igrid)
615  }
616 {end do\}
617 
618 end subroutine save_conntec
619 !=============================================================================
620 integer function nodenumbertec1d(i1,nx1,ig,igrid)
622 integer, intent(in):: i1,nx1,ig,igrid
623 !-----------------------------------------------------------------------------
624 nodenumbertec1d=i1+(ig-1)*nx1
625 
626 if(nodenumbertec1d>9999999)call mpistop("too large nodenumber")
627 end function nodenumbertec1d
628 !====================================================================================
629 integer function nodenumbertec2d(i1,i2,nx1,nx2,ig,igrid)
631 integer, intent(in):: i1,i2,nx1,nx2,ig,igrid
632 !-----------------------------------------------------------------------------
633 nodenumbertec2d=i1+i2*nx1+(ig-1)*nx1*nx2
634 
635 if(nodenumbertec2d>9999999)call mpistop("too large nodenumber")
636 end function nodenumbertec2d
637 !====================================================================================
638 integer function nodenumbertec3d(i1,i2,i3,nx1,nx2,nx3,ig,igrid)
640 integer, intent(in):: i1,i2,i3,nx1,nx2,nx3,ig,igrid
641 !-----------------------------------------------------------------------------
642 nodenumbertec3d=i1+i2*nx1+i3*nx1*nx2+(ig-1)*nx1*nx2*nx3
643 
644 if(nodenumbertec3d>9999999)call mpistop("too large nodenumber")
645 end function nodenumbertec3d
646 !=============================================================================
647 subroutine unstructuredvtk(qunit)
649 ! output for vtu format to paraview
650 ! not parallel, uses calc_grid to compute nwauxio variables
651 ! allows renormalizing using convert factors
652 
655 
656 integer, intent(in) :: qunit
657 
658 double precision :: x_VTK(1:3)
659 
660 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
661 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
662 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
663 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
664 
665 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
666 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
667 double precision, dimension(0:nw+nwauxio) :: normconv
668 integer:: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,iw
669 integer:: NumGridsOnLevel(1:nlevelshi)
670 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
671 
672 character(len=80):: filename
673 integer :: filenr
674 
675 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
676 character(len=1024) :: outfilehead
677 
678 logical :: fileopen
679 !-----------------------------------------------------------------------------
680 
681 if(npe>1)then
682  if(mype==0) print *,'unstructuredvtk not parallel, use vtumpi'
683  call mpistop('npe>1, unstructuredvtk')
684 end if
685 
686 inquire(qunit,opened=fileopen)
687 if(.not.fileopen)then
688  ! generate filename
689  filenr=snapshotini
690  if (autoconvert) filenr=snapshotnext
691  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
692  ! Open the file for the header part
693  open(qunit,file=filename,status='unknown')
694 end if
695 
696 call getheadernames(wnamei,xandwnamei,outfilehead)
697 
698 ! generate xml header
699 write(qunit,'(a)')'<?xml version="1.0"?>'
700 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
701 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
702 write(qunit,'(a)')'<UnstructuredGrid>'
703 write(qunit,'(a)')'<FieldData>'
704 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
705  'NumberOfTuples="1" format="ascii">'
706 write(qunit,*) dble(dble(global_time)*time_convert_factor)
707 write(qunit,'(a)')'</DataArray>'
708 write(qunit,'(a)')'</FieldData>'
709 
710 ! number of cells, number of corner points, per grid.
711 nx^d=ixmhi^d-ixmlo^d+1;
712 nxc^d=nx^d+1;
713 nc={nx^d*}
714 np={nxc^d*}
715 
716 ! Note: using the w_write, writelevel, writespshift
717 ! we can clip parts of the grid away, select variables, levels etc.
718 do level=levmin,levmax
719  if (writelevel(level)) then
720  do iigrid=1,igridstail; igrid=igrids(iigrid);
721  if (node(plevel_,igrid)/=level) cycle
722  block=>ps(igrid)
723  ! only output a grid when fully within clipped region selected
724  ! by writespshift array
725  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
726  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
727  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
728  call calc_x(igrid,xc,xcc)
729  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
730  ixc^l,ixcc^l,.true.)
731  select case(convert_type)
732  case('vtu')
733  ! we write out every grid as one VTK PIECE
734  write(qunit,'(a,i7,a,i7,a)') &
735  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
736  write(qunit,'(a)')'<PointData>'
737  do iw=1,nw+nwauxio
738  if(iw<=nw) then
739  if(.not.w_write(iw)) cycle
740  end if
741  write(qunit,'(a,a,a)')&
742  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
743  write(qunit,'(200(1pe14.6))') {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
744  write(qunit,'(a)')'</DataArray>'
745  end do
746  write(qunit,'(a)')'</PointData>'
747  write(qunit,'(a)')'<Points>'
748  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
749  ! write cell corner coordinates in a backward dimensional loop, always 3D output
750  {do ix^db=ixcmin^db,ixcmax^db \}
751  x_vtk(1:3)=zero;
752  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
753  write(qunit,'(3(1pe14.6))') x_vtk
754  {end do \}
755  write(qunit,'(a)')'</DataArray>'
756  write(qunit,'(a)')'</Points>'
757  case('vtuCC')
758  ! we write out every grid as one VTK PIECE
759  write(qunit,'(a,i7,a,i7,a)') &
760  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
761  write(qunit,'(a)')'<CellData>'
762  do iw=1,nw+nwauxio
763  if(iw<=nw) then
764  if(.not.w_write(iw)) cycle
765  end if
766  write(qunit,'(a,a,a)')&
767  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
768  write(qunit,'(200(1pe14.6))') {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
769  write(qunit,'(a)')'</DataArray>'
770  end do
771  write(qunit,'(a)')'</CellData>'
772  write(qunit,'(a)')'<Points>'
773  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
774  ! write cell corner coordinates in a backward dimensional loop, always 3D output
775  {do ix^db=ixcmin^db,ixcmax^db \}
776  x_vtk(1:3)=zero;
777  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
778  write(qunit,'(3(1pe14.6))') x_vtk
779  {end do \}
780  write(qunit,'(a)')'</DataArray>'
781  write(qunit,'(a)')'</Points>'
782  end select
783 
784  write(qunit,'(a)')'<Cells>'
785  ! connectivity part
786  write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
787  call save_connvtk(qunit,igrid)
788  write(qunit,'(a)')'</DataArray>'
789 
790  ! offsets data array
791  write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
792  do icel=1,nc
793  write(qunit,'(i7)') icel*(2**^nd)
794  end do
795  write(qunit,'(a)')'</DataArray>'
796 
797  ! VTK cell type data array
798  write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
799  ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
800  {^ifoned vtk_type=3 \}
801  {^iftwod vtk_type=8 \}
802  {^ifthreed vtk_type=11 \}
803  do icel=1,nc
804  write(qunit,'(i2)') vtk_type
805  end do
806  write(qunit,'(a)')'</DataArray>'
807 
808  write(qunit,'(a)')'</Cells>'
809 
810  write(qunit,'(a)')'</Piece>'
811  end if
812  end do
813  end if
814 end do
815 
816 write(qunit,'(a)')'</UnstructuredGrid>'
817 write(qunit,'(a)')'</VTKFile>'
818 close(qunit)
819 
820 end subroutine unstructuredvtk
821 !====================================================================================
822 subroutine unstructuredvtkb(qunit)
824 ! output for vtu format to paraview, binary version output
825 ! not parallel, uses calc_grid to compute nwauxio variables
826 ! allows renormalizing using convert factors
830 
831 integer, intent(in) :: qunit
832 
833 double precision :: x_VTK(1:3)
834 
835 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
836 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
837 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
838 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
839 
840 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
841 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
842 double precision :: normconv(0:nw+nwauxio)
843 
844 integer, allocatable :: intstatus(:,:)
845 integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
846 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
847 integer*8 :: offset
848 
849 integer:: k,iw
850 integer:: length,lengthcc,length_coords,length_conn,length_offsets
851 character:: buf
852 character(len=80):: filename
853 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
854 character(len=1024) :: outfilehead
855 
856 logical :: fileopen,cell_corner=.false.
857 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
858 !-----------------------------------------------------------------------------
859 
860 normconv=one
861 morton_length=morton_stop(npe-1)-morton_start(0)+1
862 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
863 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
864 morton_aim=.false.
865 morton_aim_p=.false.
866 do morton_no=morton_start(mype),morton_stop(mype)
867  igrid=sfc_to_igrid(morton_no)
868  level=node(plevel_,igrid)
869  ! we can clip parts of the grid away, select variables, levels etc.
870  if(writelevel(level)) then
871  ! only output a grid when fully within clipped region selected
872  ! by writespshift array
873  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
874  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
875  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
876  morton_aim_p(morton_no)=.true.
877  end if
878  end if
879 end do
880 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
881  icomm,ierrmpi)
882 select case(convert_type)
883  case('vtuB','vtuBmpi')
884  cell_corner=.true.
885  case('vtuBCC','vtuBCCmpi')
886  cell_corner=.false.
887 end select
888 if (mype /= 0) then
889  do morton_no=morton_start(mype),morton_stop(mype)
890  if(.not. morton_aim(morton_no)) cycle
891  igrid=sfc_to_igrid(morton_no)
892  call calc_x(igrid,xc,xcc)
893  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
894  ixc^l,ixcc^l,.true.)
895  itag=morton_no
896  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
897  if(cell_corner) then
898  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
899  else
900  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
901  end if
902  end do
903 
904 else
905  ! mype==0
906  offset=0
907  inquire(qunit,opened=fileopen)
908  if(.not.fileopen)then
909  ! generate filename
910  filenr=snapshotini
911  if (autoconvert) filenr=snapshotnext
912  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
913  ! Open the file for the header part
914  open(qunit,file=filename,status='replace')
915  end if
916  call getheadernames(wnamei,xandwnamei,outfilehead)
917  ! generate xml header
918  write(qunit,'(a)')'<?xml version="1.0"?>'
919  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
920  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
921  write(qunit,'(a)')'<UnstructuredGrid>'
922  write(qunit,'(a)')'<FieldData>'
923  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
924  'NumberOfTuples="1" format="ascii">'
925  write(qunit,*) real(global_time*time_convert_factor)
926  write(qunit,'(a)')'</DataArray>'
927  write(qunit,'(a)')'</FieldData>'
928 
929  ! number of cells, number of corner points, per grid.
930  nx^d=ixmhi^d-ixmlo^d+1;
931  nxc^d=nx^d+1;
932  nc={nx^d*}
933  np={nxc^d*}
934  length=np*size_real
935  lengthcc=nc*size_real
936  length_coords=3*length
937  length_conn=2**^nd*size_int*nc
938  length_offsets=nc*size_int
939 
940  ! Note: using the w_write, writelevel, writespshift
941  do morton_no=morton_start(0),morton_stop(0)
942  if(.not. morton_aim(morton_no)) cycle
943  if(cell_corner) then
944  ! we write out every grid as one VTK PIECE
945  write(qunit,'(a,i7,a,i7,a)') &
946  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
947  write(qunit,'(a)')'<PointData>'
948  do iw=1,nw+nwauxio
949  if(iw<=nw) then
950  if(.not.w_write(iw)) cycle
951  endif
952 
953  write(qunit,'(a,a,a,i16,a)')&
954  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
955  '" format="appended" offset="',offset,'">'
956  write(qunit,'(a)')'</DataArray>'
957  offset=offset+length+size_int
958  end do
959  write(qunit,'(a)')'</PointData>'
960  write(qunit,'(a)')'<Points>'
961  write(qunit,'(a,i16,a)') &
962  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
963  ! write cell corner coordinates in a backward dimensional loop, always 3D output
964  offset=offset+length_coords+size_int
965  write(qunit,'(a)')'</Points>'
966  else
967  ! we write out every grid as one VTK PIECE
968  write(qunit,'(a,i7,a,i7,a)') &
969  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
970  write(qunit,'(a)')'<CellData>'
971  do iw=1,nw+nwauxio
972  if(iw<=nw) then
973  if(.not.w_write(iw)) cycle
974  end if
975  write(qunit,'(a,a,a,i16,a)')&
976  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
977  '" format="appended" offset="',offset,'">'
978  write(qunit,'(a)')'</DataArray>'
979  offset=offset+lengthcc+size_int
980  end do
981  write(qunit,'(a)')'</CellData>'
982  write(qunit,'(a)')'<Points>'
983  write(qunit,'(a,i16,a)') &
984  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
985  ! write cell corner coordinates in a backward dimensional loop, always 3D output
986  offset=offset+length_coords+size_int
987  write(qunit,'(a)')'</Points>'
988  end if
989  write(qunit,'(a)')'<Cells>'
990 
991  ! connectivity part
992  write(qunit,'(a,i16,a)')&
993  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
994  offset=offset+length_conn+size_int
995 
996  ! offsets data array
997  write(qunit,'(a,i16,a)') &
998  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
999  offset=offset+length_offsets+size_int
1000 
1001  ! VTK cell type data array
1002  write(qunit,'(a,i16,a)') &
1003  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1004  offset=offset+size_int+nc*size_int
1005 
1006  write(qunit,'(a)')'</Cells>'
1007 
1008  write(qunit,'(a)')'</Piece>'
1009  end do
1010  ! write metadata communicated from other processors
1011  if(npe>1)then
1012  do ipe=1, npe-1
1013  do morton_no=morton_start(ipe),morton_stop(ipe)
1014  if(.not. morton_aim(morton_no)) cycle
1015  if(cell_corner) then
1016  ! we write out every grid as one VTK PIECE
1017  write(qunit,'(a,i7,a,i7,a)') &
1018  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1019  write(qunit,'(a)')'<PointData>'
1020  do iw=1,nw+nwauxio
1021  if(iw<=nw) then
1022  if(.not.w_write(iw)) cycle
1023  end if
1024  write(qunit,'(a,a,a,i16,a)')&
1025  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1026  '" format="appended" offset="',offset,'">'
1027  write(qunit,'(a)')'</DataArray>'
1028  offset=offset+length+size_int
1029  end do
1030  write(qunit,'(a)')'</PointData>'
1031  write(qunit,'(a)')'<Points>'
1032  write(qunit,'(a,i16,a)') &
1033  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1034  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1035  offset=offset+length_coords+size_int
1036  write(qunit,'(a)')'</Points>'
1037  else
1038  ! we write out every grid as one VTK PIECE
1039  write(qunit,'(a,i7,a,i7,a)') &
1040  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1041  write(qunit,'(a)')'<CellData>'
1042  do iw=1,nw+nwauxio
1043  if(iw<=nw) then
1044  if(.not.w_write(iw)) cycle
1045  end if
1046  write(qunit,'(a,a,a,i16,a)')&
1047  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1048  '" format="appended" offset="',offset,'">'
1049  write(qunit,'(a)')'</DataArray>'
1050  offset=offset+lengthcc+size_int
1051  end do
1052  write(qunit,'(a)')'</CellData>'
1053  write(qunit,'(a)')'<Points>'
1054  write(qunit,'(a,i16,a)') &
1055  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1056  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1057  offset=offset+length_coords+size_int
1058  write(qunit,'(a)')'</Points>'
1059  end if
1060  write(qunit,'(a)')'<Cells>'
1061  ! connectivity part
1062  write(qunit,'(a,i16,a)')&
1063  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1064  offset=offset+length_conn+size_int
1065 
1066  ! offsets data array
1067  write(qunit,'(a,i16,a)') &
1068  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1069  offset=offset+length_offsets+size_int
1070 
1071  ! VTK cell type data array
1072  write(qunit,'(a,i16,a)') &
1073  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1074  offset=offset+size_int+nc*size_int
1075 
1076  write(qunit,'(a)')'</Cells>'
1077 
1078  write(qunit,'(a)')'</Piece>'
1079  end do
1080  end do
1081  end if
1082 
1083  write(qunit,'(a)')'</UnstructuredGrid>'
1084  write(qunit,'(a)')'<AppendedData encoding="raw">'
1085  close(qunit)
1086  open(qunit,file=filename,access='stream',form='unformatted',position='append')
1087  buf='_'
1088  write(qunit) trim(buf)
1089 
1090  do morton_no=morton_start(0),morton_stop(0)
1091  if(.not. morton_aim(morton_no)) cycle
1092  igrid=sfc_to_igrid(morton_no)
1093  call calc_x(igrid,xc,xcc)
1094  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1095  ixc^l,ixcc^l,.true.)
1096  do iw=1,nw+nwauxio
1097  if(iw<=nw) then
1098  if(.not.w_write(iw)) cycle
1099  end if
1100  if(cell_corner) then
1101  write(qunit) length
1102  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
1103  else
1104  write(qunit) lengthcc
1105  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
1106  end if
1107  end do
1108 
1109  write(qunit) length_coords
1110  {do ix^db=ixcmin^db,ixcmax^db \}
1111  x_vtk(1:3)=zero;
1112  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1113  do k=1,3
1114  write(qunit) real(x_vtk(k))
1115  end do
1116  {end do \}
1117 
1118  write(qunit) length_conn
1119  {do ix^db=1,nx^db\}
1120  {^ifoned write(qunit)ix1-1,ix1 \}
1121  {^iftwod
1122  write(qunit)(ix2-1)*nxc1+ix1-1, &
1123  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1124  \}
1125  {^ifthreed
1126  write(qunit)&
1127  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1128  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1129  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1130  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1131  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1132  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1133  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1134  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1135  \}
1136  {end do\}
1137 
1138  write(qunit) length_offsets
1139  do icel=1,nc
1140  write(qunit) icel*(2**^nd)
1141  end do
1142 
1143  {^ifoned vtk_type=3 \}
1144  {^iftwod vtk_type=8 \}
1145  {^ifthreed vtk_type=11 \}
1146  write(qunit) size_int*nc
1147  do icel=1,nc
1148  write(qunit) vtk_type
1149  end do
1150  end do
1151  allocate(intstatus(mpi_status_size,1))
1152  if(npe>1)then
1153  ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1154  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1155  do ipe=1, npe-1
1156  do morton_no=morton_start(ipe),morton_stop(ipe)
1157  if(.not. morton_aim(morton_no)) cycle
1158  itag=morton_no
1159  call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1160  if(cell_corner) then
1161  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1162  else
1163  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1164  end if
1165  do iw=1,nw+nwauxio
1166  if(iw<=nw) then
1167  if(.not.w_write(iw)) cycle
1168  end if
1169  if(cell_corner) then
1170  write(qunit) length
1171  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
1172  else
1173  write(qunit) lengthcc
1174  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
1175  end if
1176  end do
1177 
1178  write(qunit) length_coords
1179  {do ix^db=ixcmin^db,ixcmax^db \}
1180  x_vtk(1:3)=zero;
1181  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1182  do k=1,3
1183  write(qunit) real(x_vtk(k))
1184  end do
1185  {end do \}
1186 
1187  write(qunit) length_conn
1188  {do ix^db=1,nx^db\}
1189  {^ifoned write(qunit)ix1-1,ix1 \}
1190  {^iftwod
1191  write(qunit)(ix2-1)*nxc1+ix1-1, &
1192  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1193  \}
1194  {^ifthreed
1195  write(qunit)&
1196  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1197  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1198  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1199  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1200  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1201  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1202  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1203  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1204  \}
1205  {end do\}
1206 
1207  write(qunit) length_offsets
1208  do icel=1,nc
1209  write(qunit) icel*(2**^nd)
1210  end do
1211  {^ifoned vtk_type=3 \}
1212  {^iftwod vtk_type=8 \}
1213  {^ifthreed vtk_type=11 \}
1214  write(qunit) size_int*nc
1215  do icel=1,nc
1216  write(qunit) vtk_type
1217  end do
1218  end do
1219  end do
1220  end if
1221  close(qunit)
1222  open(qunit,file=filename,status='unknown',form='formatted',position='append')
1223  write(qunit,'(a)')'</AppendedData>'
1224  write(qunit,'(a)')'</VTKFile>'
1225  close(qunit)
1226  deallocate(intstatus)
1227 end if
1228 
1229 deallocate(morton_aim,morton_aim_p)
1230 if (npe>1) then
1231  call mpi_barrier(icomm,ierrmpi)
1232 end if
1233 
1234 end subroutine unstructuredvtkb
1235 !====================================================================================
1236 subroutine save_connvtk(qunit,igrid)
1238 ! this saves the basic line, pixel and voxel connectivity,
1239 ! as used by VTK file outputs for unstructured grid
1240 
1242 
1243 integer, intent(in) :: qunit, igrid
1244 
1245 integer :: nx^D, nxC^D, ix^D
1246 !-----------------------------------------------------------------------------
1247 nx^d=ixmhi^d-ixmlo^d+1;
1248 nxc^d=nx^d+1;
1249 
1250 {do ix^db=1,nx^db\}
1251  {^ifoned write(qunit,'(2(i7,1x))')ix1-1,ix1 \}
1252  {^iftwod
1253  write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
1254  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1255  \}
1256  {^ifthreed
1257  write(qunit,'(8(i7,1x))')&
1258  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1259  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1260  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1261  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1262  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1263  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1264  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1265  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1266  \}
1267 {end do\}
1268 
1269 end subroutine save_connvtk
1270 !=============================================================================
1271 subroutine imagedatavtk_mpi(qunit)
1273 ! output for vti format to paraview, non-binary version output
1274 ! parallel, uses calc_grid to compute nwauxio variables
1275 ! allows renormalizing using convert factors
1276 ! allows skipping of w_write selected variables
1277 
1278 ! implementation such that length of ASCII output is identical when
1279 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1280 
1283 use mod_calculate_xw
1284 
1285 integer, intent(in) :: qunit
1286 
1287 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1288 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1289 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1290 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1291 
1292 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1293 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1294 double precision, dimension(0:nw+nwauxio) :: normconv
1295 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1296 integer:: NumGridsOnLevel(1:nlevelshi)
1297 integer :: nx^D
1298 
1299 character(len=80):: filename
1300 integer :: filenr
1301 
1302 integer, allocatable :: intstatus(:,:)
1303 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1304 
1305 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1306 character(len=1024) :: outfilehead
1307 
1308 logical :: fileopen
1309 integer :: itag,ipe,Morton_no,Morton_length
1310 integer :: ixrvC^L, ixrvCC^L, siz_ind, ind_send(5*^nd), ind_recv(5*^nd)
1311 double precision :: origin(1:3), spacing(1:3)
1312 integer :: wholeExtent(1:6), ig^D
1313 type(tree_node_ptr) :: tree
1314 !-----------------------------------------------------------------------------
1315 if(levmin/=levmax) call mpistop('ImageData can only be used when levmin=levmax')
1316 
1317 normconv(0) = length_convert_factor
1318 normconv(1:nw) = w_convert_factor
1319 siz_ind=5*^nd
1320 morton_length=morton_stop(npe-1)-morton_start(0)+1
1321 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1322 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1323 morton_aim=.false.
1324 morton_aim_p=.false.
1325 do morton_no=morton_start(mype),morton_stop(mype)
1326  igrid=sfc_to_igrid(morton_no)
1327  level=node(plevel_,igrid)
1328  ! we can clip parts of the grid away, select variables, levels etc.
1329  if(writelevel(level)) then
1330  ! only output a grid when fully within clipped region selected
1331  ! by writespshift array
1332  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1333  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1334  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1335  morton_aim_p(morton_no)=.true.
1336  end if
1337  end if
1338 end do
1339 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1340  icomm,ierrmpi)
1341 
1342 
1343 if (mype /= 0) then
1344  do morton_no=morton_start(mype),morton_stop(mype)
1345  if(.not. morton_aim(morton_no)) cycle
1346  igrid=sfc_to_igrid(morton_no)
1347  call calc_x(igrid,xc,xcc)
1348  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1349  ixc^l,ixcc^l,.true.)
1350  tree%node => igrid_to_node(igrid, mype)%node
1351  {^d& ig^d = tree%node%ig^d; }
1352  itag=morton_no
1353  ind_send=(/ ixc^l,ixcc^l, ig^d /)
1354  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1355  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1356  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1357  end do
1358 
1359 else
1360 
1361  inquire(qunit,opened=fileopen)
1362  if(.not.fileopen)then
1363  ! generate filename
1364  filenr=snapshotini
1365  if (autoconvert) filenr=snapshotnext
1366  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vti"
1367  ! Open the file for the header part
1368  open(qunit,file=filename,status='unknown',form='formatted')
1369  endif
1370 
1371 call getheadernames(wnamei,xandwnamei,outfilehead)
1372 
1373 ! number of cells per grid.
1374 nx^d=ixmhi^d-ixmlo^d+1;
1375 
1376 origin = 0
1377 {^d& origin(^d) = xprobmin^d*normconv(0); }
1378 spacing = zero
1379 {^d&spacing(^d) = dxlevel(^d)*normconv(0); }
1380 
1381 wholeextent = 0
1382 ! if we use writespshift, the whole extent has to be calculated:
1383 {^d&wholeextent(^d*2-1) = nx^d * ceiling(((xprobmax^d-xprobmin^d)*writespshift(^d,1)) &
1384  /(nx^d*dxlevel(^d))) \}
1385 {^d&wholeextent(^d*2) = nx^d * floor(((xprobmax^d-xprobmin^d)*(1.0d0-writespshift(^d,2))) &
1386  /(nx^d*dxlevel(^d))) \}
1387 
1388 ! generate xml header
1389 write(qunit,'(a)')'<?xml version="1.0"?>'
1390 write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
1391 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1392 write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
1393  origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
1394  write(qunit,'(a)')'<FieldData>'
1395  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1396  'NumberOfTuples="1" format="ascii">'
1397  write(qunit,*) real(global_time*time_convert_factor)
1398  write(qunit,'(a)')'</DataArray>'
1399  write(qunit,'(a)')'</FieldData>'
1400 
1401 ! write the data from proc 0
1402 do morton_no=morton_start(0),morton_stop(0)
1403  if(.not. morton_aim(morton_no)) cycle
1404  igrid=sfc_to_igrid(morton_no)
1405  tree%node => igrid_to_node(igrid, 0)%node
1406  {^d& ig^d = tree%node%ig^d; }
1407  call calc_x(igrid,xc,xcc)
1408  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1409  ixc^l,ixcc^l,.true.)
1410  call write_vti(qunit,ixg^ll,ixc^l,ixcc^l,ig^d,&
1411  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1412 end do
1413 
1414 if(npe>1)then
1415  allocate(intstatus(mpi_status_size,1))
1416  do ipe=1, npe-1
1417  do morton_no=morton_start(ipe),morton_stop(ipe)
1418  if(.not. morton_aim(morton_no)) cycle
1419  itag=morton_no
1420  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1421  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1422  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1423  ig^d=ind_recv(4*^nd+^d);
1424  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1425  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1426  call write_vti(qunit,ixg^ll,ixrvc^l,ixrvcc^l,ig^d,&
1427  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1428  end do
1429  end do
1430 end if
1431 
1432 write(qunit,'(a)')'</ImageData>'
1433 write(qunit,'(a)')'</VTKFile>'
1434 close(qunit)
1435 if(npe>1) deallocate(intstatus)
1436 endif
1437 
1438 deallocate(morton_aim,morton_aim_p)
1439 if (npe>1) then
1440  call mpi_barrier(icomm,ierrmpi)
1441 endif
1442 
1443 end subroutine imagedatavtk_mpi
1444 !============================================================================
1445 subroutine punstructuredvtk_mpi(qunit)
1447 ! Write one pvtu and vtu files for each processor
1448 ! Otherwise like unstructuredvtk_mpi
1449 
1452 use mod_calculate_xw
1453 
1454 integer, intent(in) :: qunit
1455 !
1456 double precision, dimension(0:nw+nwauxio) :: normconv
1457 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1458 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1459 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1460 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1461 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
1462 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1463 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1464 character(len=1024) :: outfilehead
1465 integer :: nx^D,nxC^D,nc,np, igrid,ixC^L,ixCC^L,level,Morton_no
1466 character(len=80) :: pfilename
1467 integer :: filenr
1468 logical :: fileopen,conv_grid
1469 !----------------------------------------------------------------------------
1470 
1471 ! Write pvtu-file:
1472 if (mype==0) then
1473  call write_pvtu(qunit)
1474 endif
1475 ! Now write the Source files:
1476 
1477 inquire(qunit,opened=fileopen)
1478 if(.not.fileopen)then
1479  ! generate filename
1480  filenr=snapshotini
1481  if (autoconvert) filenr=snapshotnext
1482  ! Open the file for the header part
1483  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
1484  open(qunit,file=pfilename,status='unknown',form='formatted')
1485 endif
1486 ! generate xml header
1487 write(qunit,'(a)')'<?xml version="1.0"?>'
1488 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1489 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1490 write(qunit,'(a)')' <UnstructuredGrid>'
1491 write(qunit,'(a)')'<FieldData>'
1492 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1493  'NumberOfTuples="1" format="ascii">'
1494 write(qunit,*) real(global_time*time_convert_factor)
1495 write(qunit,'(a)')'</DataArray>'
1496 write(qunit,'(a)')'</FieldData>'
1497 
1498 call getheadernames(wnamei,xandwnamei,outfilehead)
1499 
1500 ! number of cells, number of corner points, per grid.
1501 nx^d=ixmhi^d-ixmlo^d+1;
1502 nxc^d=nx^d+1;
1503 nc={nx^d*}
1504 np={nxc^d*}
1505 
1506 ! Note: using the w_write, writelevel, writespshift
1507 ! we can clip parts of the grid away, select variables, levels etc.
1508 do level=levmin,levmax
1509  if (.not.writelevel(level)) cycle
1510  do morton_no=morton_start(mype),morton_stop(mype)
1511  igrid=sfc_to_igrid(morton_no)
1512  if (node(plevel_,igrid)/=level) cycle
1513 
1514  ! only output a grid when fully within clipped region selected
1515  ! by writespshift array
1516  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1517  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1518  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1519  if (.not.conv_grid) cycle
1520 
1521  call calc_x(igrid,xc,xcc)
1522  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1523  ixc^l,ixcc^l,.true.)
1524 
1525  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1526  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1527  end do ! Morton_no loop
1528 end do ! level loop
1529 
1530  write(qunit,'(a)')' </UnstructuredGrid>'
1531  write(qunit,'(a)')'</VTKFile>'
1532  close(qunit)
1533 
1534 if (npe>1) then
1535  call mpi_barrier(icomm,ierrmpi)
1536 endif
1537 end subroutine punstructuredvtk_mpi
1538 !============================================================================
1539 subroutine unstructuredvtk_mpi(qunit)
1541 ! output for vtu format to paraview, non-binary version output
1542 ! parallel, uses calc_grid to compute nwauxio variables
1543 ! allows renormalizing using convert factors
1544 ! allows skipping of w_write selected variables
1545 
1546 ! implementation such that length of ASCII output is identical when
1547 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1548 
1551 use mod_calculate_xw
1552 
1553 integer, intent(in) :: qunit
1554 
1555 double precision :: x_VTK(1:3)
1556 
1557 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1558 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1559 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1560 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1561 
1562 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1563 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1564 double precision, dimension(0:nw+nwauxio) :: normconv
1565 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1566 integer:: NumGridsOnLevel(1:nlevelshi)
1567 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,ix^D
1568 
1569 character(len=80):: filename
1570 integer :: filenr
1571 
1572 integer, allocatable :: intstatus(:,:)
1573 
1574 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1575 character(len=1024) :: outfilehead
1576 
1577 logical :: fileopen,conv_grid,cond_grid_recv
1578 integer :: itag,ipe,Morton_no,siz_ind
1579 integer :: ind_send(4*^nd),ind_recv(4*^nd)
1580 integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvC^L,ixrvCC^L
1581 !-----------------------------------------------------------------------------
1582 if (mype==0) then
1583  inquire(qunit,opened=fileopen)
1584  if(.not.fileopen)then
1585  ! generate filename
1586  filenr=snapshotini
1587  if (autoconvert) filenr=snapshotnext
1588  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
1589  ! Open the file for the header part
1590  open(qunit,file=filename,status='unknown',form='formatted')
1591  endif
1592  ! generate xml header
1593  write(qunit,'(a)')'<?xml version="1.0"?>'
1594  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1595  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1596  write(qunit,'(a)')'<UnstructuredGrid>'
1597  write(qunit,'(a)')'<FieldData>'
1598  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1599  'NumberOfTuples="1" format="ascii">'
1600  write(qunit,*) real(global_time*time_convert_factor)
1601  write(qunit,'(a)')'</DataArray>'
1602  write(qunit,'(a)')'</FieldData>'
1603 end if
1604 
1605 call getheadernames(wnamei,xandwnamei,outfilehead)
1606 ! number of cells, number of corner points, per grid.
1607 nx^d=ixmhi^d-ixmlo^d+1;
1608 nxc^d=nx^d+1;
1609 nc={nx^d*}
1610 np={nxc^d*}
1611 
1612 ! all slave processors send their minmal/maximal levels
1613 if (mype/=0) then
1614  if (morton_stop(mype)==0) call mpistop("nultag")
1615  itag=1000*morton_stop(mype)
1616  !print *,'ype,itag for levmin=',mype,itag,levmin
1617  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
1618  itag=2000*morton_stop(mype)
1619  !print *,'mype,itag for levmax=',mype,itag,levmax
1620  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
1621 end if
1622 
1623 
1624 ! Note: using the w_write, writelevel, writespshift
1625 ! we can clip parts of the grid away, select variables, levels etc.
1626 do level=levmin,levmax
1627  if (.not.writelevel(level)) cycle
1628  do morton_no=morton_start(mype),morton_stop(mype)
1629  igrid=sfc_to_igrid(morton_no)
1630  if (mype/=0)then
1631  itag=morton_no
1632  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
1633  itag=igrid
1634  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
1635  end if
1636  if (node(plevel_,igrid)/=level) cycle
1637 
1638  ! only output a grid when fully within clipped region selected
1639  ! by writespshift array
1640  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1641  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1642  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1643  if (mype/=0)then
1644  call mpi_send(conv_grid,1,mpi_logical,0,itag,icomm,ierrmpi)
1645  end if
1646  if (.not.conv_grid) cycle
1647 
1648  call calc_x(igrid,xc,xcc)
1649  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1650  ixc^l,ixcc^l,.true.)
1651 
1652  if (mype/=0) then
1653  itag=morton_no
1654  ind_send=(/ ixc^l,ixcc^l /)
1655  siz_ind=4*^nd
1656  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1657  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
1658 
1659  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1660  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
1661  itag=igrid
1662  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1663  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
1664  else
1665  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1666  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1667  end if
1668  end do ! Morton_no loop
1669 end do ! level loop
1670 
1671 
1672 if (mype==0) then
1673  allocate(intstatus(mpi_status_size,1))
1674  if(npe>1)then
1675  do ipe=1,npe-1
1676  itag=1000*morton_stop(ipe)
1677  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1678  !!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
1679  itag=2000*morton_stop(ipe)
1680  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1681  !!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
1682  do level=levmin_recv,levmax_recv
1683  if (.not.writelevel(level)) cycle
1684  do morton_no=morton_start(ipe),morton_stop(ipe)
1685  itag=morton_no
1686  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1687  itag=igrid_recv
1688  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1689  if (level_recv/=level) cycle
1690 
1691  call mpi_recv(cond_grid_recv,1,mpi_logical, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1692  if(.not.cond_grid_recv)cycle
1693 
1694  itag=morton_no
1695  siz_ind=4*^nd
1696  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1697  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1698  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1699  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
1700  ,icomm,intstatus(:,1),ierrmpi)
1701 
1702  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1703  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1704 
1705  itag=igrid_recv
1706  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1707  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1708  call write_vtk(qunit,ixg^ll,ixrvc^l,ixrvcc^l,igrid_recv,&
1709  nc,np,nx^d,nxc^d,normconv,wnamei,&
1710  xc_tmp_recv,xcc_tmp_recv,wc_tmp_recv,wcc_tmp_recv)
1711  enddo ! Morton_no loop
1712  enddo ! level loop
1713  enddo ! processor loop
1714  endif ! multiple processors
1715  write(qunit,'(a)')'</UnstructuredGrid>'
1716  write(qunit,'(a)')'</VTKFile>'
1717  close(qunit)
1718 endif
1719 
1720 if (npe>1) then
1721  call mpi_barrier(icomm,ierrmpi)
1722  if(mype==0)deallocate(intstatus)
1723 endif
1724 
1725 end subroutine unstructuredvtk_mpi
1726 !============================================================================
1727 subroutine write_vtk(qunit,ixI^L,ixC^L,ixCC^L,igrid,nc,np,nx^D,nxC^D,&
1728  normconv,wnamei,xC,xCC,wC,wCC)
1731 
1732 integer, intent(in) :: qunit
1733 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
1734 integer, intent(in) :: igrid,nc,np,nx^D,nxC^D
1735 double precision, intent(in) :: normconv(0:nw+nwauxio)
1736 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
1737 
1738 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1739 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1740 
1741 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
1742 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
1743 
1744 double precision :: x_VTK(1:3)
1745 integer :: iw,ix^D,icel,VTK_type
1746 !----------------------------------------------------------------------------
1747 
1748 select case(convert_type)
1749  case('vtumpi','pvtumpi')
1750  ! we write out every grid as one VTK PIECE
1751  write(qunit,'(a,i7,a,i7,a)') &
1752  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1753  write(qunit,'(a)')'<PointData>'
1754  do iw=1,nw+nwauxio
1755  if(iw<=nw) then
1756  if(.not.w_write(iw)) cycle
1757  endif
1758 
1759  write(qunit,'(a,a,a)')&
1760  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1761  write(qunit,'(200(1pe14.6))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1762  write(qunit,'(a)')'</DataArray>'
1763  enddo
1764  write(qunit,'(a)')'</PointData>'
1765 
1766  write(qunit,'(a)')'<Points>'
1767  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
1768  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1769  {do ix^db=ixcmin^db,ixcmax^db \}
1770  x_vtk(1:3)=zero;
1771  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
1772  write(qunit,'(3(1pe14.6))') x_vtk
1773  {end do \}
1774  write(qunit,'(a)')'</DataArray>'
1775  write(qunit,'(a)')'</Points>'
1776 
1777  case('vtuCCmpi','pvtuCCmpi')
1778  ! we write out every grid as one VTK PIECE
1779  write(qunit,'(a,i7,a,i7,a)') &
1780  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1781  write(qunit,'(a)')'<CellData>'
1782  do iw=1,nw+nwauxio
1783  if(iw<=nw) then
1784  if(.not.w_write(iw)) cycle
1785  endif
1786  write(qunit,'(a,a,a)')&
1787  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1788  write(qunit,'(200(1pe14.6))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1789  write(qunit,'(a)')'</DataArray>'
1790  enddo
1791  write(qunit,'(a)')'</CellData>'
1792 
1793  write(qunit,'(a)')'<Points>'
1794  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
1795  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1796  {do ix^db=ixcmin^db,ixcmax^db \}
1797  x_vtk(1:3)=zero;
1798  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
1799  write(qunit,'(3(1pe14.6))') x_vtk
1800  {end do \}
1801  write(qunit,'(a)')'</DataArray>'
1802  write(qunit,'(a)')'</Points>'
1803 end select
1804 
1805 write(qunit,'(a)')'<Cells>'
1806 
1807 ! connectivity part
1808 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
1809 call save_connvtk(qunit,igrid)
1810 write(qunit,'(a)')'</DataArray>'
1811 
1812 ! offsets data array
1813 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
1814 do icel=1,nc
1815  write(qunit,'(i7)') icel*(2**^nd)
1816 end do
1817 write(qunit,'(a)')'</DataArray>'
1818 
1819 ! VTK cell type data array
1820 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
1821 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
1822 {^ifoned vtk_type=3 \}
1823 {^iftwod vtk_type=8 \}
1824 {^ifthreed vtk_type=11 \}
1825 do icel=1,nc
1826  write(qunit,'(i2)') vtk_type
1827 enddo
1828 write(qunit,'(a)')'</DataArray>'
1829 
1830 write(qunit,'(a)')'</Cells>'
1831 
1832 write(qunit,'(a)')'</Piece>'
1833 
1834 end subroutine write_vtk
1835 !============================================================================
1836 subroutine write_vti(qunit,ixI^L,ixC^L,ixCC^L,ig^D,nx^D,&
1837  normconv,wnamei,wC,wCC)
1839 
1840 integer, intent(in) :: qunit
1841 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
1842 integer, intent(in) :: ig^D,nx^D
1843 double precision, intent(in) :: normconv(0:nw+nwauxio)
1844 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
1845 
1846 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
1847 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
1848 
1849 integer :: iw,ix^D
1850 integer :: extent(1:6)
1851 !----------------------------------------------------------------------------
1852 
1853 extent = 0
1854 {^d& extent(^d*2-1) = (ig^d-1) * nx^d; }
1855 {^d& extent(^d*2) = (ig^d) * nx^d; }
1856 
1857 
1858 select case(convert_type)
1859  case('vtimpi','pvtimpi')
1860  ! we write out every grid as one VTK PIECE
1861  write(qunit,'(a,6(i10),a)') &
1862  '<Piece Extent="',extent,'">'
1863  write(qunit,'(a)')'<PointData>'
1864  do iw=1,nw+nwauxio
1865  if(iw<=nw) then
1866  if(.not.w_write(iw)) cycle
1867  endif
1868 
1869  write(qunit,'(a,a,a)')&
1870  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1871  write(qunit,'(200(1pe20.12))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1872  write(qunit,'(a)')'</DataArray>'
1873  enddo
1874  write(qunit,'(a)')'</PointData>'
1875 
1876  case('vtiCCmpi','pvtiCCmpi')
1877  ! we write out every grid as one VTK PIECE
1878  write(qunit,'(a,6(i10),a)') &
1879  '<Piece Extent="',extent,'">'
1880  write(qunit,'(a)')'<CellData>'
1881  do iw=1,nw+nwauxio
1882  if(iw<=nw) then
1883  if(.not.w_write(iw)) cycle
1884  endif
1885  write(qunit,'(a,a,a)')&
1886  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
1887  write(qunit,'(200(1pe20.12))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1888  write(qunit,'(a)')'</DataArray>'
1889  enddo
1890  write(qunit,'(a)')'</CellData>'
1891 end select
1892 
1893 write(qunit,'(a)')'</Piece>'
1894 
1895 end subroutine write_vti
1896 !=============================================================================
1897 subroutine write_pvtu(qunit)
1900 use mod_calculate_xw
1901 
1902 integer, intent(in) :: qunit
1903 
1904 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio),outtype
1905 character(len=1024) :: outfilehead
1906 character(len=80) :: filename,pfilename
1907 integer :: filenr,iw,ipe,iscalars
1908 logical :: fileopen
1909 
1910 select case(convert_type)
1911 case('pvtumpi','pvtuBmpi')
1912  outtype="PPointData"
1913 case('pvtuCCmpi','pvtuBCCmpi')
1914  outtype="PCellData"
1915 end select
1916 inquire(qunit,opened=fileopen)
1917 if(.not.fileopen)then
1918  ! generate filename
1919  filenr=snapshotini
1920  if (autoconvert) filenr=snapshotnext
1921  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".pvtu"
1922  ! Open the file
1923  open(qunit,file=filename,status='unknown',form='formatted')
1924 endif
1925 
1926 call getheadernames(wnamei,xandwnamei,outfilehead)
1927 
1928 ! Get the default selection:
1929 iscalars=1
1930 do iw=nw,1, -1
1931  if (w_write(iw)) iscalars=iw
1932 end do
1933 
1934 
1935 ! generate xml header
1936 write(qunit,'(a)')'<?xml version="1.0"?>'
1937 write(qunit,'(a)',advance='no') '<VTKFile type="PUnstructuredGrid"'
1938 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1939 write(qunit,'(a)')' <PUnstructuredGrid GhostLevel="0">'
1940 ! Either celldata or pointdata:
1941 write(qunit,'(a,a,a,a,a)')&
1942  ' <',trim(outtype),' Scalars="',trim(wnamei(iscalars))//'">'
1943 do iw=1,nw
1944  if(.not.w_write(iw))cycle
1945  write(qunit,'(a,a,a)')&
1946  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
1947 end do
1948 do iw=nw+1,nw+nwauxio
1949  write(qunit,'(a,a,a)')&
1950  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
1951 end do
1952 write(qunit,'(a,a,a)')' </',trim(outtype),'>'
1953 
1954 write(qunit,'(a)')' <PPoints>'
1955 write(qunit,'(a)')' <PDataArray type="Float32" NumberOfComponents="3"/>'
1956 write(qunit,'(a)')' </PPoints>'
1957 
1958 do ipe=0,npe-1
1959  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename(&
1960  index(base_filename, '/', back = .true.)+1:&
1961  len(base_filename))),filenr,"p",&
1962  ipe,".vtu"
1963  write(qunit,'(a,a,a)')' <Piece Source="',trim(pfilename),'"/>'
1964 end do
1965 write(qunit,'(a)')' </PUnstructuredGrid>'
1966 write(qunit,'(a)')'</VTKFile>'
1967 
1968 close(qunit)
1969 
1970 end subroutine write_pvtu
1971 !=============================================================================
1972 subroutine tecplot_mpi(qunit)
1974 ! output for tecplot (ASCII format)
1975 ! parallel, uses calc_grid to compute nwauxio variables
1976 ! allows renormalizing using convert factors
1977 
1978 ! the current implementation is such that tecplotmpi and tecplotCCmpi will
1979 ! create different length output ASCII files when used on 1 versus multiple CPUs
1980 ! in fact, on 1 CPU, there will be as many zones as there are levels
1981 ! on multiple CPUs, there will be a number of zones up to the number of
1982 ! levels times the number of CPUs (can be less, when some level not on a CPU)
1983 
1986 use mod_calculate_xw
1987 
1988 integer, intent(in) :: qunit
1989 
1990 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
1991 integer:: NumGridsOnLevel(1:nlevelshi)
1992 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
1993 integer :: nodesonlevelmype,elemsonlevelmype
1994 
1995 integer :: nodes, elems
1996 
1997 integer, allocatable :: intstatus(:,:)
1998 
1999 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
2000 
2001 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
2002 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
2003 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2004 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2005 
2006 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
2007 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
2008 double precision, dimension(0:nw+nwauxio) :: normconv
2009 logical :: fileopen,first
2010 integer :: itag,Morton_no,ipe,levmin_recv,levmax_recv,igrid_recv,level_recv
2011 integer :: ixrvC^L,ixrvCC^L
2012 integer :: ind_send(2*^nd),ind_recv(2*^nd),siz_ind,igonlevel_recv
2013 integer :: NumGridsOnLevel_mype(1:nlevelshi,0:npe-1)
2014 character(len=80) :: filename
2015 integer :: filenr
2016 character(len=1024) :: tecplothead
2017 
2018 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2019 character(len=1024) :: outfilehead
2020 !-----------------------------------------------------------------------------
2021 if(nw/=count(w_write(1:nw)))then
2022  if(mype==0) print *,'tecplot_mpi does not use w_write=F'
2023  call mpistop('w_write, tecplot')
2024 end if
2025 
2026 if(nocartesian)then
2027  if(mype==0) print *,'tecplot_mpi with nocartesian'
2028 endif
2029 
2030 master_cpu_open : if (mype == 0) then
2031  inquire(qunit,opened=fileopen)
2032  if (.not.fileopen) then
2033  ! generate filename
2034  filenr=snapshotini
2035  if (autoconvert) filenr=snapshotnext
2036  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
2037  open(qunit,file=filename,status='unknown')
2038  end if
2039 
2040  call getheadernames(wnamei,xandwnamei,outfilehead)
2041 
2042  write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
2043  write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
2044 end if master_cpu_open
2045 
2046 
2047 ! determine overall number of grids per level, and the same info per CPU
2048 numgridsonlevel(1:nlevelshi)=0
2049 do level=levmin,levmax
2050  numgridsonlevel(level)=0
2051  do morton_no=morton_start(mype),morton_stop(mype)
2052  igrid = sfc_to_igrid(morton_no)
2053  if (node(plevel_,igrid)/=level) cycle
2054  numgridsonlevel(level)=numgridsonlevel(level)+1
2055  end do
2056  numgridsonlevel_mype(level,0:npe-1)=0
2057  numgridsonlevel_mype(level,mype) = numgridsonlevel(level)
2058  call mpi_allreduce(mpi_in_place,numgridsonlevel_mype(level,0:npe-1),npe,mpi_integer,&
2059  mpi_max,icomm,ierrmpi)
2060  call mpi_allreduce(mpi_in_place,numgridsonlevel(level),1,mpi_integer,mpi_sum, &
2061  icomm,ierrmpi)
2062 end do
2063 
2064 
2065 !!do level=levmin,levmax
2066 !! print *,'mype, level en NumGridsOnLevel_mype(level,0:npe-1)=', &
2067 !! mype,level,NumGridsOnLevel_mype(level,0:npe-1)
2068 !! print *,'mype, level en NumGridsOnLevel(level)=', &
2069 !! mype,level,NumGridsOnLevel(level)
2070 !!enddo
2071 
2072 
2073 nx^d=ixmhi^d-ixmlo^d+1;
2074 nxc^d=nx^d+1;
2075 
2076 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
2077 
2078 {^ifoned
2079 if(convert_type=='teclinempi') then
2080  nodes=0
2081  elems=0
2082  do level=levmin,levmax
2083  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
2084  elems=elems + numgridsonlevel(level)*{nx^d*}
2085  enddo
2086 
2087  if (mype==0) write(qunit,"(a,i7,a,1pe12.5,a)") &
2088  'ZONE T="all levels", I=',elems, &
2089  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
2090 
2091  igonlevel=0
2092  do morton_no=morton_start(mype),morton_stop(mype)
2093  igrid = sfc_to_igrid(morton_no)
2094  call calc_x(igrid,xc,xcc)
2095  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
2096  if (mype==0) then
2097  {do ix^db=ixccmin^db,ixccmax^db\}
2098  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
2099  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2100  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2101  {end do\}
2102  else if (mype/=0) then
2103  itag=morton_no
2104  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2105  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
2106  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2107  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
2108  end if
2109  enddo
2110  if (mype==0) then
2111  do ipe=1,npe-1
2112  do morton_no=morton_start(ipe),morton_stop(ipe)
2113  itag=morton_no
2114  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2115  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,&
2116  itag,icomm,intstatus(:,1),ierrmpi)
2117  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,&
2118  icomm,intstatus(:,1),ierrmpi)
2119  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,&
2120  icomm,intstatus(:,1),ierrmpi)
2121  {do ix^db=ixccmin^db,ixccmax^db\}
2122  x_tec(1:ndim)=xcc_tmp_recv(ix^d,1:ndim)*normconv(0)
2123  w_tec(1:nw+nwauxio)=wcc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2124  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2125  {end do\}
2126  end do
2127  end do
2128  close(qunit)
2129  end if
2130 else
2131 }
2132 
2133 if (mype/=0) then
2134  itag=1000*morton_stop(mype)
2135  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
2136  itag=2000*morton_stop(mype)
2137  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
2138 end if
2139 
2140 do level=levmin,levmax
2141  nodesonlevelmype=numgridsonlevel_mype(level,mype)*{nxc^d*}
2142  elemsonlevelmype=numgridsonlevel_mype(level,mype)*{nx^d*}
2143  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2144  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2145  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
2146  ! with the AMR grid LEVEL. Other options would be
2147  ! let each grid define a zone: inefficient for TECPLOT internal workings
2148  ! hence not implemented
2149  ! let entire octree define 1 zone: no difference in interpolation
2150  ! properties across TECPLOT zones detected as yet, hence not done
2151  select case(convert_type)
2152  case('tecplotmpi')
2153  ! in this option, we store the corner coordinates, as well as the corner
2154  ! values of all variables (obtained by averaging). This allows POINT packaging,
2155  ! and thus we can save full grid info by using one call to calc_grid
2156  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2157  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2158  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2159  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2160  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2161  do morton_no=morton_start(mype),morton_stop(mype)
2162  igrid = sfc_to_igrid(morton_no)
2163  if (mype/=0)then
2164  itag=morton_no
2165  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2166  itag=igrid
2167  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2168  end if
2169  if (node(plevel_,igrid)/=level) cycle
2170  call calc_x(igrid,xc,xcc)
2171  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2172  ixc^l,ixcc^l,.true.)
2173  if (mype/=0) then
2174  itag=morton_no
2175  ind_send=(/ ixc^l /)
2176  siz_ind=2*^nd
2177  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2178  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2179 
2180  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
2181  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2182  else
2183  {do ix^db=ixcmin^db,ixcmax^db\}
2184  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
2185  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2186  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2187  {end do\}
2188  end if
2189  enddo
2190 
2191  case('tecplotCCmpi')
2192  ! in this option, we store the corner coordinates, and the cell center
2193  ! values of all variables. Due to this mix of corner/cell center, we must
2194  ! use BLOCK packaging, and thus we have enormous overhead by using
2195  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2196  ! and cell center values per dimension, per variable
2197  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2198  if(nw+nwauxio==1)then
2199  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2200  ! and just set [ndim+1]
2201  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2202  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2203  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2204  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2205  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2206  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2207  else
2208  if(ndim+nw+nwauxio<10) then
2209  ! difference only in length of integer format specification for ndim+nw+nwauxio
2210  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2211  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2212  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2213  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2214  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2215  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2216  else
2217  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2218  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2219  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2220  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2221  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2222  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2223  endif
2224  endif
2225 
2226  do idim=1,ndim
2227  first=(idim==1)
2228  do morton_no=morton_start(mype),morton_stop(mype)
2229  igrid = sfc_to_igrid(morton_no)
2230  if (mype/=0)then
2231  itag=morton_no*idim
2232  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2233  itag=igrid*idim
2234  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2235  end if
2236  if (node(plevel_,igrid)/=level) cycle
2237 
2238  call calc_x(igrid,xc,xcc)
2239  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2240  ixc^l,ixcc^l,first)
2241  if (mype/=0)then
2242  ind_send=(/ ixc^l /)
2243  siz_ind=2*^nd
2244  itag=igrid*idim
2245  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2246  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2247  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2248  else
2249  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
2250  end if
2251  enddo
2252  enddo
2253 
2254  do iw=1,nw+nwauxio
2255  do morton_no=morton_start(mype),morton_stop(mype)
2256  igrid = sfc_to_igrid(morton_no)
2257  if (mype/=0)then
2258  itag=morton_no*(ndim+iw)
2259  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2260  itag=igrid*(ndim+iw)
2261  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2262  end if
2263  if (node(plevel_,igrid)/=level) cycle
2264 
2265  call calc_x(igrid,xc,xcc)
2266  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2267  ixc^l,ixcc^l,.true.)
2268 
2269  if (mype/=0)then
2270  ind_send=(/ ixcc^l /)
2271  siz_ind=2*^nd
2272  itag=igrid*(ndim+iw)
2273  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2274  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2275  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2276  else
2277  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
2278  endif
2279  enddo
2280  enddo
2281  case default
2282  call mpistop('no such tecplot type')
2283  end select
2284 
2285 
2286  igonlevel=0
2287  do morton_no=morton_start(mype),morton_stop(mype)
2288  igrid = sfc_to_igrid(morton_no)
2289  if (mype/=0)then
2290  itag=morton_no
2291  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2292  itag=igrid
2293  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2294  end if
2295  if (node(plevel_,igrid)/=level) cycle
2296 
2297  igonlevel=igonlevel+1
2298  if (mype/=0)then
2299  itag=igrid
2300  call mpi_send(igonlevel,1,mpi_integer, 0,itag,icomm,ierrmpi)
2301  end if
2302  if(mype==0)then
2303  call save_conntec(qunit,igrid,igonlevel)
2304  endif
2305  end do
2306 end do
2307 
2308 if (mype==0) then
2309  if (npe>1) then
2310  do ipe=1,npe-1
2311  itag=1000*morton_stop(ipe)
2312  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2313  itag=2000*morton_stop(ipe)
2314  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2315  do level=levmin_recv,levmax_recv
2316  nodesonlevelmype=numgridsonlevel_mype(level,ipe)*{nxc^d*}
2317  elemsonlevelmype=numgridsonlevel_mype(level,ipe)*{nx^d*}
2318  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2319  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2320  select case(convert_type)
2321  case('tecplotmpi')
2322  ! in this option, we store the corner coordinates, as well as the corner
2323  ! values of all variables (obtained by averaging). This allows POINT packaging,
2324  ! and thus we can save full grid info by using one call to calc_grid
2325  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2326  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2327  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2328  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2329  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2330  do morton_no=morton_start(ipe),morton_stop(ipe)
2331  itag=morton_no
2332  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2333  itag=igrid_recv
2334  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2335  if (level_recv/=level) cycle
2336 
2337  itag=morton_no
2338  siz_ind=2*^nd
2339  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,&
2340  icomm,intstatus(:,1),ierrmpi)
2341  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2342  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2343  ,icomm,intstatus(:,1),ierrmpi)
2344 
2345  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,&
2346  icomm,intstatus(:,1),ierrmpi)
2347  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,&
2348  icomm,intstatus(:,1),ierrmpi)
2349 
2350  {do ix^db=ixrvcmin^db,ixrvcmax^db\}
2351  x_tec(1:ndim)=xc_tmp_recv(ix^d,1:ndim)*normconv(0)
2352  w_tec(1:nw+nwauxio)=wc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2353  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2354  {end do\}
2355  end do
2356  case('tecplotCCmpi')
2357  ! in this option, we store the corner coordinates, and the cell center
2358  ! values of all variables. Due to this mix of corner/cell center, we must
2359  ! use BLOCK packaging, and thus we have enormous overhead by using
2360  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2361  ! and cell center values per dimension, per variable
2362  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2363  if(nw+nwauxio==1)then
2364  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2365  ! and just set [ndim+1]
2366  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2367  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2368  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2369  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2370  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2371  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2372  else
2373  if(ndim+nw+nwauxio<10) then
2374  ! difference only in length of integer format specification for ndim+nw+nwauxio
2375  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2376  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2377  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2378  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2379  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2380  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2381  else
2382  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2383  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2384  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2385  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2386  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2387  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2388  endif
2389  endif
2390 
2391  do idim=1,ndim
2392  do morton_no=morton_start(ipe),morton_stop(ipe)
2393  itag=morton_no*idim
2394  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2395  itag=igrid_recv*idim
2396  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2397  if (level_recv/=level) cycle
2398 
2399  siz_ind=2*^nd
2400  itag=igrid_recv*idim
2401  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2402  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2403  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2404  ,icomm,intstatus(:,1),ierrmpi)
2405  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2406  write(qunit,fmt="(100(e14.6))") xc_tmp_recv(ixrvc^s,idim)*normconv(0)
2407  end do
2408  end do
2409 
2410  do iw=1,nw+nwauxio
2411  do morton_no=morton_start(ipe),morton_stop(ipe)
2412  itag=morton_no*(ndim+iw)
2413  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2414  itag=igrid_recv*(ndim+iw)
2415  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2416  if (level_recv/=level) cycle
2417 
2418  siz_ind=2*^nd
2419  itag=igrid_recv*(ndim+iw)
2420  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2421  ixrvccmin^d=ind_recv(^d);ixrvccmax^d=ind_recv(^nd+^d);
2422  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2423  ,icomm,intstatus(:,1),ierrmpi)
2424  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2425  write(qunit,fmt="(100(e14.6))") wcc_tmp_recv(ixrvcc^s,iw)*normconv(iw)
2426  enddo
2427  enddo
2428  case default
2429  call mpistop('no such tecplot type')
2430  end select
2431 
2432  do morton_no=morton_start(ipe),morton_stop(ipe)
2433  itag=morton_no
2434  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2435  itag=igrid_recv
2436  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2437  if (level_recv/=level) cycle
2438 
2439  itag=igrid_recv
2440  call mpi_recv(igonlevel_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2441  call save_conntec(qunit,igrid_recv,igonlevel_recv)
2442  end do ! morton loop
2443  end do ! level loop
2444  end do ! ipe loop
2445  end if ! npe>1 if
2446 end if ! mype=0 if
2447 {^ifoned endif}
2448 
2449 if (npe>1) then
2450  call mpi_barrier(icomm,ierrmpi)
2451  if(mype==0)deallocate(intstatus)
2452 endif
2453 
2454 end subroutine tecplot_mpi
2455 !=============================================================================
2456 subroutine punstructuredvtkb_mpi(qunit)
2458 ! Write one pvtu and vtu files for each processor
2459 ! Otherwise like unstructuredvtk_mpi
2460 ! output for vtu format to paraview, binary version output
2461 ! uses calc_grid to compute nwauxio variables
2462 ! allows renormalizing using convert factors
2463 
2466 use mod_calculate_xw
2467 
2468 integer, intent(in) :: qunit
2469 
2470 double precision :: x_VTK(1:3)
2471 
2472 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
2473 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
2474 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2475 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2476 
2477 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
2478 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
2479 
2480 integer :: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,Morton_no
2481 integer :: NumGridsOnLevel(1:nlevelshi)
2482 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
2483 double precision :: normconv(0:nw+nwauxio)
2484 character(len=80) :: pfilename
2485 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2486 character(len=1024) :: outfilehead
2487 
2488 integer*8 :: offset
2489 integer:: recsep,k,iw,filenr
2490 integer:: length,lengthcc,offset_points,offset_cells, &
2491  length_coords,length_conn,length_offsets
2492 character:: buf
2493 character(len=6):: bufform
2494 
2495 logical :: fileopen
2496 !-----------------------------------------------------------------------------
2497 
2498 ! Write pvtu-file:
2499 if (mype==0) then
2500  call write_pvtu(qunit)
2501 endif
2502 ! Now write the Source files:
2503 inquire(qunit,opened=fileopen)
2504 if(.not.fileopen)then
2505  ! generate filename
2506  filenr=snapshotnext-1
2507  if (autoconvert) filenr=snapshotnext
2508  ! Open the file for the header part
2509  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
2510  open(qunit,file=pfilename,status='unknown',form='formatted')
2511 endif
2512 ! generate xml header
2513 write(qunit,'(a)')'<?xml version="1.0"?>'
2514 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2515 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2516 write(qunit,'(a)')' <UnstructuredGrid>'
2517 write(qunit,'(a)')'<FieldData>'
2518 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2519  'NumberOfTuples="1" format="ascii">'
2520 write(qunit,*) real(global_time*time_convert_factor)
2521 write(qunit,'(a)')'</DataArray>'
2522 write(qunit,'(a)')'</FieldData>'
2523 
2524 
2525 offset=0
2526 recsep=4
2527 
2528 call getheadernames(wnamei,xandwnamei,outfilehead)
2529 
2530 ! number of cells, number of corner points, per grid.
2531 nx^d=ixmhi^d-ixmlo^d+1;
2532 nxc^d=nx^d+1;
2533 nc={nx^d*}
2534 np={nxc^d*}
2535 
2536 length=np*size_real
2537 lengthcc=nc*size_real
2538 
2539 length_coords=3*length
2540 length_conn=2**^nd*size_int*nc
2541 length_offsets=nc*size_int
2542 
2543 ! Note: using the w_write, writelevel, writespshift
2544 ! we can clip parts of the grid away, select variables, levels etc.
2545 do level=levmin,levmax
2546  if (writelevel(level)) then
2547  do morton_no=morton_start(mype),morton_stop(mype)
2548  igrid=sfc_to_igrid(morton_no)
2549  if (node(plevel_,igrid)/=level) cycle
2550  ! only output a grid when fully within clipped region selected
2551  ! by writespshift array
2552  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2553  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2554  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2555  select case(convert_type)
2556  case('pvtuBmpi')
2557  ! we write out every grid as one VTK PIECE
2558  write(qunit,'(a,i7,a,i7,a)') &
2559  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2560  write(qunit,'(a)')'<PointData>'
2561  do iw=1,nw
2562  if(.not.w_write(iw))cycle
2563 
2564  write(qunit,'(a,a,a,i16,a)')&
2565  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2566  '" format="appended" offset="',offset,'">'
2567  write(qunit,'(a)')'</DataArray>'
2568  offset=offset+length+size_int
2569  enddo
2570  do iw=nw+1,nw+nwauxio
2571 
2572  write(qunit,'(a,a,a,i16,a)')&
2573  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2574  '" format="appended" offset="',offset,'">'
2575  write(qunit,'(a)')'</DataArray>'
2576  offset=offset+length+size_int
2577  enddo
2578  write(qunit,'(a)')'</PointData>'
2579 
2580  write(qunit,'(a)')'<Points>'
2581  write(qunit,'(a,i16,a)') &
2582  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2583  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2584  offset=offset+length_coords+size_int
2585  write(qunit,'(a)')'</Points>'
2586  case('pvtuBCCmpi')
2587  ! we write out every grid as one VTK PIECE
2588  write(qunit,'(a,i7,a,i7,a)') &
2589  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2590  write(qunit,'(a)')'<CellData>'
2591  do iw=1,nw
2592  if(.not.w_write(iw))cycle
2593 
2594  write(qunit,'(a,a,a,i16,a)')&
2595  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2596  '" format="appended" offset="',offset,'">'
2597  write(qunit,'(a)')'</DataArray>'
2598  offset=offset+lengthcc+size_int
2599  enddo
2600  do iw=nw+1,nw+nwauxio
2601 
2602  write(qunit,'(a,a,a,i16,a)')&
2603  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2604  '" format="appended" offset="',offset,'">'
2605  write(qunit,'(a)')'</DataArray>'
2606  offset=offset+lengthcc+size_int
2607  enddo
2608  write(qunit,'(a)')'</CellData>'
2609 
2610  write(qunit,'(a)')'<Points>'
2611  write(qunit,'(a,i16,a)') &
2612  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
2613  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2614  offset=offset+length_coords+size_int
2615  write(qunit,'(a)')'</Points>'
2616  end select
2617 
2618 
2619  write(qunit,'(a)')'<Cells>'
2620 
2621  ! connectivity part
2622  write(qunit,'(a,i16,a)')&
2623  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
2624  offset=offset+length_conn+size_int
2625 
2626  ! offsets data array
2627  write(qunit,'(a,i16,a)') &
2628  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
2629  offset=offset+length_offsets+size_int
2630 
2631  ! VTK cell type data array
2632  write(qunit,'(a,i16,a)') &
2633  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
2634  offset=offset+size_int+nc*size_int
2635 
2636  write(qunit,'(a)')'</Cells>'
2637 
2638  write(qunit,'(a)')'</Piece>'
2639  endif
2640  enddo
2641  endif
2642 enddo
2643 
2644 write(qunit,'(a)')'</UnstructuredGrid>'
2645 write(qunit,'(a)')'<AppendedData encoding="raw">'
2646 
2647 close(qunit)
2648 ! next to make gfortran compiler happy, as it does not know
2649 ! form='binary' and produces error on compilation
2650 !bufform='binary'
2651 !open(qunit,file=pfilename,form=bufform,position='append')
2652 !This should in principle do also for gfortran (tested with gfortran 4.6.0 and Intel 11.1):
2653 open(qunit,file=pfilename,access='stream',form='unformatted',position='append')
2654 buf='_'
2655 write(qunit) trim(buf)
2656 
2657 do level=levmin,levmax
2658  if (writelevel(level)) then
2659  do morton_no=morton_start(mype),morton_stop(mype)
2660  igrid=sfc_to_igrid(morton_no)
2661  if (node(plevel_,igrid)/=level) cycle
2662  ! only output a grid when fully within clipped region selected
2663  ! by writespshift array
2664  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2665  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2666  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2667  call calc_x(igrid,xc,xcc)
2668  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2669  ixc^l,ixcc^l,.true.)
2670  do iw=1,nw
2671  if(.not.w_write(iw))cycle
2672  select case(convert_type)
2673  case('pvtuBmpi')
2674  write(qunit) length
2675  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
2676  case('pvtuBCCmpi')
2677  write(qunit) lengthcc
2678  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
2679  end select
2680  enddo
2681  do iw=nw+1,nw+nwauxio
2682  select case(convert_type)
2683  case('pvtuBmpi')
2684  write(qunit) length
2685  write(qunit) {(|}real(wC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixcmin^d,ixcmax^d)}
2686  case('pvtuBCCmpi')
2687  write(qunit) lengthcc
2688  write(qunit) {(|}real(wCC_TMP(ix^D,iw)*normconv(iw)),{ix^D=ixccmin^d,ixccmax^d)}
2689  end select
2690  enddo
2691 
2692  write(qunit) length_coords
2693  {do ix^db=ixcmin^db,ixcmax^db \}
2694  x_vtk(1:3)=zero;
2695  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
2696  do k=1,3
2697  write(qunit) real(x_vtk(k))
2698  end do
2699  {end do \}
2700 
2701  write(qunit) length_conn
2702  {do ix^db=1,nx^db\}
2703  {^ifoned write(qunit)ix1-1,ix1 \}
2704  {^iftwod
2705  write(qunit)(ix2-1)*nxc1+ix1-1, &
2706  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
2707  \}
2708  {^ifthreed
2709  write(qunit)&
2710  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
2711  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2712  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2713  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
2714  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
2715  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
2716  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
2717  ix3*nxc2*nxc1+ ix2*nxc1+ix1
2718  \}
2719  {end do\}
2720 
2721  write(qunit) length_offsets
2722  do icel=1,nc
2723  write(qunit) icel*(2**^nd)
2724  end do
2725 
2726 
2727  {^ifoned vtk_type=3 \}
2728  {^iftwod vtk_type=8 \}
2729  {^ifthreed vtk_type=11 \}
2730  write(qunit) size_int*nc
2731  do icel=1,nc
2732  write(qunit) vtk_type
2733  enddo
2734  endif
2735  end do
2736  endif
2737 end do
2738 
2739 close(qunit)
2740 open(qunit,file=pfilename,status='unknown',form='formatted',position='append')
2741 
2742 write(qunit,'(a)')'</AppendedData>'
2743 write(qunit,'(a)')'</VTKFile>'
2744 close(qunit)
2745 
2746 end subroutine punstructuredvtkb_mpi
2747 !=============================================================================
subroutine write_pvtu(qunit)
Definition: convert.t:1898
This module contains definitions of global parameters and variables and some generic functions/subrou...
subroutine get_euv_image(qunit)
subroutine imagedatavtk_mpi(qunit)
Definition: convert.t:1272
update ghost cells of all blocks including physical boundaries
logical, dimension(:), allocatable w_write
subroutine onegrid(qunit)
Definition: convert.t:256
logical saveprim
If true, convert from conservative to primitive variables in output.
integer max_blocks
The maximum number of grid blocks in a processor.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
subroutine oneblock(qunit)
Definition: convert.t:58
integer type_block_wc_io
MPI type for IO: cell corner (wc) or cell center (wcc) variables.
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
double precision, dimension(ndim, 2) writespshift
double precision time_convert_factor
Conversion factor for time unit.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter plevel_
subroutine write_vti(qunit, ixIL, ixCL, ixCCL, igD, nxD, normconv, wnamei, wC, wCC)
Definition: convert.t:1838
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
Module with basic grid data structures.
Definition: mod_forest.t:2
subroutine calc_grid(qunit, igrid, xC, xCC, xC_TMP, xCC_TMP, wC_TMP, wCC_TMP, normconv, ixCL, ixCCL, first)
Compute both corner as well as cell-centered values for output.
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl...
Definition: mod_physics.t:22
integer npe
The number of MPI tasks.
subroutine save_conntec(qunit, igrid, igonlevel)
Definition: convert.t:573
logical b0field
split magnetic field as background B0 field
logical, dimension(:), allocatable writelevel
subroutine unstructuredvtk(qunit)
Definition: convert.t:648
integer function nodenumbertec1d(i1, nx1, ig, igrid)
Definition: convert.t:621
Pointer to a tree_node.
Definition: mod_forest.t:6
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
subroutine calc_x(igrid, xC, xCC)
computes cell corner (xC) and cell center (xCC) coordinates
procedure(aux_output), pointer usr_aux_output
Module with all the methods that users can customize in AMRVAC.
integer, dimension(:), allocatable, save sfc_to_igrid
Go from a Morton number to an igrid index (for a single processor)
Definition: mod_forest.t:53
integer type_block_io
MPI type for IO: block excluding ghost cells.
integer, parameter unitconvert
subroutine punstructuredvtkb_mpi(qunit)
Definition: convert.t:2457
logical nocartesian
IO switches for conversion.
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
subroutine unstructuredvtk_mpi(qunit)
Definition: convert.t:1540
subroutine tecplot(qunit)
Definition: convert.t:366
logical image
output image
subroutine write_vtk(qunit, ixIL, ixCL, ixCCL, igrid, nc, np, nxD, nxCD, normconv, wnamei, xC, xCC, wC, wCC)
Definition: convert.t:1729
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
subroutine punstructuredvtk_mpi(qunit)
Definition: convert.t:1446
integer, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
subroutine unstructuredvtkb(qunit)
Definition: convert.t:823
integer, parameter unitterm
Unit for standard output.
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
integer mype
The rank of the current MPI task.
double precision global_time
The global simulation time.
subroutine resettree_convert
Force the tree to desired level(s) from level_io(_min/_max) used for conversion to vtk output...
Definition: amrgrid.t:77
subroutine get_euv_spectra(qunit)
subroutine save_connvtk(qunit, igrid)
Definition: convert.t:1237
subroutine tecplot_mpi(qunit)
Definition: convert.t:1973
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:208
integer function nodenumbertec3d(i1, i2, i3, nx1, nx2, nx3, ig, igrid)
Definition: convert.t:639
double precision length_convert_factor
double precision, dimension(ndim) dxlevel
integer snapshotini
Resume from the snapshot with this index.
type(state), pointer block
Block pointer for using one block and its previous state.
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine generate_plotfile
Definition: convert.t:3
integer nwauxio
Number of auxiliary variables that are only included in output.
integer icomm
The MPI communicator.
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
logical phys_energy
Solve energy equation or not.
Definition: mod_physics.t:25
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_convert), pointer phys_to_primitive
Definition: mod_physics.t:63
integer, dimension(:,:), allocatable node
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer type_block_xc_io
MPI type for IO: cell corner (xc) or cell center (xcc) coordinates.
logical spectrum
output spectrum
logical autoconvert
If true, already convert to output format during the run.
integer function nodenumbertec2d(i1, i2, nx1, nx2, ig, igrid)
Definition: convert.t:630
procedure(special_convert), pointer usr_special_convert
Handles computations for coordinates and variables in output.