  SUBROUTINE ts_diffh_rho(sal,temp,xev)

#ifdef key_ts_diffh_rho
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE ts_diffh_rho  ***
   !&E
   !&E ** Purpose : compute horizontal diffusion of temperature and salinity
   !&E              along isopycnal surfaces (rho=constant)
   !&E
   !&E ** Description : For very low water heights ( < hm), 
   !&E                  no horizontal diffusion is applied.
   !&E                  - The cycling indices j1,j2,j3 correspond to j-1,j,j+1
   !&E                  - The "gradient" function is a CPP function that makes the code clearer
   !&E                    and is defined at the begining of the routine.
   !&E
   !&E ** Called by : tssub_eq_wcen.h, tssub_eq_wquick.h
   !&E
   !&E ** External calls :
   !&E
   !&E ** Reference : Marchesiello, 1995, PhD
   !&E
   !&E ** History :
   !&E       !  2010-10  (T. Duhaut) Original code
   !&E
   !&E ** Remarque :
   !&E       ! Attention diffusion non nulle pour les petites profondeurs < hminkxky
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters
   USE comvars2d, ONLY : ig,id,h0,surf,hex,hx,kx,h0fond,hm,hey,hy,   &
                         dy_u,dy_v,ky,dx_u,dx_v,jbound_omp,halfdt
   USE comvars3d, ONLY : sig,dsigu,rhoref,dsigw, kz
   USE comvarp3d, ONLY : boz
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzex,hzey,hc_sig,csu_sig,dcwsds
   USE toolgeom,  ONLY : f_hzw
#endif
   USE comfluxdo, ONLY : fdifsal
   USE cometemp,  ONLY : fdiftemp

   IMPLICIT NONE

   !! * Arguments
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: sal,temp
   REAL(KIND=rsh), DIMENSION(limin:limax,ljmin:ljmax), INTENT(in)  :: xev

   !! * Local declarations
   INTEGER                                        :: j1,j2,j3,i,j,k
   INTEGER                                        :: lompjmin,lompjmax
   INTEGER                                        :: OMP_GET_THREAD_NUM
   REAL(KIND=rsh)                                 :: cffx,cffy,cff1,cff2,cff3,cff4,cff5,cff58,cff59
   REAL(KIND=rsh)                                 :: drho,EPS,hz_w,dZdr, rho_s_g, pp
   REAL(KIND=rsh),DIMENSION(0:kmax,imin:imax)     :: FXt,FXs
   REAL(KIND=rsh),DIMENSION(0:kmax,imin:imax,0:2) :: dTdr,dSdr
   REAL(KIND=rsh),DIMENSION(0:kmax,imin:imax)     :: dTdx,dSdx
   REAL(KIND=rsh),DIMENSION(0:kmax,imin:imax)     :: dRdx,FSt,FSs
   REAL(KIND=rsh),DIMENSION(0:kmax,imin:imax,0:2) :: FEt,FEs,dTde,dSde,dRde

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

   !   definition de la fonction inline gradient  
   ! macro de gradient entre deux valeurs ecartee d une distance de 1/delta
# define gradient(unsdelta,val1,val2) (unsdelta)*((val2)-(val1))

   ! define EPS value
   EPS = 0.5_rsh
   rho_s_g=-rhoref/grav
   pp = 0.0_rsh
! cthom test hex,hzex
!hex = 4000.0
!hey = 4000.0

   ! define rotating indexes j1,j2,j3
   j1=1
   j2=2
   j3=3-j1-j2

   ! LOOP ON J
   ! modify loop bounds for OMP parallelization
!lompjmin = MAX(jmin,OMP_GET_THREAD_NUM()*116-2)
!lompjmax = MIN(jmax-1,(OMP_GET_THREAD_NUM()+1)*116-1)
!$ lompjmin = MAX(jmin,jbound_omp(OMP_GET_THREAD_NUM())-2)
!$ lompjmax = MIN(jmax-1,jbound_omp(OMP_GET_THREAD_NUM()+1)-1)
   lompjmin = MAX(jmin,ljmin-2)
   lompjmax = ljmaxm1

   DO j=lompjmin,lompjmax

