!Code for generating initial field for homogeneous isotropic turbulence
!Adapted from M. Shoeybi's code

!compilation
!export PATH=/usr/bin:$PATH
!gfortran -I/usr/local/Cellar/fftw/3.3.7/include/ -o generate_init_hit generate_init_hit.f90 -L/usr/local/Cellar/fftw/3.3.7/lib -lfftw3

! gfortran -I/usr/include/ -o generate_init_hit generate_init_hit.f90 -lfftw3

!ifort -I$WORKDIR/SRC/fftw/include -o generate_init_hit generate_init_hit.f90 -L$WORKDIR/SRC/fftw/lib -lfftw3

module rand_tools_m

contains

  subroutine randset(idum)

    implicit none

    integer :: idum

    real(8) :: r(97)
    integer :: ix1,ix2,ix3
    common /random_common/ r,ix1,ix2,ix3

    real(8) :: rm1,rm2
    integer :: m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3
    parameter(m1=259200,ia1=7141,ic1=54773,rm1=1.d0/m1)
    parameter(m2=134456,ia2=8121,ic2=28441,rm2=1.d0/m2)
    parameter(m3=243000,ia3=4561,ic3=51349)

    integer :: j

    ix1 = mod(ic1 - idum , m1)
    ix1 = mod(ia1 * ix1 + ic1 , m1)
    ix2 = mod(ix1 , m2)
    ix1 = mod(ia1 * ix1 + ic1 , m1)
    ix3 = mod(ix1 , m3)

    do j = 1, 97
       ix1 = mod(ia1 * ix1 + ic1 , m1)
       ix2 = mod(ia2 * ix2 + ic2 , m2)
       r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
       if (r(j).gt.1.d0) r(j)=1.d0
       if (r(j).lt.0.d0) r(j)=0.d0
    enddo

  end subroutine randset

  subroutine randvec(ran,N)

    implicit none

    integer :: N
    real(8) :: ran(N)
    real(8) :: r(97)
    integer :: ix1,ix2,ix3
    common /random_common/ r,ix1,ix2,ix3

    real(8) :: rm1,rm2
    integer :: m1,m2,m3,ia1,ia2,ia3,ic1,ic2,ic3
    parameter(m1=259200,ia1=7141,ic1=54773,rm1=1.d0/m1)
    parameter(m2=134456,ia2=8121,ic2=28441,rm2=1.d0/m2)
    parameter(m3=243000,ia3=4561,ic3=51349)

    integer :: nombre,j

    do nombre=1,N
       ix1 = mod(ia1 * ix1 + ic1 , m1)
       ix2 = mod(ia2 * ix2 + ic2 , m2)
       ix3 = mod(ia3 * ix3 + ic3 , m3)

       j = 1 + mod(((97 * ix3) / m3),97)
       ran(nombre) = r(j)
       r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
       if (r(j).gt.1.d0) r(j)=1.d0
       if (r(j).lt.0.d0) r(j)=0.d0
    enddo

  end subroutine randvec

end module rand_tools_m

