!  SUBROUTINE parsub_eq(xe2,xe3        &
!#ifdef 1
!                                      ,cv_wat &
!#endif
!                                      ,sv_part,niter,num_ws)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE parsub_eq  ***
!&E
!&E ** Purpose : framework for the subroutine of tracers transport computation
!&E
!&E ** Description : see parsub_eq_wquick.h or parsub_eq_wcen.h
!&E                  for different possibilities of subroutine body
!&E
!&E ** Called by : stepu, stepv
!&E
!&E ** External calls : parsub_eq_hm
!&E
!&E ** Used ij-arrays : xe2, xe3, u, v
!&E
!&E ** Modified variables :  cv_wat
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-03  (P. Douillet) substances with settling velocities
!&E       !                         issued from tssub_eq
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Arguments

!! * Local declarations

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




  SUBROUTINE parsub_eq(xe2,xe3,cv_wat,setvel,niter,num_ws)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE parsub_eq_wquick  ***
!&E
!&E ** Purpose : other possibility of subroutine body for the computation
!&E              of transport of tracers with fall velocity.
!&E
!&E ** Description : A QUICK scheme is used for depths higher than hm.
!&E                  For water heights lower than hm, a special computation
!&E                  is performed in the subroutine parsub_eq_hm.
!&E                  This subroutine body replaces the default body
!&E
!&E ** Important remark about the settling velocity with the use of key_subs_part_eqsubdt cpp key
!&E              To save computation time, we assume the substances have the same settling velocity
!&E              in a group (since the group and the corresponding sub time step depend on
!&E              the diameter and the density of the particule).
!&E              This assumption is not necessary valid for the biology. If so, modify slightly
!&E              the maximum settling velocity in variable.dat file, to create a new group.
!&E
!&E ** Called by : stepu,stepv
!&E
!&E ** External calls : parsub_advz_quickupwind
!&E
!&E ** Used ij-arrays : cv_wat,h0,xe3,flx_s2w,flx_w2s,flx_w2s_sum,fz
!&E
!&E ** Modified variables : cv_wat
!&E
!&E ** Reference :
!&E
!&E ** History :
!&E       !  2010-03  (P. Douillet) substances with settling velocities
!&E       !                         issued from tssub_eq_wquick.h
!&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       !  2011-02 (B. Thouvenin) correction for key_subs_part_eqsubdt key
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2014-01  (B. Thouvenin - F. Cayocca) introduction of 1
!&E       !  2016-04  (B.Thouvenin) only vertical transport due to settling velocoty
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,     ONLY : ig,id,h0,hm,l_tssubadv_over_smalldepth,halfdt,t,fwet,h0fond  
   USE comvars3d,     ONLY : kz,dsigw,dsigu,sig,qmaxz
   USE comsubstance,  ONLY : fz,nv_ws_max,iv_ws_max,nvp,l_subs2D,nv_adv   

   USE comMUSTANG,  ONLY : flx_s2w,flx_w2s,flx_w2s_sum
   USE comsubstance, ONLY : isand1,isand2

   USE comfluxdo,   ONLY : fzeau_ws
   USE toolgeom,  ONLY   : f_dzu


   USE comsubstance,ONLY : ws3



   IMPLICIT NONE

!! * Arguments
   INTEGER,                                                              INTENT(in)    :: niter,num_ws
   REAL(KIND=rsh),                                                       INTENT(in)    :: setvel
   REAL(kind=rsh),DIMENSION(limin:limax,ljmin:ljmax),                    INTENT(in)    :: xe2
   REAL(kind=rsh),DIMENSION(liminm2:limaxp2,ljminm2:ljmaxp2),            INTENT(in)    :: xe3
   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(inout) :: cv_wat

!! * Local declarations
   INTEGER                                :: i,j,k,l,ll,l1,iv,n,niterij
   REAL(KIND=rsh)                         :: d3,unsd3,d2,htott,wssmax,epss,epssmin
   REAL(KIND=rlg),DIMENSION(kmax,nb_var)  :: y
   REAL(KIND=rlg)                         :: xk,yk,dtm
   REAL(KIND=rlg),DIMENSION(kmax)         :: a,b,c,u0,p0,q0
   
