MODULE ibmtools

#if defined key_ibm

   !!======================================================================
   !!                   ***  MODULE ibm   ***
   !! 
   !!         M. Huret 01/2011
   !!         Generic routines and functions for IBM
   !&E
   !&E ** History :
   !&E       !  2013-04 (M. Huret) Separation from IBM
   !&E       !  2014-12 (M. Honnorat) Adapt for IBM upgrade
   !&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
   USE comvars3d,  ONLY : sigw,sig
   USE comtraj,    ONLY : type_particle


   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC :: ibm_buoy, ibm_traint, ibm_nycth_mig, ibm_proftraint, ibm_profuint, &
             ibm_profvint, ibm_loc_xyz, gasdev_s, w_dens, ibm_profmean  ! routine called ibm.F90

   !! * Shared module variables
   REAL(kind=rsh),PUBLIC, dimension(imin:imax,jmin:jmax)  ::  opt_depth_part

   !! * Private variables
   REAL(kind=rsh),dimension(imin:imax,jmin:jmax)  :: depth_strat !-1 no stratif; >0 stratif

#ifdef key_ibm_unused
   REAL(KIND=rsh),  ALLOCATABLE, DIMENSION(:,:)  :: xep
   REAL(KIND=rsh),  ALLOCATABLE, DIMENSION(:,:)  :: dxe
   REAL(KIND=rsh),  ALLOCATABLE, DIMENSION(:,:)  :: dxep
   REAL(KIND=rlg),  ALLOCATABLE, DIMENSION(:,:)  :: deb_descend
   REAL(KIND=rlg),  ALLOCATABLE, DIMENSION(:,:)  :: deb_monte
#endif

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

   FUNCTION w_dens(tempw,salw)

   !&E---------------------------------------------------------------------
   !&E ** Purpose : Estimate the water density
   !&E ** 
   !&E ** Description : 
   !&E ** Called by : w_behav_part
   !&E ** External calls : 
   !&E ** Reference : 
   !&E ** History : ! M. Huret 11/1010
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : ierrorlog,icon
     USE comvars3d,  ONLY : l_stateeq_lin,rhoref

   !! * Arguments
   REAL(KIND=rsh)   , INTENT( in )  :: salw,tempw
   REAL(KIND=rsh)                   :: w_dens

   !! * Local declarations
   INTEGER                 :: IERR_MPI
   REAL(KIND=rsh)          :: TP,S,P
   REAL(KIND=rlg)          :: row,rost,c,pcc
   REAL(kind=rsh)          :: alpha,beta

   !!----------------------------------------------------------------------
   !! * Executable part
 
   IF (icon >= 1 .AND. l_stateeq_lin) THEN

    alpha=0.00016_rsh
    beta=0.0008_rsh
    w_dens = (beta*salw-alpha*tempw)*rhoref

   ELSE IF (icon >= 1 .AND. .NOT.l_stateeq_lin) THEN
! water density (Mellor, 1991, see SUBROUTINE ts_stateeq_blumberg)
     P = 10.0_rsh ! pressure
     TP = MAX(tempw,0.001_rsh)
     S=MAX(salw,0.001_rsh)

   ! density of pure water
     row = 999.842594_rlg + (  6.793952d-02 + ( - 9.095290d-03 +                &
          (1.001685d-04 - 1.120083d-06*TP + 6.536332d-09*TP*TP)*TP )*TP)*TP

   ! density of sea water at the surface
     rost = row + S * (   0.824493_rlg + (  -4.0899d-03 + ( 7.6438d-05 +        &
        (-8.2467d-07+5.3875d-09*TP)*TP )*TP  )*TP   ) +           &
        S**1.5_rlg * ( -5.72466d-03 + (1.0227d-04-1.6546d-06*TP)*TP ) +     &
        S*S*4.8314d-04

   ! density of sea water at pressure p (decibars)
     c = 1449.2_rlg + 1.34_rlg*(S-35.0_rlg) + 4.55_rlg*TP - 0.045_rlg*TP*TP +   &
          0.00821_rlg*p + 15.0d-09*p*p
     pcc = p/(c*c)
     w_dens = rost + 1.0d4*(1.0_rlg-0.2_rlg*pcc)*pcc

!------ Fin MELLOR
   ELSE

     WRITE(ierrorlog,*) 'Function IBM_BUOY'
     WRITE(ierrorlog,*) 'No state equation have been chosen'
     WRITE(ierrorlog,*) 'because icon<1'
     WRITE(ierrorlog,*) 'The simulation is stopped'
     CALL_MPI MPI_FINALIZE(IERR_MPI)
     STOP

   END IF

   END FUNCTION w_dens


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

   FUNCTION ibm_buoy(density,size,tempw,salw)

   !&E---------------------------------------------------------------------
   !&E ** Purpose : Estimate the terminal velocity from buoyancy scheme
   !&E ** 
   !&E ** Description : 
   !&E ** Called by : w_behav_part
   !&E ** External calls : 
   !&E ** Reference : 
   !&E ** History : ! M. Huret 11/1010
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE parameters, ONLY : grav

   !! * Arguments
   REAL(KIND=rsh)   , INTENT( in )  :: salw,tempw,density,size
   REAL(KIND=rsh)                   :: ibm_buoy

   !! * Local declarations
   REAL(KIND=rsh)          :: w_buoy
   REAL(KIND=rsh)          :: wsto,wdal,cor,rayoeuf,diam,difdens
   REAL(KIND=rsh)          :: TP,S,wdens,dynvisc,reyn

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

   cor=1.0_rsh
   rayoeuf=size*0.01_rsh/2.0_rsh    ! diameter to radius

   TP = MAX(tempw,0.001_rsh)
   S=MAX(salw,0.001_rsh)
   
   wdens=w_dens(tempw,salw)