!print*,'THREAD NUM : ',OMP_GET_THREAD_NUM(), ' loop on j = ',lompjmin,lompjmax,j

     j1=j2
     j2=j3
     j3=3-j1-j2

     ! fields out of domain set to 0
     FXt(:,imin:MIN(ig(j),ig(j+1)))=0.0_rsh
     FXs(:,imin:MIN(ig(j),ig(j+1)))=0.0_rsh
     FXt(:,MAX(id(j),id(j+1))-1:imax)=0.0_rsh
     FXs(:,MAX(id(j),id(j+1))-1:imax)=0.0_rsh
     FEt(:,imin:MIN(ig(j),ig(j+1)),j2)=0.0_rsh
     FEs(:,imin:MIN(ig(j),ig(j+1)),j2)=0.0_rsh
     FEt(:,MAX(id(j),id(j+1))-1:imax,j2)=0.0_rsh
     FEs(:,MAX(id(j),id(j+1))-1:imax,j2)=0.0_rsh

     ! LOOP ON I 
     ! correspond a : do i=ig(j),id(j)-1 mais plus large
     ! pour n exclure aucun cas particulier de bathymetrie
     DO i=MIN(ig(j),ig(j+1)),MAX(id(j),id(j+1))-1 ! calcul des d.dx et d.dz

       ! flux entre i et i+1
       IF ( hex(i,j) > hm .AND. (hx(i,j).NE.( -valmanq-h0fond))) THEN
         cffx=1.0_rsh/dx_u(i,j)
       ELSE
         cffx=0.0_rsh
       END IF
       ! flux entre j et j+1
       IF ( hey(i,j) > hm .AND. (hy(i,j).NE.( -valmanq-h0fond))) THEN
         cffy=1.0_rsh/dy_v(i,j) 
       ELSE
         cffy=0.0_rsh
       END IF
            
       !---------------------------------------
       ! 1- calcul des gradients de temperature et salinite
       !---------------------------------------

       DO k=1,kmax-1

         ! a- derivees suivant x a la latitude j
         !    gradient entre i et i+1
         !--------------------------------------------------

         dRdx(k,i)=rho_s_g*gradient(cffx,boz(k,i,j),boz(k,i+1,j))
         dTdx(k,i)=gradient(cffx,temp(k,i,j),temp(k,i+1,j))
         dSdx(k,i)=gradient(cffx,sal(k,i,j),sal(k,i+1,j))


         ! b- derivees suivant y a la latitude j+1/2 (j2)
         !    gradient entre j et j+1
         !--------------------------------------------------
         dRde(k,i,j2)=(-rhoref/grav)*gradient(cffy,boz(k,i,j),boz(k,i,j+1)) 
         dTde(k,i,j2)=gradient(cffy,temp(k,i,j),temp(k,i,j+1))
         dSde(k,i,j2)=gradient(cffy,sal (k,i,j),sal (k,i,j+1))
!         dRde(k,i,j2)=0.0_rsh
!         dTde(k,i,j2)=0.0_rsh
!         dSde(k,i,j2)=0.0_rsh


         ! c- derivees suivant rho la latitude j+1+1/2 (j3)
         !    gradient entre j+2 et j+1
         !--------------------------------------------------
         drho = -1.0_rsh/MAX(-rhoref/grav*(boz(k,i,j+1)-boz(k+1,i,j+1)),EPS)
         dTdr(k,i,j3)=gradient(drho, temp(k,i,j+1),temp(k+1,i,j+1))
         dSdr(k,i,j3)=gradient(drho, sal(k,i,j+1) ,sal (k+1,i,j+1))

