MODULE ibm
#if defined key_ibm
   !!======================================================================
   !!                   ***  MODULE ibm   ***
   !! 
   !!         M. Huret 01/2011
   !!         IBM (Individual Based model)
   !!         To be adapted for each case for which a CPP key has to be defined
   !&E
   !&E ** History :
   !&E       !  2011-10 (M. Huret) Introduction of MPI
   !&E       !  2012-02 (M. Huret) Introduction of NAMELIST
   !&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,    &
                           rg_valmanq_io,dg_valmanq_io,valmanq
   USE_MPI toolmpi, ONLY : iproc, nprocs

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC :: ibm_3d         ! routine called by step3dxy.F90 step3dyx.F90
   PUBLIC :: ibm_init       ! routine called by main.F90

   !! * Shared module variables

   !! * Private variables
   REAL(kind=rsh)   :: w_max, alpha_w
   INTEGER          :: max_part

CONTAINS

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

   SUBROUTINE ibm_init(xe,sal,temp,          &
#ifdef key_substance
                                   cv_wat,   &
#endif
                                   boz)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ibm_init  ***
   !&E
   !&E ** Purpose : Initialize properties of particles with call to traj_init3d
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : 
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2010-11 (M. Huret) Original code
   !&E       !  2011-10 (M. Huret) Introduction of MPI
   !&E       !  2011-11 (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,    ONLY : iscreenlog
   USE trajinitsave, ONLY : traj_init3d
#if defined key_siggen
   USE trajectools,  ONLY : ztosiggen
#endif
#ifdef key_ibm_species
   USE debmodel,     ONLY : deb_init,debuse,F_Fix,ffix,file_food,file_NBSS
!   USE ibmmove,      ONLY : fish_move_init
   USE ibmtools,     ONLY : ibm_loc_xyz,gasdev_s,ibm_traint,w_dens
#endif
   USE comtraj,      ONLY : type_particle, type_patch, patches,                 &
                            file_trajec, itypetraj, ndtz

   USE ionc4,        ONLY : ionc4_openr, ionc4_read_trajt, ionc4_read_traj,     &
                            ionc4_close, ionc4_read_time
   USE_MPI toolmpi,  ONLY : MASTER

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(     liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: sal,temp,boz
#ifdef key_substance
   REAL(KIND=rsh), DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in ) :: cv_wat
#endif

   !! * Local declarations
   CHARACTER(LEN=lchain)    :: file_inp
   LOGICAL                  :: ibm_restart, ibm_l_time
   INTEGER                  :: i, n, m, nb_part_nc, num
   INTEGER                  :: i_stage
   REAL(KIND=rsh)           :: i_w,i_size,i_density,i_Drate,i_number
   REAL(KIND=rlg)           :: dt_spawn
#ifdef key_ibm_species
   INTEGER                  :: igg,idd,jhh,jbb,kp,km
   INTEGER                  :: hlb,hlt,hrb,hrt  
   REAL(KIND=rsh)           :: px,py
   REAL(KIND=rsh)           :: sal_surf,temp_surf,dens_surf,prof,sprof
   REAL(KIND=rsh)           :: tir,ec_type,ec_type_size
   TYPE(type_particle), POINTER     :: particle
#endif
   TYPE(type_patch),    POINTER     :: patch
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)    :: flag_nc, temp_nc, super_nc
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)    :: dens_nc, size_nc, drate_nc
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)    :: dayb_nc
   INTEGER,        ALLOCATABLE, DIMENSION(:)    :: stage_nc

   NAMELIST/namibmin/      file_trajec, itypetraj
   NAMELIST/namibmphys/    ndtz
   NAMELIST/namibmrestart/ ibm_restart, ibm_l_time
   NAMELIST/namibminit/    i_stage, i_number, i_w, i_size, i_density, i_Drate
   NAMELIST/namibmbio/     w_max, alpha_w
   NAMELIST/namibmout/     dt_spawn, max_part
#ifdef key_ibm_species
   NAMELIST/namibmdeb/     debuse, F_Fix, ffix, file_NBSS, file_food
#endif

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

   PRINT_DBG*, 'enter IBM_init routine'

   ! namelists in paraibm.txt
   !--------------------------

   OPEN(50,file='paraibm.txt',status='old',form='formatted',access='sequential')
   READ(50,namibmin)
   READ(50,namibmphys)
   READ(50,namibmrestart)
   READ(50,namibminit)
   READ(50,namibmbio)
   READ(50,namibmout)
#ifdef key_ibm_species
   READ(50,namibmdeb)
#endif
   CLOSE(50)

   CALL traj_init3d(xe)

   ! save into simu.log
   !-------------------
   IF_MPI (MASTER) THEN
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) '**************************************************'
      WRITE(iscreenlog,*) '***************** IBM_INIT.F90 ****************'
      WRITE(iscreenlog,*) '**************************************************'
      WRITE(iscreenlog,*) ' '
      WRITE(iscreenlog,*) 'fichier definissant les caracteristiques des trajectoires : ',trim(file_trajec)
   ENDIF_MPI      

   patch => patches % first
   DO n = 1,patches % nb

      patch % t_spawn = patch % t_beg
      patch % dt_spawn = dt_spawn*3600.0_rlg

      ! Set patch base ouput name
      i = INDEX(patch % file_out, '.nc')-1
      IF ( i == -1 ) i = LEN_TRIM(patch % file_out)
      patch % file_out_base = patch % file_out(1:i)

      patch % init_particle % stage   = i_stage            
      patch % init_particle % super   = i_number
      patch % init_particle % w       = i_w
      patch % init_particle % size    = i_size  
      patch % init_particle % density = i_density
      patch % init_particle % temp    = 0.0_rsh   
      patch % init_particle % Drate   = i_Drate

