MODULE comtraj
  !!======================================================================
  !!                   ***  MODULE comtraj  ***
  !!
  !!======================================================================
  USE parameters, ONLY       : rsh,rlg,lchain

  IMPLICIT NONE
  PRIVATE

  !! * Accessibility
  PUBLIC    :: patch_list_append
  PUBLIC    :: get_patch
  PUBLIC    :: enlarge_patch, resize_patch
#ifdef key_MPI_2D
  PUBLIC    :: init_mpi_type_particle
#endif

  !----------------------------------------
  ! TYPE type_particle

  TYPE, PUBLIC :: type_particle

    LOGICAL :: active = .FALSE. ! .True.  if particle actually active
    
#ifdef key_MPI_2D
    ! MPI managing
    !-----------------
    INTEGER :: limitbye = 0     ! Specify the direction of the boundary crossing
                                ! 1=SW, 2=S, 3=SE, 4=E, 5=NE, 6=N, 7=NW, 8=N
#endif

    INTEGER :: itypevert        ! type of trajectory (z=cst, random walk...)
    INTEGER :: num = 0          ! to keep track of particles when save/restart 

    ! main properties
    !-----------------
    REAL(KIND=rlg)  :: date_orig   ! date of release

    ! location
    !-----------------
    REAL(KIND=rsh)  :: xpos, ypos
    REAL(KIND=rsh)  :: spos, zpos
    REAL(KIND=rsh)  :: d3, h0, xe, hc
    REAL(KIND=rsh)  :: flag = 0.0_rsh   ! wet-drying flag

#ifdef key_ibm
    ! Biological parameters
    INTEGER         :: stage
    REAL(KIND=rsh)  :: w, size, density, Drate
    REAL(KIND=rsh)  :: temp
    REAL(KIND=rsh)  :: super ! one part represent number individuals (superindividual)
#ifdef key_ibm_species
    ! DEB state variables and parameters (with default value)
    INTEGER         :: dayjuv, dayspawn, yearspawn
    INTEGER         :: season = 0
    INTEGER         :: hmove  = 0
    REAL(KIND=rsh)  :: E
    REAL(KIND=rsh)  :: L
    REAL(KIND=rsh)  :: H
    REAL(KIND=rsh)  :: R
    REAL(KIND=rsh)  :: Gam    = 0.0_rsh
    REAL(KIND=rsh)  :: X      = 0.0_rsh
    REAL(KIND=rsh)  :: Wdeb   = 0.0_rsh
    REAL(KIND=rsh)  :: Neggs  = 0.0_rsh
    REAL(KIND=rsh)  :: NRJ    = 0.0_rsh
    REAL(KIND=rsh)  :: denspawn
#endif
#endif

  END TYPE type_particle

  !----------------------------------------
  ! TYPE type_patch

  TYPE, PUBLIC :: type_patch

    INTEGER               :: id = -1
#ifdef key_ibm
    INTEGER               :: parent_id = -2
    INTEGER               :: generation = 0
#endif

    INTEGER               :: nb_part_alloc = 0      ! Size of allocated data array for particles
    INTEGER               :: nb_part_total = 0      ! Sum of particles over all MPI domains
    INTEGER               :: nb_part_batch = 10     ! Size of batch for new allocations
    INTEGER               :: nb_part_max   = -1     ! Maximum allowed number of particles
    
    REAL(KIND=rlg)        :: t_beg, t_end
    REAL(KIND=rlg)        :: t_save,  dt_save
#ifdef key_ibm
    REAL(KIND=rlg)        :: t_spawn, dt_spawn
#endif
    CHARACTER(LEN=lchain) :: file_inp               ! Input data file name
    CHARACTER(LEN=lchain) :: file_out               ! Output NetCDF file name
    LOGICAL               :: file_out_init = .FALSE.
#ifdef key_ibm
    CHARACTER(LEN=lchain) :: file_out_base          ! Output NetCDF base file name
    LOGICAL               :: has_spawned = .FALSE.
#endif
#ifdef key_ibm_species
    CHARACTER(LEN=lchain) :: species          ! name of species
#endif
    TYPE(type_particle)                             :: init_particle    ! Init values used for new particles
    TYPE(type_particle), ALLOCATABLE, DIMENSION(:)  :: particles
    TYPE(type_patch), POINTER                       :: next => NULL()   ! Next patch in the list

  END TYPE type_patch

  !----------------------------------------
  ! TYPE type_patch_list

  TYPE, PUBLIC :: type_patch_list
      INTEGER                   :: nb = 0           ! Number of elements in the list
      TYPE(type_patch), POINTER :: first => NULL()  ! First patch in the list
      TYPE(type_patch), POINTER :: last  => NULL()  ! Last  patch inserted in the list
  END TYPE type_patch_list
  
  !----------------------------------------
  !! * Shared module variables
  
  TYPE(type_patch_list), PUBLIC :: patches
                  
  CHARACTER(LEN=lchain), PUBLIC :: file_trajec      ! name of configuration file
  INTEGER,               PUBLIC :: itypetraj        ! initialisation type (circle, rectangle,netcdf)
  INTEGER,               PUBLIC :: ndtz             ! number of time step division for vertical subloop

