module handlearray_module

! BAMJ Module containing (inefficient) subroutines to shift from 
! (nx+1,ny+1,...) to ((nx+1)*(ny+1),...) structure ("reduction")
! and vice-versa ("expansion"). They are used different subroutines 
! to pass from rank 4 to 2 or to pass from rank 3 to 2 (and vice-versa).
! TODO crate a single subroutine for both 4->2 and 3->2 (and vice-versa)

interface handlearray
    module procedure handlearray_rdc43
    module procedure handlearray_rdc42
    module procedure handlearray_rdc32
    module procedure handlearray_rdc21
    module procedure handlearray_exp34
    module procedure handlearray_exp24
    module procedure handlearray_exp23
    module procedure handlearray_exp12
end interface handlearray

contains

! Rank 4-->3
subroutine handlearray_rdc43(varin,varout)
    !
    implicit none    
    real*8, dimension(:,:,:,:), intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:,:),   intent(out) :: varout   ! Output variable
    real*8, dimension(:,:,:),   allocatable :: vartemp  ! Auxiliary variable
    integer                                 :: i,j,k1,k2
    integer                                 :: nx,ny,nd,nz,pos
    !
    nx = size(varin,dim=1)
    ny = size(varin,dim=2)
    nd = size(varin,dim=3)
    nz = size(varin,dim=4)
    allocate(vartemp(1:nx*ny,1:nd,1:nz))
    do j=1,ny
       do i=1,nx
          pos = nx*(j-1)+i
          do k1=1,nd
             do k2=1,nz
                vartemp(pos,k1,k2) = varin(i,j,k1,k2)
             end do
          end do
       end do
    end do
    !vartemp(pos+1:m,k+1:m,2:m) = 0 
    varout = vartemp
    deallocate(vartemp)
    !   
end subroutine handlearray_rdc43

! Rank 4-->2
subroutine handlearray_rdc42(varin,varout)
    !
    implicit none    
    real*8, dimension(:,:,:,:), intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:),     intent(out) :: varout   ! Output variable
    real*8, dimension(:,:),     allocatable :: vartemp  ! Auxiliary variable
    integer                                 :: i,j,k
    integer                                 :: nx,ny,nz,pos
    !
    nx = size(varin,dim=1)
    ny = size(varin,dim=2)
    nz = size(varin,dim=4)
    allocate(vartemp(1:nx*ny,1:nz))
    do j=1,ny
       do i=1,nx
          pos = nx*(j-1)+i
          do k=1,nz
             vartemp(pos,k) = varin(i,j,1,k)
          end do
       end do
    end do
    !vartemp(pos+1:m,k+1:m,2:m) = 0 
    varout = vartemp
    deallocate(vartemp)
    !   
end subroutine handlearray_rdc42

! Rank 3-->2
subroutine handlearray_rdc32(varin,varout)
    !
    implicit none
    real*8, dimension(:,:,:), intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:),   intent(out) :: varout   ! Output variable
    real*8, dimension(:,:),   allocatable :: vartemp  ! Auxiliary variable
    integer                               :: i,j,k
    integer                               :: nx,ny,nz,pos
    !
    nx = size(varin,dim=1)
    ny = size(varin,dim=2)
    nz = size(varin,dim=3)
    allocate(vartemp(1:nx*ny,1:nz))
    do j=1,ny
       do i=1,nx
          pos = nx*(j-1)+i
          do k=1,nz
             vartemp(pos,k) = varin(i,j,k)
          end do
       end do
    end do
    !vartemp(pos+1:m,k+1:m,2:m) = 0 
    varout = vartemp
    deallocate(vartemp)
    !   
end subroutine handlearray_rdc32

! Rank 2-->1
subroutine handlearray_rdc21(varin,varout)
    !
    implicit none
    real*8, dimension(:,:), intent(in)  :: varin    ! Input variable
    real*8, dimension(:),   intent(out) :: varout   ! Output variable
    real*8, dimension(:),   allocatable :: vartemp  ! Auxiliary variable
    integer                             :: i,j
    integer                             :: nx,ny,pos
    !
    nx = size(varin,dim=1)
    ny = size(varin,dim=2)
    allocate(vartemp(1:nx*ny))
    do j=1,ny
       do i=1,nx
          pos = nx*(j-1)+i
          vartemp(pos) = varin(i,j)
       end do
    end do
    varout = vartemp
    deallocate(vartemp)
    !   
