!---------------------------------------------------------------------------
! pmlib_remesh.f
!   version: 2D single-core
!---------------------------------------------------------------------------
!                  __    __
!    ___  _ _     / /__ / /
!   | _ \/ ` \ _ / // // _ \
!   | ,_/_,_,_\  \_/\_/\___/
!   |_|                      Mads M. Hejlesen* and Jens H. Walther*^
!                            * Technical University of Denmark
!                            ^ ETH-Zurich
!
!---------------------------------------------------------------------------
! Copyright (c) 2014 Technical University of Denmark, 
!                    Department of Mechanical Engineering
!
! This program 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 2 of the License, or
! (at your option) any later version.
! 
! This program 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 this program; if not, see <http://www.gnu.org/licenses/>
!
! Contact: mmhej@mek.dtu.dk or jhwa@mek.dtu.dk
!---------------------------------------------------------------------------
SUBROUTINE pmlib_remesh(patch,topo,mesh,part,min_vort,min_dvort,ierr,dvort,vel)

USE pmlib_mod_mesh
USE pmlib_mod_particles
USE pmlib_mod_patch
USE pmlib_mod_topology

IMPLICIT NONE

!---------------------------------------------------------------------------
! Arguments
!---------------------------------------------------------------------------
  TYPE(class_patch)                         :: patch
  TYPE(class_topology)                      :: topo
  TYPE(class_mesh)                          :: mesh
  TYPE(class_particles)                     :: part
  REAL(MK), INTENT(IN)                      :: min_vort
  REAL(MK), INTENT(IN)                      :: min_dvort
  INTEGER, INTENT(OUT)                      :: ierr
  LOGICAL, OPTIONAL                         :: dvort
  LOGICAL, OPTIONAL                         :: vel

!---------------------------------------------------------------------------
! Local variables
!---------------------------------------------------------------------------
  REAL(MK)                                :: px, py
  REAL(MK)                                :: strength, dstrength
  INTEGER                                 :: i, j, k, l
  INTEGER                                 :: inp
  REAL(MK),DIMENSION(pmlib_ndim)          :: xmin, dx
  INTEGER,DIMENSION(pmlib_ndim)           :: ncell

!---------------------------------------------------------------------------
! Initiate subroutine
!---------------------------------------------------------------------------
  ierr = 0

  dx     = topo%dx
  ncell  = topo%ncell
  xmin   = topo%xmin

!---------------------------------------------------------------------------
! Deallocate the particles if they are already allocated
!---------------------------------------------------------------------------
  IF(ASSOCIATED( part%pos )) DEALLOCATE(part%pos,stat=ierr)

  IF(ASSOCIATED( part%vel )) DEALLOCATE(part%vel,stat=ierr)

  IF(ASSOCIATED( part%vort )) DEALLOCATE(part%vort,stat=ierr)

  IF(ASSOCIATED( part%dvort )) DEALLOCATE(part%dvort,stat=ierr)

  IF(ASSOCIATED( part%vel_rk )) DEALLOCATE(part%vel_rk,stat=ierr)

  IF(ASSOCIATED( part%dvort_rk )) DEALLOCATE(part%dvort_rk,stat=ierr)

  IF(ierr .NE. 0) THEN
    CALL pmlib_write( 'pmlib_remesh_particles', &
                    & 'Failed to deallocate existing particles.')
    GOTO 9999
  ENDIF

!---------------------------------------------------------------------------
! Find the number of particles
!---------------------------------------------------------------------------
  inp = 0
    DO j = 1, ncell(2)
      DO i = 1, ncell(1)

        strength  = ABS( mesh%vort(1,i,j) )
        dstrength = ABS( mesh%dvort(1,i,j) )

        IF( strength .GT. min_vort .OR. &
          & dstrength .GT. min_dvort) THEN

          inp = inp + 1
        END IF

      END DO
    END DO

  topo%npart = inp

!---------------------------------------------------------------------------
! Allocate particles
!---------------------------------------------------------------------------
  ALLOCATE( part%pos(pmlib_ndim,inp), part%vort(pmlib_nvort,inp), &
          & part%vel(pmlib_ndim,inp), part%dvort(pmlib_nvort,inp) )

  ALLOCATE( part%vel_rk(pmlib_ndim,inp), &
          & part%dvort_rk(pmlib_nvort,inp) )

!---------------------------------------------------------------------------
! Create particles
!---------------------------------------------------------------------------
  inp = 0

    DO j = 1, ncell(2)
      py = xmin(2) + (REAL(j - 1,MK) + 0.5_MK) *dx(2)

      DO i = 1, ncell(1)
        px = xmin(1) + (REAL(i - 1,MK) + 0.5_MK) *dx(1)

        strength  = ABS( mesh%vort(1,i,j) )
        dstrength = ABS( mesh%dvort(1,i,j) )

        IF( strength .GT. min_vort .OR. &
          & dstrength .GT. min_dvort) THEN

          inp = inp + 1

          part%pos(1,inp) = px
          part%pos(2,inp) = py

          IF( PRESENT(vel) .AND. vel )THEN
            part%vel(1:pmlib_ndim,inp) = mesh%vel(1:pmlib_ndim,i,j)
          ELSE
            part%vel(1:pmlib_ndim,inp) = 0.0_MK
          END IF

          part%vel_rk(1:pmlib_ndim,inp) = 0.0_MK

          part%vort(1:pmlib_nvort,inp) = mesh%vort(1:pmlib_nvort,i,j)

          IF( PRESENT(dvort) .AND. dvort )THEN
            part%dvort(1:pmlib_nvort,inp) = mesh%dvort(1:pmlib_nvort,i,j)
          ELSE
            part%dvort(1:pmlib_nvort,inp) = 0.0_MK
          END IF
          part%dvort_rk(1:pmlib_nvort,inp) = 0.0_MK

        END IF

      END DO
    END DO

!---------------------------------------------------------------------------
! Return
!---------------------------------------------------------------------------
 9999 CONTINUE
RETURN
END SUBROUTINE pmlib_remesh