#ifdef key_ibm_species
      !Huret et al. 2016
      IF (patch % species == 'anchovy') THEN
               ! Egg density variability ----------------
               ec_type = 0.714_rsh
               ! Egg size and variability ----------------------
               patch % init_particle % size = 0.855_rsh * 0.1_rsh ! egg diameter (mm ->cm)    
               ec_type_size=0.036_rsh * 0.1_rsh
            ENDIF

            IF (patch % species == 'sardine') THEN
               ! Egg density variability ----------------
               ec_type = 0.832_rsh
               ! Egg size and variability ----------------------
               patch % init_particle % size = 1.608_rsh * 0.1_rsh ! egg diameter (mm ->cm) 
               ec_type_size=0.095_rsh * 0.1_rsh
            ENDIF
#endif

      IF ( ibm_restart ) THEN

         file_inp = trim(patch % file_inp)

         nb_part_nc = patch % nb_part_total
         ALLOCATE( flag_nc(nb_part_nc), temp_nc(nb_part_nc),  size_nc(nb_part_nc), stage_nc(nb_part_nc) )
         ALLOCATE( dens_nc(nb_part_nc), super_nc(nb_part_nc), drate_nc(nb_part_nc), dayb_nc(nb_part_nc) )

         CALL ionc4_openr(trim(file_inp), .false.)
         CALL ionc4_read_trajt(file_inp, "flag",      flag_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "TEMP",      temp_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "SIZE",      size_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "DENSITY",   dens_nc,  1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "STAGE",     stage_nc, 1, nb_part_nc, 1)
         CALL ionc4_read_trajt(file_inp, "NUMBER",    super_nc, 1, nb_part_nc, 1)
         CALL ionc4_read_traj (file_inp, "DRATE",     drate_nc, 1, nb_part_nc)
         CALL ionc4_read_traj (file_inp, "DAYBIRTH",  dayb_nc,  1, nb_part_nc)

         DO m = 1,patch % nb_part_alloc
            IF ( .NOT. patch % particles(m) % active ) CYCLE
            num = patch % particles(m) % num
            patch % particles(m) % flag      = flag_nc(num)
            patch % particles(m) % temp      = temp_nc(num)
            patch % particles(m) % size      = size_nc(num)
            patch % particles(m) % density   = dens_nc(num)
            patch % particles(m) % stage     = stage_nc(num)
            patch % particles(m) % super     = super_nc(num)
            patch % particles(m) % Drate     = drate_nc(num)
            patch % particles(m) % date_orig = dayb_nc(num)
         END DO
              
#ifdef key_ibm_species
         CALL ionc4_read_traj(trim(file_inp), "DAYJUV",   dayb_nc, 1, nb_part_nc)
         CALL ionc4_read_traj(trim(file_inp), "DENSPAWN", dens_nc, 1, nb_part_nc)
              
         DO m=1,patch % nb_part_alloc
            IF ( .NOT. patch % particles(m) % active ) CYCLE
            
            num = patch % particles(m) % num            

            patch % particles(m) % dayjuv   = dayb_nc(num)
            patch % particles(m) % denspawn = dens_nc(num)

            IF ( patch % particles(m) % stage > 4 )     \
                 patch % particles(m) % itypevert = 0
         END DO
#endif
         DEALLOCATE( flag_nc, temp_nc, size_nc, stage_nc, dens_nc, super_nc, drate_nc, dayb_nc )

         ! update the date of restart, and savetraj is delayed not to have twice same time step in output
         IF ( ibm_l_time ) THEN
            CALL ionc4_read_time(trim(file_inp), 1, patch % t_beg)
            patch % t_save = patch % t_beg + patch % dt_save*3600.0_rlg
         END IF
         CALL ionc4_close(file_inp)

      ELSE     ! if ( .NOT. ibm_restart ) THEN

         patch % particles(:) % stage   = i_stage            
         patch % particles(:) % super   = i_number
         patch % particles(:) % w       = i_w
         patch % particles(:) % size    = i_size * 0.1_rsh ! egg diameter (mm ->cm)    
         patch % particles(:) % density = i_density
         patch % particles(:) % temp    = 0.0_rsh   
         patch % particles(:) % Drate   = i_Drate

#ifdef key_ibm_species
         !Huret et al. 2016 (egg density =f(surf.density))
         DO m = 1,patch % nb_part_alloc
            
            particle => patch % particles(m)
            IF ( .NOT. particle % active )  CYCLE
                 
            ! Egg density ----------------
                 
            prof = min(4.0_rsh, particle % d3) ! surface
#ifdef key_siggen
            prof = -prof + particle % xe ! switch from immersion to real z
            CALL ztosiggen(prof, sprof, particle % xe, particle % h0, particle % hc)
#else
            sprof = -prof / particle % d3