program generate_init_hit

  use rand_tools_m

  implicit none

  !include fftw header file
  include "fftw3.f"

  integer, parameter :: N_box = 129

  real(8), parameter :: k0 = 6.0_8

  integer, parameter :: divergence_flag = 0
  !            0: use the noraml construction which gives
  !               the divergecne free for spectral codes
  !            1: use modified wave numbers

  !constants
  real(8) :: pi, twopi

  ! grid information
  real(8) :: my_node_cc_min(3), my_node_cc_max(3), node_cc_min(3), node_cc_max(3)
  real(8) :: my_box_volume, box_volume
  real(8) :: length(3), dx(3)
  integer  :: N_total, num_node

  ! initial condition in Fourier space
  real(8) :: u0
  real(8) :: kx, ky, kz, kxyz, kmax, tolerance, kmin, kxy, kmax_m
  integer :: kx_max_i, ky_max_i, kz_max_i

  ! energy parameters
  real(8) :: Eq_n
  complex(8) :: a_ran, b_ran
  complex, parameter :: ima = (0.0_8,1.0_8)

  ! random number
  real(8), allocatable :: phi1(:,:,:), phi2(:,:,:), phi3(:,:,:)
  real(8) :: max_phi12, max_phi13, max_phi23, min_dphi

  ! Fourier ans velocity vectors
  real(8), allocatable :: u(:,:,:), u3(:,:,:,:), u3_ext(:,:,:,:)
  complex(8), allocatable :: uk(:,:,:)!,uk_div(:,:,:,:)
  integer*8 :: fftw_plan

  ! debug stuff
  real(8) :: my_max_u, my_min_u, my_mean_u(3)
  real(8) :: max_u, min_u, mean_u(3), max_divergence
  real(8) :: u_left(1:3), u_right(1:3)
  real(8), allocatable :: div(:,:,:)

  ! loop parameters
  integer :: i, j, k, i_u, ino, ierr, n_phi

  ! dump out file
  character(len=64) :: filename
  !-------------------------------------

  write(*,*) '     Generate initial condition for HIT'

  pi = 4.0_8 * atan(1.0_8)
  twopi = 2.0_8 * pi

  if ( divergence_flag.eq.0 ) then
     write(*,*) '       using spectral divergence free field'
  elseif ( divergence_flag.eq.1 ) then
     write(*,*) '       using modified wave number to get finite differecne divergence free field'
  else
     write(*,*) 'Error in in calc_init_for_hom_iso_tur: undefined divergence_flag'
     stop
  endif

  ! total number of grid points
  N_total = N_box * N_box * N_box

  ! check even number of points has been passed
  num_node = N_box - 1
  if ( mod(num_node,2).ne.0 ) then
     write(*,*) 'Error in calc_init_for_hom_iso_tur, pass odd number of points '
     stop
  endif

  ! get the grid information
  ! box length
  length(1:3) = twopi

  ! grid spacing
  dx(1:3) = length(1:3) / dble(N_box-1)

  ! write some data on the screen
  write(*,*) '*****************************************************************************'
  write(*,*) '  grid information'
  write(*,*) '     total number of grid points(cart)               =  ', N_total
  write(*,*) '     number of grid points in each direction(cart)   =  ', N_box
  write(*,*) '     length in x, y and z direction (should be 2*pi) =  ', length(1:3)
  write(*,*) '     grid spacing in x, y and z(cart)                =  ', dx(1:3)
  write(*,*) '*****************************************************************************'

  ! spectrum paprameters
  tolerance = 0.0001_8
  kmax = dble(N_box-1) / 2.0_8 - tolerance
  kmin = tolerance
  kx_max_i = idint( dble(N_box-1) / 2.0_8 + 0.1_8 ) + 1
  ky_max_i = N_box-1
  kz_max_i = N_box-1

  ! spectrum information
  ! write on the screen
  write(*,*) '*****************************************************************************'
  write(*,*) 'spectrum properties'
  write(*,*) '   wave number corresponding to most energetic wave =  ', k0
  write(*,*) '*****************************************************************************'
  write(*,*) 'wave number information'
  write(*,*) '   maximun wave number(for homogeniety) =  ', kmax
  write(*,*) '   min wave number                      =  ', kmin
  write(*,*) '   max wave number in x direction       = ',kx_max_i
  write(*,*) '   max wave number in y direction       = ',ky_max_i
  write(*,*) '   max wave number in z direction       = ',kz_max_i
  write(*,*) '*****************************************************************************'

  ! random numbers
  ! allocate memory
  allocate( phi1(kx_max_i,ky_max_i,kz_max_i) )
  allocate( phi2(kx_max_i,ky_max_i,kz_max_i) )
  allocate( phi3(kx_max_i,ky_max_i,kz_max_i) )

  ! call random number generator
  call randset(1)
  call randvec(phi1,kx_max_i*ky_max_i*kz_max_i)
  call randvec(phi2,kx_max_i*ky_max_i*kz_max_i)
  call randvec(phi3,kx_max_i*ky_max_i*kz_max_i)

  ! check the random numbers
  max_phi12 = maxval( abs( phi1(1:kx_max_i,1:ky_max_i,1:kz_max_i) -&
       phi2(1:kx_max_i,1:ky_max_i,1:kz_max_i) ) )
  max_phi13 = maxval( abs( phi1(1:kx_max_i,1:ky_max_i,1:kz_max_i) -&
       phi3(1:kx_max_i,1:ky_max_i,1:kz_max_i) ) )
  max_phi23 = maxval( abs( phi3(1:kx_max_i,1:ky_max_i,1:kz_max_i) -&
       phi2(1:kx_max_i,1:ky_max_i,1:kz_max_i) ) )

  min_dphi = min(max_phi12,max_phi13)
  min_dphi = min(min_dphi,max_phi23)

  if ( min_dphi.le.tolerance ) then
     write(*,*) 'Error in calc_init_cond_for_hom_iso_tu: random numbers are the same'
     stop
  endif

  ! allocate memory
  allocate( uk(kx_max_i,ky_max_i,kz_max_i) )
  allocate( u(ky_max_i,ky_max_i,kz_max_i) )
  allocate( u3(3,ky_max_i,ky_max_i,kz_max_i) )
  ! allocate( uk_div(3,kx_max_i,ky_max_i,kz_max_i) )

  ! modify the random numbers to get them in the
  ! interval(-pi,pi)
  do k = 1,kz_max_i
     do j = 1,ky_max_i
        do i = 1,kx_max_i
           phi1(i,j,k) = twopi * (phi1(i,j,k)-0.5_8)
           phi2(i,j,k) = twopi * (phi2(i,j,k)-0.5_8)
           phi3(i,j,k) = twopi * (phi3(i,j,k)-0.5_8)
        enddo
     enddo
  enddo

  ! main loop for the fourier transforms
  do i_u = 1,3
     uk(1:kx_max_i,1:ky_max_i,1:kz_max_i) = ( 0.0_8, 0.0_8)
     u(1:ky_max_i,1:ky_max_i,1:kz_max_i) = 0.0_8
     do k = 1,kz_max_i
        do j = 1,ky_max_i
           do i = 1,kx_max_i

              ! modify the waves to get them symmetric
              kx = dble(i-1)
              ky = dble(j-1)
              if ( j.gt.kx_max_i ) ky = -dble(ky_max_i-j+1)
              kz = dble(k-1)
              if ( k.gt.kx_max_i ) kz = -dble(kz_max_i-k+1)

              kmax_m = dsqrt( kx**2 + ky**2 + kz**2 )

              kxyz = dsqrt( kx**2 + ky**2 + kz**2 )
              kxy = dsqrt( kx**2 + ky**2 )

              ! enegry spectrum
              Eq_n = (kxyz/k0)**4 * dexp( -2.0_8*(kxyz/k0)**2 )

              a_ran = (0.0_8,0.0_8)
              b_ran = (0.0_8,0.0_8)

              if ( kxyz.gt.kmin .and. kxyz.le.kmax_m ) then
                 a_ran = dsqrt( Eq_n/(twopi* kxyz**2) ) *&
                      cdexp(ima*phi1(i,j,k)) * dcos(phi3(i,j,k))
                 b_ran = dsqrt( Eq_n/(twopi* kxyz**2) ) *&
                      cdexp(ima*phi2(i,j,k)) * dsin(phi3(i,j,k))
              endif

              if ( divergence_flag.eq.1 ) then
                 ! replace with modefied wave number
                 kx = sin( kx*twopi/dble(ky_max_i) ) / (twopi/dble(ky_max_i))
                 ky = sin( ky*twopi/dble(ky_max_i) ) / (twopi/dble(ky_max_i))
                 kz = sin( kz*twopi/dble(ky_max_i) ) / (twopi/dble(ky_max_i))
              endif

              kxyz = dsqrt( kx**2 + ky**2 + kz**2 )
              kxy = dsqrt( kx**2 + ky**2 )

              ! fourier coefficients
              uk(i,j,k) = (0.0_8,0.0_8)
              if ( kxyz.gt.kmin ) then

                 if ( i_u.eq.1 ) then

                    if ( kxy.gt.kmin ) then
                       uk(i,j,k) = (a_ran*kxyz*ky+b_ran*kx*kz)/(kxyz*kxy)
                    else
                       !uk(i,j,k) = a_ran
                       uk(i,j,k) = ( a_ran + b_ran ) / dsqrt(2.0_8)
                    endif

                 elseif ( i_u.eq.2 ) then

                    if ( kxy.gt.kmin ) then
                       uk(i,j,k) = (b_ran*ky*kz-a_ran*kxyz*kx)/(kxyz*kxy)
                    else
                       !uk(i,j,k) = b_ran
                       uk(i,j,k) = ( b_ran - a_ran ) / dsqrt(2.0_8)
                    endif

                 elseif ( i_u.eq.3 ) then

                    !if ( kxy.gt.kmin ) then
                    uk(i,j,k) = -b_ran * kxy / kxyz
                    !else
                    !   uk(i,j,k) = (0.0_8,0.0_8)
                    !endif

                 endif

              endif

           enddo
        enddo
     enddo

     ! at plane kx=0, the waves need to be complex conjugate
     do k = 2,kz_max_i
        do j = kx_max_i+1,ky_max_i
           uk(1,j,k) = dconjg(uk(1,ky_max_i-j+2,kz_max_i-k+2))
        enddo
     enddo
     do k = kx_max_i+1,kz_max_i
        uk(1,1,k) = dconjg(uk(1,1,kz_max_i-k+2))
     enddo
     do j = kx_max_i+1,ky_max_i
        uk(1,j,1) = dconjg(uk(1,ky_max_i-j+2,1))
     enddo

     ! remove odd-ball wave numbers
     uk(kx_max_i,1:ky_max_i,1:kz_max_i) = (0.0_8,0.0_8)
     uk(1:kx_max_i,kx_max_i,1:kz_max_i) = (0.0_8,0.0_8)
     uk(1:kx_max_i,1:ky_max_i,kx_max_i) = (0.0_8,0.0_8)

     ! remove mean
     uk(1,1,1) = (0.0_8,0.0_8)

     !take the fourier transform
     call dfftw_plan_dft_c2r_3d(fftw_plan, N_box-1, N_box-1, N_box-1, uk, u, 64)
     call dfftw_execute(fftw_plan)
     call dfftw_destroy_plan(fftw_plan)

     u3(i_u,1:N_box-1,1:N_box-1,1:N_box-1) = u(1:N_box-1,1:N_box-1,1:N_box-1)

  enddo

  ! deallocate memory
  deallocate(phi1); deallocate(phi2); deallocate(phi3)
  deallocate(uk); deallocate(u);

  !add the periodic points to the grid

  allocate( u3_ext(3,ky_max_i+1,ky_max_i+1,kz_max_i+1) )
  u3_ext(1:3,1:ky_max_i,1:ky_max_i,1:kz_max_i) = u3(1:3,1:ky_max_i,1:ky_max_i,1:kz_max_i)
  !periodic boundaries
  u3_ext(1:3,ky_max_i+1,1:ky_max_i,1:kz_max_i) = u3(1:3,1,1:ky_max_i,1:kz_max_i)
  u3_ext(1:3,1:ky_max_i,ky_max_i+1,1:kz_max_i) = u3(1:3,1:ky_max_i,1,1:kz_max_i)
  u3_ext(1:3,1:ky_max_i,1:ky_max_i,kz_max_i+1) = u3(1:3,1:ky_max_i,1:ky_max_i,1)

  u3_ext(1:3,ky_max_i+1,ky_max_i+1,1:kz_max_i) = u3(1:3,1,1,1:kz_max_i)
  u3_ext(1:3,ky_max_i+1,1:ky_max_i,kz_max_i+1) = u3(1:3,1,1:ky_max_i,1)
  u3_ext(1:3,1:ky_max_i,ky_max_i+1,kz_max_i+1) = u3(1:3,1:ky_max_i,1,1)
  u3_ext(1:3,ky_max_i+1,ky_max_i+1,kz_max_i+1) = u3(1:3,1,1,1)

  deallocate(u3)

  ! calculate divergence
  allocate(div(ky_max_i,ky_max_i,kz_max_i))

  do k = 1, kz_max_i
     do j = 1, ky_max_i
        do i = 1, kx_max_i

           if (i==1) then
              u_left(1) = u3_ext(1,ky_max_i,j,k)
           else
              u_left(1) = u3_ext(1,i-1,j,k)
           end if
           u_right(1) = u3_ext(1,i+1,j,k)

           if (j==1) then
              u_left(2) = u3_ext(2,i,ky_max_i,k)
           else
              u_left(2) = u3_ext(2,i,j-1,k)
           end if
           u_right(2) = u3_ext(2,i,j+1,k)

           if (k==1) then
              u_left(3) = u3_ext(3,i,j,kz_max_i)
           else
              u_left(3) = u3_ext(3,i,j,k-1)
           end if
           u_right(3) = u3_ext(3,i,j,k+1)

           div(i,j,k) = (u_right(1)-u_left(1))/2.0d0/dx(1) + &
                        (u_right(2)-u_left(2))/2.0d0/dx(2) + &
                        (u_right(3)-u_left(3))/2.0d0/dx(3)

        end do
     end do
  end do

  max_divergence = maxval(div)

  !dump the statistics on the screen
  write(*,*) '*****************************************************************************'
  write(*,*) 'initial condition'
  write(*,*) '   maximum mach number             =  ', max_u
  write(*,*) '   minimum mach number             =  ', min_u
  write(*,*) '   mean velocity                   =  ', mean_u(1:3)/box_volume
  write(*,*) '   mean velocity2                  =  ', mean_u(1:3)
  write(*,*) '   maximum divergence of velocity  =  ', max_divergence
  write(*,*) '*****************************************************************************'

  !write the velocity field
  write(filename,'(a)') './init_field_hit.dat'
  open(10, file=filename, form='formatted', action='write')
  write(10,'(a)') 'VARIABLES = "X" "Y" "Z" "U" "V" "W"'
  write(10,'(a,i4.4,a,i4.4,a,i4.4)') "ZONE I=", ky_max_i+1, ", J=", ky_max_i+1, " K=", kz_max_i+1
  do k = 1, kz_max_i+1
     do j = 1, ky_max_i+1
        do i = 1, ky_max_i+1
           write(10,'(6e24.16)') dble(i-1)*dx(1), dble(j-1)*dx(2), dble(k-1)*dx(3), u3_ext(1:3,i,j,k)
        end do
     end do
  end do
  close(10)

  !deallocate memory
  deallocate(div)
  deallocate(u3_ext)