!         drho = -1.0_rsh/MAX(rho_s_g*(boz(k,i,j)-boz(k+1,i,j)),EPS)
!         dTdr(k,i,j2)=gradient(drho, temp(k,i,j),temp(k+1,i,j))
!         dSdr(k,i,j2)=gradient(drho, sal(k,i,j) ,sal (k+1,i,j))

       END DO ! ends loop on k                

       ! conditions aux limites: fond&surface
            
       ! surface

       dRdx(kmax,i)=rho_s_g*gradient(cffx,boz(k,i,j),boz(k,i+1,j))
       dRde(kmax,i,j2)=(-rhoref/grav)*gradient(cffy,boz(k,i,j),boz(k,i,j+1))
   !    dRde(kmax,i,j2)=(-rhoref/grav)*gradient(cffy,boz(k,i,j),boz(k,i,j+1))
        

       dTdx(kmax,i)=gradient(cffx,temp(kmax,i,j),temp(kmax,i+1,j))
       dTde(kmax,i,j2)=gradient(cffy,temp(kmax,i,j),temp(kmax,i,j+1))
   !    dTde(kmax,i,j2)=gradient(cffy,temp(kmax,i,j),temp(kmax,i,j+1))
            
       dSdx(kmax,i)=gradient(cffx,sal(kmax,i,j),sal(kmax,i+1,j))
       dSde(kmax,i,j2)=gradient(cffy,sal (kmax,i,j),sal (kmax,i,j+1))
   !    dSde(kmax,i,j2)=gradient(cffy,sal (kmax,i,j),sal (kmax,i,j+1))
            

       dTdr(kmax,i,j3)=0.0_rsh
       dSdr(kmax,i,j3)=0.0_rsh

       FSt(kmax,i)=0.0_rsh
       FSs(kmax,i)=0.0_rsh
          
       ! bottom
       dTdr(0,i,j3)=0.0_rsh
       dSdr(0,i,j3)=0.0_rsh
       FSt(0,i)=0.0_rsh
       FSs(0,i)=0.0_rsh

     END DO ! ends loop on i

     IF (id(j)==imax) THEN

       DO k=1,kmax-1
         drho = -1.0_rsh/MAX(-rhoref/grav*(boz(k,i,j+1)-boz(k+1,i,j+1)),EPS)
         dTdr(k,imax,j3)=gradient(drho, temp(k,imax,j+1),temp(k+1,imax,j+1))
         dSdr(k,imax,j3)=gradient(drho, sal (k,imax,j+1),sal (k+1,imax,j+1))
    !     drho = -1.0_rsh/MAX(-rhoref/grav*(boz(k,i,j+1)-boz(k+1,i,j+1)),EPS)
    !     dTdr(k,imax,j3)=gradient(drho, temp(k,imax,j+1),temp(k+1,imax,j+1))
    !     dSdr(k,imax,j3)=gradient(drho, sal (k,imax,j+1),sal (k+1,imax,j+1))

         drho = -1.0_rsh/MAX(rho_s_g*(boz(k,i,j)-boz(k+1,i,j)),EPS)
         dTdr(k,imax,j2)=gradient(drho, temp(k,imax,j),temp(k+1,imax,j))
         dSdr(k,imax,j2)=gradient(drho, sal (k,imax,j),sal (k+1,imax,j))
       END DO ! ends loop on k

       dTdr(kmax,imax,j3)=0.0_rsh
       dSdr(kmax,imax,j3)=0.0_rsh
       dTdr(kmax,imax,j2)=0.0_rsh
       dSdr(kmax,imax,j2)=0.0_rsh
       FSt(kmax,imax)=0.0_rsh
       FSs(kmax,imax)=0.0_rsh

       dTdr(0,imax,j3)=0.0_rsh
       dSdr(0,imax,j3)=0.0_rsh
       dTdr(0,imax,j2)=0.0_rsh
       dSdr(0,imax,j2)=0.0_rsh
       FSt(0,imax)=0.0_rsh
       FSs(0,imax)=0.0_rsh

     END IF

     IF (j>(lompjmin)) THEN ! calcul de Fe entre j et j+1 possible, car derivees selon z definies en j2 (j) et j3 (j+1)
       DO i=ig(j),id(j)-1
         DO k=1,kmax

           !----------------------------------------
           ! 2- calcul de FX temperature et salinite 
           !----------------------------------------
                      ![FX]=m.C.S^-1                                              
                      ! modif de ROMS a MARS:
                      ! Kx est constant donc on ne calcule pas de Kx moyen
                      ! Hz(i,j,k)+Hz(i-1,j,k)=>disparait dans MARS, introduit plus tard
                      ! dy =>disparait dans MARS, introduit plus tard

                      ! suivant la pente du sol, les FX et FE ne sont pas calcules
                      ! avec les memes cellules pour diminuer l erreur engendree
                      ! par la pente

         IF ( hex(i,j) > hm .AND. (hx(i,j).NE.( -valmanq-h0fond))) THEN

           ! entre i et i+1
           ! temperature
           FXt(k,i)=kx*(dTdx(k,i)-                                                          &
                    0.5_rsh*(MAX(dRdx(k,i),0.0_rsh)*(dTdr(k-1,i,j2) + dTdr(k  ,i+1,j2))+    & ! pente dZdx<0
                             MIN(dRdx(k,i),0.0_rsh)*(dTdr(k  ,i,j2) + dTdr(k-1,i+1,j2)) ))    ! pente dZdx>0          

           ! salinite
           FXs(k,i)=kx*(dSdx(k,i)-                                                          &
                    0.5_rsh*(MAX(dRdx(k,i),0.0_rsh)*(dSdr(k-1,i,j2) + dSdr(k  ,i+1,j2))+    &
                             MIN(dRdx(k,i),0.0_rsh)*(dSdr(k  ,i,j2) + dSdr(k-1,i+1,j2)) )) 

         ELSE
           FXt(k,i)=0.0_rsh
           FXs(k,i)=0.0_rsh
         ENDIF

           !----------------------------------------------
           ! 3- calcul de FE temperature et salinite en j2 
           !----------------------------------------------

           ! Entre j et j+1  
           ! temperature
           FEt(k,i,j2)=ky*(dTde(k,i,j2)-                                                       &
                       0.5_rsh*(MAX(dRde(k,i,j2),0.0_rsh)*(dTdr(k-1,i,j2) + dTdr(k  ,i,j3))+   & ! pente dZdx<0
                                MIN(dRde(k,i,j2),0.0_rsh)*(dTdr(k  ,i,j2) + dTdr(k-1,i,j3))  ))  ! pente dZdx>0
     !      FEt(k,i,j2)=ky*(dTde(k,i,j2)-                                                       &
     !                  0.5_rsh*(MAX(dRde(k,i,j2),0.0_rsh)*(dTdr(k-1,i,j2) + dTdr(k  ,i,j3))+   & ! pente dZdx<0
     !                           MIN(dRde(k,i,j2),0.0_rsh)*(dTdr(k  ,i,j2) + dTdr(k-1,i,j3))  ))  ! pente dZdx>0

                  
           ! salinite
           FEs(k,i,j2)=ky*(dSde(k,i,j2)-                                                       &
                       0.5_rsh*(MIN(dRde(k,i,j2),0.0_rsh)*(dSdr(k-1,i,j2) + dSdr(k  ,i,j3))+   &
                                MAX(dRde(k,i,j2),0.0_rsh)*(dSdr(k  ,i,j2) + dSdr(k-1,i,j3))  ))
     !      FEs(k,i,j2)=ky*(dSde(k,i,j2)-                                                       &
     !                  0.5_rsh*(MIN(dRde(k,i,j2),0.0_rsh)*(dSdr(k-1,i,j2) + dSdr(k  ,i,j3))+   &
     !                           MAX(dRde(k,i,j2),0.0_rsh)*(dSdr(k  ,i,j2) + dSdr(k-1,i,j3))  ))


         END DO ! ends loop on k 
       END DO ! ends loop on i
     END IF


     IF (j>(lompjmin+1)) THEN              
       ! calcul de Fs en j possible, car derivees selon y definies en j1 (j-1) et j2 (j)
       ! possible en jminp1 aussi mais inutile
       DO i=ig(j)+1,id(j)-1
         DO k=1,kmax-1 ! borne indices k modifiee car FS9kmax)=0

           !--------------------------------------------
           ! 4- calcul de FS temperature et salinite en j2
           !--------------------------------------------