#endif
            CALL ibm_loc_xyz(xe, particle % xpos, particle % ypos, sprof, &
                             px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt,kp,km)

            ! Egg density = f(surf. density)
            sal_surf  = ibm_traint(sal, xe,sprof,kp,km,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            temp_surf = ibm_traint(temp,xe,sprof,kp,km,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            dens_surf = w_dens(temp_surf, sal_surf)

            IF (patch % species == "anchovy") THEN
               particle % denspawn = -3.434_rsh + 1.1_rsh*(dens_surf-1000.0_rsh)+1000.0_rsh ! with thermal expansion correction
               !particle % denspawn = particle % denspawn + 0.5_rsh ! cf Paul G.
            ENDIF
            IF (patch % species == "sardine") THEN
               particle % denspawn = 5.343_rsh + 0.748_rsh*(dens_surf-1000.0_rsh)+1000.0_rsh ! sardine
               !particle % denspawn = particle % denspawn + 0.78_rh ! sardine
            ENDIF
                 
            CALL gasdev_s(tir)
            particle % denspawn = particle % denspawn + tir*ec_type 
            particle % density  = particle % denspawn
            particle % size = particle % size + tir*ec_type_size
         ENDDO      ! loop on patch % nb_part_alloc
#endif

      ENDIF  ! restart or not

      patch => patch % next

   ENDDO         ! loop on patches % nb

#ifdef key_ibm_species
   IF ( debuse ) CALL deb_init(ibm_restart)
!   CALL fish_move_init
#endif

   ! save initialization
   CALL ibm_save

   PRINT_DBG*, 'leave IBM_init routine'

   END SUBROUTINE ibm_init

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

   SUBROUTINE ibm_3d(xe,uz,vz,sal,temp,    &
#ifdef key_substance
                                     cv_wat,   &
#endif
                                     boz)
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ibm_init  ***
   !&E
   !&E ** Purpose : Main routine with call to advection, behaviour, growth...
   !&E **  To be adapted by the user depending on case study
   !&E
   !&E ** Description :
   !&E ** Called by : traject3d
   !&E ** External calls :
   !&E ** Reference :
   !&E ** History :
   !&E       !  2008-12 (Marc Sourisseau)
   !&E       !  2010-11 (M. Huret)
   !&E       !  2011-10 (M. Huret) Introduction of MPI
   !&E       !  2011-11 (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d,    ONLY : t, dt
#ifdef key_substance
   USE parameters,   ONLY : nb_var
#endif
   USE traject3d,    ONLY : traj_3d
   USE trajinitsave, ONLY : traj_init3d, init_patch, indices_loc2glob
   USE ibmtools,     ONLY : ibm_loc_xyz, ibm_buoy, ibm_traint, ibm_nycth_mig,   &
                            ibm_proftraint, ibm_profuint, ibm_profvint
#if defined key_siggen
   USE trajectools,  ONLY : ztosiggen
#endif
#ifdef key_ibm_species
   USE comvars2d,    ONLY : jjulien
   USE debmodel,     ONLY : deb_init_particle, deb_cycle, debuse
   USE debmodel,     ONLY : F_Fix, readfood3d, Hp, Hj, readtemp3d, climatemp
!   USE ibmmove,      ONLY : fish_move
   USE sflxatm,      ONLY : rad
   USE_MPI toolmpi,  ONLY : ex_traj, ADD_ALL_MPI_INT
#endif
   USE comtraj,      ONLY : type_particle, type_patch, patches, patch_list_append, resize_patch
   USE_MPI comtraj,  ONLY : down_give, up_give, right_give, left_give
   USE_MPI mpi
   USE_MPI toolmpi, ONLY : MPI_COMM_MARS

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),      INTENT( in )    :: xe
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )    :: sal, temp, boz
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )    :: uz, vz
#ifdef key_substance
   REAL(KIND=rsh), DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in ) :: cv_wat
#endif

   !! * Local declarations
   INTEGER                    :: jj,mm,aaaa,hh,minu,sec,n,m,nn
#ifdef key_ibm_species
   LOGICAL                    :: update, night
   REAL(KIND=rsh)             :: depth_day, depth_night
   REAL(KIND=rsh)             :: mindepth, maxdepth
   REAL(KIND=rsh)             :: tempo, sal_part
#endif
   REAL(KIND=rsh)             :: zlag
#ifdef key_siggen
   REAL(KIND=rsh)             :: slag
#endif
   REAL(KIND=rlg)             :: dtm
   INTEGER                    :: kp, km
   INTEGER                    :: hlb, hlt, hrb, hrt  
   INTEGER                    :: igg, idd, jhh, jbb  
   INTEGER                    :: new_eggs, new_eggs_total, ierr_mpi
   INTEGER                    :: idx_s, idx_e, new_size, part_num
   LOGICAL                    :: time_to_spawn, spawn, first_spawn
   CHARACTER(len=19)          :: tool_sectodat
   CHARACTER(len=8)           :: fileout_suffix
   REAL(KIND=rsh)             :: px, py
   TYPE(type_patch),    POINTER     :: patch,    child_patch
   TYPE(type_particle), POINTER     :: particle, new_particle

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

   PRINT_DBG*, 'enter IBM_3d routine'

   ! Save particle properties (before any change for getting exact initial properties) 
   CALL ibm_save

   ! If advection, depend de itypevert passive transport
   CALL traj_3d(xe,uz,vz)

   CALL tool_decompdate(tool_sectodat(t),jj,mm,aaaa,hh,minu,sec)

#ifdef key_ibm_species
   IF ( debuse .AND. .NOT.F_Fix ) THEN
      CALL readfood3d(jjulien)
      CALL readtemp3d(jjulien)
   ENDIF
#endif

 !--------------------------------------------------------------
 ! Common part to all species
 !---------------------------------------------------------------
   dtm = dt*0.5_rlg

#ifdef key_MPI_2D
   down_give  = 0
   up_give    = 0
   right_give = 0
   left_give  = 0
#endif

   patch => patches % first
   DO n = 1,patches % nb

      IF ( (t < patch % t_beg).OR.(t > patch % t_end) ) THEN
         patch => patch % next
         CYCLE
      ENDIF

      new_eggs = 0  ! Number of spawning particles

      DO m = 1,patch % nb_part_alloc

         particle => patch % particles(m)

         IF ( .NOT. particle % active )  CYCLE
         IF ( (particle % super <= 0.0_rsh) .OR. (particle % flag == -valmanq) )  CYCLE

         CALL ibm_loc_xyz(xe, particle % xpos, particle % ypos, particle % spos,    &
                          px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt,kp,km)

       !---------------------------------------------------------------
       ! Specific models with calls to subroutines
       !---------------------------------------------------------------