! declarations  parsub_advz

   LOGICAL                               :: decent
   REAL(KIND=rsh)                        :: unsqmaxz
   REAL(KIND=rsh)                        :: yr,h3,qot,dyd,dyg,ya,yb,yc,x1,x2,x3,xx
   REAL(KIND=rsh),DIMENSION(kmax)        :: ap,bp,cp,am,bm,cm

   REAL(KIND=rsh),DIMENSION(1:kmax)      :: ws

! declaration parsub_eqhm
   REAL(KIND=rsh)                                  :: cnouv,d3c3,cmm
   REAL(KIND=rsh),DIMENSION(kmax)                  :: dz


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


!$OMP PARALLEL DEFAULT(none) PRIVATE(d2,d3,unsd3,i,j,k,l,n,            &
!$OMP& iv,ll,l1,dtm,u0,b,y,niterij,ws,h3,cnouv,d3c3,dz,     &

!$OMP& ya,yb,yc,yr,qot,decent,dyd,dyg,x1,x2,x3,xx,ap,bp,cp,am,bm,cm,  &
!$OMP&         unsqmaxz,                                   &

!$OMP& htott,wssmax,epss,epssmin,                &

!$OMP& a,c,p0,q0,xk,yk,cmm) &
!$OMP& SHARED(halfdt,l_tssubadv_over_smalldepth,xe2,xe3,hm,ig,id,h0  &
!$OMP& ,limin,limax,ljmin,ljmax,ljminp2,ljmaxm1,fwet,ws3,kz,h0fond   &
!$OMP& ,liminm1,limaxp2,ljminm1,ljmaxp2,t,nvp,qmaxz,sig,dsigw    &
!$OMP& ,cv_wat,fz,nv_ws_max,iv_ws_max,l_subs2D,fzeau_ws          &

!$OMP& ,flx_s2w,flx_w2s,flx_w2s_sum                               &

!$OMP& ,niter,num_ws,setvel,dsigu)

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljminp2,ljmaxm1
     DO i=MAX0(limin,ig(j)+1),MIN0(limax,id(j)-1)

       htott=h0(i,j)+min(xe3(i,j),xe2(i,j))
       IF (htott > h0fond) THEN
         iv=iv_ws_max(num_ws,1)

         IF (htott > hm .AND. .NOT.l_subs2D(iv) ) THEN
         
! calcul du nombre de pas fractionnaire par maille i,j
!------------------------------------------------------
! pour le groupe de variable num_ws dont les vitesses de chute sont identiques (quelque soit k)
            wssmax=-100._rsh
            epssmin=1000.0_rsh
            DO k=1,kmax
               epss=f_dzu(h0(i,j),xe3(i,j),k,i,j)
               IF( wssmax < ws3(k,iv_ws_max(num_ws,1),i,j)) wssmax=ws3(k,iv_ws_max(num_ws,1),i,j)
               IF( epssmin > epss) epssmin=epss
            ENDDO
            niterij=NINT(halfdt/MIN((0.5_rsh*epssmin/wssmax),halfdt))
            IF(niterij > 5 ) THEN
                write(*,*)i,j,t,iv,htott,epssmin,wssmax,niterij
            ENDIF
 
!------------------------------------------------------------------
! les flux horizontaux ont deja ete calcules comme pour les dissous
! flux verticaux dus a la vitesse de chute calcules au temps final (xe3)
!------------------------------------------------------------------
         

!!!! CALL parsub_advz_quickupwind(xe3,setvel,cv_wat,ll,niter,num_ws)
!! ************  include parsub_advzquickupwind   ****************************************

