 MODULE tssubeqtest
   !!======================================================================
   !!                   ***  MODULE  tssubeqtest  ***
   !! Tracers advection : no dynamics (momentum,continuity)
   !!======================================================================

#include "toolcpp.h"

   !! * Modules used
   USE parameters, ONLY : imin,imax,jmin,jmax,kmax,rsh,        &
                          ljmin,ljmax,ljmaxp1,ljmaxp2,ljminm1, &
                          limin,limax,limaxp1,limaxp2,liminm1, &
                          liminm2,ljminm2
#ifdef key_substance
   USE parameters, ONLY : nb_var
#endif

   IMPLICIT NONE

   PRIVATE

   !! * Accessibility
   PUBLIC tssub_eq_test                 ! routine called by step

   !! * Shared or public module variables
#ifdef key_tssub_casadvsmolartime
   REAL(KIND=rsh), DIMENSION(imin-1:imax,jmin-1:jmax), PUBLIC :: PSISM
#elif key_tssub_casadvlense
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(:,:,:), PUBLIC :: uzc0,vzc0,wzc0,wzc0_phys
   REAL(KIND=rsh), ALLOCATABLE, DIMENSION(  :,:), PUBLIC :: uc0,vc0
#endif

   !! * Private variables

 CONTAINS
   !!======================================================================
 
  SUBROUTINE tssub_eq_test
 
   !&E---------------------------------------------------------------------
   !&E                 ***  ROUTINE tssub_eq_test  ***
   !&E
   !&E ** Purpose : 
   !&E
   !&E ** Description :
   !&E
   !&E ** Called by :
   !&E
   !&E ** External calls : 
   !&E
   !&E ** Reference :
   !&E
   !&E ** History :
   !&E       !  2008-09  (V. Garnier)  Original code
   !&E       !  2010-04  (T. Odaka) MPI2D cyclic
   !&E
   !&E---------------------------------------------------------------------
   !! * Modules used
   USE parameters,ONLY :  ishift_obc_west,ishift_obc_east, &
                          jshift_obc_south,jshift_obc_north
   USE comvarp2d, ONLY : u,v,ssh
   USE comvarp3d, ONLY : uz,vz
   USE comvarp2d, ONLY : sal,temp
   USE comvars2d, ONLY : l_obc_cycl_x,l_obc_cycl_y,l_modele2d,hey,jb,jh,id,ig
#if defined key_tssub_casadvsmolartime || defined key_tssub_casadvlense
   USE parameters,ONLY : rsh,rlg
   USE comvars2d, ONLY : dt,nbouc,dx,dy,tfin,hex,hey,halfdt,h0
   USE comfluxdo, ONLY : fxeau,fyeau,fzeau
   USE comvars3d, ONLY : wz
#ifdef key_diag_physwz
   USE comvars3d, ONLY : wz_phys
#endif
#if defined key_siggen || defined key_gencoord
   USE comsiggen, ONLY : hzex,hzey
   USE toolgeom,  ONLY : f_hzw
#endif
#endif
#ifdef key_substance
   USE comsubstance, ONLY : cv_wat,nv_adv,nvp,nb_ws_max,itypv,iv_ws_max,sv_part,ndt_part
#endif
   USE_MPI toolmpi,   ONLY : ex_i_rsh,ex_j_rsh,ex_i_rlg,ex_j_rlg ,ex_peri_x_rsh,ex_peri_y_rsh

   !! * Declaration function

   !! * Arguments

   !! * Local declarations
   INTEGER             :: i,j,k
#if defined key_tssub_casadvsmolartime || defined key_tssub_casadvlense
   REAL(KIND=rsh)      :: phit
#endif
#if defined key_substance
   INTEGER             :: l
#endif

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

#if defined key_tssub_casadvsmolartime || defined key_tssub_casadvlense

