MODULE trajectools

#if defined key_trajec3d || defined key_trajec2d || defined key_ibm
   !!======================================================================
   !!                   ***  MODULE trajectools  ***
   !! Trajectories dynamics:  routines for trajectory estimates
   !! Several routines called by traject3d
   !!======================================================================

#include "toolcpp.h"

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

   IMPLICIT NONE
   PRIVATE

   !! * Accessibility
   PUBLIC h0int,xeint,ksupkinf,loc_h0,kzprofile,splint,uint,vint,wint,dksdzint,interpvit
#if defined key_siggen
   PUBLIC siggentoz,ztosiggen,hc_sigint,h_int_siggen
#endif

 CONTAINS

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

   FUNCTION uint(xpos,ypos,spos,u,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION uint  ***
   !&E
   !&E ** Purpose : Linear interpolation of zonal velocity (uint)
   !&E              at the location (xpos,ypos,spos)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : avance.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  08-2002 (F. Dumas, P. Lazure) 
   !&E       !  01-2007 (V. Garnier) correction bug declaration px,py
   !&E       !  01-2011 (M. Huret) Adaptation to siggen (scale factor for vertical interpolation)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,  ONLY : htx,l_modele2d
   USE comvars3d,  ONLY : sig,dsigw

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )                :: xpos,ypos,spos
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in ) :: u
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rsh)                              :: uint

   !! * Local declarations
   INTEGER                    :: igg,idd,jhh,jbb,k,kuv,kuvm
   INTEGER                    :: hlb,hlt,hrb,hrt  

   REAL(KIND=rsh)             :: uintm,uintp,px,py
   REAL(KIND=rsh)             ::f_lag_hzp,f_lag_hzm,f_lag_hz
   REAL(KIND=rsh)             ::xe_lag,hc_sig_lag,h0_lag

   !!----------------------------------------------------------------------
   !! * 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)

   IF (l_modele2d) THEN
      uint=interpvit(px,py,u(1,igg,jbb),u(1,idd,jbb),u(1,igg,jhh),u(1,idd,jhh),&
                            htx(igg,jbb),htx(idd,jbb),htx(igg,jhh),htx(idd,jhh)) 
   ELSE    ! not modele2d

     DO k=kmax,1,-1
        IF (spos <= sig(k).and.spos >= sig(k-1)) THEN
           kuv=k
           kuvm=k-1
        END IF
     END DO
     IF (spos <= sig(1)) THEN
        kuv=1
        kuvm=0
     END IF
     IF (spos >= sig(kmax)) THEN
        kuv=kmax
        kuvm=kmax
     END IF

     !-- interpolation proprement dite sur la nappe sigma au-dessus
     uintp=interpvit(px,py,u(kuv,igg,jbb),u(kuv,idd,jbb), &
                            u(kuv,igg,jhh),u(kuv,idd,jhh), &
                            htx(igg,jbb),htx(idd,jbb), &
                            htx(igg,jhh),htx(idd,jhh))

     !-- interpolation proprement dite sur la nappe sigma au-dessous
     IF (kuvm /= 0) THEN
       uintm=interpvit(px,py,u(kuvm,igg,jbb),u(kuvm,idd,jbb), &
                             u(kuvm,igg,jhh),u(kuvm,idd,jhh), &
                             htx(igg,jbb),htx(idd,jbb), &
                             htx(igg,jhh),htx(idd,jhh))
     ELSE
       uintm=0.0_rsh
     END IF

     !-- interpolation sur la verticale entre la valeur au-dessus (uintp)
     !-- et la valeur au-dessous (uintm)
     IF (spos >= sig(kmax)) THEN
       uint=uintm
     ELSE IF(kuv == 1) THEN
       uint = ((sig(kuv)-spos)*uintm+(spos+1.0_rsh)*uintp)/(sig(kuv)+1.0_rsh)
     ELSE
        f_lag_hzp=1.0_rsh
        f_lag_hzm=1.0_rsh
        f_lag_hz=1.0_rsh