end program generate_init_hit


! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.

!>
!! Provide wrapped allocate and deallocate routines with basic error handling
!! and memory usage output.
!!
module ELSI_MALLOC

   use ELSI_DATATYPE, only: elsi_basic_t
   use ELSI_MPI
   use ELSI_OUTPUT, only: elsi_say
   use ELSI_PRECISION, only: i4,i8,r4,r8

   implicit none

   private

   public :: elsi_allocate
   public :: elsi_deallocate

   interface elsi_allocate
      module procedure elsi_allocate_integer4_1d
      module procedure elsi_allocate_integer4_2d
      module procedure elsi_allocate_integer4_3d
      module procedure elsi_allocate_integer8_1d
      module procedure elsi_allocate_real8_1d
      module procedure elsi_allocate_real8_2d
      module procedure elsi_allocate_real8_3d
      module procedure elsi_allocate_real4_1d
      module procedure elsi_allocate_real4_2d
      module procedure elsi_allocate_complex16_1d
      module procedure elsi_allocate_complex16_2d
      module procedure elsi_allocate_complex16_3d
      module procedure elsi_allocate_complex8_2d
   end interface

   interface elsi_deallocate
      module procedure elsi_deallocate_integer4_1d
      module procedure elsi_deallocate_integer4_2d
      module procedure elsi_deallocate_integer4_3d
      module procedure elsi_deallocate_integer8_1d
      module procedure elsi_deallocate_real8_1d
      module procedure elsi_deallocate_real8_2d
      module procedure elsi_deallocate_real8_3d
      module procedure elsi_deallocate_real4_1d
      module procedure elsi_deallocate_real4_2d
      module procedure elsi_deallocate_complex16_1d
      module procedure elsi_deallocate_complex16_2d
      module procedure elsi_deallocate_complex16_3d
      module procedure elsi_deallocate_complex8_2d
   end interface