!$OMP PARALLEL DEFAULT(NONE) PRIVATE(phit,i,j,k)                         &
!$OMP& SHARED (dx,dy,halfdt,u,v,uz,vz,hex,hey,nbouc,dt,tfin,             &
#ifdef key_tssub_casadvsmolartime
!$OMP&         PSISM,                                                    &
#elif key_tssub_casadvlense
!$OMP&         uc0,vc0,uzc0,vzc0,wzc0,wz,wzc0_phys,fzeau,h0,ssh,         &
#endif
#ifdef key_diag_physwz
!$OMP&         wz_phys,                                                  &
#endif
!$OMP&         limin,limax,ljmin,ljmax,fxeau,fyeau)
   ! Initialization of velocities and water fluxes
   !----------------------------------------------
   ! phit should be cos ( 2 pi t / T ) to do a round trip
   ! with phit equal to cos ( pi t / T ), one do a back and forth
#ifdef key_tssub_casadvsmolartime 
   phit = COS(4.0_rlg*ATAN(1.0_rlg)*REAL(nbouc,rlg)*dt/tfin)
!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=limin,limax
       u(i,j) = -(PSISM(i,j)-PSISM(i,j-1))/dy(i,j)*phit
       v(i,j) =  (PSISM(i,j)-PSISM(i-1,j))/dx(i,j)*phit
       DO k=1,kmax
         uz(k,i,j) = u(i,j)
         vz(k,i,j) = v(i,j)
#if defined key_siggen || defined key_gencoord
         fxeau(k,i,j)=hzex(k,i,j)*uz(k,i,j)*halfdt
         fyeau(k,i,j)=hzey(k,i,j)*vz(k,i,j)*halfdt
#else
         fxeau(k,i,j)=hex(i,j)*uz(k,i,j)*halfdt
         fyeau(k,i,j)=hey(i,j)*vz(k,i,j)*halfdt
#endif
       END DO
     END DO
   END DO
!$OMP END DO

#elif key_tssub_casadvlense
   phit = SIN(4.0_rlg*ATAN(1.0_rlg)*REAL(nbouc,rlg)*dt/tfin)
!$OMP SINGLE
   u(:,:)=uc0(:,:)*phit
   v(:,:)=vc0(:,:)*phit
   uz(:,:,:)=uzc0(:,:,:)*phit
   vz(:,:,:)=vzc0(:,:,:)*phit
   wz(:,:,:)=wzc0(:,:,:)*phit
#ifdef key_diag_physwz
   wz_phys(:,:,:)=wzc0_phys(:,:,:)*phit
#endif
!$OMP END SINGLE

!$OMP DO SCHEDULE(RUNTIME)
   DO j=ljmin,ljmax
     DO i=limin,limax
#if defined key_siggen || defined key_gencoord
         fxeau(:,i,j)=hzex(:,i,j)*uz(:,i,j)*halfdt
         fyeau(:,i,j)=hzey(:,i,j)*vz(:,i,j)*halfdt
         DO k=1,kmax
           fzeau(k,i,j)=f_hzw(h0(i,j),ssh(i,j),k,i,j)*wz(k,i,j)*halfdt
         END DO
#else
         fxeau(:,i,j)=hex(i,j)*uz(:,i,j)*halfdt
         fyeau(:,i,j)=hey(i,j)*vz(:,i,j)*halfdt
         fzeau(1:kmax,i,j)=(h0(i,j)+ssh(i,j))*wz(1:kmax,i,j)*halfdt
#endif
     END DO
   END DO
!$OMP END DO
#endif

!$OMP END PARALLEL

   CALL_MPI ex_i_rlg(-1,0,kmax,liminm1,limax,ljmin,ljmax,fxeau(:,liminm1:limax,ljmin:ljmax))
   CALL_MPI ex_j_rlg(-1,0,kmax,limin,limax,ljminm1,ljmax,fyeau(:,limin:limax,ljminm1:ljmax))

#endif

   ! advection of tracers
   ! --------------------

#ifdef key_substance
   CALL_MPI ex_i_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,nb_var*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_i_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
   CALL_MPI ex_j_rsh(-1,2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))

   IF (l_modele2d) THEN
     CALL tssub2d_eq(                                             &
#ifdef key_substance
                     cv_wat,                                      &
#endif
                     sal,temp,ssh(limin:limax,ljmin:ljmax),ssh )
   ELSE
     CALL tssub_eq(ssh(limin:limax,ljmin:ljmax),                  &
                   ssh,sal,temp                                   &
#ifdef key_substance
                  ,cv_wat                                         &
#endif
                  )
     ! Solve equation of variables with settling velocity