! Species eggs and larvae
#ifdef key_ibm_species
! size in cm

         update = .TRUE.

        !---------------------------------------------------------------

         IF ( particle % stage == 1 ) THEN  ! stade 1, eggs, buoyancy  
            particle % temp = ibm_traint(temp, xe, particle % spos, kp, km, &
                                         px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
 
            IF (patch % species == 'anchovy') THEN
               ! Egg density dependent on development (Adapted from Ospina 2012 with our data)
               particle % density = particle % denspawn + 3.2029_rsh * particle % Drate - 7.4937_rsh * (particle % Drate)**2 &
                    - 3.0858_rsh * (particle % Drate)**3 + 9.0376_rsh * (particle % Drate)**4 

               ! Egg development rate (Regner, 1996)
               ! temp climato for Drate
               !            tempo = climatemp(kp, NINT(particle % xpos), NINT(particle % ypos))
               !            IF (tempo > 0.0_rsh .and. tempo < 30.0_rsh)  particle % temp = tempo
               particle % Drate = particle % Drate + dtm/86400.0_rsh * 0.000559153_rsh*(particle % temp)**2.290236_rsh
               ! Egg mortality
               !          particle % super = (1.0_rsh-0.266_rsh*dtm/86400.0_rsh)*particle % super
            ENDIF

            IF (patch % species == 'sardine') THEN
               ! Egg density dependent on development (Garcia-Garcia)
               particle % density = particle % denspawn + 3.58_rsh * (0.0012_rsh - 1.8528_rsh * particle % Drate &
                     + 15.507_rsh *(particle % Drate)**2 - 41.1312_rsh * (particle % Drate)**3 + 40.5289_rsh * (particle % Drate)**4) 

               ! Egg developement miranda Peck (dvpt prédit pour dt en jour) dtm en s
               particle % Drate = particle % Drate + (dtm/86400.0_rsh * 0.17_rsh*(particle % temp)**1.9286_rsh) / 100.0_rsh
            ENDIF

            sal_part = ibm_traint(sal, xe, particle % spos,kp,km,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            particle % w = ibm_buoy(particle % density, particle % size, particle % temp, sal_part)

            IF ( particle % Drate >= 1.0_rsh ) THEN
               particle % Drate = 0.0_rsh
               particle % stage = 2
               particle % size  = 0.35_rsh ! on force eclosion taille fixe
               particle % w     = 0.0_rsh
               update = .FALSE.
            ENDIF
         ENDIF ! end of stage 1

         !---------------------------------------------------------------

         IF ( particle % stage == 2 ) THEN  ! stade 2, yolk sac larva, buoyancy  
            ! temp climato
!            tempo = climatemp(kp,NINT(particle % xpos),NINT(particle % ypos))
!            IF (tempo > 0.0_rsh .AND. tempo < 30.0_rsh) THEN
!               particle % temp = tempo
!            ELSE
            particle % temp = ibm_traint(temp, xe, particle % spos, kp, km, &
                                            px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            !ENDIF
            
            ! Buoyancy
            ! sal_part = ibm_traint(sal, xe, particle % spos, kp, km, &
!                                  px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
            ! particle % w = ibm_buoy(particle % density, particle % size, particle % temp, sal_part)

            IF (patch % species == 'anchovy') THEN
               ! Yolk sac larva development rate (Regner, 1996)
               !   particle % Drate = particle % Drate + dtm/86400.0_rsh*0.000003702809_rsh * (particle % temp)**3.8079_rsh
              ! Yolk sac larva development rate (Boussouar et al., 2001), Emmanuelle
               particle % Drate = particle % Drate + dtm/3600_rsh * 0.0000311717 * (particle % temp)**2.1749
            ENDIF
            IF (patch % species == 'sardine') THEN
               ! Yolk sac larvae developement rate Zweifel and Lasker 1976 + Rose et al 2015 for sardinops sagax
               particle % Drate = particle % Drate + dtm/86400*0.001_rsh*exp(6.19_rsh*(1.0_rsh-exp(-0.050_rsh*particle%temp)))/log(5.15_rsh/(5.15_rsh-log(5.970_rsh/3.740_rsh)))     
            ENDIF

            ! Growth (simple as function of Drate as not handled in DEB model so far)
            particle % size = 0.35_rsh + (0.4_rsh-0.35_rsh)*particle % Drate
            IF ( particle % Drate >= 1.0_rsh ) THEN
               particle % stage = 3
               particle % size  = 0.4_rsh ! on force ouverture bouche a taille fixe (cm) tant que yolk sac pas dans DEB, Regner 1996
               update = .FALSE.
            ENDIF
         ENDIF ! end of stage 2

         !---------------------------------------------------------------

         IF ( particle % stage == 3 .AND. update ) THEN ! stage 3, larva without vertical migration
           ! temp climato
!            tempo = climatemp(kp, NINT(particle % xpos), NINT(particle % ypos))
!            IF ( tempo > 0.0_rsh .and. tempo < 30.0_rsh ) THEN
!               particle % temp = tempo
!            ELSE
            particle % temp = ibm_traint(temp, xe, particle % spos, kp, km, &
                                            px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
!            ENDIF
            depth_day   = 25.0_rsh
            depth_night = 25.0_rsh
 
!            sal_part = ibm_traint(sal, xe, particle % spos, kp, km,    &
!                                  px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
!            particle % w = ibm_buoy(particle % density, particle % size, particle % temp, sal_part)

            w_max = particle % size / 10.0_rsh / 100.0_rsh ! m.s-1
            particle % w = ibm_nycth_mig(depth_day, depth_night, w_max, w_max, particle, alpha_w)
                ! to avoid deep diving at slope, to remove later on   

            ! Growth
            IF ( debuse ) CALL deb_cycle(particle,dtm,aaaa,mm,jj)
            IF ( particle % size >= 0.6_rsh ) THEN 
               particle % stage = 4
               update = .FALSE.
            ENDIF
         ENDIF ! end of stage 3
   
         !---------------------------------------------------------------

         IF ( particle % stage == 4 .AND. update ) THEN ! stage 4, nycthemeral migration
            ! temp climato
!            tempo = climatemp(kp, NINT(particle % xpos), NINT(particle % ypos))
!            IF (tempo > 0.0_rsh .and. tempo < 30.0_rsh) THEN
!               particle % temp = tempo
!            ELSE
               particle % temp = ibm_traint(temp, xe, particle % spos, kp, km,  &
                                            px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
!            ENDIF
            depth_day   = 25.0_rsh
            depth_night = 0.0_rsh
            ! Calculate an optimal depth ?
            w_max = particle % size / 10.0_rsh / 100.0_rsh ! m.s-1
            particle % w = ibm_nycth_mig(depth_day, depth_night, w_max, w_max, particle, alpha_w)       
            ! Growth
            IF ( debuse ) CALL deb_cycle(particle,dtm,aaaa,mm,jj)
!            IF ( particle % size >= 3.5_rsh ) THEN  ! to be made more continuous
            IF (particle % H >= Hj) THEN
              particle % stage = 5
              update = .FALSE.
              particle % itypevert = 0 ! Stop vertical advection/diffusion from stage 5
              particle % zpos   = 0.0_rsh
              particle % w      = 0.0_rsh
              particle % dayjuv = jjulien
           ENDIF
        ENDIF ! end of stage 4

         !---------------------------------------------------------------

        IF ( particle % stage == 5 .AND. update ) THEN ! stage 5, Horizontal swimming, distribution from file     
           ! Vertical position
           night = .FALSE.
           IF( rad(NINT(particle % xpos), NINT(particle % ypos)) == 0.0_rsh ) THEN
              night=.TRUE.
           ENDIF
           IF ( particle % H < Hp .AND. jjulien >= particle % dayjuv .AND. jjulien < 330 ) THEN
              IF (night) THEN ! Juvenile Nuit
                 mindepth = 0.0_rsh
                 maxdepth = 30.0_rsh
              ELSE! Juvenile day
                 mindepth = min(30.0_rsh, 2.27_rsh*particle % size) ! Boyra, 2013
                 maxdepth = 30.0_rsh + (120.0_rsh-30.0_rsh)*(1.0_rsh-(330.0_rsh-REAL(jjulien,rsh)) &
                                                                    /(330.0_rsh-REAL(particle % dayjuv,rsh)))
              !  maxdepth = min(particle % d3,120.0_rsh)
              ENDIF
           ELSE
              IF (jjulien <= 90 .OR. jjulien >= 330) THEN !winter 
                 mindepth = 40.0_rsh
                 maxdepth = 120.0_rsh
              ELSE
                 IF (night) THEN ! spring/summer
                    ! Adulte Nuit
                    mindepth = 0.0_rsh
                    maxdepth = 40.0_rsh
                 ELSE ! Adult day
                    mindepth = 30.0_rsh 
                    maxdepth = 120.0_rsh
                 ENDIF
              ENDIF
           ENDIF

           ! temp climato
!           tempo = ibm_proftraint(climatemp, mindepth, maxdepth, px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
!           IF (tempo > 0.0_rsh .and. tempo < 30.0_rsh) THEN
!              particle % temp = tempo
!           ELSE
              particle % temp = ibm_proftraint(temp, mindepth, maxdepth, px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
!           ENDIF
           ! Growth
           IF ( debuse ) CALL deb_cycle(particle,dtm,aaaa,mm,jj)
           ! Movement (time step 5H, cruising speed =1.5BL ~ 4km/h pour TL=15cm) 
           ! if dt or swimming speed too big may jump over 1 cell and MPI will crash (depend also on hor. resolution)
           IF ( hh >= particle % hmove+5 ) THEN
!              CALL fish_move(particle)
              particle % hmove = hh ! update of the saved hour
              IF ( particle % hmove >= 19 ) particle % hmove = particle % hmove-24
           ENDIF
        ENDIF  ! end of stage 5 

        !---------------------------------------------------------------

        ! Mortality
!        IF ( particle % stage > 1 ) particle % super = (1.0_rsh-0.28_rsh*dtm/86400.0_rsh)*particle % super
        IF ( particle % Neggs > 0.0_rsh ) new_eggs = new_eggs+1
#endif
! key_ibm_species


    !--------------------------------------------------------------------------------
    ! Vertical movement and check for bottom and surface boundaries 
         
         zlag = -particle % zpos + particle % xe ! switch from immersion to real z
         zlag = zlag + particle % w*dtm         ! update position
         zlag = MIN(MAX(zlag, -particle % h0), particle % xe) ! check boundaries
      
         particle % zpos = -zlag + particle % xe  ! immersion
      
         IF ( particle % stage <= 2 ) particle % zpos = min(particle % zpos, 40.0_rsh) !patch to avoid going too deep...
      
#ifdef key_siggen
         CALL ztosiggen(zlag, slag, particle % xe, particle % h0, particle % hc)
         particle % spos = slag
#else
         particle % spos = (zlag - particle % xe) / particle % d3                    
#endif

      ENDDO ! particle

#ifdef key_ibm_species
    !--------------------------------------------------------------------------------
    ! Maybe some adults have spawned...

      time_to_spawn = (t >= patch % t_spawn) .AND. (t <= patch % t_end)

      IF ( time_to_spawn ) THEN
         patch % t_spawn = patch % t_spawn + patch % dt_spawn

         ! Gather eggs production information from all MPI domains
         new_eggs_total = new_eggs
#ifdef key_MPI_2D
         CALL ADD_ALL_MPI_INT(new_eggs_total)
#endif
      END IF

      spawn = time_to_spawn .AND. (new_eggs_total > 0)          ! Some particles in the patch are spawning...
      first_spawn = spawn .AND. (.NOT. patch % has_spawned)     ! ... and this is the first time.

      IF ( first_spawn ) THEN

         patch % has_spawned = .TRUE.

         ! If spawning happen for the first time, we have to create an appropriate data structure :
         ! it is a new patch (a child one).
         child_patch => patch_list_append(patches)
         child_patch % t_beg   = t
         child_patch % t_end   = patch % t_end
         child_patch % t_save  = t
         child_patch % dt_save = patch % dt_save
         child_patch % t_spawn = t
         child_patch % dt_spawn = patch % dt_spawn

         child_patch % parent_id   = patch % id                 ! keep track of the childs parent
         child_patch % generation  = patch % generation + 1
         
         ! Inherit default values for new particles
         child_patch % init_particle = patch % init_particle

         ! Set patch output file name
         WRITE( fileout_suffix, '("_G",i0,".nc")' ) child_patch % generation
         child_patch % file_out = trim(patch % file_out_base) // fileout_suffix
         child_patch % file_out_base = patch % file_out_base

         ! Allocate memory for new particles
         CALL init_patch(child_patch, 0)    ! No allocation yet

         child_patch % nb_part_max = max_part

      END IF

      IF ( spawn ) THEN
      
         ! Look for the child patch
         child_patch => patches % first
         DO WHILE ( child_patch % parent_id /= patch % id )
             child_patch => child_patch % next
             IF ( .NOT. ASSOCIATED(child_patch) ) THEN
                ! This shall never happen !
                PRINT*, " ERROR: child patch not found !!! I am ", patch % parent_id, \
                        " looking for", patch % parent_id
                STOP
             END IF
         ENDDO

         ! Compute start index for numbering new particles
         CALL indices_loc2glob( child_patch % nb_part_total, new_eggs, idx_s, idx_e )

         ! Compute new total number of particles
         child_patch % nb_part_total = child_patch % nb_part_total + new_eggs_total

         ! Check if we have not reached the limit allowed by NetCDF output
         IF ( child_patch % nb_part_total > child_patch % nb_part_max ) THEN
            PRINT*, " ERROR: Not enough space in NetCDF file to store additional particles." //         &
                    " We have", child_patch % nb_part_max, " and need", child_patch % nb_part_total,    &
                    " Please consider increasing nb_part_max."
#ifdef key_MPI_2D
            CALL MPI_FINALIZE(ierr_mpi)
#endif
            STOP
         ENDIF

         ! Increase particle array for child patch
         new_size = child_patch % nb_part_alloc + new_eggs
         CALL resize_patch(child_patch, new_size)

         ! Loop over all spawning particles of parent patch
         m = 0
         part_num = 0
         DO nn = 1,patch % nb_part_alloc

            particle => patch % particles(nn)

            ! Skip particle if inactive or has not spawn
            IF ( .NOT. particle % active )     CYCLE
            IF ( particle % Neggs <= 0.0_rsh ) CYCLE

            ! Initialize a new particle
            new_particle => NULL()
            DO WHILE ( .NOT. ASSOCIATED(new_particle) )
               m = m+1
!               IF ( m > child_patch % nb_part_alloc ) THEN
!                  PRINT*, " ERROR: max number of particles reached: ", child_patch % nb_part_alloc, \
!                          " This should never happen !"    ! <- Yes indeed !
!                  CALL resize_patch(child_patch, 1)
!               END IF
               IF ( .NOT. child_patch % particles(m) % active ) THEN 
                   new_particle => child_patch % particles(m)
                   part_num = part_num+1
               END IF
            END DO

            new_particle = patch % init_particle
            new_particle % active = .TRUE.
            new_particle % num = idx_s + part_num
            new_particle % date_orig = t
            new_particle % xpos = particle % xpos
            new_particle % ypos = particle % ypos
            new_particle % spos = particle % spos
            new_particle % zpos = particle % zpos
            new_particle % d3   = particle % d3
            new_particle % h0   = particle % h0
            new_particle % xe   = particle % xe
            new_particle % hc   = particle % hc
            new_particle % super = particle % Neggs

            ! Reset number of eggs in spawning particle
            particle % Neggs = 0.0_rsh

            ! Biological parameters
            IF ( debuse ) CALL deb_init_particle(new_particle)

         ENDDO
         
      END IF    ! ( spawn )
#endif
! key_ibm_species

      patch => patch % next

   ENDDO ! patch

#ifdef key_ibm_species
   CALL_MPI ex_traj(down_give, up_give, right_give, left_give) ! because fish_move
#endif

   END SUBROUTINE ibm_3d

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

   SUBROUTINE ibm_save
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ibm_save  ***
   !&E
   !&E ** Purpose : Save trajectories
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : ibm_init and ibm_3d (IBM/ibm.F90)
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2011-01 (M. Huret) 
   !&E       !  2011-10 (M. Huret) Introduction of MPI
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,   ONLY : t
   USE ionc4,       ONLY : ionc4_createfile_traj, ionc4_createvar_traj, &
                           ionc4_write_traj, ionc4_write_trajt,         &
                           ionc4_write_time, ionc4_sync
   USE comtraj,     ONLY : patches, type_patch, type_particle

   USE trajinitsave, ONLY : indices_loc2glob
   USE_MPI toolmpi, ONLY : ADD_ALL_MPI_INT

   !! * Arguments

   !! * Local declarations
   CHARACTER(LEN=lchain)    :: file_out
   LOGICAL                                     :: l_out_nc4par
   INTEGER                                     :: n,m,num1,num2,nb_part,nb_part_nc,p
   REAL(kind=rsh)                              :: fillval
   REAL(KIND=rlg)                              :: tool_ind2lat,tool_ind2lon
   INTEGER,        ALLOCATABLE, DIMENSION(:)   :: num_out,stage_out
   REAL(KIND=rlg), ALLOCATABLE, DIMENSION(:)   :: lat_out,lon_out,dateo_out
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: temp_out,flag_out,spos_out,zpos_out,xpos_out,ypos_out 
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: h0pos_out,size_out,nb_out,dens_out,Drate_out
#ifdef key_ibm_species
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: food_out,Wdeb_out,Denspawn_out
!   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:)   :: E_out,H_out,R_out,Neggs_out,NRJ_out,Gam_out
!   INTEGER,        ALLOCATABLE, DIMENSION(:)   :: dayjuv_out, dayspawn_out,yearspawn_out,season_out
#endif
   TYPE(type_patch),    POINTER     :: patch
   TYPE(type_particle), POINTER     :: particle
   INTEGER :: idx_s, idx_e

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

   PRINT_DBG*, 'enter IBM_save routine'

   IF(rsh==8) THEN 
      fillval = dg_valmanq_io
   ELSE
      fillval = rg_valmanq_io
   ENDIF

#ifdef key_MPI_2D
   l_out_nc4par = .true.
#else
   l_out_nc4par = .false.
#endif
   patch => patches % first

   DO n = 1,patches % nb

      IF ( (t < patch % t_save) .OR. (t > patch % t_end) ) THEN
         patch => patch % next
         CYCLE
      END IF

      file_out = trim(patch % file_out)

      IF ( .NOT. patch % file_out_init ) THEN
         ! Create output file
         patch % file_out_init = .TRUE.
         nb_part_nc = patch % nb_part_max

         CALL ionc4_createfile_traj(file_out, nb_part_nc, 0, 0, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "latitude", "degrees_north","latitude",                    &
                                             fill_value=dg_valmanq_io, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "longitude","degrees_east", "longitude",                   &
                                             fill_value=dg_valmanq_io, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DEPTH","m","depth as immersion",                          &
                                             fill_value=-fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "H0","m","h0",                                             &
                                             fill_value=-fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "NUM","","number of the particle",                         &
                                             fill_value=0,        l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "flag","nbr","flag",                                       &
                                             fill_value=fillval,  l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "TEMP","degrees_Celsius","temperature",                    &
                                             fill_value=fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "SIZE","Centimeters","Size",                               &
                                             fill_value=fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "STAGE","","Stage",                                        &
                                             fill_value=-1,       l_out_nc4par=l_out_nc4par)
  !       CALL ionc4_createvar_traj(file_out, "NUMBER","Number","Number of individuals",                 &
  !                                           fill_value=fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DENSITY","sigma","Density",                               &
                                             fill_value=fillval, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DRATE","","Development rate of egg or larva",             &
                                             fill_value=fillval,  l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DAYBIRTH","","Date of birth",                             &
                                             fill_value=dg_valmanq_io, l_out_nc4par=l_out_nc4par, ndims=1)
#ifdef key_ibm_species
         CALL ionc4_createvar_traj(file_out, "FOOD","mg/m3","food",                                     &
                                             fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "EDEB","Joules","Energy",                                  &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "HDEB","Joules","Maturity",                                &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "RDEB","Joules","Repro",                                   &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "GAM","Joules","Energy gametes",                           &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "WEIGHT","g","weight",                                     &
                                             fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "NEGGS","Number","Number of eggs",                         &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "NRJ","J/g","Energy Density",                              &
    !                                         fill_value=fillval,  l_out_nc4par=l_out_nc4par)
     
    !     CALL ionc4_createvar_traj(file_out, "YEARSPAWN","","Year authorised to spawn",                 &
    !                                         fill_value=0, ndims=1, l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "DAYSPAWN","","Julian day start spawning",                 &
    !                                         fill_value=0, ndims=1, l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "DAYJUV","","Julien day at metamorphosis",                 &
    !                                         fill_value=0, ndims=1, l_out_nc4par=l_out_nc4par)
    !     CALL ionc4_createvar_traj(file_out, "SEASON","","Wether within spawning season",               &
    !                                         fill_value=-1, ndims=1, l_out_nc4par=l_out_nc4par)
         CALL ionc4_createvar_traj(file_out, "DENSPAWN","sigma","Density of egg at spawning",           &
                                             fill_value=fillval, ndims=1, l_out_nc4par=l_out_nc4par)
#endif
      END IF    ! (.NOT. patch % file_out_init)

      nb_part = patch % nb_part_alloc
      ALLOCATE( lat_out(nb_part),  lon_out(nb_part) )
      ALLOCATE( xpos_out(nb_part), ypos_out(nb_part) )
      ALLOCATE( spos_out(nb_part), zpos_out(nb_part) )
      ALLOCATE( num_out(nb_part),  h0pos_out(nb_part) )
      ALLOCATE( flag_out(nb_part), dens_out(nb_part) )
      ALLOCATE( temp_out(nb_part), Drate_out(nb_part) )
      ALLOCATE( size_out(nb_part), dateo_out(nb_part) )
      ALLOCATE( stage_out(nb_part), nb_out(nb_part) )

      lat_out(:)  = dg_valmanq_io ; lon_out(:) = dg_valmanq_io ; dateo_out(:) = dg_valmanq_io
      xpos_out(:) = fillval ; ypos_out(:) = fillval ; spos_out(:) = fillval ; zpos_out(:) = -fillval
      h0pos_out(:)=-fillval ; flag_out(:) = fillval ; num_out(:)  = 0
      size_out(:) = fillval ; stage_out(:) = -1 ; dens_out(:) = fillval ; Drate_out(:) = fillval
      temp_out(:) = fillval ; nb_out(:) = fillval

#ifdef key_ibm_species
   !   ALLOCATE( dayjuv_out(nb_part), dayspawn_out(nb_part) )
   !   ALLOCATE( yearspawn_out(nb_part), season_out(nb_part) )
      ALLOCATE( Denspawn_out(nb_part) )
      ALLOCATE( food_out(nb_part))
   !   ALLOCATE( H_out(nb_part), E_out(nb_part), Gam_out(nb_part), R_out(nb_part))
      ALLOCATE(Wdeb_out(nb_part) )
   !   ALLOCATE( Neggs_out(nb_part), NRJ_out(nb_part) )

   !   dayjuv_out(:) = 0 ; dayspawn_out(:) = 0 ; yearspawn_out(:) = 0 ; season_out(:) = -1
      Denspawn_out(:) = fillval ; food_out(:) = fillval ; Wdeb_out(:) = fillval
   !   H_out(:) = fillval ; E_out(:) = fillval ; R_out(:) = fillval ; Gam_out(:) = fillval
   !   Neggs_out(:) = fillval ; NRJ_out(:) = fillval
#endif

      p = 0
      DO m = 1,nb_part
         particle => patch % particles(m)
         IF ( .NOT. particle % active ) CYCLE
         p = p+1
         xpos_out(p)  = particle % xpos
         ypos_out(p)  = particle % ypos
         lat_out(p)   = tool_ind2lat(particle % xpos, particle % ypos)
         lon_out(p)   = tool_ind2lon(particle % xpos, particle % ypos)
         dateo_out(p) = particle % date_orig
         spos_out(p)  = particle % spos
         zpos_out(p)  = particle % zpos ! for output as immersion
         h0pos_out(p) = particle % d3
         flag_out(p)  = particle % flag
         temp_out(p)  = particle % temp
         size_out(p)  = particle % size
         stage_out(p) = particle % stage
         dens_out(p)  = particle % density
    !     nb_out(p)    = particle % super
         num_out(p)   = particle % num
         Drate_out(p) = particle % Drate
#ifdef key_ibm_species
         food_out(p)  = particle % X
    !     E_out(p)     = particle % E
    !     H_out(p)     = particle % H
    !     R_out(p)     = particle % R
    !     Gam_out(p)   = particle % Gam
         Wdeb_out(p)  = particle % Wdeb
    !     Neggs_out(p) = particle % Neggs
    !     NRJ_out(p)   = particle % NRJ
       
    !     yearspawn_out(p)= particle % yearspawn
    !     dayspawn_out(p) = particle % dayspawn
    !     dayjuv_out(p)   = particle % dayjuv
    !     season_out(p)   = particle % season
         denspawn_out(p) = particle % denspawn
#endif
      END DO

      CALL ionc4_write_time(file_out,0,t)

      nb_part = count(patch%particles(:)%active)
      CALL indices_loc2glob(1, nb_part, idx_s, idx_e)

#ifdef key_MPI_2D
      num1 = idx_s
      num2 = idx_e
#else
      num1 = 1
      num2 = nb_part
#endif
      CALL ionc4_write_trajt(file_out, 'latitude',  lat_out(1:nb_part),num1,num2,0,dg_valmanq_io)
      CALL ionc4_write_trajt(file_out, 'longitude', lon_out(1:nb_part),num1,num2,0,dg_valmanq_io)
      CALL ionc4_write_trajt(file_out, 'DEPTH',    zpos_out(1:nb_part),num1,num2,0,-fillval)
      CALL ionc4_write_trajt(file_out, 'H0',      h0pos_out(1:nb_part),num1,num2,0,-fillval)
      CALL ionc4_write_trajt(file_out, 'NUM',       num_out(1:nb_part),num1,num2,0,0)
      CALL ionc4_write_trajt(file_out, 'flag',     flag_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_trajt(file_out, 'TEMP',     temp_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_trajt(file_out, 'SIZE',     size_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_trajt(file_out, 'STAGE',   stage_out(1:nb_part),num1,num2,0,-1)
      CALL ionc4_write_trajt(file_out, 'DENSITY',  dens_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'NUMBER',     nb_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_trajt(file_out,  'DRATE',   Drate_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_traj(file_out,  'DAYBIRTH',dateo_out(1:nb_part),num1,num2,dg_valmanq_io)
#ifdef key_ibm_species
      CALL ionc4_write_trajt(file_out, 'FOOD',          food_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'EDEB',             E_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'HDEB',             H_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'RDEB',             R_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'GAM',            Gam_out(1:nb_part),num1,num2,0,fillval)
      CALL ionc4_write_trajt(file_out, 'WEIGHT',        Wdeb_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'NEGGS',        Neggs_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_trajt(file_out, 'NRJ',            NRJ_out(1:nb_part),num1,num2,0,fillval)
   !   CALL ionc4_write_traj(file_out,  'YEARSPAWN',yearspawn_out(1:nb_part),num1,num2,0)
   !   CALL ionc4_write_traj(file_out,  'DAYSPAWN',  dayspawn_out(1:nb_part),num1,num2,0)
   !   CALL ionc4_write_traj(file_out,  'DAYJUV',      dayjuv_out(1:nb_part),num1,num2,0)
   !   CALL ionc4_write_traj(file_out,  'SEASON',      season_out(1:nb_part),num1,num2,0)
      CALL ionc4_write_traj(file_out,  'DENSPAWN',  denspawn_out(1:nb_part),num1,num2,fillval)
#endif

      ! To write the data on the disk and not loose data in case of run crash
      CALL ionc4_sync(file_out)
      DEALLOCATE(lat_out, lon_out, xpos_out, ypos_out, spos_out, zpos_out)
      DEALLOCATE(num_out, h0pos_out, flag_out, dens_out, temp_out, Drate_out)
      DEALLOCATE(size_out, dateo_out, stage_out, nb_out)
#ifdef key_ibm_species
   !   DEALLOCATE(dayjuv_out, dayspawn_out, yearspawn_out, season_out)
      DEALLOCATE(Denspawn_out, food_out, Wdeb_out)
   !   DEALLOCATE(, Gam_out, H_out, E_out, R_out, Neggs_out, NRJ_out)
#endif

      patch % t_save = patch % t_save + patch % dt_save*3600.0_rlg
      patch => patch % next
      
   END DO

   PRINT_DBG*, 'leave IBM_save routine'

   END SUBROUTINE ibm_save

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

#endif
END MODULE