end subroutine handlearray_rdc21  
    
! Rank 3-->4      
subroutine handlearray_exp34(varin,varout)
    !
    implicit none
    real*8, dimension(:,:,:),   intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:,:,:), intent(out) :: varout   ! Output variable
    real*8, dimension(:,:,:,:), allocatable :: vartemp  ! Auxiliary variable
    integer                                 :: i,j,k1,k2
    integer                                 :: nx,ny,nd,nz
    real*8                                  :: pos
    real*8                                  :: pos1
    !
    nx = size(varout,dim=1)
    ny = size(varout,dim=2)
    nd = size(varout,dim=3)
    nz = size(varout,dim=4)
    allocate(vartemp(1:nx,1:ny,1:nd,1:nz))
    do pos=1,ny*nx
          pos1 = pos/nx
          j = ceiling(pos1)
          i = pos-(nx)*(j-1)
       do k1=1,nd
          do k2=1,nz
             vartemp(i,j,k1,k2) = varin(pos,k1,k2)
          end do
       end do
    end do
    varout = vartemp
    deallocate(vartemp)
    !
end subroutine handlearray_exp34
       
! Rank 2-->4      
subroutine handlearray_exp24(varin,varout)
    !
    implicit none
    real*8, dimension(:,:),     intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:,:,:), intent(out) :: varout   ! Output variable
    real*8, dimension(:,:,:,:), allocatable :: vartemp  ! Auxiliary variable
    integer                                 :: i,j,k
    integer                                 :: nx,ny,nz
    real*8                                  :: pos
    real*8                                  :: pos1
    !
    nx = size(varout,dim=1)
    ny = size(varout,dim=2)
    nz = size(varout,dim=4)
    allocate(vartemp(1:nx,1:ny,1,1:nz))
    do pos=1,ny*nx
          pos1 = pos/nx
          j = ceiling(pos1)
          i = pos-(nx)*(j-1)
       do k=1,nz
          vartemp(i,j,1,k) = varin(pos,k)
       end do
    end do
    varout = vartemp
    deallocate(vartemp)
    !
end subroutine handlearray_exp24

! Rank 2-->3
subroutine handlearray_exp23(varin,varout)
    !
    implicit none
    real*8, dimension(:,:),   intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:,:), intent(out) :: varout   ! Output variable
    real*8, dimension(:,:,:), allocatable :: vartemp  ! Auxiliary variable
    integer                               :: i,j,k
    integer                               :: nx,ny,nz
    real*8                                :: pos
    real*8                                :: pos1
    !
    nx = size(varout,dim=1)
    ny = size(varout,dim=2)
    nz = size(varout,dim=3)
    allocate(vartemp(1:nx,1:ny,1:nz))
    do pos=1,ny*nx
          pos1 = pos/nx
          j = ceiling(pos1)
          i = pos-(nx)*(j-1)
       do k=1,nz
          vartemp(i,j,k) = varin(pos,k)
       end do
    end do
    varout = vartemp
    deallocate(vartemp)
    !
end subroutine handlearray_exp23

! Rank 1-->2
subroutine handlearray_exp12(varin,varout)
    !
    implicit none
    real*8, dimension(:),   intent(in)  :: varin    ! Input variable
    real*8, dimension(:,:), intent(out) :: varout   ! Output variable
    real*8, dimension(:,:), allocatable :: vartemp  ! Auxiliary variable
    integer                             :: i,j
    integer                             :: nx,ny
    real*8                              :: pos
    real*8                              :: pos1
    !
    nx = size(varout,dim=1)
    ny = size(varout,dim=2)
    allocate(vartemp(1:nx,1:ny))
    do pos=1,ny*nx
          pos1 = pos/nx
          j = ceiling(pos1)
          i = pos-(nx)*(j-1)
          vartemp(i,j) = varin(pos)
    end do
    varout = vartemp
    deallocate(vartemp)
    !
end subroutine handlearray_exp12

end module handlearray_module