text:
! ================================================================================================ !
!  SixTrack Array Alloc Module
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  J. Molson, K.N. Sjobak, V.K. Berglyd Olsen, BE-ABP-HSS
!  Last Modified: 2018-12-17
! ================================================================================================ !

module mod_alloc

  use, intrinsic :: iso_fortran_env
  use crcoall

  implicit none

  ! The total quantity of memory allocated (bits)
  integer(kind=int64) :: allocated_bits = 0

  ! The maximum quantity of memory allocated (bits)
  integer(kind=int64) :: maximum_bits = 0
  integer(kind=int32) :: alloc_count  = 0

  integer(kind=int64), parameter :: mbyte = 1024*1024*8
  integer(kind=int64), parameter :: kbyte = 1024*8
  integer(kind=int64), parameter :: byte  = 8

  integer, parameter :: size_of_logical = 32 ! Unless specific flags are set, a logical is stored as an int32

  interface alloc

    module procedure alloc1di16  ! 1D int16
    module procedure alloc2di16  ! 2D int16
    module procedure alloc3di16  ! 3D int16

    module procedure alloc1di32  ! 1D int32
    module procedure alloc2di32  ! 2D int32
    module procedure alloc3di32  ! 3D int32

    module procedure alloc1di64  ! 1D int64
    module procedure alloc2di64  ! 2D int64
    module procedure alloc3di64  ! 3D int64

    module procedure alloc1dr32  ! 1D real32
    module procedure alloc2dr32  ! 2D real32
    module procedure alloc3dr32  ! 3D real32
    module procedure alloc4dr32  ! 4D real32

    module procedure alloc1dr64  ! 1D real64
    module procedure alloc2dr64  ! 2D real64
    module procedure alloc3dr64  ! 3D real64
    module procedure alloc4dr64  ! 4D real64

    module procedure alloc1dr128 ! 1D real128
    module procedure alloc2dr128 ! 2D real128
    module procedure alloc3dr128 ! 3D real128
    module procedure alloc4dr128 ! 4D real128

    module procedure alloc1dc    ! 1D character
    module procedure alloc2dc    ! 2D character

    module procedure alloc1ds    ! 1D string

    module procedure alloc1dl    ! 1D logical
    module procedure alloc2dl    ! 2D logical

  end interface alloc

  interface dealloc

    module procedure dealloc1di16  ! 1D int16
    module procedure dealloc2di16  ! 2D int16
    module procedure dealloc3di16  ! 3D int16

    module procedure dealloc1di32  ! 1D int32
    module procedure dealloc2di32  ! 2D int32
    module procedure dealloc3di32  ! 3D int32

    module procedure dealloc1di64  ! 1D int64
    module procedure dealloc2di64  ! 2D int64
    module procedure dealloc3di64  ! 3D int64

    module procedure dealloc1dr32  ! 1D real32
    module procedure dealloc2dr32  ! 2D real32
    module procedure dealloc3dr32  ! 3D real32
    module procedure dealloc4dr32  ! 4D real32

    module procedure dealloc1dr64  ! 1D real64
    module procedure dealloc2dr64  ! 2D real64
    module procedure dealloc3dr64  ! 3D real64
    module procedure dealloc4dr64  ! 4D real64

    module procedure dealloc1dr128 ! 1D real128
    module procedure dealloc2dr128 ! 2D real128
    module procedure dealloc3dr128 ! 3D real128
    module procedure dealloc4dr128 ! 4D real128

    module procedure dealloc1dc    ! 1D character
    module procedure dealloc2dc    ! 2D character

    module procedure dealloc1dl    ! 1D logical
    module procedure dealloc2dl    ! 2D logical

  end interface dealloc

  private :: alloc1di16,alloc2di16,alloc3di16
  private :: alloc1di32,alloc2di32,alloc3di32
  private :: alloc1di64,alloc2di64,alloc3di64
  private :: alloc1dr32,alloc2dr32,alloc3dr32,alloc4dr32
  private :: alloc1dr64,alloc2dr64,alloc3dr64,alloc4dr64
  private :: alloc1dr128,alloc2dr128,alloc3dr128,alloc4dr128
  private :: alloc1dc,alloc2dc
  private :: alloc1ds
  private :: alloc1dl,alloc2dl

  private :: dealloc1di16,dealloc2di16,dealloc3di16
  private :: dealloc1di32,dealloc2di32,dealloc3di32
  private :: dealloc1di64,dealloc2di64,dealloc3di64
  private :: dealloc1dr32,dealloc2dr32,dealloc3dr32,dealloc4dr32
  private :: dealloc1dr64,dealloc2dr64,dealloc3dr64,dealloc4dr64
  private :: dealloc1dr128,dealloc2dr128,dealloc3dr128,dealloc4dr128
  private :: dealloc1dc,dealloc2dc
  private :: dealloc1dl,dealloc2dl

contains

subroutine alloc_init

  return

end subroutine alloc_init

subroutine alloc_error(ename, error, requested_bits)

  implicit none

  character(len=*),    intent(in) :: ename          ! The variable name we are allocating
  integer(kind=int64), intent(in) :: requested_bits ! The requested allocation size
  integer,             intent(in) :: error          ! The error code

  write(lerr,"(a)")      "ALLOC> ERROR Memory allocation error for array "//ename
  write(lerr,"(a,i0,a)") "ALLOC>       Current allocation is:    ",allocated_bits," bits"
  write(lerr,"(a,i0,a)") "ALLOC>       Requested allocation is:  ",requested_bits," bits"
  write(lerr,"(a,i0)")   "ALLOC        Allocation error code is: ",error
  write(lerr,"(a)")      "ALLOC>       Exiting!"
  stop 1

end subroutine alloc_error

subroutine print_alloc(ename, type, requested_bits)

  use numerical_constants, only : one

  character(len=*), intent(in) :: ename, type
  integer(kind=int64) :: requested_bits

  alloc_count = alloc_count + 1
  if(allocated_bits > maximum_bits) then
    maximum_bits = allocated_bits
  end if

end subroutine print_alloc

! ================================================================================================ !
!  Allocation and Resizing
! ~~~~~~~~~~~~~~~~~~~~~~~~~
!  Inputs:
!    ename   :: The name of the variable being assigned,for error printing
!    eIdx    :: The new end index (size) that the array should be set to
!    initial :: The inital value for the array to be set to
!    input   :: The input variable to allocated
!    fIdx    :: The index of the first element
! ================================================================================================ !

! ================================================================================================ !
!  INTEGERS
! ================================================================================================ !

! 1D int16 Array
subroutine alloc1di16(input, eIdx, initial, ename, fIdxIn)

  implicit none

  integer(kind=int16), allocatable, intent(inout) :: input(:)
  integer,                          intent(in)    :: eIdx
  integer(kind=int16),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn

  integer(kind=int16), allocatable :: buffer(:)   ! Buffer array
  integer                          :: fIdx        ! First index
  integer                          :: oIdx        ! Old end index
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(int16)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(int16)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D int16",request)

end subroutine alloc1di16