! Dynamic viscosity (kg.m-1.s-1) (approx. of table values, vertegg v1.0 toolbox, Adlandsvik, 2000)
   dynvisc=(1.7915_rsh - 0.0538_rsh*TP + 0.0007_rsh*TP*TP + 0.0023_rsh*S)*0.001_rsh

!  calcul de la vitesse de remontee des oeufs de poissons (Stokes)             
!  unites: m s-1
   difdens=density-REAL(wdens,rsh)
   wsto=2.0_rsh*grav*rayoeuf**2.0_rsh*difdens/(9.0_rsh*dynvisc)
   wsto=wsto*cor
   w_buoy=-wsto
              
! calcul du nombre de Reynolds (Re)
   reyn=(wdens*abs(wsto)*2.0_rsh*rayoeuf)/dynvisc

! choix de la vitesse de Dallavalle si Re > 0.5     
   if (reyn.gt.0.5_rsh) then
      diam=(9.0_rsh*(dynvisc**2)/(wdens*grav*abs(difdens)))**(1.0_rsh/3.0_rsh)      
      wdal=0.0875_rsh*abs(difdens)**(2.0_rsh/3.0_rsh)*((2.0_rsh*rayoeuf)-0.4_rsh*diam)/dynvisc**(1.0_rsh/3.0_rsh)

      w_buoy=wdal*cor
      if(difdens.ge.0.0_rsh) then
         w_buoy=-wdal*cor
      endif
   
   endif

   ibm_buoy=w_buoy

  END FUNCTION ibm_buoy

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

  FUNCTION ibm_nycth_mig(depth1,depth2,w1,w2,p,alpha)

   !&E---------------------------------------------------------------------
   !&E
   !&E ** Purpose : Nycthemeral migration between depth1 (day) and depth2 (night)
   !&E ** Migration from depth1 to depth2 starts when rad=0 with w=w2
   !&E ** Migration from depth2 to depth1 starts when rad>0 with w=w1
   !&E **
   !&E ** Description :
   !&E **
   !&E ** Called by : ibm_3d
   !&E **
   !&E ** External calls :
   !&E ** 
   !&E ** Reference :
   !&E ** 
   !&E ** History : M. Huret 11/2010
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d, ONLY : t
   USE sflxatm,      ONLY : rad

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )          ::depth1,depth2,w1,w2,alpha
!   INTEGER, INTENT( in )                 :: n,m 
   TYPE(type_particle), INTENT(inout)    :: p
   REAL(KIND=rsh)                        :: ibm_nycth_mig

   !! * Local declarations
   REAL(kind=rsh)                        ::   depth,w
   character(len=19)                     ::  tool_sectodat,date
   integer                               :: jj,mm,aaaa,hh,minu,sec

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

   depth=depth1
   w=w1
   if(rad(NINT(p%xpos),NINT(p%ypos)).eq.0.0_rsh) THEN
!   date=tool_sectodat(t)     ! Hour of the day
!   call tool_decompdate(date,jj,mm,aaaa,hh,minu,sec)
!   if(hh <= 6. .OR. hh >= 21.) THEN
      depth=depth2
      w=w2
   endif
   if (depth .lt. p%d3) THEN  
      ! p%zpos is immersion of particle
      ibm_nycth_mig= w*tanh(alpha*(p%zpos-depth)) 
   else
      ibm_nycth_mig=0.
   endif
   
  END FUNCTION ibm_nycth_mig

!!======================================================================
#ifdef key_ibm_unused

  FUNCTION ibm_tidal_mig(xe,depth1,depth2,w,p,alpha)

   !&E---------------------------------------------------------------------
   !&E
   !&E ** Purpose : Tidal migration between depth1 (flood) and depth2 (ebb) 
   !&E ** 
   !&E ** Description :
   !&E ** Called by : w_behav_part
   !&E ** External calls :
   !&E ** 
   !&E ** Reference :
   !&E ** 
   !&E ** History : M. Huret 11/2010 (A ameliorer)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E---------------------------------------------------------------------

   !! * Modules used
   USE comvars2d, ONLY : t

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT( in )   :: xe
   REAL(KIND=rsh), INTENT( in )       :: depth1,depth2,alpha,w
   REAL(KIND=rsh)                              :: ibm_tidal_mig