#ifdef key_MPI_2D
  INTEGER,               PUBLIC :: type_mpi_particle
  INTEGER,               PUBLIC :: down_give, up_give, right_give, left_give
#endif

CONTAINS

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

#ifdef key_MPI_2D

  SUBROUTINE init_mpi_type_particle
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  init_mpi_type_particle  ***
   !&E
   !&E ** Purpose : Defines the MPI type for 'type_particle'
   !&E
   !&E       !  2014-12 (M. Honnorat) Refactor MPI routines
   !&E---------------------------------------------------------------------
   !! * Modules used

   USE comvars2d, ONLY : type_mpi_rlg, type_mpi_rsh

   IMPLICIT NONE
   include 'mpif.h'

   !! * Local declarations
#ifdef key_ibm
#ifdef key_ibm_species
   INTEGER, PARAMETER   :: nb = 36     ! key_ibm_species and key_ibm
#else
   INTEGER, PARAMETER   :: nb = 21     ! key_ibm only
#endif
#else
   INTEGER, PARAMETER   :: nb = 14     ! key_MPI_2D only
#endif

   INTEGER,                        DIMENSION(nb) :: old_types, block_lengths
   INTEGER(kind=MPI_ADDRESS_KIND), DIMENSION(nb) :: displacements, addresses
   INTEGER              :: ierr_mpi, i
   TYPE(type_particle)  :: particle
   !!----------------------------------------------------------------------
   !! * Executable part

   ! Create MPI type for particles 
   old_types = (/                                                               &
       MPI_LOGICAL, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER,                      &
       type_mpi_rlg, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh,    &
       type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh     &
#ifdef key_ibm
     , MPI_INTEGER, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh,     &
       type_mpi_rsh, type_mpi_rsh                                               &
#ifdef key_ibm_species
     , MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER, MPI_INTEGER,         &
       type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh,    &
       type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh, type_mpi_rsh     &
#endif
#endif
   /)
   block_lengths = (/                                                           &
       1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1                                 &
#ifdef key_ibm
     , 1, 1, 1, 1, 1, 1, 1                                                      &
#ifdef key_ibm_species
     , 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1                              &
#endif
#endif
   /)

   i = 1
   CALL MPI_GET_ADDRESS(particle % active,    addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % limitbye,  addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % itypevert, addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % num,       addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % date_orig, addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % xpos,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % ypos,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % spos,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % zpos,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % d3,        addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % h0,        addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % xe,        addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % hc,        addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % flag,      addresses(i), ierr_mpi) ; i = i+1
#ifdef key_ibm
   CALL MPI_GET_ADDRESS(particle % stage,     addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % w,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % size,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % density,   addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % Drate,     addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % temp,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % super,     addresses(i), ierr_mpi) ; i = i+1
#ifdef key_ibm_species
   CALL MPI_GET_ADDRESS(particle % dayjuv,    addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % dayspawn,  addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % yearspawn, addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % season,    addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % hmove,     addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % E,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % L,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % H,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % R,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % Gam,       addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % X,         addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % Wdeb,      addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % Neggs,     addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % NRJ,       addresses(i), ierr_mpi) ; i = i+1
   CALL MPI_GET_ADDRESS(particle % denspawn,  addresses(i), ierr_mpi) ; i = i+1
#endif
#endif

   displacements(:) = addresses(:) - addresses(1)
    
   CALL MPI_TYPE_CREATE_STRUCT(nb, block_lengths, displacements, old_types, type_mpi_particle, ierr_mpi)
   CALL MPI_TYPE_COMMIT(type_mpi_particle, ierr_mpi)

  END SUBROUTINE init_mpi_type_particle

