!----------------------------------------------------------------------------
! $Id: tempstorage.f90,v 1.4 1997/09/12 23:37:56 mijan Exp $
!----------------------------------------------------------------------------
! Module for keeping track of temporary storage.
!
! All temporary variables are allocated, deallocated and flagged for
!  use through this module.
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
! NOTE: The tempscalar, etc. and freescalar, etc. functions do not
!   allocate or deallocate large chunks of memory (i.e. 3D arrays). They
!   only grant and remove "permission" to use certain temporary variables.
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!
! ****************Setup/clear/inquiry subroutines:******************
!
!
!  subroutine setuptemps(n,nx,ny,nz)
!     -Allocates space for n temporary gridfunctions of size (nx,ny,nz)
!
!  subroutine cleartemps()
!     -Frees all the space allocated by setuptemps()
!
!  subroutine printtemps(string)
!     -Prints "String", number of temps allocated, (total number ever used)
!
!
!
! ****************Calls for using temporary storage:******************
!
!
!  s => tempscalar(), subroutine freescalar(s)
!     -Sets up or frees a temporary scalar s
!       s declared as type(Gridfunction), pointer :: s
!
!  v => tempvector(), subroutine freevector(v)
!     -Sets up or frees a temporary vector v
!       v declared as type(Gridfunction), pointer, dimension(:) :: v
!
!  t => tempsymmtensor(), subroutine freesymmtensor(t)
!     -Sets up or frees a temporary symmetric tensor t
!       t declared as type(Gridfunction), pointer, dimension(:,:) :: t
!
!  c => tempconnection(), subroutine freeconnection(c)
!     -Sets up or frees a temporary connection c
!       c declared as type(Gridfunction), pointer, dimension(:,:,:) :: c
!
!  r => tempriemann(), subroutine freeriemann(r)
!     -Sets up or frees a temporary riemann tensor r
!       r declared as type(Gridfunction), pointer, dimension(:,:,:,:) :: r
!
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
! BEWARE:  Because Fortran allows one to pass a n-dimensional array into
!  a subroutine which expects a one-dimensional array, the compiler will
!  NOT catch bugs where the free*** routines are called with an argument of
!  the wrong type, e.g. calling freevector(t) where t is a symmetric tensor.
!  Such a bug will either:
!      -Try to free a temp that is already free, and you get an error message.
!      -Fail to free a temp that should be freed, so eventually you will
!        run out of usable temps and get an error message.
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
!----------------------------------------------------------------------------
module tempstorage
  type tempelem
    real*8, dimension(:,:,:), pointer :: d      ! Data
    integer                           :: in_use ! Is it currently used?
  end type tempelem

  type(tempelem), dimension(:), allocatable,save :: temps
  integer,save                                   :: maxtemps
  integer,save                                   :: max_used,num_used
!##  !----------------------------------------------------------------
!##  ! NaN-ification                                                 !
!##  real*8, parameter                              :: NaN = 0.d0/0.d0
!##  !----------------------------------------------------------------
  
  PRIVATE
  PUBLIC :: setuptemps,cleartemps,printtemps
  PUBLIC :: tempscalar,tempvector,tempsymmtensor,tempconnection,tempriemann
  PUBLIC :: freescalar,freevector,freesymmtensor,freeconnection,freeriemann
  PUBLIC :: tempnvector
  PUBLIC :: freenvector

  
contains

  !----------------\
  ! Print function  >
  !----------------/

  subroutine printtemps(string)
    implicit none
    character (LEN=*) :: string
    
!##    integer :: i
!##    
!##    do i = 1,maxtemps
!##       if(temps(i)%in_use == 1) then
!##	  write(*,*)string,':Temp used:',i
!##       endif
!##    enddo
    write(*,'(A,": Ntemps=",I3," (",I3,")")')string,num_used,max_used
    
  end subroutine printtemps
  
  !-------------------------\
  ! Initialization function  >
  !-------------------------/
  
  subroutine setuptemps(ntemps,nx,ny,nz)
    implicit none
    integer :: ntemps,nx,ny,nz
    
    integer :: i