contains

!>
!! Allocate a 1D array with real(kind=r4).
!!
subroutine elsi_allocate_real4_1d(bh,array,dim_1,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r4), intent(inout), allocatable :: array(:)
   integer(kind=i4), intent(in) :: dim_1
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*4

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*4
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:) = 0.0_r4

end subroutine

!>
!! Allocate a 1D array with real(kind=r8).
!!
subroutine elsi_allocate_real8_1d(bh,array,dim_1,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:)
   integer(kind=i4), intent(in) :: dim_1
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*8

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*8
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:) = 0.0_r8

end subroutine

!>
!! Allocate a 1D array with integer(kind=i4).
!!
subroutine elsi_allocate_integer4_1d(bh,array,dim_1,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:)
   integer(kind=i4), intent(in) :: dim_1
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*4

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*4
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:) = 0

end subroutine

!>
!! Allocate a 1D array with integer(kind=i8).
!!
subroutine elsi_allocate_integer8_1d(bh,array,dim_1,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i8), intent(inout), allocatable :: array(:)
   integer(kind=i4), intent(in) :: dim_1
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*8

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*8
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:) = 0

end subroutine

!>
!! Allocate a 1D array with complex(kind=r8).
!!
subroutine elsi_allocate_complex16_1d(bh,array,dim_1,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:)
   integer(kind=i4), intent(in) :: dim_1
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*16

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*16
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:) = (0.0_r8,0.0_r8)

