! $Id: slices.f90 13524 2010-03-24 08:50:28Z tavo.buk $
!
! This module produces slices for animation purposes.
!
module Slices
!
use Cdata
use Messages
use Sub
!
implicit none
!
private
!
public :: wvid, wvid_prepare, setup_slices, wslice
!
real, target, dimension (nx,ny) :: slice_xy=0.0, slice_xy2=0.0
real, target, dimension (nx,nz) :: slice_xz=0.0
real, target, dimension (ny,nz) :: slice_yz=0.0
real, target, dimension (nx,ny) :: slice_xy3=0.0, slice_xy4=0.0
!
real, public :: tvid
integer, public :: nvid
real :: tslice
!
contains
!***********************************************************************
subroutine wvid_prepare
!
! Prepare lvideo for writing slices into video file
! This is useful for visualization of scalar field (or one component
! of a vector field) on the perifery of a box.
! Can be visualized in idl using rvid_box.pro
!
! 20-oct-97/axel: coded
! 08-oct-02/tony: increased size of file to handle datadir//'/tvid.dat'
! 13-nov-02/axel: added more fields, use wslice.
! 18-mar-03/axel: added dust velocity
!
integer, save :: ifirst=0
!
character (len=5) :: ch
character (len=130) :: file
!
! Output vid-data in 'tvid' time intervals
!
file = trim(datadir)//'/tvid.dat'
if (ifirst==0) then
call read_snaptime(trim(file),tvid,nvid,dvid,t)
ifirst=1
endif
!
! This routine sets lvideo=T whenever its time to write a slice
!
call update_snaptime(file,tvid,nvid,dvid,t,lvideo,ch,ENUM=.false.)
!
! Save current time so that the time that is written out in wslice() is not
! from the next time step
!
if (lvideo) tslice = t
!
endsubroutine wvid_prepare
!***********************************************************************
subroutine wvid(f,path)
!
! Write slices for animation of scalar field
! (or one component of a vector field) on the perifery of a box.
! Can be visualized in idl using rvid_box.pro.
!
! 13-nov-02/axel: added more fields, use wslice.
! 22-sep-07/axel: changed Xy to xy2, to be compatible with Mac
!
use Density, only: get_slices_density
use EquationOfState, only: get_slices_eos
use Entropy, only: get_slices_entropy
use General
use Hydro, only: get_slices_hydro
use Magnetic, only: get_slices_magnetic
use Shock, only: get_slices_shock
!
real, dimension (mx,my,mz,mfarray) :: f
character(len=*) :: path
!
type (slice_data) :: slices
character(len=5) :: sindex
logical, save :: lfirstloop=.true.
logical :: lnewfile=.true.
logical :: lslices_legacy=.true.
integer :: inamev
!
slices%ix =ix_loc
slices%iy =iy_loc
slices%iz =iz_loc
slices%iz2=iz2_loc
slices%index=0
!
! Loop over slices.
!
inamev=1
do while (inamev <= nnamev)
slices%ready=.false.
slices%xy=>slice_xy
slices%xz=>slice_xz
slices%yz=>slice_yz
slices%xy2=>slice_xy2
!
slices%name=trim(cnamev(inamev))
!
! By default assume we're not using module hooks to get the slice contents
!
lslices_legacy=.true.
!
! Get slice information from the modules.
!
lslices_legacy=.false.
if (ldensity) call get_slices_density (f,slices)
if (lentropy) call get_slices_entropy (f,slices)
if (leos) call get_slices_eos (f,slices)
if (lhydro) call get_slices_hydro (f,slices)
if (lmagnetic) call get_slices_magnetic (f,slices)
if (lshock) call get_slices_shock (f,slices)
!
if (lslices_legacy) then
inamev=inamev+1
cycle
endif
!
! Slice, or component of slice, ready for saving.
!
if (slices%ready) then
if (slices%index==0) then ! If this wasn't a multi slice...
if (associated(slices%yz)) &
call wslice(path//trim(slices%name)//'.yz',slices%yz, &
x(slices%ix),ny,nz)
if (associated(slices%xz)) &
call wslice(path//trim(slices%name)//'.xz',slices%xz, &
y(slices%iy),nx,nz)
if (associated(slices%xy)) &
call wslice(path//trim(slices%name)//'.xy',slices%xy, &
z(slices%iz),nx,ny)
if (associated(slices%xy2)) &
call wslice(path//trim(slices%name)//'.xy2',slices%xy2, &
z(slices%iz2),nx,ny)
inamev=inamev+1
else
call chn(slices%index, sindex)
if (associated(slices%yz)) &
call wslice(path//trim(slices%name)//trim(sindex)//'.yz', &
slices%yz, x(slices%ix),ny,nz)
if (associated(slices%xz)) &
call wslice(path//trim(slices%name)//trim(sindex)//'.xz', &
slices%xz, y(slices%iy),nx,nz)
if (associated(slices%xy)) &
call wslice(path//trim(slices%name)//trim(sindex)//'.xy', &
slices%xy, z(slices%iz),nx,ny)
if (associated(slices%xy2)) &
call wslice(path//trim(slices%name)//trim(sindex)//'.xy2', &
slices%xy2, z(slices%iz2),nx,ny)
endif
else
if (slices%index==0) then ! If this wasn't a multi slice...
if (lfirstloop.and.lroot) then
if (lnewfile) then
open(1,file='video.err')
lnewfile=.false.
else
open(1,file='video.err',position='append')
endif
write(1,*) 'unknown slice: ',trim(cnamev(inamev))
close(1)
endif
else
slices%index=0
endif
inamev=inamev+1
endif
enddo
!
lfirstloop=.false.
!
endsubroutine wvid
!***********************************************************************
subroutine wslice(filename,a,pos,ndim1,ndim2)
!
! appending to an existing slice file
!
! 12-nov-02/axel: coded
! 26-jun-06/anders: moved from Slices
! 22-sep-07/axel: changed Xy to xy2, to be compatible with Mac
!
integer :: ndim1,ndim2
character (len=*) :: filename
real, dimension (ndim1,ndim2) :: a
real, intent(in) :: pos
!
! check whether we want to write a slice on this processor
!
if ( (lwrite_slice_xy .and.index(filename,'xy' )>0) .or. &
(lwrite_slice_xy2.and.index(filename,'xy2')>0) .or. &
(lwrite_slice_xz .and.index(filename,'xz' )>0) .or. &
(lwrite_slice_yz .and.index(filename,'yz' )>0) ) then
open(1,file=filename,form='unformatted',position='append')
write(1) a,tslice,pos
close(1)
endif
!
endsubroutine wslice
!***********************************************************************
subroutine setup_slices()
!
! Determine slice positions and whether slices are to be written on this
! processor.
!
! 29-may-06/tobi: wrapped code from param_io.f90 into this subroutine
!
! set slice position. The default for slice_position is 'p' for periphery,
! although setting ix, iy, iz, iz2 by hand will overwrite this.
! If slice_position is not 'p', then ix, iy, iz, iz2 are overwritten.
!
if (slice_position=='p' .or. slice_position=='S') then
ix_loc=l1; iy_loc=m1; iz_loc=n1; iz2_loc=n2
lwrite_slice_xy2=(ipz==nprocz-1)
lwrite_slice_xy=(ipz==0)
lwrite_slice_xz=(ipy==0)
lwrite_slice_yz=.true.
!
! slice position in middle of the box independ of nprocy,nprocz
! second horizontal slice is the upper most layer
!
elseif (slice_position=='m') then
ix_loc=(l1+l2)/2
if (mod(nprocy,2)==0) then; iy_loc=m1; else; iy_loc=(m1+m2)/2; endif
if (mod(nprocz,2)==0) then; iz_loc=n1; else; iz_loc=(n1+n2)/2; endif
iz2_loc=n2
lwrite_slice_xy2=(ipz==nprocz-1)
lwrite_slice_xy=(ipz==nprocz/2)
lwrite_slice_xz=(ipy==nprocy/2)
lwrite_slice_yz=.true.
!
! slice positions for spherical coordinates
! w is for "wedges" since the outputs are
! the midplane (rphi,theta=y(mpoint)) and four
! wedges in rtheta (xy)
!
elseif (slice_position=='w') then
!midplane slices
!ix_loc=nxgrid/2+nghost
iy = nygrid/2+nghost
!meridional wedges, at 4 different
!equally spaced azimuthal locations
iz = 0*nzgrid/4+1+nghost
iz2= 1*nzgrid/4+1+nghost
!yz is not needed
lwrite_slice_yz =.false.
!
! Another slice positions for spherical coordinates
! s is for "surface" meaning theta-phi sections
! keep iz_loc=n1, corresponding to a meridional slice on n=n1.
! The lwrite_slice_xz=.true. is needed for the read_video routine.
!
elseif (slice_position=='s') then
iz_loc=n1; iz2_loc=n2
call xlocation(xtop_slice,ix_loc,lwrite_slice_yz)
lwrite_slice_xy2=(ipz==nprocz/4)
lwrite_slice_xy=(ipz==0)
lwrite_slice_xz=.true.
!
! later we may also want to write other slices
!
!call xlocation(xtop_slice,ix2_loc,lwrite_slice_yz2)
!call ylocation(ytop_slice,iy2_loc,lwrite_slice_xz2)
!
! slice position when the first meshpoint in z is the equator (sphere)
! For one z-processor, iz remains n1, but iz2 is set to the middle.
!
! TH: A certain processor layout is implied here
!
elseif (slice_position=='e') then
ix_loc=(l1+l2)/2; iy_loc=m1; iz_loc=n1; iz2_loc=n2
if (nprocy==1) then; iy_loc=(m1+m2)/2; endif
if (nprocz==1) then; iz2_loc=(iz+n2)/2; endif
lwrite_slice_xy2=(ipz==nprocz/4)
lwrite_slice_xy=(ipz==0)
lwrite_slice_xz=(ipy==nprocy/2)
lwrite_slice_yz=.true.
!
! slice position similar to periphery (p), but the two z positions
! can now be given in terms of z (zbot_slice, ztop_slice).
!
elseif (slice_position=='c') then
ix_loc=l1; iy_loc=m1
call zlocation(zbot_slice,iz_loc,lwrite_slice_xy)
call zlocation(ztop_slice,iz2_loc,lwrite_slice_xy2)
lwrite_slice_xz=(ipy==0)
lwrite_slice_yz=.true.
!
! periphery of the box, but the other way around
!
elseif (slice_position=='q') then
ix_loc=l2; iy_loc=m2; iz_loc=n2; iz2_loc=n1
lwrite_slice_xy2=(ipz==0)
lwrite_slice_xy=(ipz==nprocz-1)
lwrite_slice_xz=(ipy==nprocy-1)
lwrite_slice_yz=.true.
else
if (lroot) then
call fatal_error('setup_slices', &
'No such value for slice_position: '//slice_position)
endif
endif
!
! Overwrite slice postions if any ix,iy,iz,iz2 is greater then Zero
!
if (ix>0) then
ix_loc=ix-ipx*nx
if (ix_loc>=l1.and.ix_loc<=l2) then
lwrite_slice_yz=.true.
else
lwrite_slice_yz=.false.
endif
endif
!
if (iy>0) then
iy_loc=iy-ipy*ny
if (iy_loc>=m1.and.iy_loc<=m2) then
lwrite_slice_xz=.true.
else
lwrite_slice_xz=.false.
endif
endif
!
if (iz>0) then
iz_loc=iz-ipz*nz
if (iz_loc>=n1.and.iz_loc<=n2) then
lwrite_slice_xy=.true.
else
lwrite_slice_xy=.false.
endif
endif
!
if (iz2>0) then
iz2_loc=iz2-ipz*nz
if (iz2_loc>=n1.and.iz2_loc<=n2) then
lwrite_slice_xy2=.true.
else
lwrite_slice_xy2=.false.
endif
endif
!
! write slice position to a file (for convenient post-processing)
!
if (lroot) then
open(1,file=trim(datadir)//'/slice_position.dat',STATUS='unknown')
write(1,'(a)') slice_position
close(1)
endif
!
! write the ipz processor numbers for the two slices
! The first number (=ipz) is essential, the others just for interest.
!
if (lwrite_slice_xy.and.ipy==0) then
open(1,file=trim(directory)//'/zbot_procnum.dat',STATUS='unknown')
write(1,'(2i5,e12.4)') ipz,iz_loc,z(iz_loc)
close(1)
endif
if (lwrite_slice_xy2.and.ipy==0) then
open(1,file=trim(directory)//'/ztop_procnum.dat',STATUS='unknown')
write(1,'(2i5,e12.4)') ipz,iz2_loc,z(iz2_loc)
close(1)
endif
!
! make sure ix_loc,iy_loc,iz_loc,iz2_loc
! are not outside the boundaries
!
ix_loc=min( ix_loc,l2) ; iy_loc=min( iy_loc,m2)
ix_loc=max( ix_loc,l1) ; iy_loc=max( iy_loc,m1)
iz_loc=min( iz_loc,n2) ; iz2_loc=min(iz2_loc,n2)
iz_loc=max( iz_loc,n1) ; iz2_loc=max(iz2_loc,n1)
!
if (lroot) then
write (*,*)'read_runpars: slice_position = '//slice_position
write (*,'(1x,a,4i4)') &
'read_runpars: ix,iy,iz,iz2 (video files) =',&
ix,iy,iz,iz2
endif
!
endsubroutine setup_slices
!***********************************************************************
subroutine xlocation(xpos,ixpos,lproc)
!
! if xpos lies within this processor, then lproc=T and xpos=x(ixpos).
! Otherwise lproc=F and ixpos=1.
!
! 18-nov-06/axel: coded
! 14-oct-08/ccyang: use half-closed interval and include the top-most plane
! 20-feb-10/axel: adapted from zlocation
!
real :: xpos
integer :: ixpos,l
logical :: lproc
!
! run through all x positions until we hit the right interval.
! If the right interval is found, jump out of the loop.
!
do l=l1,l2
if (x(l)<=xpos.and.x(l+1)>xpos) then
ixpos=l
lproc=.true.
goto 900
else
endif
enddo
!
! if nothing is found, we set lproc=.false. and
! and put ixpos=1
!
ixpos=1
lproc=.false.
!
900 continue
!
endsubroutine xlocation
!***********************************************************************
subroutine ylocation(ypos,iypos,lproc)
!
! if ypos lies within this processor, then lproc=T and ypos=y(iypos).
! Otherwise lproc=F and iypos=1.
!
! 18-nov-06/axel: coded
! 14-oct-08/ccyang: use half-closed interval and include the top-most plane
! 20-feb-10/axel: adapted from xlocation
!
real :: ypos
integer :: iypos,m
logical :: lproc
!
! run through all y positions until we hit the right interval.
! If the right interval is found, jump out of the loop.
!
do m=m1,m2
if (y(m)<=ypos.and.y(m+1)>ypos) then
iypos=m
lproc=.true.
goto 900
else
endif
enddo
!
! if nothing is found, we set lproc=.false. and
! and put iypos=1
!
iypos=1
lproc=.false.
!
900 continue
!
endsubroutine ylocation
!***********************************************************************
subroutine zlocation(zpos,izpos,lproc)
!
! if zpos lies within this processor, then lproc=T and zpos=z(izpos).
! Otherwise lproc=F and izpos=1.
!
! 18-nov-06/axel: coded
! 14-oct-08/ccyang: use half-closed interval and include the top-most plane
!
real :: zpos
integer :: izpos,n
logical :: lproc
!
! run through all z positions until we hit the right interval.
! If the right interval is found, jump out of the loop.
!
do n=n1,n2
if (z(n)<=zpos.and.z(n+1)>zpos) then
izpos=n
lproc=.true.
goto 900
else
endif
enddo
!
! if nothing is found, we set lproc=.false. and
! and put izpos=1
!
izpos=1
lproc=.false.
!
900 continue
!
endsubroutine zlocation
!***********************************************************************
endmodule Slices