! 2D int16 Array
subroutine alloc2di16(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  integer(kind=int16), allocatable, intent(inout) :: input(:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2
  integer(kind=int16),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2

  integer(kind=int16), allocatable :: buffer(:,:) ! Buffer array
  integer                          :: fIdx1,fIdx2 ! First indices
  integer                          :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(int16)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * storage_size(int16)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D int16",request)

end subroutine alloc2di16

! 3D int16 Array
subroutine alloc3di16(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  integer(kind=int16), allocatable, intent(inout) :: input(:,:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2,eIdx3
  integer(kind=int16),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  integer(kind=int16), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                          :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                          :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)              :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(int16)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(int16)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D int16",request)

end subroutine alloc3di16

! 1D int32 Array
subroutine alloc1di32(input, eIdx, initial, ename, fIdxIn)

  implicit none

  integer(kind=int32), allocatable, intent(inout) :: input(:)
  integer,                          intent(in)    :: eIdx
  integer(kind=int32),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn

  integer(kind=int32), allocatable :: buffer(:)   ! Buffer array
  integer                          :: fIdx        ! First index
  integer                          :: oIdx        ! Old end index
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(int32)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(int32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D int32",request)

end subroutine alloc1di32

! 2D int32 Array
subroutine alloc2di32(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  integer(kind=int32), allocatable, intent(inout) :: input(:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2
  integer(kind=int32),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2

  integer(kind=int32), allocatable :: buffer(:,:) ! Buffer array
  integer                          :: fIdx1,fIdx2 ! First indices
  integer                          :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(int32)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * storage_size(int32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D int32",request)

end subroutine alloc2di32

! 3D int32 Array
subroutine alloc3di32(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  integer(kind=int32), allocatable, intent(inout) :: input(:,:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2,eIdx3
  integer(kind=int32),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  integer(kind=int32), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                          :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                          :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)              :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(int32)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(int32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D int32",request)

end subroutine alloc3di32

! 1D int64 Array
subroutine alloc1di64(input, eIdx, initial, ename, fIdxIn)

  implicit none

  integer(kind=int64), allocatable, intent(inout) :: input(:)
  integer,                          intent(in)    :: eIdx
  integer(kind=int64),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn

  integer(kind=int64), allocatable :: buffer(:)   ! Buffer array
  integer                          :: fIdx        ! First index
  integer                          :: oIdx        ! Old end index
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(int64)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx, eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(int64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D int64",request)

end subroutine alloc1di64

! 2D int64 Array
subroutine alloc2di64(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  integer(kind=int64), allocatable, intent(inout) :: input(:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2
  integer(kind=int64),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2

  integer(kind=int64), allocatable :: buffer(:,:) ! Buffer array
  integer                          :: fIdx1,fIdx2 ! First indices
  integer                          :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)              :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(int64)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * storage_size(int64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D int64",request)

end subroutine alloc2di64

! 3D int64 Array
subroutine alloc3di64(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  integer(kind=int64), allocatable, intent(inout) :: input(:,:,:)
  integer,                          intent(in)    :: eIdx1,eIdx2,eIdx3
  integer(kind=int64),              intent(in)    :: initial
  character(len=*),                 intent(in)    :: ename
  integer,             optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  integer(kind=int64), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                          :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                          :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)              :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(int64)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(int64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D int64",request)

end subroutine alloc3di64

! ================================================================================================ !
!  REALS
! ================================================================================================ !

! 1D real32 Array
subroutine alloc1dr32(input, eIdx, initial, ename, fIdxIn)

  implicit none

  real(kind=real32), allocatable, intent(inout) :: input(:)
  integer,                        intent(in)    :: eIdx
  real(kind=real32),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn

  real(kind=real32), allocatable :: buffer(:)   ! Buffer array
  integer                        :: fIdx        ! First index
  integer                        :: oIdx        ! Old end index
  integer(kind=int64)            :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(real32)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(real32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D real32",request)

end subroutine alloc1dr32

! 2D real32 Array
subroutine alloc2dr32(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  real(kind=real32), allocatable, intent(inout) :: input(:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2
  real(kind=real32),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2

  real(kind=real32), allocatable :: buffer(:,:) ! Buffer array
  integer                        :: fIdx1,fIdx2 ! First indices
  integer                        :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)            :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(real32)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * storage_size(real32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D real32",request)

end subroutine alloc2dr32

! 3D real32 Array
subroutine alloc3dr32(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  real(kind=real32), allocatable, intent(inout) :: input(:,:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2,eIdx3
  real(kind=real32),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  real(kind=real32), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                        :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                        :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)            :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(real32)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(real32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D real32",request)

end subroutine alloc3dr32

! 4D real32 Array
subroutine alloc4dr32(input, eIdx1, eIdx2, eIdx3, eIdx4, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3, fIdxIn4)

  implicit none

  real(kind=real32), allocatable, intent(inout) :: input(:,:,:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2,eIdx3,eIdx4
  real(kind=real32),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3,fIdxIn4

  real(kind=real32), allocatable :: buffer(:,:,:,:)         ! Buffer array
  integer                        :: fIdx1,fIdx2,fIdx3,fIdx4 ! First indices
  integer                        :: oIdx1,oIdx2,oIdx3,oIdx4 ! Old end indices
  integer(kind=int64)            :: request                 ! Requested size addition

  integer i, j, k, l, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if
  if(present(fIdxIn4)) then
    fIdx4 = fIdxIn4
  else
    fIdx4 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1) * storage_size(real32)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            input(l,k,j,i) = initial
          end do
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    oIdx4   = size(input,4)+fIdx4-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)*(oIdx4-fIdx4+1)) * storage_size(real32)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            buffer(l,k,j,i) = initial
          end do
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3
    if(oIdx4 > eIdx4) oIdx4 = eIdx4

    do i=fIdx4,oIdx4
      do j=fIdx3,oIdx3
        do k=fIdx2,oIdx2
          do l=fIdx1,oIdx1
            buffer(l,k,j,i) = input(l,k,j,i)
          end do
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"4D real32",request)

end subroutine alloc4dr32

! 1D real64 Array
subroutine alloc1dr64(input, eIdx, initial, ename, fIdxIn)

  implicit none

  real(kind=real64), allocatable, intent(inout) :: input(:)
  integer,                        intent(in)    :: eIdx
  real(kind=real64),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn

  real(kind=real64), allocatable :: buffer(:)   ! Buffer array
  integer                        :: fIdx        ! First index
  integer                        :: oIdx        ! Old end index
  integer(kind=int64)            :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(real64)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx   = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(real64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D real64",request)

end subroutine alloc1dr64

! 2D real64 Array
subroutine alloc2dr64(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  real(kind=real64), allocatable, intent(inout) :: input(:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2
  real(kind=real64),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2

  real(kind=real64), allocatable :: buffer(:,:) ! Buffer array
  integer                        :: fIdx1,fIdx2 ! First indices
  integer                        :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)            :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(real64)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * storage_size(real64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D real64",request)

end subroutine alloc2dr64

! 3D real64 Array
subroutine alloc3dr64(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  real(kind=real64), allocatable, intent(inout) :: input(:,:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2,eIdx3
  real(kind=real64),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  real(kind=real64), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                        :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                        :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)            :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(real64)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(real64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D real64",request)

end subroutine alloc3dr64

! 4d real64 array
subroutine alloc4dr64(input, eIdx1, eIdx2, eIdx3, eIdx4, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3, fIdxIn4)

  implicit none

  real(kind=real64), allocatable, intent(inout) :: input(:,:,:,:)
  integer,                        intent(in)    :: eIdx1,eIdx2,eIdx3,eIdx4
  real(kind=real64),              intent(in)    :: initial
  character(len=*),               intent(in)    :: ename
  integer,           optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3,fIdxIn4

  real(kind=real64), allocatable :: buffer(:,:,:,:)         ! Buffer array
  integer                        :: fIdx1,fIdx2,fIdx3,fIdx4 ! First indices
  integer                        :: oIdx1,oIdx2,oIdx3,oIdx4 ! Old end indices
  integer(kind=int64)            :: request                 ! Requested size addition

  integer i, j, k, l, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if
  if(present(fIdxIn4)) then
    fIdx4 = fIdxIn4
  else
    fIdx4 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1) * storage_size(real64)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            input(l,k,j,i) = initial
          end do
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    oIdx4   = size(input,4)+fIdx4-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)*(oIdx4-fIdx4+1)) * storage_size(real64)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            buffer(l,k,j,i) = initial
          end do
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3
    if(oIdx4 > eIdx4) oIdx4 = eIdx4

    do i=fIdx4,oIdx4
      do j=fIdx3,oIdx3
        do k=fIdx2,oIdx2
          do l=fIdx1,oIdx1
            buffer(l,k,j,i) = input(l,k,j,i)
          end do
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"4D real64",request)

end subroutine alloc4dr64

! 1D real128 Array
subroutine alloc1dr128(input, eIdx, initial, ename, fIdxIn)

  implicit none

  real(kind=real128), allocatable, intent(inout) :: input(:)
  integer,                         intent(in)    :: eIdx
  real(kind=real128),              intent(in)    :: initial
  character(len=*),                intent(in)    :: ename
  integer,            optional,    intent(in)    :: fIdxIn

  real(kind=real128), allocatable :: buffer(:)   ! Buffer array
  integer                         :: fIdx        ! First index
  integer                         :: oIdx        ! Old end index
  integer(kind=int64)             :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * storage_size(real128)

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx   = size(input)+fIdx-1
    request = (eIdx-oIdx) * storage_size(real128)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D real128",request)

end subroutine alloc1dr128

! 2D real128 Array
subroutine alloc2dr128(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  real(kind=real128), allocatable, intent(inout) :: input(:,:)
  integer,                         intent(in)    :: eIdx1,eIdx2
  real(kind=real128),              intent(in)    :: initial
  character(len=*),                intent(in)    :: ename
  integer,            optional,    intent(in)    :: fIdxIn1,fIdxIn2

  real(kind=real128), allocatable :: buffer(:,:) ! Buffer array
  integer                         :: fIdx1,fIdx2 ! First indices
  integer                         :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)             :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * storage_size(real128)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1) * (eIdx2-fIdx2+1) - &
               (eIdx1-oIdx1)  * (eIdx2-oIdx2)) * storage_size(real128)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D real128",request)

end subroutine alloc2dr128

! 3D real128 Array
subroutine alloc3dr128(input, eIdx1, eIdx2, eIdx3, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3)

  implicit none

  real(kind=real128), allocatable, intent(inout) :: input(:,:,:)
  integer,                         intent(in)    :: eIdx1,eIdx2,eIdx3
  real(kind=real128),              intent(in)    :: initial
  character(len=*),                intent(in)    :: ename
  integer,            optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3

  real(kind=real128), allocatable :: buffer(:,:,:)     ! Buffer array
  integer                         :: fIdx1,fIdx2,fIdx3 ! First indices
  integer                         :: oIdx1,oIdx2,oIdx3 ! Old end indices
  integer(kind=int64)             :: request           ! Requested size addition

  integer i, j, k, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1) * storage_size(real128)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          input(k,j,i) = initial
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)) * storage_size(real128)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx3,eIdx3
      do j=fIdx2,eIdx2
        do k=fIdx1,eIdx1
          buffer(k,j,i) = initial
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3

    do i=fIdx3,oIdx3
      do j=fIdx2,oIdx2
        do k=fIdx1,oIdx1
          buffer(k,j,i) = input(k,j,i)
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"3D real128",request)

end subroutine alloc3dr128

! 4D real128 Array
subroutine alloc4dr128(input, eIdx1, eIdx2, eIdx3, eIdx4, initial, ename, fIdxIn1, fIdxIn2, fIdxIn3, fIdxIn4)

  implicit none

  real(kind=real128), allocatable, intent(inout) :: input(:,:,:,:)
  integer,                         intent(in)    :: eIdx1,eIdx2,eIdx3,eIdx4
  real(kind=real128),              intent(in)    :: initial
  character(len=*),                intent(in)    :: ename
  integer,            optional,    intent(in)    :: fIdxIn1,fIdxIn2,fIdxIn3,fIdxIn4

  real(kind=real128), allocatable :: buffer(:,:,:,:)         ! Buffer array
  integer                         :: fIdx1,fIdx2,fIdx3,fIdx4 ! First indices
  integer                         :: oIdx1,oIdx2,oIdx3,oIdx4 ! Old end indices
  integer(kind=int64)             :: request                 ! Requested size addition

  integer i, j, k, l, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if
  if(present(fIdxIn3)) then
    fIdx3 = fIdxIn3
  else
    fIdx3 = 1
  end if
  if(present(fIdxIn4)) then
    fIdx4 = fIdxIn4
  else
    fIdx4 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1) * storage_size(real128)

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            input(l,k,j,i) = initial
          end do
        end do
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    oIdx3   = size(input,3)+fIdx3-1
    oIdx4   = size(input,4)+fIdx4-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)*(eIdx3-fIdx3+1)*(eIdx4-fIdx4+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)*(oIdx3-fIdx3+1)*(oIdx4-fIdx4+1)) * storage_size(real128)

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2,fIdx3:eIdx3,fIdx4:eIdx4), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx4,eIdx4
      do j=fIdx3,eIdx3
        do k=fIdx2,eIdx2
          do l=fIdx1,eIdx1
            buffer(l,k,j,i) = initial
          end do
        end do
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2
    if(oIdx3 > eIdx3) oIdx3 = eIdx3
    if(oIdx4 > eIdx4) oIdx4 = eIdx4

    do i=fIdx4,oIdx4
      do j=fIdx3,oIdx3
        do k=fIdx2,oIdx2
          do l=fIdx1,oIdx1
            buffer(l,k,j,i) = input(l,k,j,i)
          end do
        end do
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"4D real128",request)

end subroutine alloc4dr128

! ================================================================================================ !
!  CHARACTERS
! ================================================================================================ !

! 1D character Array
subroutine alloc1dc(input, strlen, eIdx, initial, ename, fIdxIn)

  implicit none

  character(len=:), allocatable, intent(inout) :: input(:)
  integer,                       intent(in)    :: strlen
  integer,                       intent(in)    :: eIdx
  character(len=*),              intent(in)    :: initial
  character(len=*),              intent(in)    :: ename
  integer,          optional,    intent(in)    :: fIdxIn

  character(len=:), allocatable :: buffer(:)   ! Buffer array
  integer                       :: fIdx        ! First index
  integer                       :: oIdx        ! Old end index
  integer(kind=int64)           :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * strlen * CHARACTER_STORAGE_SIZE

    allocate(character(strlen) :: input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = (eIdx-oIdx) * strlen * CHARACTER_STORAGE_SIZE

    if(request == 0.0) then
      return
    end if

    allocate(character(strlen) :: buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D character",request)

end subroutine alloc1dc

! 2D character Array
subroutine alloc2dc(input, strlen, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  character(len=:), allocatable, intent(inout) :: input(:,:)
  integer,                       intent(in)    :: strlen
  integer,                       intent(in)    :: eIdx1,eIdx2
  character(len=*),              intent(in)    :: initial
  character(len=*),              intent(in)    :: ename
  integer,          optional,    intent(in)    :: fIdxIn1,fIdxIn2

  character(len=:), allocatable :: buffer(:,:) ! Buffer array
  integer                       :: fIdx1,fIdx2 ! First index
  integer                       :: oIdx1,oIdx2 ! Old end index
  integer(kind=int64)           :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * strlen * CHARACTER_STORAGE_SIZE

    allocate(character(strlen) :: input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1)  - &
               (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1)) * strlen * CHARACTER_STORAGE_SIZE

    if(request == 0.0) then
      return
    end if

    allocate(character(strlen) :: buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D character",request)

end subroutine alloc2dc

! ==================================================================== !
!  STRINGS
! ==================================================================== !

subroutine alloc1ds(input, eIdx, initial, ename, fIdxIn)

  use strings

  implicit none

  type(string), allocatable, intent(inout) :: input(:)
  integer,                   intent(in)    :: eIdx
  type(string),              intent(in)    :: initial
  character(len=*),          intent(in)    :: ename
  integer,      optional,    intent(in)    :: fIdxIn

  type(string), allocatable :: buffer(:)   ! Buffer array
  integer                   :: fIdx        ! First index
  integer                   :: oIdx        ! Old end index
  integer(kind=int64)       :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = eIdx-fIdx+1

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = eIdx-oIdx

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D string",request)

end subroutine alloc1ds

! ================================================================================================ !
!  LOGICALS
! ================================================================================================ !

! 1D logical Array
subroutine alloc1dl(input, eIdx, initial, ename, fIdxIn)

  implicit none

  logical, allocatable, intent(inout) :: input(:)
  integer,              intent(in)    :: eIdx
  logical,              intent(in)    :: initial
  character(len=*),     intent(in)    :: ename
  integer, optional,    intent(in)    :: fIdxIn

  logical, allocatable :: buffer(:)   ! Buffer array
  integer              :: fIdx        ! First index
  integer              :: oIdx        ! Old end index
  integer(kind=int64)  :: request     ! Requested size addition

  integer i, error

  if(present(fIdxIn)) then
    fIdx = fIdxIn
  else
    fIdx = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx-fIdx+1) * size_of_logical

    allocate(input(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      input(i) = initial
    end do

  else

    oIdx    = size(input)+fIdx-1
    request = eIdx-oIdx

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx:eIdx), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx,eIdx
      buffer(i) = initial
    end do

    if(oIdx > eIdx) oIdx = eIdx
    do i=fIdx,oIdx
      buffer(i) = input(i)
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"1D logical",request)

end subroutine alloc1dl

! 2D logical Array
subroutine alloc2dl(input, eIdx1, eIdx2, initial, ename, fIdxIn1, fIdxIn2)

  implicit none

  logical, allocatable, intent(inout) :: input(:,:)
  integer,              intent(in)    :: eIdx1,eIdx2
  logical,              intent(in)    :: initial
  character(len=*),     intent(in)    :: ename
  integer, optional,    intent(in)    :: fIdxIn1,fIdxIn2

  logical, allocatable :: buffer(:,:) ! Buffer array
  integer              :: fIdx1,fIdx2 ! First indices
  integer              :: oIdx1,oIdx2 ! Old end indices
  integer(kind=int64)  :: request     ! Requested size addition

  integer i, j, error

  if(present(fIdxIn1)) then
    fIdx1 = fIdxIn1
  else
    fIdx1 = 1
  end if
  if(present(fIdxIn2)) then
    fIdx2 = fIdxIn2
  else
    fIdx2 = 1
  end if

  if(.not. allocated(input)) then

    request = (eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) * size_of_logical

    allocate(input(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        input(j,i) = initial
      end do
    end do

  else

    oIdx1   = size(input,1)+fIdx1-1
    oIdx2   = size(input,2)+fIdx2-1
    request = ((eIdx1-fIdx1+1)*(eIdx2-fIdx2+1) - (oIdx1-fIdx1+1)*(oIdx2-fIdx2+1))

    if(request == 0.0) then
      return
    end if

    allocate(buffer(fIdx1:eIdx1,fIdx2:eIdx2), stat=error)
    if(error /= 0) call alloc_error(ename, error, request)

    do i=fIdx2,eIdx2
      do j=fIdx1,eIdx1
        buffer(j,i) = initial
      end do
    end do

    if(oIdx1 > eIdx1) oIdx1 = eIdx1
    if(oIdx2 > eIdx2) oIdx2 = eIdx2

    do i=fIdx2,oIdx2
      do j=fIdx1,oIdx1
        buffer(j,i) = input(j,i)
      end do
    end do

    call move_alloc(buffer,input)

  end if

  allocated_bits = allocated_bits + request
  call print_alloc(ename,"2D logical",request)

end subroutine alloc2dl

! ================================================================================================ !
!  DEALLOCTIONS
! ~~~~~~~~~~~~~~
!  These exist to keep track of the total allocated memory
! ================================================================================================ !

! 1D int16 Array
subroutine dealloc1di16(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int16), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(int16))
  deallocate(input)
end subroutine dealloc1di16

! 2D int16 Array
subroutine dealloc2di16(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int16), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(int16))
  deallocate(input)
end subroutine dealloc2di16

! 3D int16 Array
subroutine dealloc3di16(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int16), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(int16))
  deallocate(input)
end subroutine dealloc3di16

! 1D int32 Array
subroutine dealloc1di32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int32), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(int32))
  deallocate(input)
end subroutine dealloc1di32

! 2D int32 Array
subroutine dealloc2di32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int32), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(int32))
  deallocate(input)
end subroutine dealloc2di32

! 3D int32 Array
subroutine dealloc3di32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int32), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(int32))
  deallocate(input)