end subroutine

!>
!! Allocate a 2D array of real(kind=r4).
!!
subroutine elsi_allocate_real4_2d(bh,array,dim_1,dim_2,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r4), intent(inout), allocatable :: array(:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*4

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*4
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:) = 0.0_r4

end subroutine

!>
!! Allocate a 2D array of real(kind=r8).
!!
subroutine elsi_allocate_real8_2d(bh,array,dim_1,dim_2,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*8

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*8
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:) = 0.0_r8

end subroutine

!>
!! Allocate a 2D array of integer(kind=i4).
!!
subroutine elsi_allocate_integer4_2d(bh,array,dim_1,dim_2,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*4

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*4
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:) = 0

end subroutine

!>
!! Allocate a 2D array of complex(kind=r4).
!!
subroutine elsi_allocate_complex8_2d(bh,array,dim_1,dim_2,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r4), intent(inout), allocatable :: array(:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*8

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*8
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:) = (0.0_r4,0.0_r4)

end subroutine

!>
!! Allocate a 2D array of complex(kind=r8).
!!
subroutine elsi_allocate_complex16_2d(bh,array,dim_1,dim_2,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*16

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*16
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:) = (0.0_r8,0.0_r8)

end subroutine

!>
!! Allocate a 3D array of real(kind=r8).
!!
subroutine elsi_allocate_real8_3d(bh,array,dim_1,dim_2,dim_3,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:,:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   integer(kind=i4), intent(in) :: dim_3
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*8

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2,dim_3),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*8
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:,:) = 0.0_r8

end subroutine

!>
!! Allocate a 3D array of integer(kind=i4).
!!
subroutine elsi_allocate_integer4_3d(bh,array,dim_1,dim_2,dim_3,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:,:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   integer(kind=i4), intent(in) :: dim_3
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*4

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2,dim_3),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*4
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:,:) = 0

end subroutine

!>
!! Allocate a 3D array of complex(kind=r8).
!!
subroutine elsi_allocate_complex16_3d(bh,array,dim_1,dim_2,dim_3,label,caller)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:,:,:)
   integer(kind=i4), intent(in) :: dim_1
   integer(kind=i4), intent(in) :: dim_2
   integer(kind=i4), intent(in) :: dim_3
   character(len=*), intent(in) :: label
   character(len=*), intent(in) :: caller

   real(kind=r8) :: mem
   integer(kind=i4) :: ierr
   character(len=200) :: msg

   if(bh%print_info > 2) then
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*16

      write(msg,"(4X,A,F10.3,2A)") "Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
   end if

   allocate(array(dim_1,dim_2,dim_3),stat=ierr)

   if(ierr > 0) then
      write(msg,"(2A)") "Error in allocating ",trim(label)
      call elsi_say(bh,msg)
      mem = 1.0e-6_r8*dim_1*dim_2*dim_3*16
      write(msg,"(4X,A,F10.3,2A)") "Trying to Allocating ",mem," MB for ",trim(label)
      call elsi_say(bh,msg)
      write(msg,*) "Error number ", ierr
      call elsi_say(bh,msg)
      write(msg,"(A)") "This probably means not enough memory"
      call elsi_stop(bh,msg,caller)
   end if

   array(:,:,:) = (0.0_r8,0.0_r8)

end subroutine

