MODULE ibmmove
#if defined key_ibm
   !!======================================================================
   !!                   ***  MODULE movement   ***
   !! 
   !!         M. Huret 03/2012
   !!         Movement for fish IBM
   !!         Several movement rules are implemented
   !!        
   !&E
   !&E ** History :
   !&E
   !&E
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY :  imin,imax,jmin,jmax,kmax,ljmin,ljmax,ljmaxp1,  &
                           limaxp1,liminm1,limaxp2,ljminm1,ljmaxp2,       &
                           rsh,rlg,lchain,limin,liminm2,ljminm2,limax,    &
                           valmanq
   USE comtraj

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC :: fish_move_init         ! routine called by ibm_init
   PUBLIC :: fish_move              ! routine called by ibm

   !! * Shared module variables
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:,:), PUBLIC :: fish

   !! * Private variables
   INTEGER,        PARAMETER    :: nsize   = 36
   REAL(KIND=rsh), PARAMETER    :: sizemin = 3.5_rsh
   CHARACTER(LEN=lchain)        :: file_fish
   CHARACTER(LEN=lchain)        :: name_in_fish

CONTAINS
   
   !!======================================================================

   SUBROUTINE fish_move_init

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE fish_move_init  ***
   !&E
   !&E ** Purpose : Read fish distribution file
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : 
   !&E
   !&E ** External calls :ionc4_openr, ionc4_read_time, ionc4_read
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_read_subxyt, ionc4_read_subzxyt
   !&E                     ionc4_close
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2012-03 (M. Huret) Original code
   !&E
   !&E---------------------------------------------------------------------
   
   !! * Modules used

   USE ionc4, ONLY : ionc4_openr,ionc4_read_dimt,ionc4_read_time,&
                     ionc4_read_subxyt,ionc4_read_subzxyt,ionc4_close

   !! * Arguments
 

   !! * Local declarations

   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:) :: fish1
   INTEGER :: idimt,t

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

   file_fish = '../../inputs/probadistrib_1_2.nc'
   name_in_fish = 'Probability'

   CALL ionc4_openr(file_fish,l_out_nc4par=.true.)

   idimt = ionc4_read_dimt(file_fish) ! read number of population distribution along the year

   ALLOCATE(fish (nsize,liminm2:limaxp2,ljminm2:ljmaxp2,idimt)) 
   ALLOCATE(fish1(nsize,liminm2:limaxp2,ljminm2:ljmaxp2)) 
   
   DO t= 1,idimt
      CALL ionc4_read_subzxyt(file_fish,TRIM(name_in_fish),fish1,liminm2,limaxp2,ljminm2,ljmaxp2,1,nsize,t,0,0,1)
      fish(:,liminm2:limaxp2,ljminm2:ljmaxp2,t) = fish1(:,:,:)
   END DO

   CALL ionc4_close(file_fish)

  ! print*, 'fish', fish(nsize,limin:limax,ljmin:ljmax) 
   
   DEALLOCATE(fish1)

   END SUBROUTINE fish_move_init

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

   SUBROUTINE fish_move(particle)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE fish_move  ***
   !&E
   !&E ** Purpose : routine for random walk within a given distribution area
   !&E              following a Metropolis algorithm  
   !&E
   !&E ** Description :
   !&E ** Called by : ibm
   !&E ** External calls :
   !&E ** Reference :
   !&E ** History :
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d,   ONLY : h0_g,t,dx,dy
   USE comtraj,     ONLY : type_particle
   USE trajectools, ONLY : ksupkinf

   !! * Arguments
   TYPE(type_particle), INTENT(inout)   :: particle

   !! * Local declarations
   LOGICAL              :: move
   CHARACTER(len=19)    :: tool_sectodat
   REAL(rsh)            :: v1,v2,Pi,Pj,harvest,xpos_n,ypos_n,cell
   INTEGER              :: index,icell,jcell,ipos,jpos,saison,icells,jcells
   INTEGER              :: jj,mm,aaaa,hh,minu,sec

        INTEGER, PARAMETER :: n = 1     ! FIXME : a l'origine, c'etait le numero du patch... (mh-2015)

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

   PRINT_DBG*, 'enter FISH_MOVE routine'

   !Get the index of size in the probability distribution matrix
   index = min(max(1,NINT(particle % size*2)-INT(sizemin*2)+1),34) ! 34=20cm, au dela rien a l automne
   icells = NINT(particle % xpos)
   jcells = NINT(particle % ypos)

   !Get the seasonal distribution
   CALL tool_decompdate(tool_sectodat(t),jj,mm,aaaa,hh,minu,sec)

   ! when transition is let free, it takes 2 to 3 months to come to a reasonable R2 for comparison
   ! between obs and model distribution by size (March-April and July-August)

   if (n==1) then ! normal migration
      saison=1
      if(mm <= 2 .OR. mm >= 7) THEN
         saison=2
      endif
      if(mm <= 2  .OR. mm >= 10) then 
         if ((t-particle % date_orig)<86400.0_rlg*365.0_rlg) THEN ! to avoid remaining offshore for juveniles (still small) in winter...
            saison=1
         endif
      endif
      if (particle % size < 7.25) saison=2 ! pour cas ou juvenile avant aout, pas de size_class correspondante saison 1
   endif

   if (n==2) then ! remain in the south
      saison=1
      if (particle % size < 7.25) saison=2 ! pour cas ou juvenile avant aout, pas de size_class correspondante saison 1
   endif

   if (n==3) then ! remain in the north
      saison=2
      if(mm <= 2  .OR. mm >= 10) then 
         if ((t-particle % date_orig)<86400.0_rlg*365.0_rlg) THEN ! to avoid remaining offshore for juveniles in winter...
            saison=1
         endif
      endif
      if (particle % size < 7.25) saison=2 ! pour cas ou juvenile avant aout, pas de size_class correspondante saison 1
   endif

   Pi = fish(index,icells,jcells,saison)

   !Select a neighbouring cell j at random among the 4 neighbouring cells
   call random_number(v1) 
   !call random_number(v2)
   !icell = icells+(2*NINT(v1)-1) ! -1 or 1
   !jcell = jcells+(2*NINT(v2)-1)

   icell=icells
   jcell=jcells
   cell = NINT(v1*4.0_rsh-0.5_rsh)
   if (cell==0 .or. cell==-1) icell=icells-1
   if (cell==1) jcell=jcells-1
   if (cell==2) icell=icells+1
   if (cell==3 .or. cell==4) jcell=jcells+1

   icell = MIN(MAX(icell,imin),imax) ! check if in boundaries
   jcell = MIN(MAX(jcell,jmin),jmax)

   ! Get the probability at this new cell
   Pj = fish(index,icell,jcell,saison)

   ! Metropolis algorithm to decide on whether to move or not
   call random_number(harvest) 
   if (Pi==0.0_rsh) move=.TRUE. ! adapted
   if (Pi>0.0_rsh) move= ((Pj>=Pi) .or. (harvest < Pj/Pi))

   ! Calculate displacement (dt forced at 1h, see ibm.f90)
   if (move) then

   ! To modulate movement with fish velocity, but not really compatible with Metropolis...
 !  if (Pi==0.0_rsh .and. h0_g(icells,jcells)>1000.0_rsh .and. particle % size>8.0_rsh) then ! to have the vagrants above a given size (8) offshore back on shelf
 !     xpos_n=particle % xpos+2.0_rsh*particle % size*0.01_rsh*3600.0_rsh/dx(icells,jcells)
 !     ypos_n=particle % ypos
 !  else if (jcell==jcells) then
 !     xpos_n=particle % xpos+2.0_rsh*particle % size*0.01_rsh*3600.0_rsh/dx(icells,jcells)*(icell-icells)
 !     ypos_n=particle % ypos
 !  else if (icell==icells) then
 !     ypos_n=particle % ypos+2.0_rsh*particle % size*0.01_rsh*3600.0_rsh/dy(icells,jcells)*(jcell-jcells)
 !     xpos_n=particle % xpos
 !  else
 !     xpos_n=particle % xpos+2.0_rsh*particle % size*0.01_rsh*3600.0_rsh/dx(icells,jcells)*(icell-icells)/sqrt(2.0_rsh)
 !     ypos_n=particle % ypos+2.0_rsh*particle % size*0.01_rsh*3600.0_rsh/dy(icells,jcells)*(jcell-jcells)/sqrt(2.0_rsh)
 !  endif
 !  if (h0_g(NINT(xpos_n),NINT(ypos_n))>0.0_rsh) then ! test si on reste en mer
 !     if (NINT(xpos_n) .gt. (NINT(particle % xpos)+1) .or. NINT(xpos_n) .lt. (NINT(particle % xpos)-1) .or. NINT(ypos_n) .gt. (NINT(particle % ypos)+1) .or. NINT(ypos_n) .lt. (NINT(particle % ypos)-1)) then
 !        write(*,*)  particle % xpos,xpos_n,particle % ypos,ypos_n
 !        write(*,*) 'saut de CPU'
 !        stop
 !     endif
 !     particle % xpos=xpos_n
 !     particle % ypos=ypos_n
 !  endif

