  SUBROUTINE tssub_advdify_uqm(sal,temp   &

                                       ,c &

                                           )

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE tssub_advdify_ultimatequickestmacho  ***
!&E
!&E ** Purpose : Estimate of meridional fluxes of temperature and salinity
!&E              (advection and diffusion)
!&E              Advection uses QUICKEST scheme with ULTIMATE limiter and
!&E              its MACHO multidimensional extension.
!&E
!&E ** Description : No horizontal diffusion at lower depths ( < hminkxky)
!&E                  - case of cyclicity in y: correct for eastward propagation,
!&E                    uncorrect for southward propagation (stencil not long enough)
!&E
!&E ** Called by :
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays : htx,hex,hx,hty,hey,hy,fxeau,fyeau,sal,temp,c(cv_wat),hzey,dx,dy_v
!&E                     fsaly,ftempy,fy
!&E
!&E ** Modified variables : fsaly,ftempy,fy
!&E
!&E ** Reference : Leonard B. P., Lock A. P. and M. K. MacVean (1996) :
!&E                Conservative Explicit Unrestricted-Time-Step
!&E                Multidimensional COnstancy-Preserving Advections Schemes.
!&E                Monthly Weather Review, 124, 2588:2605
!&E
!&E ** History :
!&E           !  2008-09 (L. Debreu, T. Duhaut, M. Honnorat) original code
!&E           !  2008-09 (V. Garnier) generalized sigma, correction for southern boundary
!&E           !  2008-10 (V. Garnier) optimization
!&E           !  2009-04 (A. Mano) MPI2D
!&E           !  2010-05 (V. Garnier) OMP
!&E           !  2010-10 (T. Duhaut) introduction of diffusion along isopycnes
!&E           !  2016-12  (B.Thouvenin) version sedim_MUSTANG : only mixsed (mudsand not operationnal)
!&E
!&E---------------------------------------------------------------------





























































!! * Modules used
   USE parameters
   USE comvars2d,    ONLY : jb,jh,hty,hey,hy,h0fond,ky,hminkxky,dy_v,dx,htx,hex,    &
                            ig,id,hx,halfdt,l_obc_cycl_y,dx_v

   USE comfluxdo,    ONLY : fsaly,fyeau,fxeau
   USE cometemp,     ONLY : ftempy
   USE comvarp2d,    ONLY : sal2,temp2

   USE comsubstance, ONLY : nv_adv,nv_grav,nv_tot,fy,c2

   USE comMUSTANG,    ONLY : corfluy


   USE toolmpi,      ONLY : ex_j_rlg,ex_j_rsh,ex_peri_y_rsh

   IMPLICIT NONE

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

   REAL(KIND=rsh),DIMENSION(nb_var,kmax,liminm1:limaxp2,ljminm1:ljmaxp2),INTENT(in) :: c


!! * Local declarations
   INTEGER                   :: i,j,k
   REAL(KIND=rsh), PARAMETER :: dsmall = 1.0d-20
   REAL(KIND=rsh)            :: vv,dyds,dygs,dydt,d0,d1,vvv,courant,uu,  &
                                dygt,kyb,yas,yat,ybs,ybt,ycs,yct,thetap, &
                                psip,deltat,deltas
   REAL(KIND=rlg), DIMENSION(kmax) :: ful,fur
   REAL(KIND=rsh), DIMENSION(kmax) :: sall,salr,templ,tempr

   INTEGER                              :: iv
   REAL(KIND=rsh), DIMENSION(nv_adv)    :: dyd,dyg,ya,yb,yc,delta
   REAL(KIND=rsh), DIMENSION(1:nv_adv,kmax) :: cl,cr


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

! estimate of intermediate tracers (from upwind scheme)
! -----------------------------------------------------

! at open boundaries (north and south)