end subroutine dealloc3di32

! 1D int64 Array
subroutine dealloc1di64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int64), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(int64))
  deallocate(input)
end subroutine dealloc1di64

! 2D int64 Array
subroutine dealloc2di64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int64), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(int64))
  deallocate(input)
end subroutine dealloc2di64

! 3D int64 Array
subroutine dealloc3di64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  integer(kind=int64), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(int64))
  deallocate(input)
end subroutine dealloc3di64

! 1D real32 Array
subroutine dealloc1dr32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real32), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(real32))
  deallocate(input)
end subroutine dealloc1dr32

! 2D real32 Array
subroutine dealloc2dr32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real32), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(real32))
  deallocate(input)
end subroutine dealloc2dr32

! 3D real32 Array
subroutine dealloc3dr32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real32), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(real32))
  deallocate(input)
end subroutine dealloc3dr32

! 4D real32 Array
subroutine dealloc4dr32(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real32), allocatable, intent(inout) :: input(:,:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*size(input,4)*storage_size(real32))
  deallocate(input)
end subroutine dealloc4dr32

! 1D real64 Array
subroutine dealloc1dr64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real64), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(real64))
  deallocate(input)
end subroutine dealloc1dr64

! 2D real64 Array
subroutine dealloc2dr64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real64), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(real64))
  deallocate(input)