#ifdef key_siggen
       CALL loc_h0(xpos,ypos,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
       CALL h_int_siggen(xe_lag,h0_lag,hc_sig_lag,xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
     !  h0_lag=h0int(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
     !  xe_lag=xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
     !  hc_sig_lag=hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
       CALL f_lag_sigz_uv(spos,kuvm,xe_lag,h0_lag,hc_sig_lag,f_lag_hzp,f_lag_hzm,f_lag_hz)
#endif
        uint=(f_lag_hzp*(sig(kuv)-spos)*uintm+f_lag_hzm*(spos-sig(kuvm))*uintp) &
             /(f_lag_hz*dsigw(kuvm)) 

     END IF


   END IF   ! Ends test on l_modele2d

   END FUNCTION uint

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

   FUNCTION vint(xpos,ypos,spos,v,xe)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION vint  ***
   !&E
   !&E ** Purpose : Linear interpolation of meridional velocity (vint)
   !&E              at the location (xpos,ypos,spos)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : avance.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  08-2002 (F. Dumas, P. Lazure)
   !&E       !  01-2007 (V. Garnier) correction bug declaration px,py
   !&E       !  01-2011 (M. Huret) Adaptation to siggen (scale factor for vertical interpolation)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d,  ONLY : hty,l_modele2d
   USE comvars3d,  ONLY : sig,dsigw

   !! * Arguments
   REAL(KIND=rsh),                                                  INTENT(in) :: xpos,ypos,spos
   REAL(KIND=rsh), DIMENSION(kmax,liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: v
   REAL(KIND=rsh), DIMENSION(     liminm2:limaxp2,ljminm2:ljmaxp2), INTENT(in) :: xe
   REAL(KIND=rsh) :: vint   ! Result

   !! * Local declarations
   INTEGER                    :: idd,igg,jhh,jbb,k,kuv,kuvm
   INTEGER                    :: hlb,hlt,hrb,hrt  

   REAL(KIND=rsh)             :: vintm,vintp,px,py
   REAL(KIND=rsh)             :: f_lag_hzp,f_lag_hzm,f_lag_hz
   REAL(KIND=rsh)             :: xe_lag,hc_sig_lag,h0_lag

   !!----------------------------------------------------------------------
   !! * 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

   IF (l_modele2d) THEN
      vint=interpvit(px,py,v(1,igg,jbb),v(1,idd,jbb),v(1,igg,jhh),v(1,idd,jhh), &
                            hty(igg,jbb),hty(idd,jbb),hty(igg,jhh),hty(idd,jhh))

   ELSE   ! not l_modele2d

     DO k=kmax,1,-1
       IF (spos <= sig(k).and.spos >= sig(k-1)) THEN
         kuv=k
         kuvm=k-1
       END IF
     END DO
     IF (spos <= sig(1)) THEN
       kuv=1
       kuvm=0
     END IF
     IF (spos >= sig(kmax)) THEN
       kuv=kmax
       kuvm=kmax
     END IF

     !-- interpolation proprement dite sur la nappe sigma au-dessus
     vintp=interpvit(px,py,v(kuv,igg,jbb),v(kuv,idd,jbb), &
                          v(kuv,igg,jhh),v(kuv,idd,jhh), &
                          hty(igg,jbb),hty(idd,jbb), &
                          hty(igg,jhh),hty(idd,jhh))

     !-- interpolation proprement dite sur la nappe sigma au-dessous
     IF (kuvm /= 0) THEN
       vintm=interpvit(px,py,v(kuvm,igg,jbb),v(kuvm,idd,jbb), &
                             v(kuvm,igg,jhh),v(kuvm,idd,jhh), &
                                   hty(igg,jbb),hty(idd,jbb), &
                                   hty(igg,jhh),hty(idd,jhh))
     ELSE
       vintm=0.0_rsh
     END IF

     !-- interpolation sur la verticale entre la valeur au-dessus (vintp)
     !-- et la valeur au-dessous (vintm).
     IF (spos > sig(kmax)) THEN
       vint=vintm
     ELSE IF(kuv == 1) THEN
       vint=((sig(kuv)-spos)*vintm+(spos+1.0_rsh)*vintp)/(sig(kuv)+1.0_rsh)
     ELSE
        f_lag_hzp=1.0_rsh
        f_lag_hzm=1.0_rsh
        f_lag_hz=1.0_rsh
#ifdef key_siggen
        CALL loc_h0(xpos,ypos,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
        CALL h_int_siggen(xe_lag,h0_lag,hc_sig_lag,xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
      !  h0_lag=h0int(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
      !  xe_lag=xeint(xe,px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
      !  hc_sig_lag=hc_sigint(px,py,igg,idd,jbb,jhh,hlb,hrb,hlt,hrt)
        CALL f_lag_sigz_uv(spos,kuvm,xe_lag,h0_lag,hc_sig_lag,f_lag_hzp,f_lag_hzm,f_lag_hz)
#endif
        vint=(f_lag_hzp*(sig(kuv)-spos)*vintm+f_lag_hzm*(spos-sig(kuvm))*vintp) &
             /(f_lag_hz*dsigw(kuvm)) 

     END IF
 
   END IF   ! Ends test on l_modele2d

   END FUNCTION vint


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

   FUNCTION interpvit(x,y,v00,v10,v01,v11,rh00,rh10,rh01,rh11)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION interpvit  ***
   !&E
   !&E ** Purpose : interpolation horizontale d une composante de vitesse
   !&E              tenant compte des profondeurs aux points encadrant la
   !&E              position ou l on doit interpoler.
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : uint.F90,vint.F90
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  08-2002 (F. Dumas, P. Lazure)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )      :: x,y,v00,v10,v01,v11,rh00,rh10,rh01,rh11
   REAL(KIND=rsh)                    :: interpvit

   !! * Local declarations
   INTEGER        :: icas
   REAL(KIND=rsh) :: xref,yref

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

   IF ( rh00 > 0.0_rsh .and. rh10 > 0.0_rsh .and. rh01 > 0.0_rsh .and. rh11 > 0.0_rsh ) THEN
     interpvit= x*y*v11+x*(1.0_rsh-y)*v10+(1.0_rsh-x)*y*v01+(1.0_rsh-x)*(1.0_rsh-y)*v00
     icas = 1
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=0.0_rsh
     xref=0.5_rsh
     yref=0.5_rsh
     icas = 0
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=v11
     xref=1.0_rsh
     yref=1.0_rsh
     icas = 2
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=v10
     xref=1.0_rsh
     yref=0.0_rsh
     icas = 4
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=(1.0_rsh-y)*v10+y*v11
     xref=1.0_rsh
     yref=0.5_rsh
     icas = 3
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=v01
     xref=0.0_rsh
     yref=1.0_rsh
     icas = 5
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=((v10-v01)/2.0_rsh)*(x-y)+(v10+v01)/2.0_rsh
     xref=0.5_rsh
     yref=0.5_rsh
     icas = 6
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=(1.0_rsh-x)*v01+x*v11
     xref=0.5_rsh
     yref=1.0_rsh
     icas = 7
   ELSE IF (rh00 <= 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=(v11-v01)*x+(v11-v10)*y+v10-(v11-v01)
     xref=1.0_rsh
     yref=1.0_rsh
     icas = 8
   ELSE IF (rh00 > 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=v00
     xref=0.0_rsh
     yref=0.0_rsh
     icas = 9
   ELSE IF (rh00 > 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=((v11-v00)/2.0_rsh)*(x+y)+v00
     xref=0.5_rsh
     yref=0.5_rsh
     icas = 10
   ELSE IF (rh00 > 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=(1.0_rsh-y)*v00+y*v01
     xref=0.0_rsh
     yref=0.5_rsh
     icas = 11
   ELSE IF (rh00 > 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 <= 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=(v11-v01)*x+(v01-v00)*y+v00
     xref=0.0_rsh
     yref=1.0_rsh
     icas = 12
   ELSE IF (rh00 > 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=(1.0_rsh-x)*v00+x*v10
     xref=0.5_rsh
     yref=0.0_rsh
     icas = 13
   ELSE IF (rh00 > 0.0_rsh .and. rh01 > 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 <= 0.0_rsh) THEN
     interpvit=(v10-v00)*x+(v01-v00)*y+v00
     xref=0.0_rsh
     yref=0.0_rsh
     icas = 14
   ELSE IF (rh00 > 0.0_rsh .and. rh01 <= 0.0_rsh .and. rh10 > 0.0_rsh .and. rh11 > 0.0_rsh) THEN
     interpvit=(v10-v00)*x+(v11-v10)*y+v00
     xref=1.0_rsh
     yref=0.0_rsh
     icas = 15
   END IF
!
   END FUNCTION interpvit

    !!======================================================================
   SUBROUTINE ksupkinf(pos,refpos,sizerefpose,kw,kwm,appel)

   !&E---------------------------------------------------------------------
   !&E                 ***  SUBROUTINE ksupkinf ***
   !&E
   !&E ** Purpose : Seek neighboor index in the vertical for a particle position
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d.F90, wint, dksdzint, splint, ztosiggen
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   
   REAL(KIND=rsh),                            INTENT(in)    :: pos
   INTEGER,                                   INTENT(in)    :: sizerefpose,appel
   REAL(KIND=rsh),DIMENSION(0:sizerefpose-1), INTENT(in)    :: refpos
   INTEGER,                                   INTENT(inout) :: kw,kwm

   !! * Local declarations
   INTEGER                                     :: k
   INTEGER                               :: IERR_MPI

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

       do while (kw-kwm> 1) 
        k=int((kw+kwm)/2)  
        if(refpos(k) > pos)then
          kw=k
        else
          kwm=k
        endif
      end do

      IF ((refpos(kw)<pos) .OR. (refpos(kwm)>(pos+0.00001_rsh))) THEN
         print*,'erreur subroutine ksupkinf ',appel
         print*,'kw',kw
         print*,'kwm',kwm
         print*,'pos',pos
         print*,'refpos(1)',refpos(1)
         print*,'refpos(kw)',refpos(kw)
         print*,'refpos(kwm)',refpos(kwm)
         print*,'refpos',refpos
         print*,'refpos',appel
         CALL_MPI MPI_FINALIZE(IERR_MPI)
         stop
      END IF

   END SUBROUTINE ksupkinf


   !!======================================================================
   SUBROUTINE loc_h0(xpos,ypos,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&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,wz,kz... 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d, traj_init3d, uint, vint
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E               12/07 M. Huret
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
     USE comvars2d,  ONLY : h0_g
     USE parameters, ONLY : valmanq

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )       :: xpos,ypos
   INTEGER, INTENT( out )             ::i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh), INTENT( out )      :: px,py 

   !! * Local declarations

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

   ! indexes of the cell where the particle is
   ! reperage de la maille XE ou se situe la particule
   ! location of ssh : xpos=REAL(i), ypos=REAL(j)

   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) THEN
     rh11=0
   END IF
   IF (h0_g(i1,j2)==-valmanq) THEN
     rh12=0
   END IF
   IF (h0_g(i2,j1)==-valmanq) THEN
     rh21=0
   END IF
   IF (h0_g(i2,j2)==-valmanq) THEN
     rh22=0
   END IF

 END SUBROUTINE loc_h0

 !!======================================================================
   FUNCTION h0int(px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE h0int  ***
   !&E
   !&E ** Purpose : calcul par interpolation lineaire de la bathymetry
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (P. Lazure, F. Dumas)  Original code
   !&E       !  08-2008 (F. Dumas ?)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY  : h0_g

   !! * Arguments
   REAL(KIND=rsh),                                 INTENT( in )  :: px,py
   INTEGER,                                        INTENT( in )  :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh)                                                :: h0int

   !! * Local declarations
   
   !!----------------------------------------------------------------------
   !! * Executable part

   ! spatial interpolation

   h0int =px*(1.0_rsh-py)*h0_g(i2,j1)*real(rh21)+ px*py*h0_g(i2,j2)*real(rh22) &
         +(1.0_rsh-px)*py*h0_g(i1,j2)*real(rh12)+(1.0_rsh-px)*(1.0_rsh-py)*h0_g(i1,j1)*real(rh11)

   END FUNCTION h0int

 !!======================================================================
   FUNCTION xeint(xe,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE xe_int  ***
   !&E
   !&E ** Purpose : calcul par interpolation lineaire de l elevation hauteur d eau
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (P. Lazure, F. Dumas)  Original code
   !&E       !  08-2008 (F. Dumas ?)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(KIND=rsh),                                 INTENT( in )  :: px,py
   INTEGER,                                        INTENT( in )  :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rsh)                                                :: xeint

   !! * Local declarations
   
   !!----------------------------------------------------------------------
   !! * Executable part

   ! spatial interpolation

   xeint =px*(1.0_rsh-py)*xe(i2,j1)*REAL(rh21,rsh)+ px*py*xe(i2,j2)*REAL(rh22,rsh) &
         +(1.0_rsh-px)*py*xe(i1,j2)*REAL(rh12,rsh)+(1.0_rsh-px)*(1.0_rsh-py)*xe(i1,j1)*REAL(rh11,rsh)


   END FUNCTION xeint

   !!======================================================================
   FUNCTION interp(varint,k,i1,i2,j1,j2,px,py,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  SUBROUTINE interp  ***
   !&E
   !&E ** Purpose : Interpolate at H0 or tracer location
   !&E              with check of terrestrial boundary 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : wint
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  01-2008 (M. Huret)
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
   REAL(KIND=rsh)                 :: interp
   REAL(KIND=rsh), DIMENSION(0:kmax,liminm1:limaxp1,ljminm1:ljmaxp1),INTENT( in ) :: varint  
   REAL(KIND=rsh), INTENT( in )   :: px,py
   INTEGER, INTENT( in )          :: k
   INTEGER, INTENT( in )          :: i1,j1,i2,j2
   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
     interp=px*(1.0_rsh-py)*varint(k,i2,j1) + px*py         *varint(k,i2,j2) &
           +(1.0_rsh-px)*py*varint(k,i1,j2) + (1.0_rsh-px)*(1.0_rsh-py) *varint(k,i1,j1)

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

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
     interp=varint(k,i2,j2)

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
     interp=varint(k,i1,j2)

   ELSE IF (rh11 <= 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 > 0) THEN
     interp=px*varint(k,i2,j2)+(1.0_rsh-px)*varint(k,i1,j2)

   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
     interp=varint(k,i2,j1)

   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
     if (px > 0.5_rsh) then
       interp=varint(k,i2,j1) 
     else
       interp=varint(k,i1,j2)
     end if

   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
     interp=(1.0_rsh-py)*varint(k,i2,j1) + py*varint(k,i2,j2)

   ELSE IF (rh11 <= 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 > 0) THEN
    interp=px*(varint(k,i2,j2)-varint(k,i1,j2))+py*(varint(k,i2,j2)-varint(k,i2,j1))-varint(k,i2,j2)+varint(k,i1,j2)+varint(k,i2,j1)
   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
     interp=varint(k,i1,j1)

   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
     if (px > 0.5_rsh) then
       interp=varint(k,i2,j2)
     else
       interp=varint(k,i1,j1)
     end if

   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 <= 0) THEN
     interp=px*varint(k,i2,j1)+(1.0_rsh-px)*varint(k,i1,j1)

   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 <= 0 .and. rh22 > 0) THEN
     interp=px*(varint(k,i2,j1)-varint(k,i1,j1))+py*(varint(k,i2,j2)-varint(k,i2,j1))+varint(k,i1,j1)

   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
     interp=py*varint(k,i1,j2) + (1.0_rsh-py)*varint(k,i1,j1)

   ELSE IF (rh11 > 0 .and. rh21 > 0 .and. rh12 > 0 .and. rh22 <= 0) THEN
     interp=px*(varint(k,i2,j1)-varint(k,i1,j1))+py*(varint(k,i1,j2)-varint(k,i1,j1))+varint(k,i1,j1)

   ELSE IF (rh11 > 0 .and. rh21 <= 0 .and. rh12 > 0 .and. rh22 > 0) THEN
     interp=px*(varint(k,i2,j2)-varint(k,i1,j2))+py*(varint(k,i1,j2)-varint(k,i1,j1))+varint(k,i1,j1)

   END IF

   END FUNCTION interp


   !!======================================================================
   FUNCTION wint(spos_s,xelag,h0lag,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  FUNCTION wint  ***
   !&E
   !&E ** Purpose : Linear interpolation of vertical velocity (wint)
   !&E              at the location (spos_s) between indices i1,i2,j1,j2
   !&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       !  01-2011 (M. Huret) Adaptation to siggen 
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars3d,  ONLY : dsigu,wz,sigw

   !! * Arguments
   REAL(KIND=rsh), INTENT( in )                :: spos_s,xelag,h0lag,px,py
   INTEGER, INTENT( in )                       :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh)                              :: wint

   !! * Local declarations
   INTEGER                               :: IERR_MPI
   INTEGER                               :: kw,kwm
   REAL(KIND=rsh)                        :: wintm,wintp,f_lag_hzp,f_lag_hzm,f_lag_hz,hc_siglag
   
   !!----------------------------------------------------------------------
   !! * Executable part

   kwm=0
   kw=kmax
   CALL ksupkinf(spos_s,sigw,kmax+1,kw,kwm,5)
   
   if ((kw>kmax) .or. (kwm>kmax) .or. (kw<0) .or. (kwm<0)) then
      write(*,*) 'mauvais sigma fonction wint'
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      stop
   end if
   
  !-- interpolation proprement dite sur la nappe sigma au-dessus
   wintp=interp(wz,kw,i1,i2,j1,j2,px,py,rh11,rh21,rh12,rh22)

   !-- interpolation proprement dite sur la nappe sigma au-dessous
   IF (kwm /= 0) THEN
      wintm=interp(wz,kwm,i1,i2,j1,j2,px,py,rh11,rh21,rh12,rh22)
   ELSE
      wintm=0.0_rsh
   END IF

   !-- interpolation sur la verticale entre la valeur au-dessus (wintp)
   !-- et la valeur au-dessous (wintm).
   f_lag_hzp=1.0_rsh
   f_lag_hzm=1.0_rsh
   f_lag_hz=1.0_rsh
#ifdef key_siggen
   hc_siglag=hc_sigint(px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)
   CALL f_lag_sigz_w(spos_s,kw,xelag,h0lag,hc_siglag,f_lag_hzp,f_lag_hzm,f_lag_hz)
#endif
   wint=(f_lag_hzp*(sigw(kw)-spos_s)*wintm+f_lag_hzm*(spos_s-sigw(kwm))*wintp) &
   /(f_lag_hz*dsigu(kw)) 

   END FUNCTION wint


   !!======================================================================
   FUNCTION dksdzint(zpos_s,zpos,kzsmth,dkz2)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE dksdzint  ***
   !&E
   !&E ** Purpose : calcul de la derive de la composante verticale
   !&E              du melange turbulent (notee dksdzint) pour le profil
   !&E              kz (kzsmth) lisse calcule dans kzprofile dont derive
   !&E              seconde est dkz2
   !&E              adapte de splint
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  12/2007  (M. Huret)
   !&E      !  12/2010  (M. Huret) Switch to z
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   !! * Arguments
     REAL(KIND=rsh)                :: dksdzint
     REAL(KIND=rsh), INTENT( in )  :: zpos_s
     REAL(KIND=rsh), DIMENSION(0:kmax),INTENT( in )  :: kzsmth,dkz2,zpos

   !! * Local declarations
     INTEGER                 :: IERR_MPI
     INTEGER          :: k,khi,klo
     REAL(KIND=rsh)   :: a,b,dzw_local
     
   !!----------------------------------------------------------------------
   !! * Executable part
   
     klo=0
     khi=kmax
     CALL ksupkinf(zpos_s,zpos,kmax+1,khi,klo,3)

     IF ((khi>kmax) .OR. (klo>kmax) .OR. (khi<0) .OR. (klo<0)) THEN
       WRITE(*,*) 'mauvais sigma fonction dksdzint'
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       STOP
     END IF
     
     dzw_local=zpos(khi)-zpos(klo)
     IF (dzw_local<=0.0_rsh) THEN
       WRITE(*,*) khi,klo,zpos_s,'bad sigw input in splint'
       READ(*,*) 
     END IF
     
     a=(zpos(khi)-zpos_s)/dzw_local
     b=(zpos_s-zpos(klo))/dzw_local
     dksdzint =(-kzsmth(klo)+kzsmth(khi))/dzw_local+ &
     ((-3.0_rsh*a**2+1)*dkz2(klo)+(3.0_rsh*b**2-1)*dkz2(khi))*dzw_local/6.0_rsh
     
   END FUNCTION dksdzint

   !!======================================================================
   SUBROUTINE kzprofile(kzsmth,dkz2,zpos,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE kzprofile  ***
   !&E
   !&E ** Purpose : Calcul les derivees secondes du coeff de diffusivite 
   !&E              vertical en sigw a l aide de spline apres avoir filtre
   !&E              le profil kz
   !&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
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars3d,  ONLY : kz

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(0:kmax),INTENT( out )    :: dkz2,kzsmth
   REAL(KIND=rsh), DIMENSION(0:kmax),INTENT( in )    :: zpos
   INTEGER, INTENT( in )                              :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh),     INTENT( in )                  :: px,py

   !! * Local declarations
   INTEGER                 :: IERR_MPI
   INTEGER                            :: k
   REAL(KIND=rsh), DIMENSION(0:kmax)  :: kzint
   REAL(kind=rsh)                     :: big

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

   big=1.0E30

   ! horizontal interpolation at each sigma level
   kzint(0)=0.0_rsh
   kzint(kmax)=0.0_rsh

   DO k=1,kmax-1
      kzint(k)=interp(kz,k,i1,i2,j1,j2,px,py,rh11,rh21,rh12,rh22)
   END DO

   ! smooth of the kz profile with weighting by the distance
   IF (kmax <= 3) THEN
      PRINT*,'float interp code cant handle'
      PRINT*,'kmax.le.3'
      CALL_MPI MPI_FINALIZE(IERR_MPI)
      stop
   ENDIF
   kzsmth(0)=kzint(0) 
   kzsmth(kmax)=kzint(kmax)
   DO k=1,kmax-1
      kzsmth(k)=(4.0_rsh/6.0_rsh)*kzint(k)+(2.0_rsh/6.0_rsh)* &
           (kzint(k+1)*(zpos(k)-zpos(k-1))+kzint(k-1)*(zpos(k+1)-zpos(k)))/(zpos(k+1)-zpos(k-1))
   ENDDO

   CALL spline(zpos(0:kmax),kzsmth,kmax+1,big,big,dkz2) 

   END SUBROUTINE kzprofile

    !!======================================================================
   SUBROUTINE spline(x,y,n,yp1,ypn,y2)

  !&E---------------------------------------------------------------------
  !&E                 ***  ROUTINE spline  ***
  !&E
  !&E     From numerical recipies vol 2, but modified so that nmax=50
  !&E---------------------------------------------------------------------

  !! * Modules used

  !! * Arguments
     real(kind=rsh), INTENT(in)           :: yp1,ypn
     real(kind=rsh),dimension(n), INTENT(in) :: x,y
     real(kind=rsh),dimension(n), INTENT(out) :: y2
     integer, INTENT(IN)         :: n

  !! * Local declarations
     integer                     :: i,k
     real(kind=rsh),dimension(n) :: u
     real(kind=rsh)              :: p,qn,sig,un
!
! executable statement
!---------------------------------------------------------
!
      if (yp1 > .99e30) then
        y2(1)=0.0
        u(1)=0.0
      else
        y2(1)=-0.5
        u(1)=(3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
!
      do i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.0
        y2(i)=(sig-1.0)/p
        u(i)=(6.0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
      end do
!
      if (ypn > 0.99e30) then
        qn=0.0
        un=0.0
      else
        qn=0.5
        un=(3.0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
!
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.0)
!
      do  k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
      end do
!
      end subroutine spline

   !!======================================================================
  SUBROUTINE splint(xa,ya,y2a,n,x,y)
!
  !&E---------------------------------------------------------------------
  !&E                 ***  ROUTINE splint  ***
  !&E
  !&E     From numerical recipies vol 2
  !&E---------------------------------------------------------------------

    !! * Arguments
    integer, INTENT(IN) :: n
    real(kind=rsh),INTENT(IN)   :: x
    real(kind=rsh),INTENT(OUT)  :: y
    real(kind=rsh),dimension(0:n-1), INTENT(IN) :: xa,y2a,ya
    
    !! * Local declarations
    INTEGER                 :: IERR_MPI
    integer                     :: k,khi,klo
    real(kind=rsh)              :: a,b,h
    
    !!----------------------------------------------------------------------
    !! * Executable statment
    
    klo=0
    khi=n-1
    CALL ksupkinf(x,xa,n,khi,klo,4)

    if ((khi>kmax) .or. (klo>kmax) .or. (khi<0) .or. (klo<0)) then
       write(*,*) 'mauvais sigma fonction splint'
       write(*,*) khi,klo,kmax
       write(*,*) x
       write(*,*) xa
       CALL_MPI MPI_FINALIZE(IERR_MPI)
       stop
    end if
    
    h=xa(khi)-xa(klo)
    if (h==0.0) then
       write(*,*) 'bad xa input in splint'
       read(*,*)
    end if
    a=(xa(khi)-x)/h
    b=(x-xa(klo))/h
    y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.0
    !
  end subroutine splint

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


#ifdef key_siggen

!!======================================================================
   SUBROUTINE f_lag_sigz_uv(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 (relative to u and v location
   !&E              for vertical interpolation at particle location
   !&E
   !&E ** Description :
   !&E
   !&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)) THEN
      dcwsdsp=  (Csu_sig(k+1)-Csusig)/(sig(k+1)-spos)
   ENDIF
   IF (spos /= sig(k)) THEN
      dcwsdsm=  (Csusig-Csu_sig(k))/(spos-sig(k))
   ENDIF

   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_uv


 !!======================================================================
   SUBROUTINE h_int_siggen(xeint,h0int,hc_sigint,xe,px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE h_int  ***
   !&E
   !&E ** Purpose : calcul par interpolation lineaire de la hauteur d eau
   !&E              a la position (xpos_s,ypos_s,spos_s)
   !&E              (also return xe and h0 interpolated)
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !          (P. Lazure, F. Dumas)  Original code
   !&E       !  08-2008 (F. Dumas ?)
   !&E       !  12/2010 (M. huret)
   !&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars2d, ONLY : h0_g
   USE comsiggen, ONLY : hc_sig

   !! * Arguments
   REAL(KIND=rsh),                                 INTENT( in )  :: px,py
   INTEGER,                                        INTENT( in )  :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh), DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2), INTENT( in )  :: xe
   REAL(KIND=rsh),                   INTENT( out )                 :: xeint,h0int,hc_sigint

   !! * Local declarations


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

   ! spatial interpolation

   xeint =px*(1.0-py)*xe(i2,j1)*rh21+ px*py*xe(i2,j2)*rh22 &
         +(1.0-px)*py*xe(i1,j2)*rh12+(1.0-px)*(1.0-py)*xe(i1,j1)*rh11
   h0int =px*(1.0-py)*h0_g(i2,j1)*rh21+ px*py*h0_g(i2,j2)*rh22 &
         +(1.0-px)*py*h0_g(i1,j2)*rh12+(1.0-px)*(1.0-py)*h0_g(i1,j1)*rh11
   hc_sigint =px*(1.0-py)*hc_sig(i2,j1)*rh21+ px*py*hc_sig(i2,j2)*rh22 &
         +(1.0-px)*py*hc_sig(i1,j2)*rh12+(1.0-px)*(1.0-py)*hc_sig(i1,j1)*rh11
 

   END SUBROUTINE h_int_siggen

  !!======================================================================
   FUNCTION hc_sigint(px,py,i1,i2,j1,j2,rh11,rh21,rh12,rh22)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE hc_sigint  ***
   !&E
   !&E ** Purpose : Linear interpolation of hc_sig at the particle location 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       ! 12/2010 (M. huret)
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comsiggen, ONLY : hc_sig

   !! * Arguments
   REAL(KIND=rsh),                   INTENT( in )  :: px,py
   INTEGER,                          INTENT( in )  :: i1,i2,j1,j2,rh11,rh12,rh21,rh22
   REAL(KIND=rsh)                                  :: hc_sigint

   !! * Local declarations


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

   ! spatial interpolation

   hc_sigint =px*(1.0-py)*hc_sig(i2,j1)*rh21+ px*py*hc_sig(i2,j2)*rh22 &
         +(1.0-px)*py*hc_sig(i1,j2)*rh12+(1.0-px)*(1.0-py)*hc_sig(i1,j1)*rh11

   END FUNCTION hc_sigint

!!======================================================================
   SUBROUTINE f_lag_sigz_w(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,dcusds,Csw_sig
   USE comvars3d,   ONLY : sigw

   !! * 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)      :: Cswsig, dcusdsp,dcusdsm
   
   !!----------------------------------------------------------------------
   !! * Executable part

   Cswsig=(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)

   dcusdsp=0.0_rsh
   dcusdsm=0.0_rsh

   IF (spos /= sigw(k)) THEN
      dcusdsp=  (Csw_sig(k)-Cswsig)/(sigw(k)-spos)
   ENDIF
   IF (spos /= sigw(k-1)) THEN
      dcusdsm=  (Cswsig-Csw_sig(k-1))/(spos-sigw(k-1))
   ENDIF

   f_lag_hzp=xe_lag+hc_sig_lag+(h0_lag-hc_sig_lag)* dcusdsp
   f_lag_hzm=xe_lag+hc_sig_lag+(h0_lag-hc_sig_lag)* dcusdsm
   f_lag_hz=xe_lag+hc_sig_lag+(h0_lag-hc_sig_lag)* dcusds(k)

  END SUBROUTINE f_lag_sigz_w


!!======================================================================
   SUBROUTINE siggentoz(zpos,spos,xe_lag,h0_lag,hc_sig_lag)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE siggen_int  ***
   !&E
   !&E ** Purpose : Calculate z from sigma (for siggen case)
   !&E              
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d, ztosiggen
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used

   USE comsiggen, ONLY : b_sig,theta_sig

   !! * Arguments
   REAL(KIND=rsh),                                 INTENT( in )  :: spos,xe_lag,hc_sig_lag,h0_lag
   REAL(KIND=rsh),                                 INTENT( out ) :: zpos

   !! * Local declarations
   REAL(KIND=rsh)      :: Cssig

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

   Cssig=(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)
   zpos=xe_lag*(1.0_rsh+spos)+hc_sig_lag*spos+(h0_lag-hc_sig_lag)*Cssig

  END SUBROUTINE siggentoz

!!======================================================================
   SUBROUTINE ztosiggen(zpos_s,spos,xelag,h0lag,hc_siglag)

   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ztosiggen  ***
   !&E
   !&E ** Purpose : Calcul du sigma a partir du z (en sigma generalises)
   !&E              
   !&E ** Description :
   !&E
   !&E ** Called by : traj_3d
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference :
   !&E
   !&E ** History : M. Huret (to be improved, only linear interpolation
   !&E                should use a correction factor considering Ds/Dz
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE comvars3d,  ONLY : sigw

   !! * Arguments
   REAL(KIND=rsh),                                 INTENT( in )  :: zpos_s,xelag,h0lag,hc_siglag
   REAL(KIND=rsh),                                 INTENT( out ) :: spos

   !! * Local declarations
   INTEGER             :: kw,kwm,l
   REAL(KIND=rsh), DIMENSION(0:kmax)           :: zpos
   
   !!----------------------------------------------------------------------
   !! * Executable part

   DO l=0,kmax
      CALL siggentoz(zpos(l),sigw(l),xelag,h0lag,hc_siglag)
   END DO
   kw=kmax
   kwm=0
   CALL ksupkinf(zpos_s,zpos,kmax+1,kw,kwm,1)
   spos=((zpos_s-zpos(kwm))*sigw(kw)+(zpos(kw)-zpos_s)*sigw(kwm))/(zpos(kw)-zpos(kwm))  !!Approximation


  END SUBROUTINE ztosiggen


#endif

#endif
 END MODULE
