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