#if defined key_siggen || defined key_gencoord
           hz_w = f_hzw(h0(i,j),xev(i,j),k,i,j)
#else
           hz_w = h0(i,j)+xev(i,j)
#endif

           dZdr=(dsigw(k)*hz_w)*(-1.0_rsh)/MAX(rho_s_g*(boz(k,i,j)-boz(k+1,i,j)),EPS)

           !------------- suivant x
           cff1=MAX(dRdx(k  ,i-1  ),0.0_rsh)
           cff2=MAX(dRdx(k+1,i),0.0_rsh)
           cff3=MIN(dRdx(k+1,i-1),0.0_rsh)
           cff4=MIN(dRdx(k  ,i  ),0.0_rsh)
!           FSt(k,i)=0.5_rsh*kx*dZdr*                           &
!                    (cff1*(cff1*dTdr(k,i,j2)-dTdx(k ,i-1))+    &
!                     cff2*(cff2*dTdr(k,i,j2)-dTdx(k+1,i))+     &
!                     cff3*(cff3*dTdr(k,i,j2)-dTdx(k+1,i-1))+   &
!                     cff4*(cff4*dTdr(k,i,j2)-dTdx(k  ,i)))

           FSt(k,i)=0.5_rsh*kx*dZdr*(-cff1*dTdx(k ,i-1)         &
                   -cff2*dTdx(k+1,i)-cff3*dTdx(k+1,i-1)-cff4*dTdx(k  ,i) )