!   INTEGER, INTENT( in )                 :: n,m
   TYPE(type_particle)                   :: p

   
   !! * Local declarations
   INTEGER                 :: jj,mm,aaaa,hh,minu,sec
   REAL(kind=rsh) ::   depth

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

   dxe(n,m) = xe(NINT(p%xpos),NINT(p%ypos))-xep(n,m) 
   
   ! Si chgmt d elevation decroissant apres au moins 6h
   if(dxep(n,m) > 0.0 .AND. dxe(n,m) <= 0.0 .AND. t >= deb_monte(n,m) + 6.0*3600.0) then
      deb_monte(n,m) = t + 4.0*3600.0 + 15.0*60.0 ! chgt de courant dans 4heures 15
   endif
   
   ! Si chgmt d elevation croissant apres au moins 6h
   if(dxep(n,m) <= 0.0 .AND. dxe(n,m) > 0.0 .AND. t >= deb_descend(n,m) + 6.0*3600.0) then
      deb_descend(n,m) = t + 4.0*3600.0 + 15.0*60.0 ! chgt de courant dans 4heures 15
   endif
   
   depth=depth1 ! en surface qd maree monte
   if(t < deb_monte(n,m) .AND. t >= deb_descend(n,m)) THEN 
      depth=depth2 ! au fond qd maree descend
   endif
   
   dxep(n,m) = dxe(n,m)
   xep(n,m)=xe(NINT(p%xpos),NINT(p%ypos))
   
   ibm_tidal_mig= w*tanh(alpha*(p%zpos-depth)) 
   
  END FUNCTION ibm_tidal_mig

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

  SUBROUTINE ibm_detect_pycno(xe,bz,sal,temp)

   !&E---------------------------------------------------------------------
   !&E ** Purpose : detect pycnocl for part behaviour
   !&E ** Description :
   !&E ** Called by : w_behav_part
   !&E ** External calls :
   !&E ** Reference :
   !&E ** History : Marc Sourisseau 2008 ! non terminée mais pas encore utilisée
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : h0_g,h0fond,hm,ig,id
   USE toolgeom, ONLY : f_dzu

   !! * 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
   REAL(KIND=rsh), DIMENSION(kmax,imin:imax,jmin:jmax), INTENT( in )  :: bz

   !! * Local declarations
   INTEGER  :: k_strat_max,i,j,k
   REAL(kind=rsh) :: d,roref,dzu
   REAL(kind=rsh),dimension(kmax-1,imin:imax,jmin:jmax) ::  gradv_bz
   INTEGER,dimension(imin:imax,jmin:jmax)  :: k_stratif


   !!----------------------------------------------------------------------
   !! * Executable part
!
! calcul du gradient vertical de densite
!-----------------------
!       sig au centre des cellules et sigw position verticale entre des limites (calcul des w)

!cval Make sure you use the same value of roref than in tsstateeqblumberg.F90
     roref=1040.0_rsh

     do  j=ljmin,ljmax
!      do  i=ig(j),id(j)
      do i=MAX(limin,ig(j)+1),MIN(limax,id(j)-1)


       depth_strat(i,j)=0.0_rsh
       d=h0_g(i,j)+xe(i,j)

       if (d .gt. hm) THEN
        k_strat_max=1
        gradv_bz(1,i,j)=0.0_rsh            !grad initial ; grad > 0
        do k=2,kmax

        dzu=REAL(f_dzu(h0_g(i,j),xe(i,j),k,i,j),rsh)   !dsigu=sigw(k)-sigw(k-1) => dzu>0
        gradv_bz(k,i,j)=((roref-bz(k-1,i,j)/grav*roref) - (roref-bz(k,i,j)/grav*roref)) &
                         /dzu !grad between k-1 and k ; grad > 0 en m-1


          if (gradv_bz(k,i,j) .gt. gradv_bz(k_strat_max,i,j)) then
           k_strat_max=k           ! numero de la couche ou le gradv de flottab. (bz) est max.sur sa frontiere inf.
          END IF
        enddo !endo on k

        if ((abs(gradv_bz(k_strat_max,i,j)) .gt. 0.02_rsh) ) then 
                                 ! .and. (d*(1.0_rsh-(sigw(k_strat_max-1))) .gt. 60.0_rsh)
                                 ! 0.02 m-1 uniquement pour les panaches / stratif estivale non visible
                                 ! 0.000005 m-1 stratif obs partout de mars à sept

         k_stratif(i,j)=k_strat_max-1  !indic of sigma levels with the highest upper gradient
         depth_strat(i,j)=d*(1.0_rsh-(sigw(k_stratif(i,j))))   !depth in m of the highest buoanc. gradient
        else
         depth_strat(i,j)=-1.0_rsh
        END IF


       END IF
      enddo
     enddo

   END SUBROUTINE ibm_detect_pycno

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

   SUBROUTINE ibm_opt_depth(xe)

   !&E---------------------------------------------------------------------
   !&E
   !&E ** Purpose : Estimate the optimal depth for part. behaviour.
   !&E ** Description :
   !&E ** Called by : w_behav_part
   !&E ** External calls :
   !&E ** Reference :
   !&E ** History : marc 11/2008 ! pour l instant opt_depth fixe
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : h0_g,h0fond,hm,ig,id

   !! * Arguments
   REAL(KIND=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),INTENT( in )   :: xe

   !! * Local declarations
   INTEGER  :: k_strat_max,i,j,k
   REAL(kind=rsh) :: d
 
   !!----------------------------------------------------------------------
   !! * Executable part