!>
!! Deallocate a 1D array with real(kind=r4).
!!
subroutine elsi_deallocate_real4_1d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r4), intent(inout), allocatable :: array(:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 1D array with real(kind=r8).
!!
subroutine elsi_deallocate_real8_1d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 1D array with integer(kind=i4).
!!
subroutine elsi_deallocate_integer4_1d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 1D array with integer(kind=i8).
!!
subroutine elsi_deallocate_integer8_1d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i8), intent(inout), allocatable :: array(:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 1D array with complex(kind=r8).
!!
subroutine elsi_deallocate_complex16_1d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 2D array with real(kind=r4).
!!
subroutine elsi_deallocate_real4_2d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r4), intent(inout), allocatable :: array(:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 2D array with real(kind=r8).
!!
subroutine elsi_deallocate_real8_2d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 2D array with integer(kind=i4).
!!
subroutine elsi_deallocate_integer4_2d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 2D array with complex(kind=r4).
!!
subroutine elsi_deallocate_complex8_2d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r4), intent(inout), allocatable :: array(:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!>
!! Deallocate a 2D array with complex(kind=r8).
!!
subroutine elsi_deallocate_complex16_2d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 3D array with real(kind=r8).
!!
subroutine elsi_deallocate_real8_3d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   real(kind=r8), intent(inout), allocatable :: array(:,:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 3D array with integer(kind=i4).
!!
subroutine elsi_deallocate_integer4_3d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   integer(kind=i4), intent(inout), allocatable :: array(:,:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

!>
!! Deallocate a 3D array with complex(kind=r8).
!!
subroutine elsi_deallocate_complex16_3d(bh,array,label)

   implicit none

   type(elsi_basic_t), intent(in) :: bh
   complex(kind=r8), intent(inout), allocatable :: array(:,:,:)
   character(len=*), intent(in) :: label

   character(len=200) :: msg

   if(bh%print_info > 2) then
      write(msg,"(4X,2A)") "Deallocating ",trim(label)
      call elsi_say(bh,msg)
   end if

   deallocate(array)

end subroutine

end module ELSI_MALLOC


SUBROUTINE calc_grav_brute
  ! Subroutine calculates gravitational forces and potentials for all particles
  ! It uses brute force O(N^2) to calculate forces and potentials
  ! for particles outside the smoothing volume

  use sphdata
  use sphneighbourdata
  use sphkerneldata
  use tachedata, only:nelement

  implicit none

  real,dimension(3) :: gravi

  integer :: ielement, jpart,k,ix
  real :: meanpot,sdpot, percent,increment

  allocate(gravxyz(3,nelement))
  allocate(poten(nelement))

  gravxyz(:,:) = 0.0
  poten(:) = 0.0

  isoft = 0

  percent = 0.0
  increment = 1.0

  do ielement = 1,nelement

     call element_percent_complete(ielement,nelement,percent,increment)

     gravi(:) = 0.0
     poteni = 0.0

     ! Loop over all elements

     !$OMP PARALLEL &
     !$OMP shared(ielement, nelement,poten,gravxyz) &
     !$OMP private(jpart,poteni,gravi,ix)
     !$OMP DO SCHEDULE(runtime)
     do jpart = 1,nelement

        if(ielement==jpart) cycle

        ! Compute the forces between ielement and jpart
        call particle_forces(ielement,jpart,poteni,gravi)

        poten(ielement) = poten(ielement) + poteni

        DO ix=1,3
           gravxyz(ix,ielement) = gravxyz(ix,ielement) + gravi(ix)
        ENDDO

     enddo
     !$OMP END DO
     !$OMP END PARALLEL
     ! End of loop over jpart

     !  Add contribution to the potential from pointmasses

     do k=1,nptmass

        jpart = listpm(k)

        if(jpart==ielement) cycle

        sep = (xyzmh(1,ielement) - xyzmh(1,jpart))**2 + &
             (xyzmh(2,ielement) - xyzmh(2,jpart))**2 +&
             (xyzmh(3,ielement) - xyzmh(3,jpart))**2

        sep =sqrt(sep)

        poten(ielement) = poten(ielement) + xyzmh(4,jpart)/sep

     enddo
     !print*, 'Particle ', ielement, ': Potential - ', poten(ielement), 'Force: ',gravi(:)
  ENDDO

  ! End loop over all particles

  ! Calculate mean and standard deviation of potential

  meanpot = sum(poten)/REAL(nelement)
  sdpot = 0.0

  !$OMP PARALLEL &
  !$OMP shared(nneigh,meanneigh,nelement)&
  !$OMP private(ielement) &
  !$OMP reduction(+:sdneigh)
  !$OMP DO SCHEDULE(runtime)
  do ielement=1,nelement
     sdpot = sdpot+(poten(ielement)-meanpot)**2
  enddo
  !$OMP END DO
  !$OMP END PARALLEL

  sdpot = sqrt(sdpot/REAL(nelement))

  print*, 'Mean potential is ', meanpot
  print*, 'Standard Deviation: ', sdpot

  return
END SUBROUTINE calc_grav_brute


text:
!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2014  DNV KEMA Renewables, Inc.
!
!    This file is part of the IceFloe suite of subroutines
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!************************************************************************

!****************************************************************
! Module to handle parameter input from user's file

module iceInputParams

   use , INTRINSIC :: iso_fortran_env , only : iostat_end  ! for end of file detection
   use precision
   use NWTC_IO, only : GetNewUnit, OpenFInpFile
   use iceLog

   implicit none

!  This is used for initial read in of variables by name
!  Extraction of values into specific named variables done by getIceInput functions
   type paramType
      character(132) :: name
      real(ReKi)     :: value
   end type paramType

!  Hold some meta data for the input parameters
   type iceInputType
      integer(IntKi) :: unitNum
      integer(IntKi) :: count
      type(paramType), allocatable :: params(:)
   end type iceInputType

!  generic functions for retrieving input parameters
   interface getIceInput
      module procedure getRealInput
      module procedure getIntInput
      module procedure getLogicalInput
   end interface getIceInput

!  data structures for tower - much of this could come from FAST/ElastoDyn/SubDyn
!  but need for other codes
   type twrLegs
      real(ReKi)     :: X
      real(ReKi)     :: Y
      real(ReKi)     :: phase
      real(ReKi)     :: ks
   end type twrLegs
   type towerData
      real(ReKi)     :: diam           ! at water line - this could come from FAST but not other codes
      type(twrLegs), allocatable :: leg(:)
      real(ReKi)     :: freq           ! leg or tower fundamental frequency, Hz
      real(ReKi)     :: coneAngle      ! for flexural failure
      real(ReKi)     :: coneTopDiam    ! for flexural failure
   end type towerData

!  This is the list of input parameters
   type inputParams
      integer(IntKi) :: iceType
      real(ReKi)     :: iceThickness
      real(ReKi)     :: iceVelocity
      real(ReKi)     :: iceDirection
      real(ReKi)     :: timeStep
      real(ReKi)     :: rampTime
      integer(IntKi) :: randomSeed

      type(towerData):: twr
      real(ReKi)     :: multiLegFactor_kn ! Applied to peak load lock-in for all legs when numLegs>1
      logical        :: singleLoad        ! true if multiple leg loads will be combined into an equivalent load
      logical        :: legAutoFactor     ! code decides on leg load shelter factors if true

      real(ReKi)     :: refIceThick
      real(ReKi)     :: refIceStrength
      real(ReKi)     :: staticExponent

      real(ReKi)     :: shapeFactor_k1
      real(ReKi)     :: contactFactor_k2

      real(ReKi)     :: coeffPSD_b
      real(ReKi)     :: coeffPSD_ks
      real(ReKi)     :: crushLoadCOV
      real(ReKi)     :: stdLoadMult
      real(ReKi)     :: freqStep

      real(ReKi)     :: interPeriod
      real(ReKi)     :: riseTime
      real(ReKi)     :: fallTime

      real(ReKi)     :: minLoadFraction

      real(ReKi)     :: flexStrength
      real(ReKi)     :: ice2twrFriction
      real(ReKi)     :: iceDensity
      logical        :: includeHb
      logical        :: includeHr

      real(ReKi)     :: peakLoadCOV
      real(ReKi)     :: coeffLoadPeaks
      real(ReKi)     :: coeffLoadMin
      real(ReKi)     :: periodCOV
      real(ReKi)     :: tauMin
      real(ReKi)     :: tauMax
      real(ReKi)     :: coeffBreakLength
      real(ReKi)     :: rubbleHeight
      real(ReKi)     :: rubblePorosity
      real(ReKi)     :: rubbleCohesion
      real(ReKi)     :: rubbleAngle
      real(ReKi)     :: frictionAngle
      real(ReKi)     :: ice2iceFriction
      real(ReKi)     :: waterDensity
      real(ReKi)     :: iceModulus
      real(ReKi)     :: poissonRatio
      logical        :: includeHp
      logical        :: includeHl
      logical        :: includeHt
      logical        :: includeLc

      real(ReKi)     :: rideUpThickness
      real(ReKi)     :: freqParamK

      real(ReKi)     :: minStrength
      real(ReKi)     :: minStrengthNegVel
    end type inputParams

contains

!========================================================================
   subroutine countIceInputs (fname, iceLog, input)
      character(*), intent(in)                  :: fname    ! name of file with input parameters
      type(iceFloe_LoggingType), intent(inout)   :: iceLog   ! structure with message and error logging variables
      type(iceInputType), intent(out)              :: input

      character(1024)   :: lineIn
      INTEGER(IntKi)    :: ErrStat     ! Error status of the operation
      CHARACTER(msgLen) :: ErrMsg      ! Error message if ErrStat /= ErrID_None
      integer(IntKi)    :: ioStatus

      call GetNewUnit ( input%UnitNum, ErrStat, ErrMsg )
      if (ErrStat >= AbortErrLev) then
         call iceErrorHndlr (iceLog, ErrStat, 'CountIceInputs => GetNewUnit: '//newLine//ErrMsg, 1)
         return
      endif

      call OpenFInpFile ( input%UnitNum, trim(fname), ErrStat, ErrMsg  )
      if (ErrStat >= AbortErrLev) then
         call iceErrorHndlr (iceLog, ErrStat, 'CountIceInputs => OpenFInpFile: '//newLine//ErrMsg, 1)
         return
      endif

      input%count = 0
      readLoop: do while (.true.)
         read(input%UnitNum,'(A)',iostat=ioStatus) lineIn
         if (ioStatus == iostat_end) exit readLoop
         if(lineIn(1:1) == '!' .or.  &   ! check for a comment
            lineIn(1:1) == '#' .or.  &
            lineIn(1:1) == '$' .or.  &
            lineIn(1:1) == '%' ) then
            cycle readLoop
         else
            input%count = input%count + 1
         endif
      enddo readLoop

      allocate(input%params(input%count), stat=errStat)
      if (errStat /= 0) then
         call iceErrorHndlr (iceLog, ErrID_Fatal, 'Error in input parameter array allocation in CountIceInputs', 1)
         return  !  error in array allocation
      endif

!   Go back to the beginning of the file
      rewind (unit=input%UnitNum, iostat = ErrStat)
      if (ErrStat /= 0) then
         call iceErrorHndlr (iceLog, ErrID_Fatal, 'CountIceInputs => OpenFInpFile: '//newLine//    &
                            ' Error rewinding parameter input file', 1)
         return
      endif

   end subroutine countIceInputs

!========================================================================
   subroutine readIceInputs (iceLog, input)
      type(iceFloe_LoggingType), intent(inout)   :: iceLog   ! structure with message and error logging variables
      type(iceInputType), intent(inout)            :: input

      character(1024)   :: lineIn
      character(132)    :: varName
      real(ReKi)        :: varValue
      integer(IntKi)    :: ioStatus, nCount, n
      INTEGER(IntKi)    :: ErrStat     ! Error status of the operation
      CHARACTER(msgLen) :: ErrMsg      ! Error message if ErrStat /= ErrID_None

!  read all of the records and put into the proper variable
      nCount = 0
      readLoop: do while (nCount < input%count)
         read(input%UnitNum,'(A)',iostat=ioStatus) lineIn
         if (ioStatus == iostat_end) exit readLoop
         if(lineIn(1:1) == '!' .or.  &   ! check for a comment
            lineIn(1:1) == '$' .or.  &
            lineIn(1:1) == '%' ) cycle readLoop

!   read in the variable name and its value, convert name to upper case
         read(lineIn,*) varName, varValue
         call Conv2UC (varName)
         nCount = nCount + 1
         input%params(nCount)%name = varName
         input%params(nCount)%value = varValue

!   loop through the list of possible input variables check for duplication
         searchLoop: do n = 1, nCount-1
            if(index(varname, input%params(n)%name) > 0) then
               call iceErrorHndlr (iceLog, ErrID_Warn, 'Input parameter '//trim(varName)//  &
                                                       ' has been specified twice.', 0)
               exit searchLoop
            endif
         enddo searchloop
      enddo readLoop

      close (input%UnitNum)

   end subroutine readIceInputs

!========================================================================
!  extract the input parameter value from the list
!  and send back to calling routine
   subroutine getRealInput (input, varName, outVal, iceLog, min, max)
      type(iceInputType), intent(in)             :: input
      type(iceFloe_LoggingType), intent(inout) :: iceLog   ! structure with message and error logging variables
      character(*), intent(in)         :: varName     ! This is the variable requested
      real(ReKi), intent(out)          :: outVal
      real(ReKi), intent(in), optional :: min, max

      integer(IntKi)                :: n
      character(len_trim(varName))  :: tmpName
      logical                       :: foundParam

      tmpName = varName  ! can't change varName since it's an intent(in)
      call Conv2UC (tmpName)

      foundParam = .false.
      do n = 1, input%count
         if(index(input%params(n)%name, tmpName) > 0) then
            outVal = input%params(n)%value
            foundParam = .true.
            exit
         endif
      enddo

      if (.not. foundParam)  then
         call iceErrorHndlr (iceLog, ErrID_Fatal, 'Input parameter '//varName//' was not in the parameter input file.', 1 )
         return
      endif

!  If requeseted, check that the value is withing reasonable bounds
      if (present(min)) then
         if (outVal < min) then
            call iceErrorHndlr (iceLog, ErrID_Warn, 'Input parameter '//varName//' is below the minimum bound', 0 )
         endif
      endif
      if (present(max)) then
         if (outVal > max) then
            call iceErrorHndlr (iceLog, ErrID_Warn, 'Input parameter '//varName//' is above the maximum bound', 0 )
         endif
      endif

   end subroutine getRealInput

!------------------------------------------------------------------
   subroutine getIntInput (input, varName, outVal, iceLog, min, max)
      type(iceInputType), intent(in)             :: input
      type(iceFloe_LoggingType), intent(inout) :: iceLog   ! structure with message and error logging variables
      character(*), intent(in)                :: varName     ! This is the variable requested
      integer(IntKi), intent(out)             :: outVal
      integer(IntKi), intent(in), optional    :: min, max

      integer(IntKi)                :: n
      character(len_trim(varName))  :: tmpInName
      logical                       :: foundParam

      tmpInName = varName  ! can't use varname since Conv2UC changes it
      call Conv2UC (tmpInName)

      foundParam = .false.
      do n = 1, input%count
         if(index(input%params(n)%name, tmpInName) > 0) then
            outVal = input%params(n)%value
            foundParam = .true.
            exit
         endif
      enddo

      if (.not. foundParam)  then
         call iceErrorHndlr (iceLog, ErrID_Fatal, 'Input parameter '//varName//' was not in the parameter input file.', 1 )
         return
      endif

!  If requeseted, check that the value is withing reasonable bounds
      if (present(min)) then
         if (outVal < min) then
            call iceErrorHndlr (iceLog, ErrID_Warn, 'Input parameter '//varName//' is below the minimum bound', 0 )
         endif
      endif
      if (present(max)) then
         if (outVal > max) then
            call iceErrorHndlr (iceLog, ErrID_Warn, 'Input parameter '//varName//' is above the maximum bound', 0 )
         endif
      endif

   end subroutine getIntInput

!------------------------------------------------------------------
   subroutine getLogicalInput (input, varName, outVal, iceLog)
      type(iceInputType), intent(in)             :: input
      type(iceFloe_LoggingType), intent(inout)   :: iceLog   ! structure with message and error logging variables
      character(*), intent(in)      :: varName     ! This is the variable requested
      logical, intent(out)          :: outVal
      character(len_trim(varName))  :: tmpName
      integer(IntKi)                :: n
      logical                       :: foundParam

      tmpName = varName
      call Conv2UC (tmpName)

      foundParam = .false.
      do n = 1, input%count
         if(index(input%params(n)%name, tmpName) > 0) then
            !outVal = input%params(n)%value ! warning #6192: Fortran 2003 does not allow this data type conversion. <- real to logical
            outVal = ( NINT(input%params(n)%value) /= 0 )  !Intel: The numeric value of .TRUE. and .FALSE. can be -1 and 0 or 1 and 0 depending on compiler option fpscomp
            foundParam = .true.
            exit
         endif
      enddo

      if (.not. foundParam)  then
         call iceErrorHndlr (iceLog, ErrID_Fatal, 'Input parameter '//varName//' was not in the parameter input file.', 1 )
         return
      endif
   end subroutine getLogicalInput

end module iceInputParams

max_stars_repo_name: abhineet-gupta/openfast

max_issues_repo_name: abhineet-gupta/openfast

max_forks_repo_name: abhineet-gupta/openfast