!##    !---------------------------------------
!##    ! Randomizer                           !
!##    integer,parameter :: iseed = -5670246  !
!##    call srand(iseed)                      !
!##    !---------------------------------------
    
    maxtemps   = ntemps
    num_used   = 0
    max_used   = 0
    
    allocate(temps(ntemps))

    do i = 1, ntemps
      allocate(temps(i)%d(nx,ny,nz))
      temps(i)%in_use = 0
    enddo
  end subroutine setuptemps


  !-------------------------\
  ! Clears the storage       >
  !-------------------------/
  
  subroutine cleartemps()
    implicit none
    integer :: i

    do i = 1,maxtemps
      deallocate(temps(i)%d)
    enddo

    deallocate(temps)
    
  end subroutine cleartemps

  !-------------\\
  !--------------\\
  ! Scalars       >>
  !--------------//
  !-------------//

  
  function tempscalar() RESULT(gf)
    use tensors
    implicit none
    type(GridFunction), pointer :: gf

    allocate(gf)
    
    gf%d => setelement()

  end function tempscalar

  subroutine freescalar(gf)
    use tensors
    implicit none
    type(GridFunction), pointer :: gf

    call freeelement(gf%d)
    deallocate(gf)
  end subroutine freescalar


  !-------------\\
  !--------------\\
  ! Vectors       >>
  !--------------//
  !-------------//

  
  function tempvector() RESULT(v)
    use tensors
    implicit none
    type(GridFunction), dimension(:), pointer :: v
    
    integer :: i
    
    allocate(v(3))
    
    do i = 1,3
      v(i)%d => setelement()
    enddo
       
  end function tempvector

  subroutine freevector(v)
    use tensors
    implicit none
    type(GridFunction), dimension(:), pointer :: v

    integer :: i

    do i = 1,3
      call freeelement(v(i)%d)
    enddo
    deallocate(v)
  end subroutine freevector
  function tempnvector(n) RESULT(v)
    use tensors
    implicit none
    integer                                   :: n
    type(GridFunction), dimension(:), pointer :: v

    integer :: i

    allocate(v(n))

    do i = 1,n
      v(i)%d => setelement()
    enddo
      
  end function tempnvector

  subroutine freenvector(v)
    use tensors
    implicit none
    type(GridFunction), dimension(:), pointer :: v

    integer :: i, shp(1)

    shp = shape(v)
    do i = 1,shp(1)
      call freeelement(v(i)%d)
    enddo
    deallocate(v)
  end subroutine freenvector



  !-------------\\
  !--------------\\
  ! Symm Tensors  >>
  !--------------//
  !-------------//
  
  function tempsymmtensor() RESULT(t)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:), pointer :: t
    
    integer :: i,j
    
    allocate(t(3,3))
    
    do j = 1,3
      do i = j,3
	t(i,j)%d => setelement()
	if (i /= j) t(j,i)%d => t(i,j)%d
      enddo
    enddo
       
  end function tempsymmtensor
  
  subroutine freesymmtensor(t)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:), pointer :: t

    integer :: i,j

    do j = 1,3
      do i = j,3
	call freeelement(t(i,j)%d)
      enddo
    enddo
    deallocate(t)
  end subroutine freesymmtensor

  !-------------\\
  !--------------\\
  ! Connection    >>
  !--------------//
  !-------------//
  
  function tempconnection() RESULT(c)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:,:), pointer :: c
    
    integer :: i,j,k
    
    allocate(c(3,3,3))

    do k=1,3
      do j = 1,3
	do i = j,3
	  c(k,i,j)%d => setelement()
	  if (i /= j) c(k,j,i)%d => c(k,i,j)%d
	enddo
      enddo
    enddo
       
  end function tempconnection

  subroutine freeconnection(c)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:,:), pointer :: c

    integer :: i,j,k

    do k = 1,3
      do j = 1,3
	do i = j,3
	  call freeelement(c(k,i,j)%d)
	enddo
      enddo
    enddo
    deallocate(c)
  end subroutine freeconnection

  
  !-------------\\
  !--------------\\
  ! Riemann       >>
  !--------------//
  !-------------//
  
  function tempriemann() RESULT(r)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:,:,:), pointer :: r
    
    real*8, dimension(:,:,:), pointer :: r0,r1,r2,r3,r4,r5
    
    allocate(r(3,3,3,3))

    r0 => setelement()
    r1 => setelement()
    r2 => setelement()
    r3 => setelement()
    r4 => setelement()
    r5 => setelement()

    r(1,1,1,1)%d => r0
    r(1,1,1,2)%d => r0
    r(1,1,1,3)%d => r0
    r(1,1,2,1)%d => r0
    r(1,1,2,2)%d => r0
    r(1,1,2,3)%d => r0
    r(1,1,3,1)%d => r0
    r(1,1,3,2)%d => r0
    r(1,1,3,3)%d => r0
    r(1,2,1,1)%d => r0
    r(1,2,1,2)%d => r0
    r(1,2,1,3)%d => r1
    r(1,2,2,1)%d => r0
    r(1,2,2,2)%d => r0
    r(1,2,2,3)%d => r2
    r(1,2,3,1)%d => r1
    r(1,2,3,2)%d => r2
    r(1,2,3,3)%d => r0
    r(1,3,1,1)%d => r0
    r(1,3,1,2)%d => r1
    r(1,3,1,3)%d => r3
    r(1,3,2,1)%d => r1
    r(1,3,2,2)%d => r0
    r(1,3,2,3)%d => r4
    r(1,3,3,1)%d => r3
    r(1,3,3,2)%d => r4
    r(1,3,3,3)%d => r0
    r(2,1,1,1)%d => r0
    r(2,1,1,2)%d => r0
    r(2,1,1,3)%d => r1
    r(2,1,2,1)%d => r0
    r(2,1,2,2)%d => r0
    r(2,1,2,3)%d => r2
    r(2,1,3,1)%d => r1
    r(2,1,3,2)%d => r2
    r(2,1,3,3)%d => r0
    r(2,2,1,1)%d => r0
    r(2,2,1,2)%d => r0
    r(2,2,1,3)%d => r0
    r(2,2,2,1)%d => r0
    r(2,2,2,2)%d => r0
    r(2,2,2,3)%d => r0
    r(2,2,3,1)%d => r0
    r(2,2,3,2)%d => r0
    r(2,2,3,3)%d => r0
    r(2,3,1,1)%d => r0
    r(2,3,1,2)%d => r2
    r(2,3,1,3)%d => r4
    r(2,3,2,1)%d => r2
    r(2,3,2,2)%d => r0
    r(2,3,2,3)%d => r5
    r(2,3,3,1)%d => r4
    r(2,3,3,2)%d => r5
    r(2,3,3,3)%d => r0
    r(3,1,1,1)%d => r0
    r(3,1,1,2)%d => r1
    r(3,1,1,3)%d => r3
    r(3,1,2,1)%d => r1 
    r(3,1,2,2)%d => r0
    r(3,1,2,3)%d => r4
    r(3,1,3,1)%d => r3
    r(3,1,3,2)%d => r4
    r(3,1,3,3)%d => r0
    r(3,2,1,1)%d => r0
    r(3,2,1,2)%d => r2
    r(3,2,1,3)%d => r4
    r(3,2,2,1)%d => r2
    r(3,2,2,2)%d => r0
    r(3,2,2,3)%d => r5
    r(3,2,3,1)%d => r4
    r(3,2,3,2)%d => r5
    r(3,2,3,3)%d => r0
    r(3,3,1,1)%d => r0
    r(3,3,1,2)%d => r0
    r(3,3,1,3)%d => r0
    r(3,3,2,1)%d => r0
    r(3,3,2,2)%d => r0
    r(3,3,2,3)%d => r0
    r(3,3,3,1)%d => r0
    r(3,3,3,2)%d => r0
    r(3,3,3,3)%d => r0
       
  end function tempriemann

  
  subroutine freeriemann(r)
    use tensors
    implicit none
    type(GridFunction), dimension(:,:,:,:), pointer :: r

    call freeelement(r(1,2,1,2)%d)
    call freeelement(r(1,2,1,3)%d)
    call freeelement(r(1,2,2,3)%d)
    call freeelement(r(1,3,1,3)%d)
    call freeelement(r(1,3,2,3)%d)
    call freeelement(r(2,3,2,3)%d)
    
    deallocate(r)
  end subroutine freeriemann

  
  !---------------------------------------------\\
  !----------------------------------------------\\
  ! Private setelement and freeelement functions  >>
  !----------------------------------------------//
  !---------------------------------------------//
  
  function setelement() RESULT(tempptr)
    implicit none
    real*8, dimension(:,:,:), pointer :: tempptr

    integer :: i
    