!
! estimation of the optimal vertical depth
!-----------------------
!
     do  j=ljmin,ljmax
      do  i=MAX(limin,ig(j)+1),MIN(limax,id(j)-1)

       d= h0_g(i,j)+xe(i,j)
       if (d .gt. hm) THEN
       !  depth_strat(i,j)=15.0_rsh;

         if (depth_strat(i,j) .gt. 0.0_rsh) then
          opt_depth_part(i,j)=depth_strat(i,j)
         else
          opt_depth_part(i,j)=0.0_rsh ! dpeth in m
         END IF

       END IF
      enddo
     enddo

   END SUBROUTINE ibm_opt_depth

!!======================================================================
   FUNCTION ibm_profuint(var,xpos,ypos,minprof,maxprof)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION ibm_profuint  ***
   !&E
   !&E ** Purpose :Interpolate vertically averaged current 
   !&E          from the 4 neighbooring location 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  04-2013 M. Huret 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : htx,hx_g
     USE trajectools, ONLY : interpvit
#if defined key_siggen 
     USE comsiggen, ONLY : hcx_sig
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: var
   REAL(KIND=rsh), INTENT( in )                :: xpos,ypos
   REAL(KIND=rsh), INTENT( in )                :: minprof, maxprof
   REAL(KIND=rsh)                              :: ibm_profuint

   !! * Local declarations
   REAL(KIND=rsh), DIMENSION(2,2)              :: val4
   INTEGER                    :: igg,idd,jhh,jbb
   INTEGER                    :: hlb,hlt,hrb,hrt  
   REAL(KIND=rsh)             :: xe,px,py


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

  !---- reperage de la maille ou se situe la particule pour les U
   ! location of ssh : xpos=REAL(i), ypos=REAL(j)
   idd=nint(xpos)
   igg=idd-1
   jbb=int(ypos)
   jhh=jbb+1

   idd=min(max(idd,imin),imax)
   igg=min(max(igg,imin),imax)
   jbb=min(max(jbb,jmin),jmax)
   jhh=min(max(jhh,jmin),jmax)

   ! calcul des positions px et py dans la maille U
   px=xpos-real(igg)-0.5_rsh
   py=ypos-real(jbb)


   val4(:,:)=0.0_rsh
   ! Processing prof mean on 4 neighbooring locations
   if (htx(igg,jbb)>0.0_rsh) THEN
      xe=htx(igg,jbb)-hx_g(igg,jbb)
      val4(1,1)=ibm_profmean(var(:,igg,jbb),minprof,maxprof,hx_g(igg,jbb),xe &
#if defined key_siggen 
           ,hcx_sig(igg,jbb) &
#endif
           )
   END IF
   if (htx(idd,jbb)>0.0_rsh) THEN
      xe=htx(idd,jbb)-hx_g(idd,jbb)
      val4(2,1)=ibm_profmean(var(:,idd,jbb),minprof,maxprof,hx_g(idd,jbb),xe &
#if defined key_siggen 
      ,hcx_sig(idd,jbb) &
#endif
           )
   END IF
   if (htx(idd,jhh)>0.0_rsh) THEN
      xe=htx(idd,jhh)-hx_g(idd,jhh)
      val4(2,2)=ibm_profmean(var(:,idd,jhh),minprof,maxprof,hx_g(idd,jhh),xe &
#if defined key_siggen 
      ,hcx_sig(idd,jhh) &
#endif
           )
   END IF
   if (htx(igg,jhh)>0.0_rsh) THEN
      xe=htx(igg,jhh)-hx_g(igg,jhh)
      val4(1,2)=ibm_profmean(var(:,igg,jhh),minprof,maxprof,hx_g(igg,jhh),xe &
#if defined key_siggen 
      ,hcx_sig(igg,jhh) &
#endif
           )
   END IF

   ! Interpolation
   ibm_profuint=interpvit(px,py,val4(1,1),val4(2,1),val4(1,2),val4(2,2), &
        htx(igg,jbb),htx(idd,jbb),htx(igg,jhh),htx(idd,jhh))

   END FUNCTION ibm_profuint

!!======================================================================
   FUNCTION ibm_profvint(var,xpos,ypos,minprof,maxprof)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION ibm_profuint  ***
   !&E
   !&E ** Purpose :Interpolate vertically averaged current 
   !&E          from the 4 neighbooring location 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  04-2013 M. Huret 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : hty,hy_g
     USE trajectools, ONLY : interpvit