!!!!SUBROUTINE parsub_advz_quickupwind(xe3,setvel,c,ll,niter,num_ws)
!              modifie en .h pour inclure dans boucle i,j
 
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE parsub_advz_quickupwind  ***
!&E
!&E ** Description : A limiter switch to the QUICK or upwind scheme.
!&E                  No estimate above small depths ( < hm)
!&E                  Used if key_parsub_wquick
!&E ** Purpose : compute vertical fluxes of tracers other than temperature and salinity
!&E              using TVD scheme (quick and upwind decentred) with all-or-nothing limitator
!&E
!&E ** Description : For very low water heights ( < hm), no vertical fluxes are computed
!&E                  because the scheme is degenerated vertically and vertical fluxes are useless.
!&E                  This routine is called only if the key_parsub_wquick CPP key
!&E                  is activated (see parsub_eq_wquick.h).
!&E
!&E ** Important remark about the settling velocity with the use of key_subs_part_eqsubdt cpp key
!&E              To save computation time, we assume the substances have the same settling velocity
!&E              in a group (since the group and the corresponding sub time step depend on
!&E              the diameter and the density of the particule).
!&E              This assumption is not necessary valid for the biology. If so, modify slightly
!&E              the maximum settling velocity in variable.dat file, to create a new group.
!&E
!&E
!&E ** History :
!&E       !  20004-10-04
!&E       !  2008-12  (V. Garnier) shift of sigma values : -1 < sig < 0
!&E       !  2009-04  (A. Mano, V. Garnier) MPI2D
!&E       !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
!&E       !  2011-02 (B. Thouvenin) correction for key_subs_part_eqsubdt key + nb of iterations
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2016-04  (B.Thouvenin) only vertical transport due to settling velocoty
!&E       !  2018-12  (B.Thouvenin) converted in.h include in parsub_eq inside the i,j loop (niterij different for each i,j)
!&E
!&E---------------------------------------------------------------------


          unsqmaxz=1.0_rsh/qmaxz

! get the settling velocity
! -------------------------

! we assume all the substances of a group have the same settling velocity (iv=1 for each group)


          iv=iv_ws_max(num_ws,1)
          ws(:)=-ws3(:,iv,i,j)



!  sea surface height at t+
! --------------------------------------------
          h3=h0(i,j)+xe3(i,j)
          IF (h3 > hm) THEN

! estimate flux of particulate substance
! --------------------------------------
            DO k=1,kmax-1

              fzeau_ws(k,i,j)=(0.5_rsh*(ws(k)+ws(k+1)))*REAL(halfdt,rsh)/REAL(niterij,rsh)

            END DO