#ifdef key_substance
     IF (nvp>0) THEN
       DO l=1,nb_ws_max
         ! not for gravel
         IF(itypv(iv_ws_max(l,1)) > 1 ) THEN
              CALL parsub_eq(ssh,cv_wat,sv_part(l),ndt_part(l),l)
         ENDIF
       END DO
     END IF
#endif
   END IF

   ! open boundary
   ! -------------
   IF (l_obc_cycl_x) THEN
#ifdef key_MPI_2D
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
         CALL_MPI ex_peri_x_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
         CALL_MPI ex_peri_x_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
         CALL_MPI ex_peri_x_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#ifdef key_substance
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         CALL_MPI ex_peri_x_rsh(3,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#else
         CALL_MPI ex_peri_x_rsh(2,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#endif
#else
     DO j=ljmin,ljmax
         temp(:,ig(j),j)=temp(:,id(j)-1,j)
         temp(:,ig(j)-1,j)=temp(:,id(j)-2,j)
         temp(:,id(j),j)=temp(:,ig(j)+1,j)
         sal(:,ig(j),j)=sal(:,id(j)-1,j)
         sal(:,ig(j)-1,j)=sal(:,id(j)-2,j)
         sal(:,id(j),j)=sal(:,ig(j)+1,j)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         temp(:,ig(j)-2,j)=temp(:,id(j)-3,j)
         temp(:,id(j)+1,j)=temp(:,ig(j)+2,j)
         temp(:,id(j)+2,j)=temp(:,ig(j)+3,j)
         sal(:,ig(j)-2,j)=sal(:,id(j)-3,j)
         sal(:,id(j)+1,j)=sal(:,ig(j)+2,j)
         sal(:,id(j)+2,j)=sal(:,ig(j)+3,j)
#endif
#ifdef key_substance
         cv_wat(:,:,ig(j),j)=cv_wat(:,:,id(j)-1,j)
         cv_wat(:,:,ig(j)-1,j)=cv_wat(:,:,id(j)-2,j)
         cv_wat(:,:,id(j),j)=cv_wat(:,:,ig(j)+1,j)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         cv_wat(:,:,ig(j)-2,j)=cv_wat(:,:,id(j)-3,j)
         cv_wat(:,:,id(j)+1,j)=cv_wat(:,:,ig(j)+2,j)
         cv_wat(:,:,id(j)+2,j)=cv_wat(:,:,ig(j)+3,j)
#endif
#endif
     END DO
#endif
   ELSE
     IF_MPI (limin==imin) THEN
         temp(:,imin+1+ishift_obc_west,ljmin:ljmax)=temp(:,imin+2+ishift_obc_west,ljmin:ljmax)
         temp(:,imin+ishift_obc_west,ljmin:ljmax)=temp(:,imin+1+ishift_obc_west,ljmin:ljmax)
         sal(:,imin+1+ishift_obc_west,ljmin:ljmax)=sal(:,imin+2+ishift_obc_west,ljmin:ljmax)
         sal(:,imin+ishift_obc_west,ljmin:ljmax)=sal(:,imin+1+ishift_obc_west,ljmin:ljmax)
#ifdef key_substance
         cv_wat(:,:,imin+1+ishift_obc_west,ljmin:ljmax)=cv_wat(:,:,imin+2+ishift_obc_west,ljmin:ljmax)
         cv_wat(:,:,imin+ishift_obc_west,ljmin:ljmax)=cv_wat(:,:,imin+1+ishift_obc_west,ljmin:ljmax)
#endif
     ENDIF_MPI
     IF_MPI (limax==imax) THEN
         temp(:,imax-ishift_obc_east,ljmin:ljmax)=temp(:,imax-1-ishift_obc_east,ljmin:ljmax)
         sal(:,imax-ishift_obc_east,ljmin:ljmax)=sal(:,imax-1-ishift_obc_east,ljmin:ljmax)
#ifdef key_substance
         cv_wat(:,:,imax-ishift_obc_east,ljmin:ljmax)=cv_wat(:,:,imax-1-ishift_obc_east,ljmin:ljmax)
#endif
     ENDIF_MPI
   END IF
   IF (l_obc_cycl_y) THEN
#ifdef key_MPI_2D
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
         CALL_MPI ex_peri_y_rsh(3,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
         CALL_MPI ex_peri_y_rsh(3,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#else
         CALL_MPI ex_peri_y_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,sal(:,liminm1:limaxp2,ljminm1:ljmaxp2))
         CALL_MPI ex_peri_y_rsh(2,kmax,liminm1,limaxp2,ljminm1,ljmaxp2,temp(:,liminm1:limaxp2,ljminm1:ljmaxp2))
#ifdef key_substance
         CALL_MPI ex_peri_y_rsh(2,nv_adv*kmax,liminm1,limaxp2,ljminm1,ljmaxp2,cv_wat(:,:,liminm1:limaxp2,ljminm1:ljmaxp2))
#endif
#endif
#else
     DO i=limin,limax
         temp(:,i,jb(i))=temp(:,i,jh(i)-1)
         temp(:,i,jb(i)-1)=temp(:,i,jh(i)-2)
         temp(:,i,jh(i))=temp(:,i,jb(i)+1)
         sal(:,i,jb(i))=sal(:,i,jh(i)-1)
         sal(:,i,jb(i)-1)=sal(:,i,jh(i)-2)
         sal(:,i,jh(i))=sal(:,i,jb(i)+1)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         temp(:,i,jb(i)-2)=temp(:,i,jh(i)-3)
         temp(:,i,jh(i)+1)=temp(:,i,jb(i)+2)
         temp(:,i,jh(i)+2)=temp(:,i,jb(i)+3)
         sal(:,i,jb(i)-2)=sal(:,i,jh(i)-3)
         sal(:,i,jh(i)+1)=sal(:,i,jb(i)+2)
         sal(:,i,jh(i)+2)=sal(:,i,jb(i)+3)
#endif
#ifdef key_substance
         cv_wat(:,:,i,jb(i))=cv_wat(:,:,i,jh(i)-1)
         cv_wat(:,:,i,jb(i)-1)=cv_wat(:,:,i,jh(i)-2)
         cv_wat(:,:,i,jh(i))=cv_wat(:,:,i,jb(i)+1)
!!!#if defined key_tssub_adv_fourthorder  || defined key_tssub_adv_fifthorder
#ifdef key_shiftobc
         cv_wat(:,:,i,jb(i)-2)=cv_wat(:,:,i,jh(i)-3)
         cv_wat(:,:,i,jh(i)+1)=cv_wat(:,:,i,jb(i)+2)
         cv_wat(:,:,i,jh(i)+2)=cv_wat(:,:,i,jb(i)+3)
#endif
#endif
     END DO
#endif
   ELSE
     IF_MPI (ljmin==jmin) THEN
         temp(:,limin:limax,jmin+1+jshift_obc_south)=temp(:,limin:limax,jmin+2+jshift_obc_south)
         temp(:,limin:limax,jmin+jshift_obc_south)=temp(:,limin:limax,jmin+1+jshift_obc_south)
         sal(:,limin:limax,jmin+1+jshift_obc_south)=sal(:,limin:limax,jmin+2+jshift_obc_south)
         sal(:,limin:limax,jmin+jshift_obc_south)=sal(:,limin:limax,jmin+1+jshift_obc_south)
#ifdef key_substance
         cv_wat(:,:,limin:limax,jmin+1+jshift_obc_south)=cv_wat(:,:,limin:limax,jmin+2+jshift_obc_south)
         cv_wat(:,:,limin:limax,jmin+jshift_obc_south)=cv_wat(:,:,limin:limax,jmin+1+jshift_obc_south)
#endif
     ENDIF_MPI
     IF_MPI (ljmax==jmax) THEN
         temp(:,limin:limax,jmax-jshift_obc_north)=temp(:,limin:limax,jmax-1-jshift_obc_north)
         sal(:,limin:limax,jmax-jshift_obc_north)=sal(:,limin:limax,jmax-1-jshift_obc_north)
#ifdef key_substance
         cv_wat(:,:,limin:limax,jmax-jshift_obc_north)=cv_wat(:,:,limin:limax,jmax-1-jshift_obc_north)
#endif
     ENDIF_MPI
   END IF

  END SUBROUTINE tssub_eq_test

 END MODULE