#if defined key_siggen 
     USE comsiggen, ONLY : hcy_sig
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: var
   REAL(KIND=rsh), INTENT( in )                :: xpos,ypos
   REAL(KIND=rsh), INTENT( in )                :: minprof, maxprof
   REAL(KIND=rsh)                              :: ibm_profvint

   !! * Local declarations
   REAL(KIND=rsh), DIMENSION(2,2)              :: val4
   INTEGER                    :: igg,idd,jhh,jbb
   INTEGER                    :: hlb,hlt,hrb,hrt  
   REAL(KIND=rsh)             :: xe,px,py


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

   !---- reperage de la maille ou se situe la particule pour les V
   ! location of ssh : xpos=REAL(i), ypos=REAL(j)
   igg=int(xpos)
   idd=igg+1
   jhh=nint(ypos)
   jbb=jhh-1

   idd=min(max(idd,imin),imax)
   igg=min(max(igg,imin),imax)
   jbb=min(max(jbb,jmin),jmax)
   jhh=min(max(jhh,jmin),jmax)

   ! calcul des positions px et py dans la maille V
   px=xpos-real(igg)
   py=ypos-real(jbb)-0.5_rsh


   val4(:,:)=0.0_rsh
   ! Processing prof mean on 4 neighbooring locations
   if (hty(igg,jbb)>0.0_rsh) THEN
      xe=hty(igg,jbb)-hy_g(igg,jbb)
      val4(1,1)=ibm_profmean(var(:,igg,jbb),minprof,maxprof,hy_g(igg,jbb),xe &
#if defined key_siggen 
           ,hcy_sig(igg,jbb) &
#endif
           )
   END IF
   if (hty(idd,jbb)>0.0_rsh) THEN
      xe=hty(idd,jbb)-hy_g(idd,jbb)
      val4(2,1)=ibm_profmean(var(:,idd,jbb),minprof,maxprof,hy_g(idd,jbb),xe &
#if defined key_siggen 
      ,hcy_sig(idd,jbb) &
#endif
           )
   END IF
   if (hty(idd,jhh)>0.0_rsh) THEN
      xe=hty(idd,jhh)-hy_g(idd,jhh)
      val4(2,2)=ibm_profmean(var(:,idd,jhh),minprof,maxprof,hy_g(idd,jhh),xe &
#if defined key_siggen 
      ,hcy_sig(idd,jhh) &
#endif
           )
   END IF
   if (hty(igg,jhh)>0.0_rsh) THEN
      xe=hty(igg,jhh)-hy_g(igg,jhh)
      val4(1,2)=ibm_profmean(var(:,igg,jhh),minprof,maxprof,hy_g(igg,jhh),xe &
#if defined key_siggen 
      ,hcy_sig(igg,jhh) &
#endif
           )
   END IF

   ! Interpolation
   ibm_profvint=interpvit(px,py,val4(1,1),val4(2,1),val4(1,2),val4(2,2), &
        hty(igg,jbb),hty(idd,jbb),hty(igg,jhh),hty(idd,jhh))

   END FUNCTION ibm_profvint

!!======================================================================
   FUNCTION ibm_traint(varint,xe,spos,khi,klo,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION traint  ***
   !&E
   !&E ** Purpose : Linear interpolation of tracer values 
   !&E              at the location (xpos,ypos,zpos)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  08-2002 (F. Dumas, P. Lazure) 
   !&E       !  11-2010 M. Huret 
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : ierrorlog
   USE comvars3d,  ONLY : sig,dsigw
#ifdef key_siggen
   USE trajectools, ONLY : h_int_siggen
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: varint
   REAL(KIND=rsh), INTENT( in )                :: spos,px,py
   REAL(KIND=rsh)                              :: ibm_traint
   INTEGER, INTENT( in )                       :: khi,klo
   INTEGER, INTENT( in )                       :: i1,i2,j1,j2
   INTEGER, INTENT( in )                       :: rh11,rh21,rh12,rh22
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe


   !! * Local declarations
   REAL(KIND=rsh)                              :: traintm,traintp
   REAL(KIND=rsh)                              :: xe_lag,hc_sig_lag,h0_lag
   REAL(KIND=rsh)                              :: f_lag_hzp,f_lag_hzm,f_lag_hz
   REAL(KIND=rsh), DIMENSION(2,2)              :: val4
   INTEGER                                     :: IERR_MPI


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

   
   if ((khi>kmax+1) .or. (klo>kmax+1) .or. (khi<0) .or. (klo<0)) then
      WRITE(ierrorlog,*) 'Function IBM_TRAINT'
      WRITE(ierrorlog,*) 'bad layer from sigma'
      WRITE(ierrorlog,*) 'The simulation is stopped'
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      STOP
   end if

   !-- interpolation proprement dite sur la nappe sigma au-dessus
   IF (khi /= kmax+1) THEN
      val4(:,:)=varint(khi,i1:i2,j1:j2)
      traintp=ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)
   ELSE
      val4(:,:)=varint(klo,i1:i2,j1:j2)
      ibm_traint=ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)
   END IF
   !-- interpolation proprement dite sur la nappe sigma au-dessous
   IF (klo /= 0) THEN
      val4(:,:)=varint(klo,i1:i2,j1:j2)
      traintm=ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)
   ELSE
      val4(:,:)=varint(khi,i1:i2,j1:j2)
      ibm_traint=ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)
   END IF

   IF (khi /= kmax+1 .and. klo /= 0) THEN
   !-- interpolation sur la verticale entre la valeur au-dessus (traintp)
   !-- et la valeur au-dessous (traintm).
        f_lag_hzp=1.0_rsh
        f_lag_hzm=1.0_rsh
        f_lag_hz=1.0_rsh