! calcul des coef du schema QUICK pour w positif
!-----------------------------------------------

            DO k=2,kmax-1

              x1=sig(k-1)
              x2=sig(k)
              x3=sig(k+1)

              xx=(x1-(x2+x3)/2.0_rsh)*2.0_rsh/(x3-x2)
              ap(k)=1.0_rsh/(1.0_rsh-xx**2)
              bp(k)=(xx-xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
              cp(k)=-(xx+xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
            END DO

! cas w negatif
!-------------------

            DO k=1,kmax-2

              x1=sig(k)
              x2=sig(k+1)
              x3=sig(k+2)

              xx=(x1-(x2+x3)/2.0_rsh)*2.0_rsh/(x2-x3)
              am(k)=1.0_rsh/(1.0_rsh-xx**2)
              bm(k)=-(xx+xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
              cm(k)=(xx-xx**2)/(2.0_rsh*(1.0_rsh-xx**2))
            END DO

! pour faire simple : decentre amont pres des limites
! ---------------------------------------------------

            DO l1=1,nv_ws_max(num_ws)
              iv=iv_ws_max(num_ws,l1)
              fz(iv,0:kmax,i,j)=0.0_rsh

              k=1
              IF (fzeau_ws(k,i,j) > 0.0_rsh) THEN
                fz(iv,k,i,j)=fzeau_ws(k,i,j)*cv_wat(iv,k,i,j)
              ELSE
                fz(iv,k,i,j)=fzeau_ws(k,i,j)*cv_wat(iv,k+1,i,j)
              END IF

              k=kmax-1
              IF (fzeau_ws(k,i,j) > 0.0_rsh) THEN
               fz(iv,k,i,j)=fzeau_ws(k,i,j)*cv_wat(iv,k,i,j)
              ELSE
               fz(iv,k,i,j)=fzeau_ws(k,i,j)*cv_wat(iv,k+1,i,j)
              END IF

              DO k=2,kmax-2
                IF (fzeau_ws(k,i,j) <= 0.0_rsh) THEN
                  yb=cv_wat(iv,k+1,i,j)
                  yc=cv_wat(iv,k,i,j)
                  IF (k <= kmax-2) THEN
                    ya=cv_wat(iv,k+2,i,j)
                  ELSE
                   ya=yb
                  END IF
                  dyd=(yb-yc)/dsigw(k)
                  dyg=(yb-ya)/dsigw(k+1)

                  decent=(dyd*dyg) >= 0.0_rsh
                  IF (.NOT.decent) THEN
                    qot=ABS(dyd/dyg)
                    decent=.TRUE.
                    IF (qot <= qmaxz .AND. qot >= unsqmaxz) decent=.FALSE.
                  END IF

                  IF (decent) THEN
                    fz(iv,k,i,j)=fzeau_ws(k,i,j)*yb
                  ELSE
                    yr=cm(k)*yc+bm(k)*yb+am(k)*ya
                    fz(iv,k,i,j)=fzeau_ws(k,i,j)*yr
                  END IF

                ELSE ! test sur la vitesse cas w >0

                  yb=cv_wat(iv,k,i,j)
                  yc=cv_wat(iv,k+1,i,j)
                  IF (k >= 2) THEN
                    ya=cv_wat(iv,k-1,i,j)
                  ELSE
                    ya=yb
                  ENDIF

                  dyd=(yc-yb)/dsigw(k)
                  dyg=(yb-ya)/dsigw(k-1)

                  decent=(dyd*dyg) <= 0.0_rsh
                  IF (.NOT.decent) THEN
                    qot=ABS(dyd/dyg)
                    decent=.TRUE.
                    IF (qot <= qmaxz .AND. qot >= unsqmaxz) decent=.FALSE.
                  ENDIF

                  IF (decent) THEN
                    fz(iv,k,i,j)=fzeau_ws(k,i,j)*yb
                  ELSE
                    yr=cp(k)*yc+bp(k)*yb+ap(k)*ya
                    fz(iv,k,i,j)=fzeau_ws(k,i,j)*yr
                  END IF

                ENDIF    ! fin de test sur la vitesse
              END DO   ! fin de boucle en k

            END DO   ! fin de boucle sur les traceurs

          ENDIF     ! fin de test sur la hauteur d eau



!! ************END include parsub_advzquickupwind   ****************************************


 
         ELSE
! htott < hm & > h0fond or l_subs2D  => pas de chute (juste depot)
            niterij=1
! pas de flux verticaux pour les variables traitees en 2D (sable par exemple)
            fz(iv,:,i,j)=0.0_rsh
         
         ENDIF
         
! Boucle sur sous pas de temps fractionnaire
!---------------------------------------------
         DO ll=1,niterij

           dtm=halfdt/REAL(niterij,rlg)

           d3=h0(i,j)+REAL(ll  ,rsh)*(xe3(i,j)-xe2(i,j))/REAL(niterij,rsh)+xe2(i,j)
           d2=h0(i,j)+REAL(ll-1  ,rsh)*(xe3(i,j)-xe2(i,j))/REAL(niterij,rsh)+xe2(i,j)

           IF(d3 > hm) THEN

!  resolution 3D : hauteur d eau superieure a hm
!--------------------------------------------------
            unsd3=1.0_rsh/d3

! boucle sur chaque variable du groupe
!---------------------------------------
            DO l1=1,nv_ws_max(num_ws)         !  vitesse de chute,
              iv=iv_ws_max(num_ws,l1)

              IF( l_subs2D(iv)) THEN     ! for variables "sable" or "sand"
              
! variable traitee en 2D
!------------------------
                cv_wat(iv,2:kmax,i,j)=0.0_rsh  ! concentration in the upper layers =0
                k=1
                b(k)=d3*fwet(i,j)+dtm*flx_w2s(iv,i,j)
!y(k,iv)=d2*fwetp(i,j)*cv_wat(iv,k,i,j)+flx_s2w(iv,i,j)*dtm
                y(k,iv)=d3*fwet(i,j)*cv_wat(iv,k,i,j)+flx_s2w(iv,i,j)*dtm
                cv_wat(iv,k,i,j)=y(k,iv)/b(k)
                
              ELSE
! variable traitee en 3D
!------------------------

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

                b(k)=b(k)+dtm*flx_w2s(iv,i,j)/dsigu(k)


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

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

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

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

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

!-------------------------
! estimation des traceurs
!------------------------
                DO k=1,kmax
!y(k,iv)=d2*fwetp(i,j)*cv_wat(iv,k,i,j)-(fz(iv,k,i,j)-fz(iv,k-1,i,j))/dsigu(k)
                  y(k,iv)=d3*fwet(i,j)*cv_wat(iv,k,i,j)-(fz(iv,k,i,j)-fz(iv,k-1,i,j))/dsigu(k)
                ENDDO

                y(1,iv)=y(1,iv)+flx_s2w(iv,i,j)*dtm/dsigu(1)

                u0(1)=y(1,iv)/b(1)
                DO k=2,kmax
                 u0(k)=(y(k,iv)-a(k)*u0(k-1))/p0(k)
                ENDDO
                xk=u0(kmax)
                cv_wat(iv,kmax,i,j)=xk
                DO k=kmax-1,1,-1
                 xk=q0(k)*xk+u0(k)
                 cv_wat(iv,k,i,j)=xk
                ENDDO
              ENDIF


! integration of effectif settling => flx_w2s_sum=kg/m2
              flx_w2s_sum(iv,i,j)=flx_w2s_sum(iv,i,j)+flx_w2s(iv,i,j)*cv_wat(iv,1,i,j)*dtm


            ENDDO   ! end test l1 (iv of the same settling group)
            

           ELSE IF (l_tssubadv_over_smalldepth ) THEN
!  resolution 2D : hauteur d eau infeieure a hm
!--------------------------------------------------
         
!*****************  include parsub_eqhm    ***************************************
  
!  SUBROUTINE parsub_eq_hm(xe3,cv_wat,ll,niter,num_ws)
!  modifie en .h include in parsub_eq.h

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE parsub_eq_hm  ***
!&E
!&E ** Purpose : compute the transport of tracers with fall velocity,
!&E              for meshes with low water heights (< hm).
!&E
!&E ** Description : WARNING : The resulting vertical profiles are homogenous.
!&E
!&E ** Used ij-arrays : h0,xe3
!&E                     flx_s2w,flx_w2s,flx_w2s_sum
!&E                     cv_wat
!&E
!&E ** Modified variables : cv_wat
!&E
!&E ** History :
!&E       !  2010-03  (P. Douillet) substances with settling velocities
!&E       !                         issued from tssub_eq_hm
!&E       !  2010-06  (V.M. Faure, N. Bousserez, M. Baklouti) coupling with eco3m
!&E       !  2011-11  (V. Garnier) Spatial extension of ssh(liminm2:limaxp2,ljminm2:ljmaxp2)
!&E       !  2016-04  (B.Thouvenin) only vertical transport due to settling velocoty
!&E       !  2018-12  (B.Thouvenin) converted in.h include in parsub_eq inside the i,j loop (niterij different for each i,j)
!&E
!&E---------------------------------------------------------------------



            dz(:)=dsigu(:)   !  no unit


!-----------------------------------------------
! integration des flux de traceurs particulaires
!-----------------------------------------------

            DO l1=1,nv_ws_max(num_ws) ! number of variables in the group of settling velocity
                iv=iv_ws_max(num_ws,l1)  ! initialize index of the sustance
! according to the number of settling velocity group
! and the number of variables in this group
                IF( l_subs2D(iv)) THEN     ! for variables "sable" or "sand"


                  cmm=cv_wat(iv,1,i,j)

                ELSE
                  cmm=SUM(cv_wat(iv,:,i,j)*dz(:))
                ENDIF


                d3c3=d3*cmm


                d3c3=d3c3+dtm*flx_s2w(iv,i,j)
                cnouv=d3c3/(d3+dtm*flx_w2s(iv,i,j))
                IF(cnouv < 1.e-10) cnouv=0.0_rsh
                flx_w2s_sum(iv,i,j)=flx_w2s_sum(iv,i,j)+flx_w2s(iv,i,j)*cnouv*dtm

                IF( l_subs2D(iv)) THEN     ! for variables "sable" or "sand"
                 cv_wat(iv,1,i,j)=cnouv
                 cv_wat(iv,2:kmax,i,j)=0.0_rsh
                ELSE
                 cv_wat(iv,:,i,j)=cnouv
                ENDIF

            ENDDO


!!****************end include parsub_eqhm.h  ***************************************************************
 
           ENDIF ! fin test sur la hauteur

         END DO  ! ends loop on subdt
                  
       ENDIF ! fin test sur la h0 (terre ou non)
       
 
     ENDDO ! fin boucle en i
   ENDDO ! fin boucle en j
!$OMP END DO


!$OMP END PARALLEL

  END SUBROUTINE parsub_eq






!  END SUBROUTINE parsub_eq