!$OMP SINGLE

   IF (ljmin==jmin) THEN

     temp2(:,limin:limax,jmin:jmin+1) = temp(:,limin:limax,jmin:jmin+1)
     sal2(:,limin:limax,jmin:jmin+1) = sal(:,limin:limax,jmin:jmin+1)


     c2(nv_grav+1:nv_adv,:,limin:limax,jmin:jmin+1) = c(nv_grav+1:nv_adv,:,limin:limax,jmin:jmin+1)

   ENDIF
   IF (ljmax==jmax) THEN

     temp2(:,limin:limax,jmax) = temp(:,limin:limax,jmax)
     sal2(:,limin:limax,jmax) = sal(:,limin:limax,jmax)


     c2(nv_grav+1:nv_adv,:,limin:limax,jmax) = c(nv_grav+1:nv_adv,:,limin:limax,jmax)

   ENDIF

!$OMP END SINGLE

! inside the domain

!$OMP DO SCHEDULE(RUNTIME),PRIVATE(i,j,k,     &
!$OMP& salr,tempr,fur,sall,templ,ful,         &

!$OMP& iv,cr,cl,                              &

!$OMP& uu)

   DO j=ljmin,ljmax

     DO i=liminp1m1,limax
     IF (i==ig(j) .OR. i==liminp1m1) THEN


       DO k=1,kmax
         IF (htx(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN

           salr(k)=0.0_rsh       ! no flux advection from this cell
           tempr(k)=0.0_rsh


           cr(nv_grav+1:nv_adv,k)=0.0_rsh

           fur(k)=0.0_rlg
         ELSE

           fur(k)=fxeau(k,i,j)   ! h*uz*dt/2

! upwind scheme
           IF (fur(k)<= 0.0_rlg) THEN

             salr(k)=sal(k,i+1,j)
             tempr(k)=temp(k,i+1,j)


             cr(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i+1,j)

           ELSE

             salr(k)=sal(k,i,j)
             tempr(k)=temp(k,i,j)


             cr(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)

           END IF   ! ends test on the direction of propagation

         END IF   ! ends test on the water depth
       END DO   ! ends loop on k


     ELSE IF(i >= ig(j)+1 .AND. i <=id(j)-1 .AND. i>liminp1m1) THEN


       DO k=1,kmax

! left of the stencil

         sall(k)=salr(k)
         templ(k)=tempr(k)


         cl(nv_grav+1:nv_adv,k)=cr(nv_grav+1:nv_adv,k)

         ful(k)=fur(k)

! right of the stencil
         IF (htx(i,j) <= 0.0_rsh .OR. hx(i,j) == -valmanq-h0fond) THEN

           salr(k)=0.0_rsh       ! no flux advection from this cell
           tempr(k)=0.0_rsh


           cr(nv_grav+1:nv_adv,k)=0.0_rsh

           fur(k)=0.0_rlg
         ELSE

           fur(k)=fxeau(k,i,j)   ! h*uz*dt/2

! upwind scheme
           IF (fur(k)<= 0.0_rlg) THEN

             salr(k)=sal(k,i+1,j)
             tempr(k)=temp(k,i+1,j)


             cr(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i+1,j)

           ELSE

             salr(k)=sal(k,i,j)
             tempr(k)=temp(k,i,j)


             cr(nv_grav+1:nv_adv,k)=c(nv_grav+1:nv_adv,k,i,j)

           END IF   ! ends test of the direction of propagation

         END IF   ! ends test on the water depth

! estimate of intermediate tracers (from upwind scheme)

         IF ( (ful(k)>0.0_rlg) .AND. (fur(k)>0.0_rlg) ) THEN

           uu = ful(k)/(halfdt*hex(i-1,j))

         ELSE IF ( (ful(k)<0.0_rlg) .AND. (fur(k)<0.0_rlg) ) THEN

           uu = fur(k)/(halfdt*hex(i,j))

         ELSE
! no time evolution of the intermediate temperature
! along the boundaries
           uu = 0.0_rsh
         END IF

! apply euler (time) - upwind (x) scheme

         temp2(k,i,j) = temp(k,i,j) - (halfdt/dx(i,j))*uu*(tempr(k)-templ(k))
         sal2(k,i,j)  = sal(k,i,j)  - (halfdt/dx(i,j))*uu*(salr(k)-sall(k))  


         c2(nv_grav+1:nv_adv,k,i,j)  = c(nv_grav+1:nv_adv,k,i,j)  &
           - (halfdt/dx(i,j))*uu*(cr(nv_grav+1:nv_adv,k)-cl(nv_grav+1:nv_adv,k))  


       END DO   ! ends loop on k

     END IF

     END DO   ! ends loop on i
   END DO   ! ends loop on j
!$OMP END DO

!$OMP SINGLE
! periodic conditions along y-axis
   IF (l_obc_cycl_y) THEN



     CALL ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
     CALL ex_peri_y_rsh(1,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))

     CALL ex_peri_y_rsh(1,(nv_adv)*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(1:nv_adv,:,liminm1:limaxp2,ljminm1:ljmaxp2))




     END IF  ! ends periodic conditions along y-axis
!$OMP END SINGLE

!$OMP barrier
!$OMP master

   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal2(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp2(:,liminm1:limaxp2,ljminm1:ljmaxp2))


!to avoid temporary copy 1st dimension is not contiguous   CALL ex_j_rsh(-1,2,(nv_adv)*kmax,limin,limax,ljminm1,ljmaxp2,c2(1:nv_adv,:,limin:limax,ljminm1:ljmaxp2))
   CALL ex_j_rsh(-1,2,nv_tot*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,c2(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))

!$OMP end master
!$OMP barrier

!$OMP flush(sal2,temp2)


!$OMP flush(c2)


! estimate of fluxes (from ultimate quickest scheme)
! --------------------------------------------------

!$OMP DO schedule(runtime),private(i,j,k,dyds,dygs,dydt    &
!$OMP&,dygt,kyb,yas,yat,ybs,ybt,ycs,yct,deltat,deltas      &

!$OMP& ,iv,dyd,dyg,ya,yb,yc,delta                          &

!$OMP&,d0,d1,vv,vvv,courant,thetap,psip)

   DO j=ljmin,ljmax
     DO i=liminp2,limaxm1

       IF (j >= jb(i) .AND. j <= jh(i)-1) THEN
         DO k=1,kmax
           IF (hty(i,j) <= 0.0_rsh .OR. hy(i,j)==-valmanq-h0fond) THEN

             fsaly(k,i,j)=0.0_rlg
             ftempy(k,i,j)=0.0_rlg


             fy(nv_grav+1:nv_adv,k,i,j)=0.0_rsh

           ELSE
             kyb=ky
             IF (hty(i,j) < hminkxky) kyb=0.0_rsh ! no diffusion for small water depth
             vv=fyeau(k,i,j)

             vvv = vv/(halfdt*hey(i,j))

             courant = ABS(vvv)*halfdt/dy_v(i,j)
             d0=(1.0_rsh/6.0_rsh)*(2.0_rsh-courant)*(1.0_rsh-courant)
             d1=(1.0_rsh/6.0_rsh)*(1.0_rsh-courant)*(1.0_rsh+courant)
                           
             IF (vv <= 0.0_rsh) THEN  ! southward advection

               ybs=sal2(k,i,j+1)
               ycs=sal2(k,i,j)
               ybt=temp2(k,i,j+1)
               yct=temp2(k,i,j)


               yb(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j+1)
               yc(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j)


               IF (hty(i,j+1) > 0.0_rsh .AND. (j+1 < jh(i) .OR. l_obc_cycl_y)) THEN

                 yas=sal2(k,i,j+2)
                 yat=temp2(k,i,j+2)


                 ya(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j+2)

               ELSE

                 yas=ybs
                 yat=ybt


                 ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)

               END IF


               dyds=ycs-ybs
               dygs=ybs-yas
               dydt=yct-ybt
               dygt=ybt-yat
               deltas=ybs-ycs
               deltat=ybt-yct


               dyg(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-ya(nv_grav+1:nv_adv)
               dyd(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
               delta(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-yc(nv_grav+1:nv_adv)


             ELSE      ! northward advection


               ybs=sal2(k,i,j)
               ycs=sal2(k,i,j+1)
               ybt=temp2(k,i,j)
               yct=temp2(k,i,j+1)


               yb(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j)
               yc(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j+1)

               IF (hty(i,j-1) > 0.0_rsh .AND. (j > jb(i) .OR. l_obc_cycl_y)) THEN

                 yas=sal2(k,i,j-1)
                 yat=temp2(k,i,j-1)


                 ya(nv_grav+1:nv_adv)=c2(nv_grav+1:nv_adv,k,i,j-1)

               ELSE

                 yas=ybs
                 yat=ybt


                 ya(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)

               END IF


               dyds=ycs-ybs
               dygs=ybs-yas
               dydt=yct-ybt
               dygt=ybt-yat
               deltas=ycs-ybs
               deltat=yct-ybt


               dyg(nv_grav+1:nv_adv)=yb(nv_grav+1:nv_adv)-ya(nv_grav+1:nv_adv)
               dyd(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)
               delta(nv_grav+1:nv_adv)=yc(nv_grav+1:nv_adv)-yb(nv_grav+1:nv_adv)


             END IF       ! ends test on velocity sign



! flux of salinity advection
! --------------------------

             IF ((ABS(dyds).LE.dsmall).OR.(ABS(dygs).LE.dsmall)) THEN
! no horizontal gradient
               thetap = 0.0_rsh
               psip = 0.0_rsh
             ELSE
! ultimate limiter
               thetap = (dygs+dsmall)/(dyds+dsmall)
               psip = d0+d1*thetap
               psip=MAX(0.0_rsh,MIN(MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap))
             END IF
             fsaly(k,i,j) = vv*(ybs+psip*dyds)

! flux of temperature advection
! -----------------------------

             IF ((ABS(dydt).LE.dsmall).OR.(ABS(dygt).LE.dsmall)) THEN
! no horizontal gradient
               thetap = 0.0_rsh
               psip = 0.0_rsh
             ELSE
! ultimate limiter
               thetap = (dygt+dsmall)/(dydt+dsmall)
               psip = d0+d1*thetap
               psip=MAX(0.0_rsh,MIN(MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap))
             END IF
             ftempy(k,i,j) = vv*(ybt+psip*dydt)

! key_offline




! advection flux of other substances
! ----------------------------------

             DO iv=nv_grav+1,nv_adv
               IF ((ABS(dyd(iv)).LE.dsmall).OR.(ABS(dyg(iv)).LE.dsmall)) THEN
! no horizontal gradient
                 thetap = 0.0_rsh
                 psip = 0.0_rsh
               ELSE
! ultimate limiter
                 thetap = (dyg(iv)+dsmall)/(dyd(iv)+dsmall)
                 psip = d0+d1*thetap
                 psip=MAX(0.0_rsh,MIN(MIN(1.0_rsh,psip),((1.0_rsh-courant)/(courant+dsmall))*thetap))
               END IF
               fy(iv,k,i,j) = vv*(yb(iv)+psip*dyd(iv))
             END DO





! Advection/diffusion term (iso-sigma)
! ------------------------------------

             ftempy(k,i,j)=ftempy(k,i,j)-hey(i,j)*kyb*deltat/dy_v(i,j)*halfdt
             fsaly(k,i,j)=fsaly(k,i,j)-hey(i,j)*kyb*deltas/dy_v(i,j)*halfdt              




             DO iv=nv_grav+1,nv_adv


               fy(iv,k,i,j) = fy(iv,k,i,j) - hey(i,j)*kyb*delta(iv)/dy_v(i,j)*halfdt



               IF (k.eq.1) fy(iv,k,i,j)=fy(iv,k,i,j)*corfluy(iv,i,j)


             END DO


           END IF ! ends test on water depth
         END DO ! ends loop on k
       END IF ! for selecting right j loop
     END DO ! ends loop on i
   END DO ! ends loop on j
!$OMP END DO

!$OMP barrier
!$OMP master

   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fsaly(:,limin:limax,ljminm1:ljmax))
   CALL ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,ftempy(:,limin:limax,ljminm1:ljmax))


   CALL ex_j_rsh(-1,0,nv_adv*kmax,limin,limax,ljminm1,ljmax,fy(:,:,limin:limax,ljminm1:ljmax))

!$OMP end master
!$OMP barrier

!$OMP flush(fsaly,ftempy)


!$OMP flush(fy)


  END SUBROUTINE tssub_advdify_uqm