!           FSs(k,i)=0.5_rsh*kx*dZdr*                           &
!                    (cff1*(cff1*dSdr(k,i,j2)-dSdx(k,i-1))+     &
!                     cff2*(cff2*dSdr(k,i,j2)-dSdx(k+1,i))+     &
!                     cff3*(cff3*dSdr(k,i,j2)-dSdx(k+1,i-1))+   &
!                     cff4*(cff4*dSdr(k,i,j2)-dSdx(k,i)))

           FSs(k,i)=0.5_rsh*kx*dZdr*(-cff1*dSdx(k ,i-1)         &
                   -cff2*dSdx(k+1,i)-cff3*dSdx(k+1,i-1)-cff4*dSdx(k  ,i) )

           cff58 = 0.5_rsh*kx*dZdr*(cff1**2+cff2**2+cff3**2+cff4**2)

           FSt(k,i) = FSt(k,i) + (1-pp)*cff58*dTdr(k,i,j2)

           FSs(k,i) = FSs(k,i) + (1-pp)*cff58*dSdr(k,i,j2)
           kz(k,i,j) = kz(k,i,j) + pp * cff58 * dZdr * dsigu(k)/dsigw(k)

! Thomas impliciter les termes FS selon y Merci
                  
           !-------------  suivant y
     !      cff1=MAX(dRde(k  ,i,j1),0.0_rsh)
     !      cff2=MAX(dRde(k+1,i,j2),0.0_rsh)
     !      cff3=MIN(dRde(k+1,i,j1),0.0_rsh)
     !      cff4=MIN(dRde(k  ,i,j2),0.0_rsh)
     !      FSt(k,i)=FSt(k,i)+0.5_rsh*ky*dZdr*                  &
     !               (cff1*(cff1*dTdr(k,i,j2)-dTde(k ,i,j1))+   &
     !                cff2*(cff2*dTdr(k,i,j2)-dTde(k+1,i,j2))+  &
     !                cff3*(cff3*dTdr(k,i,j2)-dTde(k+1,i,j1))+  &
     !                cff4*(cff4*dTdr(k,i,j2)-dTde(k  ,i,j2)))
     !      FSs(k,i)=FSs(k,i)+0.5_rsh*ky*dZdr*                  &
     !               (cff1*(cff1*dSdr(k,i,j2)-dSde(k,i,j1))+    &
     !                cff2*(cff2*dSdr(k,i,j2)-dSde(k+1,i,j2))+  &
     !                cff3*(cff3*dSdr(k,i,j2)-dSde(k+1,i,j1))+  &
     !                cff4*(cff4*dSdr(k,i,j2)-dSde(k,i,j2)))
                  
         END DO ! ends loop on k
                            

         !--------------------------------------------
         ! 5- calcul de fdifsal et fdiftemp
         !--------------------------------------------

         DO k=1,kmax 
                      ! [FX]=m.C.s^-1
                      ! [FS]=m.C.S^-1
                      ! attention: calcul valable que si dy constant
                      ! si dy non constant, alors modifier dy_u(i,j) suivant les bords de la cellule        
                      ! on calcule d./dx entre i-1 et i: avant c etait fait dans tssub_eq_wcen.h !
                      ! sauf que le 1/surf(i,j) est encore dans tssub_eq_wcen.h
                      ! verification de h0+xe (i,j) > hm dans tssub_eq_wcen.h