end subroutine dealloc2dr64

! 3D real64 Array
subroutine dealloc3dr64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real64), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(real64))
  deallocate(input)
end subroutine dealloc3dr64

! 4D real64 Array
subroutine dealloc4dr64(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real64), allocatable, intent(inout) :: input(:,:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*size(input,4)*storage_size(real64))
  deallocate(input)
end subroutine dealloc4dr64

! 1D real128 Array
subroutine dealloc1dr128(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real128), allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*storage_size(real128))
  deallocate(input)
end subroutine dealloc1dr128

! 2D real128 Array
subroutine dealloc2dr128(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real128), allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*storage_size(real128))
  deallocate(input)
end subroutine dealloc2dr128

! 3D real128 Array
subroutine dealloc3dr128(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real128), allocatable, intent(inout) :: input(:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*storage_size(real128))
  deallocate(input)
end subroutine dealloc3dr128

! 4D real128 Array
subroutine dealloc4dr128(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  real(kind=real128), allocatable, intent(inout) :: input(:,:,:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*size(input,3)*size(input,4)*storage_size(real128))
  deallocate(input)
end subroutine dealloc4dr128

! 1D character Array
subroutine dealloc1dc(input, strlen, ename)
  implicit none
  character(len=*), intent(in) :: ename
  character(len=:), allocatable, intent(inout) :: input(:)
  integer, intent(in) :: strlen

  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  !TODO: Check what size(input) actually does
  allocated_bits = allocated_bits - (size(input,1)*strlen * CHARACTER_STORAGE_SIZE)
  deallocate(input)
end subroutine dealloc1dc

! 2D character Array
subroutine dealloc2dc(input, strlen, ename)
  implicit none
  character(len=*), intent(in) :: ename
  character(len=:), allocatable, intent(inout) :: input(:,:)
  integer, intent(in) :: strlen

  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  !TODO: Check what size(input) actually does
  allocated_bits = allocated_bits - (size(input,1)*size(input,2)*strlen * CHARACTER_STORAGE_SIZE)
  deallocate(input)
end subroutine dealloc2dc

! 1D logical Array
subroutine dealloc1dl(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  logical, allocatable, intent(inout) :: input(:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - size(input,1)*size_of_logical
  deallocate(input)
end subroutine dealloc1dl

! 2D logical Array
subroutine dealloc2dl(input, ename)
  implicit none
  character(len=*), intent(in) :: ename
  logical, allocatable, intent(inout) :: input(:,:)
  !Check that we are already allocated
  if(.not. allocated(input)) then
    write(lerr,"(a)") "ALLOC> ERROR Trying to deallocate a NULL pointer: '"//ename//"'"
    stop
  end if
  allocated_bits = allocated_bits - size(input,1)*size(input,2)*size_of_logical
  deallocate(input)
end subroutine dealloc2dl

end module mod_alloc

max_stars_repo_name: xsuite/xcol

max_issues_repo_name: xsuite/xcol

max_forks_repo_name: xsuite/xcol


text:
      subroutine runCP_0000(Gr,S0000,N0)
      implicit none
      include 'TRconstants.f'
      include 'pvCnames.f'
      include 'pvCv.f'
      include 'Carraydef.f'
      include 'Carrays.f'
      integer ep,N0,m,n,np
      parameter(np=2)
      double precision Gr(np,np)
      double complex S0000(-2:0),bit,pole

      do ep=-2,0
      bit=czip
      do n=1,np
      do m=1,np
      bit=bit+Gr(n,m)*Cv(czzii(z2(n,m))+N0,ep)
      enddo
      enddo

      pole=czip
      if (ep .gt. -2) pole=4d0*Cv(cc0000+N0,ep-1)

c--- note: we have simplified the recursion relation by using DD (5.9)
c--- so that S0000 = 2B00(0)+2*m0sq*C00
      Cv(cc0000+N0,ep)=
     . (pole
     . +S0000(ep)
     . -bit)/12d0
      enddo

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
  function mnmacro (molObj, iz, in)

! ======================================================================
!
!   This function calculates the Moller-Nix macroscopic mass excess
!   for (spherical) nuclei lying outside the tabulated range of Z and N.
!   The inputs are the proton and neutron numbers.
!   Uses formula (62) from At. Data Nucl. data Tables 59, 185 (1995).
!   Uses FRLDM, not FRDM.
!
!   Called by: MOLNIX
!
!    Written by A. J. Sierk  LANL  T-2  November, 1997.
!    "Last" change: 13-AUG-2003 by NVMokhov
!    Edited by A. J. Sierk, LANL T-16, October, 2003.
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection).
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: int32, real64
    use molnixParams, only: zro, hlf, one, two, thr, fiv, &
         & thrd, twthrd, pi, ato3rd

    implicit none
    class(Molnix),  intent(inout) :: molObj
    integer(int32), intent(in   ) :: in
    integer(int32), intent(in   ) :: iz
    real(real64)                  :: mnmacro

    real(real64) :: a, ai, akf, athrd, b3, bi, deln, delnp, delp, &
         & fkprp, temp, un, wt, xx, yc, yc2, ys, z
    logical      :: nodd, podd

! ======================================================================

    real(real64), parameter :: &
         &  emp = 7.289034d0,  emn = 8.071431d0, &
         &  esq = 1.4399764d0, ael = 1.433d-5, &
         &   rp = 0.80d0,       r0 = 1.16d0, &
         & aype = 0.68d0,       ay = 0.70d0, &
         & rmac = 4.80d0,        h = 6.6d0, &
         &    w = 30.d0,        av = 16.00126d0, &
         &  akv = 1.92240d0,    as = 21.18466d0, &
         &  aks = 2.345d0,      a0 = 2.615d0, &
         &   ca = 0.10289d0

    real(real64), parameter :: d1 = 3.020833333d0
    real(real64), parameter :: d2 = 0.113541666667d0
    real(real64), parameter :: d3 = 0.00126240079d0

    real(real64), parameter :: c1 = thr*esq/(fiv*r0)
    real(real64), parameter :: c4 = 1.25d0*c1*(1.5d0/pi)**twthrd
    real(real64), parameter :: c5 = (2.25d0*pi)**thrd/r0

! ======================================================================

    un = dble(in)
    z = dble(iz)
    a = z + un
    athrd = ato3rd(in + iz)
    if (a < div0Lim .and. a > -div0Lim) then
       a = div0Lim
       write(molObj%io%message,1000) "292, 323"
       call molObj%io%print(4, 3, molObj%io%message)
    end if
    ai = (un - z)/a
    deln = zro
    delnp = zro
    nodd = .false.
    podd = .false.
    if (2*(in/2).ne.in) then
       temp = ato3rd(in)
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(molObj%io%message,1000) "304"
          call molObj%io%print(4, 3, molObj%io%message)
       end if
       deln = rmac/temp
       nodd = .true.
    endif
    delp = zro
    if (2*(iz/2).ne.iz) then
       temp = ato3rd(iz)
       if (temp < div0Lim .and. temp > -div0Lim) then
          temp = div0Lim
          write(molObj%io%message,1000) "314"
          call molObj%io%print(4, 3, molObj%io%message)
       end if
       delp = rmac/temp
       podd = .true.
    endif
    wt = zro
    if (nodd .and. podd) then
       if (athrd < div0Lim .and. athrd > -div0Lim) then
          athrd = div0Lim
          write(molObj%io%message,1000) "322, 325"
          call molObj%io%print(4, 3, molObj%io%message)
       end if
       delnp = h/athrd**2
       if (in == iz) wt = one/a
    endif
    akf = c5*ato3rd(iz)/athrd
    xx = akf*rp
    fkprp = -0.125d0*esq*rp**2*(d1 - d2*xx**2 + d3*xx**4)/r0**3
    ys = r0*athrd/aype
    yc = ys*aype/ay
    if (ys < div0Lim .and. ys > -div0Lim) then
       ys = div0Lim
    end if
    if (yc < div0Lim .and. yc > -div0Lim) then
       yc = div0Lim
       write(molObj%io%message,1000) "344"
       call molObj%io%print(4, 3, molObj%io%message)
    end if
    yc2 = yc*yc

    bi = one - thr/ys**2 + (one + ys)*( two + thr * (one + one/ys)/ys) &
         & * exp(-two*ys)
    b3 = one - fiv*(one - 1.875d0*(one - 1.4d0/yc2)/yc - 0.75d0 * &
         & (one + 4.5d0/yc + 7.0d0*(one + hlf/yc)/yc2)*exp(-two*yc)) / &
         & yc2

    mnmacro = emp*z + emn*un - av*a*(one - akv*ai**2) + as*bi*  &
         (one - aks*ai**2)*athrd**2+ a0 + c1*b3*z**2/athrd -  &
         c4*ato3rd(iz)**4/athrd + fkprp*z**2/a - ca*ai*a +  &
         w*(abs(ai) + wt) + delp + deln - delnp - ael*z**2.39d0

    return
! ======================================================================
1000 format("Divide by zero error prevented in 'mnmacro.f90', line(s) ", &
          & A, ".")
! ======================================================================
  end function mnmacro

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


text:
      module rrlw_kg15

      use parkind ,only : im => kind_im, rb => kind_rb

      implicit none
      save

!-----------------------------------------------------------------
! rrtmg_lw ORIGINAL abs. coefficients for interval 15
! band 15:  2380-2600 cm-1 (low - n2o,co2; high - nothing)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
! Revised: MJIacono, AER, aug2008
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
!fracrefao: real
! kao     : real
! kao_mn2 : real
! selfrefo: real
! forrefo : real
!-----------------------------------------------------------------

      integer(kind=im), parameter :: no15 = 16

      real(kind=rb) :: fracrefao(no15,9)
      real(kind=rb) :: kao(9,5,13,no15)
      real(kind=rb) :: kao_mn2(9,19,no15)
      real(kind=rb) :: selfrefo(10,no15)
      real(kind=rb) :: forrefo(4,no15)

!-----------------------------------------------------------------
! rrtmg_lw COMBINED abs. coefficients for interval 15
! band 15:  2380-2600 cm-1 (low - n2o,co2; high - nothing)
!
! Initial version:  JJMorcrette, ECMWF, jul1998
! Revised: MJIacono, AER, jun2006
! Revised: MJIacono, AER, aug2008
!-----------------------------------------------------------------
!
!  name     type     purpose
!  ----   : ----   : ---------------------------------------------
!fracrefa : real
! ka      : real
! ka_mn2  : real
! selfref : real
! forref  : real
!
! absa    : real
!-----------------------------------------------------------------

      integer(kind=im), parameter :: ng15 = 2

      real(kind=rb) :: fracrefa(ng15,9)
      real(kind=rb) :: ka(9,5,13,ng15) ,absa(585,ng15)
      real(kind=rb) :: ka_mn2(9,19,ng15)
      real(kind=rb) :: selfref(10,ng15)
      real(kind=rb) :: forref(4,ng15)

      equivalence (ka(1,1,1,1),absa(1,1))

      end module rrlw_kg15

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
! -*- Mode: Fortran; -*-
!
!  (C) 2014 by Argonne National Laboratory.
!  See COPYRIGHT in top-level directory.
!
program main
use mpi_f08

implicit none
!include 'mpif.h'

! Fortran 2008 equivalent of src/mpi/romio/test/coll_test.c

integer FILESIZE
parameter (FILESIZE=32*32*32*4)

! A 32^3 array. For other array sizes, change FILESIZE above and
! array_of_gsizes below.

! Uses collective I/O. Writes a 3D block-distributed array to a file
! corresponding to the global array in row-major (C) order, reads it
! back, and checks that the data read is correct.

! Note that the file access pattern is noncontiguous.

integer i, ndims, array_of_gsizes(3)
integer order, intsize, nprocs, j, array_of_distribs(3)
integer array_of_dargs(3), array_of_psizes(3)
integer readbuf(FILESIZE), writebuf(FILESIZE), bufcount
integer mynod, tmpbuf(FILESIZE), array_size, argc
integer ierr
character*256 str   ! used to store the filename
integer errs, toterrs
integer(MPI_OFFSET_KIND) :: disp

type(MPI_Datatype) ::  newtype
type(MPI_Status) :: status
type(MPI_Request) :: request
type(MPI_File) :: fh

errs = 0
str = "iotest.txt"

call MTEST_INIT(ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, mynod, ierr)

! create the distributed array filetype

ndims = 3
order = MPI_ORDER_FORTRAN

array_of_gsizes(1) = 32
array_of_gsizes(2) = 32
array_of_gsizes(3) = 32

array_of_distribs(1) = MPI_DISTRIBUTE_BLOCK
array_of_distribs(2) = MPI_DISTRIBUTE_BLOCK
array_of_distribs(3) = MPI_DISTRIBUTE_BLOCK

array_of_dargs(1) = MPI_DISTRIBUTE_DFLT_DARG
array_of_dargs(2) = MPI_DISTRIBUTE_DFLT_DARG
array_of_dargs(3) = MPI_DISTRIBUTE_DFLT_DARG

do i=1, ndims
   array_of_psizes(i) = 0
end do

call MPI_DIMS_CREATE(nprocs, ndims, array_of_psizes, ierr)

call MPI_TYPE_CREATE_DARRAY(nprocs, mynod, ndims,   &
     array_of_gsizes, array_of_distribs, array_of_dargs,  &
     array_of_psizes, order, MPI_INTEGER, newtype, ierr)

call MPI_TYPE_COMMIT(newtype, ierr)

! initialize writebuf

call MPI_TYPE_SIZE(newtype, bufcount, ierr)
call MPI_TYPE_SIZE(MPI_INTEGER, intsize, ierr)
bufcount = bufcount/intsize
do i=1, bufcount
   writebuf(i) = 1
end do

do i=1, FILESIZE
   tmpbuf(i) = 0
end do

call MPI_IRECV(tmpbuf, 1, newtype, mynod, 10, MPI_COMM_WORLD, request, ierr)
call MPI_SEND(writebuf, bufcount, MPI_INTEGER, mynod, 10, MPI_COMM_WORLD, ierr)
call MPI_WAIT(request, status, ierr)

j = 1
array_size = array_of_gsizes(1) * array_of_gsizes(2) * array_of_gsizes(3)
do i=1, array_size
   if (tmpbuf(i) .ne. 0) then
      writebuf(j) = i
      j = j + 1
   end if
end do

! end of initialization

! write the array to the file

call MPI_FILE_OPEN(MPI_COMM_WORLD, str, MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

disp = 0
call MPI_FILE_SET_VIEW(fh, disp, MPI_INTEGER, newtype, "native", MPI_INFO_NULL, ierr)
call MPI_FILE_IWRITE_ALL(fh, writebuf, bufcount, MPI_INTEGER, request, ierr)
call MPI_WAIT(request, status, ierr)
call MPI_FILE_CLOSE(fh, ierr)

!now read it back

call MPI_FILE_OPEN(MPI_COMM_WORLD, str, MPI_MODE_CREATE+MPI_MODE_RDWR, MPI_INFO_NULL, fh, ierr)

call MPI_FILE_SET_VIEW(fh, disp, MPI_INTEGER, newtype, "native", MPI_INFO_NULL, ierr)
call MPI_FILE_IREAD_ALL(fh, readbuf, bufcount, MPI_INTEGER, request, ierr)
call MPI_WAIT(request, status, ierr)
call MPI_FILE_CLOSE(fh, ierr)

! check the data read
do i=1, bufcount
   if (readbuf(i) .ne. writebuf(i)) then
      errs = errs + 1
       print *, 'Node ', mynod, '  readbuf ', readbuf(i),  &
          '  writebuf ', writebuf(i), '  i', i
   end if
end do

call MPI_TYPE_FREE(newtype, ierr)

if (mynod .eq. 0) then
    call MPI_FILE_DELETE(str, MPI_INFO_NULL, ierr)
    if (ierr .ne. MPI_SUCCESS) then
        errs = errs + 1
    endif
endif

call MTEST_FINALIZE(errs)

stop
end

max_stars_repo_name: ucd-plse/mpi-error-prop

max_issues_repo_name: ucd-plse/mpi-error-prop

max_forks_repo_name: ucd-plse/mpi-error-prop


text:
*----------------------------------------------------------------------*
      subroutine mod_op_for_ge_test(label_mel,
     &     iRdef,norb,icase,icaseF,
     &     op_info,str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
*     modify operator elements (most likely the Hamiltonian)
*     for generalized extensivity test
*     iRdef(norb) defines system R, the rest is assumed system S
*     icaseF: 0 = Fock is part of E0, else = Fock is also split into R/S
*     case 1: set R and S to zero
*     case 2: only R
*     case 3: only S
*     case 4: R and S (coupling elements are 0 nontheless)
*----------------------------------------------------------------------*

      implicit none

      integer, parameter ::
     &     ntest =   00

      include 'stdunit.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_orbinf.h'
      include 'par_opnames_gen.h'
      include 'mdef_operator_info.h'
      include 'def_strmapinf.h'

      character(*), intent(in) ::
     &     label_mel
      type(operator_info), intent(in) ::
     &     op_info
      integer, intent(in) ::
     &     norb, iRdef(norb), icase, icaseF
      type(strinf), intent(in) ::
     &     str_info
      type(strmapinf) ::
     &     strmap_info
      type(orbinf), intent(in) ::
     &     orb_info

      integer ::
     &     ipri, mode, scaling, idx_mel
      real(8) ::
     &     xdum
      logical ::
     &     list_exists
      type(me_list), pointer ::
     &     mel_target

      integer, external ::
     &     idx_mel_list

      if (ntest.ge.10) then
        call write_title(lulog,wst_dbg_subr,'mod_op_for_ge_test')
        write(lulog,*) 'Rdef:  ',iRdef(1:norb)
        write(lulog,*) 'case:  ',icase
        write(lulog,*) 'icaseF:',icaseF
      end if

      if (icase.lt.1.and.icase.gt.4)
     &     call quit(1,'mod_op_for_ge_test',
     &     'illegal value for icase')

      idx_mel = idx_mel_list(label_mel,op_info)
      if (idx_mel.lt.0)
     &     call quit(1,'mod_op_for_ge_test','Label not on list: "'//
     &     trim(label_mel)//'"')

      mel_target => op_info%mel_arr(idx_mel)%mel

      call touch_file_rec(mel_target%fhand)

      call getest_mod_mel(mel_target,
     &     iRdef,norb,icase,icaseF,
     &     str_info,orb_info)

      if (ntest.ge.10.and.(.not.mel_target%op%formal)) then
        write(lulog,*)
        write(lulog,*) 'modified list: ',trim(mel_target%label)
        if (ntest.ge.10) ipri = 1
        if (ntest.ge.50) ipri = 2
        if (ntest.ge.100) ipri = 3
        if (ntest.ge.500) ipri = 4
        if (ntest.ge.1000) ipri = 5
        call wrt_mel_file(lulog,ipri,mel_target,
     &       1,mel_target%op%n_occ_cls,
     &       str_info,orb_info)
      end if

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
submodule(forlab_math) forlab_math_is_close

    use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
    implicit none

contains

    elemental module logical function is_close_rsp(a, b, rel_tol, abs_tol, equal_nan) result(close)
        real(sp), intent(in) :: a, b
        real(sp), intent(in), optional :: rel_tol, abs_tol
        logical, intent(in), optional :: equal_nan
        logical :: equal_nan_

        equal_nan_ = optval(equal_nan, .false.)

        if (ieee_is_nan(a) .or. ieee_is_nan(b)) then
            close = merge(.true., .false., equal_nan_ .and. ieee_is_nan(a) .and. ieee_is_nan(b))
        else
            close = abs(a - b) <= max(abs(optval(rel_tol, 1.0e-9_sp)*max(abs(a), abs(b))), &
                                      abs(optval(abs_tol, 0.0_sp)))
        end if

    end function is_close_rsp
    elemental module logical function is_close_rdp(a, b, rel_tol, abs_tol, equal_nan) result(close)
        real(dp), intent(in) :: a, b
        real(dp), intent(in), optional :: rel_tol, abs_tol
        logical, intent(in), optional :: equal_nan
        logical :: equal_nan_

        equal_nan_ = optval(equal_nan, .false.)

        if (ieee_is_nan(a) .or. ieee_is_nan(b)) then
            close = merge(.true., .false., equal_nan_ .and. ieee_is_nan(a) .and. ieee_is_nan(b))
        else
            close = abs(a - b) <= max(abs(optval(rel_tol, 1.0e-9_dp)*max(abs(a), abs(b))), &
                                      abs(optval(abs_tol, 0.0_dp)))
        end if

    end function is_close_rdp

    elemental module logical function is_close_csp(a, b, rel_tol, abs_tol, equal_nan) result(close)
        complex(sp), intent(in) :: a, b
        real(sp), intent(in), optional :: rel_tol, abs_tol
        logical, intent(in), optional :: equal_nan

        close = is_close_rsp(a%re, b%re, rel_tol, abs_tol, equal_nan) .and. &
                is_close_rsp(a%im, b%im, rel_tol, abs_tol, equal_nan)

    end function is_close_csp
    elemental module logical function is_close_cdp(a, b, rel_tol, abs_tol, equal_nan) result(close)
        complex(dp), intent(in) :: a, b
        real(dp), intent(in), optional :: rel_tol, abs_tol
        logical, intent(in), optional :: equal_nan

        close = is_close_rdp(a%re, b%re, rel_tol, abs_tol, equal_nan) .and. &
                is_close_rdp(a%im, b%im, rel_tol, abs_tol, equal_nan)

    end function is_close_cdp

end submodule forlab_math_is_close

max_stars_repo_name: fortran-fans/forlab

max_issues_repo_name: fortran-fans/forlab

max_forks_repo_name: fortran-fans/forlab


text:
PROGRAM msckb

! Code converted using TO_F90 by Alan Miller
! Date: 2004-06-28  Time: 12:58:12

!     ============================================================
!     Purpose: This program computes the expansion coefficients
!     of the prolate and oblate spheroidal functions,
!     c2k, using subroutine SCKB
!     Input :  m  --- Mode parameter
!     n  --- Mode parameter
!     c  --- Spheroidal parameter
!     cv --- Characteristic value
!     KD --- Function code
!     KD=1 for prolate; KD=-1 for oblate
!     Output:  CK(k) --- Expansion coefficients ck;
!     CK(1), CK(2), ... correspond to
!     c0, c2, ...
!     Example: Compute the first 13 expansion coefficients C2k for
!     KD= 1, m=2, n=3, c=3.0 and cv=14.8277782138; and
!     KD=-1, m=2, n=3, c=3.0 and cv=8.80939392077

!     Coefficients of Prolate and oblate functions

!     k         C2k(c)              C2k(-ic)
!     ---------------------------------------------
!     0     .9173213327D+01     .2489664942D+02
!     1     .4718258929D+01    -.1205287032D+02
!     2     .9841212916D+00     .2410564082D+01
!     3     .1151870224D+00    -.2735821590D+00
!     4     .8733916403D-02     .2026057157D-01
!     5     .4663888254D-03    -.1061946315D-02
!     6     .1853910398D-04     .4158091152D-04
!     7     .5708084895D-06    -.1264400411D-05
!     8     .1402786472D-07     .3074963448D-07
!     9     .2817194508D-09    -.6120579463D-09
!     10     .4712094447D-11     .1015900041D-10
!     11     .6667838485D-13    -.1427953361D-12
!     12     .8087995432D-15     .1721924955D-14
!     ============================================================

IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION ck(200),df(200),eg(200)
WRITE(*,*)'Please KD, m, n and c '
!     READ(*,*)KD,M,N,C
kd=1
m=2
n=3
c=3.0
CALL segv(m,n,c,kd,cv,eg)
WRITE(*,30)kd,m,n,c,cv
CALL sdmn(m,n,c,cv,kd,df)
CALL sckb(m,n,c,df,ck)
WRITE(*,*)
IF (kd == 1) THEN
  WRITE(*,*)'Coefficients of Prolate function'
  WRITE(*,*)
  WRITE(*,*)'   k            C2k(c)'
ELSE
  WRITE(*,*)'Coefficients of Oblate function'
  WRITE(*,*)
  WRITE(*,*)'   k           C2k(-ic)'
END IF
WRITE(*,*)'----------------------------'
nm=25+INT((n-m)/2+c)
DO  k=1,nm
  WRITE(*,20)k-1,ck(k)
END DO
20   FORMAT(2X,i3,4X,d18.10)
30   FORMAT(1X,3HKD=,i3,',  ',2HM=,i3,',  ',2HN=,i3,  &
    ',  ',2HC=,f5.1,',  ',4HCV =,f18.10)
END PROGRAM msckb

SUBROUTINE sckb(m,n,c,df,ck)

!     ======================================================
!     Purpose: Compute the expansion coefficients of the
!     prolate and oblate spheroidal functions, c2k
!     Input :  m  --- Mode parameter
!     n  --- Mode parameter
!     c  --- Spheroidal parameter
!     DF(k) --- Expansion coefficients dk
!     Output:  CK(k) --- Expansion coefficients ck;
!     CK(1), CK(2), ... correspond to
!     c0, c2, ...
!     ======================================================

INTEGER, INTENT(IN)                      :: m
INTEGER, INTENT(IN)                      :: n
DOUBLE PRECISION, INTENT(IN OUT)         :: c
DOUBLE PRECISION, INTENT(IN)             :: df(200)
DOUBLE PRECISION, INTENT(OUT)            :: ck(200)
IMPLICIT DOUBLE PRECISION (a-h,o-z)

sw=0.0
IF (c <= 1.0D-10) c=1.0D-10
nm=25+INT(0.5*(n-m)+c)
ip=1
IF (n-m == 2*INT((n-m)/2)) ip=0
fac=-0.5D0**m
reg=1.0D0
IF (m+nm > 80) reg=1.0D-200
DO  k=0,nm-1
  fac=-fac
  i1=2*k+ip+1
  r=reg
  DO  i=i1,i1+2*m-1
    r=r*i
  END DO
  i2=k+m+ip
  DO  i=i2,i2+k-1
    r=r*(i+0.5D0)
  END DO
  sum=r*df(k+1)
  DO  i=k+1,nm
    d1=2.0D0*i+ip
    d2=2.0D0*m+d1
    d3=i+m+ip-0.5D0
    r=r*d2*(d2-1.0D0)*i*(d3+k)/(d1*(d1-1.0D0)*(i-k)*d3)
    sum=sum+r*df(i+1)
    IF (DABS(sw-sum) < DABS(sum)*1.0D-14) EXIT
    sw=sum
  END DO
  25    r1=reg
  DO  i=2,m+k
    r1=r1*i
  END DO
  ck(k+1)=fac*sum/r1
END DO
RETURN
END SUBROUTINE sckb

SUBROUTINE sdmn(m,n,c,cv,kd,df)

!     =====================================================
!     Purpose: Compute the expansion coefficients of the
!     prolate and oblate spheroidal functions
!     Input :  m  --- Mode parameter
!     n  --- Mode parameter
!     c  --- Spheroidal parameter
!     cv --- Characteristic value
!     KD --- Function code
!     KD=1 for prolate; KD=-1 for oblate
!     Output:  DF(k) --- Expansion coefficients dk;
!     DF(1), DF(2), ... correspond to
!     d0, d2, ... for even n-m and d1,
!     d3, ... for odd n-m
!     =====================================================

INTEGER, INTENT(IN)                      :: m
INTEGER, INTENT(IN)                      :: n
DOUBLE PRECISION, INTENT(IN)             :: c
DOUBLE PRECISION, INTENT(IN OUT)         :: cv
INTEGER, INTENT(IN)                      :: kd
DOUBLE PRECISION, INTENT(OUT)            :: df(200)
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION a(200),d(200),g(200)

sw=0.0
fl=0.0
nm=25+INT(0.5*(n-m)+c)
IF (c < 1.0D-10) THEN
  DO  i=1,nm
    df(i)=0D0
  END DO
  df((n-m)/2+1)=1.0D0
  RETURN
END IF
cs=c*c*kd
ip=1
IF (n-m == 2*INT((n-m)/2)) ip=0
DO  i=1,nm+2
  IF (ip == 0) k=2*(i-1)
  IF (ip == 1) k=2*i-1
  dk0=m+k
  dk1=m+k+1
  dk2=2*(m+k)
  d2k=2*m+k
  a(i)=(d2k+2.0)*(d2k+1.0)/((dk2+3.0)*(dk2+5.0))*cs
  d(i)=dk0*dk1+(2.0*dk0*dk1-2.0*m*m-1.0)/((dk2-1.0) *(dk2+3.0))*cs
  g(i)=k*(k-1.0)/((dk2-3.0)*(dk2-1.0))*cs
END DO
fs=1.0D0
f1=0.0D0
f0=1.0D-100
kb=0
df(nm+1)=0.0D0
DO  k=nm,1,-1
  f=-((d(k+1)-cv)*f0+a(k+1)*f1)/g(k+1)
  IF (DABS(f) > DABS(df(k+1))) THEN
    df(k)=f
    f1=f0
    f0=f
    IF (DABS(f) > 1.0D+100) THEN
      DO  k1=k,nm
        df(k1)=df(k1)*1.0D-100
      END DO
      f1=f1*1.0D-100
      f0=f0*1.0D-100
    END IF
  ELSE
    kb=k
    fl=df(k+1)
    f1=1.0D-100
    f2=-(d(1)-cv)/a(1)*f1
    df(1)=f1
    IF (kb == 1) THEN
      fs=f2
    ELSE IF (kb == 2) THEN
      df(2)=f2
      fs=-((d(2)-cv)*f2+g(2)*f1)/a(2)
    ELSE
      df(2)=f2
      DO  j=3,kb+1
        f=-((d(j-1)-cv)*f2+g(j-1)*f1)/a(j-1)
        IF (j <= kb) df(j)=f
        IF (DABS(f) > 1.0D+100) THEN
          DO  k1=1,j
            df(k1)=df(k1)*1.0D-100
          END DO
          f=f*1.0D-100
          f2=f2*1.0D-100
        END IF
        f1=f2
        f2=f
      END DO
      fs=f
    END IF
    EXIT
  END IF
END DO
35   su1=0.0D0
r1=1.0D0
DO  j=m+ip+1,2*(m+ip)
  r1=r1*j
END DO
su1=df(1)*r1
DO  k=2,kb
  r1=-r1*(k+m+ip-1.5D0)/(k-1.0D0)
  su1=su1+r1*df(k)
END DO
su2=0.0D0
DO  k=kb+1,nm
  IF (k /= 1) r1=-r1*(k+m+ip-1.5D0)/(k-1.0D0)
  su2=su2+r1*df(k)
  IF (DABS(sw-su2) < DABS(su2)*1.0D-14) EXIT
  sw=su2
END DO
55   r3=1.0D0
DO  j=1,(m+n+ip)/2
  r3=r3*(j+0.5D0*(n+m+ip))
END DO
r4=1.0D0
DO  j=1,(n-m-ip)/2
  r4=-4.0D0*r4*j
END DO
s0=r3/(fl*(su1/fs)+su2)/r4
DO  k=1,kb
  df(k)=fl/fs*s0*df(k)
END DO
DO  k=kb+1,nm
  df(k)=s0*df(k)
END DO
RETURN
END SUBROUTINE sdmn

SUBROUTINE segv(m,n,c,kd,cv,eg)

!     =========================================================
!     Purpose: Compute the characteristic values of spheroidal
!     wave functions
!     Input :  m  --- Mode parameter
!     n  --- Mode parameter
!     c  --- Spheroidal parameter
!     KD --- Function code
!     KD=1 for Prolate; KD=-1 for Oblate
!     Output:  CV --- Characteristic value for given m, n and c
!     EG(L) --- Characteristic value for mode m and n'
!     ( L = n' - m + 1 )
!     =========================================================

INTEGER, INTENT(IN)                      :: m
INTEGER, INTENT(IN)                      :: n
DOUBLE PRECISION, INTENT(IN)             :: c
INTEGER, INTENT(IN)                      :: kd
DOUBLE PRECISION, INTENT(OUT)            :: cv
DOUBLE PRECISION, INTENT(OUT)            :: eg(200)
IMPLICIT DOUBLE PRECISION (a-h,o-z)
DIMENSION b(100),h(100),d(300),e(300),f(300),cv0(100),a(300) ,g(300)

IF (c < 1.0D-10) THEN
  DO  i=1,n
    eg(i)=(i+m)*(i+m-1.0D0)
  END DO
ELSE
  icm=INT((n-m+2)/2)
  nm=10+INT(0.5*(n-m)+c)
  cs=c*c*kd
  DO  l=0,1
    DO  i=1,nm
      IF (l == 0) k=2*(i-1)
      IF (l == 1) k=2*i-1
      dk0=m+k
      dk1=m+k+1
      dk2=2*(m+k)
      d2k=2*m+k
      a(i)=(d2k+2.0)*(d2k+1.0)/((dk2+3.0)*(dk2+5.0))*cs
      d(i)=dk0*dk1+(2.0*dk0*dk1-2.0*m*m-1.0)/((dk2-1.0) *(dk2+3.0))*cs
      g(i)=k*(k-1.0)/((dk2-3.0)*(dk2-1.0))*cs
    END DO
    DO  k=2,nm
      e(k)=DSQRT(a(k-1)*g(k))
      f(k)=e(k)*e(k)
    END DO
    f(1)=0.0D0
    e(1)=0.0D0
    xa=d(nm)+DABS(e(nm))
    xb=d(nm)-DABS(e(nm))
    nm1=nm-1
    DO  i=1,nm1
      t=DABS(e(i))+DABS(e(i+1))
      t1=d(i)+t
      IF (xa < t1) xa=t1
      t1=d(i)-t
      IF (t1 < xb) xb=t1
    END DO
    DO  i=1,icm
      b(i)=xa
      h(i)=xb
    END DO
    DO  k=1,icm
      DO  k1=k,icm
        IF (b(k1) < b(k)) THEN
          b(k)=b(k1)
          EXIT
        END IF
      END DO
      35      IF (k /= 1.AND.h(k) < h(k-1)) h(k)=h(k-1)
      DO WHILE (1==1)
        40       x1=(b(k)+h(k))/2.0D0
        cv0(k)=x1
        IF (DABS((b(k)-h(k))/x1) < 1.0D-14) EXIT
        j=0
        s=1.0D0
        DO  i=1,nm
          IF (s == 0.0D0) s=s+1.0D-30
          t=f(i)/s
          s=d(i)-t-x1
          IF (s < 0.0D0) j=j+1
        END DO
        IF (j < k) THEN
          h(k)=x1
        ELSE
          b(k)=x1
          IF (j >= icm) THEN
            b(icm)=x1
          ELSE
            IF (h(j+1) < x1) h(j+1)=x1
            IF (x1 < b(j)) b(j)=x1
          END IF
        END IF
      END DO
      50      cv0(k)=x1
      IF (l == 0) eg(2*k-1)=cv0(k)
      IF (l == 1) eg(2*k)=cv0(k)
    END DO
  END DO
END IF
70   cv=eg(n-m+1)
RETURN
END SUBROUTINE segv

max_stars_repo_name: vasaantk/bin

max_issues_repo_name: vasaantk/bin

max_forks_repo_name: vasaantk/bin