#endif
  
  !!======================================================================

  FUNCTION patch_list_append ( patch_list ) RESULT(new_patch)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  patch_list_append  ***
   !&E
   !&E ** Purpose : Append a new patch to 'patch_list'. A pointer to the
   !&E              newly created patch is returned.
   !&E
   !&E       !  2014-12 (M. Honnorat) Refactor MPI routines
   !&E---------------------------------------------------------------------
   !! * Arguments
   TYPE(type_patch_list), INTENT(inout) :: patch_list

   !! * Local declarations
   TYPE(type_patch), POINTER   :: new_patch
   INTEGER, SAVE               :: id = 0

   !!----------------------------------------------------------------------
   !! * Executable part

   ALLOCATE( new_patch )

   id = id+1
   new_patch % id = id

   IF ( ASSOCIATED(patch_list % last) ) then
     ! the list is not empty, append 'new_patch' after the last one.
     patch_list % last % next => new_patch
   ELSE
     ! the list has just been initialized. Let 'new_patch' be the first one.
     patch_list % first => new_patch
   END IF

   ! anyway, for next time 'new_patch' will be the last one. 
   patch_list % last => new_patch
   patch_list % nb = patch_list % nb + 1

  END FUNCTION patch_list_append

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

  FUNCTION get_patch ( patch_list, n ) RESULT(patch)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE  get_patch  ***
   !&E
   !&E ** Purpose : Returns a pointer to the nth patch from 'patch_list'.
   !&E
   !&E       !  2014-12 (M. Honnorat) Refactor MPI routines
   !&E---------------------------------------------------------------------
   !! * Arguments
   TYPE(type_patch_list), INTENT(in)    :: patch_list
   INTEGER,               INTENT(in)    :: n

   !! * Local declarations
   TYPE(type_patch), POINTER    :: patch
   INTEGER                      :: i

   !!----------------------------------------------------------------------
   !! * Executable part

   patch => NULL()

   IF ( (patch_list % nb <= 0) .OR. (n > patch_list % nb) )  RETURN
   
   patch => patch_list % first

   DO i = 1, patch_list % nb

      IF (i == n) RETURN
      patch => patch % next
   
   END DO

  END FUNCTION get_patch

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

  SUBROUTINE enlarge_patch(patch, n)

  !&E---------------------------------------------------------------------
  !&E                 ***  ROUTINE enlarge_patch  ***
  !&E
  !&E ** Purpose : Reallocate array patch % particles so that it can store
  !&E              n additional particles. The size n may be negative.
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : ex_traj_1d
  !&E
  !&E ** External calls :
  !&E
  !&E ** History :
  !&E       !  2015-01  (M. Honnorat) IBM upgrade
  !&E---------------------------------------------------------------------
  !! * Modules used

  !! * Arguments
  TYPE(type_patch), INTENT(inout)  :: patch
  INTEGER,          INTENT(in)     :: n     ! number of additional particles

  !! * Local declarations
  INTEGER      :: old_size, new_size
  INTEGER      :: i, j
  TYPE(type_particle), ALLOCATABLE, DIMENSION(:) :: tmp_array

  !!----------------------------------------------------------------------
  !! * Executable part

  ! Compute new_size
  new_size = patch % nb_part_alloc + patch % nb_part_batch * CEILING(REAL(n)/patch % nb_part_batch)
  old_size = patch % nb_part_alloc

  ! Save old data in resized array
  ALLOCATE(tmp_array(new_size))
  IF ( new_size > old_size ) THEN
    tmp_array(1:old_size) = patch % particles
    tmp_array(old_size+1:new_size) = patch % init_particle
  ELSE
     print*,'FIXME: enlarge_patch : NOT TESTED new_size < old_size', new_size, old_size
     i = 1
     DO j = 1,old_size
        IF ( patch % particles(j) % active ) THEN 
           tmp_array(i) = patch % particles(j)
           i = i+1
        END IF
     END DO
  END IF

  ! Update patch data with new array
  CALL MOVE_ALLOC(TO=patch % particles, FROM=tmp_array)
  patch % nb_part_alloc = new_size

  END SUBROUTINE enlarge_patch

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

  SUBROUTINE resize_patch(patch, required_size)

  !&E---------------------------------------------------------------------
  !&E                 ***  ROUTINE resize_patch  ***
  !&E
  !&E ** Purpose : Resize array patch % particles so that it can store
  !&E              #new_size# particles.
  !&E
  !&E ** Description :
  !&E
  !&E ** Called by : ibm_3d
  !&E
  !&E ** External calls :
  !&E
  !&E ** History :
  !&E       !  2015-01  (M. Honnorat) IBM upgrade
  !&E---------------------------------------------------------------------
  !! * Modules used

  !! * Arguments
  TYPE(type_patch), INTENT(inout)  :: patch
  INTEGER,          INTENT(in)     :: required_size

  !! * Local declarations
  INTEGER                                        :: old_size, new_size
  TYPE(type_particle), ALLOCATABLE, DIMENSION(:) :: tmp_array

  !!----------------------------------------------------------------------
  !! * Executable part

  old_size = patch % nb_part_alloc
  new_size = patch % nb_part_batch * CEILING(REAL(required_size)/patch % nb_part_batch)
  
  IF ( old_size >= new_size ) RETURN     ! Nothing to do
  
  ! Save old data in resized array
  ALLOCATE(tmp_array(new_size))
  tmp_array(1:old_size) = patch % particles
  tmp_array(old_size+1:new_size) = patch % init_particle

  ! Update patch data with new array
  CALL MOVE_ALLOC(TO=patch % particles, FROM=tmp_array)
  patch % nb_part_alloc = new_size

  END SUBROUTINE resize_patch

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

END MODULE
