! $Id: pdf_xyz.f90,v 1.3 2004/04/02 07:15:35 nilshau Exp $
!
!***********************************************************************
program pdf
  !
  !
  !-----------------------------------------------------------------------
  !   24-july-03/nils: code
  !
  use Cdata
  use Cparam
  use Sub
  use Deriv
  use General
  !
  implicit none
  !
  !
  !  define parameters
  !
  character (len=130) ::  file,pdf_result_file,file_ext,pdf_in_file
  character (len=4) :: ch,csnap
  integer :: cpu,nnn
  real :: dummy,ddx,ddy,ddz
  real, dimension(3) :: minarr,maxarr,ddarr
  real :: minx=1e6, miny=1e6, minz=1e6, maxx=-1e6,maxy=-1e6,maxz=-1e6
  real, dimension (mx,my,mz,mvar+maux) :: f
  real, dimension (nx,3) :: bb
  integer, dimension (3,nx) :: ppdf
  integer, dimension (nx) :: ppdfx=0,ppdfy=0,ppdfz=0
  integer :: snap,snap_start,snap_stop
  integer :: ll,i,dir_nr,dir_count,jj,ii,indeks
  !
  ! Setting some constants
  !
  iux=1 
  iuy=2
  iuz=3
  ilnrho=4
  iax=5
  iay=6
  iaz=7
  !
  ip=10
  !
  ! Reading pdf.in
  !
  call safe_character_assign(pdf_in_file,'pdf.in')
  if (ip<8) print*,'Reading '//trim(pdf_in_file)
  open(1,FILE=pdf_in_file,FORM='formatted')
  read(1,*) snap_start,snap_stop
  close(1)
  nnn=nx
  !
  !  Loop over snapshots
  !
  do snap=snap_start,snap_stop
    !
    ! Initialize
    !
    minx=1e6
    miny=1e6
    minz=1e6
    maxx=-1e6
    maxy=-1e6
    maxz=-1e6
    ppdfx=0
    ppdfy=0
    ppdfz=0
    !
    ! Read through all processors to find global max and min values
    !
    do cpu=0,ncpus-1
      print*,'Reading processor',cpu+1,'of',ncpus,'processors'
      ipy = modulo(cpu, nprocy)
      ipz = cpu/(nprocy)
      !
      ! Read var.dat 
      !
      call chn(cpu,ch)
      call chn(snap,csnap)
      if (snap_stop .eq. 0) call safe_character_assign(file,trim(datadir) &
           //'/proc'//trim(ch)//'/var.dat')
      if (snap_stop .ne. 0) call safe_character_assign(file,trim(datadir) &
           //'/proc'//trim(ch)//'/VAR'//trim(csnap))
      if (ip<8) print*,'Reading '//trim(file)
      open(1,file=file,form='unformatted')
      read(1) f(1:mx,1:my,1:mz,1:mvar)
      if (lshear) then
        read(1) t,x,y,z,dx,dy,dz,deltay
      else
        read(1) t,x,y,z,dx,dy,dz
      endif
      close(1)
      !
      ! Find min and max values
      !
      do m=m1,m2
        do n=n1,n2
          call curl(f,iax,bb)
          if (minval(bb(:,1)) .lt. minx) minx = minval(bb(:,1))
          if (minval(bb(:,2)) .lt. miny) miny = minval(bb(:,2))
          if (minval(bb(:,3)) .lt. minz) minz = minval(bb(:,3))
          if (maxval(bb(:,1)) .gt. maxx) maxx = maxval(bb(:,1))
          if (maxval(bb(:,2)) .gt. maxy) maxy = maxval(bb(:,2))
          if (maxval(bb(:,3)) .gt. maxz) maxz = maxval(bb(:,3))
        enddo
      enddo
    enddo
    !
    ! Must make max and min values slightly larger and smaller 
    ! in order for all values to fall inside the domain.
    !
    minx=minx-abs(minx)/nnn ; minarr(1)=minx
    miny=miny-abs(miny)/nnn ; minarr(2)=miny
    minz=minz-abs(minz)/nnn ; minarr(3)=minz
    maxx=maxx*(1+1./nnn) ; maxarr(1)=maxx
    maxy=maxy*(1+1./nnn) ; maxarr(2)=maxy
    maxz=maxz*(1+1./nnn) ; maxarr(3)=maxz
    !
    ! Find size of each bin
    !
    ddx=(maxx-minx)/nnn  ; ddarr(1)=ddx
    ddy=(maxy-miny)/nnn  ; ddarr(2)=ddy
    ddz=(maxz-minz)/nnn  ; ddarr(3)=ddz
    !
    ! Read through all processors to find global pdf's
    !
    do cpu=0,ncpus-1
      ppdf=0
      print*,'Reading processor',cpu+1,'of',ncpus,'processors'
      ipy = modulo(cpu, nprocy)
      ipz = cpu/(nprocy)
      !
      ! Read var.dat 
      !
      call chn(cpu,ch)
      call chn(snap,csnap)
      if (snap_stop .eq. 0) call safe_character_assign(file,trim(datadir) &
           //'/proc'//trim(ch)//'/var.dat')
      if (snap_stop .ne. 0) call safe_character_assign(file,trim(datadir) &
           //'/proc'//trim(ch)//'/VAR'//trim(csnap))
      if (ip<8) print*,'Reading '//trim(file)
      open(1,file=file,form='unformatted')
      read(1) f(1:mx,1:my,1:mz,1:mvar)
      if (lshear) then
        read(1) t,x,y,z,dx,dy,dz,deltay
      else
        read(1) t,x,y,z,dx,dy,dz
      endif
      close(1)
      !
      ! Find PDF
      !
      do m=m1,m2
        do n=n1,n2
          call curl(f,iax,bb)
          do ii=1,3
            do jj=1,nx
              indeks=ceiling((bb(jj,ii)-minarr(ii))/ddarr(ii))
              ppdf(ii,indeks)=ppdf(ii,indeks)+1
            enddo
          enddo
        enddo
      enddo
      !
      ! Collect results from all processors
      !
      ppdfx=ppdfx+ppdf(1,:)
      ppdfy=ppdfy+ppdf(2,:)
      ppdfz=ppdfz+ppdf(3,:)
      !
      !  Finnishing loop over processors
      !
    enddo
    !
    !  Saving the results
    !
    if (snap_stop .eq. 0) call safe_character_assign(pdf_result_file, &
         trim(datadir)//'/pdf_xyz.dat')
    if (snap_stop .ne. 0) call safe_character_assign(pdf_result_file, &
         trim(datadir)//'/pdf_xyz_VAR'//trim(csnap)//'.dat')
    close(1)
    if (ip<8) print*,'Writing '//trim(pdf_result_file)
    open(1,file=pdf_result_file)
    write(1,'(3e12.4)') t
    write(1,'(I)') nnn
    write(1,'(3e12.4)') minx,maxx
    write(1,'(3I)') ppdfx
    write(1,'(3e12.4)') miny,maxy
    write(1,'(3I)') ppdfy
    write(1,'(3e12.4)') minz,maxz
    write(1,'(3I)') ppdfz
    close(1)
    !
    !  Finnishing loop over snapshots (if snap_stop .ne. 0)
    !
  enddo
  !
end program pdf