#if defined key_siggen || defined key_gencoord
           fdifsal(k,i,j)= -halfdt*(                                 &
                          (FXs(k,i  )*hzex(k,i  ,j)*dy_u(i,j)-       &
                           FXs(k,i-1)*hzex(k,i-1,j)*dy_u(i-1,j)  )+  &
     !                     (FEs(k,i,j2)*hzey(k,i  ,j)*dx_v(i,j)-      &
     !                      FEs(k,i,j1)*hzey(k,i,j-1)*dx_v(i,j-1)  )+ &
                          (FSs(k,i)-FSs(k-1,i))*surf(i,j)   )
           fdiftemp(k,i,j)= -halfdt*(                                &
                          (FXt(k,i  )*hzex(k,i  ,j)*dy_u(i,j)-       &
                           FXt(k,i-1)*hzex(k,i-1,j)*dy_u(i-1,j)  )+  & 
     !                     (FEt(k,i,j2  )*hzey(k,i  ,j)*dx_v(i,j)-    &
     !                      FEt(k,i,j1)*hzey(k,i,j-1)*dx_v(i,j-1)  )+ & 
                          (FSt(k,i)-FSt(k-1,i))*surf(i,j)   )
                          if ( k == 50 .and. i==200 ) then
                            print*,'FXt(k,i,j)=',FXt(k,i),FXt(k,i-1)
                            print*,'FET(k,i,j)=',FEt(k,i,j2),FEt(k,i,j1)
                            print*,'FSt(k,i,j)=',FSt(k,i),FSt(k-1,i)
                            print*,'fdiftemp(k,i,j)=',fdiftemp(k,i,j)
                          endif
#else
           fdifsal(k,i,j)= -halfdt*(                              &
                          (FXs(k,i  )*hex(i  ,j)*dy_u(i,j)-       &
                           FXs(k,i-1)*hex(i-1,j)*dy_u(i-1,j)  )+  &
     !                     (FEs(k,i,j2)*hey(i  ,j)*dx_v(i,j)-      &
     !                      FEs(k,i,j1)*hey(i,j-1)*dx_v(i,j-1)  )+ &
                          (FSs(k,i)-FSs(k-1,i))*surf(i,j)   )
           fdiftemp(k,i,j)= -halfdt*(                             &
                          (FXt(k,i  )*hex(i  ,j)*dy_u(i,j)-       &
                           FXt(k,i-1)*hex(i-1,j)*dy_u(i-1,j)  )+  & 
     !                     (FEt(k,i,j2  )*hey(i  ,j)*dx_v(i,j)-    &
     !                      FEt(k,i,j1)*hey(i,j-1)*dx_v(i,j-1)  )+ & 
                          (FSt(k,i)-FSt(k-1,i))*surf(i,j)   )
#endif

           fdifsal(k,i,j)= 0.0
 
         END DO ! ends loop on k

       END DO ! ends loop on i
               
     END IF ! ends test if j>jmin+1
   END DO ! ends loop on j

#else

   USE parameters, ONLY : rsh,liminm1,limaxp2,ljminm1,ljmaxp2,limin,limax,ljmin,ljmax,kmax
   REAL(KIND=rsh), DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2), INTENT(in) :: sal,temp
   REAL(KIND=rsh), DIMENSION(limin:limax,ljmin:ljmax), INTENT(in)  :: xev

#endif

  END SUBROUTINE ts_diffh_rho