#ifdef key_siggen
       CALL h_int_siggen(xe_lag,h0_lag,hc_sig_lag,xe,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)
       CALL f_lag_sigz_tra(spos,klo,xe_lag,h0_lag,hc_sig_lag,f_lag_hzp,f_lag_hzm,f_lag_hz)
#endif
      ibm_traint=(f_lag_hzp*(sig(khi)-spos)*traintm+f_lag_hzm*(spos-sig(klo))*traintp) &
           /(f_lag_hz*dsigw(klo)) 
   ENDIF

   END FUNCTION ibm_traint


!!======================================================================
   FUNCTION ibm_proftraint(var,minprof,maxprof,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION traint  ***
   !&E
   !&E ** Purpose :Interpolate on the horizontal and average over minprof - maxprof
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  04-2013 M. Huret 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : h0_g
     USE comvarp2d,  ONLY : ssh
#if defined key_siggen 
     USE comsiggen, ONLY : hc_sig
#endif


   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT( in )  :: var
   REAL(KIND=rsh), INTENT( in )                :: px,py,minprof, maxprof
   REAL(KIND=rsh)                              :: ibm_proftraint
   INTEGER, INTENT( in )                       :: i1,i2,j1,j2
   INTEGER, INTENT( in )                       :: rh11,rh21,rh12,rh22

   !! * Local declarations
   REAL(KIND=rsh), DIMENSION(2,2)              :: val4
   
   !!----------------------------------------------------------------------
   !! * Executable part

   val4(:,:)=0.0_rsh

   ! Processing prof mean on 4 neighbooring locations
   if (rh11==1) THEN
      val4(1,1)=ibm_profmean(var(:,i1,j1),minprof,maxprof,h0_g(i1,j1),ssh(i1,j1) &
#if defined key_siggen 
      ,hc_sig(i1,j1) &
#endif
           )
   END IF
   if (rh21==1) THEN
      val4(2,1)=ibm_profmean(var(:,i2,j1),minprof,maxprof,h0_g(i2,j1),ssh(i2,j1) &
#if defined key_siggen 
      ,hc_sig(i2,j1) &
#endif
           )
   END IF
   if (rh22==1) THEN
      val4(2,2)=ibm_profmean(var(:,i2,j2),minprof,maxprof,h0_g(i2,j2),ssh(i2,j2) &
#if defined key_siggen 
      ,hc_sig(i2,j2) &
#endif
           )
   END IF
   if (rh12==1) THEN
      val4(1,2)=ibm_profmean(var(:,i1,j2),minprof,maxprof,h0_g(i1,j2),ssh(i1,j2) &
#if defined key_siggen 
      ,hc_sig(i1,j2) &
#endif
           )
   END IF

   ! Interpolation
   ibm_proftraint=ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)

   END FUNCTION ibm_proftraint

!!======================================================================
   FUNCTION ibm_profmean(var,minprof,maxprof,h,xe &
#if defined key_siggen
        ,hc_sig &
#endif
)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION ibm_profmean  ***
   !&E
   !&E ** Purpose : Average of values of var over a given prof range (minprof, maxprof)
   !&E               minprof and maxprof are provided as immersion values (positive)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  04-2013 M. Huret 
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars3d,  ONLY : sigw
   USE trajectools, ONLY : ksupkinf
#ifdef key_siggen
   USE trajectools, ONLY : siggentoz
#endif

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax), INTENT( in )  :: var
   REAL(KIND=rsh), INTENT( in )                :: minprof, maxprof,h,xe
   REAL(KIND=rsh)                              :: ibm_profmean
#if defined key_siggen
   REAL(KIND=rsh)                              ::hc_sig
#endif

   !! * Local declarations
   REAL(KIND=rsh)                              :: zmin, zmax, valint,d3
   REAL(KIND=rsh), DIMENSION(0:kmax)           :: zpos
   INTEGER                                     :: kw,kwm,l,ksub,ksup
   
!---------------------------------------------------------
   !! * Executable part

   valint=0.0_rsh 

   ! switch from immersion to real z
   d3=h+xe
   zmin=-minprof+xe  
   zmax=-maxprof+xe

  ! To make sure minprof is smaller than maxprof
   if (minprof>=maxprof) THEN 
      write(*,*) 'attention minprof > maxprof', minprof, maxprof, d3
   end if

   ! Convert to z because we have minprof and maxprof in z
   DO l=0,kmax
#if defined key_siggen
      CALL siggentoz(zpos(l),sigw(l),xe,h,hc_sig)
