!  SUBROUTINE tssub_eq(xe2,xe3,sal,temp        &
!#ifdef 1
!                                      ,cv_wat &
!#endif
!                                             )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_eq  ***
!&E
!&E ** Purpose : framework for the subroutine of tracers transport computation
!&E
!&E ** Description : see tssub_eq_wquick.h or tssub_eq_wcen.h
!&E                  for different possibilities of subroutine body
!&E
!&E ** Called by : stepu, stepv
!&E
!&E ** External calls : tssub_eq_hm
!&E
!&E ** Used ij-arrays : xe2, xe3, u, v
!&E
!&E ** Modified variables : sal, temp, cv_wat
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2004-10-06
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Arguments

!! * Local declarations

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



  SUBROUTINE tssub_eq(xe2,xe3,sal,temp        &

                                      ,cv_wat &

                                             )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_eq_wcompact  ***
!&E
!&E ** Purpose : other possibility of subroutine body for the computation
!&E              of transport of salinity, temperature and tracers without fall velocity.
!&E
!&E ** Description : An Euler-Compact scheme is used for depths higher than hm.
!&E                  For water heights lower than hm, a special computation
!&E                  is performed in the subroutine tssub_eq_hm.
!&E                  This subroutine body replaces the default body
!&E                  when the key 1 is actived.
!&E                  remark : substances advection is identical whatever is the variable type
!&E                          (dissolved or with settling velocity)
!&E                  The particles fall velocity is taken into account in verti_bio.
!&E                  WARNING : valid for very low fall velocities !!!
!&E
!&E ** Called by : stepu,stepv
!&E
!&E ** External calls : tssub_advdifx,tssub_advdify,ts_advz_compact,sub_advz_compact
!&E
!&E ** Used ij-arrays : sal,temp,cv_wat,h0,xe2,xe3,kz,flx_w2s,
!&E                     fsalx,fsaly,fdifsal,phisal,
!&E                     ftempx,ftempy,fdiftemp,phitemp,
!&E                     fx,fy,phicon,dy_u,dx_v,surf
!&E                     ftempz,fsalz,fz
!&E
!&E ** Modified variables : sal,temp,cv_wat
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-05-14  (V. Garnier) adaptation of quickupwind for euler compact scheme
!&E       !  2010-06     (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
!&E       !  2010-06  (F. Cayocca) coupling with the module of sedimentology
!&E       !  2010-10 (T. Duhaut) introduction of diffusion along isopycnes
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2016-04  (B.Thouvenin) applied to all substances (dissolved and particulate, include sand)
!&E       !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,   ONLY : ig,id,h0,hm,dy_u,dx_v,surf,l_tssubadv_over_smalldepth,halfdt,fwet,fwetp
   USE comvars3d,   ONLY : kz,dsigw,dsigu
   USE comfluvert,  ONLY : ftempz,fsalz

   USE comsubstance,ONLY : fx,fy,fz,nv_adv,nv_grav,nvp,phicon,f90,typdiss,l_subs2D



   USE comfluxdo,   ONLY : fsalx,fsaly
   USE cometemp,    ONLY : ftempx,ftempy,phitemp

   USE cometemp,    ONLY : phisal

!!#ifdef key_sedim_MUSTANG_mudsand
!!   USE comsubstance,ONLY : flx_s2w,flx_w2s
!!#else
   USE comMUSTANG,ONLY : flx_s2w
!!#endif


   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax),                    INTENT(in)    :: xe2
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),            INTENT(in)    :: xe3
!cval xe2 not ready for tsdiffh
   REAL(kind=rsh),DIMENSION(kmax,liminm1:limaxp2,ljminm1:ljmaxp2),       INTENT(inout) :: sal,temp

   REAL(kind=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_wat


!! * Local declarations
   INTEGER                                   :: n,i,j,k,l
   REAL(kind=rsh)                            :: dtm,tfleuve,d3,d2,unsd3

   REAL(kind=rlg),DIMENSION(kmax,nb_var)     :: y

   REAL(kind=rlg)                            :: xk,yk,b2d
   REAL(kind=rlg),DIMENSION(kmax)            :: a,b,c,ysal,u0,v0,p0,q0,ytemp

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

   !PRINT*, 'enter TSSUB_EQ_WCOMPACT.H routine'



!$OMP parallel default(none)                               &
!$OMP& PRIVATE(d2,d3,unsd3,i,j,k,l,xk,yk,                  &
!$OMP&         dtm,p0,q0,u0,v0,a,b,c,b2d,                  &

!$OMP&                               y,                    &

!$OMP&                                 ysal,ytemp)         &
!$OMP& shared(sal,temp,xe3,phisal,phitemp,hm,ig,id,h0      &
!$OMP& ,kz,dsigw,dsigu,dy_u                                &
!$OMP& ,limin,limax,ljmin,ljmax,ljminp2,ljmaxm1            &
!$OMP& ,liminm1,limaxp2,ljminm1,ljmaxp2,fwet,fwetp           &
!$OMP& ,fsalx,dx_v,fsaly,fsalz,surf,ftempx,ftempy,ftempz   &


!$OMP& ,cv_wat,nvp,nv_adv,nv_grav,fx,fy,fz,phicon,f90,typdiss,l_subs2D   &


!$OMP& ,flx_s2w                                            &


!$OMP& ,l_tssubadv_over_smalldepth,xe2,halfdt)

   dtm=halfdt

! calcul des flux horizontaux et verticaux (advectifs et diffusifs)
!------------------------------------------------------------------

! horizontal diffusion
! --------------------




! advection (+ horizontal diffusion if not key_ts_diffh_geo or key_ts_diffh_rho)
! ------------------------------------------------------------------------------






   CALL tssub_advdifx(sal,temp,cv_wat)
   CALL tssub_advdify(sal,temp,cv_wat)


   CALL ts_advz_compact(xe3,sal,temp)


   CALL sub_advz_compact(xe3,cv_wat)






!$OMP DO schedule(runtime)

   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       d3=h0(i,j)+xe3(i,j)
       d2=h0(i,j)+xe2(i,j)
       IF(d3 > hm) THEN
         unsd3=1.0_rsh/d3

! matrice tridiagonale du modele vertical
!----------------------------------------

!-------------------------
         k=1
!-------------------------
         c(k)=-dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
         b(k)=d3*fwet(i,j)+dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
!!#if defined 1 && defined key_sedim_mudsand
!!! pour mixsed :  plus besoin car pour gestion des variables dissoutes maintenant
!!! le flux w2s a deja ete ajoute au flux s2w car traite en explicite
!!         b(k)=b(k)+dtm*flx_w2s(-1,i,j)/dsigu(k) ! car flux identique pour toute
!!                                                ! variable dissoute et calcule
!!                                                ! pour temp(-1) et sal (0)
!!                                                ! flx_w2s(-1,i,j)=flx_w2s( 0,i,j)
!!#endif

         a(k)=0.0_rlg





         ysal(k)=d2*fwetp(i,j)*sal(k,i,j)-(dy_u(i,j)*fsalx(k,i,j)-dy_u(i-1,j)*fsalx(k,i-1,j) &
                +dx_v(i,j)*fsaly(k,i,j)-dx_v(i,j-1)*fsaly(k,i,j-1))/surf(i,j)     &
                +phisal(k,i,j)/dsigu(k)
         ysal(k)=ysal(k)-(fsalz(k,i,j)-fsalz(k-1,i,j))/dsigu(k)

         ysal(k)=ysal(k)+dtm*flx_s2w(0,i,j)/dsigu(k)


         ytemp(k)=d2*fwetp(i,j)*temp(k,i,j)-(dy_u(i,j)*ftempx(k,i,j)-dy_u(i-1,j)*ftempx(k,i-1,j) &
                 +dx_v(i,j)*ftempy(k,i,j)-dx_v(i,j-1)*ftempy(k,i,j-1))/surf(i,j)      &
                 +phitemp(k,i,j)/dsigu(k)
         ytemp(k)=ytemp(k)-(ftempz(k,i,j)-ftempz(k-1,i,j))/dsigu(k)

         ytemp(k)=ytemp(k)+dtm*flx_s2w(-1,i,j)/dsigu(k)  !! ATTENTION PB dx dy et dt !!!!!




! key_offline


         q0(k)=-c(k)/b(k)        ! factorisation LU


         u0(k)=ysal(k)/b(k)      ! salinite
         v0(k)=ytemp(k)/b(k)     ! temperature


!-------------------------
         DO k=2,kmax-1
!-------------------------
           c(k)=-dtm*kz(k,i,j)*unsd3/dsigw(k)/dsigu(k)
           a(k)=-dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
           b(k)=d3*fwet(i,j)+dtm*(kz(k,i,j)*unsd3/dsigw(k)+kz(k-1,i,j)*unsd3/dsigw(k-1))/dsigu(k)





           ysal(k)=d2*fwetp(i,j)*sal(k,i,j)-(dy_u(i,j)*fsalx(k,i,j)-dy_u(i-1,j)*fsalx(k,i-1,j) &
                  +dx_v(i,j)*fsaly(k,i,j)-dx_v(i,j-1)*fsaly(k,i,j-1))/surf(i,j)     &
                  +phisal(k,i,j)/dsigu(k)
           ysal(k)=ysal(k)-(fsalz(k,i,j)-fsalz(k-1,i,j))/dsigu(k)

           ytemp(k)=d2*fwetp(i,j)*temp(k,i,j)-(dy_u(i,j)*ftempx(k,i,j)-dy_u(i-1,j)*ftempx(k,i-1,j) &
                   +dx_v(i,j)*ftempy(k,i,j)-dx_v(i,j-1)*ftempy(k,i,j-1))/surf(i,j)      &
                   +phitemp(k,i,j)/dsigu(k)
           ytemp(k)=ytemp(k)-(ftempz(k,i,j)-ftempz(k-1,i,j))/dsigu(k)



! key_offline


           p0(k)=a(k)*q0(k-1)+b(k)                ! factorisation LU
           q0(k)=-c(k)/p0(k)


           u0(k)=(ysal(k)-a(k)*u0(k-1))/p0(k)     ! salinite
           v0(k)=(ytemp(k)-a(k)*v0(k-1))/p0(k)    ! temperature


!-------------------------
         ENDDO
!-------------------------

!-------------------------
         k=kmax
!-------------------------
         a(k)=-dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
         b(k)=d3*fwet(i,j)+dtm*kz(k-1,i,j)*unsd3/dsigw(k-1)/dsigu(k)
         c(k)=0.0_rlg





         ysal(k)=d2*fwetp(i,j)*sal(k,i,j)-(dy_u(i,j)*fsalx(k,i,j)-dy_u(i-1,j)*fsalx(k,i-1,j) &
                +dx_v(i,j)*fsaly(k,i,j)-dx_v(i,j-1)*fsaly(k,i,j-1))/surf(i,j)     &
                +phisal(k,i,j)/dsigu(k)
         ysal(k)=ysal(k)-(fsalz(k,i,j)-fsalz(k-1,i,j))/dsigu(k)

         ytemp(k)=d2*fwetp(i,j)*temp(k,i,j)-(dy_u(i,j)*ftempx(k,i,j)-dy_u(i-1,j)*ftempx(k,i-1,j) &
                 +dx_v(i,j)*ftempy(k,i,j)-dx_v(i,j-1)*ftempy(k,i,j-1))/surf(i,j)      &
                 +phitemp(k,i,j)/dsigu(k)
         ytemp(k)=ytemp(k)-(ftempz(k,i,j)-ftempz(k-1,i,j))/dsigu(k)



! key_offline


         p0(k)=a(k)*q0(k-1)+b(k)                ! factorisation LU


         u0(k)=(ysal(k)-a(k)*u0(k-1))/p0(k)     ! salinite
         v0(k)=(ytemp(k)-a(k)*v0(k-1))/p0(k)    ! temperature

!-------------------------------------------------
!   resolution des systemes lineaires tridiagonaux
!-------------------------------------------------

         xk=u0(kmax)
         sal(kmax,i,j)=xk
         yk=v0(kmax)
         temp(kmax,i,j)=yk
         DO k=kmax-1,1,-1
           xk=q0(k)*xk+u0(k)
           sal(k,i,j)=xk
           yk=q0(k)*yk+v0(k)
           temp(k,i,j)=yk
         ENDDO

! key_offline



! estimation des traceurs
!------------------------
! refonte avril2016 : all substances
! ========================
! all variables except sand with l_subs2D=TRUE
         DO l=nv_grav+1,nv_adv
           IF(.NOT.l_subs2D(l)) THEN
             DO k=1,kmax
               y(k,l)=d2*fwetp(i,j)*cv_wat(l,k,i,j)-(dy_u(i,j)*fx(l,k,i,j)-dy_u(i-1,j)*fx(l,k,i-1,j) &
                                       +dx_v(i,j)*fy(l,k,i,j)-dx_v(i,j-1)*fy(l,k,i,j-1) &
                                       -phicon(l,k,i,j)/dsigu(k))/surf(i,j)
               y(k,l)=y(k,l)-(fz(l,k,i,j)-fz(l,k-1,i,j))/dsigu(k)

             ENDDO
             

!!!!!!!flux a l interface eau/sediment pour les dissoutes only!!!!!!!!!
             y(1,l)=y(1,l)+typdiss(l)*dtm*flx_s2w(l,i,j)/dsigu(1)

             IF (l<=nvp) THEN ! NO DIFFUSION
              b2d=d3*fwet(i,j)
              DO k=1,kmax
                cv_wat(l,k,i,j)=f90(l)*y(k,l)/b2d
              ENDDO
             ELSE
              u0(1)=y(1,l)/b(1)
              DO k=2,kmax
               u0(k)=(y(k,l)-a(k)*u0(k-1))/p0(k)
              ENDDO
              xk=u0(kmax)
              cv_wat(l,kmax,i,j)=f90(l)*(xk)
              DO k=kmax-1,1,-1
                xk=q0(k)*xk+u0(k)
                cv_wat(l,k,i,j)=f90(l)*(xk)
              ENDDO
             ENDIF
           ELSE
             cv_wat(l,2:kmax,i,j)=0.0_rsh  ! concentration in the upper layers =0
             k=1
             b2d=d3*fwet(i,j)
             y(k,l)=d2*fwetp(i,j)*cv_wat(l,k,i,j)-((dy_u(i,j)*fx(l,k,i,j)-dy_u(i-1,j)*fx(l,k,i-1,j) &
                                  +dx_v(i,j)*fy(l,k,i,j)-dx_v(i,j-1)*fy(l,k,i,j-1)) &
                                  -phicon(l,k,i,j))/surf(i,j) 

             y(k,l)=y(k,l)+flx_s2w(l,i,j)*dtm*typdiss(l)

             u0(k)=y(k,l)/b2d
             cv_wat(l,k,i,j)=f90(l)*u0(k)   ! concentration in the row k=1 is the average concentration
           ENDIF
         ENDDO


       ENDIF ! fin test sur la hauteur
     ENDDO ! fin boucle en i
   ENDDO ! fin boucle en j

!$OMP END DO

   IF (l_tssubadv_over_smalldepth) CALL tssub_eq_hm(xe2,xe3,sal,temp &

                                     ,cv_wat           &

                                     )

!$OMP END PARALLEL



   !PRINT*, 'exit TSSUB_EQ_WCOMPACT.H routine'

  END SUBROUTINE tssub_eq





!  END SUBROUTINE tssub_eq