!      if (Pi==0.0_rsh .and. h0_g(icells,jcells)>1000.0_rsh .and. particle % size>8.0_rsh) then ! to have the vagrants above a given size (8) offshore back on shelf
      if (Pi==0.0_rsh) then ! to have the vagrants above a given size offshore back on shelf
         xpos_n=particle % xpos+1
         ypos_n=particle % ypos
      else
         call random_number(harvest) 
         xpos_n=min(real(icell,kind=rsh)+harvest-0.49_rsh, real(icell,kind=rsh)+0.49_rsh) ! to have it randomly within new cell
         call random_number(harvest) 
         ypos_n=min(real(jcell,kind=rsh)+harvest-0.49_rsh, real(jcell,kind=rsh)+0.49_rsh)
      endif

   ! detection d un franchissement limite de domaine
      IF (ypos_n <= 1.0_rsh.or.ypos_n >= real(jmax)) THEN
         particle % flag=-valmanq ! on prefere flagger la particule et la laisser ou elle est (ne bougera plus) 
         particle % super=0
      END IF
      IF (xpos_n <= 1.0_rsh.or.xpos_n >= real(imax)) THEN
         particle % flag=-valmanq ! on prefere flagger la particule et la laisser ou elle est (ne bougera plus) 
         particle % super=0
      END IF

      if (h0_g(NINT(xpos_n),NINT(ypos_n))>0.0_rsh .and. particle % flag /= -valmanq) then
         ! test si on reste en mer (possibilite si incompatibilite de grille)
         if ( (NINT(xpos_n) > icells+1) .or. (NINT(xpos_n) < icells-1) .or. &
              (NINT(ypos_n) < jcells-1) .or. (NINT(ypos_n) < jcells-1) ) then
            write(*,*) 'stop jump MPI'
            stop
         endif
         particle % xpos=xpos_n
         particle % ypos=ypos_n
      endif

   endif    ! move