!##    !-----------------------------------------------
!##    ! Randomizer                                   !
!##    integer :: whichone,j                          !
!##    real*8  :: rand                                !
!##    whichone = CEILING(rand()*(maxtemps-num_used)) !
!##    j=0                                            !
!##    !-----------------------------------------------

    do i = 1, maxtemps
      if(temps(i)%in_use==0) then
!##	!--------------------------------
!##	! Randomizer                    !
!##	j=j+1; if(j <  whichone) cycle  !
!##	!--------------------------------
	temps(i)%in_use=1
	tempptr => temps(i)%d
!##	!--------------------------------
!##	! NaN-ification                 !
!##	temps(i)%d = NaN                !
!##	!--------------------------------
	num_used = num_used + 1
	max_used = MAX(max_used,num_used)
	return
      endif
    enddo

    call system('banner "oh SHIT!"')
    write(0,*)'Tmpstorage: Attempt to allocate more than ',&
	 & maxtemps,' temps!'
    stop
  end function setelement
  
  subroutine freeelement(tempptr)
    implicit none
    real*8, dimension(:,:,:), pointer :: tempptr

    integer :: i

    do i = 1, maxtemps
      if(temps(i)%in_use==1) then
	if(associated(temps(i)%d,tempptr)) then
	  temps(i)%in_use=0
	  num_used = num_used - 1
	  return
	endif
      endif
    enddo

    call system('banner "oh SHIT!"')
    write(0,*)'Tmpstorage: Attempt to free element never allocated'
    stop
  end subroutine freeelement

end module tempstorage