#else
      zpos(l)= sigw(l)*d3+xe
#endif
   END DO

!  Calcul de la prof sigma min et max
   if (minprof >= d3) THEN
      ksup=1
      zmin=zpos(1)
   else
      kw=kmax
      kwm=0
      CALL ksupkinf(zmin,zpos,kmax+1,kw,kwm,10)
      ksup=kw
   end if
   
   if (maxprof >= d3) THEN
      ksub=0
      zmax=zpos(0)
   else
      kw=kmax
      kwm=0
      CALL ksupkinf(zmax,zpos,kmax+1,kw,kwm,11)
      ksub=kwm
   end if

   ! Cumul des valeurs entre minprof et maxprof
   if (ksup==ksub+1) THEN
      valint=var(ksup)*(zmin-zmax)
   else if (ksup==ksub+2) THEN
      valint=valint+var(ksup)*(zmin-zpos(ksup-1))
      valint=valint+var(ksup-1)*(zpos(ksup-1)-zmax)
   else
      valint=valint+var(ksup)*(zmin-zpos(ksup-1))
      DO l=ksub+2,ksup-1 
         valint=valint+var(l)*(zpos(l)-zpos(l-1))
      END DO
      valint=valint+var(ksub+1)*(zpos(ksub+1)-zmax)
   endif

   ! Average over minprof - maxprof
   ibm_profmean=valint/(zmin-zmax)

   END FUNCTION ibm_profmean

 !!======================================================================
   FUNCTION ibm_interp(val4,px,py,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  SUBROUTINE interp  ***
   !&E
   !&E ** Purpose : Horizontal interpolation with verification of boundary limits
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  01-2011 (M. Huret)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(KIND=rsh)                 :: ibm_interp
   REAL(KIND=rsh), DIMENSION(2,2), INTENT( in )  :: val4
   REAL(KIND=rsh), INTENT( in )   :: px,py
   INTEGER, INTENT( in )          :: rh11,rh21,rh12,rh22

   !! * Local declarations

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

    ! check if boundary condition
   IF (rh11 > 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 > 0) THEN
      ibm_interp=px*(1-py)*val4(2,1) + px*py  *val4(2,2) &
           +(1-px)*py*val4(1,2) + (1-px)*(1-py) *val4(1,1)

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
      ibm_interp=0

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
      ibm_interp=val4(2,2)

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
      ibm_interp=val4(1,2)

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 > 0) THEN
      ibm_interp=px*val4(2,2)+(1-px)*val4(1,2)
      
   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
      ibm_interp=val4(2,1)
      
   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
      if (px > 0.5_rsh) then
         ibm_interp=val4(2,1) 
      else
         ibm_interp=val4(1,2)
      end if
      
   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
      ibm_interp=(1-py)*val4(2,1) + py*val4(2,2)
      
   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 > 0) THEN
      ibm_interp=px*(val4(2,2)-val4(1,2))+py*(val4(2,2)-val4(2,1))-val4(2,2) &
           +val4(1,2)+val4(2,1)
      
   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
      ibm_interp=val4(1,1)
      
   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
      if (px > 0.5_rsh) then
         ibm_interp=val4(2,2)
      else
         ibm_interp=val4(1,1)
      end if
      
   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
      ibm_interp=px*val4(2,1)+(1-px)*val4(1,1)
      
   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
      ibm_interp=px*(val4(2,1)-val4(1,1))+py*(val4(2,2)-val4(2,1))+val4(1,1)
      
   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
      ibm_interp=py*val4(1,2) + (1-py)*val4(1,1)
      
   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
      ibm_interp=px*(val4(2,1)-val4(1,1))+py*(val4(1,2)-val4(1,1))+val4(1,1)
      
   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 > 0) THEN
      ibm_interp=px*(val4(2,2)-val4(1,2))+py*(val4(1,2)-val4(1,1))+val4(1,1)
      
   END IF
   
 END FUNCTION ibm_interp
 
   !!======================================================================
   SUBROUTINE ibm_loc_xyz(xe,xpos,ypos,spos,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22,khi,klo)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE location  ***
   !&E
   !&E ** Purpose : search index of the cells around particle
   !&E              and location of the particle inside the cell 
   !&E              for later interpolation 
   !&E              Valid relative to h0,xe,wkz... 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E               12/07 M. Huret
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,   ONLY : h0_g
   USE parameters,  ONLY : valmanq
   USE trajectools, ONLY : ksupkinf
     
   !! * Arguments
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rsh), INTENT( in )       :: xpos,ypos,spos
   REAL(KIND=rsh), INTENT( out )      :: px,py 
   INTEGER, INTENT( out )             :: i1,i2,j1,j2
   INTEGER, INTENT( out )             :: rh11,rh21,rh12,rh22
   INTEGER, INTENT( out )             :: khi,klo
     
   !! * Local declarations
   REAL(KIND=rsh)                         :: xe_lag,h0_lag
   REAL(KIND=rsh), DIMENSION(0:kmax+1)    :: spos_w  

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

   ! indexes of the cell where the particle is
   ! reperage de la maille XE ou se situe la particule

   i1=INT(xpos)
   i2=INT(xpos)+1
   j1=INT(ypos)
   j2=INT(ypos)+1

   i1 = MIN(MAX(i1,imin),imax)
   i2 = MIN(MAX(i2,imin),imax)
   j1 = MIN(MAX(j1,jmin),jmax)
   j2 = MIN(MAX(j2,jmin),jmax)

   ! location of the particle inside the cell
   ! calcul des positions px et py dans la maille XE

   px=xpos-INT(xpos)
   py=ypos-INT(ypos)

   ! check for next cells in water
   rh11 = 1
   rh12 = 1
   rh21 = 1
   rh22 = 1

   IF (h0_g(i1,j1)==-valmanq)  rh11 = 0
   IF (h0_g(i1,j2)==-valmanq)  rh12 = 0
   IF (h0_g(i2,j1)==-valmanq)  rh21 = 0
   IF (h0_g(i2,j2)==-valmanq)  rh22 = 0

   ! get neighbooring indices on the vertical
   spos_w(0)=-1.0_rsh
   spos_w(kmax+1)=0.0_rsh
   spos_w(1:kmax)=sig(1:kmax)

   klo=0
   khi=kmax+1
   CALL ksupkinf(spos,spos_w,kmax+2,khi,klo,6)

 END SUBROUTINE ibm_loc_xyz


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

 SUBROUTINE gasdev_s(harvest)

   REAL(rsh), INTENT(OUT) :: harvest
   !Returns in harvest a normally distributed deviate with zero mean and unit variance, using
   !ran1 as the source of uniform deviates.
   REAL(rsh) :: rsq,v1,v2
   REAL(rsh), SAVE :: g
   LOGICAL, SAVE :: gaus_stored=.false.
   
   if (gaus_stored) then ! We have an extra deviate handy,
      harvest=g          !  so return it,
      gaus_stored=.false. !and unset the flag.
   else !We dont have an extra deviate handy, so
      do
         call random_number(v1) ! pick two uniform numbers in the square extending from 
         call random_number(v2) !-1 to +1 in each direction,
         v1=2.0_rsh*v1-1.0_rsh
         v2=2.0_rsh*v2-1.0_rsh
         rsq=v1**2+v2**2 !see if they are in the unit circle,
         if (rsq > 0.0 .and. rsq < 1.0) exit
      end do !otherwise try again.
      rsq=sqrt(-2.0_rsh*log(rsq)/rsq) !Now make the Box-Muller transformation to get two normal deviates. Return one and save the other for next time.
      harvest=v1*rsq
      g=v2*rsq
      gaus_stored=.true. ! Set flag.
   end if

 END SUBROUTINE gasdev_s