#ifdef key_MPI_2D

 ! check if we need to exchange or not between CPUs
 ! -----------------------------------
   ipos=NINT(particle % xpos)
   jpos=NINT(particle % ypos)

 ! here we verify if the particle stays in your zone
 ! although it gets out of its zone; it will only move max to 1+-1 and j+-1 
 ! thus we can compute the particle s new value in this time step without pb
 !
 ! if the particle went out the zone; we count on which side it went out and
 ! note it in the down give up give r give l give 
 ! because this value is needed when all the particle finished calculated for this time step
 ! to call ex_traj
   IF ( jpos < ljmin ) THEN
      particle % limitbye = 2
      down_give = down_give+1
      IF ( ipos < limin ) THEN 
         particle % limitbye = 1
      ENDIF
      IF ( ipos > limax ) THEN
         particle % limitbye = 3
      ENDIF
   ELSE IF ( jpos > ljmax ) THEN
      particle % limitbye = 6
      up_give = up_give+1
      IF ( ipos < limin ) THEN
         particle % limitbye = 7
      ENDIF
      IF ( ipos > limax ) THEN 
         particle % limitbye = 5
      ENDIF
   ELSE IF ( ipos < limin ) THEN
      left_give = left_give+1
      particle % limitbye = 8
   ELSE IF ( ipos > limax ) THEN
      particle % limitbye = 4
      right_give = right_give+1
   END IF

#endif

   PRINT_DBG*, 'exit FISH_MOVE routine'

   END SUBROUTINE fish_move

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

#endif
    ! key_ibm

END MODULE
