! File:         submodule_ejecta_generic_constructor.f90
! Authors:      Francesco Torsello (FT)
!************************************************************************
! Copyright (C) 2020-2023 Francesco Torsello                            *
!                                                                       *
! This file is part of SPHINCS_ID                                       *
!                                                                       *
! SPHINCS_ID is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by  *
! the Free Software Foundation, either version 3 of the License, or     *
! (at your option) any later version.                                   *
!                                                                       *
! SPHINCS_ID is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of        *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
! GNU General Public License for more details.                          *
!                                                                       *
! You should have received a copy of the GNU General Public License     *
! along with SPHINCS_ID. If not, see <https://www.gnu.org/licenses/>.   *
! The copy of the GNU General Public License should be in the file      *
! 'COPYING'.                                                            *
!************************************************************************

SUBMODULE (ejecta_generic) constructor

  !*********************************************************
  !
  !# Implementation of the constructor and
  !  destructor of TYPE [[ejecta]]
  !
  !  FT xx.11.2021
  !
  !*********************************************************


  IMPLICIT NONE


  CONTAINS


  !
  !-- Implementation of the constructor of the ejecta object
  !
  MODULE PROCEDURE construct_ejecta

    !****************************************************
    !
    !# Constructs an object of TYPE [[ejecta]]
    !  @todo to be OMP parallelized
    !
    !  FT xx.11.2021
    !
    !****************************************************

    USE constants, ONLY: pi
    USE utility,   ONLY: zero, one, two, four, ten, eos$poly, eos$pwpoly, &
                         sph_path
    USE NR,        ONLY: indexx
    USE pwp_EOS,   ONLY: get_Gamma0, get_Gamma1, get_Gamma2, get_Gamma3, &
                         get_K0, get_K1, get_K2, get_K3, get_p1, &
                         get_rho_0, get_rho_1, get_rho_2, select_EOS_parameters
    USE timing,    ONLY: timer

    IMPLICIT NONE


    INTEGER, PARAMETER:: unit_pos= 2589
    INTEGER, PARAMETER:: n_matter= 1
    DOUBLE PRECISION, PARAMETER:: atmosphere_density= 1.0439859633622731D-17

    INTEGER:: header_lines= 2 ! TODO: give this as input
    INTEGER:: nlines, ntmp
    INTEGER:: i_matter, itr, i, j, k
   ! INTEGER, DIMENSION(:), ALLOCATABLE:: x_sorted, y_sorted, z_sorted
    INTEGER, DIMENSION(:), ALLOCATABLE:: mass_profile_idx

    DOUBLE PRECISION:: xtmp, ytmp, ztmp, &
                       rhotmp, epstmp, vxtmp, vytmp, vztmp, &
                       dr, dphi, dth
    DOUBLE PRECISION, DIMENSION(:,:),   ALLOCATABLE:: grid_tmp
    !DOUBLE PRECISION, DIMENSION(:,:,:), ALLOCATABLE:: rho_tmp
    DOUBLE PRECISION, DIMENSION(:,:),   ALLOCATABLE:: mass_profile

    LOGICAL:: exist

    CHARACTER(LEN=:), ALLOCATABLE:: finalnamefile

    PRINT *, " * Reading ID on Cartesian, uniform grid from formatted file " &
             // TRIM(filename), "..."

    CALL derived_type% set_n_matter(n_matter)

    CALL derived_type% set_cold_system(.FALSE.)

    derived_type% construction_timer= timer( "ejecta_construction_timer" )

    CALL derived_type% construction_timer% start_timer()

    INQUIRE( FILE= TRIM(filename), EXIST= exist )

    IF( exist )THEN
      OPEN( UNIT= unit_pos, FILE= TRIM(filename), &
            FORM= "FORMATTED", ACTION= "READ", IOSTAT= ios, &
            IOMSG= err_msg )
      IF( ios > 0 )THEN
        PRINT *, "...error when opening " // TRIM(filename), &
                ". The error message is", err_msg
        STOP
      ENDIF
    ELSE
      PRINT *, "** ERROR! Unable to find file " // TRIM(filename)
      STOP
    ENDIF

    ! Get total number of lines in the file
    nlines = 0
    DO
      READ( unit_pos, * , IOSTAT= ios )
      IF ( ios /= 0 ) EXIT
      nlines = nlines + 1
    ENDDO

    CLOSE( UNIT= unit_pos )

    derived_type% n_gridpoints= nlines - header_lines

    ! Allocate the temporary array to store data
    ALLOCATE( grid_tmp( 2*derived_type% n_gridpoints, 8 ) )
    grid_tmp= zero

    ! Read the ID
    OPEN( UNIT= unit_pos, FILE= TRIM(filename), &
          FORM= "FORMATTED", ACTION= "READ" )

    ! Skip header
    DO itr= 1, header_lines, 1
      READ( unit_pos, * )
    ENDDO

    ! Read the data into the temporary array
    ntmp= 0
    DO itr= 1, derived_type% n_gridpoints, 1

      READ( UNIT= unit_pos, FMT= *, IOSTAT = ios, IOMSG= err_msg ) &
        xtmp, ytmp, ztmp, rhotmp, epstmp, vxtmp, vytmp, vztmp

      IF( ztmp > 0 )THEN
        ntmp= ntmp + 1
        grid_tmp( ntmp, 1 )= xtmp
        grid_tmp( ntmp, 2 )= ytmp
        grid_tmp( ntmp, 3 )= ztmp
        grid_tmp( ntmp, 4 )= rhotmp
        grid_tmp( ntmp, 5 )= epstmp
        grid_tmp( ntmp, 6 )= vxtmp
        grid_tmp( ntmp, 7 )= vytmp
        grid_tmp( ntmp, 8 )= vztmp
      ENDIF

      IF( ios > 0 )THEN
        PRINT *, "...error when reading " // TRIM(filename), &
                " at particle ", itr,". The status variable is ", ios, &
                ". The error message is", err_msg
        STOP
      ENDIF

    ENDDO
    derived_type% n_gridpoints= ntmp
    grid_tmp= grid_tmp( 1:2*derived_type% n_gridpoints, : )

    CLOSE( UNIT= unit_pos )

    DO itr= 1, SIZE(grid_tmp(:,1)), 1
      IF( grid_tmp(itr,1) > grid_tmp(1,1) )THEN
        derived_type% dx_grid= grid_tmp(itr,1) - grid_tmp(1,1)
        EXIT
      ENDIF
    ENDDO
    DO itr= 1, SIZE(grid_tmp(:,2)), 1
      IF( grid_tmp(itr,2) > grid_tmp(1,2) )THEN
        derived_type% dy_grid= grid_tmp(itr,2) - grid_tmp(1,2)
        EXIT
      ENDIF
    ENDDO
    DO itr= 1, SIZE(grid_tmp(:,3)), 1
      IF( grid_tmp(itr,3) > grid_tmp(1,3) )THEN
        derived_type% dz_grid= grid_tmp(itr,3) - grid_tmp(1,3)
        EXIT
      ENDIF
    ENDDO

    derived_type% xL_grid= MINVAL(grid_tmp( :, 1 ))
    derived_type% yL_grid= MINVAL(grid_tmp( :, 2 ))
    derived_type% zL_grid= MINVAL(grid_tmp( :, 3 ), grid_tmp( :, 3 ) /= 0 )
    derived_type% xR_grid= MAXVAL(grid_tmp( :, 1 ))
    derived_type% yR_grid= MAXVAL(grid_tmp( :, 2 ))
    derived_type% zR_grid= MAXVAL(grid_tmp( :, 3 ) )

    derived_type% nx_grid= NINT((MAXVAL(grid_tmp(:,1)) - derived_type% xL_grid)&
                                /derived_type% dx_grid + 1)
    derived_type% ny_grid= NINT((MAXVAL(grid_tmp(:,2)) - derived_type% yL_grid)&
                                /derived_type% dy_grid + 1 )
    derived_type% nz_grid= NINT((MAXVAL(grid_tmp(:,3)) - derived_type% zL_grid)&
                                /derived_type% dz_grid + 1 )

    PRINT *, " * ID on Cartesian, uniform grid read."
    PRINT *
    PRINT *, "** Grid information:"
    PRINT *
    PRINT *, " * Grid size in x direction:", &
             derived_type% xL_grid, derived_type% xR_grid
    PRINT *, " * Grid size in y direction:", &
             derived_type% yL_grid, derived_type% yR_grid
    PRINT *, " * Grid size in z direction:", &
             derived_type% zL_grid, derived_type% zR_grid
    PRINT *
    PRINT *, " * Grid spacing in the x direction:", derived_type% dx_grid
    PRINT *, " * Grid spacing in the y direction:", derived_type% dy_grid
    PRINT *, " * Grid spacing in the z direction:", derived_type% dz_grid
    PRINT *
    PRINT *, " * Number of grid points in the x direction:", &
             derived_type% nx_grid
    PRINT *, " * Number of grid points in the y direction:", &
             derived_type% ny_grid
    PRINT *, " * Number of grid points in the z direction:", &
             derived_type% nz_grid
    PRINT *

    PRINT *, "** Checking that the grid dimensions are consistent with the " &
             //"number of lines in the file..."
    ! Check that the grid dimensions are consistent
    IF( derived_type% nx_grid*derived_type% ny_grid*derived_type% nz_grid &
        /= derived_type% n_gridpoints )THEN

      PRINT *, derived_type% nx_grid
      PRINT *, derived_type% ny_grid
      PRINT *, derived_type% nz_grid
      PRINT *, derived_type% nx_grid*derived_type% ny_grid*derived_type% nz_grid
      PRINT *, derived_type% n_gridpoints
      STOP

    ENDIF

    PRINT *, "** Checking that the number of grid points, the grid spacings " &
             //"and the grid sizes are consistent..."
    ! Check that nx dx and the grid extent are consistent
    ztmp= derived_type% zL_grid
    DO k= 1, derived_type% nz_grid - 1, 1
      ztmp= ztmp + derived_type% dz_grid
    ENDDO
    IF( ABS( ztmp - MAXVAL(grid_tmp( :, 3 )) ) &
        > derived_type% dz_grid/1.0D+6 )THEN
      PRINT *, "** ERROR! ztmp=", ztmp
      PRINT *, "          zR_grid=", MAXVAL(grid_tmp( :, 3 ))
      STOP
    ENDIF
    ytmp= derived_type% yL_grid
    DO j= 1, derived_type% ny_grid - 1, 1
      ytmp= ytmp + derived_type% dy_grid
    ENDDO
    IF( ABS( ytmp - MAXVAL(grid_tmp( :, 2 )) ) &
        > derived_type% dz_grid/1.0D+6 )THEN
      PRINT *, "** ERROR! ytmp=", ytmp
      PRINT *, "          yR_grid=", MAXVAL(grid_tmp( :, 2 ))
      STOP
    ENDIF
    xtmp= derived_type% xL_grid
    DO i= 1, derived_type% nx_grid - 1, 1
      xtmp= xtmp + derived_type% dx_grid
    ENDDO
    IF( ABS( xtmp - MAXVAL(grid_tmp( :, 1 )) ) &
        > derived_type% dx_grid/1.0D+6 )THEN
      PRINT *, "** ERROR! xtmp=", xtmp
      PRINT *, "          xR_grid=", MAXVAL(grid_tmp( :, 1 ))
      STOP
    ENDIF
    PRINT *

    ! Allocate and initialize member arrays
    CALL derived_type% allocate_gridid_memory( n_matter )

    derived_type% grid= zero
    derived_type% baryon_mass_density= zero
    derived_type% specific_energy= zero
    derived_type% vel= zero

    ! Store the ID into the member arrays
    DO i= 1, derived_type% nx_grid, 1
      DO j= 1, derived_type% ny_grid, 1
        DO k= 1, derived_type% nz_grid, 1

          derived_type% grid( i, j, k, 1 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 1 )
          derived_type% grid( i, j, k, 2 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 2 )
          derived_type% grid( i, j, k, 3 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 3 )

          derived_type% baryon_mass_density( i, j, k )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 4 )

          derived_type% specific_energy( i, j, k )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 5 )

          derived_type% vel( i, j, k, 1 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 6 )
          derived_type% vel( i, j, k, 2 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 7 )
          derived_type% vel( i, j, k, 3 )= &
            grid_tmp( (i-1)*(derived_type% ny_grid)*(derived_type% nz_grid) &
                      + (j-1)*(derived_type% nz_grid) + k, 8 )

        ENDDO
      ENDDO
    ENDDO

    ! Get rid of the atmosphere coming from a mesh-based simulation, if present
    DO i= 1, derived_type% nx_grid, 1
      DO j= 1, derived_type% ny_grid, 1
        DO k= 1, derived_type% nz_grid, 1

          IF( derived_type% baryon_mass_density( i, j, k ) &
              <= atmosphere_density )THEN

            derived_type% baryon_mass_density( i, j, k )= zero

            derived_type% specific_energy( i, j, k )= zero

            derived_type% vel( i, j, k, : )= zero

          ENDIF

          !derived_type% baryon_mass_density( i, j, k )= &
          !  MAX( zero, &
          !  derived_type% baryon_mass_density( i, j, k ) - atmosphere_density )
          !
          !IF( derived_type% baryon_mass_density( i, j, k ) == zero )THEN
          !
          !  derived_type% specific_energy( i, j, k )= zero
          !
          !  derived_type% vel( i, j, k, : )= zero
          !
          !ENDIF

        ENDDO
      ENDDO
    ENDDO


    ! Assign ID properties to member arrays
    DO i_matter= 1, n_matter, 1

      !derived_type% masses(i_matter)= zero
      derived_type% centers(i_matter,:)= zero
      derived_type% barycenters(i_matter,:)= zero
      derived_type% sizes(i_matter,:)= [ &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 2 )))**two ), &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 2 )))**two ), &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 2 )))**two ), &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 2 )))**two ), &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 3 )))**two ), &
                  !1.3D0*SQRT( ABS(MAXVAL(grid_tmp( :, 1 )))**two &
                  !    + ABS(MAXVAL(grid_tmp( :, 3 )))**two ) ]
                                         ABS(derived_type% xL_grid), &
                                         ABS(MAXVAL(grid_tmp( :, 1 ))), &
                                         ABS(derived_type% yL_grid), &
                                         ABS(MAXVAL(grid_tmp( :, 2 ))), &
                                         ABS(MAXVAL(grid_tmp( :, 3 ))), &
                                         ABS(MAXVAL(grid_tmp( :, 3 ))) ]

    ENDDO

    finalnamefile= TRIM(sph_path)//"pos_ejecta.dat"

    INQUIRE( FILE= TRIM(finalnamefile), EXIST= exist )

    IF( exist )THEN
        OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "REPLACE", &
              FORM= "FORMATTED", &
              POSITION= "REWIND", ACTION= "WRITE", IOSTAT= ios, &
              IOMSG= err_msg )
    ELSE
        OPEN( UNIT= 2, FILE= TRIM(finalnamefile), STATUS= "NEW", &
              FORM= "FORMATTED", &
              ACTION= "WRITE", IOSTAT= ios, IOMSG= err_msg )
    ENDIF
    IF( ios > 0 )THEN
      PRINT *, "...error when opening " // TRIM(finalnamefile), &
               ". The error message is", err_msg
      STOP
    ENDIF

    DO i= 1, derived_type% nx_grid - 1, 1
      DO j= 1, derived_type% ny_grid - 1, 1
        DO k= 1, derived_type% nz_grid - 1, 1

          WRITE( UNIT = 2, IOSTAT = ios, IOMSG = err_msg, FMT = * ) &
            derived_type% grid( i, j, k, 1 ), &
            derived_type% grid( i, j, k, 2 ), &
            derived_type% grid( i, j, k, 3 ), &
            derived_type% baryon_mass_density( i, j, k ), &
            derived_type% read_mass_density( &
              derived_type% grid( i, j, k, 1 ) + derived_type% dx_grid/two, &
              derived_type% grid( i, j, k, 2 ), &
              derived_type% grid( i, j, k, 3 ) ), &
            derived_type% grid( i, j, k, 1 ) + derived_type% dx_grid/two, &
            derived_type% specific_energy( i, j, k )
        ENDDO
      ENDDO
    ENDDO

    CLOSE( UNIT= 2 )

    ! Assign total mass to member variable
    dr             = derived_type% dx_grid/four
    dth            = pi/two/ten*ten
    dphi           = two*pi/ten*ten

    ALLOCATE( mass_profile( 3, 0:NINT(ABS(MAXVAL(grid_tmp( :, 1 )))/dr) ), &
              STAT= ios, ERRMSG= err_msg )
    ALLOCATE( mass_profile_idx( 0:NINT(ABS(MAXVAL(grid_tmp( :, 1 )))/dr) ), &
              STAT= ios, ERRMSG= err_msg )

    CALL derived_type% integrate_baryon_mass_density( &
                            derived_type% centers(1,:), &
                            MAXVAL(grid_tmp( :, 1 )), &
                            derived_type% read_mass_density( &
                              derived_type% centers(1,1), &
                              derived_type% centers(1,2), &
                              derived_type% centers(1,3) ), &
                            dr, dth, dphi, &
                            derived_type% masses(1), mass_profile, &
                            mass_profile_idx )!, &
                            !radii= [MAXVAL(derived_type% sizes(1,3:4)), &
                            !        MAXVAL(derived_type% sizes(1,5:6))] )

    ! Set the EOS parameters
    derived_type% eos_id= eos$pwpoly

    CALL select_EOS_parameters("APR4")

    derived_type% npeos  = 3
    derived_type% gamma0 = get_Gamma0()
    derived_type% gamma1 = get_Gamma1()
    derived_type% gamma2 = get_Gamma2()
    derived_type% gamma3 = get_Gamma3()
    derived_type% kappa0 = get_K0()
    derived_type% kappa1 = get_K1()
    derived_type% kappa2 = get_K2()
    derived_type% kappa3 = get_K3()
    derived_type% logP1  = LOG10(get_p1())
    derived_type% logRho0= LOG10(get_rho_0())
    derived_type% logRho1= LOG10(get_rho_1())
    derived_type% logRho2= LOG10(get_rho_2())

    derived_type% finalize_sph_id_ptr => finalize

    CALL derived_type% construction_timer% stop_timer()

  END PROCEDURE construct_ejecta


  MODULE PROCEDURE finalize

    !***********************************************
    !
    !# This SUBROUTINE is curretly just a placeholder.
    !  It could be used, for exmaple, to correct
    !  the ADM momentum at the end of the execution,
    !  or correct for residual eccentricity, etc.
    !
    !  @note Temporary implementation, to avoid warnings
    !        about unused variables.
    !
    !  FT 14.04.2022
    !
    !***********************************************

    IMPLICIT NONE

    ! Temporary implementation, to avoid warnings about unused variables

    pos  = pos
    nlrf = nlrf
    nu   = nu
    pr   = pr
    vel_u= vel_u
    theta= theta
    nstar= nstar
    u    = u

  END PROCEDURE finalize


  MODULE PROCEDURE nothing

    !***********************************************
    !
    !# Procedure that does nothing. It is used to instantiate a deferred
    !  idbase procedure which s not needed in TYPE [[ejecta_generic]].
    !  It also serves as a placeholder in case the idbase procedure
    !  will be needed in the future.
    !
    !  FT 15.09.2022
    !
    !***********************************************

    IMPLICIT NONE

  END PROCEDURE nothing


  !
  !-- Implementation of the destructor of the bns object
  !
  MODULE PROCEDURE destruct_ejecta

    !****************************************************
    !
    !# Destructs an object of TYPE [[ejecta]]
    !
    !  FT xx.11.2021
    !
    !****************************************************

    IMPLICIT NONE

    CALL this% deallocate_gridid_memory()


  END PROCEDURE destruct_ejecta


END SUBMODULE constructor