#ifdef key_siggen

!!======================================================================
   SUBROUTINE f_lag_sigz_tra(spos,k,xe_lag,h0_lag,hc_sig_lag,f_lag_hzp,f_lag_hzm,f_lag_hz)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE f_lag  ***
   !&E
   !&E ** Purpose : Scale factor from generalized sigma to z
   !&E              for vertical interpolation at particle location
   !&E
   !&E ** Description :
   !&Ex3
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History : M. Huret 12/1010
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comsiggen, ONLY : b_sig,theta_sig,Csu_sig,dcwsds
   USE comvars3d,  ONLY : sig

   !! * Arguments
   REAL(KIND=rsh), INTENT( in ) :: spos,xe_lag,h0_lag,hc_sig_lag
   REAL(KIND=rsh), INTENT(out ) :: f_lag_hzp,f_lag_hzm,f_lag_hz
   INTEGER,        INTENT( IN ) :: k

   !! * Local declarations
   REAL(KIND=rsh)      :: Csusig, dcwsdsp,dcwsdsm
   
   !!----------------------------------------------------------------------
   !! * Executable part

   Csusig=(1.0_rsh-b_sig)*SINH(theta_sig*spos)/SINH(theta_sig) &
         + 0.5_rsh*b_sig*(TANH(theta_sig*(spos+0.5_rsh))-TANH(0.5_rsh*theta_sig))/TANH(theta_sig/2.0_rsh)

   dcwsdsp=0.0_rsh
   dcwsdsm=0.0_rsh

   IF (spos /= sig(k+1)) dcwsdsp = (Csu_sig(k+1)-Csusig)/(sig(k+1)-spos)
   IF (spos /= sig(k))   dcwsdsm = (Csusig-Csu_sig(k))/(spos-sig(k))

   f_lag_hzp = xe_lag + hc_sig_lag + (h0_lag-hc_sig_lag)*dcwsdsp
   f_lag_hzm = xe_lag + hc_sig_lag + (h0_lag-hc_sig_lag)*dcwsdsm
   f_lag_hz  = xe_lag + hc_sig_lag + (h0_lag-hc_sig_lag)*dcwsds(k)

  END SUBROUTINE f_lag_sigz_tra

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

   ! key_siggen
#endif

   ! key_ibm  
#endif

 END MODULE
