MPI-AMRVAC  3.0
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
convert.t
Go to the documentation of this file.
1 !=============================================================================
2 subroutine generate_plotfile
7 use mod_convert, only: convert_all
9 !-----------------------------------------------------------------------------
10 
11 if(mype==0.and.level_io>0) write(unitterm,*)'reset tree to fixed level=',level_io
12 if(level_io>0 .or. level_io_min.ne.1 .or. level_io_max.ne.nlevelshi) then
14 else if(.not. phys_req_diagonal) then
15  call getbc(global_time,0.d0,ps,iwstart,nwgc)
16 end if
17 
18 select case(convert_type)
19  case('tecplot','tecplotCC','tecline')
20  call tecplot(unitconvert)
21  case('tecplotmpi','tecplotCCmpi','teclinempi')
23  case('vtu','vtuCC')
25  case('vtumpi','vtuCCmpi')
27  case('vtuB','vtuBCC','vtuBmpi','vtuBCCmpi')
29  case('vtuB64','vtuBCC64','vtuBmpi64','vtuBCCmpi64')
31  {^iftwod
32  case('vtuB23','vtuBCC23')
34  case('vtuBsym23','vtuBCCsym23')
36  }
37  case('pvtumpi','pvtuCCmpi')
39  case('pvtuBmpi','pvtuBCCmpi')
41  case('vtimpi','vtiCCmpi')
43  case('onegrid','onegridmpi')
44  call onegrid(unitconvert)
45  case('oneblock','oneblockB')
46  call oneblock(unitconvert)
47  case('EIvtiCCmpi','ESvtiCCmpi','SIvtiCCmpi','EIvtuCCmpi','ESvtuCCmpi','SIvtuCCmpi')
48  ! output synthetic euv emission
49  if (ndim==3 .and. slab .and. associated(phys_te_images)) then
50  call phys_te_images
51  endif
52  case('dat_generic_mpi')
53  ! it is the responsability of the physcics module to pouplate the array of
54  ! conversion methods by calling
55  call convert_all()
56  case('user','usermpi')
57  if (.not. associated(usr_special_convert)) then
58  call mpistop("usr_special_convert not defined")
59  else
61  end if
62  case default
63  call mpistop("Error in generate_plotfile: Unknown convert_type")
64 end select
65 
66 
67 
68 end subroutine generate_plotfile
69 !=============================================================================
70 subroutine oneblock(qunit)
71 ! this is for turning an AMR run into a single block
72 ! the data will be all on selected level level_io
73 
74 ! this version should work for any dimension
75 ! only writes w_write selected 1:nw variables, also nwauxio
76 ! may use saveprim to switch to primitives
77 ! this version can not work on multiple CPUs
78 ! does not renormalize variables
79 
80 ! header info differs from onegrid below
81 
82 ! ASCII or binary output
83 
87 use mod_physics
89 integer, intent(in) :: qunit
90 
91 integer :: Morton_no,igrid,ix^D,ig^D,level
92 integer, pointer :: ig_to_igrid(:^D&,:)
93 logical :: fileopen,writeblk(max_blocks)
94 character(len=80) :: filename
95 integer :: filenr,ncells^D,ncellx^D,jg^D,jig^D
96 
97 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
98 character(len=1024) :: outfilehead
99 
100 double precision :: wval1,xval1
101 double precision, dimension({^D&1:1},1:nw+nwauxio) :: wval
102 double precision, dimension({^D&1:1},1:ndim) :: xval
103 double precision:: normconv(0:nw+nwauxio)
104 
105 integer :: iw,iiw,writenw,iwrite(1:nw+nwauxio),iigrid,idim
106 logical :: patchw(ixG^T)
107 !-----------------------------------------------------------------------------
108 
109 if(level_io<1)then
110  call mpistop('please specify level_io>0 for usage with oneblock')
111 end if
112 
113 if(autoconvert)then
114  call mpistop('Set autoconvert=F and convert oneblock data manually')
115 end if
116 
117 if(npe>1)then
118  if(mype==0) print *,'ONEBLOCK as yet to be parallelized'
119  call mpistop('npe>1, oneblock')
120 end if
121 
122 ! only variables selected by w_write will be written out
123 normconv(0:nw+nwauxio)=one
124 normconv(0) = length_convert_factor
125 normconv(1:nw) = w_convert_factor
126 writenw=count(w_write(1:nw))+nwauxio
127 iiw=0
128 do iw =1,nw
129  if (.not.w_write(iw))cycle
130  iiw=iiw+1
131  iwrite(iiw)=iw
132 end do
133 if(nwauxio>0)then
134  do iw =nw+1,nw+nwauxio
135  iiw=iiw+1
136  iwrite(iiw)=iw
137  end do
138 endif
139 
140 allocate(ig_to_igrid(ng^d(level_io),0:npe-1))
141 ig_to_igrid=-1
142 writeblk=.false.
143 do morton_no=morton_start(mype),morton_stop(mype)
144  igrid=sfc_to_igrid(morton_no)
145  level=node(plevel_,igrid)
146  ig^d=igrid_to_node(igrid,mype)%node%ig^d;
147  ig_to_igrid(ig^d,mype)=igrid
148  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
149  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
150  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
151  writeblk(igrid)=.true.
152  end if
153 end do
154 
155 call getheadernames(wnamei,xandwnamei,outfilehead)
156 ncells^d=0;
157 ncellx^d=ixmhi^d-ixmlo^d+1\
158 {do ig^d=1,ng^d(level_io)\}
159  igrid=ig_to_igrid(ig^d,mype)
160  if(writeblk(igrid)) go to 20
161 {end do\}
162 20 continue
163 jg^d=ig^d;
164 {
165 jig^dd=jg^dd;
166 do ig^d=1,ng^d(level_io)
167  jig^d=ig^d
168  igrid=ig_to_igrid(jig^dd,mype)
169  if(writeblk(igrid)) ncells^d=ncells^d+ncellx^d
170 end do
171 \}
172 
173 do iigrid=1,igridstail; igrid=igrids(iigrid)
174  if(.not.writeblk(igrid)) cycle
175  block=>ps(igrid)
176  if (nwauxio > 0) then
177  if (.not. associated(usr_aux_output)) then
178  call mpistop("usr_aux_output not defined")
179  else
180  call usr_aux_output(ixg^ll,ixm^ll^ladd1, &
181  ps(igrid)%w,ps(igrid)%x,normconv)
182  end if
183  end if
184 end do
185 
186 if (saveprim) then
187  do iigrid=1,igridstail; igrid=igrids(iigrid)
188  if (.not.writeblk(igrid)) cycle
189  block=>ps(igrid)
190  call phys_to_primitive(ixg^ll,ixg^ll^lsub1,ps(igrid)%w,ps(igrid)%x)
191  if(b0field) then
192  ! add background magnetic field B0 to B
193  ps(igrid)%w(ixg^t,iw_mag(:))=ps(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
194  end if
195  end do
196 else
197  do iigrid=1,igridstail; igrid=igrids(iigrid)
198  if (.not.writeblk(igrid)) cycle
199  block=>ps(igrid)
200  if (b0field) then
201  ! add background magnetic field B0 to B
202  if(phys_energy) &
203  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) &
204  + sum(ps(igrid)%w(ixg^t,iw_mag(:))*ps(igrid)%B0(ixg^t,:,0),dim=ndim+1)
205  ps(igrid)%w(ixg^t,iw_mag(:))=ps(igrid)%w(ixg^t,iw_mag(:))+ps(igrid)%B0(ixg^t,:,0)
206  end if
207  end do
208 end if
209 
210 master_cpu_open : if (mype == 0) then
211  inquire(qunit,opened=fileopen)
212  if (.not.fileopen) then
213  ! generate filename
214  filenr=snapshotini
215  if (autoconvert) filenr=snapshotnext
216  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
217  select case(convert_type)
218  case("oneblock")
219  open(qunit,file=filename,status='unknown')
220  write(qunit,*) trim(outfilehead)
221  write(qunit,*) ncells^d
222  write(qunit,*) real(global_time*time_convert_factor)
223  case("oneblockB")
224  open(qunit,file=filename,form='unformatted',status='unknown')
225  write(qunit) outfilehead
226  write(qunit) ncells^d
227  write(qunit) real(global_time*time_convert_factor)
228  end select
229  end if
230 end if master_cpu_open
231 
232 {^ifthreed
233 do ig3=1,ng3(level_io)
234  do ix3=ixmlo3,ixmhi3}
235 
236  {^nooned
237  do ig2=1,ng2(level_io)
238  do ix2=ixmlo2,ixmhi2}
239 
240  do ig1=1,ng1(level_io)
241  igrid=ig_to_igrid(ig^d,mype)
242  if(.not.writeblk(igrid)) cycle
243  do ix1=ixmlo1,ixmhi1
244  master_write : if(mype==0) then
245  select case(convert_type)
246  case("oneblock")
247  write(qunit,fmt="(100(e14.6))") &
248  ps(igrid)%x(ix^d,1:ndim)*normconv(0),&
249  (ps(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw)),iw=1,writenw)
250  case("oneblockB")
251  write(qunit) real(ps(igrid)%x(ix^d,1:ndim)*normconv(0)),&
252  (real(ps(igrid)%w(ix^d,iwrite(iw))*normconv(iwrite(iw))),iw=1,writenw)
253  end select
254  end if master_write
255  end do
256  end do
257  {^nooned
258  end do
259  end do}
260  {^ifthreed
261  end do
262 end do}
263 
264 close(qunit)
265 
266 end subroutine oneblock
267 !=============================================================================
268 subroutine onegrid(qunit)
269 
270 ! this is for turning an AMR run into a single grid
271 ! this version should work for any dimension, can be in parallel
272 ! in 1D, should behave much like oneblock, except for header info
273 
274 ! only writes all 1:nw variables, no nwauxio
275 ! may use saveprim to switch to primitives
276 ! this version can work on multiple CPUs
277 ! does not renormalize variables
278 ! ASCII output
279 
282 use mod_physics
284 integer, intent(in) :: qunit
285 
286 integer :: itag,Morton_no,igrid,ix^D,iw
287 logical :: fileopen
288 character(len=80) :: filename
289 integer :: filenr
290 
291 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
292 character(len=1024) :: outfilehead
293 
294 !.. MPI variables ..
295 integer :: igrid_recv,ipe
296 double precision :: w_recv(ixG^T,1:nw),x_recv(ixG^T,1:ndim)
297 integer, allocatable :: intstatus(:,:)
298 
299 !-----------------------------------------------------------------------------
300 
301 if(nwauxio>0)then
302  if(mype==0) print *,'ONEGRID to be used without nwauxio'
303  call mpistop('nwauxio>0, onegrid')
304 end if
305 
306 if(saveprim)then
307  if(mype==0.and.nwaux>0) print *,'warning: ONEGRID used with saveprim, check auxiliaries'
308 end if
309 
310 
311 
312 master_cpu_open : if (mype == 0) then
313  call getheadernames(wnamei,xandwnamei,outfilehead)
314  write(outfilehead,'(a)') "#"//" "//trim(outfilehead)
315  inquire(qunit,opened=fileopen)
316  if (.not.fileopen) then
317  ! generate filename
318  filenr=snapshotini
319  if (autoconvert) filenr=snapshotnext
320  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".blk"
321  open(qunit,file=filename,status='unknown')
322  end if
323  write(qunit,"(a)")outfilehead
324  write(qunit,"(i7)") ( {^d&(ixmhi^d-ixmlo^d+1)*} )*(morton_stop(npe-1)-morton_start(0)+1)
325 end if master_cpu_open
326 
327 do morton_no=morton_start(mype),morton_stop(mype)
328  igrid=sfc_to_igrid(morton_no)
329  if(saveprim) call phys_to_primitive(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x)
330  if (mype/=0)then
331  itag=morton_no
332  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
333  call mpi_send(ps(igrid)%x,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
334  itag=igrid
335  call mpi_send(ps(igrid)%w,1,type_block_io, 0,itag,icomm,ierrmpi)
336  else
337  {do ix^db=ixmlo^db,ixmhi^db\}
338  do iw=1,nw
339  if( dabs(ps(igrid)%w(ix^d,iw)) < 1.0d-32 ) ps(igrid)%w(ix^d,iw) = zero
340  enddo
341  write(qunit,fmt="(100(e14.6))") ps(igrid)%x(ix^d,1:ndim)&
342  ,ps(igrid)%w(ix^d,1:nw)
343  {end do\}
344  end if
345 end do
346 
347 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
348 
349 manycpu : if (npe>1) then
350  if (mype==0) then
351  loop_cpu : do ipe =1, npe-1
352  loop_morton : do morton_no=morton_start(ipe),morton_stop(ipe)
353  itag=morton_no
354  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
355  call mpi_recv(x_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
356  itag=igrid_recv
357  call mpi_recv(w_recv,1,type_block_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
358  {do ix^db=ixmlo^db,ixmhi^db\}
359  do iw=1,nw
360  if( dabs(ps(igrid)%w(ix^d,iw)) < smalldouble ) ps(igrid)%w(ix^d,iw) = zero
361  enddo
362  write(qunit,fmt="(100(e14.6))") x_recv(ix^d,1:ndim)&
363  ,w_recv(ix^d,1:nw)
364  {end do\}
365  end do loop_morton
366  end do loop_cpu
367  end if
368 end if manycpu
369 
370 if (npe>1) then
371  call mpi_barrier(icomm,ierrmpi)
372  if(mype==0)deallocate(intstatus)
373 endif
374 
375 if(mype==0) close(qunit)
376 end subroutine onegrid
377 !============================================================================
378 subroutine tecplot(qunit)
379 
380 ! output for tecplot (ASCII format)
381 ! not parallel, uses calc_grid to compute nwauxio variables
382 ! allows renormalizing using convert factors
383 
386 
387 integer, intent(in) :: qunit
388 
389 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
390 integer:: NumGridsOnLevel(1:nlevelshi)
391 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
392 
393 integer :: nodes, elems
394 
395 
396 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
397 
398 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
399 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
400 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
401 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
402 
403 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
404 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
405 double precision, dimension(0:nw+nwauxio) :: normconv
406 logical :: fileopen,first
407 character(len=80) :: filename
408 integer :: filenr
409 
410 !!! possible length conflict
411 character(len=1024) :: tecplothead
412 
413 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
414 character(len=1024) :: outfilehead
415 !-----------------------------------------------------------------------------
416 
417 if(npe>1)then
418  if(mype==0) print *,'tecplot not parallel, use tecplotmpi'
419  call mpistop('npe>1, tecplot')
420 end if
421 
422 if(nw/=count(w_write(1:nw)))then
423  if(mype==0) print *,'tecplot does not use w_write=F'
424  call mpistop('w_write, tecplot')
425 end if
426 
427 if(nocartesian)then
428  if(mype==0) print *,'tecplot with nocartesian'
429 endif
430 
431 inquire(qunit,opened=fileopen)
432 if (.not.fileopen) then
433  ! generate filename
434  filenr=snapshotini
435  if (autoconvert) filenr=snapshotnext
436  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
437  open(qunit,file=filename,status='unknown')
438 end if
439 
440 call getheadernames(wnamei,xandwnamei,outfilehead)
441 
442 write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
443 
444 write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
445 
446 numgridsonlevel(1:nlevelshi)=0
447 do level=levmin,levmax
448  numgridsonlevel(level)=0
449  do iigrid=1,igridstail; igrid=igrids(iigrid);
450  if (node(plevel_,igrid)/=level) cycle
451  numgridsonlevel(level)=numgridsonlevel(level)+1
452  end do
453 end do
454 
455 nx^d=ixmhi^d-ixmlo^d+1;
456 nxc^d=nx^d+1;
457 
458 {^ifoned
459 if(convert_type=='tecline') then
460  nodes=0
461  elems=0
462  do level=levmin,levmax
463  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
464  elems=elems + numgridsonlevel(level)*{nx^d*}
465  enddo
466 
467  write(qunit,"(a,i7,a,1pe12.5,a)") &
468  'ZONE T="all levels", I=',elems, &
469  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
470 
471  igonlevel=0
472  do iigrid=1,igridstail; igrid=igrids(iigrid);
473  block=>ps(igrid)
474  call calc_x(igrid,xc,xcc)
475  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
476  {do ix^db=ixccmin^db,ixccmax^db\}
477  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
478  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
479  !write(qunit,fmt="(100(e14.6))") x_TEC, w_TEC
480  write(qunit,fmt="(100(e24.16))") x_tec, w_tec
481  {end do\}
482  enddo
483  close(qunit)
484 else
485 }
486 do level=levmin,levmax
487  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
488  elemsonlevel=numgridsonlevel(level)*{nx^d*}
489  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
490  ! with the AMR grid LEVEL. Other options would be
491  ! let each grid define a zone: inefficient for TECPLOT internal workings
492  ! hence not implemented
493  ! let entire octree define 1 zone: no difference in interpolation
494  ! properties across TECPLOT zones detected as yet, hence not done
495  select case(convert_type)
496  case('tecplot')
497  ! in this option, we store the corner coordinates, as well as the corner
498  ! values of all variables (obtained by averaging). This allows POINT packaging,
499  ! and thus we can save full grid info by using one call to calc_grid
500  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
501  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
502  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
503  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
504  do iigrid=1,igridstail; igrid=igrids(iigrid);
505  if (node(plevel_,igrid)/=level) cycle
506  block=>ps(igrid)
507  call calc_x(igrid,xc,xcc)
508  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
509  ixc^l,ixcc^l,.true.)
510  {do ix^db=ixcmin^db,ixcmax^db\}
511  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
512  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
513  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
514  {end do\}
515  enddo
516  case('tecplotCC')
517  ! in this option, we store the corner coordinates, and the cell center
518  ! values of all variables. Due to this mix of corner/cell center, we must
519  ! use BLOCK packaging, and thus we have enormous overhead by using
520  ! calc_grid repeatedly to merely fill values of cell corner coordinates
521  ! and cell center values per dimension, per variable
522  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
523  if(nw+nwauxio==1)then
524  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
525  ! and just set [ndim+1]
526  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
527  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
528  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
529  ndim+1,']=CELLCENTERED), ZONETYPE=', &
530  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
531  else
532  if(ndim+nw+nwauxio<10) then
533  ! difference only in length of integer format specification for ndim+nw+nwauxio
534  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
535  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
536  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
537  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
538  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
539  else
540  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
541  'ZONE T="',level,'"',', N=',nodesonlevel,', E=',elemsonlevel, &
542  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
543  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
544  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
545  endif
546  endif
547  do idim=1,ndim
548  first=(idim==1)
549  do iigrid=1,igridstail; igrid=igrids(iigrid);
550  if (node(plevel_,igrid)/=level) cycle
551  block=>ps(igrid)
552  call calc_x(igrid,xc,xcc)
553  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
554  ixc^l,ixcc^l,first)
555  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
556  enddo
557  enddo
558  do iw=1,nw+nwauxio
559  do iigrid=1,igridstail; igrid=igrids(iigrid);
560  if (node(plevel_,igrid)/=level) cycle
561  block=>ps(igrid)
562  call calc_x(igrid,xc,xcc)
563  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
564  ixc^l,ixcc^l,.true.)
565  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
566  enddo
567  enddo
568  case default
569  call mpistop('no such tecplot type')
570  end select
571  igonlevel=0
572  do iigrid=1,igridstail; igrid=igrids(iigrid);
573  if (node(plevel_,igrid)/=level) cycle
574  block=>ps(igrid)
575  igonlevel=igonlevel+1
576  call save_conntec(qunit,igrid,igonlevel)
577  end do
578 end do
579 {^ifoned endif}
580 
581 close(qunit)
582 
583 end subroutine tecplot
584 !=============================================================================
585 subroutine save_conntec(qunit,igrid,igonlevel)
586 
587 ! this saves the basic line, quad and brick connectivity,
588 ! as used by TECPLOT file outputs for unstructured grid
589 
591 
592 integer, intent(in) :: qunit, igrid, igonlevel
593 
594 integer :: nx^D, nxC^D, ix^D
595 {^ifoned integer, external:: nodenumbertec1D \}
596 {^iftwod integer, external:: nodenumbertec2D \}
597 {^ifthreed integer, external:: nodenumbertec3D \}
598 !-----------------------------------------------------------------------------
599 nx^d=ixmhi^d-ixmlo^d+1;
600 nxc^d=nx^d+1;
601 
602 ! connectivity list
603 {do ix^db=1,nx^db\}
604  {^ifthreed
605  ! basic brick connectivity
606  write(qunit,'(8(i7,1x))') &
607  nodenumbertec3d(ix1, ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
608  nodenumbertec3d(ix1+1,ix2-1,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
609  nodenumbertec3d(ix1+1,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
610  nodenumbertec3d(ix1 ,ix2 ,ix3-1,nxc1,nxc2,nxc3,igonlevel,igrid),&
611  nodenumbertec3d(ix1 ,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
612  nodenumbertec3d(ix1+1,ix2-1,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
613  nodenumbertec3d(ix1+1,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid),&
614  nodenumbertec3d(ix1 ,ix2 ,ix3 ,nxc1,nxc2,nxc3,igonlevel,igrid)
615  }
616  {^iftwod
617  ! basic quadrilateral connectivity
618  write(qunit,'(4(i7,1x))') &
619  nodenumbertec2d(ix1, ix2-1,nxc1,nxc2,igonlevel,igrid),&
620  nodenumbertec2d(ix1+1,ix2-1,nxc1,nxc2,igonlevel,igrid),&
621  nodenumbertec2d(ix1+1,ix2 ,nxc1,nxc2,igonlevel,igrid),&
622  nodenumbertec2d(ix1 ,ix2 ,nxc1,nxc2,igonlevel,igrid)
623  }
624  {^ifoned
625  ! basic line connectivity
626  write(qunit,'(2(i7,1x))') nodenumbertec1d(ix1,nxc1,igonlevel,igrid),&
627  nodenumbertec1d(ix1+1,nxc1,igonlevel,igrid)
628  }
629 {end do\}
630 
631 end subroutine save_conntec
632 !=============================================================================
633 integer function nodenumbertec1d(i1,nx1,ig,igrid)
634 
635 integer, intent(in):: i1,nx1,ig,igrid
636 !-----------------------------------------------------------------------------
637 nodenumbertec1d=i1+(ig-1)*nx1
638 
639 if(nodenumbertec1d>9999999)call mpistop("too large nodenumber")
640 end function nodenumbertec1d
641 !====================================================================================
642 integer function nodenumbertec2d(i1,i2,nx1,nx2,ig,igrid)
643 
644 integer, intent(in):: i1,i2,nx1,nx2,ig,igrid
645 !-----------------------------------------------------------------------------
646 nodenumbertec2d=i1+i2*nx1+(ig-1)*nx1*nx2
647 
648 if(nodenumbertec2d>9999999)call mpistop("too large nodenumber")
649 end function nodenumbertec2d
650 !====================================================================================
651 integer function nodenumbertec3d(i1,i2,i3,nx1,nx2,nx3,ig,igrid)
652 
653 integer, intent(in):: i1,i2,i3,nx1,nx2,nx3,ig,igrid
654 !-----------------------------------------------------------------------------
655 nodenumbertec3d=i1+i2*nx1+i3*nx1*nx2+(ig-1)*nx1*nx2*nx3
656 
657 if(nodenumbertec3d>9999999)call mpistop("too large nodenumber")
658 end function nodenumbertec3d
659 !=============================================================================
660 subroutine unstructuredvtk(qunit)
661 
662 ! output for vtu format to paraview
663 ! not parallel, uses calc_grid to compute nwauxio variables
664 ! allows renormalizing using convert factors
665 
668 
669 integer, intent(in) :: qunit
670 
671 double precision :: x_VTK(1:3)
672 
673 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
674 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
675 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
676 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
677 
678 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
679 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
680 double precision, dimension(0:nw+nwauxio) :: normconv
681 integer:: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,iw
682 integer:: NumGridsOnLevel(1:nlevelshi)
683 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
684 
685 character(len=80):: filename
686 integer :: filenr
687 
688 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
689 character(len=1024) :: outfilehead
690 
691 logical :: fileopen
692 !-----------------------------------------------------------------------------
693 
694 if(npe>1)then
695  if(mype==0) print *,'unstructuredvtk not parallel, use vtumpi'
696  call mpistop('npe>1, unstructuredvtk')
697 end if
698 
699 inquire(qunit,opened=fileopen)
700 if(.not.fileopen)then
701  ! generate filename
702  filenr=snapshotini
703  if (autoconvert) filenr=snapshotnext
704  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
705  ! Open the file for the header part
706  open(qunit,file=filename,status='unknown')
707 end if
708 
709 call getheadernames(wnamei,xandwnamei,outfilehead)
710 
711 ! generate xml header
712 write(qunit,'(a)')'<?xml version="1.0"?>'
713 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
714 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
715 write(qunit,'(a)')'<UnstructuredGrid>'
716 write(qunit,'(a)')'<FieldData>'
717 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
718  'NumberOfTuples="1" format="ascii">'
719 write(qunit,*) dble(dble(global_time)*time_convert_factor)
720 write(qunit,'(a)')'</DataArray>'
721 write(qunit,'(a)')'</FieldData>'
722 
723 ! number of cells, number of corner points, per grid.
724 nx^d=ixmhi^d-ixmlo^d+1;
725 nxc^d=nx^d+1;
726 nc={nx^d*}
727 np={nxc^d*}
728 
729 ! Note: using the w_write, writelevel, writespshift
730 ! we can clip parts of the grid away, select variables, levels etc.
731 do level=levmin,levmax
732  if (writelevel(level)) then
733  do iigrid=1,igridstail; igrid=igrids(iigrid);
734  if (node(plevel_,igrid)/=level) cycle
735  block=>ps(igrid)
736  ! only output a grid when fully within clipped region selected
737  ! by writespshift array
738  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
739  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
740  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
741  call calc_x(igrid,xc,xcc)
742  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
743  ixc^l,ixcc^l,.true.)
744  select case(convert_type)
745  case('vtu')
746  ! we write out every grid as one VTK PIECE
747  write(qunit,'(a,i7,a,i7,a)') &
748  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
749  write(qunit,'(a)')'<PointData>'
750  do iw=1,nw+nwauxio
751  if(iw<=nw) then
752  if(.not.w_write(iw)) cycle
753  end if
754  write(qunit,'(a,a,a)')&
755  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
756  write(qunit,'(200(1pe14.6))') {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
757  write(qunit,'(a)')'</DataArray>'
758  end do
759  write(qunit,'(a)')'</PointData>'
760  write(qunit,'(a)')'<Points>'
761  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
762  ! write cell corner coordinates in a backward dimensional loop, always 3D output
763  {do ix^db=ixcmin^db,ixcmax^db \}
764  x_vtk(1:3)=zero;
765  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
766  write(qunit,'(3(1pe14.6))') x_vtk
767  {end do \}
768  write(qunit,'(a)')'</DataArray>'
769  write(qunit,'(a)')'</Points>'
770  case('vtuCC')
771  ! we write out every grid as one VTK PIECE
772  write(qunit,'(a,i7,a,i7,a)') &
773  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
774  write(qunit,'(a)')'<CellData>'
775  do iw=1,nw+nwauxio
776  if(iw<=nw) then
777  if(.not.w_write(iw)) cycle
778  end if
779  write(qunit,'(a,a,a)')&
780  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
781  write(qunit,'(200(1pe14.6))') {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
782  write(qunit,'(a)')'</DataArray>'
783  end do
784  write(qunit,'(a)')'</CellData>'
785  write(qunit,'(a)')'<Points>'
786  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
787  ! write cell corner coordinates in a backward dimensional loop, always 3D output
788  {do ix^db=ixcmin^db,ixcmax^db \}
789  x_vtk(1:3)=zero;
790  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
791  write(qunit,'(3(1pe14.6))') x_vtk
792  {end do \}
793  write(qunit,'(a)')'</DataArray>'
794  write(qunit,'(a)')'</Points>'
795  end select
796 
797  write(qunit,'(a)')'<Cells>'
798  ! connectivity part
799  write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
800  call save_connvtk(qunit,igrid)
801  write(qunit,'(a)')'</DataArray>'
802 
803  ! offsets data array
804  write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
805  do icel=1,nc
806  write(qunit,'(i7)') icel*(2**^nd)
807  end do
808  write(qunit,'(a)')'</DataArray>'
809 
810  ! VTK cell type data array
811  write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
812  ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
813  {^ifoned vtk_type=3 \}
814  {^iftwod vtk_type=8 \}
815  {^ifthreed vtk_type=11 \}
816  do icel=1,nc
817  write(qunit,'(i2)') vtk_type
818  end do
819  write(qunit,'(a)')'</DataArray>'
820 
821  write(qunit,'(a)')'</Cells>'
822 
823  write(qunit,'(a)')'</Piece>'
824  end if
825  end do
826  end if
827 end do
828 
829 write(qunit,'(a)')'</UnstructuredGrid>'
830 write(qunit,'(a)')'</VTKFile>'
831 close(qunit)
832 
833 end subroutine unstructuredvtk
834 !====================================================================================
835 subroutine unstructuredvtkb(qunit)
836 
837 ! output for vtu format to paraview, binary version output
838 ! not parallel, uses calc_grid to compute nwauxio variables
839 ! allows renormalizing using convert factors
843 
844 integer, intent(in) :: qunit
845 
846 double precision :: x_VTK(1:3)
847 
848 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
849 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
850 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
851 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
852 
853 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
854 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
855 double precision :: normconv(0:nw+nwauxio)
856 
857 integer, allocatable :: intstatus(:,:)
858 integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
859 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
860 integer*8 :: offset
861 
862 integer:: k,iw
863 integer:: length,lengthcc,length_coords,length_conn,length_offsets
864 character:: buf
865 character(len=80):: filename
866 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
867 character(len=1024) :: outfilehead
868 
869 logical :: fileopen,cell_corner=.false.
870 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
871 !-----------------------------------------------------------------------------
872 
873 normconv=one
874 morton_length=morton_stop(npe-1)-morton_start(0)+1
875 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
876 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
877 morton_aim=.false.
878 morton_aim_p=.false.
879 do morton_no=morton_start(mype),morton_stop(mype)
880  igrid=sfc_to_igrid(morton_no)
881  level=node(plevel_,igrid)
882  ! we can clip parts of the grid away, select variables, levels etc.
883  if(writelevel(level)) then
884  ! only output a grid when fully within clipped region selected
885  ! by writespshift array
886  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
887  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
888  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
889  morton_aim_p(morton_no)=.true.
890  end if
891  end if
892 end do
893 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
894  icomm,ierrmpi)
895 select case(convert_type)
896  case('vtuB','vtuBmpi')
897  cell_corner=.true.
898  case('vtuBCC','vtuBCCmpi')
899  cell_corner=.false.
900 end select
901 if (mype /= 0) then
902  do morton_no=morton_start(mype),morton_stop(mype)
903  if(.not. morton_aim(morton_no)) cycle
904  igrid=sfc_to_igrid(morton_no)
905  call calc_x(igrid,xc,xcc)
906  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
907  ixc^l,ixcc^l,.true.)
908  itag=morton_no
909  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
910  if(cell_corner) then
911  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
912  else
913  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
914  end if
915  end do
916 
917 else
918  ! mype==0
919  offset=0
920  inquire(qunit,opened=fileopen)
921  if(.not.fileopen)then
922  ! generate filename
923  filenr=snapshotini
924  if (autoconvert) filenr=snapshotnext
925  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
926  ! Open the file for the header part
927  open(qunit,file=filename,status='replace')
928  end if
929  call getheadernames(wnamei,xandwnamei,outfilehead)
930  ! generate xml header
931  write(qunit,'(a)')'<?xml version="1.0"?>'
932  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
933  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
934  write(qunit,'(a)')'<UnstructuredGrid>'
935  write(qunit,'(a)')'<FieldData>'
936  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
937  'NumberOfTuples="1" format="ascii">'
938  write(qunit,*) real(global_time*time_convert_factor)
939  write(qunit,'(a)')'</DataArray>'
940  write(qunit,'(a)')'</FieldData>'
941 
942  ! number of cells, number of corner points, per grid.
943  nx^d=ixmhi^d-ixmlo^d+1;
944  nxc^d=nx^d+1;
945  nc={nx^d*}
946  np={nxc^d*}
947  length=np*size_real
948  lengthcc=nc*size_real
949  length_coords=3*length
950  length_conn=2**^nd*size_int*nc
951  length_offsets=nc*size_int
952 
953  ! Note: using the w_write, writelevel, writespshift
954  do morton_no=morton_start(0),morton_stop(0)
955  if(.not. morton_aim(morton_no)) cycle
956  if(cell_corner) then
957  ! we write out every grid as one VTK PIECE
958  write(qunit,'(a,i7,a,i7,a)') &
959  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
960  write(qunit,'(a)')'<PointData>'
961  do iw=1,nw+nwauxio
962  if(iw<=nw) then
963  if(.not.w_write(iw)) cycle
964  endif
965 
966  write(qunit,'(a,a,a,i16,a)')&
967  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
968  '" format="appended" offset="',offset,'">'
969  write(qunit,'(a)')'</DataArray>'
970  offset=offset+length+size_int
971  end do
972  write(qunit,'(a)')'</PointData>'
973  write(qunit,'(a)')'<Points>'
974  write(qunit,'(a,i16,a)') &
975  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
976  ! write cell corner coordinates in a backward dimensional loop, always 3D output
977  offset=offset+length_coords+size_int
978  write(qunit,'(a)')'</Points>'
979  else
980  ! we write out every grid as one VTK PIECE
981  write(qunit,'(a,i7,a,i7,a)') &
982  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
983  write(qunit,'(a)')'<CellData>'
984  do iw=1,nw+nwauxio
985  if(iw<=nw) then
986  if(.not.w_write(iw)) cycle
987  end if
988  write(qunit,'(a,a,a,i16,a)')&
989  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
990  '" format="appended" offset="',offset,'">'
991  write(qunit,'(a)')'</DataArray>'
992  offset=offset+lengthcc+size_int
993  end do
994  write(qunit,'(a)')'</CellData>'
995  write(qunit,'(a)')'<Points>'
996  write(qunit,'(a,i16,a)') &
997  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
998  ! write cell corner coordinates in a backward dimensional loop, always 3D output
999  offset=offset+length_coords+size_int
1000  write(qunit,'(a)')'</Points>'
1001  end if
1002  write(qunit,'(a)')'<Cells>'
1003 
1004  ! connectivity part
1005  write(qunit,'(a,i16,a)')&
1006  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1007  offset=offset+length_conn+size_int
1008 
1009  ! offsets data array
1010  write(qunit,'(a,i16,a)') &
1011  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1012  offset=offset+length_offsets+size_int
1013 
1014  ! VTK cell type data array
1015  write(qunit,'(a,i16,a)') &
1016  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1017  offset=offset+size_int+nc*size_int
1018 
1019  write(qunit,'(a)')'</Cells>'
1020 
1021  write(qunit,'(a)')'</Piece>'
1022  end do
1023  ! write metadata communicated from other processors
1024  if(npe>1)then
1025  do ipe=1, npe-1
1026  do morton_no=morton_start(ipe),morton_stop(ipe)
1027  if(.not. morton_aim(morton_no)) cycle
1028  if(cell_corner) then
1029  ! we write out every grid as one VTK PIECE
1030  write(qunit,'(a,i7,a,i7,a)') &
1031  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1032  write(qunit,'(a)')'<PointData>'
1033  do iw=1,nw+nwauxio
1034  if(iw<=nw) then
1035  if(.not.w_write(iw)) cycle
1036  end if
1037  write(qunit,'(a,a,a,i16,a)')&
1038  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1039  '" format="appended" offset="',offset,'">'
1040  write(qunit,'(a)')'</DataArray>'
1041  offset=offset+length+size_int
1042  end do
1043  write(qunit,'(a)')'</PointData>'
1044  write(qunit,'(a)')'<Points>'
1045  write(qunit,'(a,i16,a)') &
1046  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1047  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1048  offset=offset+length_coords+size_int
1049  write(qunit,'(a)')'</Points>'
1050  else
1051  ! we write out every grid as one VTK PIECE
1052  write(qunit,'(a,i7,a,i7,a)') &
1053  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1054  write(qunit,'(a)')'<CellData>'
1055  do iw=1,nw+nwauxio
1056  if(iw<=nw) then
1057  if(.not.w_write(iw)) cycle
1058  end if
1059  write(qunit,'(a,a,a,i16,a)')&
1060  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
1061  '" format="appended" offset="',offset,'">'
1062  write(qunit,'(a)')'</DataArray>'
1063  offset=offset+lengthcc+size_int
1064  end do
1065  write(qunit,'(a)')'</CellData>'
1066  write(qunit,'(a)')'<Points>'
1067  write(qunit,'(a,i16,a)') &
1068  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1069  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1070  offset=offset+length_coords+size_int
1071  write(qunit,'(a)')'</Points>'
1072  end if
1073  write(qunit,'(a)')'<Cells>'
1074  ! connectivity part
1075  write(qunit,'(a,i16,a)')&
1076  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1077  offset=offset+length_conn+size_int
1078 
1079  ! offsets data array
1080  write(qunit,'(a,i16,a)') &
1081  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1082  offset=offset+length_offsets+size_int
1083 
1084  ! VTK cell type data array
1085  write(qunit,'(a,i16,a)') &
1086  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1087  offset=offset+size_int+nc*size_int
1088 
1089  write(qunit,'(a)')'</Cells>'
1090 
1091  write(qunit,'(a)')'</Piece>'
1092  end do
1093  end do
1094  end if
1095 
1096  write(qunit,'(a)')'</UnstructuredGrid>'
1097  write(qunit,'(a)')'<AppendedData encoding="raw">'
1098  close(qunit)
1099  open(qunit,file=filename,access='stream',form='unformatted',position='append')
1100  buf='_'
1101  write(qunit) trim(buf)
1102 
1103  do morton_no=morton_start(0),morton_stop(0)
1104  if(.not. morton_aim(morton_no)) cycle
1105  igrid=sfc_to_igrid(morton_no)
1106  call calc_x(igrid,xc,xcc)
1107  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1108  ixc^l,ixcc^l,.true.)
1109  do iw=1,nw+nwauxio
1110  if(iw<=nw) then
1111  if(.not.w_write(iw)) cycle
1112  end if
1113  if(cell_corner) then
1114  write(qunit) length
1115  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
1116  else
1117  write(qunit) lengthcc
1118  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
1119  end if
1120  end do
1121 
1122  write(qunit) length_coords
1123  {do ix^db=ixcmin^db,ixcmax^db \}
1124  x_vtk(1:3)=zero;
1125  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1126  do k=1,3
1127  write(qunit) real(x_vtk(k))
1128  end do
1129  {end do \}
1130 
1131  write(qunit) length_conn
1132  {do ix^db=1,nx^db\}
1133  {^ifoned write(qunit)ix1-1,ix1 \}
1134  {^iftwod
1135  write(qunit)(ix2-1)*nxc1+ix1-1, &
1136  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1137  \}
1138  {^ifthreed
1139  write(qunit)&
1140  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1141  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1142  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1143  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1144  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1145  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1146  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1147  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1148  \}
1149  {end do\}
1150 
1151  write(qunit) length_offsets
1152  do icel=1,nc
1153  write(qunit) icel*(2**^nd)
1154  end do
1155 
1156  {^ifoned vtk_type=3 \}
1157  {^iftwod vtk_type=8 \}
1158  {^ifthreed vtk_type=11 \}
1159  write(qunit) size_int*nc
1160  do icel=1,nc
1161  write(qunit) vtk_type
1162  end do
1163  end do
1164  allocate(intstatus(mpi_status_size,1))
1165  if(npe>1)then
1166  ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1167  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1168  do ipe=1, npe-1
1169  do morton_no=morton_start(ipe),morton_stop(ipe)
1170  if(.not. morton_aim(morton_no)) cycle
1171  itag=morton_no
1172  call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1173  if(cell_corner) then
1174  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1175  else
1176  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1177  end if
1178  do iw=1,nw+nwauxio
1179  if(iw<=nw) then
1180  if(.not.w_write(iw)) cycle
1181  end if
1182  if(cell_corner) then
1183  write(qunit) length
1184  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
1185  else
1186  write(qunit) lengthcc
1187  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
1188  end if
1189  end do
1190 
1191  write(qunit) length_coords
1192  {do ix^db=ixcmin^db,ixcmax^db \}
1193  x_vtk(1:3)=zero;
1194  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1195  do k=1,3
1196  write(qunit) real(x_vtk(k))
1197  end do
1198  {end do \}
1199 
1200  write(qunit) length_conn
1201  {do ix^db=1,nx^db\}
1202  {^ifoned write(qunit)ix1-1,ix1 \}
1203  {^iftwod
1204  write(qunit)(ix2-1)*nxc1+ix1-1, &
1205  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1206  \}
1207  {^ifthreed
1208  write(qunit)&
1209  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1210  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1211  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1212  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1213  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1214  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1215  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1216  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1217  \}
1218  {end do\}
1219 
1220  write(qunit) length_offsets
1221  do icel=1,nc
1222  write(qunit) icel*(2**^nd)
1223  end do
1224  {^ifoned vtk_type=3 \}
1225  {^iftwod vtk_type=8 \}
1226  {^ifthreed vtk_type=11 \}
1227  write(qunit) size_int*nc
1228  do icel=1,nc
1229  write(qunit) vtk_type
1230  end do
1231  end do
1232  end do
1233  end if
1234  close(qunit)
1235  open(qunit,file=filename,status='unknown',form='formatted',position='append')
1236  write(qunit,'(a)')'</AppendedData>'
1237  write(qunit,'(a)')'</VTKFile>'
1238  close(qunit)
1239  deallocate(intstatus)
1240 end if
1241 
1242 deallocate(morton_aim,morton_aim_p)
1243 if (npe>1) then
1244  call mpi_barrier(icomm,ierrmpi)
1245 end if
1246 
1247 end subroutine unstructuredvtkb
1248 !====================================================================================
1249 subroutine unstructuredvtkb64(qunit)
1250 
1251 ! output for vtu format to paraview, binary version output
1252 ! not parallel, uses calc_grid to compute nwauxio variables
1253 ! allows renormalizing using convert factors
1256 use mod_calculate_xw
1257 
1258 integer, intent(in) :: qunit
1259 
1260 double precision :: x_VTK(1:3)
1261 
1262 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1263 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1264 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1265 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1266 
1267 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio):: wC_TMP
1268 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1269 double precision :: normconv(0:nw+nwauxio)
1270 
1271 integer, allocatable :: intstatus(:,:)
1272 integer :: itag,ipe,igrid,level,icel,ixC^L,ixCC^L,Morton_no,Morton_length
1273 integer :: nx^D,nxC^D,nc,np,VTK_type,ix^D,filenr
1274 integer*8 :: offset
1275 
1276 integer:: k,iw
1277 integer:: length,lengthcc,length_coords,length_conn,length_offsets
1278 character:: buf
1279 character(len=80):: filename
1280 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1281 character(len=1024) :: outfilehead
1282 
1283 logical :: fileopen,cell_corner=.false.
1284 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1285 !-----------------------------------------------------------------------------
1286 
1287 normconv=one
1288 morton_length=morton_stop(npe-1)-morton_start(0)+1
1289 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1290 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1291 morton_aim=.false.
1292 morton_aim_p=.false.
1293 do morton_no=morton_start(mype),morton_stop(mype)
1294  igrid=sfc_to_igrid(morton_no)
1295  level=node(plevel_,igrid)
1296  ! we can clip parts of the grid away, select variables, levels etc.
1297  if(writelevel(level)) then
1298  ! only output a grid when fully within clipped region selected
1299  ! by writespshift array
1300  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1301  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1302  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1303  morton_aim_p(morton_no)=.true.
1304  end if
1305  end if
1306 end do
1307 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1308  icomm,ierrmpi)
1309 select case(convert_type)
1310  case('vtuB64','vtuBmpi64')
1311  cell_corner=.true.
1312  case('vtuBCC64','vtuBCCmpi64')
1313  cell_corner=.false.
1314 end select
1315 if (mype /= 0) then
1316  do morton_no=morton_start(mype),morton_stop(mype)
1317  if(.not. morton_aim(morton_no)) cycle
1318  igrid=sfc_to_igrid(morton_no)
1319  call calc_x(igrid,xc,xcc)
1320  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1321  ixc^l,ixcc^l,.true.)
1322  itag=morton_no
1323  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
1324  if(cell_corner) then
1325  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1326  else
1327  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1328  end if
1329  end do
1330 
1331 else
1332  ! mype==0
1333  offset=0
1334  inquire(qunit,opened=fileopen)
1335  if(.not.fileopen)then
1336  ! generate filename
1337  filenr=snapshotini
1338  if (autoconvert) filenr=snapshotnext
1339  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
1340  ! Open the file for the header part
1341  open(qunit,file=filename,status='replace')
1342  end if
1343  call getheadernames(wnamei,xandwnamei,outfilehead)
1344  ! generate xml header
1345  write(qunit,'(a)')'<?xml version="1.0"?>'
1346  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1347  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1348  write(qunit,'(a)')'<UnstructuredGrid>'
1349  write(qunit,'(a)')'<FieldData>'
1350  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1351  'NumberOfTuples="1" format="ascii">'
1352  write(qunit,*) real(global_time*time_convert_factor)
1353  write(qunit,'(a)')'</DataArray>'
1354  write(qunit,'(a)')'</FieldData>'
1355 
1356  ! number of cells, number of corner points, per grid.
1357  nx^d=ixmhi^d-ixmlo^d+1;
1358  nxc^d=nx^d+1;
1359  nc={nx^d*}
1360  np={nxc^d*}
1361  length=np*size_double
1362  lengthcc=nc*size_double
1363  length_coords=3*length
1364  length_conn=2**^nd*size_int*nc
1365  length_offsets=nc*size_int
1366 
1367  ! Note: using the w_write, writelevel, writespshift
1368  do morton_no=morton_start(0),morton_stop(0)
1369  if(.not. morton_aim(morton_no)) cycle
1370  if(cell_corner) then
1371  ! we write out every grid as one VTK PIECE
1372  write(qunit,'(a,i7,a,i7,a)') &
1373  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1374  write(qunit,'(a)')'<PointData>'
1375  do iw=1,nw+nwauxio
1376  if(iw<=nw) then
1377  if(.not.w_write(iw)) cycle
1378  endif
1379 
1380  write(qunit,'(a,a,a,i16,a)')&
1381  '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1382  '" format="appended" offset="',offset,'">'
1383  write(qunit,'(a)')'</DataArray>'
1384  offset=offset+length+size_int
1385  end do
1386  write(qunit,'(a)')'</PointData>'
1387  write(qunit,'(a)')'<Points>'
1388  write(qunit,'(a,i16,a)') &
1389  '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1390  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1391  offset=offset+length_coords+size_int
1392  write(qunit,'(a)')'</Points>'
1393  else
1394  ! we write out every grid as one VTK PIECE
1395  write(qunit,'(a,i7,a,i7,a)') &
1396  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1397  write(qunit,'(a)')'<CellData>'
1398  do iw=1,nw+nwauxio
1399  if(iw<=nw) then
1400  if(.not.w_write(iw)) cycle
1401  end if
1402  write(qunit,'(a,a,a,i16,a)')&
1403  '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1404  '" format="appended" offset="',offset,'">'
1405  write(qunit,'(a)')'</DataArray>'
1406  offset=offset+lengthcc+size_int
1407  end do
1408  write(qunit,'(a)')'</CellData>'
1409  write(qunit,'(a)')'<Points>'
1410  write(qunit,'(a,i16,a)') &
1411  '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1412  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1413  offset=offset+length_coords+size_int
1414  write(qunit,'(a)')'</Points>'
1415  end if
1416  write(qunit,'(a)')'<Cells>'
1417 
1418  ! connectivity part
1419  write(qunit,'(a,i16,a)')&
1420  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1421  offset=offset+length_conn+size_int
1422 
1423  ! offsets data array
1424  write(qunit,'(a,i16,a)') &
1425  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1426  offset=offset+length_offsets+size_int
1427 
1428  ! VTK cell type data array
1429  write(qunit,'(a,i16,a)') &
1430  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1431  offset=offset+size_int+nc*size_int
1432 
1433  write(qunit,'(a)')'</Cells>'
1434 
1435  write(qunit,'(a)')'</Piece>'
1436  end do
1437  ! write metadata communicated from other processors
1438  if(npe>1)then
1439  do ipe=1, npe-1
1440  do morton_no=morton_start(ipe),morton_stop(ipe)
1441  if(.not. morton_aim(morton_no)) cycle
1442  if(cell_corner) then
1443  ! we write out every grid as one VTK PIECE
1444  write(qunit,'(a,i7,a,i7,a)') &
1445  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1446  write(qunit,'(a)')'<PointData>'
1447  do iw=1,nw+nwauxio
1448  if(iw<=nw) then
1449  if(.not.w_write(iw)) cycle
1450  end if
1451  write(qunit,'(a,a,a,i16,a)')&
1452  '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1453  '" format="appended" offset="',offset,'">'
1454  write(qunit,'(a)')'</DataArray>'
1455  offset=offset+length+size_int
1456  end do
1457  write(qunit,'(a)')'</PointData>'
1458  write(qunit,'(a)')'<Points>'
1459  write(qunit,'(a,i16,a)') &
1460  '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1461  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1462  offset=offset+length_coords+size_int
1463  write(qunit,'(a)')'</Points>'
1464  else
1465  ! we write out every grid as one VTK PIECE
1466  write(qunit,'(a,i7,a,i7,a)') &
1467  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
1468  write(qunit,'(a)')'<CellData>'
1469  do iw=1,nw+nwauxio
1470  if(iw<=nw) then
1471  if(.not.w_write(iw)) cycle
1472  end if
1473  write(qunit,'(a,a,a,i16,a)')&
1474  '<DataArray type="Float64" Name="',trim(wnamei(iw)), &
1475  '" format="appended" offset="',offset,'">'
1476  write(qunit,'(a)')'</DataArray>'
1477  offset=offset+lengthcc+size_int
1478  end do
1479  write(qunit,'(a)')'</CellData>'
1480  write(qunit,'(a)')'<Points>'
1481  write(qunit,'(a,i16,a)') &
1482  '<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
1483  ! write cell corner coordinates in a backward dimensional loop, always 3D output
1484  offset=offset+length_coords+size_int
1485  write(qunit,'(a)')'</Points>'
1486  end if
1487  write(qunit,'(a)')'<Cells>'
1488  ! connectivity part
1489  write(qunit,'(a,i16,a)')&
1490  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
1491  offset=offset+length_conn+size_int
1492 
1493  ! offsets data array
1494  write(qunit,'(a,i16,a)') &
1495  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
1496  offset=offset+length_offsets+size_int
1497 
1498  ! VTK cell type data array
1499  write(qunit,'(a,i16,a)') &
1500  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
1501  offset=offset+size_int+nc*size_int
1502 
1503  write(qunit,'(a)')'</Cells>'
1504 
1505  write(qunit,'(a)')'</Piece>'
1506  end do
1507  end do
1508  end if
1509 
1510  write(qunit,'(a)')'</UnstructuredGrid>'
1511  write(qunit,'(a)')'<AppendedData encoding="raw">'
1512  close(qunit)
1513  open(qunit,file=filename,access='stream',form='unformatted',position='append')
1514  buf='_'
1515  write(qunit) trim(buf)
1516 
1517  do morton_no=morton_start(0),morton_stop(0)
1518  if(.not. morton_aim(morton_no)) cycle
1519  igrid=sfc_to_igrid(morton_no)
1520  call calc_x(igrid,xc,xcc)
1521  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1522  ixc^l,ixcc^l,.true.)
1523  do iw=1,nw+nwauxio
1524  if(iw<=nw) then
1525  if(.not.w_write(iw)) cycle
1526  end if
1527  if(cell_corner) then
1528  write(qunit) length
1529  write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1530  else
1531  write(qunit) lengthcc
1532  write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1533  end if
1534  end do
1535 
1536  write(qunit) length_coords
1537  {do ix^db=ixcmin^db,ixcmax^db \}
1538  x_vtk(1:3)=zero;
1539  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1540  do k=1,3
1541  write(qunit) x_vtk(k)
1542  end do
1543  {end do \}
1544 
1545  write(qunit) length_conn
1546  {do ix^db=1,nx^db\}
1547  {^ifoned write(qunit)ix1-1,ix1 \}
1548  {^iftwod
1549  write(qunit)(ix2-1)*nxc1+ix1-1, &
1550  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1551  \}
1552  {^ifthreed
1553  write(qunit)&
1554  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1555  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1556  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1557  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1558  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1559  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1560  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1561  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1562  \}
1563  {end do\}
1564 
1565  write(qunit) length_offsets
1566  do icel=1,nc
1567  write(qunit) icel*(2**^nd)
1568  end do
1569 
1570  {^ifoned vtk_type=3 \}
1571  {^iftwod vtk_type=8 \}
1572  {^ifthreed vtk_type=11 \}
1573  write(qunit) size_int*nc
1574  do icel=1,nc
1575  write(qunit) vtk_type
1576  end do
1577  end do
1578  allocate(intstatus(mpi_status_size,1))
1579  if(npe>1)then
1580  ixccmin^d=ixmlo^d; ixccmax^d=ixmhi^d;
1581  ixcmin^d=ixmlo^d-1; ixcmax^d=ixmhi^d;
1582  do ipe=1, npe-1
1583  do morton_no=morton_start(ipe),morton_stop(ipe)
1584  if(.not. morton_aim(morton_no)) cycle
1585  itag=morton_no
1586  call mpi_recv(xc_tmp,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1587  if(cell_corner) then
1588  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1589  else
1590  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1591  end if
1592  do iw=1,nw+nwauxio
1593  if(iw<=nw) then
1594  if(.not.w_write(iw)) cycle
1595  end if
1596  if(cell_corner) then
1597  write(qunit) length
1598  write(qunit) {(|}wc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
1599  else
1600  write(qunit) lengthcc
1601  write(qunit) {(|}wcc_tmp(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
1602  end if
1603  end do
1604 
1605  write(qunit) length_coords
1606  {do ix^db=ixcmin^db,ixcmax^db \}
1607  x_vtk(1:3)=zero;
1608  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
1609  do k=1,3
1610  write(qunit) x_vtk(k)
1611  end do
1612  {end do \}
1613 
1614  write(qunit) length_conn
1615  {do ix^db=1,nx^db\}
1616  {^ifoned write(qunit)ix1-1,ix1 \}
1617  {^iftwod
1618  write(qunit)(ix2-1)*nxc1+ix1-1, &
1619  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1620  \}
1621  {^ifthreed
1622  write(qunit)&
1623  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1624  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1625  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1626  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1627  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1628  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1629  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1630  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1631  \}
1632  {end do\}
1633 
1634  write(qunit) length_offsets
1635  do icel=1,nc
1636  write(qunit) icel*(2**^nd)
1637  end do
1638  {^ifoned vtk_type=3 \}
1639  {^iftwod vtk_type=8 \}
1640  {^ifthreed vtk_type=11 \}
1641  write(qunit) size_int*nc
1642  do icel=1,nc
1643  write(qunit) vtk_type
1644  end do
1645  end do
1646  end do
1647  end if
1648  close(qunit)
1649  open(qunit,file=filename,status='unknown',form='formatted',position='append')
1650  write(qunit,'(a)')'</AppendedData>'
1651  write(qunit,'(a)')'</VTKFile>'
1652  close(qunit)
1653  deallocate(intstatus)
1654 end if
1655 
1656 deallocate(morton_aim,morton_aim_p)
1657 if (npe>1) then
1658  call mpi_barrier(icomm,ierrmpi)
1659 end if
1660 
1661 end subroutine unstructuredvtkb64
1662 !====================================================================================
1663 subroutine save_connvtk(qunit,igrid)
1664 
1665 ! this saves the basic line, pixel and voxel connectivity,
1666 ! as used by VTK file outputs for unstructured grid
1667 
1669 
1670 integer, intent(in) :: qunit, igrid
1671 
1672 integer :: nx^D, nxC^D, ix^D
1673 !-----------------------------------------------------------------------------
1674 nx^d=ixmhi^d-ixmlo^d+1;
1675 nxc^d=nx^d+1;
1676 
1677 {do ix^db=1,nx^db\}
1678  {^ifoned write(qunit,'(2(i7,1x))')ix1-1,ix1 \}
1679  {^iftwod
1680  write(qunit,'(4(i7,1x))')(ix2-1)*nxc1+ix1-1, &
1681  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
1682  \}
1683  {^ifthreed
1684  write(qunit,'(8(i7,1x))')&
1685  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
1686  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1687  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1688  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
1689  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
1690  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
1691  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
1692  ix3*nxc2*nxc1+ ix2*nxc1+ix1
1693  \}
1694 {end do\}
1695 
1696 end subroutine save_connvtk
1697 !=============================================================================
1698 subroutine imagedatavtk_mpi(qunit)
1699 
1700 ! output for vti format to paraview, non-binary version output
1701 ! parallel, uses calc_grid to compute nwauxio variables
1702 ! allows renormalizing using convert factors
1703 ! allows skipping of w_write selected variables
1704 
1705 ! implementation such that length of ASCII output is identical when
1706 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1707 
1710 use mod_calculate_xw
1711 
1712 integer, intent(in) :: qunit
1713 
1714 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1715 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1716 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1717 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1718 
1719 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1720 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1721 double precision, dimension(0:nw+nwauxio) :: normconv
1722 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1723 integer:: NumGridsOnLevel(1:nlevelshi)
1724 integer :: nx^D
1725 
1726 character(len=80):: filename
1727 integer :: filenr
1728 
1729 integer, allocatable :: intstatus(:,:)
1730 logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
1731 
1732 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1733 character(len=1024) :: outfilehead
1734 
1735 logical :: fileopen
1736 integer :: itag,ipe,Morton_no,Morton_length
1737 integer :: ixrvC^L, ixrvCC^L, siz_ind, ind_send(5*^ND), ind_recv(5*^ND)
1738 double precision :: origin(1:3), spacing(1:3)
1739 integer :: wholeExtent(1:6), ig^D
1740 type(tree_node_ptr) :: tree
1741 !-----------------------------------------------------------------------------
1742 if(levmin/=levmax) call mpistop('ImageData can only be used when levmin=levmax')
1743 
1744 normconv(0) = length_convert_factor
1745 normconv(1:nw) = w_convert_factor
1746 siz_ind=5*^nd
1747 morton_length=morton_stop(npe-1)-morton_start(0)+1
1748 allocate(morton_aim(morton_start(0):morton_stop(npe-1)))
1749 allocate(morton_aim_p(morton_start(0):morton_stop(npe-1)))
1750 morton_aim=.false.
1751 morton_aim_p=.false.
1752 do morton_no=morton_start(mype),morton_stop(mype)
1753  igrid=sfc_to_igrid(morton_no)
1754  level=node(plevel_,igrid)
1755  ! we can clip parts of the grid away, select variables, levels etc.
1756  if(writelevel(level)) then
1757  ! only output a grid when fully within clipped region selected
1758  ! by writespshift array
1759  if(({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1760  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1761  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
1762  morton_aim_p(morton_no)=.true.
1763  end if
1764  end if
1765 end do
1766 call mpi_allreduce(morton_aim_p,morton_aim,morton_length,mpi_logical,mpi_lor,&
1767  icomm,ierrmpi)
1768 
1769 
1770 if (mype /= 0) then
1771  do morton_no=morton_start(mype),morton_stop(mype)
1772  if(.not. morton_aim(morton_no)) cycle
1773  igrid=sfc_to_igrid(morton_no)
1774  call calc_x(igrid,xc,xcc)
1775  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1776  ixc^l,ixcc^l,.true.)
1777  tree%node => igrid_to_node(igrid, mype)%node
1778  {^d& ig^d = tree%node%ig^d; }
1779  itag=morton_no
1780  ind_send=(/ ixc^l,ixcc^l, ig^d /)
1781  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
1782  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
1783  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
1784  end do
1785 
1786 else
1787 
1788  inquire(qunit,opened=fileopen)
1789  if(.not.fileopen)then
1790  ! generate filename
1791  filenr=snapshotini
1792  if (autoconvert) filenr=snapshotnext
1793  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vti"
1794  ! Open the file for the header part
1795  open(qunit,file=filename,status='unknown',form='formatted')
1796  endif
1797 
1798 call getheadernames(wnamei,xandwnamei,outfilehead)
1799 
1800 ! number of cells per grid.
1801 nx^d=ixmhi^d-ixmlo^d+1;
1802 
1803 origin = 0
1804 {^d& origin(^d) = xprobmin^d*normconv(0); }
1805 spacing = zero
1806 {^d&spacing(^d) = dxlevel(^d)*normconv(0); }
1807 
1808 wholeextent = 0
1809 ! if we use writespshift, the whole extent has to be calculated:
1810 {^d&wholeextent(^d*2-1) = nx^d * ceiling(((xprobmax^d-xprobmin^d)*writespshift(^d,1)) &
1811  /(nx^d*dxlevel(^d))) \}
1812 {^d&wholeextent(^d*2) = nx^d * floor(((xprobmax^d-xprobmin^d)*(1.0d0-writespshift(^d,2))) &
1813  /(nx^d*dxlevel(^d))) \}
1814 
1815 ! generate xml header
1816 write(qunit,'(a)')'<?xml version="1.0"?>'
1817 write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
1818 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1819 write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')' <ImageData Origin="',&
1820  origin,'" WholeExtent="',wholeextent,'" Spacing="',spacing,'">'
1821  write(qunit,'(a)')'<FieldData>'
1822  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1823  'NumberOfTuples="1" format="ascii">'
1824  write(qunit,*) real(global_time*time_convert_factor)
1825  write(qunit,'(a)')'</DataArray>'
1826  write(qunit,'(a)')'</FieldData>'
1827 
1828 ! write the data from proc 0
1829 do morton_no=morton_start(0),morton_stop(0)
1830  if(.not. morton_aim(morton_no)) cycle
1831  igrid=sfc_to_igrid(morton_no)
1832  tree%node => igrid_to_node(igrid, 0)%node
1833  {^d& ig^d = tree%node%ig^d; }
1834  call calc_x(igrid,xc,xcc)
1835  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1836  ixc^l,ixcc^l,.true.)
1837  call write_vti(qunit,ixg^ll,ixc^l,ixcc^l,ig^d,&
1838  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1839 end do
1840 
1841 if(npe>1)then
1842  allocate(intstatus(mpi_status_size,1))
1843  do ipe=1, npe-1
1844  do morton_no=morton_start(ipe),morton_stop(ipe)
1845  if(.not. morton_aim(morton_no)) cycle
1846  itag=morton_no
1847  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1848  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
1849  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
1850  ig^d=ind_recv(4*^nd+^d);
1851  call mpi_recv(wc_tmp,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1852  call mpi_recv(wcc_tmp,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
1853  call write_vti(qunit,ixg^ll,ixrvc^l,ixrvcc^l,ig^d,&
1854  nx^d,normconv,wnamei,wc_tmp,wcc_tmp)
1855  end do
1856  end do
1857 end if
1858 
1859 write(qunit,'(a)')'</ImageData>'
1860 write(qunit,'(a)')'</VTKFile>'
1861 close(qunit)
1862 if(npe>1) deallocate(intstatus)
1863 endif
1864 
1865 deallocate(morton_aim,morton_aim_p)
1866 if (npe>1) then
1867  call mpi_barrier(icomm,ierrmpi)
1868 endif
1869 
1870 end subroutine imagedatavtk_mpi
1871 !============================================================================
1872 subroutine punstructuredvtk_mpi(qunit)
1873 
1874 ! Write one pvtu and vtu files for each processor
1875 ! Otherwise like unstructuredvtk_mpi
1876 
1879 use mod_calculate_xw
1880 
1881 integer, intent(in) :: qunit
1882 !
1883 double precision, dimension(0:nw+nwauxio) :: normconv
1884 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
1885 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
1886 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1887 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1888 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
1889 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
1890 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
1891 character(len=1024) :: outfilehead
1892 integer :: nx^D,nxC^D,nc,np, igrid,ixC^L,ixCC^L,level,Morton_no
1893 character(len=80) :: pfilename
1894 integer :: filenr
1895 logical :: fileopen,conv_grid
1896 !----------------------------------------------------------------------------
1897 
1898 ! Write pvtu-file:
1899 if (mype==0) then
1900  call write_pvtu(qunit)
1901 endif
1902 ! Now write the Source files:
1903 
1904 inquire(qunit,opened=fileopen)
1905 if(.not.fileopen)then
1906  ! generate filename
1907  filenr=snapshotini
1908  if (autoconvert) filenr=snapshotnext
1909  ! Open the file for the header part
1910  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
1911  open(qunit,file=pfilename,status='unknown',form='formatted')
1912 endif
1913 ! generate xml header
1914 write(qunit,'(a)')'<?xml version="1.0"?>'
1915 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
1916 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
1917 write(qunit,'(a)')' <UnstructuredGrid>'
1918 write(qunit,'(a)')'<FieldData>'
1919 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
1920  'NumberOfTuples="1" format="ascii">'
1921 write(qunit,*) real(global_time*time_convert_factor)
1922 write(qunit,'(a)')'</DataArray>'
1923 write(qunit,'(a)')'</FieldData>'
1924 
1925 call getheadernames(wnamei,xandwnamei,outfilehead)
1926 
1927 ! number of cells, number of corner points, per grid.
1928 nx^d=ixmhi^d-ixmlo^d+1;
1929 nxc^d=nx^d+1;
1930 nc={nx^d*}
1931 np={nxc^d*}
1932 
1933 ! Note: using the w_write, writelevel, writespshift
1934 ! we can clip parts of the grid away, select variables, levels etc.
1935 do level=levmin,levmax
1936  if (.not.writelevel(level)) cycle
1937  do morton_no=morton_start(mype),morton_stop(mype)
1938  igrid=sfc_to_igrid(morton_no)
1939  if (node(plevel_,igrid)/=level) cycle
1940 
1941  ! only output a grid when fully within clipped region selected
1942  ! by writespshift array
1943  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
1944  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
1945  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
1946  if (.not.conv_grid) cycle
1947 
1948  call calc_x(igrid,xc,xcc)
1949  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
1950  ixc^l,ixcc^l,.true.)
1951 
1952  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
1953  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
1954  end do ! Morton_no loop
1955 end do ! level loop
1956 
1957  write(qunit,'(a)')' </UnstructuredGrid>'
1958  write(qunit,'(a)')'</VTKFile>'
1959  close(qunit)
1960 
1961 if (npe>1) then
1962  call mpi_barrier(icomm,ierrmpi)
1963 endif
1964 end subroutine punstructuredvtk_mpi
1965 !============================================================================
1966 subroutine unstructuredvtk_mpi(qunit)
1967 
1968 ! output for vtu format to paraview, non-binary version output
1969 ! parallel, uses calc_grid to compute nwauxio variables
1970 ! allows renormalizing using convert factors
1971 ! allows skipping of w_write selected variables
1972 
1973 ! implementation such that length of ASCII output is identical when
1974 ! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
1975 
1978 use mod_calculate_xw
1979 
1980 integer, intent(in) :: qunit
1981 
1982 double precision :: x_VTK(1:3)
1983 
1984 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
1985 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
1986 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
1987 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
1988 
1989 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
1990 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
1991 double precision, dimension(0:nw+nwauxio) :: normconv
1992 integer:: igrid,iigrid,level,ixC^L,ixCC^L
1993 integer:: NumGridsOnLevel(1:nlevelshi)
1994 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,ix^D
1995 
1996 character(len=80):: filename
1997 integer :: filenr
1998 
1999 integer, allocatable :: intstatus(:,:)
2000 
2001 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2002 character(len=1024) :: outfilehead
2003 
2004 logical :: fileopen,conv_grid,cond_grid_recv
2005 integer :: itag,ipe,Morton_no,siz_ind
2006 integer :: ind_send(4*^ND),ind_recv(4*^ND)
2007 integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvC^L,ixrvCC^L
2008 !-----------------------------------------------------------------------------
2009 if (mype==0) then
2010  inquire(qunit,opened=fileopen)
2011  if(.not.fileopen)then
2012  ! generate filename
2013  filenr=snapshotini
2014  if (autoconvert) filenr=snapshotnext
2015  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".vtu"
2016  ! Open the file for the header part
2017  open(qunit,file=filename,status='unknown',form='formatted')
2018  endif
2019  ! generate xml header
2020  write(qunit,'(a)')'<?xml version="1.0"?>'
2021  write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2022  write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2023  write(qunit,'(a)')'<UnstructuredGrid>'
2024  write(qunit,'(a)')'<FieldData>'
2025  write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2026  'NumberOfTuples="1" format="ascii">'
2027  write(qunit,*) real(global_time*time_convert_factor)
2028  write(qunit,'(a)')'</DataArray>'
2029  write(qunit,'(a)')'</FieldData>'
2030 end if
2031 
2032 call getheadernames(wnamei,xandwnamei,outfilehead)
2033 ! number of cells, number of corner points, per grid.
2034 nx^d=ixmhi^d-ixmlo^d+1;
2035 nxc^d=nx^d+1;
2036 nc={nx^d*}
2037 np={nxc^d*}
2038 
2039 ! all slave processors send their minmal/maximal levels
2040 if (mype/=0) then
2041  if (morton_stop(mype)==0) call mpistop("nultag")
2042  itag=1000*morton_stop(mype)
2043  !print *,'ype,itag for levmin=',mype,itag,levmin
2044  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
2045  itag=2000*morton_stop(mype)
2046  !print *,'mype,itag for levmax=',mype,itag,levmax
2047  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
2048 end if
2049 
2050 
2051 ! Note: using the w_write, writelevel, writespshift
2052 ! we can clip parts of the grid away, select variables, levels etc.
2053 do level=levmin,levmax
2054  if (.not.writelevel(level)) cycle
2055  do morton_no=morton_start(mype),morton_stop(mype)
2056  igrid=sfc_to_igrid(morton_no)
2057  if (mype/=0)then
2058  itag=morton_no
2059  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2060  itag=igrid
2061  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2062  end if
2063  if (node(plevel_,igrid)/=level) cycle
2064 
2065  ! only output a grid when fully within clipped region selected
2066  ! by writespshift array
2067  conv_grid=({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2068  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2069  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})
2070  if (mype/=0)then
2071  call mpi_send(conv_grid,1,mpi_logical,0,itag,icomm,ierrmpi)
2072  end if
2073  if (.not.conv_grid) cycle
2074 
2075  call calc_x(igrid,xc,xcc)
2076  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2077  ixc^l,ixcc^l,.true.)
2078 
2079  if (mype/=0) then
2080  itag=morton_no
2081  ind_send=(/ ixc^l,ixcc^l /)
2082  siz_ind=4*^nd
2083  call mpi_send(ind_send,siz_ind,mpi_integer, 0,itag,icomm,ierrmpi)
2084  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2085 
2086  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
2087  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2088  itag=igrid
2089  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2090  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
2091  else
2092  call write_vtk(qunit,ixg^ll,ixc^l,ixcc^l,igrid,nc,np,nx^d,nxc^d,&
2093  normconv,wnamei,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp)
2094  end if
2095  end do ! Morton_no loop
2096 end do ! level loop
2097 
2098 
2099 if (mype==0) then
2100  allocate(intstatus(mpi_status_size,1))
2101  if(npe>1)then
2102  do ipe=1,npe-1
2103  itag=1000*morton_stop(ipe)
2104  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2105  !!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
2106  itag=2000*morton_stop(ipe)
2107  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2108  !!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
2109  do level=levmin_recv,levmax_recv
2110  if (.not.writelevel(level)) cycle
2111  do morton_no=morton_start(ipe),morton_stop(ipe)
2112  itag=morton_no
2113  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2114  itag=igrid_recv
2115  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2116  if (level_recv/=level) cycle
2117 
2118  call mpi_recv(cond_grid_recv,1,mpi_logical, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2119  if(.not.cond_grid_recv)cycle
2120 
2121  itag=morton_no
2122  siz_ind=4*^nd
2123  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2124  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2125  ixrvccmin^d=ind_recv(2*^nd+^d);ixrvccmax^d=ind_recv(3*^nd+^d);
2126  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2127  ,icomm,intstatus(:,1),ierrmpi)
2128 
2129  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2130  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2131 
2132  itag=igrid_recv
2133  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2134  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2135  call write_vtk(qunit,ixg^ll,ixrvc^l,ixrvcc^l,igrid_recv,&
2136  nc,np,nx^d,nxc^d,normconv,wnamei,&
2137  xc_tmp_recv,xcc_tmp_recv,wc_tmp_recv,wcc_tmp_recv)
2138  enddo ! Morton_no loop
2139  enddo ! level loop
2140  enddo ! processor loop
2141  endif ! multiple processors
2142  write(qunit,'(a)')'</UnstructuredGrid>'
2143  write(qunit,'(a)')'</VTKFile>'
2144  close(qunit)
2145 endif
2146 
2147 if (npe>1) then
2148  call mpi_barrier(icomm,ierrmpi)
2149  if(mype==0)deallocate(intstatus)
2150 endif
2151 
2152 end subroutine unstructuredvtk_mpi
2153 !============================================================================
2154 subroutine write_vtk(qunit,ixI^L,ixC^L,ixCC^L,igrid,nc,np,nx^D,nxC^D,&
2155  normconv,wnamei,xC,xCC,wC,wCC)
2156 
2158 
2159 integer, intent(in) :: qunit
2160 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
2161 integer, intent(in) :: igrid,nc,np,nx^D,nxC^D
2162 double precision, intent(in) :: normconv(0:nw+nwauxio)
2163 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
2164 
2165 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2166 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2167 
2168 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
2169 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
2170 
2171 double precision :: x_VTK(1:3)
2172 integer :: iw,ix^D,icel,VTK_type
2173 !----------------------------------------------------------------------------
2174 
2175 select case(convert_type)
2176  case('vtumpi','pvtumpi')
2177  ! we write out every grid as one VTK PIECE
2178  write(qunit,'(a,i7,a,i7,a)') &
2179  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2180  write(qunit,'(a)')'<PointData>'
2181  do iw=1,nw+nwauxio
2182  if(iw<=nw) then
2183  if(.not.w_write(iw)) cycle
2184  endif
2185 
2186  write(qunit,'(a,a,a)')&
2187  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2188  write(qunit,'(200(1pe14.6))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
2189  write(qunit,'(a)')'</DataArray>'
2190  enddo
2191  write(qunit,'(a)')'</PointData>'
2192 
2193  write(qunit,'(a)')'<Points>'
2194  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2195  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2196  {do ix^db=ixcmin^db,ixcmax^db \}
2197  x_vtk(1:3)=zero;
2198  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
2199  write(qunit,'(3(1pe14.6))') x_vtk
2200  {end do \}
2201  write(qunit,'(a)')'</DataArray>'
2202  write(qunit,'(a)')'</Points>'
2203 
2204  case('vtuCCmpi','pvtuCCmpi')
2205  ! we write out every grid as one VTK PIECE
2206  write(qunit,'(a,i7,a,i7,a)') &
2207  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2208  write(qunit,'(a)')'<CellData>'
2209  do iw=1,nw+nwauxio
2210  if(iw<=nw) then
2211  if(.not.w_write(iw)) cycle
2212  endif
2213  write(qunit,'(a,a,a)')&
2214  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2215  write(qunit,'(200(1pe14.6))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
2216  write(qunit,'(a)')'</DataArray>'
2217  enddo
2218  write(qunit,'(a)')'</CellData>'
2219 
2220  write(qunit,'(a)')'<Points>'
2221  write(qunit,'(a)')'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
2222  ! write cell corner coordinates in a backward dimensional loop, always 3D output
2223  {do ix^db=ixcmin^db,ixcmax^db \}
2224  x_vtk(1:3)=zero;
2225  x_vtk(1:ndim)=xc(ix^d,1:ndim)*normconv(0);
2226  write(qunit,'(3(1pe14.6))') x_vtk
2227  {end do \}
2228  write(qunit,'(a)')'</DataArray>'
2229  write(qunit,'(a)')'</Points>'
2230 end select
2231 
2232 write(qunit,'(a)')'<Cells>'
2233 
2234 ! connectivity part
2235 write(qunit,'(a)')'<DataArray type="Int32" Name="connectivity" format="ascii">'
2236 call save_connvtk(qunit,igrid)
2237 write(qunit,'(a)')'</DataArray>'
2238 
2239 ! offsets data array
2240 write(qunit,'(a)')'<DataArray type="Int32" Name="offsets" format="ascii">'
2241 do icel=1,nc
2242  write(qunit,'(i7)') icel*(2**^nd)
2243 end do
2244 write(qunit,'(a)')'</DataArray>'
2245 
2246 ! VTK cell type data array
2247 write(qunit,'(a)')'<DataArray type="Int32" Name="types" format="ascii">'
2248 ! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
2249 {^ifoned vtk_type=3 \}
2250 {^iftwod vtk_type=8 \}
2251 {^ifthreed vtk_type=11 \}
2252 do icel=1,nc
2253  write(qunit,'(i2)') vtk_type
2254 enddo
2255 write(qunit,'(a)')'</DataArray>'
2256 
2257 write(qunit,'(a)')'</Cells>'
2258 
2259 write(qunit,'(a)')'</Piece>'
2260 
2261 end subroutine write_vtk
2262 !============================================================================
2263 subroutine write_vti(qunit,ixI^L,ixC^L,ixCC^L,ig^D,nx^D,&
2264  normconv,wnamei,wC,wCC)
2266 
2267 integer, intent(in) :: qunit
2268 integer, intent(in) :: ixI^L,ixC^L,ixCC^L
2269 integer, intent(in) :: ig^D,nx^D
2270 double precision, intent(in) :: normconv(0:nw+nwauxio)
2271 character(len=name_len), intent(in):: wnamei(1:nw+nwauxio)
2272 
2273 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC
2274 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC
2275 
2276 integer :: iw,ix^D
2277 integer :: extent(1:6)
2278 !----------------------------------------------------------------------------
2279 
2280 extent = 0
2281 {^d& extent(^d*2-1) = (ig^d-1) * nx^d; }
2282 {^d& extent(^d*2) = (ig^d) * nx^d; }
2283 
2284 
2285 select case(convert_type)
2286  case('vtimpi','pvtimpi')
2287  ! we write out every grid as one VTK PIECE
2288  write(qunit,'(a,6(i10),a)') &
2289  '<Piece Extent="',extent,'">'
2290  write(qunit,'(a)')'<PointData>'
2291  do iw=1,nw+nwauxio
2292  if(iw<=nw) then
2293  if(.not.w_write(iw)) cycle
2294  endif
2295 
2296  write(qunit,'(a,a,a)')&
2297  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2298  write(qunit,'(200(1pe20.12))') {(|}wc(ix^d,iw)*normconv(iw),{ix^d=ixcmin^d,ixcmax^d)}
2299  write(qunit,'(a)')'</DataArray>'
2300  enddo
2301  write(qunit,'(a)')'</PointData>'
2302 
2303  case('vtiCCmpi','pvtiCCmpi')
2304  ! we write out every grid as one VTK PIECE
2305  write(qunit,'(a,6(i10),a)') &
2306  '<Piece Extent="',extent,'">'
2307  write(qunit,'(a)')'<CellData>'
2308  do iw=1,nw+nwauxio
2309  if(iw<=nw) then
2310  if(.not.w_write(iw)) cycle
2311  endif
2312  write(qunit,'(a,a,a)')&
2313  '<DataArray type="Float64" Name="',trim(wnamei(iw)),'" format="ascii">'
2314  write(qunit,'(200(1pe20.12))') {(|}wcc(ix^d,iw)*normconv(iw),{ix^d=ixccmin^d,ixccmax^d)}
2315  write(qunit,'(a)')'</DataArray>'
2316  enddo
2317  write(qunit,'(a)')'</CellData>'
2318 end select
2319 
2320 write(qunit,'(a)')'</Piece>'
2321 
2322 end subroutine write_vti
2323 !=============================================================================
2324 subroutine write_pvtu(qunit)
2325 
2327 use mod_calculate_xw
2328 
2329 integer, intent(in) :: qunit
2330 
2331 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio),outtype
2332 character(len=1024) :: outfilehead
2333 character(len=80) :: filename,pfilename
2334 integer :: filenr,iw,ipe,iscalars
2335 logical :: fileopen
2336 
2337 select case(convert_type)
2338 case('pvtumpi','pvtuBmpi')
2339  outtype="PPointData"
2340 case('pvtuCCmpi','pvtuBCCmpi')
2341  outtype="PCellData"
2342 end select
2343 inquire(qunit,opened=fileopen)
2344 if(.not.fileopen)then
2345  ! generate filename
2346  filenr=snapshotini
2347  if (autoconvert) filenr=snapshotnext
2348  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".pvtu"
2349  ! Open the file
2350  open(qunit,file=filename,status='unknown',form='formatted')
2351 endif
2352 
2353 call getheadernames(wnamei,xandwnamei,outfilehead)
2354 
2355 ! Get the default selection:
2356 iscalars=1
2357 do iw=nw,1, -1
2358  if (w_write(iw)) iscalars=iw
2359 end do
2360 
2361 
2362 ! generate xml header
2363 write(qunit,'(a)')'<?xml version="1.0"?>'
2364 write(qunit,'(a)',advance='no') '<VTKFile type="PUnstructuredGrid"'
2365 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2366 write(qunit,'(a)')' <PUnstructuredGrid GhostLevel="0">'
2367 ! Either celldata or pointdata:
2368 write(qunit,'(a,a,a,a,a)')&
2369  ' <',trim(outtype),' Scalars="',trim(wnamei(iscalars))//'">'
2370 do iw=1,nw
2371  if(.not.w_write(iw))cycle
2372  write(qunit,'(a,a,a)')&
2373  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
2374 end do
2375 do iw=nw+1,nw+nwauxio
2376  write(qunit,'(a,a,a)')&
2377  ' <PDataArray type="Float32" Name="',trim(wnamei(iw)),'"/>'
2378 end do
2379 write(qunit,'(a,a,a)')' </',trim(outtype),'>'
2380 
2381 write(qunit,'(a)')' <PPoints>'
2382 write(qunit,'(a)')' <PDataArray type="Float32" NumberOfComponents="3"/>'
2383 write(qunit,'(a)')' </PPoints>'
2384 
2385 do ipe=0,npe-1
2386  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename(&
2387  index(base_filename, '/', back = .true.)+1:&
2388  len(base_filename))),filenr,"p",&
2389  ipe,".vtu"
2390  write(qunit,'(a,a,a)')' <Piece Source="',trim(pfilename),'"/>'
2391 end do
2392 write(qunit,'(a)')' </PUnstructuredGrid>'
2393 write(qunit,'(a)')'</VTKFile>'
2394 
2395 close(qunit)
2396 
2397 end subroutine write_pvtu
2398 !=============================================================================
2399 subroutine tecplot_mpi(qunit)
2400 
2401 ! output for tecplot (ASCII format)
2402 ! parallel, uses calc_grid to compute nwauxio variables
2403 ! allows renormalizing using convert factors
2404 
2405 ! the current implementation is such that tecplotmpi and tecplotCCmpi will
2406 ! create different length output ASCII files when used on 1 versus multiple CPUs
2407 ! in fact, on 1 CPU, there will be as many zones as there are levels
2408 ! on multiple CPUs, there will be a number of zones up to the number of
2409 ! levels times the number of CPUs (can be less, when some level not on a CPU)
2410 
2413 use mod_calculate_xw
2414 
2415 integer, intent(in) :: qunit
2416 
2417 integer:: igrid,iigrid,level,igonlevel,iw,idim,ix^D
2418 integer:: NumGridsOnLevel(1:nlevelshi)
2419 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,ixC^L,ixCC^L
2420 integer :: nodesonlevelmype,elemsonlevelmype
2421 
2422 integer :: nodes, elems
2423 
2424 integer, allocatable :: intstatus(:,:)
2425 
2426 double precision :: x_TEC(ndim), w_TEC(nw+nwauxio)
2427 
2428 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP,xC_TMP_recv
2429 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP,xCC_TMP_recv
2430 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2431 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2432 
2433 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP,wC_TMP_recv
2434 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
2435 double precision, dimension(0:nw+nwauxio) :: normconv
2436 logical :: fileopen,first
2437 integer :: itag,Morton_no,ipe,levmin_recv,levmax_recv,igrid_recv,level_recv
2438 integer :: ixrvC^L,ixrvCC^L
2439 integer :: ind_send(2*^ND),ind_recv(2*^ND),siz_ind,igonlevel_recv
2440 integer :: NumGridsOnLevel_mype(1:nlevelshi,0:npe-1)
2441 character(len=80) :: filename
2442 integer :: filenr
2443 character(len=1024) :: tecplothead
2444 
2445 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2446 character(len=1024) :: outfilehead
2447 !-----------------------------------------------------------------------------
2448 if(nw/=count(w_write(1:nw)))then
2449  if(mype==0) print *,'tecplot_mpi does not use w_write=F'
2450  call mpistop('w_write, tecplot')
2451 end if
2452 
2453 if(nocartesian)then
2454  if(mype==0) print *,'tecplot_mpi with nocartesian'
2455 endif
2456 
2457 master_cpu_open : if (mype == 0) then
2458  inquire(qunit,opened=fileopen)
2459  if (.not.fileopen) then
2460  ! generate filename
2461  filenr=snapshotini
2462  if (autoconvert) filenr=snapshotnext
2463  write(filename,'(a,i4.4,a)') trim(base_filename),filenr,".plt"
2464  open(qunit,file=filename,status='unknown')
2465  end if
2466 
2467  call getheadernames(wnamei,xandwnamei,outfilehead)
2468 
2469  write(tecplothead,'(a)') "VARIABLES = "//trim(outfilehead)
2470  write(qunit,'(a)') tecplothead(1:len_trim(tecplothead))
2471 end if master_cpu_open
2472 
2473 
2474 ! determine overall number of grids per level, and the same info per CPU
2475 numgridsonlevel(1:nlevelshi)=0
2476 do level=levmin,levmax
2477  numgridsonlevel(level)=0
2478  do morton_no=morton_start(mype),morton_stop(mype)
2479  igrid = sfc_to_igrid(morton_no)
2480  if (node(plevel_,igrid)/=level) cycle
2481  numgridsonlevel(level)=numgridsonlevel(level)+1
2482  end do
2483  numgridsonlevel_mype(level,0:npe-1)=0
2484  numgridsonlevel_mype(level,mype) = numgridsonlevel(level)
2485  call mpi_allreduce(mpi_in_place,numgridsonlevel_mype(level,0:npe-1),npe,mpi_integer,&
2486  mpi_max,icomm,ierrmpi)
2487  call mpi_allreduce(mpi_in_place,numgridsonlevel(level),1,mpi_integer,mpi_sum, &
2488  icomm,ierrmpi)
2489 end do
2490 
2491 
2492 !!do level=levmin,levmax
2493 !! print *,'mype, level en NumGridsOnLevel_mype(level,0:npe-1)=', &
2494 !! mype,level,NumGridsOnLevel_mype(level,0:npe-1)
2495 !! print *,'mype, level en NumGridsOnLevel(level)=', &
2496 !! mype,level,NumGridsOnLevel(level)
2497 !!enddo
2498 
2499 
2500 nx^d=ixmhi^d-ixmlo^d+1;
2501 nxc^d=nx^d+1;
2502 
2503 if(mype==0.and.npe>1) allocate(intstatus(mpi_status_size,1))
2504 
2505 {^ifoned
2506 if(convert_type=='teclinempi') then
2507  nodes=0
2508  elems=0
2509  do level=levmin,levmax
2510  nodes=nodes + numgridsonlevel(level)*{nxc^d*}
2511  elems=elems + numgridsonlevel(level)*{nx^d*}
2512  enddo
2513 
2514  if (mype==0) write(qunit,"(a,i7,a,1pe12.5,a)") &
2515  'ZONE T="all levels", I=',elems, &
2516  ', SOLUTIONTIME=',global_time*time_convert_factor,', F=POINT'
2517 
2518  igonlevel=0
2519  do morton_no=morton_start(mype),morton_stop(mype)
2520  igrid = sfc_to_igrid(morton_no)
2521  call calc_x(igrid,xc,xcc)
2522  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,ixc^l,ixcc^l,.true.)
2523  if (mype==0) then
2524  {do ix^db=ixccmin^db,ixccmax^db\}
2525  x_tec(1:ndim)=xcc_tmp(ix^d,1:ndim)*normconv(0)
2526  w_tec(1:nw+nwauxio)=wcc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2527  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2528  {end do\}
2529  else if (mype/=0) then
2530  itag=morton_no
2531  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2532  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision,0,itag,icomm,ierrmpi)
2533  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2534  call mpi_send(xcc_tmp,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
2535  end if
2536  enddo
2537  if (mype==0) then
2538  do ipe=1,npe-1
2539  do morton_no=morton_start(ipe),morton_stop(ipe)
2540  itag=morton_no
2541  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2542  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,&
2543  itag,icomm,intstatus(:,1),ierrmpi)
2544  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,&
2545  icomm,intstatus(:,1),ierrmpi)
2546  call mpi_recv(xcc_tmp_recv,1,type_block_xcc_io, ipe,itag,&
2547  icomm,intstatus(:,1),ierrmpi)
2548  {do ix^db=ixccmin^db,ixccmax^db\}
2549  x_tec(1:ndim)=xcc_tmp_recv(ix^d,1:ndim)*normconv(0)
2550  w_tec(1:nw+nwauxio)=wcc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2551  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2552  {end do\}
2553  end do
2554  end do
2555  close(qunit)
2556  end if
2557 else
2558 }
2559 
2560 if (mype/=0) then
2561  itag=1000*morton_stop(mype)
2562  call mpi_send(levmin,1,mpi_integer, 0,itag,icomm,ierrmpi)
2563  itag=2000*morton_stop(mype)
2564  call mpi_send(levmax,1,mpi_integer, 0,itag,icomm,ierrmpi)
2565 end if
2566 
2567 do level=levmin,levmax
2568  nodesonlevelmype=numgridsonlevel_mype(level,mype)*{nxc^d*}
2569  elemsonlevelmype=numgridsonlevel_mype(level,mype)*{nx^d*}
2570  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2571  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2572  ! for all tecplot variants coded up here, we let the TECPLOT ZONES coincide
2573  ! with the AMR grid LEVEL. Other options would be
2574  ! let each grid define a zone: inefficient for TECPLOT internal workings
2575  ! hence not implemented
2576  ! let entire octree define 1 zone: no difference in interpolation
2577  ! properties across TECPLOT zones detected as yet, hence not done
2578  select case(convert_type)
2579  case('tecplotmpi')
2580  ! in this option, we store the corner coordinates, as well as the corner
2581  ! values of all variables (obtained by averaging). This allows POINT packaging,
2582  ! and thus we can save full grid info by using one call to calc_grid
2583  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2584  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2585  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2586  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2587  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2588  do morton_no=morton_start(mype),morton_stop(mype)
2589  igrid = sfc_to_igrid(morton_no)
2590  if (mype/=0)then
2591  itag=morton_no
2592  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2593  itag=igrid
2594  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2595  end if
2596  if (node(plevel_,igrid)/=level) cycle
2597  call calc_x(igrid,xc,xcc)
2598  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2599  ixc^l,ixcc^l,.true.)
2600  if (mype/=0) then
2601  itag=morton_no
2602  ind_send=(/ ixc^l /)
2603  siz_ind=2*^nd
2604  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2605  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2606 
2607  call mpi_send(wc_tmp,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
2608  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2609  else
2610  {do ix^db=ixcmin^db,ixcmax^db\}
2611  x_tec(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0)
2612  w_tec(1:nw+nwauxio)=wc_tmp(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2613  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2614  {end do\}
2615  end if
2616  enddo
2617 
2618  case('tecplotCCmpi')
2619  ! in this option, we store the corner coordinates, and the cell center
2620  ! values of all variables. Due to this mix of corner/cell center, we must
2621  ! use BLOCK packaging, and thus we have enormous overhead by using
2622  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2623  ! and cell center values per dimension, per variable
2624  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2625  if(nw+nwauxio==1)then
2626  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2627  ! and just set [ndim+1]
2628  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2629  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2630  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2631  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2632  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2633  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2634  else
2635  if(ndim+nw+nwauxio<10) then
2636  ! difference only in length of integer format specification for ndim+nw+nwauxio
2637  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2638  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2639  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2640  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2641  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2642  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2643  else
2644  if (mype==0.and.(nodesonlevelmype>0.and.elemsonlevelmype>0))&
2645  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2646  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2647  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2648  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2649  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2650  endif
2651  endif
2652 
2653  do idim=1,ndim
2654  first=(idim==1)
2655  do morton_no=morton_start(mype),morton_stop(mype)
2656  igrid = sfc_to_igrid(morton_no)
2657  if (mype/=0)then
2658  itag=morton_no*idim
2659  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2660  itag=igrid*idim
2661  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2662  end if
2663  if (node(plevel_,igrid)/=level) cycle
2664 
2665  call calc_x(igrid,xc,xcc)
2666  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2667  ixc^l,ixcc^l,first)
2668  if (mype/=0)then
2669  ind_send=(/ ixc^l /)
2670  siz_ind=2*^nd
2671  itag=igrid*idim
2672  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2673  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2674  call mpi_send(xc_tmp,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
2675  else
2676  write(qunit,fmt="(100(e14.6))") xc_tmp(ixc^s,idim)*normconv(0)
2677  end if
2678  enddo
2679  enddo
2680 
2681  do iw=1,nw+nwauxio
2682  do morton_no=morton_start(mype),morton_stop(mype)
2683  igrid = sfc_to_igrid(morton_no)
2684  if (mype/=0)then
2685  itag=morton_no*(ndim+iw)
2686  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2687  itag=igrid*(ndim+iw)
2688  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2689  end if
2690  if (node(plevel_,igrid)/=level) cycle
2691 
2692  call calc_x(igrid,xc,xcc)
2693  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
2694  ixc^l,ixcc^l,.true.)
2695 
2696  if (mype/=0)then
2697  ind_send=(/ ixcc^l /)
2698  siz_ind=2*^nd
2699  itag=igrid*(ndim+iw)
2700  call mpi_send(ind_send,siz_ind, mpi_integer, 0,itag,icomm,ierrmpi)
2701  call mpi_send(normconv,nw+nwauxio+1,mpi_double_precision, 0,itag,icomm,ierrmpi)
2702  call mpi_send(wcc_tmp,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
2703  else
2704  write(qunit,fmt="(100(e14.6))") wcc_tmp(ixcc^s,iw)*normconv(iw)
2705  endif
2706  enddo
2707  enddo
2708  case default
2709  call mpistop('no such tecplot type')
2710  end select
2711 
2712 
2713  igonlevel=0
2714  do morton_no=morton_start(mype),morton_stop(mype)
2715  igrid = sfc_to_igrid(morton_no)
2716  if (mype/=0)then
2717  itag=morton_no
2718  call mpi_send(igrid,1,mpi_integer, 0,itag,icomm,ierrmpi)
2719  itag=igrid
2720  call mpi_send(node(plevel_,igrid),1,mpi_integer, 0,itag,icomm,ierrmpi)
2721  end if
2722  if (node(plevel_,igrid)/=level) cycle
2723 
2724  igonlevel=igonlevel+1
2725  if (mype/=0)then
2726  itag=igrid
2727  call mpi_send(igonlevel,1,mpi_integer, 0,itag,icomm,ierrmpi)
2728  end if
2729  if(mype==0)then
2730  call save_conntec(qunit,igrid,igonlevel)
2731  endif
2732  end do
2733 end do
2734 
2735 if (mype==0) then
2736  if (npe>1) then
2737  do ipe=1,npe-1
2738  itag=1000*morton_stop(ipe)
2739  call mpi_recv(levmin_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2740  itag=2000*morton_stop(ipe)
2741  call mpi_recv(levmax_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2742  do level=levmin_recv,levmax_recv
2743  nodesonlevelmype=numgridsonlevel_mype(level,ipe)*{nxc^d*}
2744  elemsonlevelmype=numgridsonlevel_mype(level,ipe)*{nx^d*}
2745  nodesonlevel=numgridsonlevel(level)*{nxc^d*}
2746  elemsonlevel=numgridsonlevel(level)*{nx^d*}
2747  select case(convert_type)
2748  case('tecplotmpi')
2749  ! in this option, we store the corner coordinates, as well as the corner
2750  ! values of all variables (obtained by averaging). This allows POINT packaging,
2751  ! and thus we can save full grid info by using one call to calc_grid
2752  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2753  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,a)") &
2754  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2755  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=POINT, ZONETYPE=', &
2756  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2757  do morton_no=morton_start(ipe),morton_stop(ipe)
2758  itag=morton_no
2759  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2760  itag=igrid_recv
2761  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2762  if (level_recv/=level) cycle
2763 
2764  itag=morton_no
2765  siz_ind=2*^nd
2766  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,&
2767  icomm,intstatus(:,1),ierrmpi)
2768  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2769  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2770  ,icomm,intstatus(:,1),ierrmpi)
2771 
2772  call mpi_recv(wc_tmp_recv,1,type_block_wc_io, ipe,itag,&
2773  icomm,intstatus(:,1),ierrmpi)
2774  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,&
2775  icomm,intstatus(:,1),ierrmpi)
2776 
2777  {do ix^db=ixrvcmin^db,ixrvcmax^db\}
2778  x_tec(1:ndim)=xc_tmp_recv(ix^d,1:ndim)*normconv(0)
2779  w_tec(1:nw+nwauxio)=wc_tmp_recv(ix^d,1:nw+nwauxio)*normconv(1:nw+nwauxio)
2780  write(qunit,fmt="(100(e14.6))") x_tec, w_tec
2781  {end do\}
2782  end do
2783  case('tecplotCCmpi')
2784  ! in this option, we store the corner coordinates, and the cell center
2785  ! values of all variables. Due to this mix of corner/cell center, we must
2786  ! use BLOCK packaging, and thus we have enormous overhead by using
2787  ! calc_grid repeatedly to merely fill values of cell corner coordinates
2788  ! and cell center values per dimension, per variable
2789  if(ndim+nw+nwauxio>99) call mpistop("adjust format specification in writeout")
2790  if(nw+nwauxio==1)then
2791  ! to make tecplot happy: avoid [ndim+1-ndim+1] in varlocation varset
2792  ! and just set [ndim+1]
2793  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2794  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,a)") &
2795  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2796  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2797  ndim+1,']=CELLCENTERED), ZONETYPE=', &
2798  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2799  else
2800  if(ndim+nw+nwauxio<10) then
2801  ! difference only in length of integer format specification for ndim+nw+nwauxio
2802  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2803  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i1,a,a)") &
2804  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2805  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2806  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2807  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2808  else
2809  if(nodesonlevelmype>0.and.elemsonlevelmype>0) &
2810  write(qunit,"(a,i7,a,a,i7,a,i7,a,f25.16,a,i1,a,i2,a,a)") &
2811  'ZONE T="',level,'"',', N=',nodesonlevelmype,', E=',elemsonlevelmype, &
2812  ', SOLUTIONTIME=',global_time*time_convert_factor,', DATAPACKING=BLOCK, VARLOCATION=([', &
2813  ndim+1,'-',ndim+nw+nwauxio,']=CELLCENTERED), ZONETYPE=', &
2814  {^ifoned 'FELINESEG'}{^iftwod 'FEQUADRILATERAL'}{^ifthreed 'FEBRICK'}
2815  endif
2816  endif
2817 
2818  do idim=1,ndim
2819  do morton_no=morton_start(ipe),morton_stop(ipe)
2820  itag=morton_no*idim
2821  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2822  itag=igrid_recv*idim
2823  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2824  if (level_recv/=level) cycle
2825 
2826  siz_ind=2*^nd
2827  itag=igrid_recv*idim
2828  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2829  ixrvcmin^d=ind_recv(^d);ixrvcmax^d=ind_recv(^nd+^d);
2830  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2831  ,icomm,intstatus(:,1),ierrmpi)
2832  call mpi_recv(xc_tmp_recv,1,type_block_xc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2833  write(qunit,fmt="(100(e14.6))") xc_tmp_recv(ixrvc^s,idim)*normconv(0)
2834  end do
2835  end do
2836 
2837  do iw=1,nw+nwauxio
2838  do morton_no=morton_start(ipe),morton_stop(ipe)
2839  itag=morton_no*(ndim+iw)
2840  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2841  itag=igrid_recv*(ndim+iw)
2842  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2843  if (level_recv/=level) cycle
2844 
2845  siz_ind=2*^nd
2846  itag=igrid_recv*(ndim+iw)
2847  call mpi_recv(ind_recv,siz_ind, mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2848  ixrvccmin^d=ind_recv(^d);ixrvccmax^d=ind_recv(^nd+^d);
2849  call mpi_recv(normconv,nw+nwauxio+1, mpi_double_precision,ipe,itag&
2850  ,icomm,intstatus(:,1),ierrmpi)
2851  call mpi_recv(wcc_tmp_recv,1,type_block_wcc_io, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2852  write(qunit,fmt="(100(e14.6))") wcc_tmp_recv(ixrvcc^s,iw)*normconv(iw)
2853  enddo
2854  enddo
2855  case default
2856  call mpistop('no such tecplot type')
2857  end select
2858 
2859  do morton_no=morton_start(ipe),morton_stop(ipe)
2860  itag=morton_no
2861  call mpi_recv(igrid_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2862  itag=igrid_recv
2863  call mpi_recv(level_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2864  if (level_recv/=level) cycle
2865 
2866  itag=igrid_recv
2867  call mpi_recv(igonlevel_recv,1,mpi_integer, ipe,itag,icomm,intstatus(:,1),ierrmpi)
2868  call save_conntec(qunit,igrid_recv,igonlevel_recv)
2869  end do ! morton loop
2870  end do ! level loop
2871  end do ! ipe loop
2872  end if ! npe>1 if
2873 end if ! mype=0 if
2874 {^ifoned endif}
2875 
2876 if (npe>1) then
2877  call mpi_barrier(icomm,ierrmpi)
2878  if(mype==0)deallocate(intstatus)
2879 endif
2880 
2881 end subroutine tecplot_mpi
2882 !=============================================================================
2883 subroutine punstructuredvtkb_mpi(qunit)
2884 
2885 ! Write one pvtu and vtu files for each processor
2886 ! Otherwise like unstructuredvtk_mpi
2887 ! output for vtu format to paraview, binary version output
2888 ! uses calc_grid to compute nwauxio variables
2889 ! allows renormalizing using convert factors
2890 
2893 use mod_calculate_xw
2894 
2895 integer, intent(in) :: qunit
2896 
2897 double precision :: x_VTK(1:3)
2898 
2899 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC_TMP
2900 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC_TMP
2901 double precision, dimension(ixMlo^D-1:ixMhi^D,ndim) :: xC
2902 double precision, dimension(ixMlo^D:ixMhi^D,ndim) :: xCC
2903 
2904 double precision, dimension(ixMlo^D-1:ixMhi^D,nw+nwauxio) :: wC_TMP
2905 double precision, dimension(ixMlo^D:ixMhi^D,nw+nwauxio) :: wCC_TMP
2906 
2907 integer :: igrid,iigrid,level,igonlevel,icel,ixC^L,ixCC^L,Morton_no
2908 integer :: NumGridsOnLevel(1:nlevelshi)
2909 integer :: nx^D,nxC^D,nodesonlevel,elemsonlevel,nc,np,VTK_type,ix^D
2910 double precision :: normconv(0:nw+nwauxio)
2911 character(len=80) :: pfilename
2912 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:ndim+nw+nwauxio)
2913 character(len=1024) :: outfilehead
2914 
2915 integer*8 :: offset
2916 integer:: recsep,k,iw,filenr
2917 integer:: length,lengthcc,offset_points,offset_cells, &
2918  length_coords,length_conn,length_offsets
2919 character:: buf
2920 character(len=6):: bufform
2921 
2922 logical :: fileopen
2923 !-----------------------------------------------------------------------------
2924 
2925 ! Write pvtu-file:
2926 if (mype==0) then
2927  call write_pvtu(qunit)
2928 endif
2929 ! Now write the Source files:
2930 inquire(qunit,opened=fileopen)
2931 if(.not.fileopen)then
2932  ! generate filename
2933  filenr=snapshotnext-1
2934  if (autoconvert) filenr=snapshotnext
2935  ! Open the file for the header part
2936  write(pfilename,'(a,i4.4,a,i4.4,a)') trim(base_filename),filenr,"p",mype,".vtu"
2937  open(qunit,file=pfilename,status='unknown',form='formatted')
2938 endif
2939 ! generate xml header
2940 write(qunit,'(a)')'<?xml version="1.0"?>'
2941 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
2942 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
2943 write(qunit,'(a)')' <UnstructuredGrid>'
2944 write(qunit,'(a)')'<FieldData>'
2945 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
2946  'NumberOfTuples="1" format="ascii">'
2947 write(qunit,*) real(global_time*time_convert_factor)
2948 write(qunit,'(a)')'</DataArray>'
2949 write(qunit,'(a)')'</FieldData>'
2950 
2951 
2952 offset=0
2953 recsep=4
2954 
2955 call getheadernames(wnamei,xandwnamei,outfilehead)
2956 
2957 ! number of cells, number of corner points, per grid.
2958 nx^d=ixmhi^d-ixmlo^d+1;
2959 nxc^d=nx^d+1;
2960 nc={nx^d*}
2961 np={nxc^d*}
2962 
2963 length=np*size_real
2964 lengthcc=nc*size_real
2965 
2966 length_coords=3*length
2967 length_conn=2**^nd*size_int*nc
2968 length_offsets=nc*size_int
2969 
2970 ! Note: using the w_write, writelevel, writespshift
2971 ! we can clip parts of the grid away, select variables, levels etc.
2972 do level=levmin,levmax
2973  if (writelevel(level)) then
2974  do morton_no=morton_start(mype),morton_stop(mype)
2975  igrid=sfc_to_igrid(morton_no)
2976  if (node(plevel_,igrid)/=level) cycle
2977  ! only output a grid when fully within clipped region selected
2978  ! by writespshift array
2979  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
2980  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
2981  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
2982  select case(convert_type)
2983  case('pvtuBmpi')
2984  ! we write out every grid as one VTK PIECE
2985  write(qunit,'(a,i7,a,i7,a)') &
2986  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
2987  write(qunit,'(a)')'<PointData>'
2988  do iw=1,nw
2989  if(.not.w_write(iw))cycle
2990 
2991  write(qunit,'(a,a,a,i16,a)')&
2992  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
2993  '" format="appended" offset="',offset,'">'
2994  write(qunit,'(a)')'</DataArray>'
2995  offset=offset+length+size_int
2996  enddo
2997  do iw=nw+1,nw+nwauxio
2998 
2999  write(qunit,'(a,a,a,i16,a)')&
3000  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
3001  '" format="appended" offset="',offset,'">'
3002  write(qunit,'(a)')'</DataArray>'
3003  offset=offset+length+size_int
3004  enddo
3005  write(qunit,'(a)')'</PointData>'
3006 
3007  write(qunit,'(a)')'<Points>'
3008  write(qunit,'(a,i16,a)') &
3009  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
3010  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3011  offset=offset+length_coords+size_int
3012  write(qunit,'(a)')'</Points>'
3013  case('pvtuBCCmpi')
3014  ! we write out every grid as one VTK PIECE
3015  write(qunit,'(a,i7,a,i7,a)') &
3016  '<Piece NumberOfPoints="',np,'" NumberOfCells="',nc,'">'
3017  write(qunit,'(a)')'<CellData>'
3018  do iw=1,nw
3019  if(.not.w_write(iw))cycle
3020 
3021  write(qunit,'(a,a,a,i16,a)')&
3022  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
3023  '" format="appended" offset="',offset,'">'
3024  write(qunit,'(a)')'</DataArray>'
3025  offset=offset+lengthcc+size_int
3026  enddo
3027  do iw=nw+1,nw+nwauxio
3028 
3029  write(qunit,'(a,a,a,i16,a)')&
3030  '<DataArray type="Float32" Name="',trim(wnamei(iw)), &
3031  '" format="appended" offset="',offset,'">'
3032  write(qunit,'(a)')'</DataArray>'
3033  offset=offset+lengthcc+size_int
3034  enddo
3035  write(qunit,'(a)')'</CellData>'
3036 
3037  write(qunit,'(a)')'<Points>'
3038  write(qunit,'(a,i16,a)') &
3039  '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',offset,'"/>'
3040  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3041  offset=offset+length_coords+size_int
3042  write(qunit,'(a)')'</Points>'
3043  end select
3044 
3045 
3046  write(qunit,'(a)')'<Cells>'
3047 
3048  ! connectivity part
3049  write(qunit,'(a,i16,a)')&
3050  '<DataArray type="Int32" Name="connectivity" format="appended" offset="',offset,'"/>'
3051  offset=offset+length_conn+size_int
3052 
3053  ! offsets data array
3054  write(qunit,'(a,i16,a)') &
3055  '<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,'"/>'
3056  offset=offset+length_offsets+size_int
3057 
3058  ! VTK cell type data array
3059  write(qunit,'(a,i16,a)') &
3060  '<DataArray type="Int32" Name="types" format="appended" offset="',offset,'"/>'
3061  offset=offset+size_int+nc*size_int
3062 
3063  write(qunit,'(a)')'</Cells>'
3064 
3065  write(qunit,'(a)')'</Piece>'
3066  endif
3067  enddo
3068  endif
3069 enddo
3070 
3071 write(qunit,'(a)')'</UnstructuredGrid>'
3072 write(qunit,'(a)')'<AppendedData encoding="raw">'
3073 
3074 close(qunit)
3075 ! next to make gfortran compiler happy, as it does not know
3076 ! form='binary' and produces error on compilation
3077 !bufform='binary'
3078 !open(qunit,file=pfilename,form=bufform,position='append')
3079 !This should in principle do also for gfortran (tested with gfortran 4.6.0 and Intel 11.1):
3080 open(qunit,file=pfilename,access='stream',form='unformatted',position='append')
3081 buf='_'
3082 write(qunit) trim(buf)
3083 
3084 do level=levmin,levmax
3085  if (writelevel(level)) then
3086  do morton_no=morton_start(mype),morton_stop(mype)
3087  igrid=sfc_to_igrid(morton_no)
3088  if (node(plevel_,igrid)/=level) cycle
3089  ! only output a grid when fully within clipped region selected
3090  ! by writespshift array
3091  if (({rnode(rpxmin^d_,igrid)>=xprobmin^d+(xprobmax^d-xprobmin^d)&
3092  *writespshift(^d,1)|.and.}).and.({rnode(rpxmax^d_,igrid)&
3093  <=xprobmax^d-(xprobmax^d-xprobmin^d)*writespshift(^d,2)|.and.})) then
3094  call calc_x(igrid,xc,xcc)
3095  call calc_grid(qunit,igrid,xc,xcc,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
3096  ixc^l,ixcc^l,.true.)
3097  do iw=1,nw
3098  if(.not.w_write(iw))cycle
3099  select case(convert_type)
3100  case('pvtuBmpi')
3101  write(qunit) length
3102  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
3103  case('pvtuBCCmpi')
3104  write(qunit) lengthcc
3105  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
3106  end select
3107  enddo
3108  do iw=nw+1,nw+nwauxio
3109  select case(convert_type)
3110  case('pvtuBmpi')
3111  write(qunit) length
3112  write(qunit) {(|}real(wc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixcmin^d,ixcmax^d)}
3113  case('pvtuBCCmpi')
3114  write(qunit) lengthcc
3115  write(qunit) {(|}real(wcc_tmp(ix^d,iw)*normconv(iw)),{ix^d=ixccmin^d,ixccmax^d)}
3116  end select
3117  enddo
3118 
3119  write(qunit) length_coords
3120  {do ix^db=ixcmin^db,ixcmax^db \}
3121  x_vtk(1:3)=zero;
3122  x_vtk(1:ndim)=xc_tmp(ix^d,1:ndim)*normconv(0);
3123  do k=1,3
3124  write(qunit) real(x_vtk(k))
3125  end do
3126  {end do \}
3127 
3128  write(qunit) length_conn
3129  {do ix^db=1,nx^db\}
3130  {^ifoned write(qunit)ix1-1,ix1 \}
3131  {^iftwod
3132  write(qunit)(ix2-1)*nxc1+ix1-1, &
3133  (ix2-1)*nxc1+ix1,ix2*nxc1+ix1-1,ix2*nxc1+ix1
3134  \}
3135  {^ifthreed
3136  write(qunit)&
3137  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3138  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3139  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3140  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3141  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3142  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3143  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3144  ix3*nxc2*nxc1+ ix2*nxc1+ix1
3145  \}
3146  {end do\}
3147 
3148  write(qunit) length_offsets
3149  do icel=1,nc
3150  write(qunit) icel*(2**^nd)
3151  end do
3152 
3153 
3154  {^ifoned vtk_type=3 \}
3155  {^iftwod vtk_type=8 \}
3156  {^ifthreed vtk_type=11 \}
3157  write(qunit) size_int*nc
3158  do icel=1,nc
3159  write(qunit) vtk_type
3160  enddo
3161  endif
3162  end do
3163  endif
3164 end do
3165 
3166 close(qunit)
3167 open(qunit,file=pfilename,status='unknown',form='formatted',position='append')
3168 
3169 write(qunit,'(a)')'</AppendedData>'
3170 write(qunit,'(a)')'</VTKFile>'
3171 close(qunit)
3172 
3173 end subroutine punstructuredvtkb_mpi
3174 {^iftwod
3175 ! subroutines to convert 2.5D data to 3D data
3176 !====================================================================================
3177 subroutine unstructuredvtkb23(qunit)
3178 
3179 ! output for vtu format to paraview, binary version output
3180 ! not parallel, uses calc_grid to compute nwauxio variables
3181 
3183 use mod_physics
3184 use mod_calculate_xw
3185 
3186 integer, intent(in) :: qunit
3187 
3188 double precision :: x_VTK(1:3)
3189 
3190 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
3191 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
3192 
3193 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
3194 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
3195 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
3196 integer:: igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
3197  ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
3198  ixCCmax2,ixCCmax3
3199 integer:: NumGridsOnLevel(1:nlevelshi)
3200 integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
3201  VTK_type,ix1,ix2,ix3
3202 double precision :: normconv(0:nw+nwauxio)
3203 character(len=80):: filename
3204 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
3205 character(len=1024) :: outfilehead
3206 integer*8 :: offset
3207 integer :: size_length,recsep,k,iw
3208 integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
3209  length_conn,length_offsets
3210 integer :: i3grid,n3grid
3211 double precision ::d3grid,zlengsc,zgridsc
3212 character:: buffer
3213 character(len=6):: bufform
3214 
3215 double precision :: zlength
3216 
3217 logical :: fileopen
3218 !-----------------------------------------------------------------------------
3219 if(npe>1)then
3220  if(mype==0) print *,'unstructuredvtkB23 not parallel, use vtumpi'
3221  call mpistop('npe>1, unstructuredvtkB23')
3222 end if
3223 
3224 offset=0
3225 recsep=4
3226 size_length=4
3227 
3228 inquire(qunit,opened=fileopen)
3229 if(.not.fileopen)then
3230  ! generate filename
3231  write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
3232  ! Open the file for the header part
3233  open(qunit,file=filename,status='replace')
3234 endif
3235 
3236 call getheadernames(wnamei,xandwnamei,outfilehead)
3237 ! generate xml header
3238 write(qunit,'(a)')'<?xml version="1.0"?>'
3239 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
3240 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3241 write(qunit,'(a)')'<UnstructuredGrid>'
3242 write(qunit,'(a)')'<FieldData>'
3243 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3244  'NumberOfTuples="1" format="ascii">'
3245 write(qunit,'(f10.2)') real(global_time*time_convert_factor)
3246 write(qunit,'(a)')'</DataArray>'
3247 write(qunit,'(a)')'</FieldData>'
3248 
3249 ! number of cells, number of corner points, per grid.
3250 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3251 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3252 nc=nx1*nx2*nx3
3253 np=nxc1*nxc2*nxc3
3254 
3255 length=np*size_real
3256 lengthcc=nc*size_real
3257 
3258 length_coords=3*length
3259 length_conn=2**3*size_int*nc
3260 length_offsets=nc*size_int
3261 
3262 ! Note: using the w_write, writelevel, writespshift
3263 ! we can clip parts of the grid away, select variables, levels etc.
3264 zgridsc=2.d0
3265 zlengsc=2.d0*zgridsc
3266 zlength=zlengsc*(xprobmax1-xprobmin1)
3267 do level=levmin,levmax
3268  if (writelevel(level)) then
3269  do iigrid=1,igridstail; igrid=igrids(iigrid);
3270  if (node(plevel_,igrid)/=level) cycle
3271  block=>ps(igrid)
3272  ! only output a grid when fully within clipped region selected
3273  ! by writespshift array
3274  if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3275  *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3276  +(xprobmax2-xprobmin2)*writespshift(2,1))&
3277  .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3278  *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3279  -xprobmin2)*writespshift(2,2))) then
3280  d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3281  n3grid=nint(zlength/d3grid)
3282  do i3grid=1,n3grid !subcycles
3283  select case(convert_type)
3284  case('vtuB23')
3285  ! we write out every grid as one VTK PIECE
3286  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3287  '" NumberOfCells="',nc,'">'
3288  write(qunit,'(a)')'<PointData>'
3289  do iw=1,nw
3290  if(.not.w_write(iw))cycle
3291  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3292  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3293  write(qunit,'(a)')'</DataArray>'
3294  offset=offset+length+size_int
3295  enddo
3296  if(nwauxio>0)then
3297  do iw=nw+1,nw+nwauxio
3298  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3299  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3300  write(qunit,'(a)')'</DataArray>'
3301  offset=offset+length+size_int
3302  enddo
3303  endif
3304  write(qunit,'(a)')'</PointData>'
3305 
3306  write(qunit,'(a)')'<Points>'
3307  write(qunit,'(a,i16,a)') &
3308 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3309  offset,'"/>'
3310  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3311  offset=offset+length_coords+size_int
3312  write(qunit,'(a)')'</Points>'
3313  case('vtuBCC23')
3314  ! we write out every grid as one VTK PIECE
3315  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3316  '" NumberOfCells="',nc,'">'
3317  write(qunit,'(a)')'<CellData>'
3318  do iw=1,nw
3319  if(.not.w_write(iw))cycle
3320  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3321  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3322  write(qunit,'(a)')'</DataArray>'
3323  offset=offset+lengthcc+size_int
3324  enddo
3325  if(nwauxio>0)then
3326  do iw=nw+1,nw+nwauxio
3327  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3328  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3329  write(qunit,'(a)')'</DataArray>'
3330  offset=offset+lengthcc+size_int
3331  enddo
3332  endif
3333  write(qunit,'(a)')'</CellData>'
3334 
3335  write(qunit,'(a)')'<Points>'
3336  write(qunit,'(a,i16,a)') &
3337 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3338  offset,'"/>'
3339  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3340  offset=offset+length_coords+size_int
3341  write(qunit,'(a)')'</Points>'
3342  end select
3343 
3344  write(qunit,'(a)')'<Cells>'
3345 
3346  ! connectivity part
3347  write(qunit,'(a,i16,a)')&
3348 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3349  offset,'"/>'
3350  offset=offset+length_conn+size_int
3351 
3352  ! offsets data array
3353  write(qunit,'(a,i16,a)') &
3354  '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3355  offset,'"/>'
3356  offset=offset+length_offsets+size_int
3357 
3358  ! VTK cell type data array
3359  write(qunit,'(a,i16,a)') &
3360  '<DataArray type="Int32" Name="types" format="appended" offset="',&
3361  offset,'"/>'
3362  offset=offset+size_length+nc*size_int
3363 
3364  write(qunit,'(a)')'</Cells>'
3365 
3366  write(qunit,'(a)')'</Piece>'
3367  enddo !subcycles
3368  endif
3369  enddo
3370  endif
3371 enddo
3372 
3373 write(qunit,'(a)')'</UnstructuredGrid>'
3374 write(qunit,'(a)')'<AppendedData encoding="raw">'
3375 
3376 close(qunit)
3377 open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
3378 buffer='_'
3379 write(qunit) trim(buffer)
3380 
3381 do level=levmin,levmax
3382  if (writelevel(level)) then
3383  do iigrid=1,igridstail; igrid=igrids(iigrid);
3384  if (node(plevel_,igrid)/=level) cycle
3385  block=>ps(igrid)
3386  ! only output a grid when fully within clipped region selected
3387  ! by writespshift array
3388  if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3389  *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3390  +(xprobmax2-xprobmin2)*writespshift(2,1))&
3391  .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3392  *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3393  -xprobmin2)*writespshift(2,2))) then
3394  d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3395  n3grid=nint(zlength/d3grid)
3396  ! In case primitives to be saved: use primitive subroutine
3397  ! extra layer around mesh only needed when storing corner values and averaging
3398  if(saveprim) then
3399  call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
3400  ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
3401  endif
3402  ! using array w so that new output auxiliaries can be calculated by the user
3403  ! extend 2D data to 3D insuring variables are independent on the third coordinate
3404  do ix3=ixglo1,ixghi1
3405  w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
3406  ixglo2:ixghi2,1:nw)
3407  end do
3408  do i3grid=1,n3grid !subcycles
3409  call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
3410  ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
3411  ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
3412  do iw=1,nw
3413  if(.not.w_write(iw))cycle
3414  select case(convert_type)
3415  case('vtuB23')
3416  write(qunit) length
3417  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3418  =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3419  case('vtuBCC23')
3420  write(qunit) lengthcc
3421  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3422  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3423  =ixccmin3,ixccmax3)
3424  end select
3425  enddo
3426  if(nwauxio>0)then
3427  do iw=nw+1,nw+nwauxio
3428  select case(convert_type)
3429  case('vtuB23')
3430  write(qunit) length
3431  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3432  =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3433  case('vtuBCC23')
3434  write(qunit) lengthcc
3435  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3436  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3437  =ixccmin3,ixccmax3)
3438  end select
3439  end do
3440  end if
3441  write(qunit) length_coords
3442  do ix3=ixcmin3,ixcmax3
3443  do ix2=ixcmin2,ixcmax2
3444  do ix1=ixcmin1,ixcmax1
3445  x_vtk(1:3)=zero;
3446  x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3447  do k=1,3
3448  write(qunit) real(x_vtk(k))
3449  end do
3450  end do
3451  end do
3452  end do
3453 
3454  write(qunit) length_conn
3455  do ix3=1,nx3
3456  do ix2=1,nx2
3457  do ix1=1,nx1
3458  write(qunit)&
3459  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3460  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3461  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3462  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3463  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3464  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3465  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3466  ix3*nxc2*nxc1+ ix2*nxc1+ix1
3467 
3468  end do
3469  end do
3470  end do
3471  write(qunit) length_offsets
3472  do icel=1,nc
3473  write(qunit) icel*(2**3)
3474  end do
3475 
3476  vtk_type=11
3477  write(qunit) size_int*nc
3478  do icel=1,nc
3479  write(qunit) vtk_type
3480  end do
3481  end do !subcycles
3482  end if
3483  end do
3484  end if
3485 end do
3486 
3487 close(qunit)
3488 open(qunit,file=filename,status='unknown',form='formatted',position='append')
3489 
3490 write(qunit,'(a)')'</AppendedData>'
3491 write(qunit,'(a)')'</VTKFile>'
3492 close(qunit)
3493 end subroutine unstructuredvtkb23
3494 !====================================================================================
3495 subroutine unstructuredvtkbsym23(qunit)
3496 
3497 ! output for vtu format to paraview, binary version output
3498 ! not parallel, uses calc_grid to compute nwauxio variables
3499 ! use this subroutine when the physical domain is symmetric/asymmetric about (0,y,z)
3500 ! plane, xprobmin1=0 and the computational domain is a half of the physical domain
3501 
3503 use mod_physics
3504 use mod_calculate_xw
3505 
3506 integer, intent(in) :: qunit
3507 
3508 double precision :: x_VTK(1:3)
3509 
3510 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
3511 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
3512 
3513 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
3514 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
3515 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
3516 integer:: igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
3517  ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
3518  ixCCmax2,ixCCmax3
3519 integer:: NumGridsOnLevel(1:nlevelshi)
3520 integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
3521  VTK_type,ix1,ix2,ix3
3522 double precision :: normconv(0:nw+nwauxio)
3523 character(len=80):: filename
3524 character(len=name_len) :: wnamei(1:nw+nwauxio),xandwnamei(1:3+nw+nwauxio)
3525 character(len=1024) :: outfilehead
3526 integer*8 :: offset
3527 integer :: size_length,recsep,k,iw
3528 integer :: length,lengthcc,offset_points,offset_cells, length_coords,&
3529  length_conn,length_offsets
3530 integer :: i3grid,n3grid
3531 double precision ::d3grid,zlengsc,zgridsc
3532 character:: buffer
3533 character(len=6):: bufform
3534 
3535 double precision :: zlength
3536 
3537 logical :: fileopen
3538 !-----------------------------------------------------------------------------
3539 if(npe>1)then
3540  if(mype==0) print *,'unstructuredvtkBsym23 not parallel, use vtumpi'
3541  call mpistop('npe>1, unstructuredvtkBsym23')
3542 end if
3543 
3544 offset=0
3545 recsep=4
3546 size_length=4
3547 
3548 inquire(qunit,opened=fileopen)
3549 if(.not.fileopen)then
3550  ! generate filename
3551  write(filename,'(a,a,i4.4,a)') trim(base_filename),"3D",snapshotini,".vtu"
3552  ! Open the file for the header part
3553  open(qunit,file=filename,status='unknown')
3554 endif
3555 
3556 call getheadernames(wnamei,xandwnamei,outfilehead)
3557 ! generate xml header
3558 write(qunit,'(a)')'<?xml version="1.0"?>'
3559 write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
3560 write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
3561 write(qunit,'(a)')'<UnstructuredGrid>'
3562 write(qunit,'(a)')'<FieldData>'
3563 write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
3564  'NumberOfTuples="1" format="ascii">'
3565 write(qunit,'(f10.2)') real(global_time*time_convert_factor)
3566 write(qunit,'(a)')'</DataArray>'
3567 write(qunit,'(a)')'</FieldData>'
3568 
3569 ! number of cells, number of corner points, per grid.
3570 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
3571 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
3572 nc=nx1*nx2*nx3
3573 np=nxc1*nxc2*nxc3
3574 
3575 length=np*size_real
3576 lengthcc=nc*size_real
3577 
3578 length_coords=3*length
3579 length_conn=2**3*size_int*nc
3580 length_offsets=nc*size_int
3581 
3582 ! Note: using the w_write, writelevel, writespshift
3583 ! we can clip parts of the grid away, select variables, levels etc.
3584 zlengsc=4.d0
3585 zgridsc=2.d0
3586 zlength=zlengsc*(xprobmax1-xprobmin1)
3587 do level=levmin,levmax
3588  if (writelevel(level)) then
3589  do iigrid=1,igridstail; igrid=igrids(iigrid);
3590  if (node(plevel_,igrid)/=level) cycle
3591  block=>ps(igrid)
3592  ! only output a grid when fully within clipped region selected
3593  ! by writespshift array
3594  if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3595  *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3596  +(xprobmax2-xprobmin2)*writespshift(2,1))&
3597  .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3598  *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3599  -xprobmin2)*writespshift(2,2))) then
3600  d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3601  n3grid=nint(zlength/d3grid)
3602  do i3grid=1,n3grid !subcycles
3603  !! original domain ----------------------------------start
3604  select case(convert_type)
3605  case('vtuBsym23')
3606  ! we write out every grid as one VTK PIECE
3607  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3608  '" NumberOfCells="',nc,'">'
3609  write(qunit,'(a)')'<PointData>'
3610  do iw=1,nw
3611  if(.not.w_write(iw))cycle
3612  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3613  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3614  write(qunit,'(a)')'</DataArray>'
3615  offset=offset+length+size_length
3616  enddo
3617  if(nwauxio>0)then
3618  do iw=nw+1,nw+nwauxio
3619  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3620  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3621  write(qunit,'(a)')'</DataArray>'
3622  offset=offset+length+size_length
3623  enddo
3624  endif
3625  write(qunit,'(a)')'</PointData>'
3626 
3627  write(qunit,'(a)')'<Points>'
3628  write(qunit,'(a,i16,a)') &
3629 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3630  offset,'"/>'
3631  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3632  offset=offset+length_coords+size_length
3633  write(qunit,'(a)')'</Points>'
3634  case('vtuBCCsym23')
3635  ! we write out every grid as one VTK PIECE
3636  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3637  '" NumberOfCells="',nc,'">'
3638  write(qunit,'(a)')'<CellData>'
3639  do iw=1,nw
3640  if(.not.w_write(iw))cycle
3641  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3642  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3643  write(qunit,'(a)')'</DataArray>'
3644  offset=offset+lengthcc+size_length
3645  enddo
3646  if(nwauxio>0)then
3647  do iw=nw+1,nw+nwauxio
3648  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3649  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3650  write(qunit,'(a)')'</DataArray>'
3651  offset=offset+lengthcc+size_length
3652  enddo
3653  endif
3654  write(qunit,'(a)')'</CellData>'
3655 
3656  write(qunit,'(a)')'<Points>'
3657  write(qunit,'(a,i16,a)') &
3658 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3659  offset,'"/>'
3660  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3661  offset=offset+length_coords+size_length
3662  write(qunit,'(a)')'</Points>'
3663  end select
3664 
3665 
3666  write(qunit,'(a)')'<Cells>'
3667 
3668  ! connectivity part
3669  write(qunit,'(a,i16,a)')&
3670 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3671  offset,'"/>'
3672  offset=offset+length_conn+size_length
3673 
3674  ! offsets data array
3675  write(qunit,'(a,i16,a)') &
3676  '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3677  offset,'"/>'
3678  offset=offset+length_offsets+size_length
3679 
3680  ! VTK cell type data array
3681  write(qunit,'(a,i16,a)') &
3682  '<DataArray type="Int32" Name="types" format="appended" offset="',&
3683  offset,'"/>'
3684  offset=offset+size_length+nc*size_int
3685 
3686  write(qunit,'(a)')'</Cells>'
3687 
3688  write(qunit,'(a)')'</Piece>'
3689  !! original domain ----------------------------------end
3690  !! symetric/asymetric mirror domain -----------------start
3691  select case(convert_type)
3692  case('vtuBsym23')
3693  ! we write out every grid as one VTK PIECE
3694  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3695  '" NumberOfCells="',nc,'">'
3696  write(qunit,'(a)')'<PointData>'
3697  do iw=1,nw
3698  if(.not.w_write(iw))cycle
3699  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3700  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3701  write(qunit,'(a)')'</DataArray>'
3702  offset=offset+length+size_length
3703  enddo
3704  if(nwauxio>0)then
3705  do iw=nw+1,nw+nwauxio
3706  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3707  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3708  write(qunit,'(a)')'</DataArray>'
3709  offset=offset+length+size_length
3710  enddo
3711  endif
3712  write(qunit,'(a)')'</PointData>'
3713 
3714  write(qunit,'(a)')'<Points>'
3715  write(qunit,'(a,i16,a)') &
3716 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3717  offset,'"/>'
3718  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3719  offset=offset+length_coords+size_length
3720  write(qunit,'(a)')'</Points>'
3721  case('vtuBCCsym23')
3722  ! we write out every grid as one VTK PIECE
3723  write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
3724  '" NumberOfCells="',nc,'">'
3725  write(qunit,'(a)')'<CellData>'
3726  do iw=1,nw
3727  if(.not.w_write(iw))cycle
3728  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3729  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3730  write(qunit,'(a)')'</DataArray>'
3731  offset=offset+lengthcc+size_length
3732  enddo
3733  if(nwauxio>0)then
3734  do iw=nw+1,nw+nwauxio
3735  write(qunit,'(a,a,a,i16,a)')'<DataArray type="Float32" Name="',&
3736  trim(wnamei(iw)), '" format="appended" offset="',offset,'">'
3737  write(qunit,'(a)')'</DataArray>'
3738  offset=offset+lengthcc+size_length
3739  enddo
3740  endif
3741  write(qunit,'(a)')'</CellData>'
3742 
3743  write(qunit,'(a)')'<Points>'
3744  write(qunit,'(a,i16,a)') &
3745 '<DataArray type="Float32" NumberOfComponents="3" format="appended" offset="',&
3746  offset,'"/>'
3747  ! write cell corner coordinates in a backward dimensional loop, always 3D output
3748  offset=offset+length_coords+size_length
3749  write(qunit,'(a)')'</Points>'
3750  end select
3751 
3752 
3753  write(qunit,'(a)')'<Cells>'
3754 
3755  ! connectivity part
3756  write(qunit,'(a,i16,a)')&
3757 '<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
3758  offset,'"/>'
3759  offset=offset+length_conn+size_length
3760 
3761  ! offsets data array
3762  write(qunit,'(a,i16,a)') &
3763  '<DataArray type="Int32" Name="offsets" format="appended" offset="',&
3764  offset,'"/>'
3765  offset=offset+length_offsets+size_length
3766 
3767  ! VTK cell type data array
3768  write(qunit,'(a,i16,a)') &
3769  '<DataArray type="Int32" Name="types" format="appended" offset="',&
3770  offset,'"/>'
3771  offset=offset+size_length+nc*size_int
3772 
3773  write(qunit,'(a)')'</Cells>'
3774 
3775  write(qunit,'(a)')'</Piece>'
3776  !! symetric/asymetric mirror domain -----------------end
3777  enddo !subcycles
3778  endif
3779  enddo
3780  endif
3781 enddo
3782 
3783 write(qunit,'(a)')'</UnstructuredGrid>'
3784 write(qunit,'(a)')'<AppendedData encoding="raw">'
3785 
3786 close(qunit)
3787 open(qunit,file=filename,form='unformatted',access='stream',status='old',position='append')
3788 buffer='_'
3789 write(qunit) trim(buffer)
3790 do level=levmin,levmax
3791  if (writelevel(level)) then
3792  do iigrid=1,igridstail; igrid=igrids(iigrid);
3793  if (node(plevel_,igrid)/=level) cycle
3794  block=>ps(igrid)
3795  ! only output a grid when fully within clipped region selected
3796  ! by writespshift array
3797  if ((rnode(rpxmin1_,igrid)>=xprobmin1+(xprobmax1-xprobmin1)&
3798  *writespshift(1,1).and.rnode(rpxmin2_,igrid)>=xprobmin2&
3799  +(xprobmax2-xprobmin2)*writespshift(2,1))&
3800  .and.(rnode(rpxmax1_,igrid)<=xprobmax1-(xprobmax1-xprobmin1)&
3801  *writespshift(1,2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2&
3802  -xprobmin2)*writespshift(2,2))) then
3803  d3grid=zgridsc*(rnode(rpxmax1_,igrid)-rnode(rpxmin1_,igrid))
3804  n3grid=nint(zlength/d3grid)
3805  ! In case primitives to be saved: use primitive subroutine
3806  ! extra layer around mesh only needed when storing corner values and averaging
3807  if(saveprim) then
3808  call phys_to_primitive(ixglo1,ixglo2,ixghi1,ixghi2,&
3809  ixglo1,ixglo2,ixghi1,ixghi2,ps(igrid)%w,ps(igrid)%x)
3810  endif
3811  ! using array w so that new output auxiliaries can be calculated by the user
3812  ! extend 2D data to 3D insuring variables are independent on the third coordinate
3813  do ix3=ixglo1,ixghi1
3814  w(ixglo1:ixghi1,ixglo2:ixghi2,ix3,1:nw)=ps(igrid)%w(ixglo1:ixghi1,&
3815  ixglo2:ixghi2,1:nw)
3816  end do
3817  do i3grid=1,n3grid !subcycles
3818  call calc_grid23(qunit,igrid,xc_tmp,xcc_tmp,wc_tmp,wcc_tmp,normconv,&
3819  ixcmin1,ixcmin2,ixcmin3,ixcmax1,ixcmax2,ixcmax3,ixccmin1,ixccmin2,&
3820  ixccmin3,ixccmax1,ixccmax2,ixccmax3,.true.,i3grid,d3grid,w,zlength,zgridsc)
3821  !! original domain ----------------------------------start
3822  do iw=1,nw
3823  if(.not.w_write(iw))cycle
3824  select case(convert_type)
3825  case('vtuBsym23')
3826  write(qunit) length
3827  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3828  =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3829  case('vtuBCCsym23')
3830  write(qunit) lengthcc
3831  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3832  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3833  =ixccmin3,ixccmax3)
3834  end select
3835  enddo
3836  if(nwauxio>0)then
3837  do iw=nw+1,nw+nwauxio
3838  select case(convert_type)
3839  case('vtuBsym23')
3840  write(qunit) length
3841  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3842  =ixcmin1,ixcmax1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3843  case('vtuBCCsym23')
3844  write(qunit) lengthcc
3845  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3846  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3847  =ixccmin3,ixccmax3)
3848  end select
3849  enddo
3850  endif
3851  write(qunit) length_coords
3852  do ix3=ixcmin3,ixcmax3
3853  do ix2=ixcmin2,ixcmax2
3854  do ix1=ixcmin1,ixcmax1
3855  x_vtk(1:3)=zero;
3856  x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3857  do k=1,3
3858  write(qunit) real(x_vtk(k))
3859  end do
3860  end do
3861  end do
3862  end do
3863 
3864  write(qunit) length_conn
3865  do ix3=1,nx3
3866  do ix2=1,nx2
3867  do ix1=1,nx1
3868  write(qunit)&
3869  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3870  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3871  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3872  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3873  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3874  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3875  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3876  ix3*nxc2*nxc1+ ix2*nxc1+ix1
3877  end do
3878  end do
3879  end do
3880  write(qunit) length_offsets
3881  do icel=1,nc
3882  write(qunit) icel*(2**3)
3883  end do
3884 
3885  vtk_type=11
3886  write(qunit) size_int*nc
3887  do icel=1,nc
3888  write(qunit) vtk_type
3889  end do
3890  !! original domain ----------------------------------end
3891  !! symetric/asymetric mirror domain -----------------start
3892  do iw=1,nw
3893  if(.not.w_write(iw))cycle
3894  if(iw==2 .or. iw==4 .or. iw==7) then
3895  wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)=&
3896  -wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,iw)
3897  wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)=&
3898  -wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,iw)
3899  end if
3900  select case(convert_type)
3901  case('vtuBsym23')
3902  write(qunit) length
3903  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3904  =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3905  case('vtuBCCsym23')
3906  write(qunit) lengthcc
3907  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3908  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3909  =ixccmin3,ixccmax3)
3910  end select
3911  enddo
3912  if(nwauxio>0)then
3913  do iw=nw+1,nw+nwauxio
3914  select case(convert_type)
3915  case('vtuBsym23')
3916  write(qunit) length
3917  write(qunit) (((real(wc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3918  =ixcmax1,ixcmin1,-1),ix2=ixcmin2,ixcmax2),ix3=ixcmin3,ixcmax3)
3919  case('vtuBCCsym23')
3920  write(qunit) lengthcc
3921  write(qunit) (((real(wcc_tmp(ix1,ix2,ix3,iw)*normconv(iw)),ix1&
3922  =ixccmin1,ixccmax1),ix2=ixccmin2,ixccmax2),ix3&
3923  =ixccmin3,ixccmax3)
3924  end select
3925  end do
3926  end if
3927  write(qunit) length_coords
3928  do ix3=ixcmin3,ixcmax3
3929  do ix2=ixcmin2,ixcmax2
3930  do ix1=ixcmax1,ixcmin1,-1
3931  x_vtk(1:3)=zero;
3932  x_vtk(1:3)=xc_tmp(ix1,ix2,ix3,1:3)*normconv(0);
3933  x_vtk(1)=-x_vtk(1)
3934  do k=1,3
3935  write(qunit) real(x_vtk(k))
3936  end do
3937  end do
3938  end do
3939  end do
3940 
3941  write(qunit) length_conn
3942  do ix3=1,nx3
3943  do ix2=1,nx2
3944  do ix1=nx1,1,-1
3945  write(qunit)&
3946  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
3947  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3948  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3949  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
3950  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
3951  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
3952  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
3953  ix3*nxc2*nxc1+ ix2*nxc1+ix1
3954  end do
3955  end do
3956  end do
3957  write(qunit) length_offsets
3958  do icel=1,nc
3959  write(qunit) icel*(2**3)
3960  end do
3961 
3962  vtk_type=11
3963  write(qunit) size_int*nc
3964  do icel=1,nc
3965  write(qunit) vtk_type
3966  end do
3967  !! symetric/asymetric mirror domain -----------------end
3968  end do !subcycles
3969  end if
3970  end do
3971  end if
3972 end do
3973 
3974 close(qunit)
3975 open(qunit,file=filename,status='unknown',form='formatted',position='append')
3976 
3977 write(qunit,'(a)')'</AppendedData>'
3978 write(qunit,'(a)')'</VTKFile>'
3979 close(qunit)
3980 end subroutine unstructuredvtkbsym23
3981 !=============================================================================
3982 subroutine calc_grid23(qunit,igrid,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,normconv,&
3983 ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,&
3984  ixCCmax1,ixCCmax2,ixCCmax3,first,i3grid,d3grid,w,zlength,zgridsc)
3985 
3986 ! this subroutine computes both corner as well as cell-centered values
3987 ! it handles how we do the center to corner averaging, as well as
3988 ! whether we switch to cartesian or want primitive or conservative output,
3989 ! handling the addition of B0 in B0+B1 cases, ...
3990 !
3991 ! the normconv is passed on to specialvar_output for extending with
3992 ! possible normalization values for the nw+1:nw+nwauxio entries
3993 
3995 integer, intent(in) :: qunit, igrid,i3grid
3996 logical, intent(in) :: first
3997 
3998 integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3, ix, iw, level, idir
3999 integer :: ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
4000  ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,nxCC1,nxCC2,nxCC3
4001 double precision :: dx1,dx2,dx3,d3grid,zlength,zgridsc
4002 
4003 integer :: idims,jxCmin1,jxCmin2,jxCmin3,jxCmax1,jxCmax2,jxCmax3
4004 double precision :: ldw(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1),&
4005  dwC(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1)
4006 
4007 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC
4008 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC
4009 
4010 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC
4011 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC
4012 
4013 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,3) :: xC_TMP
4014 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,& 3) :: xCC_TMP
4015 
4016 double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,ixMlo1& -1:ixMhi1,nw+nwauxio) :: wC_TMP
4017 double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo1:ixMhi1,nw& +nwauxio) :: wCC_TMP
4018 double precision, dimension(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:nw& +nwauxio) :: w
4019 
4020 double precision,dimension(0:nw+nwauxio) :: normconv
4021 logical, save :: subfirst=.true.
4022 !-----------------------------------------------------------------------------
4023 ! following only for allowing compiler to go through with debug on
4024 
4025 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
4026 level=node(plevel_,igrid)
4027 dx1=dx(1,level);dx2=dx(2,level);dx3=zgridsc*dx(1,level);
4028 
4029 ! for normalization within the code
4030 if(saveprim) then
4031  normconv(0) = length_convert_factor
4032  normconv(1:nw) = w_convert_factor
4033 else
4034  normconv(0)=length_convert_factor
4035  ! assuming density
4036  normconv(1)=w_convert_factor(1)
4037  ! assuming momentum=density*velocity
4038  if (nw>=2) normconv(2:2+3)=w_convert_factor(1)*w_convert_factor(2:2+3)
4039  ! assuming energy/pressure and magnetic field
4040  if (nw>=2+3) normconv(2+3:nw)=w_convert_factor(2+3:nw)
4041 end if
4042 
4043 ! coordinates of cell centers
4044 nxcc1=nx1;nxcc2=nx2;nxcc3=nx3;
4045 ixccmin1=ixmlo1;ixccmin2=ixmlo2;ixccmin3=ixmlo1; ixccmax1=ixmhi1
4046 ixccmax2=ixmhi2;ixccmax3=ixmhi1;
4047 do ix=ixccmin1,ixccmax1
4048  xcc(ix,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1)=rnode(rpxmin1_,igrid)&
4049  +(dble(ix-ixccmin1)+half)*dx1
4050 end do
4051 do ix=ixccmin2,ixccmax2
4052  xcc(ixccmin1:ixccmax1,ix,ixccmin3:ixccmax3,2)=rnode(rpxmin2_,igrid)&
4053  +(dble(ix-ixccmin2)+half)*dx2
4054 end do
4055 do ix=ixccmin3,ixccmax3
4056  xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ix,3)=-zlength/two+&
4057  dble(i3grid-1)*d3grid+(dble(ix-ixccmin3)+half)*dx3
4058 end do
4059 
4060 ! coordinates of cell corners
4061 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
4062 ixcmin1=ixmlo1-1;ixcmin2=ixmlo2-1;ixcmin3=ixmlo1-1; ixcmax1=ixmhi1
4063 ixcmax2=ixmhi2;ixcmax3=ixmhi1;
4064 do ix=ixcmin1,ixcmax1
4065  xc(ix,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1)=rnode(rpxmin1_,igrid)&
4066  +dble(ix-ixcmin1)*dx1
4067 end do
4068 do ix=ixcmin2,ixcmax2
4069  xc(ixcmin1:ixcmax1,ix,ixcmin3:ixcmax3,2)=rnode(rpxmin2_,igrid)&
4070  +dble(ix-ixcmin2)*dx2
4071 end do
4072 do ix=ixcmin3,ixcmax3
4073  xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ix,3)=-zlength/two+&
4074  dble(i3grid-1)*d3grid+dble(ix-ixcmin3)*dx3
4075 end do
4076 
4077 
4078 if (nwextra>0) then
4079  ! here we actually fill the ghost layers for the nwextra variables using
4080  ! continuous extrapolation (as these values do not exist normally in ghost
4081  ! cells)
4082  do idims=1,3
4083  select case(idims)
4084  case(1)
4085  jxcmin1=ixghi1+1-nghostcells;jxcmin2=ixglo2;jxcmin3=ixglo1;
4086  jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
4087  do ix1=jxcmin1,jxcmax1
4088  w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmin1&
4089  -1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
4090  end do
4091  jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
4092  jxcmax1=ixglo1-1+nghostcells;jxcmax2=ixghi2;jxcmax3=ixghi1;
4093  do ix1=jxcmin1,jxcmax1
4094  w(ix1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw) = w(jxcmax1&
4095  +1,jxcmin2:jxcmax2,jxcmin3:jxcmax3,nw-nwextra+1:nw)
4096  end do
4097  case(2)
4098  jxcmin1=ixglo1;jxcmin2=ixghi2+1-nghostcells;jxcmin3=ixglo1;
4099  jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
4100  do ix2=jxcmin2,jxcmax2
4101  w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
4102  = w(jxcmin1:jxcmax1,jxcmin2-1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
4103  end do
4104  jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
4105  jxcmax1=ixghi1;jxcmax2=ixglo2-1+nghostcells;jxcmax3=ixghi1;
4106  do ix2=jxcmin2,jxcmax2
4107  w(jxcmin1:jxcmax1,ix2,jxcmin3:jxcmax3,nw-nwextra+1:nw) &
4108  = w(jxcmin1:jxcmax1,jxcmax2+1,jxcmin3:jxcmax3,nw-nwextra+1:nw)
4109  end do
4110  case(3)
4111  jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixghi1+1-nghostcells;
4112  jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixghi1;
4113  do ix3=jxcmin3,jxcmax3
4114  w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
4115  = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmin3-1,nw-nwextra+1:nw)
4116  end do
4117  jxcmin1=ixglo1;jxcmin2=ixglo2;jxcmin3=ixglo1;
4118  jxcmax1=ixghi1;jxcmax2=ixghi2;jxcmax3=ixglo1-1+nghostcells;
4119  do ix3=jxcmin3,jxcmax3
4120  w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,ix3,nw-nwextra+1:nw) &
4121  = w(jxcmin1:jxcmax1,jxcmin2:jxcmax2,jxcmax3+1,nw-nwextra+1:nw)
4122  end do
4123  end select
4124  end do
4125 end if
4126 ! next lines needed when specialvar_output uses gradients
4127 ! and later on when dwlimiter2 is used
4128 if(nwauxio>0)then
4129  ! auxiliary io variables can be computed and added by user
4130  ! next few lines ensure correct usage of routines like divvector etc
4131  dxlevel(1)=rnode(rpdx1_,igrid);dxlevel(2)=rnode(rpdx2_,igrid)
4132  ! default (no) normalization for auxiliary variables
4133  normconv(nw+1:nw+nwauxio)=one
4134  ! maybe need for restriction to ixG^LL^LSUB1
4135  call specialvar_output23(ixglo1,ixglo2,ixglo1,ixghi1,ixghi2,ixghi1,ixglo1&
4136  +1,ixglo2+1,ixglo1+1,ixghi1-1,ixghi2-1,ixghi1-1,w,xcc,normconv)
4137 endif
4138 
4139 ! compute the cell-center values for w first
4140 !===========================================
4141 ! cell center values obtained from mere copy, while B0+B1 split handled here
4142 
4143 wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)=w(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,:)
4144 if(b0field) then
4145  do ix3=ixccmin3,ixccmax3
4146  do ix2=ixccmin2,ixccmax2
4147  do ix1=ixccmin1,ixccmax1
4148  wcc(ix1,ix2,ix3,iw_mag(:))=wcc(ix1,ix2,ix3,iw_mag(:))+ps(igrid)%B0(ix1,ix2,&
4149  :,0)
4150  end do
4151  end do
4152  end do
4153 end if
4154 
4155 if(.not.saveprim .and. b0field .and. iw_e>0) then
4156  do ix3=ixccmin3,ixccmax3
4157  do ix2=ixccmin2,ixccmax2
4158  do ix1=ixccmin1,ixccmax1
4159  wcc(ix1,ix2,ix3,iw_e)=w(ix1,ix2,ix3,iw_e) +half*sum(ps(igrid)%B0(ix1,&
4160  ix2,:,0)**2 ) + sum(w(ix1,ix2,ix3,&
4161  iw_mag(:))*ps(igrid)%B0(ix1,ix2,:,0))
4162  end do
4163  end do
4164  end do
4165 end if
4166 ! compute the corner values for w now by averaging
4167 !=================================================
4168 if(slab_uniform)then
4169  ! for slab symmetry: no geometrical info required
4170  do iw=1,nw+nwauxio
4171  if (b0field.and.iw>iw_mag(1)-1.and.iw<=iw_mag(ndir)) then
4172  idir=iw-iw_mag(1)+1
4173  do ix3=ixcmin3,ixcmax3
4174  do ix2=ixcmin2,ixcmax2
4175  do ix1=ixcmin1,ixcmax1
4176  wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3,iw) &
4177  +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4178  ,idir,0))/dble(2**3)+&
4179  sum(w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw) &
4180  +ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4181  ,idir,0))/dble(2**3)
4182  end do
4183  end do
4184  end do
4185  else
4186  do ix3=ixcmin3,ixcmax3
4187  do ix2=ixcmin2,ixcmax2
4188  do ix1=ixcmin1,ixcmax1
4189  wc(ix1,ix2,ix3,iw)=sum(w(ix1:ix1+1,ix2:ix2+1,ix3:ix3&
4190  +1,iw))/dble(2**3)
4191  end do
4192  end do
4193  end do
4194  end if
4195  end do
4196  if(.not.saveprim .and. b0field .and. iw_e>0) then
4197  do ix3=ixcmin3,ixcmax3
4198  do ix2=ixcmin2,ixcmax2
4199  do ix1=ixcmin1,ixcmax1
4200  wc(ix1,ix2,ix3,iw_e)=sum( w(ix1:ix1+1,ix2:ix2+1,ix3,iw_e) &
4201  +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4202  ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
4203  ,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4204  ,:,0),dim=ndim+1) ) /dble(2**3)+&
4205  sum( w(ix1:ix1+1,ix2:ix2+1,ix3+1,iw_e) &
4206  +half*sum(ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4207  ,:,0)**2,dim=ndim+1) + sum( w(ix1:ix1+1,ix2:ix2+1,ix3&
4208  +1,iw_mag(:))*ps(igrid)%B0(ix1:ix1+1,ix2:ix2+1&
4209  ,:,0),dim=ndim+1) ) /dble(2**3)
4210  end do
4211  end do
4212  end do
4213  end if
4214 end if
4215 ! keep the coordinate and vector components
4216 xc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3) &
4217  = xc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:3)
4218 wc_tmp(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
4219  +nwauxio) = wc(ixcmin1:ixcmax1,ixcmin2:ixcmax2,ixcmin3:ixcmax3,1:nw&
4220  +nwauxio)
4221 xcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
4222  1:3) = xcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,&
4223  ixccmin3:ixccmax3,1:3)
4224 wcc_tmp(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,1:nw&
4225  +nwauxio) = wcc(ixccmin1:ixccmax1,ixccmin2:ixccmax2,ixccmin3:ixccmax3,&
4226  1:nw+nwauxio)
4227 end subroutine calc_grid23
4228 !=============================================================================
4229 subroutine save_connvtk23(qunit,igrid)
4230 
4231 ! this saves the basic line, pixel and voxel connectivity,
4232 ! as used by VTK file outputs for unstructured grid
4233 
4235 
4236 integer, intent(in) :: qunit, igrid
4237 
4238 integer :: nx1,nx2,nx3, nxC1,nxC2,nxC3, ix1,ix2,ix3
4239 !-----------------------------------------------------------------------------
4240 nx1=ixmhi1-ixmlo1+1;nx2=ixmhi2-ixmlo2+1;nx3=ixmhi1-ixmlo1+1;
4241 nxc1=nx1+1;nxc2=nx2+1;nxc3=nx3+1;
4242 
4243 do ix3=1,nx3
4244 do ix2=1,nx2
4245 do ix1=1,nx1
4246 
4247 
4249  write(qunit,'(8(i7,1x))')&
4250  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1-1, &
4251  (ix3-1)*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
4252  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1-1,&
4253  (ix3-1)*nxc2*nxc1+ ix2*nxc1+ix1,&
4254  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1-1,&
4255  ix3*nxc2*nxc1+(ix2-1)*nxc1+ix1,&
4256  ix3*nxc2*nxc1+ ix2*nxc1+ix1-1,&
4257  ix3*nxc2*nxc1+ ix2*nxc1+ix1
4258 
4259 end do
4260 end do
4261 end do
4262 end subroutine save_connvtk23
4263 !=============================================================================
4264 subroutine specialvar_output23(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
4265  ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,normconv)
4266 
4267 ! this subroutine can be used in convert, to add auxiliary variables to the
4268 ! converted output file, for further analysis using tecplot, paraview, ....
4269 ! these auxiliary values need to be stored in the nw+1:nw+nwauxio slots
4270 !
4271 ! the array normconv can be filled in the (nw+1:nw+nwauxio) range with
4272 ! corresponding normalization values (default value 1)
4273 
4275 
4276 integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
4277  ixImax3,ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3
4278 double precision, intent(in) :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
4279  ixImin3:ixImax3,1:3)
4280 double precision :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
4281  ixImin3:ixImax3,nw+nwauxio)
4282 double precision :: normconv(0:nw+nwauxio)
4284 double precision :: qvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir),&
4285  curlvec(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo1:ixGhi1,1:ndir)
4286 integer :: idirmin
4287 !-----------------------------------------------------------------------------
4288 ! output Te
4289 !if(saveprim)then
4290 ! w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+1)=w(ixOmin1:ixOmax1,&
4291 ! ixOmin2:ixOmax2,ixOmin3:ixOmax3,iw_e)/w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
4292 ! ixOmin3:ixOmax3,iw_rho)
4293 !endif
4294 !!! store current
4295 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,1)=w(ixImin1:ixImax1,&
4296 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(1))
4297 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,2)=w(ixImin1:ixImax1,&
4298 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(2))
4299 ! qvec(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3,3)=w(ixImin1:ixImax1,&
4300 ! ixImin2:ixImax2,ixImin3:ixImax3,mag(3));
4301 !call curlvector3D(qvec,ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,ixOmin1,&
4302 ! ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,curlvec,idirmin,1,ndir)
4303 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+2)=curlvec&
4304 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1)
4305 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+3)=curlvec&
4306 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,2)
4307 !w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,nw+4)=curlvec&
4308 ! (ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,3);
4309 end subroutine specialvar_output23
4310 \}
4311 !=============================================================================
4312 
subroutine resettree_convert
Force the tree to desired level(s) from level_io(_min/_max) used for conversion to vtk output.
Definition: amrgrid.t:69
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:194
integer function nodenumbertec3d(i1, i2, i3, nx1, nx2, nx3, ig, igrid)
Definition: convert.t:652
subroutine unstructuredvtkbsym23(qunit)
Definition: convert.t:3501
integer function nodenumbertec2d(i1, i2, nx1, nx2, ig, igrid)
Definition: convert.t:643
subroutine unstructuredvtk_mpi(qunit)
Definition: convert.t:1967
subroutine tecplot(qunit)
Definition: convert.t:379
subroutine write_vti(qunit, ixIL, ixCL, ixCCL, igD, nxD, normconv, wnamei, wC, wCC)
Definition: convert.t:2265
subroutine imagedatavtk_mpi(qunit)
Definition: convert.t:1699
subroutine unstructuredvtk(qunit)
Definition: convert.t:661
subroutine save_connvtk(qunit, igrid)
Definition: convert.t:1664
integer function nodenumbertec1d(i1, nx1, ig, igrid)
Definition: convert.t:634
subroutine onegrid(qunit)
Definition: convert.t:269
subroutine unstructuredvtkb23(qunit)
Definition: convert.t:3178
subroutine specialvar_output23(ixImin1, ixImin2, ixImin3, ixImax1, ixImax2, ixImax3, ixOmin1, ixOmin2, ixOmin3, ixOmax1, ixOmax2, ixOmax3, w, x, normconv)
Definition: convert.t:4285
subroutine oneblock(qunit)
Definition: convert.t:71
subroutine unstructuredvtkb64(qunit)
Definition: convert.t:1250
subroutine unstructuredvtkb(qunit)
Definition: convert.t:836
subroutine save_connvtk23(qunit, igrid)
Definition: convert.t:4249
subroutine punstructuredvtkb_mpi(qunit)
Definition: convert.t:2884
subroutine generate_plotfile
Definition: convert.t:3
subroutine punstructuredvtk_mpi(qunit)
Definition: convert.t:1873
subroutine tecplot_mpi(qunit)
Definition: convert.t:2400
subroutine write_vtk(qunit, ixIL, ixCL, ixCCL, igrid, nc, np, nxD, nxCD, normconv, wnamei, xC, xCC, wC, wCC)
Definition: convert.t:2156
subroutine write_pvtu(qunit)
Definition: convert.t:2325
subroutine save_conntec(qunit, igrid, igonlevel)
Definition: convert.t:586
subroutine calc_grid23(qunit, igrid, xC_TMP, xCC_TMP, wC_TMP, wCC_TMP, normconv, ixCmin1, ixCmin2, ixCmin3, ixCmax1, ixCmax2, ixCmax3, ixCCmin1, ixCCmin2, ixCCmin3, ixCCmax1, ixCCmax2, ixCCmax3, first, i3grid, d3grid, w, zlength, zgridsc)
Definition: convert.t:3995
Handles computations for coordinates and variables in output.
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.
subroutine getheadernames(wnamei, xandwnamei, outfilehead)
get all variables names
subroutine calc_x(igrid, xC, xCC)
computes cell corner (xC) and cell center (xCC) coordinates
subroutine convert_all()
Definition: mod_convert.t:120
Module with basic grid data structures.
Definition: mod_forest.t:2
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, dimension(:), allocatable, save morton_start
First Morton number per processor.
Definition: mod_forest.t:62
integer, dimension(:), allocatable, save morton_stop
Last Morton number per processor.
Definition: mod_forest.t:65
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
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(:), allocatable w_convert_factor
Conversion factors the primitive variables.
type(state), pointer block
Block pointer for using one block and its previous state.
logical nocartesian
IO switches for conversion.
double precision global_time
The global simulation time.
integer type_block_xc_io
MPI type for IO: cell corner (xc) or cell center (xcc) coordinates.
integer snapshotini
Resume from the snapshot with this index.
logical saveprim
If true, convert from conservative to primitive variables in output.
character(len=std_len) convert_type
Which format to use when converting.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
double precision time_convert_factor
Conversion factor for time unit.
integer icomm
The MPI communicator.
integer, dimension(:), allocatable ng
number of grid blocks in domain per dimension, in array over levels
integer mype
The rank of the current MPI task.
integer type_block_io
MPI type for IO: block excluding ghost cells.
integer, parameter plevel_
double precision length_convert_factor
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with ghost cells)
integer ierrmpi
A global MPI error return code.
logical autoconvert
If true, already convert to output format during the run.
logical slab
Cartesian geometry or not.
integer type_block_wc_io
MPI type for IO: cell corner (wc) or cell center (wcc) variables.
integer snapshotnext
IO: snapshot and collapsed views output numbers/labels.
integer npe
The number of MPI tasks.
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision, dimension(ndim, 2) writespshift
integer, parameter unitterm
Unit for standard output.
logical, dimension(:), allocatable w_write
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
logical, dimension(:), allocatable writelevel
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer, parameter unitconvert
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
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:59
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:27
procedure(sub_check_params), pointer phys_te_images
Definition: mod_physics.t:91
Module with all the methods that users can customize in AMRVAC.
procedure(aux_output), pointer usr_aux_output
procedure(special_convert), pointer usr_special_convert
Pointer to a tree_node.
Definition: mod_forest.t:6